/*  Convolve an image with the derivatives of Gaussians; part of detect-lines.
    Copyright (C) 1996-1998 Carsten Steger

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2, or (at your option)
    any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */

/* Changes Made by R. Balasubramanian for incorporating the the detect lines code to incorporate
   within GRASP (May 10th 1999) 

   The type of the pointer image has been changed to float from unsigned char in routines
   convolve_rows_gauss() and convolve_gauss()

   */

#include "lines.h"



/* Constants */

/* 1/sqrt(2*PI) */
#define SQRT_2_PI_INV 0.398942280401432677939946059935



/* Local function prototypes */

static void convolve_rows_gauss __P((float *image, double *mask, long n,
                                     float *h, long width, long height));
static void convolve_cols_gauss __P((float *h, double *mask, long n, float *k,
                                     long width, long height));
static double *compute_gauss_mask_0 __P((long *num, double sigma));
static double *compute_gauss_mask_1 __P((long *num, double sigma));
static double *compute_gauss_mask_2 __P((long *num, double sigma));



/* Functions to compute the integral, and the 0th and 1st derivative of the
   Gaussian function 1/(sqrt(2*PI)*sigma)*exp(-0.5*x^2/sigma^2) */

/* Integral of the Gaussian function */
double phi0(x,sigma)
  double x;
  double sigma;
{
  return normal(x/sigma);
}



/* The Gaussian function */
double phi1(x,sigma)
  double x;
  double sigma;
{
  double t;

  t = x/sigma;
  return SQRT_2_PI_INV/sigma*exp(-0.5*t*t);
}



/* First derivative of the Gaussian function */
double phi2(x,sigma)
  double x;
  double sigma;
{
  double t;

  t = x/sigma;
  return -x*SQRT_2_PI_INV/pow(sigma,3.0)*exp(-0.5*t*t);
}



/* Functions to compute the one-dimensional convolution masks of the 0th, 1st,
   and 2nd derivative of the Gaussian kernel for a certain smoothing level
   given by sigma.  The mask is allocated by the function and given as the
   return value.  The caller must ensure that this memory is freed.  The
   output is intended to be used as an array with range [-num:num].  Therefore,
   the caller should add num to the return value.  Examples for the calling
   sequence can be found in convolve_gauss.  Examples for the usage of the
   masks are given in convolve_rows_gauss and convolve_cols_gauss. */

/* Gaussian smoothing mask */
static double *compute_gauss_mask_0(num,sigma)
  long   *num;
  double sigma;
{
  long   i, n;
  double limit;
  double *h, *mask;

  limit = MASK_SIZE(MAX_SIZE_MASK_0,sigma); /* Error < 0.001 on each side */
  n = (long)limit;
  h = xcalloc(2*n+1,sizeof(double));
  mask = h + n;
  for (i=-n+1;i<=n-1;i++)
    mask[i] = phi0(-i+0.5,sigma) - phi0(-i-0.5,sigma);
  mask[-n] = 1.0 - phi0(n-0.5,sigma);
  mask[n] = phi0(-n+0.5,sigma);
  *num = n;
  return h;
}



/* First derivative of Gaussian smoothing mask */
static double *compute_gauss_mask_1(num,sigma)
  long   *num;
  double sigma;
{
  long   i, n;
  double limit;
  double *h, *mask;

  limit = MASK_SIZE(MAX_SIZE_MASK_1,sigma); /* Error < 0.001 on each side */
  n = (long)limit;
  h = xcalloc(2*n+1,sizeof(double));
  mask = h + n;
  for (i=-n+1;i<=n-1;i++)
    mask[i] = phi1(-i+0.5,sigma) - phi1(-i-0.5,sigma);
  mask[-n] = -phi1(n-0.5,sigma);
  mask[n] = phi1(-n+0.5,sigma);
  *num = n;
  return h;
}



/* Second derivative of Gaussian smoothing mask */
static double *compute_gauss_mask_2(num,sigma)
  long   *num;
  double sigma;
{
  long   i, n;
  double limit;
  double *h, *mask;

  limit = MASK_SIZE(MAX_SIZE_MASK_2,sigma); /* Error < 0.001 on each side */
  n = (long)limit;
  h = xcalloc(2*n+1,sizeof(double));
  mask = h + n;
  for (i=-n+1;i<=n-1;i++)
    mask[i] = phi2(-i+0.5,sigma) - phi2(-i-0.5,sigma);
  mask[-n] = -phi2(n-0.5,sigma);
  mask[n] = phi2(-n+0.5,sigma);
  *num = n;
  return h;
}





