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.lib.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 ij.IJ; 019 020/** 021 * This class implements the Mahalanobis distance using the Apache commons math library. 022 * See the numerical example in Theodoridis/Koutumbras, "Pattern recognition", p. 27. 023 */ 024public class MahalanobisDistance extends VectorNorm { 025 026 public static final double DEFAULT_DIAGONAL_INCREMENT = 10E-6; 027 public static final boolean BIAS_CORRECTION = false; // we use no bias correction here (factor is 1/n) 028 029 private final int n; // number of samples 030 private final int m; // feature dimension 031 private final double[] mean; // the distribution's mean vector (\mu) 032 private final double[][] Cov; // covariance matrix of size n x n 033 private final double[][] iCov; // inverse covariance matrix size n x n 034 035 /** 036 * Create a new instance from an array of m-dimensional samples, e.g. 037 * samples = {{x1,y1}, {x2,y2},...,{xn,yn}} for a vector of n two-dimensional 038 * samples. 039 * @param samples A vector of length n with m-dimensional samples, i.e. 040 * samples[k][i] represents the i-th component of the k-th sample. 041 * @param diagonalIncrement Quantity added to the diagonal values of the 042 * covariance matrix to avoid singularity. 043 */ 044 public MahalanobisDistance(double[][] samples, double diagonalIncrement) { 045 n = samples.length; 046 m = samples[0].length; 047 if (n < 1 || m < 1) { 048 throw new IllegalArgumentException("dimension less than 1"); 049 } 050 if (diagonalIncrement < 0) { 051 throw new IllegalArgumentException("diagonal increment must be positive"); 052 } 053 054 mean = makeMeanVector(samples); 055 056 Covariance cov = new Covariance(samples, BIAS_CORRECTION); 057 RealMatrix S = cov.getCovarianceMatrix(); 058 059 // condition the covariance matrix to avoid singularity: 060 S = conditionCovarianceMatrix(S); 061 062 Cov = S.getData(); 063 // get the inverse covariance matrix 064 iCov = MatrixUtils.inverse(S).getData(); // not always symmetric? 065 } 066 067 /** 068 * Create a new instance from an array of m-dimensional samples, e.g. 069 * samples = {{x1,y1}, {x2,y2},...,{xn,yn}} for a vector of n two-dimensional 070 * samples. 071 * @param samples A vector of length n with m-dimensional samples, i.e. 072 * samples[k][i] represents the i-th component of the k-th sample. 073 */ 074 public MahalanobisDistance(double[][] samples) { 075 this(samples, DEFAULT_DIAGONAL_INCREMENT); 076 } 077 078 // 079 /** 080 * Conditions the supplied covariance matrix by enforcing 081 * positive eigenvalues along its main diagonal. 082 * @param S original covariance matrix 083 * @return modified covariance matrix 084 */ 085 private RealMatrix conditionCovarianceMatrix(RealMatrix S) { 086 EigenDecomposition ed = new EigenDecomposition(S); // S -> V . D . V^T 087 RealMatrix V = ed.getV(); 088 RealMatrix D = ed.getD(); // diagonal matrix of eigenvalues 089 RealMatrix VT = ed.getVT(); 090 for (int i = 0; i < D.getRowDimension(); i++) { 091 D.setEntry(i, i, Math.max(D.getEntry(i, i), 10E-6)); // setting eigenvalues to zero is not enough! 092 } 093 return V.multiply(D).multiply(VT); 094 } 095 096 // -------------------------------------------------------------- 097 098 private double[] makeMeanVector(double[][] samples) { 099 int n = samples.length; 100 int m = samples[0].length; 101 double[] mu = new double[m]; 102 for (int k = 0; k < n; k++) { 103 for (int i = 0; i < m; i++) { 104 mu[i] = mu[i] + samples[k][i]; 105 } 106 } 107 for (int i = 0; i < m; i++) { 108 mu[i] = mu[i] / n; 109 } 110 return mu; 111 } 112 113 public int getNumberOfSamples() { 114 return n; 115 } 116 117 public int getSampleDimension() { 118 return m; 119 } 120 121 //------------------------------------------------------------------------------------ 122 123 public double[][] getCovarianceMatrix() { 124 return Cov; 125 } 126 127 public double[][] getInverseCovarianceMatrix() { 128 return iCov; 129 } 130 131 public double[] getMeanVector() { 132 return mean; 133 } 134 135 //------------------------------------------------------------------------------------ 136 137 /** 138 * Returns the 'root' (U) of the inverse covariance matrix S^{-1}, 139 * such that S^{-1} = U^T . U 140 * This matrix can be used to pre-transform the original sample 141 * vectors X (by X → U . X) to a space where distance measurement (in the Mahalanobis 142 * sense) can be calculated with the usual Euclidean norm. 143 * The matrix U is invertible in case the reverse mapping is required. 144 * 145 * @return The matrix for pre-transforming the original sample vectors. 146 */ 147 public double[][] getWhiteningTransformation() { 148 IJ.log("in whitening"); 149 double relativeSymmetryThreshold = 1.0E-6; // CholeskyDecomposition.DEFAULT_RELATIVE_SYMMETRY_THRESHOLD == 1.0E-15; too small! 150 double absolutePositivityThreshold = 1.0E-10; // CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD == 1.0E-10; 151 CholeskyDecomposition cd = 152 new CholeskyDecomposition(MatrixUtils.createRealMatrix(iCov), 153 relativeSymmetryThreshold, absolutePositivityThreshold); 154 RealMatrix U = cd.getLT(); 155 return U.getData(); 156 } 157 158 //------------------------------------------------------------------------------------ 159 160 public double distance(double[] X, double[] Y) { 161 return Math.sqrt(distance2(X, Y)); 162 } 163 164 public double distance2(double[] X, double[] Y) { 165 if (X.length != m || Y.length != m) { 166 throw new IllegalArgumentException("vectors must be of length " + m); 167 } 168 // TODO: implement with methods from Matrix class 169 // d^2(X,Y) = (X-Y)^T . S^{-1} . (X-Y) 170 double[] dXY = new double[m]; // = (X-Y) 171 for (int k = 0; k < m; k++) { 172 dXY[k] = X[k] - Y[k]; 173 } 174 double[] iSdXY = new double[m]; // = S^{-1} . (X-Y) 175 for (int k = 0; k < m; k++) { 176 for (int i = 0; i < m; i++) { 177 iSdXY[k] += iCov[i][k] * dXY[i]; 178 } 179 } 180 double d2 = 0; // final sum 181 for (int k = 0; k < m; k++) { 182 d2 += dXY[k] * iSdXY[k]; 183 } // d = (X-Y)^T . S^{-1} . (X-Y) 184 return d2; 185 } 186 187 /** 188 * Calculates the Mahalanobis distance between point X and 189 * the mean of the reference distribution. 190 * @param X an arbitrary K-dimensional vector 191 * @return Mahalanobis distance between the point X and 192 * the mean of the reference distribution. 193 */ 194 public double distance(double[] X) { 195 return distance(X, mean); 196 } 197 198 /** 199 * Calculates the sqared Mahalanobis distance between point X and 200 * the mean of the reference distribution. 201 * @param X an arbitrary K-dimensional vector 202 * @return squared Mahalanobis distance between the point X and 203 * the mean of the reference distribution. 204 */ 205 public double distance2(double[] X) { 206 return distance2(X, mean); 207 } 208 209 //----------- unsupported (but required) stuff ------------------------------------------------- 210 211 public double distance(int[] a, int[] b) { 212 throw new UnsupportedOperationException(); 213 } 214 215 public double distance2(int[] a, int[] b) { 216 throw new UnsupportedOperationException(); 217 } 218 219 public double magnitude(double[] X) { 220 throw new UnsupportedOperationException(); 221 } 222 223 public double magnitude(int[] X) { 224 throw new UnsupportedOperationException(); 225 } 226 227 public double getScale(int n) { 228 return 1.0; 229 } 230 231 @Override 232 public double distance(float[] a, float[] b) { 233 throw new UnsupportedOperationException(); 234 } 235 236 @Override 237 public double distance2(float[] a, float[] b) { 238 throw new UnsupportedOperationException(); 239 } 240 241 // ---------------------------------------------------------------- 242 243 /** 244 * Test example from Burger/Burge UTICS-C Appendix: 245 * N = 4 samples, K = 3 dimensions 246 * 247 * @param args ignored 248 */ 249 public static void main(String[] args) { 250 251 double[] X0 = {0, 0, 0}; 252 double[] X1 = {75, 37, 12}; 253 double[] X2 = {41, 27, 20}; 254 double[] X3 = {93, 81, 11}; 255 double[] X4 = {12, 48, 52}; 256 double[] X5 = {12, 48, 100}; 257 258 259 System.out.println("X0 = " + Matrix.toString(X0)); 260 System.out.println("X1 = " + Matrix.toString(X1)); 261 System.out.println("X2 = " + Matrix.toString(X2)); 262 System.out.println("X3 = " + Matrix.toString(X3)); 263 System.out.println("X4 = " + Matrix.toString(X4)); 264 System.out.println("X5 = " + Matrix.toString(X5)); 265 System.out.println(); 266 267 double[][] samples = {X1, X2, X3, X4}; 268 269 MahalanobisDistance mhd = new MahalanobisDistance(samples); 270 271 double[] mu = mhd.getMeanVector(); 272 System.out.println("mu = " + Matrix.toString(mu)); 273 274 System.out.println(); 275 276 // covariance matrix Cov (3x3) 277 double[][] cov = mhd.getCovarianceMatrix(); 278 System.out.println("cov = " + Matrix.toString(cov)); 279 280 System.out.println(); 281 282 double[][] icov = mhd.getInverseCovarianceMatrix(); 283 System.out.println("icov = " + Matrix.toString(icov)); System.out.println(); 284 285 System.out.format("MH-dist(X1,X1) = %.3f\n", mhd.distance(X1, X1)); 286 System.out.format("MH-dist(X1,X2) = %.3f\n", mhd.distance(X1, X2)); 287 System.out.format("MH-dist(X2,X1) = %.3f\n", mhd.distance(X2, X1)); 288 System.out.format("MH-dist(X1,X3) = %.3f\n", mhd.distance(X1, X3)); 289 System.out.format("MH-dist(X1,X4) = %.3f\n", mhd.distance(X1, X4)); 290 System.out.format("MH-dist(X1,X5) = %.3f\n", mhd.distance(X1, X5)); 291 System.out.println(); 292 293 System.out.format("MH-dist(X1,mu) = %.3f\n", mhd.distance(X1, mu)); 294 System.out.format("MH-dist(X2,mu) = %.3f\n", mhd.distance(X2, mu)); 295 System.out.format("MH-dist(X3,mu) = %.3f\n", mhd.distance(X3, mu)); 296 System.out.format("MH-dist(X4,mu) = %.3f\n", mhd.distance(X4, mu)); 297 System.out.println(); 298 299 System.out.format("MH-dist(X5,mu) = %.3f\n", mhd.distance(X5, mu)); 300 System.out.format("MH-dist(X0,mu) = %.3f\n", mhd.distance(X0, mu)); 301 System.out.println(); 302 303 VectorNorm L2 = new VectorNorm.L2(); 304 System.out.format("L2-distance(X1,X2) = %.3f\n", L2.distance(X1, X2)); 305 System.out.format("L2-distance(X2,X1) = %.3f\n", L2.distance(X2, X1)); 306 307 // ------------------------------------------------------ 308 309 System.out.println(); 310 System.out.println("Testing pre-transformed Mahalanobis distances:"); 311 RealMatrix U = MatrixUtils.createRealMatrix(mhd.getWhiteningTransformation()); 312 System.out.println("U = " + Matrix.toString(U.getData())); 313 double[] Y0 = U.operate(X0); 314 double[] Y1 = U.operate(X1); 315 double[] Y2 = U.operate(X2); 316 double[] Y3 = U.operate(X3); 317 double[] Y4 = U.operate(X4); 318 double[] Y5 = U.operate(X5); 319 320 System.out.println("Y0 = " + Matrix.toString(Y0)); 321 System.out.println("Y1 = " + Matrix.toString(Y1)); 322 System.out.println("Y2 = " + Matrix.toString(Y2)); 323 System.out.println("Y3 = " + Matrix.toString(Y3)); 324 System.out.println("Y4 = " + Matrix.toString(Y4)); 325 System.out.println("Y5 = " + Matrix.toString(Y5)); 326 327 System.out.format("pre-transformed MH-distance(X1,X2) = %.3f\n", L2.distance(Y1, Y2)); 328 System.out.format("pre-transformed MH-distance(X1,X3) = %.3f\n", L2.distance(Y1, Y3)); 329 System.out.format("pre-transformed MH-distance(X1,X4) = %.3f\n", L2.distance(Y1, Y4)); 330 System.out.format("pre-transformed MH-distance(X1,X5) = %.3f\n", L2.distance(Y1, Y5)); 331 332 } 333 334/* Results: verified with Mathematica (bias-corrected): 335X0 = {0.000, 0.000, 0.000} 336X1 = {75.000, 37.000, 12.000} 337X2 = {41.000, 27.000, 20.000} 338X3 = {93.000, 81.000, 11.000} 339X4 = {12.000, 48.000, 52.000} 340X5 = {12.000, 48.000, 100.000} 341 342mu = {55.250, 48.250, 23.750} 343 344cov = {{1296.250, 442.583, -627.250}, 345{442.583, 550.250, -70.917}, 346{-627.250, -70.917, 370.917}} 347 348icov = {{0.024, -0.014, 0.038}, 349{-0.014, 0.011, -0.022}, 350{0.038, -0.022, 0.063}} 351 352MH-dist(X1,X1) = 0.000 353MH-dist(X1,X2) = 2.449 354MH-dist(X2,X1) = 2.449 355MH-dist(X1,X3) = 2.449 356MH-dist(X1,X4) = 2.449 357MH-dist(X1,X5) = 11.714 358 359MH-dist(X1,mu) = 1.500 360MH-dist(X2,mu) = 1.500 361MH-dist(X3,mu) = 1.500 362MH-dist(X4,mu) = 1.500 363 364MH-dist(X5,mu) = 12.613 365MH-dist(X0,mu) = 10.214 366 367L2-distance(X1,X2) = 36.332 368L2-distance(X2,X1) = 36.332 369 370Testing pre-transformed Mahalanobis distances: 371U = {{0.155, -0.093, 0.245}, 372{0.000, 0.043, 0.008}, 373{0.000, 0.000, 0.052}} 374pre-transformed MH-distance(X1,X2) = 2.449 375pre-transformed MH-distance(X1,X3) = 2.449 376pre-transformed MH-distance(X1,X4) = 2.449 377pre-transformed MH-distance(X1,X5) = 11.714 378*/ 379 380/* Results non-bias corrected: 381X0 = {0.000, 0.000, 0.000} 382X1 = {75.000, 37.000, 12.000} 383X2 = {41.000, 27.000, 20.000} 384X3 = {93.000, 81.000, 11.000} 385X4 = {12.000, 48.000, 52.000} 386X5 = {12.000, 48.000, 100.000} 387 388mu = {55.250, 48.250, 23.750} 389 390cov = {{972.188, 331.938, -470.438}, 391{331.938, 412.687, -53.188}, 392{-470.438, -53.188, 278.188}} 393 394icov = {{0.032, -0.019, 0.051}, 395{-0.019, 0.014, -0.030}, 396{0.051, -0.030, 0.083}} 397 398MH-dist(X1,X1) = 0.000 399MH-dist(X1,X2) = 2.828 400MH-dist(X2,X1) = 2.828 401MH-dist(X1,X3) = 2.828 402MH-dist(X1,X4) = 2.828 403MH-dist(X1,X5) = 13.527 404 405MH-dist(X1,mu) = 1.732 406MH-dist(X2,mu) = 1.732 407MH-dist(X3,mu) = 1.732 408MH-dist(X4,mu) = 1.732 409 410MH-dist(X5,mu) = 14.564 411MH-dist(X0,mu) = 11.794 412 413L2-distance(X1,X2) = 36.332 414L2-distance(X2,X1) = 36.332 415 416Testing pre-transformed Mahalanobis distances: 417U = {{0.179, -0.108, 0.282}, 418{0.000, 0.050, 0.010}, 419{0.000, 0.000, 0.060}} 420Y0 = {0.000, 0.000, 0.000} 421Y1 = {12.840, 1.959, 0.719} 422Y2 = {10.085, 1.536, 1.199} 423Y3 = {11.044, 4.142, 0.660} 424Y4 = {11.664, 2.888, 3.118} 425Y5 = {25.218, 3.345, 5.996} 426pre-transformed MH-distance(X1,X2) = 2.828 427pre-transformed MH-distance(X1,X3) = 2.828 428pre-transformed MH-distance(X1,X4) = 2.828 429pre-transformed MH-distance(X1,X5) = 13.527 430*/ 431}