home *** CD-ROM | disk | FTP | other *** search
/ BCI NET / BCI NET Dec 94.iso / archives / programming / source / fbm12s.lha / flconv.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-07-18  |  8.4 KB  |  335 lines
  1. /*****************************************************************
  2.  * flconv.c: FBM Release 1.2 18-Apr-93 Michael Mauldin
  3.  *
  4.  * Copyright (C) 1993 by Michael Mauldin.  Permission is granted
  5.  * to use this file in whole or in part for any purpose, educational,
  6.  * recreational or commercial, provided that this copyright notice
  7.  * is retained unchanged.  This software is available to all free of
  8.  * charge by anonymous FTP and in the UUNET archives.
  9.  *
  10.  * flconv.c: 
  11.  *    Apply arbitary convolutions to an image.  For multi-plane
  12.  *    images, the convolution is applied to each image separately.
  13.  *
  14.  *    Only odd kernel sizes are allowed (3x3, 5x5, etc.).  Of course,
  15.  *    you can always use a kernel one size larger with the rightmost
  16.  *    and bottommost values set to 0 to implement the smaller kernel.
  17.  *
  18.  *     Edges are set to the product of the edge pixel value and
  19.  *    the sum of the kernel (so for kernels that sum to zero, the
  20.  *    edges are set to black).
  21.  *
  22.  *    If absval is non-zero, then the absolute value of the weighted
  23.  *    sum is used, else it is set to 0.
  24.  *
  25.  *    Internally, the kernel is converted to fixed binary arithmetic,
  26.  *    with 8 bits used to represent the fraction, 23 for the integer,
  27.  *    and 1 for the sign.
  28.  *
  29.  * CONTENTS
  30.  *    convolve_fbm (input, output, kernel, n, absval)
  31.  *    FBM *input, *output;
  32.  *    double **kernel;
  33.  *
  34.  * EDITLOG
  35.  *    LastEditDate = Mon Jun 25 00:18:02 1990 - Michael Mauldin
  36.  *    LastFileName = /usr2/mlm/src/misc/fbm/flconv.c
  37.  *
  38.  * HISTORY
  39.  * 18-Apr-93  Michael Mauldin (mlm) at Carnegie-Mellon University
  40.  *    Created from flshrp.c
  41.  *****************************************************************/
  42.  
  43. # include <stdio.h>
  44. # include <math.h>
  45. # include <ctype.h>
  46. # include "fbm.h"
  47.  
  48. /****************************************************************
  49.  * convolve_fbm: apply a convolution
  50.  ****************************************************************/
  51.  
  52. #ifndef lint
  53. static char *fbmid =
  54. "$FBM flconv.c <1.2> 18-Apr-93  (C) 1993 by Michael Mauldin, source \
  55. code available free from MLM@CS.CMU.EDU and from UUNET archives$";
  56. #endif
  57.  
  58. convolve_fbm (input, output, kernel, n, absval)
  59. FBM *input, *output;
  60. int n, absval;
  61. double **kernel;
  62. { register int i, m, *ko, *kw, sum;
  63.   register unsigned char *ibm, *obm;
  64.   int j, k, n2, nsq, rowlen, plnlen, w, h, *koff, *kwgt;
  65.   double msum = 0.0, floor ();
  66.   
  67.   if ((n&1) == 0)
  68.   { fprintf (stderr, "convolve_fbm: only odd kernel sizes are allowed.\n");
  69.     return (0);
  70.   }
  71.  
  72.   fprintf (stderr, "\nIn flconv, kernel is:\n");
  73.  
  74.     /* Print the kernel */
  75.     for (j=0; j<n; j++)
  76.     { fprintf (stderr, "\t");
  77.       for (i=0; i<n; i++)
  78.       { fprintf (stderr, "% 6.2lf", kernel[j][i]); }
  79.       fprintf (stderr, "\n");
  80.     }
  81.   
  82.  
  83.   /* Allocate output */
  84.   output->hdr = input->hdr;
  85.   alloc_fbm (output);
  86.  
  87.   /* Cache important size values */
  88.   w = input->hdr.cols;
  89.   h = input->hdr.rows;
  90.   rowlen = input->hdr.rowlen;
  91.   plnlen = input->hdr.plnlen;
  92.   
  93.   /* Preprocess kernel, convert to fixed arithmetic */
  94.   n2 = n/2;
  95.   nsq = n * n;
  96.  
  97.   kwgt = (int *) malloc (nsq * sizeof (int *));
  98.   koff = (int *) malloc (nsq * sizeof (int *));
  99.   
  100.   for (j=0; j<n; j++)
  101.   { for (i=0; i<n; i++)
  102.     { msum += kernel[j][i];
  103.       kwgt[j*n + i] = floor (kernel[j][i] * 256.0 + 0.5);
  104.       koff[j*n + i] = (j-n2) * rowlen + (i-n2);
  105.     }  
  106.   }
  107.   
  108.   fprintf (stderr, "Kernel sums to %8.2lf\n", msum);
  109.   
  110.   /* Loop through each plane */
  111.   for (k=0; k<input->hdr.planes; k++)
  112.   {
  113.     /* Copy edges directly */
  114.     for (j=0; j<h; j++)
  115.     { ibm = &input->bm[k*plnlen + j*rowlen];
  116.       obm = &output->bm[k*plnlen + j*rowlen];
  117.  
  118.       for (m=0; m < n2; m++)
  119.       { obm[m] = ibm[m] * msum + 0.5;
  120.         obm[w-(m+1)] = ibm[w-(m+1)] * msum + 0.5;
  121.       }
  122.     }
  123.   
  124.     ibm = &input->bm[k*plnlen];
  125.     obm = &output->bm[k*plnlen]; 
  126.       
  127.     for (i=0; i<w; i++)
  128.     { for (m=0; m < n2; m++)
  129.       { obm[i + m*rowlen] = ibm[i + m* rowlen] * msum + 0.5;
  130.         obm[i + (h-(m+1))*rowlen] = ibm[i + (h-(m+1))*rowlen] * msum + 0.5;
  131.       }
  132.     }
  133.  
  134.     /* Now do the convolution */
  135.     for (j = n2; j < h-n2; j++)
  136.     { ibm = &(input->bm[k*plnlen + j*rowlen]);
  137.       obm = &(output->bm[k*plnlen + j*rowlen]);
  138.       
  139.       for (i=n2; i < w-n2; i++)
  140.       { sum = 0;
  141.  
  142.     kw = kwgt; ko = koff;
  143.  
  144.     for (m = nsq; --m >= 0; )
  145.         { sum += *kw++ * ibm[i + *ko++]; }
  146.     
  147.     sum >>= 8; /* Divide by 256 */
  148.     
  149.     if (sum < BLACK) sum = absval ? -sum : BLACK;
  150.     if (sum > WHITE) sum = WHITE;
  151.     
  152.     obm[i] = sum;
  153.       }
  154.     }
  155.   }
  156.   
  157.   return (1);
  158. }
  159.  
  160. sobel_fbm (input, output)
  161. FBM *input, *output;
  162. { register int i, j, sum1, sum2;
  163.   register unsigned char *ibm, *obm;
  164.   int k, rowlen, plnlen, w, h;
  165.   
  166.   /* Allocate output */
  167.   output->hdr = input->hdr;
  168.   alloc_fbm (output);
  169.  
  170.   /* Cache important size values */
  171.   w = input->hdr.cols;
  172.   h = input->hdr.rows;
  173.   rowlen = input->hdr.rowlen;
  174.   plnlen = input->hdr.plnlen;
  175.   
  176.   /* Loop through each plane */
  177.   for (k=0; k<input->hdr.planes; k++)
  178.   {
  179.     /* Clear edges */
  180.     for (j=0; j<h; j++)
  181.     { output->bm[k*plnlen + j*rowlen] = 0;
  182.       output->bm[k*plnlen + j*rowlen + w-1] = 0;
  183.     }
  184.   
  185.     ibm = &input->bm[k*plnlen];
  186.     obm = &output->bm[k*plnlen]; 
  187.       
  188.     for (i=0; i<w; i++)
  189.     { obm[i] = 0;
  190.       obm[i + (h-1)*rowlen] = 0;
  191.     }
  192.  
  193.     /* Now do the convolution */
  194.     for (j = 1; j < h-1; j++)
  195.     { ibm = &(input->bm[k*plnlen + j*rowlen]);
  196.       obm = &(output->bm[k*plnlen + j*rowlen]);
  197.       
  198.       for (i=1; i < w-1; i++)
  199.       { /* Compute horizontal gradient */
  200.     sum1 = ibm[i-rowlen-1] + 2 * ibm[i-rowlen] + ibm[i-rowlen+1] +
  201.           -ibm[i+rowlen-1] - 2 * ibm[i+rowlen] - ibm[i+rowlen+1];
  202.  
  203.     /* Compute vertical gradient */
  204.     if (sum1 < 0) sum1 = -sum1;
  205.  
  206.         sum2 =      ibm[i-rowlen-1] - ibm[i-rowlen+1] +
  207.            2 * (ibm[i-1]        - ibm[i+1]) +
  208.                   ibm[i+rowlen-1] - ibm[i+rowlen+1];
  209.  
  210.     if (sum2 < 0) sum2 = -sum2;
  211.  
  212.     sum1 += sum2;
  213.     if (sum1 > WHITE) sum1 = WHITE;
  214.     obm[i] = sum1;
  215.       }
  216.     }
  217.   }
  218.   
  219.   return (1);
  220. }
  221.  
  222. sobel_5x5_fbm (input, output)
  223. FBM *input, *output;
  224. { register int i, j, sum1, sum2;
  225.   register unsigned char *ibm, *obm;
  226.   int k, rowlen, plnlen, w, h;
  227.   int max1=0, min1=0, max2=0, min2=0;
  228.   
  229.   /* Allocate output */
  230.   output->hdr = input->hdr;
  231.   alloc_fbm (output);
  232.  
  233.   /* Cache important size values */
  234.   w = input->hdr.cols;
  235.   h = input->hdr.rows;
  236.   rowlen = input->hdr.rowlen;
  237.   plnlen = input->hdr.plnlen;
  238.   
  239.   /* Loop through each plane */
  240.   for (k=0; k<input->hdr.planes; k++)
  241.   {
  242.     /* Clear edges */
  243.     for (j=0; j<h; j++)
  244.     { output->bm[k*plnlen + j*rowlen] = 0;
  245.       output->bm[k*plnlen + j*rowlen+1] = 0;
  246.       output->bm[k*plnlen + j*rowlen + w-2] = 0;
  247.       output->bm[k*plnlen + j*rowlen + w-1] = 0;
  248.     }
  249.   
  250.     ibm = &input->bm[k*plnlen];
  251.     obm = &output->bm[k*plnlen]; 
  252.       
  253.     for (i=0; i<w; i++)
  254.     { obm[i] = 0;
  255.       obm[i+rowlen] = 0;
  256.       obm[i + (h-2)*rowlen] = 0;
  257.       obm[i + (h-1)*rowlen] = 0;
  258.       
  259.     }
  260.  
  261.     /* Now do the convolution */
  262.     for (j = 2; j < h-2; j++)
  263.     { ibm = &(input->bm[k*plnlen + j*rowlen]);
  264.       obm = &(output->bm[k*plnlen + j*rowlen]);
  265.       
  266.       for (i=2; i < w-2; i++)
  267.       { /* Compute horizontal gradient */
  268.     sum1 =       -ibm[i-2*rowlen-1] +
  269.              -2 * ibm[i-2*rowlen] +
  270.                  -ibm[i-2*rowlen+1] +
  271.  
  272.                  -ibm[i-rowlen-2] +
  273.              -2 * ibm[i-rowlen-1] +
  274.              -4 * ibm[i-rowlen] +
  275.              -2 * ibm[i-rowlen+1] +
  276.                  -ibm[i-rowlen+2] +
  277.  
  278.                   ibm[i+rowlen-2] +
  279.               2 * ibm[i+rowlen-1] +
  280.               4 * ibm[i+rowlen] +
  281.               2 * ibm[i+rowlen+1] +
  282.                   ibm[i+rowlen+2] +
  283.  
  284.                   ibm[i+2*rowlen-1] +
  285.               2 * ibm[i+2*rowlen] +
  286.                   ibm[i+2*rowlen+1];
  287.  
  288.     if (sum1 < min1) min1 = sum1;
  289.     if (sum1 > max1) max1 = sum1;
  290.  
  291.     if (sum1 < 0) sum1 = -sum1;
  292.  
  293.     /* Compute vertical gradient */
  294.         sum2 =       -ibm[i-2*rowlen-1] +
  295.                   ibm[i-2*rowlen+1] +
  296.            
  297.                  -ibm[i-rowlen-2] +
  298.              -2 * ibm[i-rowlen-1] +
  299.               2 * ibm[i-rowlen+1] +
  300.                   ibm[i-rowlen+2] +
  301.            
  302.              -2 * ibm[i-2] +
  303.              -4 * ibm[i-1] +
  304.               4 * ibm[i+1] +
  305.               2 * ibm[i+2] +
  306.  
  307.                  -ibm[i+rowlen-2] +
  308.              -2 * ibm[i+rowlen-1] +
  309.               2 * ibm[i+rowlen+1] +
  310.                   ibm[i+rowlen+2] +
  311.            
  312.              -ibm[i+2*rowlen-1] +
  313.                   ibm[i+2*rowlen+1];
  314.  
  315.     if (sum2 < min2) min2 = sum2;
  316.     if (sum2 > max2) max2 = sum2;
  317.            
  318.     if (sum2 < 0) sum2 = -sum2;
  319.  
  320.     sum1 += sum2;
  321.     sum1 >>= 2;
  322.     if (sum1 > WHITE) sum1 = WHITE;
  323.  
  324.     obm[i] = sum1;
  325.       }
  326.     }
  327.   }
  328.   
  329.   fprintf (stderr, "Vertical range %d..%d, horizontal range %d..%d\n",
  330.        min1, max1, min2, max2);
  331.   
  332.   return (1);
  333. }
  334.  
  335.