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