home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_4.10 / HOUGH / HOUGH.C next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  7.0 KB  |  237 lines
  1. /* 
  2.  * hough.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /* HOUGH:       program performs Hough transform on binary image
  10.  *                    and yields an image of the transform space as well
  11.  *                      as the line corresponding to the highest peak
  12.  */
  13.  
  14. #include <stdio.h>
  15. #include <string.h>
  16. #include <stdlib.h>
  17. #include <tiffimage.h>          /* tiff info on images */
  18. #include <images.h>
  19. #include <math.h>
  20. extern void print_sos_lic ();
  21.  
  22. #define ON 0
  23. #define DFLT_BORDER 2           /* default border of image */
  24. #if defined(WIN32)
  25. #define PI 3.1415926536
  26. #endif
  27. #define PID2 1.5707963268
  28. #define TAN_TOO_BIG 263.0       /* tan so big that angle is PI/2 */
  29. #define INBOUNDX(X) ((X) >= 0 && (X) < width)
  30. #define INBOUNDY(Y) ((Y) >= 0 && (Y) < height)
  31.  
  32. int usage (short);
  33. int input (int, char **, long *, short *);
  34.  
  35. main (argc, argv)
  36.      int argc;
  37.      char *argv[];
  38. {
  39.   Image *imgI, *imgH;           /* I/O image structures */
  40.   unsigned char **imgIn, **hough;  /* input/output images */
  41.   long width, height;           /* input image size */
  42.   long border;                  /* zero out border of image */
  43.   short dFlag;                  /* if = 0, display peak in Hough space */
  44.   long xEnd, yEnd;              /* end of image within borders */
  45.   long thetaHt, rhoWid;         /* width and height of Hough space */
  46.   long rhoWidM1;                /* rho width minus 1 */
  47.   double rho, theta;            /* radius and angle in Hough space */
  48.   double tanTheta;              /* tan of theta angle */
  49.   double denom;                 /* denominator */
  50.   double rhoNorm;               /* normalization factor for rho */
  51.   long max, xMax, yMax;         /* peak point in Hough space */
  52.   double x1, y1;
  53.   long i, j;
  54.   long x, y;
  55.  
  56.   if ((input (argc, argv, &border, &dFlag)) < 0)
  57.     return (-1);
  58.  
  59. /* open input and output images */
  60.   imgI = ImageIn (argv[1]);
  61.   imgIn = imgI->img;
  62.   height = ImageGetHeight (imgI);
  63.   width = ImageGetWidth (imgI);
  64.   printf ("image size is %dx%d.\n", width, height);
  65.  
  66.   rhoWid = width;
  67.   thetaHt = height;
  68.   rhoWidM1 = rhoWid - 1;
  69.   imgH = ImageAlloc (thetaHt, rhoWid, 8);
  70.   hough = ImageGetPtr (imgH);
  71.  
  72. /* initialize accumulator bins */
  73.   for (i = 0; i < thetaHt; i++)
  74.     for (j = 0; j < rhoWid; j++)
  75.       hough[i][j] = 0;
  76.  
  77. /* perform Hough transform */
  78.  
  79.   rhoNorm = rhoWidM1 / sqrt ((double) (width * width) + (double) height * height);
  80.  
  81.   yEnd = height - border;
  82.   xEnd = width - border;
  83.   for (y = border + 1; y < yEnd; y++) {
  84.     for (x = border + 1; x < xEnd; x++) {
  85.       if (imgIn[y][x] == ON) {
  86.         for (i = 0; i < thetaHt; i++) {
  87.           theta = (double) i *PI / (double) thetaHt;
  88.           tanTheta = tan (theta);
  89.           if (tanTheta > TAN_TOO_BIG) {
  90.             /*printf ("too big: x = %d, i = %d, j = %d\n", x, i, j);*/
  91.             rho = (double) x;
  92.           }
  93.           else {
  94.             denom = tanTheta * tanTheta + 1.0;
  95.             y1 = ((double) y - (double) x * tanTheta) / denom;
  96.             x1 = ((double) x * tanTheta * tanTheta - (double) y * tanTheta) / denom;
  97.             rho = sqrt (x1 * x1 + y1 * y1);
  98.           }
  99.           j = (long) (rho * rhoNorm + 0.5);
  100.           if (hough[i][j] < 255)
  101.             hough[i][j]++;
  102.         }
  103.       }
  104.     }
  105.   }
  106.  
  107. /* determine peak in Hough space and (theta, rho) coords */
  108.   max = 0;
  109.   for (i = 0; i < thetaHt; i++) {
  110.     for (j = 0; j < rhoWid; j++) {
  111.       if (hough[i][j] > max) {
  112.         max = hough[i][j];
  113.         xMax = j;
  114.         yMax = i;
  115.       }
  116.     }
  117.   }
  118.   if (max > 0) {
  119.     rho = xMax / rhoNorm;
  120.     theta = yMax * PI / (double) thetaHt;
  121.     printf ("Line located with parameters:\n");
  122.     printf ("\tpolar coords: (%5.2f, %5.2f) (rho [pixels], theta [radians])\n",
  123.             rho, theta);
  124.     if (dFlag) {                /* display peak if flag set */
  125.       printf ("Peak located at: (%d,%d)\n", xMax, yMax);
  126.       for (i = yMax - 3; i <= yMax + 3; i++) {
  127.         if (INBOUNDY (i)) {
  128.           if (INBOUNDX (xMax - 3))
  129.             hough[i][xMax - 3] = (hough[i][xMax - 3] > 128) ? 0 : 255;
  130.           if (INBOUNDX (xMax + 4))
  131.             hough[i][xMax + 3] = (hough[i][xMax + 3] > 128) ? 0 : 255;
  132.         }
  133.       }
  134.       for (j = xMax - 3; j <= xMax + 3; j++) {
  135.         if (INBOUNDX (j)) {
  136.           if (INBOUNDY (yMax - 3))
  137.             hough[yMax - 3][j] = (hough[yMax - 3][j] > 128) ? 0 : 255;
  138.           if (INBOUNDY (yMax + 4))
  139.             hough[yMax + 3][j] = (hough[yMax + 3][j] > 128) ? 0 : 255;
  140.         }
  141.       }
  142.     }
  143.   }
  144.   else
  145.     printf ("No line found.\n");
  146.  
  147. /* invert Hough image for display */
  148.   for (i = 0; i < thetaHt; i++)
  149.     for (j = 0; j < rhoWid; j++)
  150.       hough[i][j] = 255 - hough[i][j];
  151.  
  152.   ImageOut (argv[2], imgH);
  153.  
  154.   return (0);
  155. }
  156.  
  157.  
  158.  
  159. /* USAGE:       function gives instructions on usage of program
  160.  *                    usage: usage (flag)
  161.  *              When flag is 1, the long message is given, 0 gives short.
  162.  */
  163.  
  164. int
  165. usage (flag)
  166.      short flag;                /* flag =1 for long message; =0 for short message */
  167. {
  168.  
  169. /* print short usage message or long */
  170.   printf ("USAGE: hough inimg outimg [-b BORDER] [-d] [-L]\n");
  171.   if (flag == 0)
  172.     return (-1);
  173.  
  174.   printf ("\nhough transforms binary image from\n");
  175.   printf ("spatial domain ((x,y) coordinates) to polar coordinate\n");
  176.   printf ("domain ((rho, theta) coordinates); peak in the polar\n");
  177.   printf ("(\"Hough\")domain indicates a dominant line in the spatial domain;\n");
  178.   printf ("NOTE: origin (0,0) of Hough transform\n");
  179.   printf ("      image is in top-left corner;\n");
  180.   printf ("      rho increases along horizontal axis to\n");
  181.   printf ("      maximum rho equal to image diagonal length;\n");
  182.   printf ("      theta increases downward from 0 radians to PI radians.\n\n");
  183.   printf ("ARGUMENTS:\n");
  184.   printf ("    inimg: input image filename (TIF)\n");
  185.   printf ("   outimg: output image filename (TIF)\n\n");
  186.   printf ("OPTIONS:\n");
  187.   printf ("  -b BORDER: to remove noise at borders of image by omitting\n");
  188.   printf ("             image for this number of rows/cols [default = %d].\n", DFLT_BORDER);
  189.   printf ("         -d: display peak in Hough transform space.\n");
  190.   printf ("         -L: print Software License for this module\n");
  191.  
  192.   return (-1);
  193. }
  194.  
  195.  
  196. /* INPUT:       function reads input parameters
  197.  *                  usage: input (argc, argv, &dFlag)
  198.  */
  199.  
  200. #define USAGE_EXIT(VALUE) {usage (VALUE); return (-1);}
  201.  
  202. int
  203. input (argc, argv, border, dFlag)
  204.      int argc;
  205.      char *argv[];
  206.      long *border;              /* zero out border of image */
  207.      short *dFlag;              /* if = 1, display peak */
  208. {
  209.   long n;
  210.  
  211.   if (argc == 1)
  212.     USAGE_EXIT (1);
  213.   if (argc == 2)
  214.     USAGE_EXIT (0);
  215.  
  216.   *dFlag = 0;
  217.   *border = DFLT_BORDER;
  218.  
  219.   for (n = 3; n < argc; n++) {
  220.     if (strcmp (argv[n], "-b") == 0) {
  221.       if (++n == argc || argv[n][0] == '-')
  222.         USAGE_EXIT (0);
  223.       *border = atol (argv[n]);
  224.     }
  225.     else if (strcmp (argv[n], "-d") == 0)
  226.       *dFlag = 1;
  227.     else if (strcmp (argv[n], "-L") == 0) {
  228.       print_sos_lic ();
  229.       exit (0);
  230.     }
  231.     else
  232.       USAGE_EXIT (0);
  233.   }
  234.  
  235.   return (0);
  236. }
  237.