Leolo's guassian2 patch
[imager.git] / gaussian.im
1 #define IMAGER_NO_CONTEXT
2 #include "imager.h"
3 #include <math.h>
4
5 static double
6 gauss(int x, double std) {
7   return 1.0/(sqrt(2.0*PI)*std)*exp(-(double)(x)*(double)(x)/(2*std*std));
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) {
22   return i_gaussian2( im, stddev, stddev );
23 }
24
25 typedef struct s_gauss_coeff {
26     int diameter;
27     int radius;
28     double *coeff;
29 } t_gauss_coeff;
30
31  
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;
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);
50   pc=0.0;
51   for(i=0; i < diameter; i++)
52     pc+=coeff[i];
53   for(i=0;i < diameter;i++) {
54     coeff[i] /= pc;
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   }
96
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   }
121
122 #code im->bits <= 8
123   IM_COLOR rcolor;
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
132   for(y = 0; y < im->ysize; y++) {
133     for(x = 0; x < im->xsize; x++) {
134       pc=0.0;
135       for(ch=0;ch<im->channels;ch++) 
136         res[ch]=0; 
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];
142         }
143       for(ch=0;ch<im->channels;ch++) {
144         double value = res[ch] / pc;
145         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
146       }
147       IM_PPIX(timg, x, y, &rcolor);
148     }
149   }
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   }
159   
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));
175   for(x = 0;x < im->xsize; x++) {
176     for(y = 0; y < im->ysize; y++) {
177       pc=0.0;
178       for(ch=0; ch<im->channels; ch++)
179         res[ch]=0; 
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];
185         }
186         for(ch=0;ch<yin->channels;ch++) {
187         double value = res[ch]/pc;
188         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
189       }
190         IM_PPIX(yout, x, y, &rcolor);
191       }
192     }
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 );
196   }
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
213 #/code
214   if( co != NULL )
215     free_coeff(co);
216
217   i_img_destroy(timg);
218   
219   return 1;
220 }
221
222
223
224
225
226
227