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.edgepreservingfilters;
011
012import imagingbook.lib.filters.GenericFilter;
013import imagingbook.lib.image.ImageAccessor;
014import imagingbook.lib.math.Arithmetic;
015import imagingbook.lib.math.VectorNorm;
016import imagingbook.lib.math.VectorNorm.NormType;
017
018/**
019 * This class implements a bilateral filter as proposed in
020 * C. Tomasi and R. Manduchi, "Bilateral Filtering for Gray and Color Images",
021 * Proceedings of the 1998 IEEE International Conference on Computer Vision,
022 * Bombay, India.
023 * The filter uses Gaussian domain and range kernels and can be applied to all 
024 * image types.
025 * @author W. Burger
026 * @version 2013/05/30
027 */
028public class BilateralFilter extends GenericFilter {
029        
030        public static class Parameters {
031                /** Sigma (width) of domain filter */
032                public double sigmaD = 2;               
033                /** Sigma (width) of range filter */
034                public double sigmaR = 50;
035                /** Specify which distance norm to use */
036                public NormType colorNormType = NormType.L2;
037                
038                /**
039                 * Create a default parameter object.
040                 * @param sigmaD Sigma (width) of domain filter
041                 * @param sigmaR Sigma (width) of range filter
042                 * @return a new parameter object
043                 */
044                static Parameters create(double sigmaD, double sigmaR) {
045                        Parameters p = new Parameters();
046                        p.sigmaD = sigmaD;
047                        p.sigmaR = sigmaR;
048                        return p;
049                }
050        }
051        
052        protected final Parameters params;
053        
054        private float[][] Hd;   // the domain kernel
055        protected final int K;
056        protected final float[] rgb = {0,0,0};
057        protected final double sigmaR2;
058        protected final VectorNorm colorNorm;
059        protected final double colorScale;
060        
061        public BilateralFilter() {
062                this(new Parameters());
063        }
064        
065        // only for convenience / book compatibility:
066        public BilateralFilter(double sigmaD, double sigmaR) {
067                this(Parameters.create(sigmaD, sigmaR));
068        }
069        
070        public BilateralFilter(Parameters params) {
071                this.params = params;
072                K = (int) Math.max(1, 3.5 * params.sigmaD);
073                sigmaR2 = params.sigmaR * params.sigmaR;
074                colorNorm = params.colorNormType.create();
075                colorScale = Arithmetic.sqr(colorNorm.getScale(3));
076                initialize();
077        }
078        
079        private void initialize() {
080                Hd = makeDomainKernel2D(params.sigmaD, K);
081        }
082        
083        public float filterPixel(ImageAccessor.Scalar I, int u, int v) {
084                float S = 0;                    // sum of weighted pixel values
085                float W = 0;                    // sum of weights
086                
087                float a = I.getVal(u, v); // value of the current center pixel
088                
089                for (int m = -K; m <= K; m++) {
090                        for (int n = -K; n <= K; n++) {
091                                float b = I.getVal(u + m, v + n);
092                                float wd = Hd[m + K][n + K];
093                                float wr = similarityGauss(a, b);
094                                float w = wd * wr;
095                                S = S + w * b;
096                                W = W + w;
097                        }
098                }
099                return S / W;
100        }
101        
102        public float[] filterPixel(ImageAccessor.Rgb I, int u, int v) {
103                float[] S = new float[3];       // sum of weighted RGB values
104                float W = 0;                            // sum of weights
105                //int[] a = new int[3];
106                //int[] b = new int[3];
107                
108                float[] a = I.getPix(u, v);                     // value of the current center pixel
109                
110                for (int m = -K; m <= K; m++) {
111                        for (int n = -K; n <= K; n++) {
112                                float[] b = I.getPix(u + m, v + n);
113                                float wd = Hd[m + K][n + K];
114                                float wr = similarityGauss(a, b);
115                                float w = wd * wr;
116                                S[0] = S[0] + w * b[0];
117                                S[1] = S[1] + w * b[1];
118                                S[2] = S[2] + w * b[2];
119                                W = W + w;
120                        }
121                }
122                rgb[0] = Math.round(S[0] / W);
123                rgb[1] = Math.round(S[1] / W);
124                rgb[2] = Math.round(S[2] / W);
125                return rgb;
126        }
127        
128        // ------------------------------------------------------
129        // This returns the weights for a Gaussian range kernel (scalar version):
130        protected float similarityGauss(float a, float b) {
131                double dI = a - b;
132                return (float) Math.exp(-(dI * dI) / (2 * sigmaR2));
133        }
134        
135        // This returns the weights for a Gaussian range kernel (color vector version):
136        protected float similarityGauss(float[] a, float[] b) {
137                double d2 = colorScale * colorNorm.distance2(a, b);
138                return (float) Math.exp(-d2 / (2 * sigmaR2));
139        }
140        
141        // Color distances need to be scaled to yield the same range of
142        // values. This method returns the scale factor for squared 
143        // distances, since this is what we use in the Gaussian.
144//      static double getColorDistanceScale(NormType norm) {
145//              double s = 1.0;
146//              switch (norm) {
147//              case L1 : s = 1/3.0; break;                             // L1-dist is in [0,...,3*255*3]
148//              case L2:  s = Math.sqrt(1/3.0); break;  // L2-dist is in [0,...,sqrt(3*255^2)] = [0,...,sqrt(3)*255]
149//              case Linf: s = 1.0; break;                              // Linf-dist is in [0,...,255]
150//              default: break;
151//              }
152//              return s * s;   // scale factor for the squared distance
153//      }
154
155        // ------------------------------------------------------
156
157        protected float[][] makeDomainKernel2D(double sigma, int K) {
158                int size = K + 1 + K;
159                float[][] domainKernel = new float[size][size]; //center cell = kernel[K][K]
160                double sigma2 = sigma * sigma;
161                double scale = 1.0 / (2 * Math.PI * sigma2);
162                for (int i = 0; i < size; i++) {
163                        double x = K - i;
164                        for (int j = 0; j < size; j++) {
165                                double y = K - j;
166                                domainKernel[i][j] =  (float) (scale * Math.exp(-0.5 * (x*x + y*y) / sigma2));
167                        }
168                }
169                return domainKernel;
170        }
171        
172        protected float[] makeRangeKernel(double sigma, int K) {
173                int size = K + 1 + K;
174                float[] rangeKernel = new float[size]; //center cell = kernel[K]
175                double sigma2 = sigma * sigma;
176                for (int i = 0; i < size; i++) {
177                        double x = K - i;
178                        rangeKernel[i] =  (float) Math.exp(-0.5 * (x*x) / sigma2);
179                }
180                return rangeKernel;
181        }
182
183}