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.circle.algebraic; 010 011import imagingbook.common.geometry.basic.Pnt2d; 012import imagingbook.common.geometry.circle.AlgebraicCircle; 013import imagingbook.common.geometry.circle.GeometricCircle; 014import org.apache.commons.math3.linear.MatrixUtils; 015import org.apache.commons.math3.linear.RealMatrix; 016 017import static imagingbook.common.math.Arithmetic.isZero; 018import static imagingbook.common.math.Arithmetic.sqr; 019 020/** 021 * Common interface for all algebraic circle fits. 022 * 023 * @author WB 024 * 025 */ 026public interface CircleFitAlgebraic { 027 028 public enum FitType { 029 KasaA, 030 KasaB, 031 KasaC, 032 Pratt, 033 Taubin, 034 HyperSimple, 035 HyperStable 036 } 037 038 /** 039 * Creates and returns a new circle fit instance of the specified type for the given sample points. 040 * 041 * @param type the circle fit type 042 * @param points an array of 2D sample points 043 * @return a new circle fit instance 044 * @see FitType 045 */ 046 public static CircleFitAlgebraic getFit(FitType type, Pnt2d[] points) { 047 switch (type) { 048 case KasaA: return new CircleFitKasaA(points); 049 case KasaB: return new CircleFitKasaB(points); 050 case KasaC: return new CircleFitKasaC(points); 051 case Pratt: return new CircleFitPratt(points); 052 case Taubin: return new CircleFitTaubin(points); 053 case HyperSimple: return new CircleFitHyperSimple(points); 054 case HyperStable: return new CircleFitHyperStable(points); 055 } 056 throw new IllegalArgumentException("unknown algebraic fit type: " + type); 057 } 058 059 /** 060 * Returns parameters (A, B, C, D) of the {@link AlgebraicCircle} or {@code null} if the fit was unsuccessful. 061 * Parameters are not normalized. 062 * 063 * @return the algebraic circle parameters or {@code null} 064 */ 065 public double[] getParameters(); 066 067 068 /** 069 * Returns a {@link AlgebraicCircle} instance for this fit or {@code null} if the fit was unsuccessful. 070 * 071 * @return a {@link AlgebraicCircle} instance or null 072 */ 073 public default AlgebraicCircle getAlgebraicCircle() { 074 double[] p = getParameters(); 075 return (p == null) ? null : new AlgebraicCircle(p); 076 } 077 078 /** 079 * Returns a {@link GeometricCircle} instance for this fit or {@code null} if the fit was unsuccessful. 080 * 081 * @return a {@link GeometricCircle} instance or null 082 */ 083 public default GeometricCircle getGeometricCircle() { 084 double[] q = this.getParameters(); // assumed to be (A, B, C, D) 085 if (q == null || isZero(q[0])) { // (abs(2 * A / s) < (1 / Rmax)) 086 return null; // return a straight line? 087 } 088 else { 089 return new GeometricCircle(getAlgebraicCircle()); 090 } 091 } 092 093 /** 094 * Used to transform the algebraic parameter vector for sample data shifted to the reference point specified. If qq 095 * is the parameter vector for the data centered at (xr, yr), the original parameters q are obtained as 096 * <pre>q = M * qq</pre> 097 * 098 * @param xr reference point (x) 099 * @param yr reference point (y) 100 * @return the transformation matrix 101 */ 102 static RealMatrix getDecenteringMatrix(double xr, double yr) { 103 return MatrixUtils.createRealMatrix(new double[][] 104 {{ 1, 0, 0, 0 }, 105 {-2*xr, 1, 0, 0 }, 106 {-2*yr, 0, 1, 0 }, 107 {sqr(xr) + sqr(yr), -xr, -yr, 1}}); 108 } 109 110 /** 111 * Normalize parameter vector q=(A,B,C,D) by enforcing the constraint B^2 + C^2 - 4 A D = 1. 112 * 113 * @param q original parameter vector 114 * @return normalized parameter vector 115 */ 116 public default double[] normalizeP(double[] q) { 117 final double A = q[0]; 118 final double B = q[1]; 119 final double C = q[2]; 120 final double D = q[3]; 121 double s = Math.sqrt(sqr(B) + sqr(C) - 4 * A * D); 122 return new double[] {A/s, B/s, C/s, D/s}; 123 } 124}