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-2023 Wilhelm Burger, Mark J. Burge. All rights reserved. 007 * Visit https://imagingbook.com for additional details. 008 ******************************************************************************/ 009package imagingbook.common.geometry.fitting.line; 010 011import imagingbook.common.geometry.basic.Pnt2d; 012import imagingbook.common.math.Matrix; 013import org.apache.commons.math3.linear.Array2DRowRealMatrix; 014import org.apache.commons.math3.linear.MatrixUtils; 015import org.apache.commons.math3.linear.QRDecomposition; 016import org.apache.commons.math3.linear.RealMatrix; 017import org.apache.commons.math3.linear.SingularValueDecomposition; 018 019/** 020 * <p> 021 * This class implements line fitting by orthogonal regression to a set of 2D points by solving a homogeneous linear 022 * system. See Sec. 6.7.1 of [1] and Sec. 10.2.2 (Alg. 10.2) of [2] for additional details. 023 * </p> 024 * <p> 025 * [1] W. Gander, M. J. Gander, and F. Kwok. "Scientific Computing – An Introduction using Maple and MATLAB". Springer 026 * (2014). <br> [2] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd 027 * ed, Springer (2022). 028 * </p> 029 * 030 * @author WB 031 * @version 2022/09/29 032 * @see OrthogonalLineFitEigen 033 * @see OrthogonalLineFitSvd 034 */ 035public class OrthogonalLineFitHomogeneous implements LineFit { 036 037 private final Pnt2d[] points; 038 private final double[] p; // algebraic line parameters A,B,C 039 040 /** 041 * Constructor, performs a orthogonal regression fit to the specified points. At least two different points are 042 * required. 043 * 044 * @param points an array with at least 2 points 045 */ 046 public OrthogonalLineFitHomogeneous(Pnt2d[] points) { 047 if (points.length < 2) { 048 throw new IllegalArgumentException("line fit requires at least 2 points"); 049 } 050 this.points = points; 051 this.p = fit(); 052 } 053 054 @Override 055 public int getSize() { 056 return points.length; 057 } 058 059 @Override 060 public double[] getLineParameters() { 061 return p; 062 } 063 064 private double[] fit() { 065 final int n = points.length; //.max(3, pts.length); 066 final double[][] X = new double[n][3]; 067 068 for (int i = 0; i < points.length; i++) { 069 X[i][0] = 1; 070 X[i][1] = points[i].getX(); 071 X[i][2] = points[i].getY(); 072 } 073 074 QRDecomposition qr = new QRDecomposition(MatrixUtils.createRealMatrix(X)); 075 RealMatrix R = qr.getR(); 076 077// RealMatrix RR = R.getSubMatrix(1, 2, 1, 2); 078 RealMatrix RR = new Array2DRowRealMatrix(2,2); 079 RR.setEntry(0, 0, R.getEntry(1, 1)); 080 RR.setEntry(0, 1, R.getEntry(1, 2)); 081 if (R.getRowDimension() >= 3) { 082 RR.setEntry(1, 1, R.getEntry(2, 2)); 083 } 084 085 SingularValueDecomposition svd = new SingularValueDecomposition(RR); 086 087 RealMatrix V = svd.getV(); 088 double[] s = svd.getSingularValues(); 089 int k = Matrix.idxMin(s); 090 double[] e = V.getColumn(k); 091 092 double R01 = R.getEntry(0, 1); 093 double R02 = R.getEntry(0, 2); 094 double R00 = R.getEntry(0, 0); 095 096 double A = e[0]; 097 double B = e[1]; 098 double C = ( -A * R01 - B * R02) / R00; 099 return new double[] {A,B,C}; 100 } 101 102}