/*****************************************************************
 * flconv.c: FBM Release 1.2 18-Apr-93 Michael Mauldin
 *
 * Copyright (C) 1993 by Michael Mauldin.  Permission is granted
 * to use this file in whole or in part for any purpose, educational,
 * recreational or commercial, provided that this copyright notice
 * is retained unchanged.  This software is available to all free of
 * charge by anonymous FTP and in the UUNET archives.
 *
 * flconv.c: 
 *	Apply arbitary convolutions to an image.  For multi-plane
 *	images, the convolution is applied to each image separately.
 *
 *	Only odd kernel sizes are allowed (3x3, 5x5, etc.).  Of course,
 *	you can always use a kernel one size larger with the rightmost
 *	and bottommost values set to 0 to implement the smaller kernel.
 *
 * 	Edges are set to the product of the edge pixel value and
 *	the sum of the kernel (so for kernels that sum to zero, the
 *	edges are set to black).
 *
 *	If absval is non-zero, then the absolute value of the weighted
 *	sum is used, else it is set to 0.
 *
 *	Internally, the kernel is converted to fixed binary arithmetic,
 *	with 8 bits used to represent the fraction, 23 for the integer,
 *	and 1 for the sign.
 *
 * CONTENTS
 *	convolve_fbm (input, output, kernel, n, absval)
 *	FBM *input, *output;
 *	double **kernel;
 *
 * EDITLOG
 *	LastEditDate = Mon Jun 25 00:18:02 1990 - Michael Mauldin
 *	LastFileName = /usr2/mlm/src/misc/fbm/flconv.c
 *
 * HISTORY
 * 18-Apr-93  Michael Mauldin (mlm) at Carnegie-Mellon University
 *	Created from flshrp.c
 *****************************************************************/

# include <stdio.h>
# include <math.h>
# include <ctype.h>
# include "fbm.h"

/****************************************************************
 * convolve_fbm: apply a convolution
 ****************************************************************/

#ifndef lint
static char *fbmid =
"$FBM flconv.c <1.2> 18-Apr-93  (C) 1993 by Michael Mauldin, source \
code available free from MLM@CS.CMU.EDU and from UUNET archives$";
#endif

convolve_fbm (input, output, kernel, n, absval)
FBM *input, *output;
int n, absval;
double **kernel;
{ register int i, m, *ko, *kw, sum;
  register unsigned char *ibm, *obm;
  int j, k, n2, nsq, rowlen, plnlen, w, h, *koff, *kwgt;
  double msum = 0.0, floor ();
  
  if ((n&1) == 0)
  { fprintf (stderr, "convolve_fbm: only odd kernel sizes are allowed.\n");
    return (0);
  }

  fprintf (stderr, "\nIn flconv, kernel is:\n");

    /* Print the kernel */
    for (j=0; j<n; j++)
    { fprintf (stderr, "\t");
      for (i=0; i<n; i++)
      { fprintf (stderr, "% 6.2lf", kernel[j][i]); }
      fprintf (stderr, "\n");
    }
  

  /* Allocate output */
  output->hdr = input->hdr;
  alloc_fbm (output);

  /* Cache important size values */
  w = input->hdr.cols;
  h = input->hdr.rows;
  rowlen = input->hdr.rowlen;
  plnlen = input->hdr.plnlen;
  
  /* Preprocess kernel, convert to fixed arithmetic */
  n2 = n/2;
  nsq = n * n;

  kwgt = (int *) malloc (nsq * sizeof (int *));
  koff = (int *) malloc (nsq * sizeof (int *));
  
  for (j=0; j<n; j++)
  { for (i=0; i<n; i++)
    { msum += kernel[j][i];
      kwgt[j*n + i] = floor (kernel[j][i] * 256.0 + 0.5);
      koff[j*n + i] = (j-n2) * rowlen + (i-n2);
    }  
  }
  
  fprintf (stderr, "Kernel sums to %8.2lf\n", msum);
  
  /* Loop through each plane */
  for (k=0; k<input->hdr.planes; k++)
  {
    /* Copy edges directly */
    for (j=0; j<h; j++)
    { ibm = &input->bm[k*plnlen + j*rowlen];
      obm = &output->bm[k*plnlen + j*rowlen];

      for (m=0; m < n2; m++)
      { obm[m] = ibm[m] * msum + 0.5;
        obm[w-(m+1)] = ibm[w-(m+1)] * msum + 0.5;
      }
    }
  
    ibm = &input->bm[k*plnlen];
    obm = &output->bm[k*plnlen]; 
      
    for (i=0; i<w; i++)
    { for (m=0; m < n2; m++)
      { obm[i + m*rowlen] = ibm[i + m* rowlen] * msum + 0.5;
        obm[i + (h-(m+1))*rowlen] = ibm[i + (h-(m+1))*rowlen] * msum + 0.5;
      }
    }

    /* Now do the convolution */
    for (j = n2; j < h-n2; j++)
    { ibm = &(input->bm[k*plnlen + j*rowlen]);
      obm = &(output->bm[k*plnlen + j*rowlen]);
      
      for (i=n2; i < w-n2; i++)
      { sum = 0;

	kw = kwgt; ko = koff;

	for (m = nsq; --m >= 0; )
        { sum += *kw++ * ibm[i + *ko++]; }
	
	sum >>= 8; /* Divide by 256 */
	
	if (sum < BLACK) sum = absval ? -sum : BLACK;
	if (sum > WHITE) sum = WHITE;
	
	obm[i] = sum;
      }
    }
  }
  
  return (1);
}

