---

cluster.c

Go to the documentation of this file.

#include        "sadie.h"

typedef struct _MERGE
{
  short cls1;
  short cls2;
  PIXEL dist;
  struct _MERGE *next;
} MERGE;

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1990 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function determines the intrinsic groupings of pixels in feature    */
/*   space, using an iterative K-means clustering algorithm. Clusters with    */
/*   too few pixels are discarded, and clusters that are too close are merged */
/*   in the iteration process. After the clusters have been determined, the   */
/*   entire input image is classified in a final pass.                        */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Background Information-----------------------------------------------------*/
/*                                                                            */
/*   Duda, R.D., and Hart, P.E.:                                              */
/*   "Pattern Classification and Scene Analysis,"                             */
/*   J. Wiley and Sons, New York, 1973.                                       */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/* Interface Information------------------------------------------------------*/
void
CLUSTER (IMAGE * in,            /*  I   Pointer to the input image.                           */
         short niter,           /*  I   Number of iterations.                                 */
         short nclass,          /*  I   Number of seed classes.                               */
         short csize,           /*  |   Minimum number of pixels in each class.               */
         PIXEL mthr,            /*  I   Threshold for merging clusters.                       */
         short iinc,            /*  I   Band increment.                                       */
         short jinc,            /*  I   Line increment.                                       */
         short kinc,            /*  I   Pixels/line increment.                                */
         PIXEL othr,            /*  I   Outlier threshold.                                    */
         IMAGE ** out           /*  O   Address of a pointer to the output image.             */
/*----------------------------------------------------------------------------*/
  )
{
  register short i, j, k, m, n, class, nsize, *mrgd = 0;
  char msg[SLEN];
  long *cnt = 0;
  double pinc, psum;
  PIXEL *mean = 0, *nwmn = 0, *dev =
    0, gmin, gmax, ginc, gsum, gmig, gdiff, gval;
  MERGE *head = 0, *tail, *ptr, *newClass;

  if (TIMES)
    TIMING (T_CLUSTER);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "CLUSTER");
      MESSAGE ('I', "");
      sprintf (msg, " Input image:                   %s", in->text);
      MESSAGE ('I', msg);
      sprintf (msg, " Number of iterations:          %d", niter);
      MESSAGE ('I', msg);
      sprintf (msg, " Number of seed classes:        %d", nclass);
      MESSAGE ('I', msg);
      sprintf (msg, " Minimum class size:            %d", csize);
      MESSAGE ('I', msg);
      sprintf (msg, " Class merging threshold:     %12.4e", mthr);
      MESSAGE ('I', msg);
      sprintf (msg, " Subsample increment: bands:    %d", iinc);
      MESSAGE ('I', msg);
      sprintf (msg, "                      lines:    %d", jinc);
      MESSAGE ('I', msg);
      sprintf (msg, "                      pixels:   %d", kinc);
      MESSAGE ('I', msg);
      sprintf (msg, " Outlier threshold:           %12.4e", othr);
      MESSAGE ('I', msg);
      MESSAGE ('I', " ...............");
    }

  /* check input */
  if (!CHECKIMG (in))
    {
      MESSAGE ('E', " Can't identify image.");
      goto the_end;
    }
  else if (iinc <= 0)
    {
      MESSAGE ('E', " Band subsample increment must be greater than zero.");
      goto the_end;
    }
  else if (jinc <= 0)
    {
      MESSAGE ('E', " Line subsample increment must be greater than zero.");
      goto the_end;
    }
  else if (kinc <= 0)
    {
      MESSAGE ('E', " Pixel subsample increment must be greater than zero.");
      goto the_end;
    }
  else if (niter <= 0)
    {
      MESSAGE ('E', " Number of iterations must be greater than zero.");
      goto the_end;
    }
  else if (nclass <= 0)
    {
      MESSAGE ('E', " Number of classes must be greater than zero.");
      goto the_end;
    }
  else if (csize <= 0)
    {
      MESSAGE ('E', " Minimum class size must be greater than zero.");
      goto the_end;
    }
  else if (mthr < 0.0)
    {
      MESSAGE ('E',
               " Distance threshold must be greater than or equal to zero.");
      goto the_end;
    }
  else if (othr < 0.0)
    {
      MESSAGE ('E',
               " Distance threshold must be greater than or equal to 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 */
  nsize = nclass + 1;
  if (!(mrgd = (short *) malloc (nsize * sizeof (short))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(cnt = (long *) malloc (nsize * sizeof (long))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(mean = (PIXEL *) malloc (nsize * in->nbnd * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(nwmn = (PIXEL *) malloc (nsize * in->nbnd * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(dev = (PIXEL *) malloc (nsize * in->nbnd * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!
      (head =
       (MERGE *) malloc ((((nclass - 1) * nclass) / 2 + 2) * sizeof (MERGE))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }

  /* initialize progress indicator */
  if (LINES && !PROGRESS (psum = 0.0))
    goto the_end;
  pinc =
    0.5 / (double) (niter * (in->nlin / jinc + (in->nlin % jinc ? 1 : 0)) +
                    in->nlin);

  /* initialize */
  for (i = 0; i < in->nbnd; i++)
    {
      gmax = gmin = in->data[i][0][0];
      for (j = 0; j < in->nlin; j++)
        {
          for (k = 0; k < in->npix; k++)
            {
              gmin = min (gmin, in->data[i][j][k]);
              gmax = max (gmax, in->data[i][j][k]);
            }
        }
      ginc = (gmax - gmin) / (PIXEL) nclass;
      for (m = 1; m < nsize; m++)
        {
          mean[i * nsize + m] =
            gmin + (PIXEL) (((double) m - 0.5) * (double) ginc);
        }
    }
  tail = head + 1;
  tail->cls1 = 0;
  tail->cls2 = 0;
  tail->dist = mthr * (PIXEL) in->nbnd;
  tail->next = 0;

  /* iterate */
  for (m = 0; m < niter; m++)
    {

      /* initialize */
      for (n = 0; n < nsize; n++)
        {
          cnt[n] = 0L;
          mrgd[n] = TRUE;
          for (i = 0; i < in->nbnd; i++)
            {
              nwmn[i * nsize + n] = (PIXEL) 0;
              dev[i * nsize + n] = (PIXEL) 0;
            }
        }

      /* cluster */
      for (class = 1, j = 0; j < in->nlin; j += jinc)
        {
          for (k = 0; k < in->npix; k += kinc)
            {
              gmin = (PIXEL) 0;
              for (i = 0; i < in->nbnd; i += iinc)
                {
                  gmin +=
                    (PIXEL)
                    fabs ((double)
                          (in->data[i][j][k] - mean[i * nsize + class]));
                }
              for (n = 1; n <= nclass; n++)
                {
                  if (n != class)
                    {
                      gsum = (PIXEL) 0;
                      for (i = 0; i < in->nbnd && gmin > gsum; i += iinc)
                        {
                          gsum +=
                            (PIXEL)
                            fabs ((double)
                                  (in->data[i][j][k] - mean[i * nsize + n]));
                        }
                      if (gmin > gsum)
                        {
                          gmin = gsum;
                          class = n;
                        }
                    }
                }
              (*out)->data[0][j][k] = (PIXEL) class;
              for (i = 0; i < in->nbnd; i += iinc)
                {
                  nwmn[i * nsize + class] += in->data[i][j][k];
                }
              cnt[class] += 1L;
            }
          if (LINES && !PROGRESS (psum += pinc))
            goto the_end;
        }

      /* compute new mean vectors and average migration */
      gmig = (PIXEL) 0;
      for (i = 0; i < in->nbnd; i += iinc)
        {
          for (n = 1; n <= nclass; n++)
            {
              if (cnt[n] > 0L)
                {
                  nwmn[i * nsize + n] /= (PIXEL) cnt[n];
                  gmig +=
                    (PIXEL)
                    fabs ((double)
                          (nwmn[i * nsize + n] - mean[i * nsize + n]));
                  mean[i * nsize + n] = nwmn[i * nsize + n];
                }
            }
        }
      gmig /= (double) (in->nbnd * nclass);

      /* compute new standard deviation vectors */
      for (j = 0; j < in->nlin; j += jinc)
        {
          for (k = 0; k < in->npix; k += kinc)
            {
              for (i = 0; i < in->nbnd; i += iinc)
                {
                  class = (short) (*out)->data[0][j][k];
                  gdiff = in->data[i][j][k] - mean[i * nsize + class];
                  dev[i * nsize + class] += gdiff * gdiff;
                }
            }
          if (LINES && !PROGRESS (psum += pinc))
            goto the_end;
        }
      for (n = 1; n <= nclass; n++)
        {
          for (i = 0; i < in->nbnd; i += iinc)
            {
              if (cnt[n] > 0L)
                {
                  dev[i * nsize + n] =
                    (PIXEL) sqrt ((double) dev[i * nsize + n] /
                                  (double) cnt[n]);
                }
            }
        }

      /* eliminate classes that are too small */
      for (n = 1; n <= nclass; n++)
        {
          if (cnt[n] < csize)
            {
              for (j = n; j < nclass; j++)
                {
                  for (i = 0; i < in->nbnd; i += iinc)
                    {
                      mean[i * nsize + j] = mean[i * nsize + j + 1];
                      dev[i * nsize + j] = dev[i * nsize + j + 1];
                    }
                  cnt[j] = cnt[j + 1];
                }
              nclass -= 1;
              n -= 1;
            }
        }

      /* find and merge classes that are too close together */
      newClass = head + 2;
      head->next = tail;
      for (n = 1; n <= nclass - 1; n++)
        {
          for (k = n + 1; k <= nclass; k++)
            {
              gsum = (PIXEL) 0;
              for (i = 0; i < in->nbnd; i += iinc)
                {
                  gval = dev[i * nsize + n] + dev[i * nsize + k];
                  if (gval != (PIXEL) 0
                      && (gval =
                          (PIXEL)
                          fabs ((double)
                                (mean[i * nsize + n] -
                                 mean[i * nsize + k])) / gval) < mthr)
                    {
                      gsum += gval;
                      if (i + iinc >= in->nbnd)
                        {
                          for (ptr = head; ptr->next->dist <= gsum;
                               ptr = ptr->next);
                          newClass->cls1 = n;
                          newClass->cls2 = k;
                          newClass->dist = gsum;
                          newClass->next = ptr->next;
                          ptr->next = newClass++;
                        }
                    }
                  else
                    {
                      break;
                    }
                }
            }
        }
      for (ptr = head->next; ptr != tail; ptr = ptr->next)
        {
          if (mrgd[ptr->cls1] == TRUE && mrgd[ptr->cls2] == TRUE)
            {
              for (i = 0; i < in->nbnd; i += iinc)
                {
                  gval =
                    mean[i * nsize + ptr->cls1] * cnt[ptr->cls1] +
                    mean[i * nsize + ptr->cls2] * cnt[ptr->cls2];
                  mean[i * nsize + ptr->cls1] =
                    gval / (PIXEL) (cnt[ptr->cls1] + cnt[ptr->cls2]);
                }
              for (j = ptr->cls2; j <= nclass - 1; j++)
                {
                  for (i = 0; i < in->nbnd; i++)
                    {
                      mean[i * nsize + j] = mean[i * nsize + j + 1];
                    }
                }
              mrgd[ptr->cls1] = FALSE;
              mrgd[ptr->cls2] = FALSE;
              nclass -= 1;
            }
        }

      /* output iteration statistics */
      if (m == 0)
        {
          MESSAGE ('I', "");
          MESSAGE ('I', " Mean vector migration averaged over all classes:");
          MESSAGE ('I',
                   " Iteration     Migration     Number of classes after merging");
        }
      sprintf (msg, "%6d%18.4e%20d", m + 1, (double) gmig, nclass);
      MESSAGE ('I', msg);
    }

  /* initialize */
  for (n = 0; n < nsize; n++)
    {
      cnt[n] = 0L;
      for (i = 0; i < in->nbnd; i++)
        {
          nwmn[i * nsize + n] = (PIXEL) 0;
        }
    }

  /* classify */
  for (class = 1, j = 0; j < in->nlin; j++)
    {
      for (k = 0; k < in->npix; k++)
        {
          class = max (1, class);
          gmin = (PIXEL) 0;
          for (i = 0; i < in->nbnd; i++)
            {
              gmin +=
                (PIXEL)
                fabs ((double) (in->data[i][j][k] - mean[i * nsize + class]));
            }
          for (n = 1; n <= nclass; n++)
            {
              if (n != class)
                {
                  gsum = (PIXEL) 0;
                  for (i = 0; i < in->nbnd && gmin > gsum; i++)
                    {
                      gsum +=
                        (PIXEL)
                        fabs ((double)
                              (in->data[i][j][k] - mean[i * nsize + n]));
                    }
                  if (gmin > gsum)
                    {
                      gmin = gsum;
                      class = n;
                    }
                }
            }
          for (i = 0; i < in->nbnd; i++)
            {
              if ((PIXEL)
                  fabs ((double)
                        (in->data[i][j][k] - mean[i * nsize + class])) >
                  othr * dev[i * nsize + class])
                {
                  class = 0;
                  break;
                }
            }
          (*out)->data[0][j][k] = class;
          if (class != 0)
            {
              for (i = 0; i < in->nbnd; i++)
                {
                  nwmn[i * nsize + class] += in->data[i][j][k];
                }
            }
          cnt[class] += 1L;
        }
      if (LINES && !PROGRESS (psum += pinc))
        goto the_end;
    }

  /* compute new mean vectors and average migration */
  gmig = (PIXEL) 0;
  for (i = 0; i < in->nbnd; i++)
    {
      for (n = 1; n <= nclass; n++)
        {
          if (cnt[n] > 0L)
            {
              nwmn[i * nsize + n] /= (PIXEL) cnt[n];
              gmig +=
                (PIXEL)
                fabs ((double) (nwmn[i * nsize + n] - mean[i * nsize + n]));
              mean[i * nsize + n] = nwmn[i * nsize + n];
            }
        }
    }
  gmig /= (double) (in->nbnd * nclass);

  /* compute new standard deviation vectors */
  for (i = 0; i < nsize * in->nbnd; i++)
    {
      dev[i] = (PIXEL) 0;
    }
  for (j = 0; j < in->nlin; j++)
    {
      for (k = 0; k < in->npix; k++)
        {
          for (i = 0; i < in->nbnd; i++)
            {
              class = (short) (*out)->data[0][j][k];
              if (class != 0)
                {
                  gdiff = in->data[i][j][k] - mean[i * nsize + class];
                  dev[i * nsize + class] += gdiff * gdiff;
                }
            }
        }
      if (LINES && !PROGRESS (psum += pinc))
        goto the_end;
    }
  for (n = 1; n <= nclass; n++)
    {
      for (i = 0; i < in->nbnd; i++)
        {
          if (cnt[n] > 0L)
            {
              dev[i * nsize + n] =
                (PIXEL) sqrt ((double) dev[i * nsize + n] / (double) cnt[n]);
            }
        }
    }

  /* output iteration statistics */
  sprintf (msg, "     *%18.4e%20d", (double) gmig, nclass);
  MESSAGE ('I', msg);

  /* output classification statistics */
  MESSAGE ('I', "");
  MESSAGE ('I', " Mean and standard deviation vectors of each class:");
  sprintf (msg, " Class  ");
  for (i = 0; i < in->nbnd; i++)
    {
      sprintf (msg + strlen (msg), "           Band %-11d", i + 1);
    }
  sprintf (msg + strlen (msg), "# Pixels     %%");
  MESSAGE ('I', msg);
  sprintf (msg, "%6d", 0);
  for (i = 0; i < in->nbnd; i++)
    {
      sprintf (msg + strlen (msg), "              N/A          ", "");
    }
  sprintf (msg + strlen (msg), "%10ld%7.1f   (Threshold class)", cnt[0],
           100.0 * (double) cnt[0] / ((double) in->nlin * (double) in->npix));
  MESSAGE ('I', msg);
  for (n = 1; n <= nclass; n++)
    {
      sprintf (msg, "%6d", n);
      for (i = 0; i < in->nbnd; i++)
        {
          sprintf (msg + strlen (msg), "%13.4e +/-%10.4e",
                   (double) mean[i * nsize + n], (double) dev[i * nsize + n]);
        }
      sprintf (msg + strlen (msg), "%10ld%7.1f", cnt[n],
               100.0 * (double) cnt[n] / ((double) in->nlin *
                                          (double) in->npix));
      MESSAGE ('I', msg);
    }

the_end:
  if (mrgd)
    free (mrgd);
  if (cnt)
    free (cnt);
  if (mean)
    free (mean);
  if (nwmn)
    free (nwmn);
  if (dev)
    free (dev);
  if (head)
    free (head);
  if (TIMES)
    TIMING (T_EXIT);
}

---