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);
}