/* Convolve an image with the derivatives of a Gaussian smoothing kernel.
   Since all of the masks are separable, this is done in two steps in the
   function convolve_gauss.  Firstly, the rows of the image are convolved by
   an appropriate one-dimensional mask in convolve_rows_gauss, yielding an
   intermediate float-image h.  Then the columns of this image are convolved
   by another appropriate mask in convolve_cols_gauss to yield the final
   result k.  At the border of the image the gray values are mirrored. */

/* Convolve the rows of an image with the derivatives of a Gaussian. */
static void convolve_rows_gauss(image,mask,n,h,width,height)
  float    *image;
  double  *mask;
  long    n;
  float   *h;
  long    width;
  long    height;
{
  long      j, r, c, l;
  double    sum;

  /* Inner region */
  for (r=n; r<height-n; r++) {
    for (c=0; c<width; c++) {
      l = LINCOOR(r,c,width);
      sum = 0.0;
      for (j=-n;j<=n;j++)
        sum += (double)(image[l+j*width])*mask[j];
      h[l] = sum;
    }
  }
  /* Border regions */
  for (r=0; r<n; r++) {
    for (c=0; c<width; c++) {
      l = LINCOOR(r,c,width);
      sum = 0.0;
      for (j=-n;j<=n;j++)
        sum += (double)(image[LINCOOR(BR(r+j),c,width)])*mask[j];
      h[l] = sum;
    }
  }
  for (r=height-n; r<height; r++) {
    for (c=0; c<width; c++) {
      l = LINCOOR(r,c,width);
      sum = 0.0;
      for (j=-n;j<=n;j++)
        sum += (double)(image[LINCOOR(BR(r+j),c,width)])*mask[j];
      h[l] = sum;
    }
  }
}



/* Convolve the columns of an image with the derivatives of a Gaussian. */
static void convolve_cols_gauss(h,mask,n,k,width,height)
  float   *h;
  double  *mask;
  long    n;
  float   *k;
  long    width;
  long    height;
{
  long      j, r, c, l;
  double    sum;

  /* Inner region */
  for (r=0; r<height; r++) {
    for (c=n; c<width-n; c++) {
      l = LINCOOR(r,c,width);
      sum = 0.0;
      for (j=-n;j<=n;j++)
        sum += h[l+j]*mask[j];
      k[l] = (float)sum;
    }
  }
  /* Border regions */
  for (r=0; r<height; r++) {
    for (c=0; c<n; c++) {
      l = LINCOOR(r,c,width);
      sum = 0.0;
      for (j=-n;j<=n;j++)
        sum += h[LINCOOR(r,BC(c+j),width)]*mask[j];
      k[l] = (float)sum;
    }
  }
  for (r=0; r<height; r++) {
    for (c=width-n; c<width; c++) {
      l = LINCOOR(r,c,width);
      sum = 0.0;
      for (j=-n;j<=n;j++)
        sum += h[LINCOOR(r,BC(c+j),width)]*mask[j];
      k[l] = (float)sum;
    }
  }
}



/* Convolve an image with a derivative of the Gaussian. */
void convolve_gauss(image,k,width,height,sigma,deriv_type)
  float    *image;
  float   *k;
  long    width;
  long    height;
  double  sigma;
  long    deriv_type;
{
  double  *hr, *hc;
  double  *maskr, *maskc;
  long    nr, nc;
  float   *h;

  h = xcalloc(width*height,sizeof(float));

  switch (deriv_type) {
    case DERIV_R:
      hr = compute_gauss_mask_1(&nr,sigma);
      hc = compute_gauss_mask_0(&nc,sigma);
      break;
    case DERIV_C:
      hr = compute_gauss_mask_0(&nr,sigma);
      hc = compute_gauss_mask_1(&nc,sigma);
      break;
    case DERIV_RR:
      hr = compute_gauss_mask_2(&nr,sigma);
      hc = compute_gauss_mask_0(&nc,sigma);
      break;
    case DERIV_RC:
      hr = compute_gauss_mask_1(&nr,sigma);
      hc = compute_gauss_mask_1(&nc,sigma);
      break;
    case DERIV_CC:
      hr = compute_gauss_mask_0(&nr,sigma);
      hc = compute_gauss_mask_2(&nc,sigma);
      break;
  }

  maskr = hr + nr;
  maskc = hc + nc;

  convolve_rows_gauss(image,maskr,nr,h,width,height);
  convolve_cols_gauss(h,maskc,nc,k,width,height);

  free(h);
  free(hr);
  free(hc);
}