---

vstretch.c

Go to the documentation of this file.

#include        "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function performs a spatially-variable graylevel transformation     */
/*   on a single-band image.                                                  */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Background Information-----------------------------------------------------*/
/*                                                                            */
/*   Fahnestock, J.D. and Schowengerdt, R.A.:                                 */
/*   "Spatially variant contrast enhancement using local range modification." */
/*   Optical Engineering, Vol. 22, No. 3, May/June 1983                       */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
VSTRETCH (IMAGE * in,           /*  I   Pointer to the input image.                           */
          short size,           /*  I   Number of lines and pixels in image blocks.           */
          IMAGE ** out          /*  O   Address of a pointer to the output image.             */
/*----------------------------------------------------------------------------*/
  )
{
  register short j, k, u, v;
  char msg[SLEN];
  long total, x, y, idx1, idx2, blkx, blky;
  double rat1, rat2, rat3, rat4, pinc, psum;
  PIXEL *hi = 0, *lo = 0, *mp = 0, *np = 0, m, n;

  if (TIMES)
    TIMING (T_VSTRETCH);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "VSTRETCH");
      MESSAGE ('I', "");
      sprintf (msg, " Input image:                    %s", in->text);
      MESSAGE ('I', msg);
      sprintf (msg, " Block size: lines and pixels:   %d", size);
      MESSAGE ('I', msg);
      MESSAGE ('I', " ...............");
    }

  /* check input */
  if (!CHECKIMG (in))
    {
      MESSAGE ('E', " Can't identify image.");
      goto the_end;
    }
  else if (in->nbnd > 1)
    {
      MESSAGE ('E', " Image must be single-band.");
      goto the_end;
    }
  else if (size < 1)
    {
      MESSAGE ('E',
               " Size of subblock must be equal to or greater than one.");
      goto the_end;
    }
  else if (size > in->nlin || size > in->npix)
    {
      MESSAGE ('E',
               " Image size must be equal to or greater than subblock size.");
      goto the_end;
    }

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

  /* allocate space for arrays */
  total = blkx * blky;
  if (!(hi = (PIXEL *) malloc (total * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(lo = (PIXEL *) malloc (total * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  total = (blkx + 1L) * (blky + 1L);
  if (!(mp = (PIXEL *) malloc (total * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(np = (PIXEL *) malloc (total * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }

  /* initialize progress indicator */
  if (LINES && !PROGRESS (psum = 0.0))
    goto the_end;
  pinc = 0.5 / (double) blkx / (double) blky;

  /* first pass: obtain transformation parameters from subblocks */
  /* search for maximum and minimum */
  for (y = 0L; y < blky; y++)
    {
      for (x = 0L; x < blkx; x++)
        {
          for (m = n = in->data[0][y * size][x * size], j = y * size;
               j < (y + 1L) * size && j < in->nlin; j++)
            {
              for (k = x * size; k < (x + 1L) * size && k < in->npix; k++)
                {
                  if (in->data[0][j][k] > m)
                    m = in->data[0][j][k];
                  if (in->data[0][j][k] < n)
                    n = in->data[0][j][k];
                }
            }
          hi[y * blkx + x] = m;
          lo[y * blkx + x] = n;
          if (LINES && !PROGRESS (psum += pinc))
            goto the_end;
        }
    }

  /* obtain values for parameter arrays */
  /* four corners */
  mp[0] = hi[0];
  np[0] = lo[0];
  mp[blkx] = hi[blkx - 1L];
  np[blkx] = lo[blkx - 1L];
  mp[blky * (blkx + 1L)] = hi[(blky - 1L) * blkx];
  np[blky * (blkx + 1L)] = lo[(blky - 1L) * blkx];
  mp[(blkx + 1L) * (blky + 1L) - 1L] = hi[blkx * blky - 1L];
  np[(blkx + 1L) * (blky + 1L) - 1L] = lo[blkx * blky - 1L];

  /* four sides */
  for (y = 1L; y < blky; y++)
    {
      idx1 = (y - 1L) * blkx;
      idx2 = y * blkx;
      mp[y * (blkx + 1L)] = max (hi[idx1], hi[idx2]);
      np[y * (blkx + 1L)] = min (lo[idx1], lo[idx2]);
      idx1 = (y + 1L) * blkx - 1L;
      mp[(y + 1L) * (blkx + 1L) - 1L] = max (hi[idx2 - 1L], hi[idx1]);
      np[(y + 1L) * (blkx + 1L) - 1L] = min (lo[idx2 - 1L], lo[idx1]);
    }
  for (x = 1L; x < blkx; x++)
    {
      mp[x] = max (hi[x - 1L], hi[x]);
      np[x] = min (lo[x - 1L], lo[x]);
      idx1 = blky * (blkx + 1L) + x;
      idx2 = (blky - 1L) * blkx + x;
      mp[idx1] = max (hi[idx2 - 1L], hi[idx2]);
      np[idx1] = min (lo[idx2 - 1L], lo[idx2]);
    }

  /* inside */
  for (y = 1L; y < blky; y++)
    {
      for (x = 1L; x < blkx; x++)
        {
          idx1 = (y - 1L) * blkx + x;
          idx2 = y * blkx + x;
          mp[y * (blkx + 1L) + x] =
            max (max (hi[idx1], hi[idx1 - 1L]),
                 max (hi[idx2], hi[idx2 - 1L]));
          np[y * (blkx + 1L) + x] =
            min (min (lo[idx1], lo[idx1 - 1L]),
                 min (lo[idx2], lo[idx2 - 1L]));
        }
    }

  /* second pass: perform transformation */
  for (x = 0L; x < blkx; x++)
    {
      for (y = 0L; y < blky; y++)
        {
          for (u = 0, j = y * size; j < (y + 1L) * size && j < in->nlin;
               u++, j++)
            {
              for (v = 0, k = x * size; k < (x + 1L) * size && k < in->npix;
                   v++, k++)
                {
                  idx1 = (y + 1L) * (blkx + 1L) + x;
                  idx2 = y * (blkx + 1L) + x;
                  rat1 = (double) u / (double) size;
                  rat2 = (double) (size - u) / (double) size;
                  rat3 = (double) (size - v) / (double) size;
                  rat4 = (double) v / (double) size;
                  m =
                    (rat1 * (double) mp[idx1] +
                     rat2 * (double) mp[idx2]) * rat3 +
                    (rat1 * (double) mp[idx1 + 1L] +
                     rat2 * (double) mp[idx2 + 1L]) * rat4;
                  n =
                    (rat1 * (double) np[idx1] +
                     rat2 * (double) np[idx2]) * rat3 +
                    (rat1 * (double) np[idx1 + 1L] +
                     rat2 * (double) np[idx2 + 1L]) * rat4;
                  (*out)->data[0][j][k] =
                    m ==
                    n ? in->
                    data[0][j][k] : (PIXEL) ((double) (in->gmax - in->gmin) /
                                             (double) (m -
                                                       n) *
                                             ((double) in->data[0][j][k] -
                                              n)) + in->gmin;
                }
            }
          if (LINES && !PROGRESS (psum += pinc))
            goto the_end;
        }
    }

the_end:
  if (hi)
    free (hi);
  if (lo)
    free (lo);
  if (mp)
    free (mp);
  if (np)
    free (np);
  if (TIMES)
    TIMING (T_EXIT);
}

---