home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_6.4 / XKNN / XKNN.C < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  8.6 KB  |  378 lines
  1. /* 
  2.  * xknn.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * XKNN(K-Nearest Neighbors) 
  11.  *
  12.  * test knn.c module
  13.  * 
  14.  */
  15.  
  16. #include "xknn.h"
  17.  
  18. #define BUFSZ 1024
  19. #define     ON        1
  20. #define     OFF        0
  21. extern short tiffInput;         /* flag=0 if no ImageIn to set tags; else =1 */
  22.  
  23. /*
  24.  * Globals
  25.  */
  26.  
  27. int maxX;
  28. int maxY;
  29. int minX;
  30. int minY;
  31. int npoints;
  32. extern char *optarg;
  33. extern int optind, opterr;
  34.  
  35. /*
  36.  * usage of routine
  37.  */
  38. void
  39. usage (char *progname)
  40. {
  41.   progname = last_bs (progname);
  42.   printf ("USAGE: %s infile knn outimg [-l][-L]\n", progname);
  43.   printf ("\n%s finds K-nearest neighbors for each point\n", progname);
  44.   printf ("in a list of x-y points\n\n");
  45.   printf ("ARGUMENTS:\n");
  46.   printf ("   infile: input file with list of (x,y) pairs (ASCII)\n");
  47.   printf ("      knn: number of nearest neighbors to calculate (int)\n");
  48.   printf ("   outimg: output image filename (TIF)\n\n");
  49.   printf ("OPTIONS:\n");
  50.   printf ("       -l: prints (x,y) point with neighbor list\n");
  51.   printf ("       -L: print Software License for this module\n");
  52.   exit (1);
  53. }
  54.  
  55. struct Point *
  56. readPoints (char *filename, int nknn)
  57. {
  58.   FILE *stream;
  59.   struct Point *pointStruct, *p;
  60.   int line_no, i, got_first_point;
  61.   char inBuf[BUFSZ];
  62.  
  63.   line_no = 0;
  64.   got_first_point = 0;
  65.   if ((stream = fopen (filename, "r")) != NULL) {
  66.     for (;;) {
  67.       line_no++;
  68.       if (fgets (inBuf, BUFSZ, stream) == (char *) NULL) {
  69.         printf ("Premature EOF encountered on line %d while reading points file %s\n", line_no, filename);
  70.         exit (1);
  71.       }
  72.       /*
  73.        * Throw away comments
  74.        */
  75.       if (inBuf[0] == '#')
  76.         continue;
  77.       else {
  78.         if (sscanf (inBuf, "%d \n", &npoints) != 1) {
  79.           printf ("error getting number of points\n");
  80.           exit (1);
  81.         }
  82.         if ((pointStruct = (struct Point *) calloc (npoints, sizeof (struct Point))) == NULL) {
  83.           printf ("\n...mem allocation for retMatrix failed\n");
  84.           exit (1);
  85.         }
  86.         break;
  87.       }
  88.     }
  89. /*
  90.  * Now read the (x,y) points in
  91.  */
  92.     i = 0;
  93.     p = pointStruct;
  94.     for (;;) {
  95.       line_no++;
  96.       if (fgets (inBuf, BUFSZ, stream) == (char *) NULL) {
  97.         printf ("Premature EOF encountered on line %d while reading points file %s\n", line_no, filename);
  98.         exit (1);
  99.       }
  100.       /*
  101.        * Throw away comments
  102.        */
  103.       if (inBuf[0] == '#')
  104.         continue;
  105.       if (sscanf (inBuf, "%d %d\n", &(p->x_pos), &(p->y_pos)) == 2) {
  106.         if (!got_first_point) {
  107.           maxX = p->x_pos;
  108.           minX = p->x_pos;
  109.           maxY = p->y_pos;
  110.           minY = p->y_pos;
  111.           got_first_point = 1;
  112.         }
  113.         else {
  114.           maxX = (p->x_pos > maxX) ? p->x_pos : maxX;
  115.           minX = (p->x_pos < minX) ? p->x_pos : minX;
  116.           maxY = (p->y_pos > maxY) ? p->y_pos : maxY;
  117.           minY = (p->y_pos < minY) ? p->y_pos : minY;
  118.         }
  119.         if ((p->Knn = (struct KNN *) calloc (nknn, sizeof (struct KNN))) == NULL) {
  120.           printf ("\n...mem allocation for Knn struct failed\n");
  121.           exit (1);
  122.         }
  123.         p++;
  124.         i++;
  125.         if (i >= npoints)
  126.           break;
  127.       }
  128.       else {
  129.         printf ("error getting (x,y) value on line %d of file %s\n", line_no, filename);
  130.         fclose (stream);
  131.         exit (1);
  132.       }
  133.     }
  134.     fclose (stream);
  135.     return (pointStruct);
  136.   }
  137.   else {
  138.     printf ("Can not open points file: %s\n", filename);
  139.     exit (1);
  140.   }
  141. }
  142.  
  143. /*
  144.  * qsort comparison function for x-direction
  145.  */
  146.  
  147. int
  148. sortX (p1, p2)
  149.      struct Point *p1;
  150.      struct Point *p2;
  151. {
  152.   if (p1->x_pos == p2->x_pos)
  153.     return (0);
  154.   else
  155.     return ((p1->x_pos < p2->x_pos) ? -1 : 1);
  156. }
  157.  
  158. /*
  159.  * qsort comparison function for y-direction
  160.  */
  161.  
  162. int
  163. sortY (p1, p2)
  164.      struct Point *p1;
  165.      struct Point *p2;
  166. {
  167.   if (p1->y_pos == p2->y_pos)
  168.     return (0);
  169.   else
  170.     return ((p1->y_pos < p2->y_pos) ? -1 : 1);
  171. }
  172.  
  173. int
  174. dump_point (struct Point *p1, int knn)
  175. {
  176.   struct KNN *pk;
  177.   int i;
  178.  
  179.   pk = p1->Knn;
  180.   printf ("struct for Point (%d, %d), nNN = %d\n", p1->x_pos, p1->y_pos, p1->nNN);
  181.   for (i = 0; i < knn; i++) {
  182.     printf ("pk[%d] = (%d, %d), distSqd = %d\n", i, pk[i].x_pos, pk[i].y_pos, pk[i].distSqd);
  183.   }
  184.   return (0);
  185. }
  186.  
  187. /*
  188.  * find_knn()
  189.  * finds k-nearest neighbors for point array
  190.  * if maxRange = 0, array is sorted in x-direction, otherwise y-direction
  191.  */
  192.  
  193. void
  194. find_knn (struct Point *points, int npoints, int knn, int maxRange)
  195. {
  196.   int i, j;
  197.   long xDist, yDist, distSqd;
  198.   struct Point *p1, *p2;
  199.  
  200.   /* outer loop through all points */
  201.   for (i = 0, p1 = points; i < npoints; i++, p1++) {
  202.     for (j = i + 1, p2 = p1 + 1; j < npoints; j++, p2++) {
  203.       /*for(j = 0, p2 = points; j < npoints; j++, p2++) { */
  204.       if (p2 == p1)
  205.         break;
  206.       xDist = p2->x_pos - p1->x_pos;
  207.       yDist = p2->y_pos - p1->y_pos;
  208.       distSqd = xDist * xDist + yDist * yDist;
  209.       if (!maxRange && ((xDist * xDist) > p1->Knn[knn].distSqd))
  210.         break;
  211.       if (maxRange && (p1->Knn[knn - 1].distSqd > 0) && ((yDist * yDist) > (p1->Knn[knn - 1].distSqd))) {
  212.         //dump_point(p1, knn);
  213.         break;
  214.       }
  215.       insertNN (p1, p2, distSqd, knn);
  216.       insertNN (p2, p1, distSqd, knn);
  217.     }
  218.   }
  219. }
  220.  
  221. /*
  222.  * insertNN()
  223.  * p2 is inserted into p1's knn list if the
  224.  * distance is less than any of the distances
  225.  * currently in p1's knn list
  226.  * RETURN VALUE:
  227.  *  1 if p2 was inserted
  228.  *  0 otherwise
  229.  */
  230.  
  231. int
  232. insertNN (struct Point *p1, struct Point *p2, long distSqd, int knn)
  233. {
  234.   int i, j;
  235.   int nnInsert = 0;
  236.   struct KNN *pk;
  237.  
  238.   pk = p1->Knn;
  239.   for (i = 0; i < knn; i++) {
  240.     if (distSqd < pk[i].distSqd) {
  241.       for (j = knn - 1; j > i; --j)
  242.         pk[j] = pk[j - 1];
  243.       pk[i].x_pos = p2->x_pos;
  244.       pk[i].y_pos = p2->y_pos;
  245.       pk[i].distSqd = distSqd;
  246.       if (p1->nNN < knn)
  247.         p1->nNN++;
  248.       nnInsert = 1;
  249.       break;
  250.     }
  251.   }
  252.   if (!nnInsert && p1->nNN < knn) {
  253.     pk[p1->nNN].x_pos = p2->x_pos;
  254.     pk[p1->nNN].y_pos = p2->y_pos;
  255.     pk[p1->nNN].distSqd = distSqd;
  256.     p1->nNN++;
  257.     nnInsert = 1;
  258.   }
  259.   //dump_point(p1, knn);
  260.   return (nnInsert);
  261. }
  262.  
  263.  
  264. /*
  265.  * main() for xsgll
  266.  */
  267.  
  268. int
  269. main (int argc, char *argv[])
  270. {
  271.  
  272.   Image *imgOut;
  273.   struct Point *points;
  274.   struct Point *p;
  275.   int rangeX;
  276.   int rangeY;
  277.   int maxRange;
  278.   int knn;
  279.   int i, j;
  280.   int i_arg;
  281.   int print_list = 0;
  282.  
  283. /*
  284.  * cmd line options
  285.  */
  286.   static char *optstring = "lL";
  287.  
  288.  
  289. /*
  290.  * parse command line
  291.  */
  292.   optind = 4;                   /* set getopt to point to the 4th arg */
  293.   opterr = ON;                  /* give error messages */
  294.  
  295.  
  296.   if (argc < 4)
  297.     usage (argv[0]);
  298.  
  299.   while ((i_arg = getopt (argc, argv, optstring)) != EOF) {
  300.     switch (i_arg) {
  301.     case 'l':
  302.       printf ("\n...option %c: print list of points with K-nn\n", i_arg);
  303.       print_list = 1;
  304.       break;
  305.     case 'L':
  306.       print_sos_lic ();
  307.       exit (0);
  308.     default:
  309.       printf ("\ngetopt: unknown condition encountered\n");
  310.       exit (1);
  311.       break;
  312.     }
  313.   }
  314.  
  315.   if (sscanf (argv[2], "%d", &knn) != 1) {
  316.     printf ("Error getting value for nearnbr\n");
  317.     usage (argv[0]);
  318.   }
  319.  
  320.   /* reset tiffInput so that we write a grayscale file (i.e tags are not copied) */
  321.   tiffInput = 0;
  322.  
  323. /*
  324.  * Read in points
  325.  */
  326.   points = readPoints (argv[1], knn);
  327.   /* determine the maximum range of values, whether over x or y */
  328.   rangeX = maxX - minX;
  329.   rangeY = maxY - minY;
  330.   maxRange = (rangeX > rangeY) ? 0 : 1;  /* 0 or 1 for max X or Y range */
  331. /* order elements along x or y axes, whichever is maxRange */
  332.   if (maxRange == 0)
  333. #if defined(LINUX)
  334.     qsort (points, (size_t) npoints, sizeof (struct Point), (__compar_fn_t) sortX);
  335.   else
  336.     qsort (points, (size_t) npoints, sizeof (struct Point), (__compar_fn_t) sortY);
  337. #else
  338.     qsort (points, (size_t) npoints, sizeof (struct Point), sortX);
  339.   else
  340.     qsort (points, (size_t) npoints, sizeof (struct Point), sortY);
  341. #endif
  342.  
  343. /*
  344.  * Allocate memory for output image
  345.  */
  346.   imgOut = ImageAlloc ((long) maxY + 10, (long) maxX + 10, (long) BPS);
  347.  
  348.   if (print_list) {
  349.     printf ("Sorted array is:\n");
  350.     for (i = 0, p = points; i < npoints; i++, p++)
  351.       printf ("Point %d = (%d, %d)\n", i, p->x_pos, p->y_pos);
  352.   }
  353.  
  354. /*
  355.  * Find the nearest neighbors for all points
  356.  */
  357.   find_knn (points, npoints, knn, maxRange);
  358.  
  359. /*
  360.  * Plot the lines
  361.  */
  362.   for (i = 0, p = points; i < npoints; i++, p++) {
  363.     draw_cross (p->x_pos, p->y_pos, 10, imgOut, WHITE);
  364.     if (print_list)
  365.       printf ("Nearest neighbors for Point %d = (%d, %d), %d\n", i, p->x_pos, p->y_pos, p->nNN);
  366.     for (j = 0; j < p->nNN; j++) {
  367.       if (print_list)
  368.         printf ("\t(%d, %d)\n", p->Knn[j].x_pos, p->Knn[j].y_pos);
  369.       draw_line (p->x_pos, p->y_pos, p->Knn[j].x_pos, p->Knn[j].y_pos, imgOut, WHITE);
  370.     }
  371.   }
  372. /*
  373.  * Write the output image
  374.  */
  375.   ImageOut (argv[3], imgOut);
  376.   return (1);
  377. }
  378.