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)*(float)(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,l,c,ch;
  22   double pc;
  23   double *coeff;
  24   double res[MAXCHANNELS];
  25   i_img *timg;
  26   int radius, diameter;
  27
  28   mm_log((1,"i_gaussian(im %p, stdev %.2f)\n",im,stddev));
  29   i_clear_error();
  30
  31   if (stddev <= 0) {
  32     i_push_error(0, "stddev must be positive");
  33     return 0;
  34   }
  35   /* totally silly cutoff */
  36   if (stddev > 1000) {
  37     stddev = 1000;
  38   }
  39
  40   timg = i_sametype(im, im->xsize, im->ysize);
  41
  42   if (im->bits <= 8)
  43     radius = ceil(2 * stddev);
  44   else
  45     radius = ceil(3 * stddev);
  46
  47   diameter = 1 + radius * 2;
  48
  49   coeff = mymalloc(sizeof(double) * diameter);
  50
  51   for(i=0;i <= radius;i++)
  52     coeff[radius + i]=coeff[radius - i]=gauss(i, stddev);
  53   pc=0;
  54   for(i=0; i < diameter; i++)
  55     pc+=coeff[i];
  56   for(i=0;i < diameter;i++)
  57     coeff[i] /= pc;
  58
  59
  60 #code im->bits <= 8
  61   IM_COLOR rcolor;
  62   for(l=0;l<im->ysize;l++) {
  63     for(i=0;i<im->xsize;i++) {
  64       pc=0.0;
  65       for(ch=0;ch<im->channels;ch++)
  66         res[ch]=0;
  67       for(c = 0;c < diameter; c++)
  68         if (IM_GPIX(im,i+c-radius,l,&rcolor)!=-1) {
  69           for(ch=0;ch<im->channels;ch++)
  70             res[ch]+= rcolor.channel[ch] * coeff[c];
  71           pc+=coeff[c];
  72         }
  73       for(ch=0;ch<im->channels;ch++) {
  74         double value = res[ch] / pc;
  75         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
  76       }
  77       IM_PPIX(timg,i,l,&rcolor);
  78     }
  79   }
  80
  81   for(l=0;l<im->xsize;l++) {
  82     for(i=0;i<im->ysize;i++) {
  83       pc=0.0;
  84       for(ch=0; ch<im->channels; ch++)
  85         res[ch]=0;
  86       for(c=0; c < diameter; c++)
  87         if (IM_GPIX(timg,l,i+c-radius,&rcolor)!=-1) {
  88           for(ch=0;ch<im->channels;ch++)
  89             res[ch]+= rcolor.channel[ch] * coeff[c];
  90           pc+=coeff[c];
  91         }
  92       for(ch=0;ch<im->channels;ch++) {
  93         double value = res[ch]/pc;
  94         rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : IM_ROUND(value);
  95       }
  96       IM_PPIX(im,l,i,&rcolor);
  97     }
  98   }
  99 #/code
 100   myfree(coeff);
 101   i_img_destroy(timg);
 102
 103   return 1;
 104 }
 105
 106
 107
 108
 109
 110
 111