170e303692f3308905209adbc54cf298d441b03f4nicolasroard/* 270e303692f3308905209adbc54cf298d441b03f4nicolasroard * Copyright (C) 2013 The Android Open Source Project 370e303692f3308905209adbc54cf298d441b03f4nicolasroard * 470e303692f3308905209adbc54cf298d441b03f4nicolasroard * Licensed under the Apache License, Version 2.0 (the "License"); 570e303692f3308905209adbc54cf298d441b03f4nicolasroard * you may not use this file except in compliance with the License. 670e303692f3308905209adbc54cf298d441b03f4nicolasroard * You may obtain a copy of the License at 770e303692f3308905209adbc54cf298d441b03f4nicolasroard * 870e303692f3308905209adbc54cf298d441b03f4nicolasroard * http://www.apache.org/licenses/LICENSE-2.0 970e303692f3308905209adbc54cf298d441b03f4nicolasroard * 1070e303692f3308905209adbc54cf298d441b03f4nicolasroard * Unless required by applicable law or agreed to in writing, software 1170e303692f3308905209adbc54cf298d441b03f4nicolasroard * distributed under the License is distributed on an "AS IS" BASIS, 1270e303692f3308905209adbc54cf298d441b03f4nicolasroard * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 1370e303692f3308905209adbc54cf298d441b03f4nicolasroard * See the License for the specific language governing permissions and 1470e303692f3308905209adbc54cf298d441b03f4nicolasroard * limitations under the License. 1570e303692f3308905209adbc54cf298d441b03f4nicolasroard */ 1670e303692f3308905209adbc54cf298d441b03f4nicolasroard 1770e303692f3308905209adbc54cf298d441b03f4nicolasroardpackage com.android.gallery3d.filtershow.tools; 1870e303692f3308905209adbc54cf298d441b03f4nicolasroard 1970e303692f3308905209adbc54cf298d441b03f4nicolasroardimport android.util.Log; 2070e303692f3308905209adbc54cf298d441b03f4nicolasroard 2170e303692f3308905209adbc54cf298d441b03f4nicolasroardpublic class MatrixFit { 2270e303692f3308905209adbc54cf298d441b03f4nicolasroard // Simple implementation of a matrix fit in N dimensions. 2370e303692f3308905209adbc54cf298d441b03f4nicolasroard 2470e303692f3308905209adbc54cf298d441b03f4nicolasroard private static final String LOGTAG = "MatrixFit"; 2570e303692f3308905209adbc54cf298d441b03f4nicolasroard 2670e303692f3308905209adbc54cf298d441b03f4nicolasroard private double[][] mMatrix; 2770e303692f3308905209adbc54cf298d441b03f4nicolasroard private int mDimension; 2870e303692f3308905209adbc54cf298d441b03f4nicolasroard private boolean mValid = false; 2970e303692f3308905209adbc54cf298d441b03f4nicolasroard private static double sEPS = 1.0f/10000000000.0f; 3070e303692f3308905209adbc54cf298d441b03f4nicolasroard 3170e303692f3308905209adbc54cf298d441b03f4nicolasroard public MatrixFit(double[][] from, double[][] to) { 3270e303692f3308905209adbc54cf298d441b03f4nicolasroard mValid = fit(from, to); 3370e303692f3308905209adbc54cf298d441b03f4nicolasroard } 3470e303692f3308905209adbc54cf298d441b03f4nicolasroard 3570e303692f3308905209adbc54cf298d441b03f4nicolasroard public int getDimension() { 3670e303692f3308905209adbc54cf298d441b03f4nicolasroard return mDimension; 3770e303692f3308905209adbc54cf298d441b03f4nicolasroard } 3870e303692f3308905209adbc54cf298d441b03f4nicolasroard 3970e303692f3308905209adbc54cf298d441b03f4nicolasroard public boolean isValid() { 4070e303692f3308905209adbc54cf298d441b03f4nicolasroard return mValid; 4170e303692f3308905209adbc54cf298d441b03f4nicolasroard } 4270e303692f3308905209adbc54cf298d441b03f4nicolasroard 4370e303692f3308905209adbc54cf298d441b03f4nicolasroard public double[][] getMatrix() { 4470e303692f3308905209adbc54cf298d441b03f4nicolasroard return mMatrix; 4570e303692f3308905209adbc54cf298d441b03f4nicolasroard } 4670e303692f3308905209adbc54cf298d441b03f4nicolasroard 4770e303692f3308905209adbc54cf298d441b03f4nicolasroard public boolean fit(double[][] from, double[][] to) { 4870e303692f3308905209adbc54cf298d441b03f4nicolasroard if ((from.length != to.length) || (from.length < 1)) { 4970e303692f3308905209adbc54cf298d441b03f4nicolasroard Log.e(LOGTAG, "from and to must be of same size"); 5070e303692f3308905209adbc54cf298d441b03f4nicolasroard return false; 5170e303692f3308905209adbc54cf298d441b03f4nicolasroard } 5270e303692f3308905209adbc54cf298d441b03f4nicolasroard 5370e303692f3308905209adbc54cf298d441b03f4nicolasroard mDimension = from[0].length; 5470e303692f3308905209adbc54cf298d441b03f4nicolasroard mMatrix = new double[mDimension +1][mDimension + mDimension +1]; 5570e303692f3308905209adbc54cf298d441b03f4nicolasroard 5670e303692f3308905209adbc54cf298d441b03f4nicolasroard if (from.length < mDimension) { 5770e303692f3308905209adbc54cf298d441b03f4nicolasroard Log.e(LOGTAG, "Too few points => under-determined system"); 5870e303692f3308905209adbc54cf298d441b03f4nicolasroard return false; 5970e303692f3308905209adbc54cf298d441b03f4nicolasroard } 6070e303692f3308905209adbc54cf298d441b03f4nicolasroard 6170e303692f3308905209adbc54cf298d441b03f4nicolasroard double[][] q = new double[from.length][mDimension]; 6270e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < from.length; i++) { 6370e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int j = 0; j < mDimension; j++) { 6470e303692f3308905209adbc54cf298d441b03f4nicolasroard q[i][j] = from[i][j]; 6570e303692f3308905209adbc54cf298d441b03f4nicolasroard } 6670e303692f3308905209adbc54cf298d441b03f4nicolasroard } 6770e303692f3308905209adbc54cf298d441b03f4nicolasroard 6870e303692f3308905209adbc54cf298d441b03f4nicolasroard double[][] p = new double[to.length][mDimension]; 6970e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < to.length; i++) { 7070e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int j = 0; j < mDimension; j++) { 7170e303692f3308905209adbc54cf298d441b03f4nicolasroard p[i][j] = to[i][j]; 7270e303692f3308905209adbc54cf298d441b03f4nicolasroard } 7370e303692f3308905209adbc54cf298d441b03f4nicolasroard } 7470e303692f3308905209adbc54cf298d441b03f4nicolasroard 7570e303692f3308905209adbc54cf298d441b03f4nicolasroard // Make an empty (dim) x (dim + 1) matrix and fill it 7670e303692f3308905209adbc54cf298d441b03f4nicolasroard double[][] c = new double[mDimension+1][mDimension]; 7770e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int j = 0; j < mDimension; j++) { 7870e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int k = 0; k < mDimension + 1; k++) { 7970e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < q.length; i++) { 8070e303692f3308905209adbc54cf298d441b03f4nicolasroard double qt = 1; 8170e303692f3308905209adbc54cf298d441b03f4nicolasroard if (k < mDimension) { 8270e303692f3308905209adbc54cf298d441b03f4nicolasroard qt = q[i][k]; 8370e303692f3308905209adbc54cf298d441b03f4nicolasroard } 8470e303692f3308905209adbc54cf298d441b03f4nicolasroard c[k][j] += qt * p[i][j]; 8570e303692f3308905209adbc54cf298d441b03f4nicolasroard } 8670e303692f3308905209adbc54cf298d441b03f4nicolasroard } 8770e303692f3308905209adbc54cf298d441b03f4nicolasroard } 8870e303692f3308905209adbc54cf298d441b03f4nicolasroard 8970e303692f3308905209adbc54cf298d441b03f4nicolasroard // Make an empty (dim+1) x (dim+1) matrix and fill it 9070e303692f3308905209adbc54cf298d441b03f4nicolasroard double[][] Q = new double[mDimension+1][mDimension+1]; 9170e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int qi = 0; qi < q.length; qi++) { 9270e303692f3308905209adbc54cf298d441b03f4nicolasroard double[] qt = new double[mDimension + 1]; 9370e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < mDimension; i++) { 9470e303692f3308905209adbc54cf298d441b03f4nicolasroard qt[i] = q[qi][i]; 9570e303692f3308905209adbc54cf298d441b03f4nicolasroard } 9670e303692f3308905209adbc54cf298d441b03f4nicolasroard qt[mDimension] = 1; 9770e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < mDimension + 1; i++) { 9870e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int j = 0; j < mDimension + 1; j++) { 9970e303692f3308905209adbc54cf298d441b03f4nicolasroard Q[i][j] += qt[i] * qt[j]; 10070e303692f3308905209adbc54cf298d441b03f4nicolasroard } 10170e303692f3308905209adbc54cf298d441b03f4nicolasroard } 10270e303692f3308905209adbc54cf298d441b03f4nicolasroard } 10370e303692f3308905209adbc54cf298d441b03f4nicolasroard 10470e303692f3308905209adbc54cf298d441b03f4nicolasroard // Use a gaussian elimination to solve the linear system 10570e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < mDimension + 1; i++) { 10670e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int j = 0; j < mDimension + 1; j++) { 10770e303692f3308905209adbc54cf298d441b03f4nicolasroard mMatrix[i][j] = Q[i][j]; 10870e303692f3308905209adbc54cf298d441b03f4nicolasroard } 10970e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int j = 0; j < mDimension; j++) { 11070e303692f3308905209adbc54cf298d441b03f4nicolasroard mMatrix[i][mDimension + 1 + j] = c[i][j]; 11170e303692f3308905209adbc54cf298d441b03f4nicolasroard } 11270e303692f3308905209adbc54cf298d441b03f4nicolasroard } 11370e303692f3308905209adbc54cf298d441b03f4nicolasroard if (!gaussianElimination(mMatrix)) { 11470e303692f3308905209adbc54cf298d441b03f4nicolasroard return false; 11570e303692f3308905209adbc54cf298d441b03f4nicolasroard } 11670e303692f3308905209adbc54cf298d441b03f4nicolasroard return true; 11770e303692f3308905209adbc54cf298d441b03f4nicolasroard } 11870e303692f3308905209adbc54cf298d441b03f4nicolasroard 11970e303692f3308905209adbc54cf298d441b03f4nicolasroard public double[] apply(double[] point) { 12070e303692f3308905209adbc54cf298d441b03f4nicolasroard if (mDimension != point.length) { 12170e303692f3308905209adbc54cf298d441b03f4nicolasroard return null; 12270e303692f3308905209adbc54cf298d441b03f4nicolasroard } 12370e303692f3308905209adbc54cf298d441b03f4nicolasroard double[] res = new double[mDimension]; 12470e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int j = 0; j < mDimension; j++) { 12570e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < mDimension; i++) { 12670e303692f3308905209adbc54cf298d441b03f4nicolasroard res[j] += point[i] * mMatrix[i][j+ mDimension +1]; 12770e303692f3308905209adbc54cf298d441b03f4nicolasroard } 12870e303692f3308905209adbc54cf298d441b03f4nicolasroard res[j] += mMatrix[mDimension][j+ mDimension +1]; 12970e303692f3308905209adbc54cf298d441b03f4nicolasroard } 13070e303692f3308905209adbc54cf298d441b03f4nicolasroard return res; 13170e303692f3308905209adbc54cf298d441b03f4nicolasroard } 13270e303692f3308905209adbc54cf298d441b03f4nicolasroard 13370e303692f3308905209adbc54cf298d441b03f4nicolasroard public void printEquation() { 13470e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int j = 0; j < mDimension; j++) { 13570e303692f3308905209adbc54cf298d441b03f4nicolasroard String str = "x" + j + "' = "; 13670e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < mDimension; i++) { 13770e303692f3308905209adbc54cf298d441b03f4nicolasroard str += "x" + i + " * " + mMatrix[i][j+mDimension+1] + " + "; 13870e303692f3308905209adbc54cf298d441b03f4nicolasroard } 13970e303692f3308905209adbc54cf298d441b03f4nicolasroard str += mMatrix[mDimension][j+mDimension+1]; 14070e303692f3308905209adbc54cf298d441b03f4nicolasroard Log.v(LOGTAG, str); 14170e303692f3308905209adbc54cf298d441b03f4nicolasroard } 14270e303692f3308905209adbc54cf298d441b03f4nicolasroard } 14370e303692f3308905209adbc54cf298d441b03f4nicolasroard 14470e303692f3308905209adbc54cf298d441b03f4nicolasroard private void printMatrix(String name, double[][] matrix) { 14570e303692f3308905209adbc54cf298d441b03f4nicolasroard Log.v(LOGTAG, "name: " + name); 14670e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < matrix.length; i++) { 14770e303692f3308905209adbc54cf298d441b03f4nicolasroard String str = ""; 14870e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int j = 0; j < matrix[0].length; j++) { 14970e303692f3308905209adbc54cf298d441b03f4nicolasroard str += "" + matrix[i][j] + " "; 15070e303692f3308905209adbc54cf298d441b03f4nicolasroard } 15170e303692f3308905209adbc54cf298d441b03f4nicolasroard Log.v(LOGTAG, str); 15270e303692f3308905209adbc54cf298d441b03f4nicolasroard } 15370e303692f3308905209adbc54cf298d441b03f4nicolasroard } 15470e303692f3308905209adbc54cf298d441b03f4nicolasroard 15570e303692f3308905209adbc54cf298d441b03f4nicolasroard /* 15670e303692f3308905209adbc54cf298d441b03f4nicolasroard * Transforms the given matrix into a row echelon matrix 15770e303692f3308905209adbc54cf298d441b03f4nicolasroard */ 15870e303692f3308905209adbc54cf298d441b03f4nicolasroard private boolean gaussianElimination(double[][] m) { 15970e303692f3308905209adbc54cf298d441b03f4nicolasroard int h = m.length; 16070e303692f3308905209adbc54cf298d441b03f4nicolasroard int w = m[0].length; 16170e303692f3308905209adbc54cf298d441b03f4nicolasroard 16270e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int y = 0; y < h; y++) { 16370e303692f3308905209adbc54cf298d441b03f4nicolasroard int maxrow = y; 16470e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int y2 = y + 1; y2 < h; y2++) { // Find max pivot 16570e303692f3308905209adbc54cf298d441b03f4nicolasroard if (Math.abs(m[y2][y]) > Math.abs(m[maxrow][y])) { 16670e303692f3308905209adbc54cf298d441b03f4nicolasroard maxrow = y2; 16770e303692f3308905209adbc54cf298d441b03f4nicolasroard } 16870e303692f3308905209adbc54cf298d441b03f4nicolasroard } 16970e303692f3308905209adbc54cf298d441b03f4nicolasroard // swap 17070e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int i = 0; i < mDimension; i++) { 17170e303692f3308905209adbc54cf298d441b03f4nicolasroard double t = m[y][i]; 17270e303692f3308905209adbc54cf298d441b03f4nicolasroard m[y][i] = m[maxrow][i]; 17370e303692f3308905209adbc54cf298d441b03f4nicolasroard m[maxrow][i] = t; 17470e303692f3308905209adbc54cf298d441b03f4nicolasroard } 17570e303692f3308905209adbc54cf298d441b03f4nicolasroard 17670e303692f3308905209adbc54cf298d441b03f4nicolasroard if (Math.abs(m[y][y]) <= sEPS) { // Singular Matrix 17770e303692f3308905209adbc54cf298d441b03f4nicolasroard return false; 17870e303692f3308905209adbc54cf298d441b03f4nicolasroard } 17970e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int y2 = y + 1; y2 < h; y2++) { // Eliminate column y 18070e303692f3308905209adbc54cf298d441b03f4nicolasroard double c = m[y2][y] / m[y][y]; 18170e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int x = y; x < w; x++) { 18270e303692f3308905209adbc54cf298d441b03f4nicolasroard m[y2][x] -= m[y][x] * c; 18370e303692f3308905209adbc54cf298d441b03f4nicolasroard } 18470e303692f3308905209adbc54cf298d441b03f4nicolasroard } 18570e303692f3308905209adbc54cf298d441b03f4nicolasroard } 18670e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int y = h -1; y > -1; y--) { // Back substitution 18770e303692f3308905209adbc54cf298d441b03f4nicolasroard double c = m[y][y]; 18870e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int y2 = 0; y2 < y; y2++) { 18970e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int x = w - 1; x > y - 1; x--) { 19070e303692f3308905209adbc54cf298d441b03f4nicolasroard m[y2][x] -= m[y][x] * m[y2][y] / c; 19170e303692f3308905209adbc54cf298d441b03f4nicolasroard } 19270e303692f3308905209adbc54cf298d441b03f4nicolasroard } 19370e303692f3308905209adbc54cf298d441b03f4nicolasroard m[y][y] /= c; 19470e303692f3308905209adbc54cf298d441b03f4nicolasroard for (int x = h; x < w; x++) { // Normalize row y 19570e303692f3308905209adbc54cf298d441b03f4nicolasroard m[y][x] /= c; 19670e303692f3308905209adbc54cf298d441b03f4nicolasroard } 19770e303692f3308905209adbc54cf298d441b03f4nicolasroard } 19870e303692f3308905209adbc54cf298d441b03f4nicolasroard return true; 19970e303692f3308905209adbc54cf298d441b03f4nicolasroard } 20070e303692f3308905209adbc54cf298d441b03f4nicolasroard} 201