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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
|
/*
//
// Lynkeos
// $Id$
//
// Created by Jean-Etienne LAMIAUD on Sun Nov 4 2007.
// Copyright (c) 2007-2023. Jean-Etienne LAMIAUD
//
// 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 of the License, 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., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
//
*/
#include "ProcessingUtilities.h"
/* Cut the coordinate to the authorized range */
inline short range( long i, short l )
{
if ( i < 0 )
i = 0;
else if ( i >= l )
i = l-1;
return( i );
}
#define K_CUTOFF 2 /* The Gaussian is 1/K_CUTOFF at argument radius */
void MakeGaussian( REAL * const * const planes,
u_short width, u_short height, u_short nPlanes,
u_short lineWidth, double radius )
{
const double k = log(K_CUTOFF)/radius/radius;
u_short x, y, c;
// Fill the buffer with the gaussian
for( y = 0; y <= (height+1)/2; y++ )
{
for( x = 0; x <= (width+1)/2; x++ )
{
double d2, v;
d2 = x*x + y*y;
v = exp( -k*d2 );
for( c = 0; c < nPlanes; c++ )
{
planes[c][x+lineWidth*y] = v;
if ( x != 0 )
planes[c][width-x+lineWidth*y] = v;
if ( y != 0 )
planes[c][x+lineWidth*(height-y)] = v;
if ( x != 0 && y != 0 )
planes[c][width-x+lineWidth*(height-y)] = v;
}
}
}
}
void MakeGaussianSpectrum( LNKCOMPLEX * const * const planes,
u_short halfW, u_short height, u_short nPlanes,
u_short lineWidth, double radius )
{
const double k = M_PI*M_PI*radius*radius/4.0/M_LN2;
#warning Optimise by calculating line 0 alone, and reusing its value multiplied by the y factor on each other lines
// Fill the buffer with the gaussian
for( u_short y = 0; y < height; y++ )
{
double dy = y < height/2 ? y : height - y;
double fy2 = dy*dy/(double)height/(double)height;
for( u_short x = 0; x < halfW; x++ )
{
double f2, v;
f2 = (double)x*(double)x/(double)halfW/(double)halfW/4 + fy2;
v = exp( -k*f2 );
for( u_short c = 0; c < nPlanes; c++ )
planes[c][x+lineWidth*y] = v;
}
}
}
|
