1 /*
2  * Copyright (C) 2013 The Android Open Source Project
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  *      http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 package com.android.gallery3d.filtershow.tools;
18 
19 import android.util.Log;
20 
21 public class MatrixFit {
22     // Simple implementation of a matrix fit in N dimensions.
23 
24     private static final String LOGTAG = "MatrixFit";
25 
26     private double[][] mMatrix;
27     private int mDimension;
28     private boolean mValid = false;
29     private static double sEPS = 1.0f/10000000000.0f;
30 
MatrixFit(double[][] from, double[][] to)31     public MatrixFit(double[][] from, double[][] to) {
32         mValid = fit(from, to);
33     }
34 
getDimension()35     public int getDimension() {
36         return mDimension;
37     }
38 
isValid()39     public boolean isValid() {
40         return mValid;
41     }
42 
getMatrix()43     public double[][] getMatrix() {
44         return mMatrix;
45     }
46 
fit(double[][] from, double[][] to)47     public boolean fit(double[][] from, double[][] to) {
48         if ((from.length != to.length) || (from.length < 1)) {
49             Log.e(LOGTAG, "from and to must be of same size");
50             return false;
51         }
52 
53         mDimension = from[0].length;
54         mMatrix = new double[mDimension +1][mDimension + mDimension +1];
55 
56         if (from.length < mDimension) {
57             Log.e(LOGTAG, "Too few points => under-determined system");
58             return false;
59         }
60 
61         double[][] q = new double[from.length][mDimension];
62         for (int i = 0; i < from.length; i++) {
63             for (int j = 0; j < mDimension; j++) {
64                 q[i][j] = from[i][j];
65             }
66         }
67 
68         double[][] p = new double[to.length][mDimension];
69         for (int i = 0; i < to.length; i++) {
70             for (int j = 0; j < mDimension; j++) {
71                 p[i][j] = to[i][j];
72             }
73         }
74 
75         // Make an empty (dim) x (dim + 1) matrix and fill it
76         double[][] c = new double[mDimension+1][mDimension];
77         for (int j = 0; j < mDimension; j++) {
78             for (int k = 0; k < mDimension + 1; k++) {
79                 for (int i = 0; i < q.length; i++) {
80                     double qt = 1;
81                     if (k < mDimension) {
82                         qt = q[i][k];
83                     }
84                     c[k][j] += qt * p[i][j];
85                 }
86             }
87         }
88 
89         // Make an empty (dim+1) x (dim+1) matrix and fill it
90         double[][] Q = new double[mDimension+1][mDimension+1];
91         for (int qi = 0; qi < q.length; qi++) {
92             double[] qt = new double[mDimension + 1];
93             for (int i = 0; i < mDimension; i++) {
94                 qt[i] = q[qi][i];
95             }
96             qt[mDimension] = 1;
97             for (int i = 0; i < mDimension + 1; i++) {
98                 for (int j = 0; j < mDimension + 1; j++) {
99                     Q[i][j] += qt[i] * qt[j];
100                 }
101             }
102         }
103 
104         // Use a gaussian elimination to solve the linear system
105         for (int i = 0; i < mDimension + 1; i++) {
106             for (int j = 0; j < mDimension + 1; j++) {
107                 mMatrix[i][j] = Q[i][j];
108             }
109             for (int j = 0; j < mDimension; j++) {
110                 mMatrix[i][mDimension + 1 + j] = c[i][j];
111             }
112         }
113         if (!gaussianElimination(mMatrix)) {
114             return false;
115         }
116         return true;
117     }
118 
apply(double[] point)119     public double[] apply(double[] point) {
120         if (mDimension != point.length) {
121             return null;
122         }
123         double[] res = new double[mDimension];
124         for (int j = 0; j < mDimension; j++) {
125             for (int i = 0; i < mDimension; i++) {
126                 res[j] += point[i] * mMatrix[i][j+ mDimension +1];
127             }
128             res[j] += mMatrix[mDimension][j+ mDimension +1];
129         }
130         return res;
131     }
132 
printEquation()133     public void printEquation() {
134         for (int j = 0; j < mDimension; j++) {
135             String str = "x" + j + "' = ";
136             for (int i = 0; i < mDimension; i++) {
137                 str += "x" + i + " * " + mMatrix[i][j+mDimension+1] + " + ";
138             }
139             str += mMatrix[mDimension][j+mDimension+1];
140             Log.v(LOGTAG, str);
141         }
142     }
143 
printMatrix(String name, double[][] matrix)144     private void printMatrix(String name, double[][] matrix) {
145         Log.v(LOGTAG, "name: " + name);
146         for (int i = 0; i < matrix.length; i++) {
147             String str = "";
148             for (int j = 0; j < matrix[0].length; j++) {
149                 str += "" + matrix[i][j] + " ";
150             }
151             Log.v(LOGTAG, str);
152         }
153     }
154 
155     /*
156      * Transforms the given matrix into a row echelon matrix
157      */
gaussianElimination(double[][] m)158     private boolean gaussianElimination(double[][] m) {
159         int h = m.length;
160         int w = m[0].length;
161 
162         for (int y = 0; y < h; y++) {
163             int maxrow = y;
164             for (int y2 = y + 1; y2 < h; y2++) { // Find max pivot
165                 if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y])) {
166                     maxrow = y2;
167                 }
168             }
169             // swap
170             for (int i = 0; i < mDimension; i++) {
171                 double t = m[y][i];
172                 m[y][i] = m[maxrow][i];
173                 m[maxrow][i] = t;
174             }
175 
176             if (Math.abs(m[y][y]) <= sEPS) { // Singular Matrix
177                 return false;
178             }
179             for (int y2 = y + 1; y2 < h; y2++) { // Eliminate column y
180                 double c = m[y2][y] / m[y][y];
181                 for (int x = y; x < w; x++) {
182                     m[y2][x] -= m[y][x] * c;
183                 }
184             }
185         }
186         for (int y = h -1; y > -1; y--) { // Back substitution
187             double c = m[y][y];
188             for (int y2 = 0; y2 < y; y2++) {
189                 for (int x = w - 1; x > y - 1; x--) {
190                     m[y2][x] -= m[y][x] * m[y2][y] / c;
191                 }
192             }
193             m[y][y] /= c;
194             for (int x = h; x < w; x++) { // Normalize row y
195                 m[y][x] /= c;
196             }
197         }
198         return true;
199     }
200 }
201