---

destripe.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--------------------------------------------------------*/
/*                                                                            */
/*   There are three methods to use in destriping an image.                   */
/*   All use either a single detector or the average of all detectors as      */
/*   a reference.  The first method adjusts the means of all the detectors    */
/*   to match that of the reference.  The second method adjusts the means and */
/*   standard deviations to match those of the reference.  The third method   */
/*   computes cumulative distribution functions using a specified number of   */
/*   levels for the reference and each of the detectors.  These are then used */
/*   in matching the contrast of the detectors to that of the reference.      */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
DESTRIPE (IMAGE * in,           /*  I  Pointer to the image.                                 */
          short option,         /*  I  Switch, indicating type of line matching              */
          /*     DESMNE - Destripe using mean only levels              */
          /*     DESSTD - Destripe using mean and standard dev.        */
          /*     DESCDF - Destripe using cumulative dist. function     */
          short roption,        /*  I  Switch, indicating type of reference used             */
          /*     DESSNG - Use a single line (detector) as reference    */
          /*     DESAVG - Average all detectors for reference          */
          short lines,          /*  I  Number of lines to destripe - ie, number of           */
          /*     detectors.  Must be a factor of total image line num  */
          short refline,        /*  I  The line (detector) to be used for destriping.        */
          /*     This value does not matter if DESAVG is selected.     */
          short nlev,           /*  I  The number of levels used in the cumulative           */
          /*     distribution function for the reference and           */
          /*     transformation lines.                                 */
          IMAGE ** out          /*  O  Address of a pointer to the output image.             */
          /*     The input image may be overwritten (out = &in1).      */
/*----------------------------------------------------------------------------*/
  )
{
  register short i, j, k, l;
  short skip, index, index2;
  double mean, stddev, oldmean, oldstddev, factor, *rcdf = 0, *tcdf = 0;
  PIXEL *table = 0;

  if (TIMES)
    TIMING (T_DESTRIPE);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "DESTRIPE");
    }

  /* check input */
  if (!CHECKIMG (in))
    {
      MESSAGE ('E', " Can't identify image.");
      goto the_end;
    }
  else if (option != DESMNE && option != DESSTD && option != DESCDF)
    {
      MESSAGE ('E', " Option must be DESMNE, DESSTD, or DESCDF.");
      goto the_end;
    }
  else if (roption != DESSNG && roption != DESAVG)
    {
      MESSAGE ('E', " Reference option must be either DESSNG or DESAVG.");
      goto the_end;
    }

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

  /* Set up image hopping parameters */
  if (roption == DESAVG)
    {
      refline = 1;
      skip = 1;
    }
  else
    {
      skip = lines;
    }

  if (option == DESMNE || option == DESSTD)
    {

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

          /* compute reference mean */
          mean = 0.0;
          for (j = refline - 1; j < in->nlin; j += skip)
            {
              for (k = 0; k < in->npix; k++)
                {
                  mean += (double) in->data[i][j][k];
                }
            }
          if (roption == DESSNG)
            {
              mean /= (double) in->npix * (double) (j - refline +
                                                    1) / (double) lines;
            }
          else
            {
              mean /= (double) in->npix * (double) in->nlin;
            }

          /* compute reference standard deviation */
          if (option == DESSTD)
            {
              stddev = 0.0;
              for (j = refline - 1; j < in->nlin; j += skip)
                {
                  for (k = 0; k < in->npix; k++)
                    {
                      stddev +=
                        ((double) in->data[i][j][k] -
                         mean) * ((double) in->data[i][j][k] - mean);
                    }
                }
              if (roption == DESSNG)
                {
                  stddev /= (double) in->npix * (double) (j - refline +
                                                          1) / (double) lines;
                }
              else
                {
                  stddev /= (double) in->npix * (double) in->nlin;
                }
              stddev = sqrt (stddev);
            }

          /* Calculate mean and/or standard deviation of each line to be transformed. Copy reference line to output. */
          for (j = 0; j < lines; j++)
            {

              if (roption == DESSNG && j == refline - 1)
                {
                  /* copy reference lines (detector) directly to output image */
                  for (k = j; k < in->nlin; k += lines)
                    {
                      for (l = 0; l < in->npix; l++)
                        {
                          (*out)->data[i][k][l] = in->data[i][k][l];
                        }
                    }
                }
              else
                {
                  /* compute mean of one set of lines (one detector) */
                  oldmean = 0.0;
                  for (k = j; k < in->nlin; k += lines)
                    {
                      for (l = 0; l < in->npix; l++)
                        {
                          oldmean += (double) in->data[i][k][l];
                        }
                    }
                  oldmean /= (double) in->npix * (double) (k -
                                                           j) /
                    (double) lines;

                  /* compute standard deviation of one set of lines (one detector) */
                  if (option == DESSTD)
                    {
                      oldstddev = 0.0;
                      for (k = j; k < in->nlin; k += lines)
                        {
                          for (l = 0; l < in->npix; l++)
                            {
                              oldstddev +=
                                ((double) in->data[i][k][l] -
                                 oldmean) * ((double) in->data[i][k][l] -
                                             oldmean);
                            }
                        }
                      oldstddev /= (double) in->npix * (double) (k -
                                                                 j) /
                        (double) lines;
                      oldstddev = sqrt (oldstddev);
                    }

                  /* copy transformed line to output, adjusting gray levels to achieve mean and/or standard deviation of reference */
                  for (k = j; k < in->nlin; k += lines)
                    {
                      for (l = 0; l < in->npix; l++)
                        {
                          if (option == DESMNE)
                            {
                              (*out)->data[i][k][l] =
                                (PIXEL) ((double) in->data[i][k][l] -
                                         oldmean + mean);
                            }
                          else
                            {
                              (*out)->data[i][k][l] =
                                (PIXEL) (((double) in->data[i][k][l] -
                                          oldmean) * stddev / oldstddev +
                                         mean);
                            }
                        }
                    }
                }
            }
        }

    }
  else
    {

      if (!(rcdf = (double *) malloc (nlev * sizeof (double))))
        {
          MESSAGE ('E', " Memory request failed.");
          goto the_end;
        }
      if (!(tcdf = (double *) malloc (nlev * sizeof (double))))
        {
          MESSAGE ('E', " Memory request falied.");
          goto the_end;
        }
      if (!(table = (PIXEL *) malloc (nlev * sizeof (PIXEL))))
        {
          MESSAGE ('E', " Memory request failed.");
          goto the_end;
        }

      RANGE (in);
      factor = (double) (nlev - 1) / (double) (in->gmax - in->gmin);

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

          /* compute histogram of reference */
          for (j = 0; j < nlev; j++)
            {
              rcdf[j] = 0.0;
            }
          for (j = refline - 1; j < in->nlin; j += skip)
            {
              for (k = 0; k < in->npix; k++)
                {
                  rcdf[(short)
                       ((double) (in->data[i][j][k] - in->gmin) * factor +
                        0.5)] += 1.0;
                }
            }

          /* convert histogram to cumulative distribution function and scale to [0,1] */
          for (j = 1; j < nlev; j++)
            {
              rcdf[j] += rcdf[j - 1];
            }
          for (j = 0; j < nlev; j++)
            {
              rcdf[j] /= rcdf[nlev - 1];
            }

          /* loop through each line, calculate cdf, and use it to transform the line */
          for (j = 0; j < lines; j++)
            {

              /* copy reference lines (detector) directly to output image */
              if (roption == DESSNG && j == refline - 1)
                {
                  for (k = j; k < in->nlin; k += lines)
                    {
                      for (l = 0; l < in->npix; l++)
                        {
                          (*out)->data[i][k][l] = in->data[i][k][l];
                        }
                    }
                }
              else
                {

                  /* compute histogram of one set of lines (one detector) */
                  for (k = 0; k < nlev; k++)
                    {
                      tcdf[k] = 0.0;
                    }
                  for (k = j; k < in->nlin; k += lines)
                    {
                      for (l = 0; l < in->npix; l++)
                        {
                          tcdf[(short)
                               ((double) (in->data[i][k][l] - in->gmin) *
                                factor + 0.5)] += 1;
                        }
                    }

                  /* convert histogram to cumulative distribution function amd scale to [0,1] */
                  for (k = 1; k < nlev; k++)
                    {
                      tcdf[k] += tcdf[k - 1];
                    }
                  for (k = 0; k < nlev; k++)
                    {
                      tcdf[k] /= tcdf[nlev - 1];
                    }

                  /* loop through histograms making lookup table entries */
                  for (index = 0; index < nlev; index++)
                    {
                      for (index2 = 0; index2 < nlev; index2++)
                        {
                          if (rcdf[index2] == tcdf[index])
                            {
                              table[index] =
                                (PIXEL) index2 / (PIXEL) factor + in->gmin;
                              break;
                            }
                          else if (rcdf[index2] > tcdf[index])
                            {
                              if (index2 == 0)
                                {
                                  table[index] = in->gmin;
                                }
                              else
                                {
                                  if (rcdf[index2] - tcdf[index] <=
                                      tcdf[index] - rcdf[index2 - 1])
                                    {
                                      table[index] =
                                        (PIXEL) index2 / (PIXEL) factor +
                                        in->gmin;
                                    }
                                  else
                                    {
                                      table[index] =
                                        (PIXEL) (index2 -
                                                 1) / (PIXEL) factor +
                                        in->gmin;
                                    }
                                }
                              break;
                            }
                        }
                    }

                  /* create new set of lines (detector) for output image */
                  for (k = j; k < in->nlin; k += lines)
                    {
                      for (l = 0; l < in->npix; l++)
                        {
                          index =
                            (short) ((double) (in->data[i][k][l] - in->gmin) *
                                     factor + 0.5);
                          (*out)->data[i][k][l] = table[index];
                        }
                    }
                }
            }
        }
    }

the_end:
  if (rcdf)
    free (rcdf);
  if (tcdf)
    free (tcdf);
  if (table)
    free (table);
  if (TIMES)
    TIMING (T_EXIT);
}

---