sobel_fbm (input, output)
FBM *input, *output;
{ register int i, j, sum1, sum2;
  register unsigned char *ibm, *obm;
  int k, rowlen, plnlen, w, h;
  
  /* Allocate output */
  output->hdr = input->hdr;
  alloc_fbm (output);

  /* Cache important size values */
  w = input->hdr.cols;
  h = input->hdr.rows;
  rowlen = input->hdr.rowlen;
  plnlen = input->hdr.plnlen;
  
  /* Loop through each plane */
  for (k=0; k<input->hdr.planes; k++)
  {
    /* Clear edges */
    for (j=0; j<h; j++)
    { output->bm[k*plnlen + j*rowlen] = 0;
      output->bm[k*plnlen + j*rowlen + w-1] = 0;
    }
  
    ibm = &input->bm[k*plnlen];
    obm = &output->bm[k*plnlen]; 
      
    for (i=0; i<w; i++)
    { obm[i] = 0;
      obm[i + (h-1)*rowlen] = 0;
    }

    /* Now do the convolution */
    for (j = 1; j < h-1; j++)
    { ibm = &(input->bm[k*plnlen + j*rowlen]);
      obm = &(output->bm[k*plnlen + j*rowlen]);
      
      for (i=1; i < w-1; i++)
      { /* Compute horizontal gradient */
	sum1 = ibm[i-rowlen-1] + 2 * ibm[i-rowlen] + ibm[i-rowlen+1] +
	      -ibm[i+rowlen-1] - 2 * ibm[i+rowlen] - ibm[i+rowlen+1];

	/* Compute vertical gradient */
	if (sum1 < 0) sum1 = -sum1;

        sum2 =      ibm[i-rowlen-1] - ibm[i-rowlen+1] +
	       2 * (ibm[i-1]        - ibm[i+1]) +
  	            ibm[i+rowlen-1] - ibm[i+rowlen+1];

	if (sum2 < 0) sum2 = -sum2;

	sum1 += sum2;
	if (sum1 > WHITE) sum1 = WHITE;
	obm[i] = sum1;
      }
    }
  }
  
  return (1);
}

sobel_5x5_fbm (input, output)
FBM *input, *output;
{ register int i, j, sum1, sum2;
  register unsigned char *ibm, *obm;
  int k, rowlen, plnlen, w, h;
  int max1=0, min1=0, max2=0, min2=0;
  
  /* Allocate output */
  output->hdr = input->hdr;
  alloc_fbm (output);

  /* Cache important size values */
  w = input->hdr.cols;
  h = input->hdr.rows;
  rowlen = input->hdr.rowlen;
  plnlen = input->hdr.plnlen;
  
  /* Loop through each plane */
  for (k=0; k<input->hdr.planes; k++)
  {
    /* Clear edges */
    for (j=0; j<h; j++)
    { output->bm[k*plnlen + j*rowlen] = 0;
      output->bm[k*plnlen + j*rowlen+1] = 0;
      output->bm[k*plnlen + j*rowlen + w-2] = 0;
      output->bm[k*plnlen + j*rowlen + w-1] = 0;
    }
  
    ibm = &input->bm[k*plnlen];
    obm = &output->bm[k*plnlen]; 
      
    for (i=0; i<w; i++)
    { obm[i] = 0;
      obm[i+rowlen] = 0;
      obm[i + (h-2)*rowlen] = 0;
      obm[i + (h-1)*rowlen] = 0;
      
    }

    /* Now do the convolution */
    for (j = 2; j < h-2; j++)
    { ibm = &(input->bm[k*plnlen + j*rowlen]);
      obm = &(output->bm[k*plnlen + j*rowlen]);
      
      for (i=2; i < w-2; i++)
      { /* Compute horizontal gradient */
	sum1 =       -ibm[i-2*rowlen-1] +
	         -2 * ibm[i-2*rowlen] +
	             -ibm[i-2*rowlen+1] +

	             -ibm[i-rowlen-2] +
	         -2 * ibm[i-rowlen-1] +
	         -4 * ibm[i-rowlen] +
	         -2 * ibm[i-rowlen+1] +
	             -ibm[i-rowlen+2] +

	              ibm[i+rowlen-2] +
	          2 * ibm[i+rowlen-1] +
	          4 * ibm[i+rowlen] +
	          2 * ibm[i+rowlen+1] +
	              ibm[i+rowlen+2] +

	              ibm[i+2*rowlen-1] +
	          2 * ibm[i+2*rowlen] +
	              ibm[i+2*rowlen+1];

	if (sum1 < min1) min1 = sum1;
	if (sum1 > max1) max1 = sum1;

	if (sum1 < 0) sum1 = -sum1;

	/* Compute vertical gradient */
        sum2 =       -ibm[i-2*rowlen-1] +
	              ibm[i-2*rowlen+1] +
	       
	             -ibm[i-rowlen-2] +
	         -2 * ibm[i-rowlen-1] +
	          2 * ibm[i-rowlen+1] +
	              ibm[i-rowlen+2] +
	       
	         -2 * ibm[i-2] +
	         -4 * ibm[i-1] +
	          4 * ibm[i+1] +
	          2 * ibm[i+2] +

	             -ibm[i+rowlen-2] +
	         -2 * ibm[i+rowlen-1] +
	          2 * ibm[i+rowlen+1] +
	              ibm[i+rowlen+2] +
	       
		     -ibm[i+2*rowlen-1] +
	              ibm[i+2*rowlen+1];

	if (sum2 < min2) min2 = sum2;
	if (sum2 > max2) max2 = sum2;
	       
	if (sum2 < 0) sum2 = -sum2;

	sum1 += sum2;
	sum1 >>= 2;
	if (sum1 > WHITE) sum1 = WHITE;

	obm[i] = sum1;
      }
    }
  }
  
  fprintf (stderr, "Vertical range %d..%d, horizontal range %d..%d\n",
	   min1, max1, min2, max2);
  
  return (1);
}