use a convolution kernel size based on the stddev rather than a

fixed size when performing a gaussian blur

also uses the i_fcolor interfaces when working with >8 bit/sample
images
RT 25645

diff --git a/Changes b/Changes
index 08e56eb13ac9071d1dce24d7957f5b56db59947f..3a69225664a090c0c791b9518c313d5ebfdce365 100644 (file)
--- a/Changes
+++ b/Changes
@@ -9,6 +9,10 @@ Bug fixes:
    Imager::Color checks when deciding whether to skip testing it
    http://rt.cpan.org/Ticket/Display.html?id=26064
 
+ - use a convolution kernel size based on the stddev rather than a
+   fixed size when performing a gaussian blur
+   http://rt.cpan.org/Ticket/Display.html?id=25645
+
 Imager 0.56 - 1 Apr 2007
 ===========
 

diff --git a/MANIFEST b/MANIFEST
index c6b78b09468840dc72fdfff60744ea494540d471..17c63250ddfa6918df9f53100f830701d66e4ae8 100644 (file)
--- a/MANIFEST
+++ b/MANIFEST
@@ -99,7 +99,7 @@ fontfiles/dcr10.afm
 fontfiles/dcr10.pfb
 fontfiles/dodge.ttf
 freetyp2.c     Implements freetype2 font support
-gaussian.c
+gaussian.im
 gif.c
 hlines.c       Manage sets of horizontal line segments
 image.c

diff --git a/gaussian.c b/gaussian.c
deleted file mode 100644 (file)
index ce0c295..0000000
--- a/gaussian.c
+++ /dev/null
@@ -1,72 +0,0 @@
-#include "imager.h"
-
-static float
-gauss(int x,float std) {
-  return 1.0/(sqrt(2.0*PI)*std)*exp(-(float)(x)*(float)(x)/(2*std*std));
-}
-
-/* Counters are as follows
- l:  lines
- i:  columns
- c:  filter coeffs
- ch: channels
- pc: coeff equalization
-*/
-
-
-
-void
-i_gaussian(i_img *im,float stdev) {
-  int i,l,c,ch;
-  float pc;
-  float coeff[21];
-  i_color rcolor;
-  float res[11];
-  i_img timg;
-
-  mm_log((1,"i_gaussian(im %p, stdev %.2f)\n",im,stdev));
-  
-  i_img_empty_ch(&timg,im->xsize,im->ysize,im->channels);
-             
-  for(i=0;i<11;i++) coeff[10+i]=coeff[10-i]=gauss(i,stdev);
-  pc=0;
-  for(i=0;i<21;i++) pc+=coeff[i];
-  for(i=0;i<21;i++) coeff[i]/=pc;
-
-
-  for(l=0;l<im->ysize;l++) {
-    for(i=0;i<im->xsize;i++) {
-      pc=0.0;
-      for(ch=0;ch<im->channels;ch++) res[ch]=0; 
-      for(c=0;c<21;c++)
-       if (i_gpix(im,i+c-10,l,&rcolor)!=-1) {
-         for(ch=0;ch<im->channels;ch++) res[ch]+=(float)(rcolor.channel[ch])*coeff[c];
-         pc+=coeff[c];
-       }
-      for(ch=0;ch<im->channels;ch++) rcolor.channel[ch]=(unsigned char)(((res[ch]/(float)(pc)>255.0)?255.0:res[ch]/(float)(pc)));
-      i_ppix(&timg,i,l,&rcolor);
-    }
-  }
-  
-  for(l=0;l<im->xsize;l++) {
-    for(i=0;i<im->ysize;i++) {
-      pc=0.0;
-      for(ch=0;ch<im->channels;ch++) res[ch]=0; 
-      for(c=0;c<21;c++)
-       if (i_gpix(&timg,l,i+c-10,&rcolor)!=-1) {
-         for(ch=0;ch<im->channels;ch++) res[ch]+=(float)(rcolor.channel[ch])*coeff[c];
-         pc+=coeff[c];
-       }
-      for(ch=0;ch<im->channels;ch++) rcolor.channel[ch]=(unsigned char)(((res[ch]/(float)(pc)>255.0)?255.0:res[ch]/(float)(pc)));
-      i_ppix(im,l,i,&rcolor);
-    }
-  }
-  i_img_exorcise(&timg);
-}
-
-
-
-
-
-
-

