revert threading changes, they aren't ready for the mainline yet
[imager.git] / gaussian.im
1 #include "imager.h"
2 #include <math.h>
3
4 static double
5 gauss(int x, double std) {
6   return 1.0/(sqrt(2.0*PI)*std)*exp(-(double)(x)*(double)(x)/(2*std*std));
7 }
8
9 /* Counters are as follows
10  l:  lines
11  i:  columns
12  c:  filter coeffs
13  ch: channels
14  pc: coeff equalization
15 */
16
17
18
19 int
20 i_gaussian(i_img *im, double stddev) {
21   int i, c, ch;
22   i_img_dim x, y;
23   double pc;
24   double *coeff;
25   double res[MAXCHANNELS];
26   i_img *timg;
27   int radius, diameter;
28
29   mm_log((1,"i_gaussian(im %p, stdev %.2f)\n",im,stddev));
30   i_clear_error();
31
32   if (stddev <= 0) {
33     i_push_error(0, "stddev must be positive");
34     return 0;
35   }
36   /* totally silly cutoff */
37   if (stddev > 1000) {
38     stddev = 1000;
39   }
40  
41   timg = i_sametype(im, im->xsize, im->ysize);
42
43   if (im->bits <= 8)
44     radius = ceil(2 * stddev);
45   else
46     radius = ceil(3 * stddev);
47
48   diameter = 1 + radius * 2;
49
50   coeff = mymalloc(sizeof(double) * diameter);
51
52   for(i=0;i <= radius;i++) 
53     coeff[radius + i]=coeff[radius - i]=gauss(i, stddev);
54   pc=0;
55   for(i=0; i < diameter; i++)
56     pc+=coeff[i];
57   for(i=0;i < diameter;i++)
58     coeff[i] /= pc;
59
60
61 #code im->bits <= 8
62   IM_COLOR rcolor;
63   for(y = 0; y < im->ysize; y++) {
64     for(x = 0; x < im->xsize; x++) {
65       pc=0.0;
66       for(ch=0;ch<im->channels;ch++) 
67         res[ch]=0; 
68       for(c = 0;c < diameter; c++)
69         if (IM_GPIX(im,x+c-radius,y,&rcolor)!=-1) {
70           for(ch=0;ch<im->channels;ch++) 
71             res[ch]+= rcolor.channel[ch] * coeff[c];
72           pc+=coeff[c];
73         }
74       for(ch=0;ch<im->channels;ch++) {
75         double value = res[ch] / pc;
76         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
77       }
78       IM_PPIX(timg, x, y, &rcolor);
79     }
80   }
81   
82   for(x = 0;x < im->xsize; x++) {
83     for(y = 0; y < im->ysize; y++) {
84       pc=0.0;
85       for(ch=0; ch<im->channels; ch++)
86         res[ch]=0; 
87       for(c=0; c < diameter; c++)
88         if (IM_GPIX(timg, x, y+c-radius, &rcolor)!=-1) {
89           for(ch=0;ch<im->channels;ch++) 
90             res[ch]+= rcolor.channel[ch] * coeff[c];
91           pc+=coeff[c];
92         }
93       for(ch=0;ch<im->channels;ch++) {
94         double value = res[ch]/pc;
95         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
96       }
97       IM_PPIX(im, x, y, &rcolor);
98     }
99   }
100 #/code
101   myfree(coeff);
102   i_img_destroy(timg);
103   
104   return 1;
105 }
106
107
108
109
110
111
112