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