c_sep.c

Go to the documentation of this file.

/*!
   \file cluster/c_sep.c
 
   \brief Cluster library - Separation
 
   (C) 2001-2009 by the GRASS Development Team
 
   This program is free software under the GNU General Public License
   (>=v2). Read the file COPYING that comes with GRASS for details.
 
   \author Original author CERL
 */
 
#include <math.h>
#include <grass/cluster.h>
 
#define FAR ((double)-1.0)
 
/*!
   \brief ?
 
   \param C pointer to Cluster structure
   \param class1 1st class
   \param class2 2nd class
 */
double I_cluster_separation(struct Cluster *C, int class1, int class2)
{
    int band;
    double q;
    double d;
    double var;
    double a1, a2;
    double n1, n2;
    double m1, m2;
    double s1, s2;
 
    if (C->count[class1] < 2)
        return FAR;
    if (C->count[class2] < 2)
        return FAR;
    n1 = (double)C->count[class1];
    n2 = (double)C->count[class2];
 
    d = 0.0;
    a1 = a2 = 0.0;
    for (band = 0; band < C->nbands; band++) {
        s1 = C->sum[band][class1];
        s2 = C->sum[band][class2];
        m1 = s1 / n1;
        m2 = s2 / n2;
        q = m1 - m2;
        q = q * q;
        d += q;
 
        var = C->sum2[band][class1] - (s1 * m1);
        var /= n1 - 1;
        if (var)
            a1 += q / var;
 
        var = C->sum2[band][class2] - (s2 * m2);
        var /= n2 - 1;
        if (var)
            a2 += q / var;
    }
    if (d == 0.0)
        return d;
 
    if (a1 < 0 || a2 < 0)
        return FAR;
    if (a1)
        a1 = sqrt(6 * d / a1);
    if (a2)
        a2 = sqrt(6 * d / a2);
    q = a1 + a2;
    if (q == 0.0)
        return FAR;
 
    return (sqrt(d) / q);
}