1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // imagery_tools //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // textural_features.cpp //
14 // //
15 // Copyright (C) 2016 by //
16 // Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.org //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "textural_features.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
59 enum
60 {
61 ASM = 0,
62 CONTRAST,
63 CORRELATION,
64 VARIANCE,
65 IDM,
66 SUM_AVERAGE,
67 SUM_ENTROPY,
68 SUM_VARIANCE,
69 ENTROPY,
70 DIF_VARIANCE,
71 DIF_ENTROPY,
72 MOC_1,
73 MOC_2,
74 g_nFeatures
75 };
76
77 //---------------------------------------------------------
78 const CSG_String g_Features[g_nFeatures][2] =
79 {
80 { "ASM" , _TL("Angular Second Moment" ) },
81 { "CONTRAST" , _TL("Contrast" ) },
82 { "CORRELATION" , _TL("Correlation" ) },
83 { "VARIANCE" , _TL("Variance" ) },
84 { "IDM" , _TL("Inverse Diff Moment" ) },
85 { "SUM_AVERAGE" , _TL("Sum Average" ) },
86 { "SUM_ENTROPY" , _TL("Sum Entropy" ) },
87 { "SUM_VARIANCE", _TL("Sum Variance" ) },
88 { "ENTROPY" , _TL("Entropy" ) },
89 { "DIF_VARIANCE", _TL("Difference Variance" ) },
90 { "DIF_ENTROPY" , _TL("Difference Entropy" ) },
91 { "MOC_1" , _TL("Measure of Correlation-1") },
92 { "MOC_2" , _TL("Measure of Correlation-2") }
93 };
94
95
96 ///////////////////////////////////////////////////////////
97 // //
98 // //
99 // //
100 ///////////////////////////////////////////////////////////
101
102 //---------------------------------------------------------
CTextural_Features(void)103 CTextural_Features::CTextural_Features(void)
104 {
105 Set_Name (_TL("Textural Features"));
106
107 Set_Author ("O.Conrad (c) 2016");
108
109 Set_Description (_TW(
110 "Textural features. This tool is based on the GRASS GIS implementation by Carmine Basco (r.texture). "
111 ));
112
113 Add_Reference("Haralick, R.M.; Shanmugam, K.; Dinstein, I.", "1973",
114 "Textural Features for Image Classification",
115 "IEEE Transactions on Systems, Man, and Cybernetics. SMC-3 (6): 610�621.",
116 SG_T("http://haralick.org/journals/TexturalFeatures.pdf"), SG_T("online")
117 );
118
119 Parameters.Add_Grid("",
120 "GRID" , _TL("Grid"),
121 _TL(""),
122 PARAMETER_INPUT
123 );
124
125 for(int i=0; i<g_nFeatures; i++)
126 {
127 Parameters.Add_Grid("", g_Features[i][0], g_Features[i][1], _TL(""), PARAMETER_OUTPUT_OPTIONAL);
128 }
129
130 Parameters.Add_Choice("",
131 "DIRECTION" , _TL("Direction"),
132 _TL(""),
133 CSG_String::Format("%s|%s|%s|%s|%s",
134 _TL("all"),
135 _TL("N-S"),
136 _TL("NE-SW"),
137 _TL("E-W"),
138 _TL("SE-NW")
139 ), 0
140 );
141
142 Parameters.Add_Int("",
143 "RADIUS" , _TL("Radius"),
144 _TL("kernel radius in cells"),
145 1, 1, true
146 );
147
148 Parameters.Add_Int("",
149 "DISTANCE" , _TL("Distance"),
150 _TL("The distance between two samples."),
151 1, 1, true
152 );
153
154 Parameters.Add_Int("",
155 "MAX_CATS" , _TL("Maximum Number of Categories"),
156 _TL(""),
157 256, 2, true
158 );
159 }
160
161
162 ///////////////////////////////////////////////////////////
163 // //
164 ///////////////////////////////////////////////////////////
165
166 //---------------------------------------------------------
On_Execute(void)167 bool CTextural_Features::On_Execute(void)
168 {
169 //-----------------------------------------------------
170 CSG_Grid *pFeatures[g_nFeatures];
171
172 {
173 int i, n = 0;
174
175 for(i=0; i<g_nFeatures; i++)
176 {
177 if( (pFeatures[i] = Parameters(g_Features[i][0])->asGrid()) != NULL )
178 {
179 n++;
180 }
181 }
182
183 if( n == 0 )
184 {
185 Error_Set(_TL("Nothing to do. No feature has been selected."));
186
187 return( false );
188 }
189 }
190
191 //-----------------------------------------------------
192 m_pGrid = Parameters("GRID")->asGrid();
193
194 if( m_pGrid->Get_Range() <= 0.0 )
195 {
196 Error_Set(_TL("Nothing to do. No variation in input grid."));
197
198 return( false );
199 }
200
201 //-----------------------------------------------------
202 m_Radius = Parameters("RADIUS" )->asInt();
203 m_MaxCats = Parameters("MAX_CATS")->asInt();
204
205 int Distance = Parameters("DISTANCE" )->asInt();
206 int Direction = Parameters("DIRECTION")->asInt();
207
208 //-----------------------------------------------------
209 for(int y=0; y<Get_NY() && Set_Progress(y); y++)
210 {
211 #pragma omp parallel for
212 for(int x=0; x<Get_NX(); x++)
213 {
214 CSG_Matrix P[4];
215
216 if( !Get_Matrices(x, y, Distance, P) )
217 {
218 for(int i=0; i<g_nFeatures; i++)
219 {
220 if( pFeatures[i] )
221 {
222 pFeatures[i]->Set_NoData(x, y);
223 }
224 }
225 }
226 else
227 {
228 CSG_Vector Features(g_nFeatures);
229
230 switch( Direction )
231 {
232 default: // all
233 Get_Features(Features, P[0]);
234 Get_Features(Features, P[1]);
235 Get_Features(Features, P[2]);
236 Get_Features(Features, P[3]);
237 break;
238
239 case 1: // N-S
240 Get_Features(Features, P[0]);
241 break;
242
243 case 2: // NE-SW
244 Get_Features(Features, P[1]);
245 break;
246
247 case 3: // E-W
248 Get_Features(Features, P[2]);
249 break;
250
251 case 4: // SE-NW
252 Get_Features(Features, P[3]);
253 break;
254 }
255
256 for(int i=0; i<g_nFeatures; i++)
257 {
258 if( pFeatures[i] )
259 {
260 pFeatures[i]->Set_Value(x, y, Direction ? Features[i] : Features[i] / 4);
261 }
262 }
263 }
264 }
265 }
266
267 //-----------------------------------------------------
268 return( true );
269 }
270
271
272 ///////////////////////////////////////////////////////////
273 // //
274 ///////////////////////////////////////////////////////////
275
276 //---------------------------------------------------------
Get_Value(int x,int y)277 inline int CTextural_Features::Get_Value(int x, int y)
278 {
279 if( m_pGrid->is_InGrid(x, y) )
280 {
281 return( (int)((m_pGrid->asDouble(x, y) - m_pGrid->Get_Min()) * (m_MaxCats - 1) / m_pGrid->Get_Range()) );
282 }
283
284 return( -1 );
285 }
286
287 //---------------------------------------------------------
Get_Matrices(int x,int y,int d,CSG_Matrix P[4])288 bool CTextural_Features::Get_Matrices(int x, int y, int d, CSG_Matrix P[4])
289 {
290 if( m_pGrid->is_NoData(x, y) )
291 {
292 return( false );
293 }
294
295 size_t iTone, nTones;
296
297 int ix, iy;
298
299 //-----------------------------------------------------
300 CSG_Array_Int Tones(m_MaxCats);
301
302 for(iTone=0; iTone<Tones.Get_Size(); iTone++)
303 {
304 Tones[iTone] = -1;
305 }
306
307 // Determine the number of different gray scales (not maxval)
308 for(iy=y-m_Radius; iy<=y+m_Radius; iy++)
309 {
310 for(ix=x-m_Radius; ix<=x+m_Radius; ix++)
311 {
312 int Value = Get_Value(ix, iy);
313
314 if( Value < 0 )
315 {
316 return( false );
317 }
318
319 Tones[Value] = Value;
320 }
321 }
322
323 // Collapse array, taking out all zero values
324 for(iTone=0, nTones=0; iTone<Tones.Get_Size(); iTone++)
325 {
326 if( Tones[iTone] >= 0 )
327 {
328 Tones[nTones++] = Tones[iTone];
329 }
330 }
331
332 Tones.Set_Array(nTones);
333
334 //-----------------------------------------------------
335 P[0].Create((int)nTones, (int)nTones);
336 P[1].Create((int)nTones, (int)nTones);
337 P[2].Create((int)nTones, (int)nTones);
338 P[3].Create((int)nTones, (int)nTones);
339
340 //-----------------------------------------------------
341 // Find gray-Tones spatial dependence matrix
342 for(iy=y-m_Radius; iy<=y+m_Radius; iy++)
343 {
344 for(ix=x-m_Radius; ix<=x+m_Radius; ix++)
345 {
346 int j = 0; while(Tones[j] != Get_Value(ix, iy)) j++;
347
348 if( ix + d <= x + m_Radius )
349 {
350 int i = 0; while(Tones[i] != Get_Value(ix + d, iy )) i++;
351
352 P[0][j][i]++;
353 P[0][i][j]++;
354 }
355
356 if( iy + d <= y + m_Radius )
357 {
358 int i = 0; while(Tones[i] != Get_Value(ix , iy + d)) i++;
359
360 P[2][j][i]++;
361 P[2][i][j]++;
362 }
363
364 if( iy + d <= y + m_Radius && ix - d >= x - m_Radius )
365 {
366 int i = 0; while(Tones[i] != Get_Value(ix - d, iy + d)) i++;
367
368 P[1][j][i]++;
369 P[1][i][j]++;
370 }
371
372 if( iy + d <= y + m_Radius && ix + d <= x + m_Radius )
373 {
374 int i = 0; while(Tones[i] != Get_Value(ix + d, iy + d)) i++;
375
376 P[3][j][i]++;
377 P[3][i][j]++;
378 }
379 }
380 }
381
382 //-----------------------------------------------------
383 // Normalize gray-Tones spatial dependence matrix
384 int n = 1 + 2 * m_Radius;
385
386 P[0] *= 1. / (2. * (n ) * (n - 1));
387 P[1] *= 1. / (2. * (n - 1) * (n - 1));
388 P[2] *= 1. / (2. * (n - 1) * (n ));
389 P[3] *= 1. / (2. * (n - 1) * (n - 1));
390
391 //-----------------------------------------------------
392 return( true );
393 }
394
395
396 ///////////////////////////////////////////////////////////
397 // //
398 // //
399 // //
400 ///////////////////////////////////////////////////////////
401
402 //---------------------------------------------------------
403 // in the following are those parts of the original grass implementation
404 // (r.texture/h_measure.c) responsible for the calculation of the 'measures' from the
405 // occurrence/co-occurrence matrices:
406 //
407 // MODULE: r.texture
408 // AUTHOR(S): Carmine Basco - basco@unisannio.it
409 // with hints from:
410 // prof. Giulio Antoniol - antoniol@ieee.org
411 // prof. Michele Ceccarelli - ceccarelli@unisannio.it
412 //
413 //---------------------------------------------------------
414
415 #define EPSILON 0.000000001
416
f1_asm(const double ** P,int Ng)417 double f1_asm(const double **P, int Ng)
418 {
419 double sum = 0.0;
420
421 for(int i=0; i<Ng; i++)
422 for(int j=0; j<Ng; j++)
423 sum += P[i][j] * P[i][j];
424
425 return sum;
426 }
427
f2_contrast(const double ** P,int Ng)428 double f2_contrast(const double **P, int Ng)
429 {
430 int i, j, n;
431 double sum, bigsum = 0;
432
433 for (n = 0; n < Ng; n++) {
434 sum = 0;
435 for (i = 0; i < Ng; i++) {
436 for (j = 0; j < Ng; j++) {
437 if ((i - j) == n || (j - i) == n) {
438 sum += P[i][j];
439 }
440 }
441 }
442 bigsum += n * n * sum;
443 }
444 return bigsum;
445 }
446
f3_corr(const double ** P,int Ng,const double * px)447 double f3_corr(const double **P, int Ng, const double *px)
448 {
449 int i, j;
450 double sum_sqrx = 0, sum_sqry = 0, tmp = 0;
451 double meanx = 0, meany = 0, stddevx, stddevy;
452
453
454 /* Now calculate the means and standard deviations of px and py */
455
456 /*- fix supplied by J. Michael Christensen, 21 Jun 1991 */
457
458 /*- further modified by James Darrell McCauley, 16 Aug 1991
459 * after realizing that meanx=meany and stddevx=stddevy
460 */
461 for (i = 0; i < Ng; i++) {
462 meanx += px[i] * i;
463 sum_sqrx += px[i] * i * i;
464
465 for (j = 0; j < Ng; j++)
466 tmp += i * j * P[i][j];
467 }
468 meany = meanx;
469 sum_sqry = sum_sqrx;
470 stddevx = sqrt(sum_sqrx - (meanx * meanx));
471 stddevy = stddevx;
472
473 return (tmp - meanx * meany) / (stddevx * stddevy);
474 }
475
476 /* Sum of Squares: Variance */
f4_var(const double ** P,int Ng)477 double f4_var(const double **P, int Ng)
478 {
479 int i, j;
480 double mean = 0, var = 0;
481
482 /*- Corrected by James Darrell McCauley, 16 Aug 1991
483 * calculates the mean intensity level instead of the mean of
484 * cooccurrence matrix elements
485 */
486 for (i = 0; i < Ng; i++)
487 for (j = 0; j < Ng; j++)
488 mean += i * P[i][j];
489
490 for (i = 0; i < Ng; i++)
491 for (j = 0; j < Ng; j++)
492 var += (i + 1 - mean) * (i + 1 - mean) * P[i][j];
493
494 return var;
495 }
496
497 /* Inverse Difference Moment */
f5_idm(const double ** P,int Ng)498 double f5_idm(const double **P, int Ng)
499 {
500 int i, j;
501 double idm = 0;
502
503 for (i = 0; i < Ng; i++)
504 for (j = 0; j < Ng; j++)
505 idm += P[i][j] / (1 + (i - j) * (i - j));
506
507 return idm;
508 }
509
510 /* Sum Average */
f6_savg(const double ** P,int Ng,const double * Pxpys)511 double f6_savg(const double **P, int Ng, const double *Pxpys)
512 {
513 int i;
514 double savg = 0;
515
516 for (i = 0; i < 2 * Ng - 1; i++)
517 savg += (i + 2) * Pxpys[i];
518
519 return savg;
520 }
521
522 /* Sum Variance */
f7_svar(const double ** P,int Ng,double S,const double * Pxpys)523 double f7_svar(const double **P, int Ng, double S, const double *Pxpys)
524 {
525 int i;
526 double var = 0;
527
528 for (i = 0; i < 2 * Ng - 1; i++)
529 var += (i + 2 - S) * (i + 2 - S) * Pxpys[i];
530
531 return var;
532 }
533
534 /* Sum Entropy */
f8_sentropy(const double ** P,int Ng,const double * Pxpys)535 double f8_sentropy(const double **P, int Ng, const double *Pxpys)
536 {
537 int i;
538 double sentr = 0;
539
540 for (i = 0; i < 2 * Ng - 1; i++)
541 sentr -= Pxpys[i] * log10(Pxpys[i] + EPSILON);
542
543 return sentr;
544 }
545
546 /* Entropy */
f9_entropy(const double ** P,int Ng)547 double f9_entropy(const double **P, int Ng)
548 {
549 int i, j;
550 double entropy = 0;
551
552 for (i = 0; i < Ng; i++) {
553 for (j = 0; j < Ng; j++) {
554 entropy += P[i][j] * log10(P[i][j] + EPSILON);
555 }
556 }
557
558 return -entropy;
559 }
560
561 /* Difference Variance */
f10_dvar(const double ** P,int Ng,const double * Pxpyd)562 double f10_dvar(const double **P, int Ng, const double *Pxpyd)
563 {
564 int i, tmp;
565 double sum = 0, sum_sqr = 0, var = 0;
566
567 /* Now calculate the variance of Pxpy (Px-y) */
568 for (i = 0; i < Ng; i++) {
569 sum += Pxpyd[i];
570 sum_sqr += Pxpyd[i] * Pxpyd[i];
571 }
572 tmp = Ng * Ng;
573 var = ((tmp * sum_sqr) - (sum * sum)) / (tmp * tmp);
574
575 return var;
576 }
577
578 /* Difference Entropy */
f11_dentropy(const double ** P,int Ng,const double * Pxpyd)579 double f11_dentropy(const double **P, int Ng, const double *Pxpyd)
580 {
581 int i;
582 double sum = 0;
583
584 for (i = 0; i < Ng; i++)
585 sum += Pxpyd[i] * log10(Pxpyd[i] + EPSILON);
586
587 return -sum;
588 }
589
590 /* Information Measures of Correlation */
f12_icorr(const double ** P,int Ng,const double * px,const double * py)591 double f12_icorr(const double **P, int Ng, const double *px, const double *py)
592 {
593 int i, j;
594 double hx = 0, hy = 0, hxy = 0, hxy1 = 0;
595
596 for (i = 0; i < Ng; i++)
597 for (j = 0; j < Ng; j++) {
598 hxy1 -= P[i][j] * log10(px[i] * py[j] + EPSILON);
599 hxy -= P[i][j] * log10(P[i][j] + EPSILON);
600 }
601
602 /* Calculate entropies of px and py - is this right? */
603 for (i = 0; i < Ng; i++) {
604 hx -= px[i] * log10(px[i] + EPSILON);
605 hy -= py[i] * log10(py[i] + EPSILON);
606 }
607
608 /* fprintf(stderr,"hxy1=%f\thxy=%f\thx=%f\thy=%f\n",hxy1,hxy,hx,hy); */
609 return ((hxy - hxy1) / (hx > hy ? hx : hy));
610 }
611
612 /* Information Measures of Correlation */
f13_icorr(const double ** P,int Ng,const double * px,const double * py)613 double f13_icorr(const double **P, int Ng, const double *px, const double *py)
614 {
615 int i, j;
616 double hxy = 0, hxy2 = 0;
617
618 for (i = 0; i < Ng; i++) {
619 for (j = 0; j < Ng; j++) {
620 hxy2 -= px[i] * py[j] * log10(px[i] * py[j] + EPSILON);
621 hxy -= P[i][j] * log10(P[i][j] + EPSILON);
622 }
623 }
624
625 /* fprintf(stderr,"hx=%f\thxy2=%f\n",hx,hxy2); */
626 return (sqrt(fabs(1 - exp(-2.0 * (hxy2 - hxy)))));
627 }
628
629
630 ///////////////////////////////////////////////////////////
631 // //
632 ///////////////////////////////////////////////////////////
633
634 //---------------------------------------------------------
Get_Features(CSG_Vector & Features,const CSG_Matrix & P)635 bool CTextural_Features::Get_Features(CSG_Vector &Features, const CSG_Matrix &P)
636 {
637 //-----------------------------------------------------
638 int nTones = P.Get_NCols();
639
640 CSG_Vector px (nTones );
641 CSG_Vector py (nTones );
642 CSG_Vector Pxpys(nTones * 2);
643 CSG_Vector Pxpyd(nTones * 2);
644
645 for(int i=0; i<nTones; i++)
646 {
647 for(int j=0; j<nTones; j++)
648 {
649 px[i] += P[i][j];
650 py[j] += P[i][j];
651 Pxpys[i + j] += P[i][j];
652 Pxpyd[abs(i - j)] += P[i][j];
653 }
654 }
655
656 //-----------------------------------------------------
657 double Sum_Entropy;
658
659 Features[ASM ] += f1_asm (P, nTones);
660 Features[CONTRAST ] += f2_contrast (P, nTones);
661 Features[CORRELATION ] += f3_corr (P, nTones, px);
662 Features[VARIANCE ] += f4_var (P, nTones);
663 Features[IDM ] += f5_idm (P, nTones);
664 Features[SUM_AVERAGE ] += f6_savg (P, nTones, Pxpys);
665 Features[SUM_ENTROPY ] += Sum_Entropy = f8_sentropy (P, nTones, Pxpys);
666 Features[SUM_VARIANCE] += f7_svar (P, nTones, Sum_Entropy, Pxpys);
667 Features[ENTROPY ] += f9_entropy (P, nTones);
668 Features[DIF_VARIANCE] += f10_dvar (P, nTones, Pxpyd);
669 Features[DIF_ENTROPY ] += f11_dentropy(P, nTones, Pxpyd);
670 Features[MOC_1 ] += f12_icorr (P, nTones, px, py);
671 Features[MOC_2 ] += f13_icorr (P, nTones, px, py);
672
673 //-----------------------------------------------------
674 return( true );
675 }
676
677
678 ///////////////////////////////////////////////////////////
679 // //
680 // //
681 // //
682 ///////////////////////////////////////////////////////////
683
684 //---------------------------------------------------------
685