add github action for CI on OS X
[imager.git] / gaussian.im
CommitLineData
e215ad76 1#define IMAGER_NO_CONTEXT
73962964
TC
2#include "imager.h"
3#include <math.h>
4
5static double
6gauss(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
20int
21i_gaussian(i_img *im, double stddev) {
3d3b6bed
TC
22 return i_gaussian2( im, stddev, stddev );
23}
73962964 24
3d3b6bed
TC
25typedef struct s_gauss_coeff {
26 int diameter;
27 int radius;
28 double *coeff;
29} t_gauss_coeff;
73962964 30
73962964 31
3d3b6bed
TC
32static 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
64static 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
75int
76i_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