---

label.c

Go to the documentation of this file.

#include        "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1990 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function segments an image into regions of similar graylevels       */
/*   and creates a label map. Two regions with mean values mean1 and mean2,   */
/*   and a combined standard deviation dev, are merged if                     */
/*                                                                            */
/*           abs(mean1 - mean2)  <=  thresh * (1 - dev/sigma)                 */
/*                                                                            */
/*   where thresh and sigma are input parameters. Assuming the image is       */
/*   segmented into N regions, the resulting label map consists of N          */
/*   segments, numbered consecutively from 0 to N-1.                          */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Background Information-----------------------------------------------------*/
/*                                                                            */
/*   Mehldau, G.:                                                             */
/*   "An Image Segmentation Algorithm."                                       */
/*   University of Arizona, DIAL, Technical Report No. 90-3                   */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
LABEL (IMAGE * in,              /*  I   Pointer to the input image.                           */
       double sigma,            /*  I   Greylevel variance threshold.                         */
       PIXEL thresh,            /*  I   Graylevel difference threshold.                       */
       IMAGE ** out             /*  O   Address of a pointer to the output image.             */
/*----------------------------------------------------------------------------*/
  )
{
  register short t, l, i, j, k;
  register double p, q, r, s, topn, lftn, topsum, lftsum, pinc, psum;
  char msg[SLEN];
  long total, *cnt = 0, *eqv = 0, top, lft, newLabel, label = 0L;
  float *mean = 0, *dev = 0;
  double *topmean = 0, *topdev = 0, *lftmean = 0, *lftdev = 0;

  if (TIMES)
    TIMING (T_LABEL);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "LABEL");
      MESSAGE ('I', "");
      sprintf (msg, " Input image:                      %s", in->text);
      MESSAGE ('I', msg);
      sprintf (msg, " Graylevel difference threshold: %12.4e", thresh);
      MESSAGE ('I', msg);
      sprintf (msg, " Graylevel variance threshold:   %12.4e", sigma);
      MESSAGE ('I', msg);
      MESSAGE ('I', " ...............");
    }

  /* check input */
  if (!CHECKIMG (in))
    {
      MESSAGE ('E', " Can't identify image.");
      goto the_end;
    }
  else if (sigma < 0.0)
    {
      MESSAGE ('E', " Threshold  must be greater than or equal to zero.");
      goto the_end;
    }
  else if (thresh < 0.0)
    {
      MESSAGE ('E', " Threshold  must be greater than or equal to zero.");
      goto the_end;
    }

  /* allocate space */
  total = (long) in->nlin * (long) in->npix;
  if (!(cnt = (long *) malloc (total * sizeof (long))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(eqv = (long *) malloc (total * sizeof (long))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  total *= (long) in->nbnd;
  if (!(mean = (float *) malloc (total * sizeof (float))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(dev = (float *) malloc (total * sizeof (float))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(topmean = (double *) malloc (in->nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(topdev = (double *) malloc (in->nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(lftmean = (double *) malloc (in->nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(lftdev = (double *) malloc (in->nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }

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

  /* initialize progress indicator */
  if (LINES && !PROGRESS (psum = 0.0))
    goto the_end;
  pinc = 1.0 / (double) in->nlin;

  for (j = 0; j < in->nlin; j++)
    {
      for (k = 0; k < in->npix; k++)
        {

          /* look at adjacent segments */
          if (t = j)
            {
              for (top = eqv[(long) (*out)->data[0][j - 1][k]];
                   top != eqv[top]; top = eqv[top]);
              eqv[(long) (*out)->data[0][j - 1][k]] = top;
              topn = (double) cnt[top];
              top *= in->nbnd;
              topsum = 0.0;
            }
          if (l = k)
            {
              for (lft = eqv[(long) (*out)->data[0][j][k - 1]];
                   lft != eqv[lft]; lft = eqv[lft]);
              eqv[(long) (*out)->data[0][j][k - 1]] = lft;
              lftn = (double) cnt[lft];
              lft *= in->nbnd;
              lftsum = 0.0;
            }

          /* merge pixel with adjacent segment? */
          for (i = 0; (t || l) && i < in->nbnd; i++)
            {
              p = in->data[i][j][k];
              if (t)
                {
                  q = mean[top + i];
                  r = dev[top + i];
                  topmean[i] = (p + topn * q) / (topn + 1.0);
                  topdev[i] =
                    sqrt (topn * (r * r + (p - q) * (p - q) / (topn + 1.0)) /
                          (topn + 1.0));
                  t = fabs (p - q) <= thresh * (1.0 - topdev[i] / sigma);
                  topsum += (p - q) * (p - q);
                }
              if (l)
                {
                  q = mean[lft + i];
                  r = dev[lft + i];
                  lftmean[i] = (p + lftn * q) / (lftn + 1.0);
                  lftdev[i] =
                    sqrt (lftn * (r * r + (p - q) * (p - q) / (lftn + 1.0)) /
                          (lftn + 1.0));
                  l = fabs (p - q) <= thresh * (1.0 - lftdev[i] / sigma);
                  lftsum += (p - q) * (p - q);
                }
            }

          /* assign pixel label */
          if (t && l)
            {
              l = !(t = topsum < lftsum);
            }
          if (t)
            {
              for (i = 0; i < in->nbnd; i++)
                {
                  mean[top + i] = topmean[i];
                  dev[top + i] = topdev[i];
                }
              cnt[top / in->nbnd] += 1L;
              (*out)->data[0][j][k] = (*out)->data[0][j - 1][k];
            }
          else if (l)
            {
              for (i = 0; i < in->nbnd; i++)
                {
                  mean[lft + i] = lftmean[i];
                  dev[lft + i] = lftdev[i];
                }
              cnt[lft / in->nbnd] += 1L;
              (*out)->data[0][j][k] = (*out)->data[0][j][k - 1];
            }
          else
            {
              newLabel = label * in->nbnd;
              for (i = 0; i < in->nbnd; i++)
                {
                  mean[newLabel + i] = in->data[i][j][k];
                  dev[newLabel + i] = 0.0;
                }
              cnt[label] = 1L;
              eqv[label] = label;
              (*out)->data[0][j][k] = (PIXEL) label++;
            }

          /* merge adjacent segments? */
          if ((t || l) && top != lft && j && k)
            {
              for (t = TRUE, i = 0; t && i < in->nbnd; i++)
                {
                  p = mean[top + i];
                  q = mean[lft + i];
                  r = dev[top + i];
                  s = dev[lft + i];
                  topmean[i] = (topn * p + lftn * q) / (topn + lftn);
                  topdev[i] =
                    sqrt ((topn * (r * r + p * p) +
                           lftn * (q * q + s * s)) / (topn + lftn) -
                          topmean[i] * topmean[i]);
                  t = fabs (p - q) <= thresh * (1.0 - topdev[i] / sigma);
                }

              /* merge adjacent segments */
              if (t)
                {
                  if (top < lft)
                    {
                      for (i = 0; i < in->nbnd; i++)
                        {
                          mean[top + i] = topmean[i];
                          dev[top + i] = topdev[i];
                        }
                      top /= in->nbnd;
                      lft /= in->nbnd;
                      cnt[top] += cnt[lft];
                      eqv[lft] = top;
                    }
                  else
                    {
                      for (i = 0; i < in->nbnd; i++)
                        {
                          mean[lft + i] = topmean[i];
                          dev[lft + i] = topdev[i];
                        }
                      top /= in->nbnd;
                      lft /= in->nbnd;
                      cnt[lft] += cnt[top];
                      eqv[top] = lft;
                    }
                }
            }

        }
      if (LINES && !PROGRESS (psum += pinc))
        goto the_end;
    }

  /* assign consecutive labels */
  for (top = 0L; top < label; top++)
    {
      if (eqv[top] != top)
        {
          for (newLabel = eqv[top]; newLabel != eqv[newLabel];
               newLabel = eqv[newLabel]);
          eqv[top] = newLabel;
        }
    }
  for (newLabel = top = 0L; top < label; top++)
    {
      eqv[top] = (eqv[top] == top ? newLabel++ : eqv[eqv[top]]);
    }
  for (j = 0; j < in->nlin; j++)
    {
      for (k = 0; k < in->npix; k++)
        {
          (*out)->data[0][j][k] = (PIXEL) eqv[(long) (*out)->data[0][j][k]];
        }
    }
  (*out)->gmin = (PIXEL) 0L;
  (*out)->gmax = (PIXEL) (newLabel - 1L);
  MESSAGE ('I', "");
  sprintf (msg, " Number of segments = %ld", newLabel);
  MESSAGE ('I', msg);

the_end:
  if (cnt)
    free (cnt);
  if (eqv)
    free (eqv);
  if (mean)
    free (mean);
  if (dev)
    free (dev);
  if (topmean)
    free (topmean);
  if (topdev)
    free (topdev);
  if (lftmean)
    free (lftmean);
  if (lftdev)
    free (lftdev);
  if (TIMES)
    TIMING (T_EXIT);
}

---