/*      Name:   getg.c

Created:        Thu May 29 00:37:44 1986
Last modified:  Sat May 31 20:34:30 1986

Purpose:        Get the laplacian of a Gaussian (not normalized).

Author:         Bill Hoff,2-114C,8645,3563478 (hoff) at uicsl
*/

#include <stdio.h>
#include <math.h>
#include "numerical.h"
#include <grass/gmath.h>

int 
getg (double w, double *g[2], int size)
{
  long i, j, totsize, n, g_row;
  float rsq, sigma, two_ssq, val, sum = 0.0;

  totsize = size*size;
  n = size/2;
  for (i=0; i < totsize; i++)
    { *(g[0]+i) = 0.0; *(g[1]+i) = 0.0;}

  sigma = w / (2.0*sqrt((double) 2.0));
  two_ssq = 2.0*sigma*sigma;
  for (i=0; i < n; i++) {
    g_row = i*size;             /* start of row */
    for (j=0; j < n; j++)       {
      rsq = i*i + j*j;
      val = (rsq/two_ssq - 1)*exp(-rsq/two_ssq);
      *(g[0]+g_row+j) = val;
      sum += val;
      /* reflect into other quadrants */
      if (j > 0) {
        *(g[0]+g_row+(size-j)) = val;
        sum += val;
      }
      if (i > 0)        {
        *(g[0]+(size-i)*size+j) = val;
        sum += val;
      }
      if (i > 0 && j > 0)       {
        *(g[0]+(size-i)*size+(size-j)) = val;
        sum += val;
      }
    }
  }

  *(g[0]+0) -= sum;             /* make sure sum of all values is zero */

    return 0;
}