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.lib.filters;
011
012public class GaussianFilter extends LinearFilter {
013
014        public GaussianFilter(double sigma) {
015                super(makeGaussKernel2d(sigma));
016        }
017        
018        public GaussianFilter(double sigmaX, double sigmaY) {
019                super(makeGaussKernel2d(sigmaX, sigmaY));
020        }
021        
022        public static float[] makeGaussKernel1d(double sigma){
023                // make 1D Gauss filter kernel large enough
024                // to avoid truncation effects (too small in ImageJ!) 
025                int rad = (int) (3.5 * sigma);
026                float[] kernel = new float[rad + 1 + rad]; // odd size
027                double sigma2 = sigma * sigma;
028                
029                for (int i = 0; i < kernel.length; i++) {
030                        double r = rad - i;
031                        kernel[i] =  (float) Math.exp(-0.5 * (r*r) / sigma2);
032                }
033                return kernel;
034        }
035
036        
037        public static float[][] makeGaussKernel2d(double sigma) {
038                int rad = (int) (3.5 * sigma);
039                int size = rad+rad+1;
040                float[][] kernel = new float[size][size]; //center cell = kernel[rad][rad]
041                double sigma2 = sigma * sigma;
042                double scale = 1.0 / (2 * Math.PI * sigma * sigma);
043                double sum = 0;
044                
045                for (int i = 0; i < size; i++) {
046                        double x = rad - i;
047                        for (int j = 0; j < size; j++) {
048                                double y = rad - j;
049                                kernel[i][j] = (float) (scale * Math.exp(-0.5 * (x * x + y * y)
050                                                / sigma2));
051                                sum = sum + kernel[i][j];
052                        }
053                }
054                
055                // normalize the kernel:
056                for (int i = 0; i < size; i++) {
057                        for (int j = 0; j < size; j++) {
058                                kernel[i][j] = (float) (kernel[i][j] / sum);
059                        }
060                }
061                
062                return kernel;
063        }
064        
065        public static float[][] makeGaussKernel2d(double sigmaX, double sigmaY){
066                int radX = (int) (3.5 * sigmaX);
067                int radY = (int) (3.5 * sigmaY);
068                int sizeX = radX + radX + 1;
069                int sizeY = radY + radY + 1;
070
071                float[][] kernel = new float[sizeX][sizeY]; //center cell = kernel[rad][rad]
072                double sigmaX2 = (sigmaX > 0.1) ? sigmaX * sigmaX : 0.1;
073                double sigmaY2 = (sigmaY > 0.1) ? sigmaY * sigmaY : 0.1;
074                
075                double sum = 0;
076                for (int i = 0; i < sizeX; i++) {
077                        double x = radX - i;
078                        for (int j = 0; j < sizeY; j++) {
079                                double y = radY - j;
080                                // IJ.log("x = " + x + " / " + "y = " + y);
081                                double g = (float) Math.exp(-((x * x) / (2 * sigmaX2) + (y * y) / (2 * sigmaY2)));
082                                // IJ.log("g = " + g);
083                                kernel[i][j] = (float) g;
084                                sum = sum + g;
085                        }
086                }
087
088                // normalize the kernel to sum 1
089                double scale = 1 / sum;
090                for (int i = 0; i < sizeX; i++) {
091                        for (int j = 0; j < sizeY; j++) {
092                                kernel[i][j] = (float) (kernel[i][j] * scale);
093                        }
094                }
095                return kernel;
096        }
097
098}