001/*******************************************************************************
002 * This software is provided as a supplement to the authors' textbooks on digital
003 *  image processing published by Springer-Verlag in various languages and editions.
004 * Permission to use and distribute this software is granted under the BSD 2-Clause 
005 * "Simplified" License (see http://opensource.org/licenses/BSD-2-Clause). 
006 * Copyright (c) 2006-2016 Wilhelm Burger, Mark J. Burge. All rights reserved. 
007 * Visit http://imagingbook.com for additional details.
008 *******************************************************************************/
009
010package imagingbook.pub.sift.filters;
011
012/*
013 * TODO: Should use the GenericFilter class!
014 */
015
016import ij.IJ;
017import ij.plugin.filter.Convolver;
018import ij.process.FloatProcessor;
019
020public class GaussianFilter {
021        
022        static final double kernelSizeFactor = 3.5;   // times sigma
023        private float[] kernel1D;
024        
025        public GaussianFilter(double sigma) {
026                kernel1D = makeGaussKernel1d(sigma);
027        }
028        
029        public void applyTo(FloatProcessor fp) {
030                Convolver Conv = new Convolver();
031                Conv.setNormalize(true);        // this is important!
032                Conv.convolve(fp, kernel1D, 1, kernel1D.length);
033                Conv.convolve(fp, kernel1D, kernel1D.length, 1);
034        }
035        
036        public static float[] makeGaussKernel1d(double sigma){
037                //make 1D Gauss filter kernel large enough
038                //to avoid truncation effects (too small in ImageJ) 
039                int rad = (int) (kernelSizeFactor * sigma);
040                if (rad < 1) rad = 1;
041                int size = rad+rad+1;
042                float[] kernel = new float[size]; //center cell = kernel[rad]
043                double sigma2 = sigma * sigma;
044                double scale = 1.0 / (Math.sqrt(2 * Math.PI) * sigma);
045                
046                for (int i = 0; i < size; i++) {
047                        double x = rad - i;
048                        kernel[i] = (float) (scale * Math.exp(-0.5 * (x * x) / sigma2));
049                }
050                
051                return kernel;
052        }
053        
054        static float[][] makeGaussKernel2d(double sigma) {
055                int rad = (int) (kernelSizeFactor * sigma);
056                int size = rad+rad+1;
057                float[][] kernel = new float[size][size]; //center cell = kernel[rad][rad]
058                double sigma2 = sigma * sigma;
059                double scale = 1.0 / (2 * Math.PI * sigma * sigma);
060                double sum = 0;
061                
062                for (int i = 0; i < size; i++) {
063                        double x = rad - i;
064                        for (int j = 0; j < size; j++) {
065                                double y = rad - j;
066                                kernel[i][j] = (float) (scale * Math.exp(-0.5 * (x * x + y * y)
067                                                / sigma2));
068                                sum = sum + kernel[i][j];
069                        }
070                }
071                
072                // normalize the kernel:
073                for (int i = 0; i < size; i++) {
074                        for (int j = 0; j < size; j++) {
075                                kernel[i][j] = (float) (kernel[i][j] / sum);
076                        }
077                }
078                
079                return kernel;
080        }
081        
082        void printKernel(float[] kernel) {
083                System.out.println("****** Gaussian kernel ******* ");
084                for (int i=0; i<kernel.length; i++) {
085                        System.out.println(i+": "+kernel[i]);
086                }
087        }
088        
089        void printKernel(float[][] kernel) {
090                System.out.println("****** Gaussian kernel ******* ");
091                for (int i=0; i<kernel.length; i++) {
092                        for (int j=0; j<kernel[i].length; j++) {
093                                System.out.print(" "+kernel[i][j]);
094                        }
095                        System.out.println();
096                }
097        }
098        
099        void printKernel(float[] kernel, String title) {
100                IJ.log("****** " + title + " ******* ");
101                for (int i=0; i<kernel.length; i++) {
102                        IJ.log(i+": "+kernel[i]);
103                }
104        }
105        
106        static float sumKernel(float[] kernel) {
107                float sum = 0;
108                for (int i=0; i<kernel.length; i++) {
109                        sum = sum +kernel[i];
110                }
111                return sum;
112        }
113        
114        static float sumKernel(float[][] kernel) {
115                float sum = 0;
116                for (int i=0; i<kernel.length; i++) {
117                        for (int j=0; j<kernel[i].length; j++) {
118                                sum = sum +kernel[i][j];
119                        }
120                }
121                return sum;
122        }
123
124}