1 /*
2 * ThreePointCircleFitter.java
3 *
4 * Created on March 23, 2006, 3:41 PM
5 *
6 *
7 */
8
9 package org.lcsim.fit.threepointcircle;
10
11 import static java.lang.Math.pow;
12 import static java.lang.Math.sqrt;
13 /**
14 * A class to return the circle which passes through three points.
15 * @author Norman Graf
16 */
17 public class ThreePointCircleFitter
18 {
19 private CircleFit _fit;
20 private static double[][] a = {{0, 0, 0},{0, 0, 0},{0, 0, 0}};
21 private static double[][] m = {{0, 0, 0},{0, 0, 0},{0, 0, 0}};
22
23 /** Creates a new instance of ThreePointCircleFitter */
24 public ThreePointCircleFitter()
25 {
26 }
27
28 /**
29 * Find the circle which passes through three points
30 * @param p1 The first point.
31 * @param p2 The second point.
32 * @param p3 The third point.
33 * @return true if a circle can be found which passes through the points.
34 */
35 public boolean fit(double[] p1, double[] p2, double[] p3)
36 {
37 int i;
38 double r, m11, m12, m13, m14;
39
40 double[][] P = {p1, p2, p3};
41
42 for (i = 0; i < 3; i++) // find minor 11
43 {
44 a[i][0] = P[i][0];
45 a[i][1] = P[i][1];
46 a[i][2] = 1;
47 }
48 m11 = determinant( a, 3 );
49 if(m11==0)
50 {
51 _fit = null;
52 return false;
53 }
54
55 for (i = 0; i < 3; i++) // find minor 12
56 {
57 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1];
58 a[i][1] = P[i][1];
59 a[i][2] = 1;
60 }
61 m12 = determinant( a, 3 );
62
63 for (i = 0; i < 3; i++) // find minor 13
64 {
65 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1];
66 a[i][1] = P[i][0];
67 a[i][2] = 1;
68 }
69 m13 = determinant( a, 3 );
70
71 for (i = 0; i < 3; i++) // find minor 14
72 {
73 a[i][0] = P[i][0]*P[i][0] + P[i][1]*P[i][1];
74 a[i][1] = P[i][0];
75 a[i][2] = P[i][1];
76 }
77 m14 = determinant( a, 3 );
78
79 double x0 = 0.5 * m12 / m11; // center of circle
80 double y0 = -0.5 * m13 / m11;
81 r = sqrt( x0*x0 + y0*y0 + m14/m11 );
82
83 _fit = new CircleFit(x0, y0, r);
84
85 return true;
86 }
87
88 /**
89 * Get the result of the fit.
90 * @return The Circle which passes through the points.
91 */
92 public CircleFit getFit()
93 {
94 return _fit;
95 }
96
97 // Recursive definition of determinant using expansion by minors.
98 private double determinant(double[][] a, int n)
99 {
100 int i, j, j1, j2;
101 double d = 0.;
102
103 if (n == 2) // terminate recursion
104 {
105 d = a[0][0]*a[1][1] - a[1][0]*a[0][1];
106 }
107 else
108 {
109 d = 0;
110 for (j1 = 0; j1 < n; ++j1 ) // do each column
111 {
112 for (i = 1; i < n; ++i) // create minor
113 {
114 j2 = 0;
115 for (j = 0; j < n; ++j)
116 {
117 if (j == j1) continue;
118 m[i-1][j2] = a[i][j];
119 j2++;
120 }
121 }
122 // sum (+/-)cofactor * minor
123 d += pow(-1.0, j1)*a[0][j1]*determinant( m, n-1 );
124 }
125 }
126 return d;
127 }
128 }