diff --git a/gaussian.im b/gaussian.im
new file mode 100644 (file)
index 0000000..34fe813
--- /dev/null
+++ b/gaussian.im
@@ -0,0 +1,111 @@
+#include "imager.h"
+#include <math.h>
+
+static double
+gauss(int x, double std) {
+  return 1.0/(sqrt(2.0*PI)*std)*exp(-(double)(x)*(float)(x)/(2*std*std));
+}
+
+/* Counters are as follows
+ l:  lines
+ i:  columns
+ c:  filter coeffs
+ ch: channels
+ pc: coeff equalization
+*/
+
+
+
+int
+i_gaussian(i_img *im, double stddev) {
+  int i,l,c,ch;
+  double pc;
+  double *coeff;
+  double res[MAXCHANNELS];
+  i_img *timg;
+  int radius, diameter;
+
+  mm_log((1,"i_gaussian(im %p, stdev %.2f)\n",im,stddev));
+  i_clear_error();
+
+  if (stddev <= 0) {
+    i_push_error(0, "stddev must be positive");
+    return 0;
+  }
+  /* totally silly cutoff */
+  if (stddev > 1000) {
+    stddev = 1000;
+  }
+ 
+  timg = i_sametype(im, im->xsize, im->ysize);
+
+  if (im->bits <= 8)
+    radius = ceil(2 * stddev);
+  else
+    radius = ceil(3 * stddev);
+
+  diameter = 1 + radius * 2;
+
+  coeff = mymalloc(sizeof(double) * diameter);
+
+  for(i=0;i <= radius;i++) 
+    coeff[radius + i]=coeff[radius - i]=gauss(i, stddev);
+  pc=0;
+  for(i=0; i < diameter; i++)
+    pc+=coeff[i];
+  for(i=0;i < diameter;i++)
+    coeff[i] /= pc;
+
+
+#code im->bits <= 8
+  IM_COLOR rcolor;
+  for(l=0;l<im->ysize;l++) {
+    for(i=0;i<im->xsize;i++) {
+      pc=0.0;
+      for(ch=0;ch<im->channels;ch++) 
+       res[ch]=0; 
+      for(c = 0;c < diameter; c++)
+       if (IM_GPIX(im,i+c-radius,l,&rcolor)!=-1) {
+         for(ch=0;ch<im->channels;ch++) 
+           res[ch]+= rcolor.channel[ch] * coeff[c];
+         pc+=coeff[c];
+       }
+      for(ch=0;ch<im->channels;ch++) {
+       double value = res[ch] / pc;
+       rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : value;
+      }
+      IM_PPIX(timg,i,l,&rcolor);
+    }
+  }
+  
+  for(l=0;l<im->xsize;l++) {
+    for(i=0;i<im->ysize;i++) {
+      pc=0.0;
+      for(ch=0; ch<im->channels; ch++)
+       res[ch]=0; 
+      for(c=0; c < diameter; c++)
+       if (IM_GPIX(timg,l,i+c-radius,&rcolor)!=-1) {
+         for(ch=0;ch<im->channels;ch++) 
+           res[ch]+= rcolor.channel[ch] * coeff[c];
+         pc+=coeff[c];
+       }
+      for(ch=0;ch<im->channels;ch++) {
+       double value = res[ch]/pc;
+       rcolor.channel[ch] = value > IM_SAMPLE_MAX ? IM_SAMPLE_MAX : value;
+      }
+      IM_PPIX(im,l,i,&rcolor);
+    }
+  }
+#/code
+  myfree(coeff);
+  i_img_destroy(timg);
+  
+  return 1;
+}
+
+
+
+
+
+
+

diff --git a/t/t61filters.t b/t/t61filters.t
index dcbd431d9455a9feeba95c8611725ef7484dc171..f04c7fe2cbbbf2176e910e2763516898288aef26 100644 (file)
--- a/t/t61filters.t
+++ b/t/t61filters.t
@@ -1,8 +1,9 @@
 #!perl -w
 use strict;
 use Imager qw(:handy);
-use Test::More tests => 66;
+use Test::More tests => 69;
 Imager::init_log("testout/t61filters.log", 1);
+use Imager::Test qw(is_image_similar);
 # meant for testing the filters themselves
 my $imbase = Imager->new;
 $imbase->open(file=>'testout/t104.ppm') or die;
@@ -18,8 +19,16 @@ test($imbase, {type=>'contrast', intensity=>0.5},
 test($imbase, {type=>'conv', coef=>[ -0.5, 1, -0.5, ], },
      'testout/t61_conv.ppm');
 
-test($imbase, {type=>'gaussian', stddev=>5 },
-     'testout/t61_gaussian.ppm');
+{
+  my $gauss = test($imbase, {type=>'gaussian', stddev=>5 },
+                  'testout/t61_gaussian.ppm');
+
+  my $imbase16 = $imbase->to_rgb16;
+  my $gauss16 = test($imbase16,  {type=>'gaussian', stddev=>5 },
+                    'testout/t61_gaussian16.ppm');
+  is_image_similar($gauss, $gauss16, 200000, "8 and 16 gaussian match");
+}
+
 
 test($imbase, { type=>'gradgen', dist=>1,
                    xo=>[ 10,  10, 120 ],
@@ -244,6 +253,7 @@ sub test {
       skip("couldn't filter", 1);
     }
   }
+  $copy;
 }
 
 sub color_close {