]>
Commit | Line | Data |
---|---|---|
e215ad76 | 1 | #define IMAGER_NO_CONTEXT |
73962964 TC |
2 | #include "imager.h" |
3 | #include <math.h> | |
4 | ||
5 | static double | |
6 | gauss(int x, double std) { | |
8d14daab | 7 | return 1.0/(sqrt(2.0*PI)*std)*exp(-(double)(x)*(double)(x)/(2*std*std)); |
73962964 TC |
8 | } |
9 | ||
10 | /* Counters are as follows | |
11 | l: lines | |
12 | i: columns | |
13 | c: filter coeffs | |
14 | ch: channels | |
15 | pc: coeff equalization | |
16 | */ | |
17 | ||
18 | ||
19 | ||
20 | int | |
21 | i_gaussian(i_img *im, double stddev) { | |
3d3b6bed TC |
22 | return i_gaussian2( im, stddev, stddev ); |
23 | } | |
73962964 | 24 | |
3d3b6bed TC |
25 | typedef struct s_gauss_coeff { |
26 | int diameter; | |
27 | int radius; | |
28 | double *coeff; | |
29 | } t_gauss_coeff; | |
73962964 | 30 | |
73962964 | 31 | |
3d3b6bed TC |
32 | static t_gauss_coeff *build_coeff( i_img *im, double stddev ) { |
33 | double *coeff = NULL; | |
34 | double pc; | |
35 | int radius, diameter, i; | |
36 | t_gauss_coeff *ret = mymalloc(sizeof(struct s_gauss_coeff)); | |
37 | ret->coeff = NULL; | |
73962964 TC |
38 | |
39 | if (im->bits <= 8) | |
40 | radius = ceil(2 * stddev); | |
41 | else | |
42 | radius = ceil(3 * stddev); | |
43 | ||
44 | diameter = 1 + radius * 2; | |
45 | ||
46 | coeff = mymalloc(sizeof(double) * diameter); | |
47 | ||
48 | for(i=0;i <= radius;i++) | |
49 | coeff[radius + i]=coeff[radius - i]=gauss(i, stddev); | |
3d3b6bed | 50 | pc=0.0; |
73962964 TC |
51 | for(i=0; i < diameter; i++) |
52 | pc+=coeff[i]; | |
3d3b6bed | 53 | for(i=0;i < diameter;i++) { |
73962964 | 54 | coeff[i] /= pc; |
3d3b6bed TC |
55 | // im_log((aIMCTX, 1, "i_gaussian2 Y i=%i coeff=%.2f\n", i, coeff[i] )); |
56 | } | |
57 | ||
58 | ret->diameter = diameter; | |
59 | ret->radius = radius; | |
60 | ret->coeff = coeff; | |
61 | return ret; | |
62 | } | |
63 | ||
64 | static void free_coeff(t_gauss_coeff *co ) { | |
65 | ||
66 | if( co->coeff != NULL ) | |
67 | myfree( co->coeff ); | |
68 | myfree( co ); | |
69 | } | |
70 | ||
71 | #define img_copy(dest, src) i_copyto( (dest), (src), 0,0, (src)->xsize,(src)->ysize, 0,0); | |
72 | ||
73 | ||
74 | ||
75 | int | |
76 | i_gaussian2(i_img *im, double stddevX, double stddevY) { | |
77 | int c, ch; | |
78 | i_img_dim x, y; | |
79 | double pc; | |
80 | t_gauss_coeff *co = NULL; | |
81 | double res[MAXCHANNELS]; | |
82 | i_img *timg; | |
83 | dIMCTXim(im); | |
84 | ||
85 | im_log((aIMCTX, 1,"i_gaussian2(im %p, stddev %.2f,%.2f)\n",im,stddevX,stddevY)); | |
86 | i_clear_error(); | |
87 | ||
88 | if (stddevX < 0) { | |
89 | i_push_error(0, "stddevX must be positive"); | |
90 | return 0; | |
91 | } | |
92 | if (stddevY < 0) { | |
93 | i_push_error(0, "stddevY must be positive"); | |
94 | return 0; | |
95 | } | |
73962964 | 96 | |
3d3b6bed TC |
97 | if( stddevX == stddevY && stddevY == 0 ) { |
98 | i_push_error(0, "stddevX or stddevY must be positive"); | |
99 | return 0; | |
100 | } | |
101 | ||
102 | ||
103 | /* totally silly cutoff */ | |
104 | if (stddevX > 1000) { | |
105 | stddevX = 1000; | |
106 | } | |
107 | if (stddevY > 1000) { | |
108 | stddevY = 1000; | |
109 | } | |
110 | ||
111 | timg = i_sametype(im, im->xsize, im->ysize); | |
112 | ||
113 | if( stddevX > 0 ) { | |
114 | /* Build Y coefficient matrix */ | |
115 | co = build_coeff( im, stddevX ); | |
116 | im_log((aIMCTX, 1, "i_gaussian2 X coeff radius=%i diamter=%i coeff=%p\n", co->radius, co->diameter, co->coeff)); | |
117 | } | |
118 | else { | |
119 | im_log((aIMCTX, 1, "i_gaussian2 X coeff is unity\n")); | |
120 | } | |
73962964 TC |
121 | |
122 | #code im->bits <= 8 | |
123 | IM_COLOR rcolor; | |
3d3b6bed TC |
124 | i_img *yin; |
125 | i_img *yout; | |
126 | ||
127 | if( stddevX > 0 ) { | |
128 | /******************/ | |
129 | /* Process X blur */ | |
130 | im_log((aIMCTX, 1, "i_gaussian2 X blur from im=%p to timg=%p\n", im, timg)); | |
131 | ||
8d14daab TC |
132 | for(y = 0; y < im->ysize; y++) { |
133 | for(x = 0; x < im->xsize; x++) { | |
73962964 TC |
134 | pc=0.0; |
135 | for(ch=0;ch<im->channels;ch++) | |
136 | res[ch]=0; | |
3d3b6bed TC |
137 | for(c = 0;c < co->diameter; c++) |
138 | if (IM_GPIX(im,x+c-co->radius,y,&rcolor)!=-1) { | |
139 | for(ch=0;ch<im->channels;ch++) | |
140 | res[ch]+= rcolor.channel[ch] * co->coeff[c]; | |
141 | pc+=co->coeff[c]; | |
73962964 TC |
142 | } |
143 | for(ch=0;ch<im->channels;ch++) { | |
144 | double value = res[ch] / pc; | |
bd8052a6 | 145 | rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value); |
73962964 | 146 | } |
8d14daab | 147 | IM_PPIX(timg, x, y, &rcolor); |
73962964 TC |
148 | } |
149 | } | |
3d3b6bed TC |
150 | /* processing is im -> timg=yin -> im=yout */ |
151 | yin = timg; | |
152 | yout = im; | |
153 | } | |
154 | else { | |
155 | /* processing is im=yin -> timg=yout -> im */ | |
156 | yin = im; | |
157 | yout = timg; | |
158 | } | |
73962964 | 159 | |
3d3b6bed TC |
160 | if( stddevY > 0 ) { |
161 | if( stddevX != stddevY ) { | |
162 | if( co != NULL ) { | |
163 | free_coeff(co); | |
164 | co = NULL; | |
165 | } | |
166 | ||
167 | /* Build Y coefficient matrix */ | |
168 | co = build_coeff( im, stddevY ); | |
169 | im_log((aIMCTX, 1, "i_gaussian2 Y coeff radius=%i diamter=%i coeff=%p\n", co->radius, co->diameter, co->coeff)); | |
170 | } | |
171 | ||
172 | /******************/ | |
173 | /* Process Y blur */ | |
174 | im_log((aIMCTX, 1, "i_gaussian2 Y blur from yin=%p to yout=%p\n", yin, yout)); | |
8d14daab TC |
175 | for(x = 0;x < im->xsize; x++) { |
176 | for(y = 0; y < im->ysize; y++) { | |
73962964 TC |
177 | pc=0.0; |
178 | for(ch=0; ch<im->channels; ch++) | |
179 | res[ch]=0; | |
3d3b6bed TC |
180 | for(c=0; c < co->diameter; c++) |
181 | if (IM_GPIX(yin, x, y+c-co->radius, &rcolor)!=-1) { | |
182 | for(ch=0;ch<yin->channels;ch++) | |
183 | res[ch]+= rcolor.channel[ch] * co->coeff[c]; | |
184 | pc+=co->coeff[c]; | |
73962964 | 185 | } |
3d3b6bed | 186 | for(ch=0;ch<yin->channels;ch++) { |
73962964 | 187 | double value = res[ch]/pc; |
bd8052a6 | 188 | rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value); |
73962964 | 189 | } |
3d3b6bed TC |
190 | IM_PPIX(yout, x, y, &rcolor); |
191 | } | |
73962964 | 192 | } |
3d3b6bed TC |
193 | if( im != yout ) { |
194 | im_log((aIMCTX, 1, "i_gaussian2 copying yout=%p to im=%p\n", yout, im)); | |
195 | img_copy( im, yout ); | |
73962964 | 196 | } |
3d3b6bed TC |
197 | } |
198 | else { | |
199 | im_log((aIMCTX, 1, "i_gaussian2 Y coeff is unity\n")); | |
200 | if( yin==timg ) { | |
201 | im_log((aIMCTX, 1, "i_gaussian2 copying timg=%p to im=%p\n", timg, im)); | |
202 | img_copy( im, timg ); | |
203 | } | |
204 | } | |
205 | ||
206 | im_log((aIMCTX, 1, "i_gaussian2 im=%p\n", im)); | |
207 | im_log((aIMCTX, 1, "i_gaussian2 timg=%p\n", timg)); | |
208 | im_log((aIMCTX, 1, "i_gaussian2 yin=%p\n", yin)); | |
209 | im_log((aIMCTX, 1, "i_gaussian2 yout=%p\n", yout)); | |
210 | ||
211 | ||
212 | ||
73962964 | 213 | #/code |
3d3b6bed TC |
214 | if( co != NULL ) |
215 | free_coeff(co); | |
216 | ||
73962964 TC |
217 | i_img_destroy(timg); |
218 | ||
219 | return 1; | |
220 | } | |
221 | ||
222 | ||
223 | ||
224 | ||
225 | ||
226 | ||
227 |