home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_6.4 / XKNN / xknn.c.orig < prev    next >
Encoding:
Text File  |  1999-09-11  |  7.7 KB  |  317 lines
  1. /* 
  2. ** xknn.c
  3. ** 
  4. ** Practical Algorithms for Image Analysis
  5. ** 
  6. ** Copyright (c) 1997 ????
  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. extern short tiffInput;        /* flag=0 if no ImageIn to set tags; else =1 */
  20.  
  21. /*
  22.  * Globals
  23.  */
  24.  
  25. int maxX;
  26. int maxY;
  27. int minX;
  28. int minY;
  29. int npoints;
  30.  
  31. /*
  32.  * usage of routine
  33.  */
  34. void 
  35. usage (char *progname)
  36. {
  37.   progname = last_bs (progname);
  38.   printf ("\nusage:\n");
  39.   printf ("%s nearnbr input_points output_image [-l]\n", progname);
  40.   printf ("   nearnbr - number of nearest neighbors to calculate (int)\n");
  41.   printf ("   input_points - input file with list of (x,y) pairs\n");
  42.   printf ("   output_image - output image filename (TIFF format)\n");
  43.   printf ("\noptions:\n");
  44.   printf ("   -l prints (x,y) point with neighbor list");
  45.   exit (1);
  46. }
  47.  
  48. struct Point *
  49. readPoints(char *filename, int nknn)
  50. {
  51.     FILE *stream;
  52.     struct Point *pointStruct, *p;
  53.     int line_no, i, got_first_point;
  54.     char inBuf[BUFSZ];
  55.  
  56.     line_no = 0;
  57.     got_first_point = 0;
  58.     if( (stream = fopen( filename, "r" )) != NULL ) {
  59.         for(;;) {
  60.             line_no++;
  61.             if(fgets(inBuf, BUFSZ, stream) == (char *)NULL) {
  62.                 printf("Premature EOF encountered on line %d while reading points file %s\n", line_no, filename);
  63.                 exit(1);
  64.             }
  65.             /*
  66.              * Throw away comments
  67.              */
  68.             if(inBuf[0] == '#')
  69.                 continue;
  70.             else {
  71.                 if( sscanf( inBuf, "%d \n", &npoints ) != 1) {
  72.                     printf( "error getting number of points\n" );
  73.                     exit(1);
  74.                 }
  75.                 if( (pointStruct = (struct Point *)calloc(npoints, sizeof(struct Point)) ) == NULL) {
  76.                     printf("\n...mem allocation for retMatrix failed\n");
  77.                     exit(1);
  78.                 }
  79.                 break;
  80.             }
  81.         }
  82. /*
  83. ** Now read the (x,y) points in
  84. */
  85.         i =0;
  86.         p = pointStruct;
  87.         for(;;) {
  88.             line_no++;
  89.             if(fgets(inBuf, BUFSZ, stream) == (char *)NULL) {
  90.                 printf("Premature EOF encountered on line %d while reading points file %s\n", line_no, filename);
  91.                 exit(1);
  92.             }
  93.             /*
  94.              * Throw away comments
  95.              */
  96.             if(inBuf[0] == '#')
  97.                 continue;
  98.             if( sscanf( inBuf, "%d %d\n", &(p->x_pos), &(p->y_pos) ) == 2) {
  99.                 if(!got_first_point) {
  100.                     maxX = p->x_pos;
  101.                     minX = p->x_pos;
  102.                     maxY = p->y_pos;
  103.                     minY = p->y_pos;
  104.                     got_first_point = 1;
  105.                 }
  106.                 else {
  107.                     maxX = (p->x_pos > maxX) ? p->x_pos : maxX;
  108.                     minX = (p->x_pos < minX) ? p->x_pos : minX;
  109.                     maxY = (p->y_pos > maxY) ? p->y_pos : maxY;
  110.                     minY = (p->y_pos < minY) ? p->y_pos : minY;
  111.                 }
  112.                 if( (p->Knn = (struct KNN *)calloc(nknn, sizeof(struct KNN)) ) == NULL) {
  113.                     printf("\n...mem allocation for Knn struct failed\n");
  114.                     exit(1);
  115.                 }
  116.                 p++;
  117.                 i++;
  118.                 if(i >= npoints)
  119.                     break;
  120.             }
  121.             else {
  122.                 printf( "error getting (x,y) value on line %d of file %s\n", line_no, filename);
  123.                 fclose( stream );
  124.                 exit(1);
  125.             }
  126.         }
  127.         fclose( stream );
  128.         return(pointStruct); 
  129.    }
  130.    else {
  131.        printf("Can not open points file: %s\n", filename);
  132.        exit(1);
  133.    }
  134. }
  135.  
  136. /*
  137.  * qsort comparison function for x-direction
  138.  */
  139.  
  140. int 
  141. sortX(p1, p2)
  142. struct Point *p1;
  143. struct Point *p2;
  144. {
  145.     if(p1->x_pos == p2->x_pos)
  146.         return(0);
  147.     else
  148.         return((p1->x_pos < p2->x_pos) ? -1 : 1);
  149. }
  150.  
  151. /*
  152.  * qsort comparison function for y-direction
  153.  */
  154.  
  155. int
  156. sortY(p1, p2)
  157. struct Point *p1;
  158. struct Point *p2;
  159. {
  160.     if(p1->y_pos == p2->y_pos)
  161.         return(0);
  162.     else
  163.         return((p1->y_pos < p2->y_pos) ? -1 : 1);
  164. }
  165.  
  166. /*
  167.  * find_knn()
  168.  * finds k-nearest neighbors for point array
  169.  * if maxRange = 0, array is sorted in x-direction, otherwise y-direction
  170.  */
  171.  
  172. void
  173. find_knn(struct Point *points, int npoints, int knn, int maxRange)
  174. {
  175.     int i, j;
  176.     long xDist, yDist, distSqd;
  177.     struct Point *p1, *p2;
  178.  
  179.     /* outer loop through all points */
  180.     for(i = 0, p1 = points; i < npoints; i++, p1++) {
  181.         for(j = i + 1, p2 = p1 + 1; j < npoints; j++, p2++) {
  182.             xDist = p2->x_pos - p1->x_pos;
  183.             yDist = p2->y_pos - p1->y_pos;
  184.             distSqd = xDist * xDist + yDist * yDist;
  185.             if(!maxRange && (xDist * xDist > p1->Knn[knn].distSqd))
  186.                 break;
  187.             if(maxRange && (yDist * yDist > p1->Knn[knn].distSqd))
  188.                 break;
  189.             insertNN(p1, p2, distSqd, knn);
  190.             insertNN(p2, p1, distSqd, knn);
  191.         }
  192.     }
  193. }
  194.  
  195. /*
  196.  * insertNN()
  197.  * p2 is inserted into p1's knn list if the
  198.  * distance is less than any of the distances
  199.  * currently in p1's knn list
  200.  * RETURN VALUE:
  201.  *  1 if p2 was inserted
  202.  *  0 otherwise
  203.  */
  204.  
  205. int
  206. insertNN(struct Point *p1, struct Point *p2, long distSqd, int knn)
  207. {
  208.     int i, j;
  209.     int nnInsert = 0;
  210.     struct KNN *pk;
  211.  
  212.     pk = p1->Knn;
  213.     for(i = 0; i < knn; i++) {
  214.         if(distSqd < pk[i].distSqd) {
  215.             for(j = knn - 1; j > i; --j)
  216.                 pk[j] = pk[j-1];
  217.             pk[i].x_pos = p2->x_pos;
  218.             pk[i].y_pos = p2->y_pos;
  219.             pk[i].distSqd = distSqd;
  220.             if(p1->nNN < knn)
  221.                 p1->nNN++;
  222.             nnInsert = 1;
  223.             break;
  224.         }
  225.     }
  226.     if(!nnInsert && p1->nNN < knn) {
  227.         pk[p1->nNN].x_pos = p2->x_pos;
  228.         pk[p1->nNN].y_pos = p2->y_pos;
  229.         pk[p1->nNN].distSqd = distSqd;
  230.         p1->nNN++;
  231.         nnInsert = 1;
  232.     }
  233.     return(nnInsert);
  234. }
  235.  
  236.  
  237. /*
  238.  * main() for xsgll
  239.  */
  240.  
  241. int
  242. main(int argc, char *argv[])
  243. {
  244.  
  245.     Image *imgOut;
  246.     struct Point *points;
  247.     struct Point *p;
  248.     int rangeX;
  249.     int rangeY;
  250.     int maxRange;
  251.     int knn;
  252.     int i, j;
  253.  
  254.     if( argc < 4 )        
  255.         usage(argv[0]);
  256.  
  257.     if( sscanf( argv[1], "%d", &knn) != 1) {
  258.         printf("Error getting value for nearnbr\n");
  259.         usage(argv[0]);
  260.     }
  261.  
  262.     /* reset tiffInput so that we write a grayscale file (i.e tags are not copied) */
  263.     tiffInput = 0;
  264.  
  265. /*
  266. ** Read in points
  267. */
  268.     points = readPoints(argv[2], knn);
  269.     /* determine the maximum range of values, whether over x or y */
  270.     rangeX = maxX - minX;
  271.     rangeY = maxY - minY;
  272.     maxRange = (rangeX > rangeY) ? 0: 1;  /* 0 or 1 for max X or Y range */
  273. /* order elements along x or y axes, whichever is maxRange */
  274.     if (maxRange == 0) 
  275. #if defined(LINUX)
  276.         qsort(points, (size_t)npoints, sizeof(struct Point), (__compar_fn_t)sortX);
  277.     else
  278.         qsort(points, (size_t)npoints, sizeof(struct Point), (__compar_fn_t)sortY);
  279. #else
  280.         qsort(points, (size_t)npoints, sizeof(struct Point), sortX);
  281.     else
  282.         qsort(points, (size_t)npoints, sizeof(struct Point), sortY);
  283. #endif
  284.  
  285. /*
  286. ** Allocate memory for output image
  287. */
  288.     imgOut = ImageAlloc((long)maxY + 10, (long)maxX + 10, (long)BPS);
  289.  
  290.     printf("Sorted array is:\n");
  291.     for(i = 0, p = points; i < npoints; i++, p++)
  292.         printf("Point %d = (%d, %d)\n", i, p->x_pos, p->y_pos);
  293.  
  294. /*
  295. ** Find the nearest neighbors for all points
  296. */
  297.     find_knn(points, npoints, knn, maxRange);
  298.  
  299. /*
  300. ** Plot the lines
  301. */
  302.     for(i = 0, p = points; i < npoints; i++, p++) {
  303.         draw_cross(p->x_pos, p->y_pos, 10, imgOut, WHITE);
  304.         printf("Nearest neighbors for Point %d = (%d, %d)\n", i, p->x_pos, p->y_pos);
  305.         for(j = 0; j < p->nNN; j++) {
  306.             printf("\t(%d, %d)\n", p->Knn[j].x_pos, p->Knn[j].y_pos);
  307.             draw_line(p->x_pos, p->y_pos, p->Knn[j].x_pos, p->Knn[j].y_pos, imgOut, WHITE);
  308.         }
  309.     }
  310. /*
  311. ** Write the output image
  312. */
  313.     ImageOut (argv[3], imgOut);
  314.     return(1);
  315. }
  316.  
  317.