OrthogonalRegressionPlane.java

/*-
 *  Copyright 2020 Beuth Hochschule für Technik Berlin, Hochschule für Technik Stuttgart
 * 
 *  This file is part of CityDoctor2.
 *
 *  CityDoctor2 is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  CityDoctor2 is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with CityDoctor2.  If not, see <https://www.gnu.org/licenses/>.
 */
package de.hft.stuttgart.citydoctor2.math;

import java.util.List;

import Jama.EigenvalueDecomposition;
import Jama.Matrix;

/**
 * Utility class to calculate regression planes
 * 
 * @author Matthias Betz
 *
 */
public class OrthogonalRegressionPlane {

	private OrthogonalRegressionPlane() {

	}

	/**
	 * calculates a regression plane using the PCA algorithm
	 * 
	 * @param vertices the vertices forming a plane
	 * @return the regression plane
	 */
	public static Plane calculateOrthogonalRegressionPlane(List<? extends Vector3d> vertices) {
		// get the center of all points, the plane has to go through this point
		Vector3d centroid = CovarianceMatrix.getCentroid(vertices);
		EigenvalueDecomposition ed = decompose(vertices, centroid);
		Vector3d eigenvalues = getEigenvalues(ed);
		return calculatePlane(centroid, ed, eigenvalues);
	}

	public static Plane calculatePlane(Vector3d centroid, EigenvalueDecomposition ed, Vector3d eigenvalues) {
		double eig0 = eigenvalues.getX();
		double eig1 = eigenvalues.getY();
		double eig2 = eigenvalues.getZ();
		Matrix v = ed.getV();
		Matrix eigenVector;
		if (eig0 <= eig1 && eig0 <= eig2) {
			// the first eigenvalue is the lowest
			eigenVector = v.getMatrix(0, 2, 0, 0);
		} else if (eig1 <= eig0 && eig1 <= eig2) {
			// the second eigenvalue is the lowest
			eigenVector = v.getMatrix(0, 2, 1, 1);
		} else {
			// the third eigenvalue is the lowest
			eigenVector = v.getMatrix(0, 2, 2, 2);
		}
		// the eigenvector of the lowest eigenvalue is the normal of the plane
		Vector3d normal = new Vector3d(eigenVector.get(0, 0), eigenVector.get(1, 0), eigenVector.get(2, 0));

		return new Plane(normal, centroid);
	}

	public static Vector3d getEigenvalues(EigenvalueDecomposition ed) {
		Matrix eigenValues = ed.getD();
		// get the lowest eigenvalue and the corresponding eigenvector
		double eig0 = eigenValues.get(0, 0);
		double eig1 = eigenValues.get(1, 1);
		double eig2 = eigenValues.get(2, 2);
		return new Vector3d(eig0, eig1, eig2);
	}

	public static EigenvalueDecomposition decompose(List<? extends Vector3d> vertices, Vector3d centroid) {
		// create covariance matrix (3 x 3) out of all used points
		double[][] covMatrixValues = CovarianceMatrix.calculateCovarianceMatrix(vertices, centroid);
		Matrix covMatrix = new Matrix(covMatrixValues);
		// calculate eigenvalues and eigenvectors
		return covMatrix.eig();
	}

}