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;
014
015/**
016 * This class implements a 5x5 Nagao-Matsuyama filter, as described in
017 * NagaoMatsuyama (1979).
018 * 
019 * @author W. Burger
020 * @version 2013/05/30
021 */
022
023public class NagaoMatsuyamaFilter extends GenericFilter {
024        
025        public static class Parameters {
026                /** Variance threshold */
027                public double varThreshold = 0.0;       // 0,...,10
028        }
029        
030        private final Parameters params;
031        
032        private static final int[][] R1 =
033        {{-1,-1}, {0,-1}, {1,-1},
034         {-1, 0}, {0, 0}, {1, 0},
035         {-1, 1}, {0, 1}, {1, 1}};
036        
037        private static final int[][] R2 =
038        {{-2,-1}, {-1,-1},
039         {-2, 0}, {-1, 0}, {0, 0},
040         {-2, 1}, {-1, 1}};
041        
042        private static final int[][] R3 =
043        {{-2,-2}, {-1,-2},
044         {-2,-1}, {-1,-1}, {0,-1},
045                  {-1, 0}, {0, 0}};
046        
047        private static final int[][] R4 =
048        {{-1,-2}, {0,-2}, {1,-2}, 
049         {-1,-1}, {0,-1}, {1,-1},
050                  {0, 0}};
051
052        private static final int[][] R5 =
053        {        {1,-2}, {2,-2},
054         {0,-1}, {1,-1}, {2,-1},
055         {0, 0}, {1, 0}};
056        
057        private static final int[][] R6 =
058        {        {1,-1}, {2,-1},
059         {0, 0}, {1, 0}, {2, 0},
060                 {1, 1}, {2, 1}};
061        
062        private static final int[][] R7 =
063        {{0,0}, {1,0},
064         {0,1}, {1,1}, {2,1},
065                {1,2}, {2,2}};
066
067        private static final int[][] R8 =
068        {        {0,0},
069         {-1,1}, {0,1}, {1,1},
070         {-1,2}, {0,2}, {1,2}};
071        
072        private static final int[][] R9 =
073        {        {-1,0}, {0,0},
074         {-2,1}, {-1,1}, {0,1},
075         {-2,2}, {-1,2}};
076        
077        private static final int[][][] subRegions =
078                {R2, R3, R4, R5, R6, R7, R8, R9};
079        
080        // ------------------------------------------------------
081        
082        public NagaoMatsuyamaFilter() {
083                this(new Parameters());
084        }
085        
086        public NagaoMatsuyamaFilter(Parameters params) {
087                this.params = params;
088        }
089        
090        private float minVariance;
091        private float minMean;
092        private float minMeanR;
093        private float minMeanG;
094        private float minMeanB;
095        
096        // ------------------------------------------------------
097
098        public float filterPixel(ImageAccessor.Scalar image, int u, int v) {
099                minVariance = Float.MAX_VALUE;
100                evalSubregion(image, R1, u, v);
101                minVariance = minVariance - (float) params.varThreshold;
102                for (int[][] Rk : subRegions) {
103                        evalSubregion(image, Rk, u, v);
104                }
105                return minMean;
106        }
107        
108        void evalSubregion(ImageAccessor.Scalar ia, int[][] R, int u, int v) {
109                float sum1 = 0; 
110                float sum2 = 0;
111                int n = 0;
112                for (int[] p : R) {
113                        float a = ia.getVal(u+p[0], v+p[1]);
114                        sum1 = sum1 + a;
115                        sum2 = sum2 + a * a;
116                        n = n + 1;
117                }
118                float nr = n;
119                float var = (sum2 - sum1 * sum1 / nr) / nr;     // = sigma^2
120                if (var < minVariance) {
121                        minVariance = var;
122                        minMean = sum1 / nr; // mean
123                }
124        }
125        
126        // ------------------------------------------------------
127        
128        final float[] rgb = {0,0,0};
129        
130        public float[] filterPixel(ImageAccessor.Rgb ia, int u, int v) {
131                minVariance = Float.MAX_VALUE;
132                evalSubregionColor(ia, R1, u, v);
133                minVariance = minVariance - (3 * (float) params.varThreshold);
134                for (int[][] Rk : subRegions) {
135                        evalSubregionColor(ia, Rk, u, v);
136                }
137                rgb[0] = (int) Math.rint(minMeanR);
138                rgb[1] = (int) Math.rint(minMeanG);
139                rgb[2] = (int) Math.rint(minMeanB);
140                return rgb;
141        }
142        
143        void evalSubregionColor(ImageAccessor.Rgb ia, int[][] R, int u, int v) {
144                //final int[] cpix = {0,0,0};
145                int sum1R = 0; int sum2R = 0;
146                int sum1G = 0; int sum2G = 0;
147                int sum1B = 0; int sum2B = 0;
148                int n = 0;
149                for (int[] p : R) {
150                        final float[] cpix = ia.getPix(u + p[0], v + p[1]);
151                        int red = (int) cpix[0];
152                        int grn = (int) cpix[1];
153                        int blu = (int) cpix[2];
154                        sum1R = sum1R + red;
155                        sum1G = sum1G + grn;
156                        sum1B = sum1B + blu;
157                        sum2R = sum2R + red * red;
158                        sum2G = sum2G + grn * grn;
159                        sum2B = sum2B + blu * blu;
160                        n = n + 1;
161                }
162                float nr = n;
163                // calculate variance for this subregion:
164                float varR = (sum2R - sum1R * sum1R / nr) / nr;
165                float varG = (sum2G - sum1G * sum1G / nr) / nr;
166                float varB = (sum2B - sum1B * sum1B / nr) / nr;
167                // total variance:
168                float totalVar = varR + varG + varB;    
169                if (totalVar < minVariance) {
170                        minVariance = totalVar;
171                        minMeanR = sum1R / nr;
172                        minMeanG = sum1G / nr;
173                        minMeanB = sum1B / nr;
174                }
175        }
176
177}