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-2023 Wilhelm Burger, Mark J. Burge. All rights reserved.
007 * Visit https://imagingbook.com for additional details.
008 ******************************************************************************/
009package imagingbook.common.image;
010
011import ij.process.ByteProcessor;
012
013
014/**
015 * <p>
016 * This class represents an 'integral image' or 'summed area table' as proposed in [1], See Sec. 2.8 of [2] for a
017 * detailed description. Currently only implemented for images of type {@link ByteProcessor}.
018 * </p>
019 * <p>
020 * [1] F. C. Crow. Summed-area tables for texture mapping. SIGGRAPH, Computer Graphics 18(3), 207–212 (1984).<br> [2] W.
021 * Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer (2022).
022 * </p>
023 *
024 * @author WB
025 * @version 2022/07/11
026 */
027public class IntegralImage {
028        
029        private final int M, N; 
030        private final long[][] S1, S2;
031
032        /**
033         * Constructor, creates a new {@link IntegralImage} instance from pixel values in a 2D int-array I[x][y].
034         *
035         * @param I pixel values
036         */
037        public IntegralImage(int[][] I) {
038                M = I.length;                   // image width
039                N = I[0].length;                // image height
040                S1 = new long[M][N];
041                S2 = new long[M][N];
042                
043                // initialize top-left corner (0,0)
044                S1[0][0] = I[0][0];
045                S2[0][0] = I[0][0] * I[0][0];
046                // do line v = 0:
047                for (int u = 1; u < M; u++) {
048                        S1[u][0] = S1[u-1][0] + I[u][0];
049                        S2[u][0] = S2[u-1][0] + I[u][0] * I[u][0];
050                }
051                // do lines v = 1,...,h-1
052                for (int v = 1; v < N; v++) {
053                        S1[0][v] = S1[0][v-1] + I[0][v];
054                        S2[0][v] = S2[0][v-1] + I[0][v] * I[0][v];
055                        for (int u = 1; u < M; u++) {
056                                S1[u][v] = S1[u-1][v] + S1[u][v-1] - S1[u-1][v-1] + I[u][v];
057                                S2[u][v] = S2[u-1][v] + S2[u][v-1] - S2[u-1][v-1] + I[u][v] * I[u][v];
058                        }
059                }
060        }
061
062        /**
063         * Constructor, creates a new {@link IntegralImage} instance from pixel values in the specified
064         * {@link ByteProcessor}.
065         *
066         * @param I input image
067         */
068        public IntegralImage(ByteProcessor I) {
069                this(I.getIntArray());
070        }
071        
072        // -------------------------------------------------------
073
074        /**
075         * Returns the summed area table of pixel values (&Sigma;<sub>1</sub>).
076         *
077         * @return array of &Sigma;<sub>1</sub> values
078         */
079        public long[][] getS1() {
080                return S1;
081        }
082
083        /**
084         * Returns the summed area table of squared pixel values (&Sigma;<sub>2</sub>).
085         *
086         * @return array of &Sigma;<sub>2</sub> values
087         */
088        public long[][] getS2() {
089                return S2;
090        }
091        
092        // -------------------------------------------------------
093
094        /**
095         * Calculates the sum of the pixel values in the rectangle R, specified by the corner points a = (ua, va) and b =
096         * (ub, vb).
097         *
098         * @param ua leftmost position in R
099         * @param va top position in R
100         * @param ub rightmost position in R
101         * @param vb bottom position in R
102         * @return the first-order block sum (S1(R)) inside the specified rectangle or zero if the rectangle is empty.
103         */
104        public long getBlockSum1(int ua, int va, int ub, int vb) {      // TODO: allow coordinates outside image boundaries
105                if (ub < ua || vb < va) {
106                        return 0;
107                }
108                final long saa = (ua > 0  && va > 0)  ? S1[ua - 1][va - 1] : 0;
109                final long sba = (ub >= 0 && va > 0)  ? S1[ub][va - 1] : 0;
110                final long sab = (ua > 0  && vb >= 0) ? S1[ua - 1][vb] : 0;
111                final long sbb = (ub >= 0 && vb >= 0) ? S1[ub][vb] : 0;
112                return sbb + saa - sba - sab;
113        }
114
115        /**
116         * Calculates the sum of the squared pixel values in the rectangle R, specified by the corner points a = (ua, va)
117         * and b = (ub, vb).
118         *
119         * @param ua leftmost position in R
120         * @param va top position in R
121         * @param ub rightmost position in R
122         * @param vb bottom position in R
123         * @return the second-order block sum (S2(R)) inside the specified rectangle or zero if the rectangle is empty.
124         */
125        public long getBlockSum2(int ua, int va, int ub, int vb) {      // TODO: allow coordinates outside image boundaries
126                if (ub < ua || vb < va) {
127                        return 0;
128                }
129                final long saa = (ua > 0  && va > 0)  ? S2[ua - 1][va - 1] : 0;
130                final long sba = (ub >= 0 && va > 0)  ? S2[ub][va - 1] : 0;
131                final long sab = (ua > 0  && vb >= 0) ? S2[ua - 1][vb] : 0;
132                final long sbb = (ub >= 0 && vb >= 0) ? S2[ub][vb] : 0;
133                return sbb + saa - sba - sab;
134        }
135
136        /**
137         * Returns the size of (number of pixels in) the specified rectangle.
138         *
139         * @param ua leftmost position in R
140         * @param va top position in R
141         * @param ub rightmost position in R {@literal (ub >= ua)}
142         * @param vb bottom position in R {@literal (vb >= va)}
143         * @return the size of the specified rectangle
144         */
145        public int getSize(int ua, int va, int ub, int vb) {
146                return (ub - ua + 1) * (vb - va + 1);
147        }
148
149        /**
150         * Calculates the mean of the image values in the specified rectangle.
151         *
152         * @param ua leftmost position in R
153         * @param va top position in R
154         * @param ub rightmost position in R {@literal (ub >= ua)}
155         * @param vb bottom position in R {@literal (vb >= va)}
156         * @return the mean value for the specified rectangle
157         */
158        public double getMean(int ua, int va, int ub, int vb) {
159                int size = getSize(ua, va, ub, vb);
160                if (size <= 0) {
161                        throw new IllegalArgumentException("region size must be positive");
162                }
163                return getBlockSum1(ua, va, ub, vb) / size;
164        }
165
166        /**
167         * Calculates the variance of the image values in the specified rectangle. parameters.
168         *
169         * @param ua leftmost position in R
170         * @param va top position in R
171         * @param ub rightmost position in R {@literal (u1 >= u0)}
172         * @param vb bottom position in R {@literal (v1 >= v0)}
173         * @return the variance for the specified rectangle
174         */
175        public double getVariance(int ua, int va, int ub, int vb) {
176                int size = getSize(ua, va, ub, vb);
177                if (size <= 0) {
178                        throw new IllegalArgumentException("region size must be positive");
179                }
180                double s1 = getBlockSum1(ua, va, ub, vb);
181                double s2 = getBlockSum2(ua, va, ub, vb);
182                return (s2 - (s1 * s1) / size) / size;
183        }
184
185}