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.matching;
011
012import ij.process.FloatProcessor;
013
014/**
015 * @author W. Burger
016 * @version 2014-04-20
017 */
018public class CorrCoeffMatcher {
019        private final FloatProcessor I; // search image
020        private final int MI, NI;               // width/height of image
021
022        private FloatProcessor R;               // reference image
023        private int MR, NR;                     // width/height of reference
024        private int K;
025        private double meanR;                   // mean value of template
026        private double varR;                    // square root of template variance
027
028        public CorrCoeffMatcher(FloatProcessor I) {
029                this.I = I;
030                this.MI = this.I.getWidth();
031                this.NI = this.I.getHeight();
032        }
033        
034        public float[][] getMatch(FloatProcessor R) {
035                this.R = R;
036                this.MR = R.getWidth();
037                this.NR = R.getHeight();
038                this.K = MR * NR;
039
040                // calculate the mean and variance of template
041                double sumR = 0;
042                double sumR2 = 0;
043                for (int j = 0; j < NR; j++) {
044                        for (int i = 0; i < MR; i++) {
045                                float aR = R.getf(i,j);
046                                sumR  += aR;
047                                sumR2 += aR * aR;
048                        }
049                }
050                
051                this.meanR = sumR / K;
052                this.varR = Math.sqrt(sumR2 - K * meanR * meanR);
053                
054                float[][] C = new float[MI - MR + 1][NI - NR + 1];
055                for (int r = 0; r <= MI - MR; r++) {
056                        for (int s = 0; s <= NI - NR; s++) {
057                                float d = (float) getMatchValue(r, s);
058                                C[r][s] = d;
059                        }       
060                }
061                this.R = null;
062                return C;
063        }
064        
065        private double getMatchValue(int r, int s) {
066                double sumI = 0, sumI2 = 0, sumIR = 0;
067                for (int j = 0; j < NR; j++) {
068                        for (int i = 0; i < MR; i++) {
069                                float aI = I.getf(r + i, s + j);
070                                float aR = R.getf(i, j);
071                                sumI  += aI;
072                                sumI2 += aI * aI;
073                                sumIR += aI * aR;
074                        }
075                }
076                double meanI = sumI / K;
077                return (sumIR - K * meanI * meanR) / 
078                           (1 + Math.sqrt(sumI2 - K * meanI * meanI) * varR);
079                                // WB: added 1 in denominator to handle flat image regions (w. zero variance)
080        }  
081        
082}