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.math;
011
012import org.apache.commons.math3.linear.CholeskyDecomposition;
013import org.apache.commons.math3.linear.EigenDecomposition;
014import org.apache.commons.math3.linear.MatrixUtils;
015import org.apache.commons.math3.linear.RealMatrix;
016import org.apache.commons.math3.stat.correlation.Covariance;
017
018import static imagingbook.common.math.Statistics.conditionCovarianceMatrix;
019import static imagingbook.common.math.Statistics.covarianceMatrix;
020
021/**
022 * <p>
023 * This class implements the Mahalanobis distance using the Apache Commons Math library. No statistical bias correction
024 * is used. See the Appendix G (Sections G.2-G.3) of [1] for additional details and examples.
025 * </p>
026 * <p>
027 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing &ndash; An Algorithmic Introduction</em>, 3rd ed, Springer
028 * (2022).
029 * </p>
030 *
031 * @author WB
032 * @version 2023/01/01
033 */
034public class MahalanobisDistance {
035
036        /**
037         * The default minimum diagonal value used to condition the covariance matrix.
038         */
039        public static final double DefaultMinimumDiagonalValue = 1.0E-15;
040
041        /**
042         * Tolerance used for validating symmetry in the Cholesky decomposition. See
043         * {@link CholeskyDecomposition#DEFAULT_RELATIVE_SYMMETRY_THRESHOLD} (1.0E-15), which seems too small).
044         */
045        public static final double DefaultRelativeSymmetryThreshold = 1.0E-6;
046
047        /**
048         * Tolerance used for validating positvity in the Cholesky decomposition. See
049         * {@link CholeskyDecomposition#DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD}.
050         */
051        public static final double DefaultAbsolutePositivityThreshold =
052                        CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD;
053
054        private final int m;                            // feature dimension (length of sample vectors)
055        private final double[][] cov;           // covariance matrix of size m x m
056        private final double[][] icov;          // inverse covariance matrix of size m x m
057
058        /**
059         * Constructor, creates a new {@link MahalanobisDistance} instance from the m x m covariance matrix and the
060         * m-dimensional mean vector of a distribution of m-dimensional data. The covariance matrix is assumed to be
061         * properly conditioned, i.e., has all-positive diagonal values.
062         *
063         * @param cov the covariance matrix of size m x m
064         */
065        public MahalanobisDistance(double[][] cov) {
066                this.m = cov.length;
067                if (!Matrix.isSquare(cov)) {
068                        throw new IllegalArgumentException("covariance matrix must be square");
069                }
070                if (!Matrix.isSymmetric(cov)) {
071                        throw new IllegalArgumentException("covariance matrix must be symmetric");
072                }
073                this.cov = cov;
074                this.icov = Matrix.inverse(cov);
075        }
076
077        /**
078         * Creates a new {@link MahalanobisDistance} instance from an array of m-dimensional samples, e.g., samples = {{x1,y1},
079         * {x2,y2},...,{xn,yn}} for a vector of n two-dimensional samples. Uses {@link #DefaultMinimumDiagonalValue} to
080         * condition the covariance matrix.
081         *
082         * @param samples a sequence of m-dimensional samples, i.e., samples[k][i] represents the i-th component of the k-th
083         * sample
084         * @return a {@link MahalanobisDistance} instance
085         */
086        public static MahalanobisDistance fromSamples(double[][] samples) {
087                return fromSamples(samples, DefaultMinimumDiagonalValue);
088        }
089
090        /**
091         * Creates a new {@link MahalanobisDistance} instance from an array of m-dimensional samples, e.g. samples =
092         * {{x1,y1}, {x2,y2},...,{xn,yn}} for a vector of two-dimensional samples (n = 2).
093         *
094         * @param samples a vector of length n with m-dimensional samples, i.e., samples[k][i] represents the i-th component
095         * of the k-th sample
096         * @param minDiagVal the minimum diagonal value used to condition the covariance matrix to avoid singularity (see
097         * {@link #DefaultMinimumDiagonalValue})
098         * @return a {@link MahalanobisDistance} instance
099         */
100        public static MahalanobisDistance fromSamples(double[][] samples, double minDiagVal) {
101                if (samples.length < 2)
102                        throw new IllegalArgumentException("number of samples must be 2 or more");
103                if (samples[0].length < 1)
104                        throw new IllegalArgumentException("sample dimension must be at least 1");
105                double[][] cov = covarianceMatrix(samples);
106                return new MahalanobisDistance(conditionCovarianceMatrix(cov, minDiagVal));
107        }
108
109        //------------------------------------------------------------------------------------
110
111        /**
112         * Returns the sample dimension (M) for this instance.
113         * @return the sample dimension
114         */
115        public int getDimension() {
116                return m;
117        }
118
119        /**
120         * Returns the covariance matrix used for distance calculations.
121         * @return the conditioned covariance matrix
122         */
123        public double[][] getCovarianceMatrix() {
124                return cov;
125        }
126
127        /**
128         * Returns the inverse of the conditioned covariance matrix.
129         * @return the inverse of the conditioned covariance matrix
130         */
131        public double[][] getInverseCovarianceMatrix() {
132                return icov;
133        }
134
135        //------------------------------------------------------------------------------------
136
137        /**
138         * Returns the 'root' (U) of the inverse covariance matrix S^{-1}, such that S^{-1} = U^T . U This matrix can be
139         * used to pre-transform the original sample vectors X (by X &#8594; U . X) to a space where distance (in the
140         * Mahalanobis sense) can be measured with the usual Euclidean norm. The matrix U is invertible in case the reverse
141         * transformation is required.
142         *
143         * @return the m x m matrix for pre-transforming the original (m-dimensional) sample vectors
144         */
145        public double[][] getWhiteningTransformation() {
146                return getWhiteningTransformation(DefaultRelativeSymmetryThreshold, DefaultAbsolutePositivityThreshold);
147        }
148
149        /**
150         * Returns the 'root' (U) of the inverse covariance matrix S^{-1}, such that S^{-1} = U^T . U This matrix can be *
151         * used to pre-transform the original sample vectors X (by X &#8594; U . X) to a space where distance (in the *
152         * Mahalanobis sense) can be measured with the usual Euclidean norm. The matrix U is invertible in case the reverse
153         * * transformation is required.
154         *
155         * @param relativeSymmetryThreshold maximum deviation from symmetry (see {@link #DefaultRelativeSymmetryThreshold})
156         * @param absolutePositivityThreshold maximum deviation from positivity (see {@link #DefaultAbsolutePositivityThreshold})
157         * @return
158         */
159        public double[][] getWhiteningTransformation(double relativeSymmetryThreshold, double absolutePositivityThreshold) {
160                CholeskyDecomposition cd =
161                                new CholeskyDecomposition(MatrixUtils.createRealMatrix(icov),
162                                                relativeSymmetryThreshold, absolutePositivityThreshold);
163                RealMatrix U = cd.getLT();
164                return U.getData();
165        }
166
167        //------------------------------------------------------------------------------------
168
169        /**
170         * Returns the Mahalanobis distance between the given points.
171         * @param a first point
172         * @param b second point
173         * @return the Mahalanobis distance
174         */
175        public double distance(double[] a, double[] b) {
176                return Math.sqrt(distance2(a, b));
177        }
178
179        /**
180         * Returns the squared Mahalanobis distance between the given points a, b: d^2(a,b) = (a-b)^T . S^{-1} . (a-b)
181         *
182         * @param a first point
183         * @param b second point
184         * @return the squared Mahalanobis distance
185         */
186        public double distance2(double[] a, double[] b) {
187                if (a.length != m || b.length != m) {
188                        throw new IllegalArgumentException("vectors a, b must be of length " + m);
189                }       
190                double[] amb = Matrix.subtract(a, b);   // amb = a - b
191                return Matrix.dotProduct(Matrix.multiply(amb, icov), amb);
192                
193        }
194        
195        // original version
196//      public double distance2(double[] X, double[] Y) {
197//              if (X.length != m || Y.length != m) {
198//                      throw new IllegalArgumentException("vectors must be of length " + m);
199//              }
200//              // d^2(X,Y) = (X-Y)^T . S^{-1} . (X-Y)
201//              double[] dXY = new double[m]; // = (X-Y)
202//              for (int k = 0; k < m; k++) {
203//                      dXY[k] = X[k] - Y[k];
204//              }
205//              double[] iSdXY = new double[m]; // = S^{-1} . (X-Y)
206//              for (int k = 0; k < m; k++) {
207//                      for (int i = 0; i < m; i++) {
208//                              iSdXY[k] += icov[i][k] * dXY[i];
209//                      }
210//              }
211//              double d2 = 0;  // final sum
212//              for (int k = 0; k < m; k++) {
213//                      d2 += dXY[k] * iSdXY[k];
214//              }                                               // d = (X-Y)^T . S^{-1} . (X-Y)
215//              return d2;
216//      }
217
218}