1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
|
/* 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 <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;
}
|
