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 *******************************************************************************/
009package imagingbook.lib.image;
010
011import ij.IJ;
012import ij.process.ByteProcessor;
013
014
015/**
016 * This class represents an 'integral image' or 'summed area table' [Crow, 1984], 
017 * as described in the book (see 2nd English ed. 2016, Sec. 3.8).
018 * Currently only implemented for images of type {@link ByteProcessor}.
019 * 
020 * @author W. Burger
021 * @version 2015/11/15
022 *
023 */
024public class IntegralImage {
025        
026        private final int M, N;
027        private final long[][] S1, S2;
028        
029        /**
030         * Creates a new integral image from pixel values in a 2D int-array.
031         * @param I pixel values
032         */
033        public IntegralImage(int[][] I) {
034                M = I.length;
035                N = I[0].length;
036                S1 = new long[M][N];
037                S2 = new long[M][N];
038                
039                // initialize top-left corner (0,0)
040                S1[0][0] = I[0][0];
041                S2[0][0] = I[0][0] * I[0][0];
042                // do line v = 0:
043                for (int u = 1; u < M; u++) {
044                        S1[u][0] = S1[u-1][0] + I[u][0];
045                        S2[u][0] = S2[u-1][0] + I[u][0] * I[u][0];
046                }
047                // do lines v = 1,...,h-1
048                for (int v = 1; v < N; v++) {
049                        S1[0][v] = S1[0][v-1] + I[0][v];
050                        S2[0][v] = S2[0][v-1] + I[0][v] * I[0][v];
051                        for (int u = 1; u < M; u++) {
052                                S1[u][v] = S1[u-1][v] + S1[u][v-1] - S1[u-1][v-1] + I[u][v];
053                                S2[u][v] = S2[u-1][v] + S2[u][v-1] - S2[u-1][v-1] + I[u][v] * I[u][v];
054                        }
055                }
056        }
057        
058        /**
059         * Creates a new integral image from pixel values in a {@link ByteProcessor}.
060         * @param I input image
061         */
062        public IntegralImage(ByteProcessor I) {
063                this(I.getIntArray());
064        }
065        
066        
067        // -------------------------------------------------------
068        
069        /**
070         * Returns the summed area table of pixel values (Sigma_1).
071         * @return Array of Sigma_1 values
072         */
073        public long[][] getS1() {
074                return S1;
075        }
076        
077        /**
078         * Returns the summed area table of squared pixel values (Sigma_2).
079         * @return Array of Sigma_2 values
080         */
081        public long[][] getS2() {
082                return S2;
083        }
084        
085        // -------------------------------------------------------
086        
087        /**
088         * Calculates the sum of the pixel values in the rectangle
089         * R, specified by the corner points a = (ua, va) and b = (b1, vb).
090         * @param ua leftmost position in R
091         * @param va top position in R
092         * @param ub rightmost position in R
093         * @param vb bottom position in R
094         * @return the first-order block sum (S1(R)) inside the specified rectangle 
095         * or zero if the rectangle is empty.
096         */
097        public long getBlockSum1(int ua, int va, int ub, int vb) {
098                if (ub < ua || vb < va) 
099                        return 0;
100                final long saa = (ua > 0  && va > 0)  ? S1[ua - 1][va - 1] : 0;
101                final long sba = (ub >= 0 && va > 0)  ? S1[ub][va - 1] : 0;
102                final long sab = (ua > 0  && vb >= 0) ? S1[ua - 1][vb] : 0;
103                final long sbb = (ub >= 0 && vb >= 0) ? S1[ub][vb] : 0;
104                return sbb + saa - sba - sab;
105        }
106        
107        /**
108         * Calculates the sum of the squared pixel values in the rectangle
109         * R, specified by the corner points a = (ua, va) and b = (b1, vb).
110         * @param ua leftmost position in R
111         * @param va top position in R
112         * @param ub rightmost position in R
113         * @param vb bottom position in R
114         * @return the second-order block sum (S2(R)) inside the specified rectangle
115         * or zero if the rectangle is empty.
116         */
117        public long getBlockSum2(int ua, int va, int ub, int vb) {
118                if (ub < ua || vb < va) 
119                        return 0;
120                final long saa = (ua > 0  && va > 0)  ? S2[ua - 1][va - 1] : 0;
121                final long sba = (ub >= 0 && va > 0)  ? S2[ub][va - 1] : 0;
122                final long sab = (ua > 0  && vb >= 0) ? S2[ua - 1][vb] : 0;
123                final long sbb = (ub >= 0 && vb >= 0) ? S2[ub][vb] : 0;
124                return sbb + saa - sba - sab;
125        }
126        
127        public int getSize(int u0, int v0, int u1, int v1) {
128                return (1 + u1 - u0) * (1 + v1 - v0);
129        }
130        
131        /**
132         * Calculates the mean of the image values in the specified rectangle.
133         * 
134         * @param ua leftmost position in R
135         * @param va top position in R
136         * @param ub rightmost position in R {@literal (u1 >= u0)}
137         * @param vb bottom position in R {@literal (v1 >= v0)}
138         * @return the mean value for the specified rectangle
139         */
140        public double getMean(int ua, int va, int ub, int vb) {
141                int N = getSize(ua, va, ub, vb);
142                if (N <= 0)
143                        throw new IllegalArgumentException("region size must be positive");
144                double S1 = getBlockSum1(ua, va, ub, vb);
145                IJ.log("u0 = " + ua); IJ.log("v0 = " + va);
146                IJ.log("u1 = " + ub); IJ.log("v1 = " + vb);
147                IJ.log("S1 = " + S1);
148                return S1 / N;
149        }
150        
151        /**
152         * Calculates the variance of the image values in the specified rectangle.
153         * parameters.
154         * @param ua leftmost position in R
155         * @param va top position in R
156         * @param ub rightmost position in R {@literal (u1 >= u0)}
157         * @param vb bottom position in R {@literal (v1 >= v0)}
158         * @return the variance for the specified rectangle
159         */
160        public double getVariance(int ua, int va, int ub, int vb) {
161                int N = getSize(ua, va, ub, vb);
162                if (N <= 0)
163                        throw new IllegalArgumentException("region size must be positive");
164                double S1 = getBlockSum1(ua, va, ub, vb);
165                double S2 = getBlockSum2(ua, va, ub, vb);
166                return (S2 - (S1 * S1) / N) / N;
167        }
168
169}