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.ellipse.algebraic; 010 011import imagingbook.common.geometry.basic.Pnt2d; 012import imagingbook.common.geometry.basic.PntUtils; 013import imagingbook.common.math.Matrix; 014import imagingbook.common.math.eigen.GeneralizedSymmetricEigenDecomposition; 015import org.apache.commons.math3.linear.MatrixUtils; 016import org.apache.commons.math3.linear.RealMatrix; 017 018import static imagingbook.common.math.Arithmetic.sqr; 019 020/** 021 * <p> 022 * Algebraic ellipse fit based on Taubin's method [1]. Version 1 uses the full scatter and constraint matrix (6x6), the 023 * solution is found by a generalized symmetric eigendecomposition. Note that the constraint matrix (C) is not positive 024 * definite. See [3, Sec. 11.2.1] for a detailed description (Alg. 11.7). This implementation performs data centering 025 * or, alternatively, accepts a specific reference point. 026 * <p> 027 * [1] G. Taubin, G. Taubin. "Estimation of planar curves, surfaces, and nonplanar space curves defined by implicit 028 * equations with applications to edge and range image segmentation", IEEE Transactions on Pattern Analysis and Machine 029 * Intelligence 13(11), 1115–1138 (1991). <br> [2] W. Burger, M.J. Burge, <em>Digital Image Processing – An 030 * Algorithmic Introduction</em>, 3rd ed, Springer (2022). 031 * </p> 032 * 033 * @author WB 034 * @version 2021/11/09 035 * @see EllipseFitTaubin2 036 */ 037public class EllipseFitTaubin1 implements EllipseFitAlgebraic { 038 039 private final double[] q; // = (A,B,C,D,E,F) ellipse parameters 040 041 public EllipseFitTaubin1(Pnt2d[] points, Pnt2d xref) { 042 if (points.length < 5) { 043 throw new IllegalArgumentException("at least 5 points must be supplied for ellipse fitting"); 044 } 045 this.q = fit(points, xref); 046 047 } 048 049 public EllipseFitTaubin1(Pnt2d[] points) { 050 this(points, PntUtils.centroid(points)); 051 } 052 053 @Override 054 public double[] getParameters() { 055 return this.q; 056 } 057 058 // -------------------------------------------------------------------------- 059 060 private double[] fit(Pnt2d[] points, Pnt2d xref) { 061 final int n = points.length; 062 063 // reference point 064 final double xr = xref.getX(); 065 final double yr = xref.getY(); 066 067 // design matrix (same as Fitzgibbon): 068 RealMatrix X = MatrixUtils.createRealMatrix(n, 6); 069 for (int i = 0; i < n; i++) { 070 final double x = points[i].getX() - xr; // center data set 071 final double y = points[i].getY() - yr; 072 double[] fi = {sqr(x), x*y, sqr(y), x, y, 1}; 073 X.setRow(i, fi); 074 } 075 076 // scatter matrix S (normalized by 1/n) 077 RealMatrix S = X.transpose().multiply(X).scalarMultiply(1.0/n); 078 079 double a = S.getEntry(0, 5); 080 double b = S.getEntry(2, 5); 081 double c = S.getEntry(1, 5); //2*s[1][5]; 082 double d = S.getEntry(3, 5); 083 double e = S.getEntry(4, 5); 084 085 // constraint matrix 086 RealMatrix C = Matrix.makeRealMatrix(6, 6, 087 4*a, 2*c, 0, 2*d, 0, 0 , 088 2*c, a+b, 2*c, e, d, 0 , 089 0, 2*c, 4*b, 0, 2*e, 0 , 090 2*d, e, 0, 1, 0, 0 , 091 0, d, 2*e, 0, 1, 0 , 092 0, 0, 0, 0, 0, 0 ); 093 094 // solve C*p = lambda*S*p (note: C is not positive definite!) 095 GeneralizedSymmetricEigenDecomposition eigen = 096 new GeneralizedSymmetricEigenDecomposition(C, S, 1e-15, 1e-15); 097 098 double[] evals = eigen.getRealEigenvalues(); 099 100 int k = Matrix.idxMax(evals); // index of largest eigenvalue 101 double[] qopt = eigen.getEigenvector(k).toArray(); 102 103 RealMatrix U = getDataOffsetCorrectionMatrix(xr, yr); 104 return U.operate(qopt); 105 } 106}