home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_7.1 / powspec / powspec.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  13.4 KB  |  391 lines
  1. /* 
  2.  * powspec.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /* POWSPEC: program performs 2-dimensional FFT on image and yields
  10.  *          image containing power spectrum
  11.  *              usage: powspec inimg outimg [-l] [-p] [-w] [-c] [-L]
  12.  *        Note: image sidelength must be power of 2; if not the case,
  13.  *              this program changes the size to be so, either the
  14.  *              next higher or lower power of 2 as user-chosen (<-p>).
  15.  *
  16.  */
  17.  
  18. #include <stdio.h>
  19. #include <math.h>
  20. #include <stdlib.h>
  21. #include <malloc.h>
  22. #include <string.h>
  23.  
  24. #include "ip.h"
  25.  
  26. #if defined(WIN32) || defined(LINUX)
  27. #define M_PI    3.14159265358979323846
  28. #define log2(a) (log(a)/ 0.6931471805599)
  29. #define exp2(a) (pow((double)2, (double)a))
  30. #endif
  31. #define CROSSSEC_HEIGHT 128     /* height of cross section display */
  32.  
  33. void input (int argc, char *argv[], short *linear, short *zeroPad, short *noWndw, short *crossSec);
  34. void window (float **image, long nRow, long nCol, long nRowS, long nColS);
  35. void usage (short flag);
  36.  
  37. int
  38. main (argc, argv)
  39.      int argc;
  40.      char *argv[];
  41.  
  42. {
  43.   register long x, y, nCol, nRow,  /* image input sidelengths */
  44.     nCol2, nRow2,               /* power of 2 spectrum sidelengths */
  45.     nColD2, nRowD2;             /* middle row/column of spectrum image */
  46.   long xStart, yStart;          /* beginning coord.s of pow spec in image */
  47.   Image *imgI, *imgO;           /* I/O image structures */
  48.   unsigned char **imgIn,        /* image array */
  49.   **imgOut;                     /* output image */
  50.   short linear;                 /* if 1, linear power scale; 0 => log */
  51.   short zeroPad;                /* if 1, upsize image to pow of 2; 0 => down */
  52.   short noWndw;                 /* if 1, no smoothing window; 0 otherwise */
  53.   short crossSec;               /* if 1, display cross section; 0 otherwise */
  54.   float **imgReal, **imgImag;   /* complex image arrays */
  55.   //float *wndwX, *wndwY;     /* row and column window vectors */
  56.   double min, max;              /* maximum value in power spectrum */
  57.   double temp;
  58.   long sum, nSum;               /* sum of intensities */
  59.   long avg;                     /* average image intensity (dc) */
  60.   long nRow3, nRow4;            /* image heights with cross section displays */
  61.   long *radialP;                /* radial power for cross section display */
  62.   long *nRadialP;               /* no. of power measurements at each radius */
  63.   long radius;                  /* radius from (0,0) for radial power */
  64.   long *angleP;                 /* angular power for cross section display */
  65.   long *nAngleP;                /* no. of power measurements at each angle */
  66.   long angle;                   /* angle (0 to 180) degrees */
  67.   long maxR, maxA;              /* max of powers for normalization */
  68.   long minR, minA;              /* min of powers for normalization */
  69.   double xRad, yRad;            /* x,y radii from (0,0) */
  70.  
  71.   input (argc, argv, &linear, &zeroPad, &noWndw, &crossSec);
  72.  
  73. /* read input image */
  74.   imgI = ImageIn (argv[1]);
  75.   imgIn = imgI->img;
  76.   nRow = ImageGetHeight (imgI);
  77.   nCol = ImageGetWidth (imgI);
  78.   printf ("Image is %3dx%3d.\n", nCol, nRow);
  79.  
  80. /* check for image sidelengths powers of 2 */
  81.   nRow2 = nRow;
  82.   temp = log2 ((double) nRow);
  83.   if ((temp - (long) temp) != 0.0)
  84.     nRow2 = (zeroPad == 0) ?
  85.       (long) exp2 ((double) ((long) temp)) : (long) exp2 ((double) ((long) temp + 1));
  86.   nCol2 = nCol;
  87.   temp = log2 ((double) nCol);
  88.   if ((temp - (long) temp) != 0.0)
  89.     nCol2 = (zeroPad == 0) ?
  90.       (long) exp2 ((double) ((long) temp)) : (long) exp2 ((double) ((long) temp + 1));
  91.   if (nRow2 != nRow || nCol2 != nCol)
  92.     printf ("Image size changed to %dx%d so power of 2 (required for FFT).\n",
  93.             nRow2, nCol2);
  94.  
  95. /* allocate output image */
  96.   nRow3 = (crossSec) ? nRow2 + CROSSSEC_HEIGHT : nRow2;
  97.   nRow4 = (crossSec) ? nRow2 + 2 * CROSSSEC_HEIGHT : nRow2;
  98.   imgO = ImageAlloc (nRow4, nCol2, ImageGetDepth (imgI));
  99.   imgOut = ImageGetPtr (imgO);
  100.  
  101. /* allocate memory for complex array */
  102.   if ((imgReal = (float **) calloc (nRow2, sizeof (float *))) == NULL)
  103.       exit (2);
  104.   for (y = 0; y < nRow2; y++)
  105.     if ((imgReal[y] = (float *) calloc (nCol2, sizeof (float))) == NULL)
  106.         exit (3);
  107.  
  108.   if ((imgImag = (float **) calloc (nRow2, sizeof (float *))) == NULL)
  109.       exit (4);
  110.   for (y = 0; y < nRow2; y++)
  111.     if ((imgImag[y] = (float *) calloc (nCol2, sizeof (float))) == NULL)
  112.         exit (5);
  113.  
  114. /* find image average intensity */
  115.   for (y = 0, sum = 0, nSum = 0; y < nRow; y++) {
  116.     for (x = 0; x < nCol; x++) {
  117.       sum += (long) imgIn[y][x];
  118.       nSum++;
  119.     }
  120.   }
  121.   avg = sum / nSum;
  122.  
  123.  
  124. /* center the image, remove average, and zero-pad or reduce size if required */
  125.   printf ("Forward 2-D FFT being performed.\n");
  126.   if (zeroPad == 0) {
  127.     yStart = (nRow - nRow2) / 2;
  128.     xStart = (nCol - nCol2) / 2;
  129.     for (y = 0; y < nRow2; y++) {
  130.       for (x = 0; x < nCol2; x++) {
  131.         imgReal[y][x] = (float) imgIn[y + yStart][x + xStart] - (float) avg;
  132.         imgImag[y][x] = (float) 0.0;
  133.       }
  134.     }
  135.   }
  136.   else {
  137.     yStart = (nRow2 - nRow) / 2;
  138.     xStart = (nCol2 - nCol) / 2;
  139.     for (y = 0; y < nRow; y++) {
  140.       for (x = 0; x < nCol; x++) {
  141.         imgReal[y + yStart][x + xStart] = (float) imgIn[y][x];
  142.         imgImag[y][x] = (float) 0.0;
  143.       }
  144.     }
  145.   }
  146. /* perform windowing */
  147.   if (zeroPad == 0)
  148.     window (imgReal, nRow2, nCol2, nRow2, nCol2);
  149.   else
  150.     window (imgReal, nRow2, nCol2, nRow, nCol);
  151.  
  152. /* perform 2-dimensional FFT */
  153.   fft2d (imgReal, imgImag, nRow2, nCol2, -1);
  154.  
  155. /* construct power spectrum */
  156.   printf ("Power spectrum being constructed.\n");
  157.   for (y = 0, max = 0.0, min = 100000.0; y < nRow2; y++) {
  158.     for (x = 0; x < nCol2; x++) {
  159.       imgReal[y][x] = (float) sqrt ((double) (imgReal[y][x] * imgReal[y][x]
  160.                                           + imgImag[y][x] * imgImag[y][x]));
  161.       if (!(y == 0 && x == 0))
  162.         max = (max > imgReal[y][x]) ? max : imgReal[y][x];
  163.       min = (min < imgReal[y][x]) ? min : imgReal[y][x];
  164.     }
  165.   }
  166.   if (imgReal[0][0] > max)
  167.     imgReal[0][0] = (float) max;
  168.  
  169. /* scale, center origin at middle, and write output bytes */
  170.   max = (linear == 0) ? log10 (max - min + 1.0) : max - min;
  171.   nColD2 = nCol2 / 2;
  172.   nRowD2 = nRow2 / 2;
  173.   for (y = 0; y < nRow2; y++)
  174.     for (x = 0; x < nCol2; x++) {
  175.       temp = (linear == 0) ?
  176.         log10 ((double) imgReal[y][x] - min + 1.0) :
  177.         (double) (imgReal[y][x] - min);
  178.       imgOut[(y + nRowD2) % nRow2][(x + nColD2) % nCol2] =
  179.         (unsigned char) (temp * 255.0 / max);
  180.     }
  181.  
  182. /* display cross section of horizontal line through (0,0) if flag set */
  183.   if (crossSec) {
  184.     if ((radialP = (long *) calloc (nColD2, sizeof (long))) == NULL)
  185.         exit (6);
  186.     if ((nRadialP = (long *) calloc (nColD2, sizeof (long))) == NULL)
  187.         exit (7);
  188.     if ((angleP = (long *) calloc (nColD2, sizeof (long))) == NULL)
  189.         exit (8);
  190.     if ((nAngleP = (long *) calloc (nColD2, sizeof (long))) == NULL)
  191.         exit (9);
  192.     for (x = 0; x < nColD2; x++) {
  193.       radialP[x] = 0;
  194.       nRadialP[x] = 1;
  195.       angleP[x] = 0;
  196.       nAngleP[x] = 1;
  197.     }
  198.     for (y = 0; y < nRow2; y++) {
  199.       for (x = 0; x < nCol2; x++) {
  200.         yRad = (double) (nRowD2 - y);
  201.         xRad = (double) (nColD2 - x);
  202.         radius = (long) (sqrt (yRad * yRad + xRad * xRad) + 0.5);
  203.         if (radius < nColD2)
  204.           radialP[radius] += (long) imgOut[y][x];
  205.         nRadialP[radius]++;
  206.         if (yRad == 0.0 && xRad == 0.0)
  207.           angle = 0;
  208.         else
  209.           angle = (long) (atan2 (-yRad, xRad) *
  210.                           (double) (nColD2 - 1) / M_PI + 0.5);
  211.         if (angle < 0)
  212.           angle += nColD2;
  213.         angleP[angle] += (long) imgOut[y][x];
  214.         nAngleP[angle]++;
  215.       }
  216.     }
  217.     for (x = 0, maxR = maxA = 0; x < nColD2; x++) {
  218.       if (x != 0)
  219.         radialP[x] /= nRadialP[x];
  220.       if (radialP[x] > maxR)
  221.         maxR = radialP[x];
  222.       angleP[x] /= nAngleP[x];
  223.       if (angleP[x] > maxA)
  224.         maxA = angleP[x];
  225.     }
  226.     for (x = 0, minR = maxR, minA = maxA; x < nColD2; x++) {
  227.       if (radialP[x] < minR)
  228.         minR = radialP[x];
  229.       if (angleP[x] < minA)
  230.         minA = angleP[x];
  231.     }
  232.     for (y = nRow2 + 1; y < nRow4; y++)
  233.       for (x = 0; x < nCol2; x++)
  234.         imgOut[y][x] = 255;
  235.     for (x = 0; x < nCol2; x++)
  236.       imgOut[nRow2][x] = 0;
  237.     maxR -= minR;
  238.     for (radius = 0; radius < nColD2; radius++) {
  239.       radialP[radius] = (radialP[radius] - minR) * CROSSSEC_HEIGHT / maxR;
  240.       for (y = nRow3 - 1 - radialP[radius]; y < nRow3; y++)
  241.         imgOut[y][nColD2 + radius] = imgOut[y][nColD2 - radius] = 0;
  242.     }
  243.     maxA -= minA;
  244.     for (angle = 0; angle < nColD2; angle++) {
  245.       angleP[angle] = (angleP[angle] - minA) * CROSSSEC_HEIGHT / maxA;
  246.       for (y = nRow4 - 1 - angleP[angle]; y < nRow4; y++)
  247.         imgOut[y][nColD2 + angle] = imgOut[y][nColD2 - angle] = 0;
  248.     }
  249.   }
  250.  
  251. /* output power spectrum */
  252.   ImageOut (argv[2], imgO);
  253.   return (0);
  254. }
  255.  
  256.  
  257. /* WINDOW:      function multiplies vector by smoothing window
  258.  *                    usage: window (image, nRow, nCol, nRowS, nColS)
  259.  *              Note: There are 2 sets of image sizes for the following
  260.  *              reason. The (nRow, nCol) size is the image size to be
  261.  *              windowed; that is, this is a power of 2 for the FFT.
  262.  *              The (nRowS, nColS) size is a smaller size that is the original
  263.  *              image size before zero-padding. This smaller image has
  264.  *              been placed in the middle of the zero-padded image. It is
  265.  *              the SMALLER IMAGE, NOT THE FINAL IMAGE that should be windowed.
  266.  */
  267.  
  268. void
  269. window (image, nRow, nCol, nRowS, nColS)
  270.      float **image;             /* image to be smoothed */
  271.      long nRow, nCol;           /* size of image to be windowed */
  272.      long nRowS, nColS;         /* size of smaller image b4 padding */
  273. {
  274.   float *hammingX, *hammingY;   /* horiz. and vert. hamming windows */
  275.   double nM1;                   /* no. rows/cols minus 1 */
  276.   long xStart, yStart;          /* offset smaller image in larger */
  277.   long x, y;
  278.  
  279. /* allocate space for Hamming windows */
  280.   if ((hammingX = (float *) calloc (nCol, sizeof (*hammingX))) == NULL)
  281.     exit (1);
  282.   if ((hammingY = (float *) calloc (nRow, sizeof (*hammingY))) == NULL)
  283.     exit (2);
  284.  
  285. /* put window in middle of zeroed vector if image has been padded */
  286.   for (x = 0; x < nCol; x++)
  287.     hammingX[x] = (float) 0.0;
  288.   for (y = 0; y < nRow; y++)
  289.     hammingY[y] = (float) 0.0;
  290.   xStart = (nCol - nColS) / 2;
  291.   yStart = (nRow - nRowS) / 2;
  292.  
  293. /* calculate Hamming window */
  294.   nM1 = (double) (nColS - 1);
  295.   for (x = 0; x < nColS; x++)
  296.     hammingX[x + xStart] = (float) (0.54 - 0.46 * cos ((M_PI * 2.0 * (double) x) / nM1));
  297.   nM1 = (double) (nRowS - 1);
  298.   for (y = 0; y < nRowS; y++)
  299.     hammingY[y + yStart] = (float) (0.54 - 0.46 * cos ((M_PI * 2.0 * (double) y) / nM1));
  300.  
  301. /* apply window to image */
  302.   for (y = 0; y < nRow; y++)
  303.     for (x = 0; x < nCol; x++)
  304.       image[y][x] = image[y][x] * hammingX[x] * hammingY[y];
  305.  
  306. }
  307.  
  308.  
  309. /* USAGE:       function gives instructions on usage of program
  310.  *                      usage: usage (flag)
  311.  *              When flag is 1, the long message is given, 0 gives short.
  312.  */
  313.  
  314. void
  315. usage (flag)
  316.      short flag;                /* flag =1 for long message; =0 for short message */
  317. {
  318.  
  319. /* print short usage message or long */
  320.   printf ("USAGE: powspec inimg outimg [-l] [-p] [-w] [-c] [-L]\n");
  321.  
  322.   if (flag == 0)
  323.     exit (1);
  324.  
  325.   printf ("\npowspec performs 2-D FFT to produce\n");
  326.   printf ("image of power spectrum.\n");
  327.   printf ("ARGUMENTS:\n");
  328.   printf ("    inimg: input image filename (TIF)\n");
  329.   printf ("   outimg: output image filename (TIF)\n\n");
  330.   printf ("OPTIONS:\n");
  331.   printf ("    -l: if set, displays linear scale for power in\n");
  332.   printf ("        power spectrum, instead of logarithmic default.\n");
  333.   printf ("    -p: if set, and image row or column is\n");
  334.   printf ("        not a power of 2, image size is increased by zero-.\n");
  335.   printf ("        padding; otherwise image size is decreased (default).\n");
  336.   printf ("    -w: if set, no smoothing window is applied \n");
  337.   printf ("        to initial image; default is to apply window.\n");
  338.   printf ("    -c: display 2 plots: top plot is summation of power over all angles\n");
  339.   printf ("        shown as a function of frequency or distance from center;\n");
  340.   printf ("        bottom plot is summation of power over all frequencies\n");
  341.   printf ("        shown as a function of angle, where 0 degrees is in middle.\n");
  342.   printf ("    -L: print Software License for this module\n");
  343.  
  344.   exit (1);
  345. }
  346.  
  347.  
  348. /* INPUT:       function reads input parameters
  349.  *                   usage: input (argc, argv, &linear, &zeroPad, &noWndw)
  350.  */
  351.  
  352. void
  353. input (argc, argv, linear, zeroPad, noWndw, crossSec)
  354.      int argc;
  355.      char *argv[];
  356.      short *linear;             /* if 1, linear power display; 0 =>for log display */
  357.      short *zeroPad;            /* if 1, increase img size to 2-power; 0 => decrease */
  358.      short *noWndw;             /* if 1, do not use window; 0 => smoothing window */
  359.      short *crossSec;           /* if 1, display cross section; 0 otherwise */
  360. {
  361.   long n;
  362.   double atof ();
  363.  
  364.   if (argc < 3)
  365.     usage (1);
  366.  
  367.   *linear = 0;
  368.   *zeroPad = 0;
  369.   *noWndw = 0;
  370.   *crossSec = 0;
  371.  
  372.   for (n = 3; n < argc; n++) {
  373.     if (strcmp (argv[n], "-l") == 0)
  374.       *linear = 1;
  375.     else if (strcmp (argv[n], "-p") == 0)
  376.       *zeroPad = 1;
  377.     else if (strcmp (argv[n], "-w") == 0)
  378.       *noWndw = 1;
  379.     else if (strcmp (argv[n], "-c") == 0)
  380.       *crossSec = 1;
  381.     else if (strcmp (argv[n], "-L") == 0) {
  382.       print_sos_lic ();
  383.       exit (0);
  384.     }
  385.     else
  386.       usage (0);
  387.   }
  388.  
  389.   return;
  390. }
  391.