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 ******************************************************************************/
009
010package imagingbook.common.image.matching;
011
012import ij.process.ImageProcessor;
013
014import static imagingbook.common.math.Arithmetic.sqr;
015import static java.lang.Math.sqrt;
016
017/**
018 * <p>
019 * Instances of this class perform matching on scalar-valued images based on the correlation coefficient. The "search"
020 * image I (to be searched for matches of the "reference" image R) is initially associated with the
021 * {@link CorrCoeffMatcher}. The assumption is, that the search image I is fixed and the matcher tries to match multiple
022 * reference images R. See Sec. 23.1.1 (Alg. 23.1) of [1] for additional details.
023 * </p>
024 * <p>
025 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
026 * (2022).
027 * </p>
028 *
029 * @author WB
030 * @version 2022/09/16
031 */
032public class CorrCoeffMatcher {
033
034        private final float[][] fI;             // search image (copied to a float array)
035        private final int wI, hI;               // width/height of image
036
037        private float[][] fR;                   // reference image (copied to a float array)
038        private int wR, hR;                     // width/height of reference
039
040        // modified variables:
041        private int K;                                  // number of pixels in reference image
042        private float meanR;                    // mean value of reference
043        private float varR;                     // square root of reference variance
044
045        /**
046         * Constructor, accepts an {@link ImageProcessor} instance. Color images are converted to grayscale.
047         * @param I the search image (to be matched to)
048         */
049        public CorrCoeffMatcher(ImageProcessor I) {
050                this(I.getFloatArray());
051        }
052
053        /**
054         * Constructor, accepts a 2D float array.
055         * @param fI the search image (to be matched to)
056         */
057        public CorrCoeffMatcher(float[][] fI) {
058                this.fI = fI;
059                this.wI = fI.length;
060                this.hI = fI[0].length;
061        }
062
063        /**
064         * Matches the specified reference image R to the (fixed) search image I. Resulting score values are in [-1, 1], the
065         * score for the optimal match is +1. The returned score matrix has the size of the search image I reduced by the
066         * size of the reference image R.
067         *
068         * @param R a scalar-valued reference image
069         * @return a 2D array Q[r][s] of match scores in [-1,1]
070         */
071        public float[][] getMatch(ImageProcessor R) {
072                return getMatch(R.getFloatArray());
073        }
074
075        /**
076         * Matches the specified reference image R to the (fixed) search image I. Resulting score values are in [-1, 1], the
077         * score for the optimal match is +1. The returned score matrix has the size of the search image I reduced by the
078         * size of the reference image R.
079         *
080         * @param fR a scalar-valued reference image
081         * @return a 2D array Q[r][s] of match scores in [-1,1]
082         */
083        public float[][] getMatch(float[][] fR) {
084                this.fR = fR;
085                this.wR = fR.length;
086                this.hR = fR[0].length;
087                this.K = wR * hR;
088
089                // calculate the mean and variance of R
090                float sumR = 0;
091                float sumR2 = 0;
092                for (int j = 0; j < hR; j++) {
093                        for (int i = 0; i < wR; i++) {
094                                float b = fR[i][j];
095                                sumR  += b;
096                                sumR2 += sqr(b);
097                        }
098                }
099                
100                this.meanR = sumR / K;
101                this.varR = (float) sqrt(sumR2 - K * sqr(meanR));
102                
103                float[][] C = new float[wI - wR + 1][hI - hR + 1];
104                for (int r = 0; r < C.length; r++) {
105                        for (int s = 0; s < C[r].length; s++) {
106                                C[r][s] = getMatchValue(r, s);
107                        }       
108                }
109                this.fR = null;
110                return C;
111        }
112        
113        private float getMatchValue(int r, int s) {
114                float sumI = 0, sumI2 = 0, sumIR = 0;
115                for (int j = 0; j < hR; j++) {
116                        for (int i = 0; i < wR; i++) {
117                                float a = fI[r + i][s + j];
118                                float b = fR[i][j];
119                                sumI  += a;
120                                sumI2 += sqr(a);
121                                sumIR += a * b;
122                        }
123                }
124                float meanI = sumI / K;
125                return (sumIR - K * meanI * meanR) / 
126                           (1f + (float) sqrt(sumI2 - K * sqr(meanI)) * varR);
127                        // added 1 in denominator to handle flat image regions (w. zero variance)
128        }  
129        
130}