---

match.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 transforms the graylevels of one image to match those      */
/*   of a second image.  The first option uses a least-squares match, and     */
/*   the second option uses the cumulative distribution functions of both     */
/*   images.                                                                  */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Background Information-----------------------------------------------------*/
/*                                                                            */
/*   Dallas, W.J. and Mauser, W.:                                             */
/*   "Preparing Pictures for Visual Comparison."                              */
/*   Applied Optics, Vol. 19, No. 21, November 1980                           */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
MATCH (IMAGE * in1,             /*  I  Pointer to the transformation image.                  */
       IMAGE * in2,             /*  I  Pointer to the reference image.                       */
       short nlev,              /*  I  Number of graylevels to be used in the transformation. */
       short option,            /*  I  Switch, indicating type of image matching             */
       /*     LSQ   -   Images matched using gray levels            */
       /*     CDF   -   Images matched using cumulative dist fncns  */
       IMAGE ** out             /*  O  Address of a pointer to the output image.             */
       /*     The input image may be overwritten (out = &in1).      */
/*----------------------------------------------------------------------------*/
  )
{
  register short i, j, k, index, index2;
  char msg[SLEN];
  long *count = 0, *hist1 = 0, *hist2 = 0;
  double offset, offset2, factor, factor2, pinc, psum, *cdf1 = 0, *cdf2 = 0;
  PIXEL *table = 0;
  char minlbl[10], maxlbl[10];

  if (TIMES)
    TIMING (T_MATCH);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "MATCH");
      MESSAGE ('I', "");
      sprintf (msg, " Transformation input image:   %s", in1->text);
      MESSAGE ('I', msg);
      sprintf (msg, " Reference input image:        %s", in2->text);
      MESSAGE ('I', msg);
      sprintf (msg, " Lookup table entries:         %d", nlev);
      MESSAGE ('I', msg);
      if (option == LSQ)
        {
          MESSAGE ('I', " Least squares transformation");
        }
      else
        {
          MESSAGE ('I', " CDF transformation");
        }
      MESSAGE ('I', " ...............");
    }

  /* check input */
  if (!CHECKIMG (in1))
    {
      MESSAGE ('E', " Can't identify transformation image.");
      goto the_end;
    }
  else if (!CHECKIMG (in2))
    {
      MESSAGE ('E', " Can't identify reference image.");
      goto the_end;
    }
  else if (in1->nbnd != in2->nbnd)
    {
      MESSAGE ('E', " Number of bands must be identical in both images.");
      goto the_end;
    }
  else if (in1->nlin != in2->nlin)
    {
      MESSAGE ('E', " Number of lines must be identical in both images.");
      goto the_end;
    }
  else if (in1->npix != in2->npix)
    {
      MESSAGE ('E',
               " Number of pixels/line must be identical in both images.");
      goto the_end;
    }
  else if (nlev <= 0)
    {
      MESSAGE ('E', " Number of graylevel bins must be greater than zero.");
      goto the_end;
    }
  else if (option != LSQ && option != CDF)
    {
      MESSAGE ('E', " Option must be either LSQ or CDF.");
      goto the_end;
    }

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

  /* initialize progress indicator */
  if (LINES && !PROGRESS (psum = 0.0))
    goto the_end;
  pinc =
    (option == LSQ ? 0.5 : 1.0) / (double) in1->nbnd / (double) in1->nlin;

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

  /* initialize */
  RANGE (in1);
  offset = (double) in1->gmin;
  factor = (double) (nlev - 1) / (double) (in1->gmax - in1->gmin);

  if (option == LSQ)
    {

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

      for (i = 0; i < in1->nbnd; i++)
        {

          /* build look-up table */
          for (index = 0; index < nlev; index++)
            {
              table[index] = (PIXEL) 0;
              count[index] = 0L;
            }
          for (j = 0; j < in1->nlin; j++)
            {
              for (k = 0; k < in1->npix; k++)
                {
                  index =
                    (short) (((double) in1->data[i][j][k] - offset) * factor);
                  table[index] += in2->data[i][j][k];
                  count[index] += 1L;
                }
              if (LINES && !PROGRESS (psum += pinc))
                goto the_end;
            }

          /* average the table with the count array */
          for (index = 0; index < nlev; index++)
            {
              if (count[index] != 0)
                {
                  table[index] /= (double) count[index];
                }
            }

          /* create output image */
          for (j = 0; j < in1->nlin; j++)
            {
              for (k = 0; k < in1->npix; k++)
                {
                  index =
                    (short) (((double) in1->data[i][j][k] - offset) * factor);
                  (*out)->data[i][j][k] = table[index];
                }
              if (LINES && !PROGRESS (psum += pinc))
                goto the_end;
            }

          /* output look-up table */
          sprintf (minlbl, "0");
          sprintf (maxlbl, "%d", nlev);
          PLOT (&table, 1, nlev, minlbl, maxlbl, NORM);
        }

    }
  else
    {

      /* allocate space for histograms of reference and transform images and cdf's */
      if (!(hist1 = (long *) malloc (in1->nbnd * nlev * sizeof (long))))
        {
          MESSAGE ('E', " Memory request failed.");
          goto the_end;
        }
      if (!(hist2 = (long *) malloc (in2->nbnd * nlev * sizeof (long))))
        {
          MESSAGE ('E', " Memory request failed.");
          goto the_end;
        }
      if (!(cdf1 = (double *) malloc (in1->nbnd * nlev * sizeof (double))))
        {
          MESSAGE ('E', " Memory request failed.");
          goto the_end;
        }
      if (!(cdf2 = (double *) malloc (in2->nbnd * nlev * sizeof (double))))
        {
          MESSAGE ('E', " Memory request failed.");
          goto the_end;
        }

      /* compute cumulative distribution function for transformation image */
      RANGE (in1);
      HISTOGRM (in1, 1, in1->gmin, in1->gmax, nlev, hist1);
      for (i = 0; i < in1->nbnd; i++)
        {
          for (j = 1; j < nlev; j++)
            {
              hist1[i * nlev + j] += hist1[i * nlev + j - 1];
            }
        }

      /* scale cumulative distribution function to [0,1] */
      for (i = 0; i < in1->nbnd; i++)
        {
          for (j = 0; j < nlev; j++)
            {
              cdf1[i * nlev + j] =
                (double) (hist1[i * nlev + j]) / (double) (hist1[nlev - 1]);
            }
        }

      /* compute cumulative distribution function for reference image */
      RANGE (in2);
      HISTOGRM (in2, 1, in2->gmin, in2->gmax, nlev, hist2);
      for (i = 0; i < in2->nbnd; i++)
        {
          for (j = 1; j < nlev; j++)
            {
              hist2[i * nlev + j] += hist2[i * nlev + j - 1];
            }
        }

      /* scale cumulative distribution function to [0,1] */
      for (i = 0; i < in1->nbnd; i++)
        {
          for (j = 0; j < nlev; j++)
            {
              cdf2[i * nlev + j] =
                (double) (hist2[i * nlev + j]) / (double) (hist2[nlev - 1]);
            }
        }

      offset2 = (double) in2->gmin;
      factor2 = (double) (nlev - 1) / (double) (in2->gmax - in2->gmin);

      /* loop through histograms making the lookup table entries table[index] = cdf2 <inverse> [cdf1[index]] */
      for (i = 0; i < in1->nbnd; i++)
        {
          for (index = 0; index < nlev; index++)
            {
              for (index2 = 0; index2 < nlev; index2++)
                {
                  if (cdf2[i * nlev + index2] == cdf1[i * nlev + index])
                    {
                      table[index] =
                        (PIXEL) index2 / (PIXEL) factor2 + offset2;
                      break;
                    }
                  else if (cdf2[i * nlev + index2] > cdf1[i * nlev + index])
                    {
                      if (index2 == 0)
                        {
                          table[index] = offset2;
                        }
                      else
                        {
                          if ((cdf2[i * nlev + index2] -
                               cdf1[i * nlev + index]) <=
                              (cdf1[i * nlev + index] -
                               cdf2[i * nlev + index2 - 1]))
                            {
                              table[index] =
                                (PIXEL) index2 / (PIXEL) factor2 + offset2;
                            }
                          else
                            {
                              table[index] =
                                (PIXEL) (index2 - 1) / (PIXEL) factor2 +
                                offset2;
                            }
                        }
                      break;
                    }
                }
            }

          /* create output image */
          for (j = 0; j < in1->nlin; j++)
            {
              for (k = 0; k < in1->npix; k++)
                {
                  index =
                    (short) (((double) (in1->data[i][j][k] - offset) *
                              factor) + 0.5);
                  (*out)->data[i][j][k] = table[index];
                }
              if (LINES && !PROGRESS (psum += pinc))
                goto the_end;
            }

          /* output the lookup table */
          sprintf (minlbl, "0");
          sprintf (maxlbl, "%d", nlev);
          PLOT (&table, 1, nlev, minlbl, maxlbl, NORM);
        }
    }

the_end:
  if (table)
    free (table);
  if (count)
    free (count);
  if (hist1)
    free (hist1);
  if (hist2)
    free (hist2);
  if (cdf1)
    free (cdf1);
  if (cdf2)
    free (cdf2);
  if (TIMES)
    TIMING (T_EXIT);
}

---