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.points; 010 011import imagingbook.common.geometry.basic.Pnt2d; 012import org.apache.commons.math3.linear.DecompositionSolver; 013import org.apache.commons.math3.linear.MatrixUtils; 014import org.apache.commons.math3.linear.RealMatrix; 015import org.apache.commons.math3.linear.RealVector; 016import org.apache.commons.math3.linear.SingularValueDecomposition; 017 018/** 019 * <p> 020 * This class implements 2D point fitting under projective (4-point) transformations (exact and least-squares). See Sec. 021 * 21.1.4 of [1] for details. 022 * </p> 023 * <p> 024 * [1] W. Burger, M.J. Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer 025 * (2022). 026 * </p> 027 * 028 * @author WB 029 * @version 2022/08/01 030 */ 031public class ProjectiveFit2d implements LinearFit2d { 032 033 private final RealMatrix A; // the calculated transformation matrix 034 private final double err; // the calculated error 035 036 /** 037 * Constructor. Fits two sequences of 2D points using a projective transformation model. At least 4 point pairs are 038 * required. For 4 point pairs, the solution is an exact fit, otherwise a least-squares fit is found. 039 * 040 * @param P the source points 041 * @param Q the target points 042 */ 043 public ProjectiveFit2d(Pnt2d[] P, Pnt2d[] Q) { 044 checkSize(P, Q); 045 final int n = Math.min(P.length, Q.length); 046 047 double[] ba = new double[2 * n]; 048 double[][] Ma = new double[2 * n][]; 049 for (int i = 0; i < n; i++) { 050 double px = P[i].getX(); 051 double py = P[i].getY(); 052 double qx = Q[i].getX(); 053 double qy = Q[i].getY(); 054 ba[2 * i + 0] = qx; 055 ba[2 * i + 1] = qy; 056 Ma[2 * i + 0] = new double[] { px, py, 1, 0, 0, 0, -qx * px, -qx * py }; 057 Ma[2 * i + 1] = new double[] { 0, 0, 0, px, py, 1, -qy * px, -qy * py }; 058 } 059 060 RealMatrix M = MatrixUtils.createRealMatrix(Ma); 061 RealVector b = MatrixUtils.createRealVector(ba); 062 DecompositionSolver solver = new SingularValueDecomposition(M).getSolver(); 063 RealVector h = solver.solve(b); 064 A = MatrixUtils.createRealMatrix(3, 3); 065 A.setEntry(0, 0, h.getEntry(0)); 066 A.setEntry(0, 1, h.getEntry(1)); 067 A.setEntry(0, 2, h.getEntry(2)); 068 A.setEntry(1, 0, h.getEntry(3)); 069 A.setEntry(1, 1, h.getEntry(4)); 070 A.setEntry(1, 2, h.getEntry(5)); 071 A.setEntry(2, 0, h.getEntry(6)); 072 A.setEntry(2, 1, h.getEntry(7)); 073 A.setEntry(2, 2, 1.0); 074 075 err = Math.sqrt(LinearFit2d.getSquaredError(P, Q, A.getData())); 076 } 077 078 @Override 079 public double[][] getTransformationMatrix() { 080 return A.getData(); 081 } 082 083 @Override 084 public double getError() { 085 return err; 086 } 087 088 // ------------------------------------------------------------------------------ 089 090 private void checkSize(Pnt2d[] P, Pnt2d[] Q) { 091 if (P.length < 4 || Q.length < 4) { 092 throw new IllegalArgumentException("4 or more point pairs are required to calculate this fit"); 093 } 094 } 095 096}