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
012
013/**
014 * The this class computes the eigenvalues/eigenvectors for a 2x2
015 * matrix
016 *   | A B |
017 *   | C D |
018 * The implementation was inspired by Blinn, Jim: Consider the lowly 2x2 matrix. 
019 * IEEE Computer Graphics and Applications, 16(2):82-88, 1996.
020 * @author W. Burger
021 * @version 2013-07-29
022 */
023public class Eigensolver2x2 {
024        
025        private final double A, B, C, D;
026        
027        public Eigensolver2x2(double[][] A) {
028                this(A[0][0], A[0][1], A[1][0], A[1][1]);
029        }
030        
031        public Eigensolver2x2(double A, double B, double C, double D) {
032                this.A = A;
033                this.B = B;
034                this.C = C;
035                this.D = D;
036        }
037
038        public EigenPair[] realEigenValues2x2() {
039                final double R = (A + D) / 2;
040                final double S = (A - D) / 2;
041                if ((S * S + B * C) < 0) // no real eigenvalues
042                        return null;
043                else {
044                        double T = Math.sqrt(S * S + B * C);
045                        double eVal1 = R + T;
046                        double eVal2 = R - T;
047                        double[] eVec1, eVec2;
048                        if ((A - D) >= 0) {
049                                eVec1 = new double[] {S + T, C};
050                                eVec2 = new double[] {B, -S - T};
051                        } else {
052                                eVec1 = new double[] {B, -S + T};
053                                eVec2 = new double[] {S - T, C};
054                        }
055
056                        EigenPair[] e = new EigenPair[2];
057                        // put eigenpair with larger eigenvalue up front
058                        if (Math.abs(eVal1) >= Math.abs(eVal2)) {
059                                e[0] = new EigenPair(eVal1, eVec1);
060                                e[1] = new EigenPair(eVal2, eVec2);
061                        } else {
062                                e[1] = new EigenPair(eVal1, eVec1);
063                                e[0] = new EigenPair(eVal2, eVec2);
064                        }
065                        return e;
066                }
067        }
068
069        /**
070         * EigenPair is a tuple (eigenvalue, eigenvector) and represents
071         * the solution to an eigen problem.
072         */
073        public static class EigenPair {
074                final double eival;
075                final double[] eivec;
076
077                public EigenPair(double eival, double[] eivec) {
078                        this.eival = eival;
079                        this.eivec = eivec.clone();
080                }
081
082                public double getEigenvalue() {
083                        return this.eival;
084                }
085
086                public double[] getEigenvector() {
087                        return this.eivec;
088                }
089                
090                public String toString() {
091                        if (eivec == null)
092                                return "no eigenvalue / eigenvector";
093                        else {
094                                return String.format("eigenvalue: %.5f | eigenvector: %s", eival, Matrix.toString(eivec)) ;
095                        }
096                }
097        }
098        
099        
100        // for Testing:
101//      public static void main(String[] args) {
102//              
103//              double[][] A = {
104//                              {-0.009562, 0.011933}, 
105//                              {0.011933, -0.021158}
106//                              };
107//              /*
108//               * eigenvalue: -0.02863 | eigenvector: {0.012, -0.019}
109//               * eigenvalue: -0.00209 | eigenvector: {0.019, 0.012}
110//               */
111//              
112//              double[][] B = {
113//                              {-0.004710, -0.006970},
114//                              {-0.006970, -0.029195}};
115//              /*
116//               * eigenvalue: -0.03104 | eigenvector: {-0.007, -0.026}
117//               * eigenvalue: -0.00286 | eigenvector: {0.026, -0.007}
118//               */
119//              
120//              EigenPair[] eigenvals; 
121//              eigenvals = new Eigensolver2x2(A).realEigenValues2x2();
122//              System.out.println(eigenvals[0].toString());
123//              System.out.println(eigenvals[1].toString());
124//              
125//              eigenvals = new Eigensolver2x2(B).realEigenValues2x2();
126//              System.out.println(eigenvals[0].toString());
127//              System.out.println(eigenvals[1].toString());
128//      }
129        
130
131}