---

mindist.c

Go to the documentation of this file.

#include        "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1992 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function performs a minimum distance classification on an image.    */
/*                                                                            */
/*   The discriminant function for classification of a pixel to a class is:   */
/*                                                                            */
/*   City-block distance classifier:                                          */
/*      DC(X) = |X - M|                                                       */
/*                                                                            */
/*   Euclidean distance classifier:                                           */
/*      DE(X) = (X - M)**t * (X - M)                                          */
/*                                                                            */
/*   Mahalanobis distance classifier:                                         */
/*      DM(X) = (X - M)**t * C**-1 * (X - M)                                  */
/*                                                                            */
/*   where X is the feature vector, M is the class mean vector, and C**-1     */
/*   is the inverse of the covariance matrix (equal for all classes).         */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
MINDIST (IMAGE * in,            /*  I  Pointer to the input image.                             */
         short option,          /*  I  Switch, indicating the type of classifier:              */
         /*     CITYBLOCK     -   DC(X)                                 */
         /*     EUCLIDEAN     -   DE(X)                                 */
         /*     MAHALANOBIS   -   DM(X)                                 */
         short nlev,            /*  I  Number of classes                                       */
         /*     (and number of graylevels in the output image).         */
         double *mean,          /*  I  Pointer to the matrix[nlev][nbnd] of mean graylevels.   */
         double *cov,           /*  I  Pointer to the covariance matrix[nbnd][nbnd].           */
         IMAGE ** out           /*  O  Address of a pointer to the output image.               */
/*--------------------------------------------------------------------------- */
  )
{
  register short i, j, k, m, n, l1, l2, ll, nl, np;
  char msg[SLEN];
  PIXEL *dist = 0, mindist;
  double *inv = 0;

  if (TIMES)
    TIMING (T_MINDIST);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "MINDIST");
      MESSAGE ('I', "");
      sprintf (msg, " Input image:          %s", in->text);
      MESSAGE ('I', msg);
      if (option == CITYBLOCK)
        {
          MESSAGE ('I', " Classifier:           City-block");
        }
      else if (option == EUCLIDEAN)
        {
          MESSAGE ('I', " Classifier:           Euclidean");
        }
      else if (option == MAHALANOBIS)
        {
          MESSAGE ('I', " Classifier:           Mahalanobis");
        }
      sprintf (msg, " Number of classes:    %d", nlev);
      MESSAGE ('I', msg);
      sprintf (msg, " Mean graylevels:");
      MESSAGE ('I', msg);
      for (n = j = 0; j < nlev; j++)
        {
          sprintf (msg, "  Class%3d:          ", j + 1);
          for (k = 0, m = strlen (msg); k < in->nbnd; k++, m = strlen (msg))
            {
              sprintf (msg + m, "%12.4e", mean[n++]);
            }
          MESSAGE ('I', msg);
        }
      MESSAGE ('I', " Covariance matrix:");
      for (n = j = 0; j < in->nbnd; j++)
        {
          sprintf (msg, "  Band%4d:          ", j + 1);
          for (k = 0, m = strlen (msg); k < in->nbnd; k++, m = strlen (msg))
            {
              sprintf (msg + m, "%12.4e", cov[n++]);
            }
          MESSAGE ('I', msg);
        }
      MESSAGE ('I', " ...............");
    }

  /* check input  */
  if (!CHECKIMG (in))
    {
      MESSAGE ('E', "Can't identify image.");
      goto the_end;
    }
  else if (option != CITYBLOCK && option != EUCLIDEAN
           && option != MAHALANOBIS)
    {
      MESSAGE ('E',
               " Option must be one of CITYBLOCK, EUCLIDEAN, or MAHALANOBIS.");
      goto the_end;
    }
  else if (nlev <= 0)
    {
      MESSAGE ('E', "Number of classes must be greater than zero.");
      goto the_end;
    }

  /* create image of appropriate size */
  if (!CHECKIMG (*out))
    GETMEM (1, in->nlin, in->npix, out);
  if (!*out)
    goto the_end;

  /* allocate space for  distance array */
  if (!(dist = (PIXEL *) malloc (nlev * sizeof (PIXEL))))
    {
      MESSAGE ('E', "Memory request failed.");
      goto the_end;
    }

  /* City-block distance classifier */
  if (option == CITYBLOCK)
    {
      for (nl = 0; nl < in->nlin; nl++)
        {
          for (np = 0; np < in->npix; np++)
            {
              for (k = 0; k < nlev; k++)
                {
                  dist[k] = 0.0;
                  for (i = 0; i < in->nbnd; i++)
                    {
                      l1 = k * in->nbnd + i;
                      dist[k] += fabs (in->data[i][nl][np] - mean[l1]);
                    }
                }
              mindist = dist[0];
              for (i = 0; i < nlev; i++)
                {
                  if (dist[i] <= mindist)
                    {
                      (*out)->data[0][nl][np] = (PIXEL) (i + 1);
                    }
                }
            }
        }

      /* Euclidean distance classifier */
    }
  else if (option == EUCLIDEAN)
    {
      for (nl = 0; nl < in->nlin; nl++)
        {
          for (np = 0; np < in->npix; np++)
            {
              for (k = 0; k < nlev; k++)
                {
                  dist[k] = 0.0;
                  for (i = 0; i < in->nbnd; i++)
                    {
                      l1 = k * in->nbnd + i;
                      dist[k] +=
                        (in->data[i][nl][np] -
                         mean[l1]) * (in->data[i][nl][np] - mean[l1]);
                    }
                }
              mindist = dist[0];
              for (i = 0; i < nlev; i++)
                {
                  if (dist[i] <= mindist)
                    {
                      (*out)->data[0][nl][np] = (PIXEL) (i + 1);
                    }
                }
            }
        }

      /* Mahalanobis distance classifier */
    }
  else if (option == MAHALANOBIS)
    {

      /* allocate space for inverse of covariance matrix and distance array */
      if (!(inv = (double *) malloc (in->nbnd * in->nbnd * sizeof (double))))
        {
          MESSAGE ('E', "Memory request failed.");
          goto the_end;
        }
      if (!(dist = (PIXEL *) malloc (nlev * sizeof (PIXEL))))
        {
          MESSAGE ('E', "Memory request failed.");
          goto the_end;
        }
      MATRIX_INVERT (cov, in->nbnd, inv);

      for (nl = 0; nl < in->nlin; nl++)
        {
          for (np = 0; np < in->npix; np++)
            {
              for (k = 0; k < nlev; k++)
                {
                  dist[k] = 0.0;
                  for (i = 0; i < in->nbnd; i++)
                    {
                      for (j = 0; j < in->nbnd; j++)
                        {
                          l1 = k * in->nbnd + i;
                          l2 = k * in->nbnd + j;
                          ll = i * in->nbnd + j;
                          dist[k] +=
                            (in->data[i][nl][np] -
                             mean[l1]) * inv[ll] * (in->data[j][nl][np] -
                                                    mean[l2]);
                        }
                    }
                }
              mindist = dist[0];
              for (i = 0; i < nlev; i++)
                {
                  if (dist[i] <= mindist)
                    {
                      (*out)->data[0][nl][np] = (PIXEL) (i + 1);
                    }
                }
            }
        }
    }

the_end:
  if (inv)
    free (inv);
  if (dist)
    free (dist);
  if (TIMES)
    TIMING (T_EXIT);
}

---