1#define IMAGER_NO_CONTEXT
2#include "imager.h"
3#include <math.h>
4
5static double
6gauss(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
20int
21i_gaussian(i_img *im, double stddev) {
22  return i_gaussian2( im, stddev, stddev );
23}
24
25typedef struct s_gauss_coeff {
26    int diameter;
27    int radius;
28    double *coeff;
29} t_gauss_coeff;
30
31
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;
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
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  }
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
228