home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / CH_4.2 / KFILL / KFILL.C next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  26.3 KB  |  847 lines
  1. /* 
  2.  * kfill.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /* KFILL:       program performs kxk salt-and-pepper noise reduction
  10.  *            on binary image using table look-up and run-length coding
  11.  *              for speed.
  12.  *              usage: kfill inimg outimg [-f FDIFF] [-c] [-n NITER]
  13.  *                                        [-r R%][-k K] [-l] [-d] [-I] [-L]
  14.  *              For options, type "kfill".
  15.  *              NOTE: Because this program uses modified run-length code
  16.  *                    on only the ON pixels, the OFF and ON lengths are
  17.  *                    not treated equivalently. Specifically, while the
  18.  *                    calculations for the ON pixels are correct as 
  19.  *                    expected, those for the OFF pixels may not be. 
  20.  *                    For example, if an OFF run-length is longer than
  21.  *                    MIN0RUN, the pixels within that run are not inspected.
  22.  *                    If a scan line is omitted, then pixels on that line
  23.  *                    are never inspected, and they can never become ON.
  24.  *                    Oh, there is an exception to this on the first 
  25.  *                    iteration, since KFILL inspects all pixels to
  26.  *                    perform runlength coding on this iteration.
  27.  */
  28.  
  29. #include <stdio.h>
  30. #include <string.h>
  31. #include <stdlib.h>
  32. #include <math.h>
  33. #include "images.h"
  34. #include "tiffimage.h"
  35. extern void print_sos_lic ();
  36.  
  37. #define DFLTCFLAG 1             /* flag =1 if connectivity matters; or 0 */
  38. #define DFLTEFLAG 1             /* flag =1 if endpoint matters; or 0 */
  39. #define DFLTKMAX 3              /* dflt maximum thinning kernel size */
  40. #define MAXKMAX 21              /* max of the maximum thinning kernel size */
  41. #define DFLTITER 20             /* default number of iterations */
  42. #define DFLTDIRTY 0             /* default pct of max noise left on iter */
  43. #define DFLTDFLAG 0             /* flag =1 if display each iteration; or 0 */
  44.  
  45. #define MIN0RUN 5               /* minimum run of zeros to run-length code */
  46. #define FILLINITIAL 0           /* initial fill type */
  47. #define MAXDISPLAY 40           /* maximum length of line for display */
  48.  
  49. long OFF, ON;                   /* initial/final values of pixels */
  50. unsigned char FILL0, FILL1;     /* values of OFF fill, ON fill */
  51. long fDiff,                     /* fNum difference from default */
  52.   cFlag,                        /* if =1, retain connectivity; if 0, don't */
  53.   eFlag;                        /* if =1, retain endpoint; if 0, don't */
  54.  
  55. Image *imgIO;                   /* structure for I/O image */
  56. unsigned char **image;          /* input/output image */
  57. struct point imgSize;           /* image size */
  58. long ySizeM1, xSizeM1;          /* y,x lengths minus 1 */
  59. long **xRun;                    /* no., then x locns of 1/0 runs for each y */
  60.  
  61. int filltest (long, long, long);
  62. int fill0 (unsigned char **, unsigned char **, long, long,
  63.            long *, long *, long *);
  64. int fill (unsigned char **, unsigned char **, long, long,
  65.           long *, long *);
  66. long ksize (long, long, long, long);
  67. int getring (long, long, long, long, long *);
  68. int fillsqr (long, long, long, long, long *);
  69. int display (long);
  70. int stall (char *);
  71. long usage (short);
  72. long input (int, char **, long *, long *, long *, long *, long *,
  73.             long *, long *, long *, long *);
  74.  
  75. main (argc, argv)
  76.      int argc;
  77.      char *argv[];
  78. {
  79.   register long x, y,           /* image coordinates */
  80.     i, k;                       /* sidelength of current thinning kernel */
  81.   unsigned char **f0Table,      /* table of FILL0 values for each k */
  82.   **f1Table;                    /* table of FILL1 values for each k */
  83.   long kMax,                    /* max. sidelength of thinning kernel */
  84.     maxIter,                    /* maximum no. of iterations */
  85.     change[MAXKMAX],            /* no. erasures for each mask size */
  86.     nTable,                     /* size of fTable for each k */
  87.     nIter,                      /* no. iterations */
  88.     pctDirty,                   /* pct. of max noise left on an iteration */
  89.     nChange,                    /* no. thinning operations on iteration */
  90.     nChangeB4,                  /* no. changes on iteration before */
  91.     nChangeMax,                 /* no. max. changes in one iteration */
  92.     nChangeThresh,              /* change threshold of max pct removed */
  93.     displayFlag,                /* display results after each iter if =1 */
  94.     invertFlag,                 /* invert input image before processing */
  95.     nONs,                       /* total ONs in original image */
  96.     fillFlag,                   /* ON fill (if =1), or OFF erase (if =0) */
  97.     nFill;                      /* cumulative no. filled in image */
  98.   double pow ();
  99.  
  100.   if (input (argc, argv, &fDiff, &cFlag, &eFlag, &maxIter, &pctDirty, &kMax,
  101.              &fillFlag, &displayFlag, &invertFlag) < 0)
  102.     return (-1);
  103.  
  104.   imgIO = ImageIn (argv[1]);
  105.   image = imgIO->img;
  106.   imgSize.x = ImageGetWidth (imgIO);
  107.   imgSize.y = ImageGetHeight (imgIO);
  108.   if (imgSize.y > MAXDISPLAY)
  109.     displayFlag = 0;
  110.   ySizeM1 = imgSize.y - 1;
  111.   xSizeM1 = imgSize.x - 1;
  112.  
  113. /* invert image */
  114.   if (invertFlag) {
  115.     for (y = 0; y < imgSize.y; y++)
  116.       for (x = 0; x < imgSize.x; x++)
  117.       image[y][x] = 255 - image[y][x];
  118.   }
  119.   xRun = (long **) calloc (imgSize.y, sizeof (long));
  120.  
  121.   OFF = 0;
  122.   FILL0 = (unsigned char) OFF + 1;
  123.   ON = 255;
  124.   FILL1 = (unsigned char) (ON - 1);
  125.  
  126. /* make table of fill-values for FILL0 and FILL1 */
  127.   printf ("constructing fill tables:\n");
  128.   f0Table = (unsigned char **) calloc (kMax - 2, sizeof (unsigned char *));
  129.   f1Table = (unsigned char **) calloc (kMax - 2, sizeof (unsigned char *));
  130.   for (k = 3; k <= kMax; k++) {
  131.     nTable = (long) pow (2.0, 4.0 * ((double) k - 1.0));
  132.     f0Table[k - 3] = (unsigned char *) calloc (nTable, sizeof (unsigned char));
  133.     f1Table[k - 3] = (unsigned char *) calloc (nTable, sizeof (unsigned char));
  134.     for (i = 0; i < nTable; i++) {
  135.       f0Table[k - 3][i] = filltest (i, k, (long) FILL0);
  136.       f1Table[k - 3][i] = filltest (i, k, (long) FILL1);
  137.  
  138.     }
  139.   }
  140.  
  141. /* zero image borders */
  142.   for (y = 0; y < imgSize.y; y++)
  143.     image[y][0] = image[y][imgSize.x - 1] = (unsigned char) OFF;
  144.   for (x = 0; x < imgSize.x; x++)
  145.     image[0][x] = image[ySizeM1][x] = (unsigned char) OFF;
  146.  
  147.   printf ("performing filling:\n");
  148.   for (k = 0; k <= kMax; k++)
  149.     change[k] = 0;
  150.  
  151. /* iteratively convolve through image until filled */
  152.  
  153.   if (displayFlag != 0) {
  154.     display (fillFlag);
  155.     stall ("For next iteration");
  156.   }
  157.  
  158. /* on first iteration, perform filling and accumulate x-run info */
  159.   nChange = fill0 (f0Table, f1Table, fillFlag, kMax, change, &nONs, &nFill);
  160.   nChangeMax = nChange;
  161.   nChangeThresh = 0;
  162.   nChangeB4 = nChangeThresh + 1;
  163.   printf ("iteration  0: [%1d]\t", fillFlag);
  164.   for (i = 3; i <= kMax; i++) {
  165.     printf (" %d) %3d;  ", i, change[i]);
  166.     change[i] = 0;
  167.   }
  168.   printf ("\n");
  169.  
  170.   if (displayFlag != 0 && nChange > 0) {
  171.     display (fillFlag);
  172.     stall ("For next iteration");
  173.   }
  174.   if (fillFlag == 1) {
  175.     fillFlag = 0;
  176.     FILL0++;
  177.   }
  178.   else {
  179.     fillFlag = 1;
  180.     --FILL1;
  181.   }
  182. /* on subsequent iterations, perform filling */
  183.   for (nIter = 1; nIter < maxIter &&
  184.        (nChange > nChangeThresh || nChangeB4 > nChangeThresh); nIter++) {
  185.     nChangeB4 = nChange;
  186.     nChange = fill (f0Table, f1Table, fillFlag, kMax, change, &nFill);
  187.     printf ("iteration %2d: [%1d]\t", nIter + 1, fillFlag);
  188.     for (i = 3; i <= kMax; i++) {
  189.       printf (" %d) %3d;  ", i, change[i]);
  190.       change[i] = 0;
  191.     }
  192.     printf ("\n");
  193.  
  194.     if (displayFlag != 0 && nChange > 0) {
  195.       display (fillFlag);
  196.       stall ("For next iteration");
  197.     }
  198.  
  199.     if (fillFlag == 1) {
  200.       FILL0++;
  201.       if (nChangeB4 != 0)
  202.         fillFlag = 0;
  203.     }
  204.     else {
  205.       --FILL1;
  206.       if (nChangeB4 != 0)
  207.         fillFlag = 1;
  208.     }
  209.     if (nChange > nChangeMax)
  210.       nChangeMax = nChange;
  211.     nChangeThresh = (pctDirty * nChangeMax) / 100;
  212.   }
  213.   if (nIter >= maxIter && (nChange != 0 || nChangeB4 != 0))
  214.     printf ("maximum iterations reached\n");
  215.  
  216.   for (y = 1; y < ySizeM1; y++) {
  217.     for (x = 1; x < imgSize.x - 1; x++) {
  218.       if (image[y][x] <= FILL0)
  219.         image[y][x] = (unsigned char) OFF;
  220.       else if (image[y][x] >= FILL1)
  221.         image[y][x] = (unsigned char) ON;
  222.     }
  223.   }
  224.   if (displayFlag != 0)
  225.     display (fillFlag);
  226.  
  227. /* un-invert image */
  228.   if (invertFlag) {
  229.     for (y = 0; y < imgSize.y; y++)
  230.       for (x = 0; x < imgSize.x; x++)
  231.       image[y][x] = 255 - image[y][x];
  232.   }
  233.   ImageOut (argv[2], imgIO);
  234.  
  235.   return (0);
  236. }
  237.  
  238.  
  239.  
  240. /* FILLTEST:    function makes decision to fill or not based on CNUM
  241.  *            and FNUM in perimeter ring
  242.  *                      usage: fillFlag = filltest (pack, k, fill01)
  243.  *              fillFlag = 0 if no fill, 1 if OFF fill, 2 if ON fill,
  244.  *                              3 if ON and OFF fill
  245.  */
  246.  
  247. int
  248. filltest (pack, k, fill01)
  249.      register long pack,        /* packed ring of perimeter pixelsr */
  250.        k,                       /* square sidelength of ring */
  251.        fill01;                  /* fill value: FILL0 or FILL1 */
  252. {
  253.   register nRing,               /* no. pixels in ring */
  254.     n, i;
  255.   unsigned char *ring;          /* neighborhood pixels array */
  256.   long fNum,                    /* no. 1s on ring */
  257.     cNum,                       /* connectivity number */
  258.     m;
  259.   long lower, upper;            /* adjacent ring elements for cNum calc */
  260.   long nCornerOn,               /* no. corners ON */
  261.     fNumThresh;                 /* threshold of fNum */
  262.  
  263.   ring = (unsigned char *) malloc (4 * (k - 1));
  264.  
  265. /* unpack ring from word to array */
  266.   nRing = 4 * k - 4;
  267.   for (i = 0; i < nRing; i++)
  268.     ring[i] = (pack >> i) & 01;
  269.  
  270. /* calculate CNUM, first skipping corners */
  271.   for (i = 2, cNum = 0; i < nRing; i++) {
  272.     lower = (long) ring[i - 1];
  273.     if ((i % (k - 1)) == 0)
  274.       i++;                      /* skip the corner pixels */
  275.     upper = (long) ring[i];
  276.     if (upper != 0 && lower == 0)
  277.       cNum++;
  278.   }
  279.   if (ring[1] != 0 && ring[nRing - 1] == 0)
  280.     cNum++;
  281.  
  282. /* CNUM at corners */
  283.   for (n = 1, nCornerOn = 0; n < 4; n++) {
  284.     m = n * (k - 1);
  285.     if (ring[m] != 0) {
  286.       if (ring[m - 1] == 0 && ring[m + 1] == 0)
  287.         cNum++;
  288.       nCornerOn++;
  289.     }
  290.   }
  291.   if (ring[0] != 0) {
  292.     if (ring[1] == 0 && ring[nRing - 1] == 0)
  293.       cNum++;
  294.     nCornerOn++;
  295.   }
  296.  
  297. /* calculate FNUM */
  298.   if (fill01 == (long) FILL1) {
  299.     for (i = 0, fNum = 0; i < nRing; i++)
  300.       if (ring[i] != 0)
  301.         fNum++;
  302.   }
  303.   else {
  304.     for (i = 0, fNum = 0; i < nRing; i++)
  305.       if (ring[i] == 0)
  306.         fNum++;
  307.   }
  308.  
  309. /* to fill or not to fill */
  310.   if (cFlag == 0 || (cFlag != 0 && cNum <= 1)) {
  311.     fNumThresh = 3 * (k - 1) - 1 + fDiff;
  312.     if (fill01 == (long) FILL1 || eFlag == 0) {
  313.       if (fNum > fNumThresh)
  314.         return (1);
  315.       if (fNum == fNumThresh && nCornerOn == 2)
  316.         return (1);
  317.     }
  318.     else {
  319.       if (fNum == nRing)
  320.         return (1);
  321.       if (fNum == fNumThresh && nCornerOn == 2)
  322.         return (1);
  323.     }
  324.   }
  325.   return (0);
  326. }
  327.  
  328.  
  329.  
  330. /* FILL0:       function performs first iteration of thinning, as well
  331.  *            as compiling run-lengths
  332.  *                  usage: nChange = fill0 (f0Table, f1Table, fillFlag, kMax,
  333.  *                                              change, &nONs, &nFill)
  334.  */
  335.  
  336. int
  337. fill0 (f0Table, f1Table, fillFlag, kMax, change, nONs, nFill)
  338.      unsigned char **f0Table,   /* table of FILL0 values for each k */
  339.      **f1Table;                 /* table of FILL1 values for each k */
  340.      long fillFlag,             /* ON fill (if =1), or OFF erase (if =0) */
  341.        kMax,                    /* max sidelength of thinning kernel */
  342.       *change,                  /* no. erasures for each mask size */
  343.       *nONs,                    /* no. original ON pixels in image */
  344.       *nFill;                   /* cumulative no. filled this iteration */
  345. {
  346.   register long x, y,           /* image coordinates */
  347.     iXRun,                      /* index of runs in x */
  348.     k,                          /* sidelength of current thinning kernel */
  349.     kM1;                        /* k minus 1 */
  350.   long ring,                    /* word containing neighborhood pixels */
  351.     nChange,                    /* no. thinning operations on iteration */
  352.     fillValue,                  /* fillValue = 1 to fill or 0 */
  353.     onRun;                      /* flag = 1 for run of 1s; 0 for 0s */
  354.  
  355.   *nONs = nChange = 0;
  356.   for (y = 1; y < ySizeM1; y++) {
  357.     xRun[y] = (long *) calloc (imgSize.x + 1, sizeof (long));
  358.     xRun[y][0] = -MIN0RUN;
  359.  
  360.     for (x = 1, iXRun = 1, onRun = 0; x < xSizeM1; x++) {
  361.       if (image[y][x] <= FILL0) {
  362.         if (onRun == 1) {
  363.           onRun = 0;
  364.           xRun[y][iXRun++] =
  365.             (x - 1 >= imgSize.x) ? xSizeM1 : x - 1;
  366.         }
  367.       }
  368.       else {
  369.         if (onRun == 0) {
  370.           onRun = 1;
  371.           if ((x - xRun[y][iXRun - 1]) < MIN0RUN)
  372.             --iXRun;
  373.           else
  374.             xRun[y][iXRun++] = (x < 0) ? 1 : x;
  375.         }
  376.         (*nONs)++;
  377.       }
  378.       k = ksize (x, y, kMax, fillFlag);
  379.       kM1 = (k > 3) ? k - 1 : 3;
  380.       while (k >= kM1) {
  381.         getring (x, y, k, fillFlag, &ring);
  382.         fillValue = (fillFlag == 0) ?
  383.           f0Table[k - 3][ring] : f1Table[k - 3][ring];
  384.         if (fillValue == 1) {
  385.           nChange++;
  386.           (change[k])++;
  387.           fillsqr (x, y, k, fillFlag, nFill);
  388.           break;
  389.         }
  390.         --k;
  391.       }
  392.     }
  393.     --iXRun;
  394.     if (iXRun % 2 != 0)
  395.       xRun[y][++iXRun] = x;
  396.     xRun[y][0] = iXRun;
  397.     xRun[y] = (long *) realloc (xRun[y], (sizeof (long)) * (iXRun + 1));
  398.   }
  399.  
  400.   return (nChange);
  401. }
  402.  
  403.  
  404.  
  405. /* FILL:        function performs an iteration of thinning
  406.  *                usage: nChange = fill (f0Table, f1Table, fillFlag, kMax,
  407.  *                                                      change, &nFill)
  408.  */
  409.  
  410. int
  411. fill (f0Table, f1Table, fillFlag, kMax, change, nFill)
  412.      unsigned char **f0Table,   /* table of FILL0 values for each k */
  413.      **f1Table;                 /* table of FILL1 values for each k */
  414.      long fillFlag,             /* ON fill (if =1), or OFF erase (if =0) */
  415.        kMax,                    /* max sidelength of thinning kernel */
  416.       *change,                  /* no. erasures for each mask size */
  417.       *nFill;                   /* cumulative no. filled this iteration */
  418. {
  419.   register long x, y,           /* image coordinates */
  420.     xStart, xEnd,               /* start and end of x-run */
  421.     iXRun,                      /* index of runs in x */
  422.     k,                          /* sidelength of current thinning kernel */
  423.     kM1;                        /* k minus 1 */
  424.   long ring,                    /* word containing nbrhood pixels */
  425.     nChange,                    /* no. thinning operations on iteration */
  426.     fillValue;                  /* fillValue = 1 to fill or 0 */
  427.  
  428.   nChange = 0;
  429.   for (y = 1; y < ySizeM1; y++) {
  430.     for (iXRun = 1, x = 1; iXRun <= xRun[y][0]; iXRun += 2) {
  431.       xStart = xRun[y][iXRun] - kMax + 2;
  432.       xStart = (xStart > x) ? xStart : x;
  433.       xEnd = xRun[y][iXRun + 1] + kMax - 2;
  434.       if (xEnd > xSizeM1)
  435.         xEnd = xSizeM1;
  436.       for (x = xStart; x <= xEnd; x++) {
  437.         k = ksize (x, y, kMax, fillFlag);
  438.         kM1 = (k > 3) ? k - 1 : 3;
  439.         while (k >= kM1) {
  440.           getring (x, y, k, fillFlag, &ring);
  441.           fillValue = (fillFlag == 0) ?
  442.             f0Table[k - 3][ring] : f1Table[k - 3][ring];
  443.           if (fillValue == 1) {
  444.             nChange++;
  445.             (change[k])++;
  446.             fillsqr (x, y, k, fillFlag, nFill);
  447.             break;
  448.           }
  449.           --k;
  450.         }
  451.       }
  452.     }
  453.   }
  454.  
  455.   return (nChange);
  456. }
  457.  
  458.  
  459.  
  460. /* KSIZE:       function determines k, where kxk is largest square
  461.  *            around (x,y) which contains all OFF if fillFlag=0,
  462.  *              or all ON for fillFlag=1
  463.  *                      usage: k = ksize (x, y, kMax, fillFlag)
  464.  */
  465.  
  466. long
  467. ksize (x, y, kMax, fillFlag)
  468.      long x, y,                 /* image coordinates */
  469.        kMax,                    /* maximum k value */
  470.        fillFlag;                /* ON fill (if =1), or OFF erase (if =0) */
  471. {
  472.   register long xMask, yMask,   /* x,y mask coordinates */
  473.     xEnd, yEnd,                 /* end coord.s of square */
  474.     k;                          /* mask size */
  475.   long upHalf, downHalf,        /* half of mask below and above center */
  476.     xStart,                     /* x- start and end of square */
  477.     yStart;                     /* y- start and end of square */
  478.  
  479.   if (fillFlag == 0) {
  480.     if (image[y][x] <= FILL0)
  481.       return (0);
  482.     else if (kMax == 3)
  483.       return (3);
  484.   }
  485.   else {
  486.     if (image[y][x] >= FILL1)
  487.       return (0);
  488.     else if (kMax == 3)
  489.       return (3);
  490.   }
  491.  
  492.   for (k = 4; k <= kMax; k++) {
  493.     if (k % 2 == 1)
  494.       downHalf = upHalf = (k - 3) / 2;
  495.     else {
  496.       upHalf = (k - 2) / 2;
  497.       downHalf = (k - 4) / 2;
  498.     }
  499.     xStart = x - downHalf;
  500.     xEnd = x + upHalf;
  501.     yStart = y - downHalf;
  502.     yEnd = y + upHalf;
  503.     if (xStart <= 0 || yStart <= 0
  504.         || xEnd >= (imgSize.x - 1) || yEnd >= ySizeM1)
  505.       return (k - 1);
  506.     for (yMask = yStart; yMask <= yEnd; yMask++)
  507.       for (xMask = xStart; xMask <= xEnd; xMask++) {
  508.         if (fillFlag == 0) {
  509.           if (image[yMask][xMask] < FILL0)
  510.             return (k - 1);
  511.         }
  512.         else if (image[yMask][xMask] > FILL1)
  513.           return (k - 1);
  514.       }
  515.   }
  516.   return (kMax);
  517. }
  518.  
  519.  
  520.  
  521. /* GETRING:     function gets ring of pixels on perimeter of k-size square
  522.  *                usage:  getring (x, y, k, fillFlag, *ring)
  523.  */
  524.  
  525. int
  526. getring (x, y, k, fillFlag, ring)
  527.      register long x, y;        /* image coordinages */
  528.      long k,                    /* square sidelength of ring */
  529.        fillFlag;                /* ON fill (if =1), or OFF erase (if =0) */
  530.      long *ring;                /* ring of pixels on perimeter of kxk sqr */
  531. {
  532.   register xEnd, yEnd,          /* x,y ends of square */
  533.     i, xStart, yStart;          /* start and end of square */
  534.   long upHalf, downHalf;        /* half of mask below and above center */
  535.  
  536.   if (k % 2 == 1)
  537.     downHalf = upHalf = (k - 1) / 2;
  538.   else {
  539.     upHalf = k / 2;
  540.     downHalf = (k - 2) / 2;
  541.   }
  542.   xStart = x - downHalf;
  543.   xEnd = x + upHalf;
  544.   yStart = y - downHalf;
  545.   yEnd = y + upHalf;
  546.  
  547.   i = 0;
  548.   *ring = 0;
  549.  
  550.   if (fillFlag == 0) {
  551.     for (x = xStart, y = yStart; x <= xEnd; x++, i++)
  552.       if (image[y][x] >= FILL0)
  553.         *ring = *ring | (01 << i);
  554.     for (y = yStart + 1, x = xEnd; y <= yEnd; y++, i++)
  555.       if (image[y][x] >= FILL0)
  556.         *ring = *ring | (01 << i);
  557.     for (x = xEnd - 1, y = yEnd; x >= xStart; --x, i++)
  558.       if (image[y][x] >= FILL0)
  559.         *ring = *ring | (01 << i);
  560.     for (y = yEnd - 1, x = xStart; y > yStart; --y, i++)
  561.       if (image[y][x] >= FILL0)
  562.         *ring = *ring | (01 << i);
  563.   }
  564.   else {
  565.  
  566.     for (x = xStart, y = yStart; x <= xEnd; x++, i++)
  567.       if (image[y][x] > FILL1)
  568.         *ring = *ring | (01 << i);
  569.     for (y = yStart + 1, x = xEnd; y <= yEnd; y++, i++)
  570.       if (image[y][x] > FILL1)
  571.         *ring = *ring | (01 << i);
  572.     for (x = xEnd - 1, y = yEnd; x >= xStart; --x, i++)
  573.       if (image[y][x] > FILL1)
  574.         *ring = *ring | (01 << i);
  575.     for (y = yEnd - 1, x = xStart; y > yStart; --y, i++)
  576.       if (image[y][x] > FILL1)
  577.         *ring = *ring | (01 << i);
  578.   }
  579.  
  580.   return (0);
  581. }
  582.  
  583.  
  584.  
  585. /* FILLSQR:     function fills square with OFFs or ONs
  586.  *                    usage: fillsqr (x, y, k, fillFlag, &nFill)
  587.  */
  588.  
  589. int
  590. fillsqr (x, y, k, fillFlag, nFill)
  591.      register long x, y;        /* image coordinages */
  592.      long k;                    /* square sidelength of ring */
  593.      register long fillFlag;    /* ON fill (if =1), or OFF erase (if =0) */
  594.      long *nFill;               /* no. of filled */
  595. {
  596.   register long xEnd, yEnd;     /* upper bounds of center fill area */
  597.   long upHalf, downHalf,        /* half of mask below and above center */
  598.     yStart,                     /* bounds of center fill area */
  599.     xStart;
  600.  
  601. /* fill for 3x3 */
  602.   if (k == 3) {
  603.     if (fillFlag == 0) {
  604.       if (image[y][x] > FILL0) {
  605.         (*nFill)++;
  606.         image[y][x] = FILL0;
  607.       }
  608.     }
  609.     else {
  610.       if (image[y][x] < FILL1) {
  611.         (*nFill)++;
  612.         image[y][x] = FILL1;
  613.       }
  614.     }
  615.   }
  616.  
  617. /* fill for kxk > 3x3 */
  618.   else {
  619.     if (k % 2 == 1)
  620.       downHalf = upHalf = (k - 3) / 2;
  621.     else {
  622.       upHalf = (k - 2) / 2;
  623.       downHalf = (k - 4) / 2;
  624.     }
  625.     xStart = x - downHalf;
  626.     xEnd = x + upHalf;
  627.     yStart = y - downHalf;
  628.     yEnd = y + upHalf;
  629.     for (y = yStart; y <= yEnd; y++) {
  630.       for (x = xStart; x <= xEnd; x++) {
  631.         if (fillFlag == 0) {
  632.           if (image[y][x] > FILL0) {
  633.             (*nFill)++;
  634.             image[y][x] = FILL0;
  635.           }
  636.         }
  637.         else if (image[y][x] < FILL1) {
  638.           (*nFill)++;
  639.           image[y][x] = FILL1;
  640.         }
  641.       }
  642.     }
  643.   }
  644.   return (0);
  645. }
  646.  
  647.  
  648.  
  649. /* DISPLAY:     function displays results after each iteration of thinning
  650.  *                    usage: display (fillFlag)
  651.  */
  652.  
  653. int
  654. display (fillFlag)
  655.      long fillFlag;             /* ON fill (if =1), or OFF erase (if =0) */
  656. {
  657.   long y, x,                    /* image coordinates */
  658.     onThresh, offThresh;        /* on/off threshold for diff. fillFlag */
  659.  
  660.   if (fillFlag == 1) {
  661.     onThresh = (long) FILL1;
  662.     offThresh = (long) FILL0 + 1;
  663.   }
  664.   else {
  665.     onThresh = (long) FILL1 - 1;
  666.     offThresh = (long) FILL0;
  667.   }
  668.  
  669.   for (y = 0; y < imgSize.y; y++) {
  670.     for (x = 0; x < imgSize.x; x++) {
  671.       if (image[y][x] == onThresh)
  672.         printf ("* ");
  673.       else if ((long) image[y][x] > onThresh)
  674.         printf ("X ");
  675.       else if ((long) image[y][x] == offThresh)
  676.         printf ("- ");
  677.       else if ((long) image[y][x] < offThresh
  678.                && image[y][x] > (unsigned char) OFF)
  679.         printf (". ");
  680.       else
  681.         printf ("  ");
  682.     }
  683.     printf ("\n");
  684.   }
  685.   return (0);
  686. }
  687.  
  688.  
  689.  
  690. /* STALL: stall (or pause) until keyboard input
  691.  * log  12 Nov 81       */
  692.  
  693. int
  694. stall (prompt)
  695.      char *prompt;
  696. {
  697.   char c;
  698.  
  699.   printf (prompt);
  700.   printf (" -- enter <CR> to continue\n");
  701.   scanf ("%c", &c);
  702.   return (0);
  703. }
  704.  
  705.  
  706. /* USAGE:       function gives instructions on usage of program
  707.  *                    usage: usage (flag)
  708.  *              When flag is 1, the long message is given, 0 gives short.
  709.  */
  710.  
  711. long
  712. usage (flag)
  713.      short flag;                /* flag =1 for long message; =0 for short message */
  714. {
  715.  
  716. /* print short usage message or long */
  717.   printf ("USAGE: kfill inimg outimg [-f FDIFF] [-c] [-n NITER]\n");
  718.   printf ("                          [-r R%][-k K] [-l] [-d] [-I] [-L]\n");
  719.   if (flag == 0)
  720.     return (-1);
  721.  
  722.   printf ("\nkfill removes spatial noise from binary input image\n\n");
  723.   printf ("ARGUMENTS:\n");
  724.   printf ("    inimg: input image filename (TIF)\n");
  725.   printf ("   outimg: output image filename (TIF)\n\n");
  726.   printf ("OPTIONS:\n");
  727.   printf ("    -f FDIFF: amount of increase or decrease (negative) of\n");
  728.   printf ("              factor no. from default\n");
  729.   printf ("          -c: when set, retain connectivity (default is otherwise)\n");
  730.   printf ("          -e: when set, do not retain endpoints (default is otherwise)\n");
  731.   printf ("       -r R%: percentage of max noise to leave remaining on any\n");
  732.   printf ("              iteration before stopping (default is 0%)\n");
  733.   printf ("    -n NITER: maximum number of iterations (default max = %d)\n", DFLTITER);
  734.   printf ("        -k K: window size for kxk mask (k >= 3, default = %d)\n", DFLTKMAX);
  735.   printf ("          -l: first iteration fill-value opposite from default\n");
  736.   printf ("              (default is %d value)\n", FILLINITIAL);
  737.   printf ("          -d: to display results of each iteration (< 40x40 image)\n");
  738.   printf ("          -I: invert input image before processing\n");
  739.   printf ("          -L: print Software License for this module\n");
  740.  
  741.   return (-1);
  742. }
  743.  
  744.  
  745. /* INPUT:       function reads input parameters
  746.  *                  usage: input (argc, argv, &fDiff, &cFlag, &eFlag,
  747.  *                                      &maxIter, &pctDirty, &kMax,
  748.  *                                      &fillFlag, &displayFlag)
  749.  */
  750.  
  751. #define USAGE_EXIT(VALUE) {usage (VALUE); return (-1);}
  752.  
  753. long input
  754.   (argc, argv, fDiff, cFlag, eFlag, maxIter, pctDirty, kMax, fillFlag,
  755.    displayFlag, invertFlag)
  756.      int argc;
  757.      char *argv[];
  758.      long *fDiff,               /* factor no. difference */
  759.       *cFlag,                   /* if 1, retain connectivity; if 0, don't */
  760.       *eFlag,                   /* if 1, retain endpoint; if 0, don't */
  761.       *maxIter,                 /* max. no. iterations */
  762.       *pctDirty,                /* pct. of max noise left on an iteration */
  763.       *kMax,                    /* maximum value of k, mask sidelengths */
  764.       *fillFlag,                /* if =0, start with 0-fill; or =1 */
  765.       *displayFlag,             /* if 1, display each iteration, or 0 */
  766.       *invertFlag;              /* invert input image before processing */
  767. {
  768.   long n;
  769.   static char *fString = "-f";
  770.   static char *cString = "-c";
  771.   static char *eString = "-e";
  772.   static char *nString = "-n";
  773.   static char *rString = "-r";
  774.   static char *kString = "-k";
  775.   static char *lString = "-l";
  776.   static char *dString = "-d";
  777.   static char *LString = "-L";
  778.   static char *IString = "-I";
  779.   static char *dashString = "-";
  780.  
  781.  
  782.   if (argc == 1)
  783.     USAGE_EXIT (1);
  784.   if (argc == 2)
  785.     USAGE_EXIT (0);
  786.  
  787.   *fDiff = 0;
  788.   *cFlag = DFLTCFLAG;
  789.   *eFlag = DFLTEFLAG;
  790.   *maxIter = DFLTITER;
  791.   *pctDirty = DFLTDIRTY;
  792.   *kMax = DFLTKMAX;
  793.   *fillFlag = FILLINITIAL;
  794.   *displayFlag = DFLTDFLAG;
  795.   *invertFlag = 0;
  796.  
  797.   for (n = 3; n < argc; n++) {
  798.     if (strcmp (argv[n], fString) == 0) {
  799.       if (++n == argc)
  800.         usage (0);
  801.       *fDiff = atol (argv[n]);
  802.     }
  803.     else if (strcmp (argv[n], cString) == 0)
  804.       *cFlag = 0;
  805.     else if (strcmp (argv[n], eString) == 0)
  806.       *eFlag = 0;
  807.     else if (strcmp (argv[n], nString) == 0) {
  808.       if (++n == argc || argv[n][0] == '-')
  809.         USAGE_EXIT (0);
  810.       *maxIter = atol (argv[n]);
  811.     }
  812.     else if (strcmp (argv[n], rString) == 0) {
  813.       if (++n == argc || argv[n][0] == '-')
  814.         USAGE_EXIT (0);
  815.       *pctDirty = atol (argv[n]);
  816.     }
  817.     else if (strcmp (argv[n], kString) == 0) {
  818.       if (++n == argc || argv[n][0] == '-')
  819.         USAGE_EXIT (0);
  820.       *kMax = atol (argv[n]);
  821.     }
  822.     else if (strcmp (argv[n], lString) == 0)
  823.       *fillFlag = 1;
  824.     else if (strcmp (argv[n], dString) == 0)
  825.       *displayFlag = 1;
  826.     else if (strcmp (argv[n], IString) == 0)
  827.       *invertFlag = 1;
  828.     else if (strcmp (argv[n], LString) == 0) {
  829.       print_sos_lic ();
  830.       exit (0);
  831.     }
  832.     else
  833.       USAGE_EXIT (0);
  834.   }
  835.   if (*cFlag != 0)
  836.     *cFlag = 1;
  837.   if (*eFlag != 0)
  838.     *eFlag = 1;
  839.   if (*kMax < 2)
  840.     usage (1);
  841.   if (*fillFlag != 0)
  842.     *fillFlag = 1;
  843.   if (*displayFlag != 0)
  844.     *displayFlag = 1;
  845.   return (0);
  846. }
  847.