GRASS Programmer's Manual
6.4.2(2012)
|
00001 #include <math.h> 00002 #include <grass/cluster.h> 00003 00004 #define FAR ((double) -1.0) 00005 00006 double I_cluster_separation(struct Cluster *C, int class1, int class2) 00007 { 00008 int band; 00009 double q; 00010 double d; 00011 double var; 00012 double a1, a2; 00013 double n1, n2; 00014 double m1, m2; 00015 double s1, s2; 00016 00017 if (C->count[class1] < 2) 00018 return FAR; 00019 if (C->count[class2] < 2) 00020 return FAR; 00021 n1 = (double)C->count[class1]; 00022 n2 = (double)C->count[class2]; 00023 00024 d = 0.0; 00025 a1 = a2 = 0.0; 00026 for (band = 0; band < C->nbands; band++) { 00027 s1 = C->sum[band][class1]; 00028 s2 = C->sum[band][class2]; 00029 m1 = s1 / n1; 00030 m2 = s2 / n2; 00031 q = m1 - m2; 00032 q = q * q; 00033 d += q; 00034 00035 00036 var = C->sum2[band][class1] - (s1 * m1); 00037 var /= n1 - 1; 00038 if (var) 00039 a1 += q / var; 00040 00041 var = C->sum2[band][class2] - (s2 * m2); 00042 var /= n2 - 1; 00043 if (var) 00044 a2 += q / var; 00045 } 00046 if (d == 0.0) 00047 return d; 00048 00049 if (a1 < 0 || a2 < 0) 00050 return FAR; 00051 if (a1) 00052 a1 = sqrt(6 * d / a1); 00053 if (a2) 00054 a2 = sqrt(6 * d / a2); 00055 q = a1 + a2; 00056 if (q == 0.0) 00057 return FAR; 00058 00059 return (sqrt(d) / q); 00060 }