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 ij.IJ; 012import imagingbook.common.geometry.basic.Pnt2d; 013import imagingbook.common.geometry.basic.PntUtils; 014import imagingbook.common.math.Matrix; 015import imagingbook.common.math.eigen.EigenDecompositionJama; 016import org.apache.commons.math3.linear.MatrixUtils; 017import org.apache.commons.math3.linear.RealMatrix; 018import org.apache.commons.math3.linear.RealVector; 019 020import java.util.Arrays; 021 022import static imagingbook.common.math.Arithmetic.sqr; 023 024/** 025 * <p> 026 * Algebraic ellipse fit based on Fitzgibbon's method [1], numerically improved as suggested by Halir and Flusser [2]. 027 * See [3, Sec. 11.2.1] for a detailed description. Note: This implementation performs data centering or, alternatively, 028 * accepts a specific reference point. Capable of performing an (exact) 5-point fit! 029 * </p> 030 * <p> 031 * [1] A. W. Fitzgibbon, M. Pilu, and R. B. Fisher. Direct least- squares fitting of ellipses. IEEE Transactions on 032 * Pattern Analysis and Machine Intelligence 21(5), 476-480 (1999). <br> [2] R. Halíř and J. Flusser. Numerically stable 033 * direct least squares fitting of ellipses. In "Proceedings of the 6th International Conference in Central Europe on 034 * Computer Graphics and Visualization (WSCG’98)", pp. 125-132, Plzeň, CZ (February 1998). <br> [3] W. Burger, M.J. 035 * Burge, <em>Digital Image Processing – An Algorithmic Introduction</em>, 3rd ed, Springer (2022). 036 * </p> 037 * 038 * @author WB 039 * @version 2021/11/06 040 */ 041public class EllipseFitFitzgibbonStable implements EllipseFitAlgebraic { 042 043 private static final RealMatrix C1 = 044 Matrix.makeRealMatrix(3, 3, 045 0.0, 0.0, 2, 046 0.0, -1.0, 0.0, 047 2, 0.0, 0.0); 048 049 private static final RealMatrix C1i = 050 Matrix.makeRealMatrix(3, 3, 051 0.0, 0.0, 0.5, 052 0.0, -1.0, 0.0, 053 0.5, 0.0, 0.0); 054 055 private final double[] q; // = (A,B,C,D,E,F) ellipse parameters 056 057 public EllipseFitFitzgibbonStable(Pnt2d[] points, Pnt2d xref) { 058 this.q = fit(points, xref); 059 } 060 061 public EllipseFitFitzgibbonStable(Pnt2d[] points) { 062 this(points, PntUtils.centroid(points)); 063 } 064 065 @Override 066 public double[] getParameters() { 067 return this.q; 068 } 069 070 private double[] fit(Pnt2d[] points, Pnt2d xref) { 071 final int n = points.length; 072 if (n < 5) { 073 throw new IllegalArgumentException("fitter requires at least 5 sample points instead of " + points.length); 074 } 075 076 // reference point 077 final double xr = xref.getX(); 078 final double yr = xref.getY(); 079 080 RealMatrix X1 = MatrixUtils.createRealMatrix(n, 3); 081 RealMatrix X2 = MatrixUtils.createRealMatrix(n, 3); 082 083 for (int i = 0; i < n; i++) { 084 final double x = points[i].getX() - xr; // center data set 085 final double y = points[i].getY() - yr; 086 double[] f1 = {sqr(x), x*y, sqr(y)}; 087 double[] f2 = {x, y, 1}; 088 X1.setRow(i, f1); 089 X2.setRow(i, f2); 090 } 091 092 // build reduced scatter matrices: 093 RealMatrix S1 = X1.transpose().multiply(X1); 094 RealMatrix S2 = X1.transpose().multiply(X2); 095 RealMatrix S3 = X2.transpose().multiply(X2); 096 RealMatrix S3i = MatrixUtils.inverse(S3); 097 098 RealMatrix T = S3i.scalarMultiply(-1).multiply(S2.transpose()); 099 RealMatrix Z = C1i.multiply(S1.add(S2.multiply(T))); 100 101 // find the eigenvector of Z which satisfies the ellipse constraint: 102// EigenDecomposition ed = new EigenDecomposition(Z); 103 EigenDecompositionJama ed = new EigenDecompositionJama(Z); 104 double[] p1 = null; 105 for (int i = 0; i < 3; i++) { 106 RealVector e = ed.getEigenvector(i); 107 if (e.dotProduct(C1.operate(e)) > 0) { 108 p1 = e.toArray(); 109 break; 110 } 111 } 112 113 if (p1 == null) { 114 IJ.log("p1 is null! " + Arrays.toString(points)); 115 return null; 116 } 117 118 double[] p2 = T.operate(p1); 119 120 RealMatrix U = getDataOffsetCorrectionMatrix(xr, yr); 121 122 // assemble q 123 return U.operate(Matrix.join(p1, p2)); 124 } 125 126}