Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
fitEllipseTGraphDLSF.cxx
Go to the documentation of this file.
1//
2// The "fitEllipseTGraphDLSF" macro uses the "Direct Least Squares Fitting"
3// algorithm for fitting an ellipse to a set of data points from a TGraph
4//
5// To try this macro, in a ROOT prompt, do:
6// .L fitEllipseTGraphDLSF.cxx // or ".L fitEllipseTGraphDLSF.cxx++"
7// fitEllipseTGraphDLSF(TestGraphDLSF());
8// for (int i=0; i<10; i++) { fitEllipseTGraphDLSF(); gSystem->Sleep(333); }
9//
10// Last update: Thu Jul 31 18:00:00 UTC 2014
11//
12// Changes:
13// 2014.07.31 - (initial version)
14//
15
16#include "TROOT.h"
17#include "TMath.h"
18#include "TRandom.h"
19#include "TGraph.h"
20#include "TVectorD.h"
21#include "TMatrixD.h"
22#include "TMatrixDEigen.h"
23#include "TCanvas.h"
24#include "TEllipse.h"
25
26#include <cmath>
27#include <iostream>
28
29//
30// "NUMERICALLY STABLE DIRECT LEAST SQUARES FITTING OF ELLIPSES"
31// Radim Halir, Jan Flusser
32// http://autotrace.sourceforge.net/WSCG98.pdf
33//
34// http://en.wikipedia.org/wiki/Ellipse
35//
36// An "algebraic distance" of a point "(x, y)" to a given conic:
37// F(x, y) = A * (x - X0)^2 + B * (x - X0) * (y - Y0) + C * (y - Y0)^2
38// + D * (x - X0) + E * (y - Y0) + F
39//
40// Ellipse-specific constraints:
41// F(x, y) = 0
42// B^2 - 4 * A * C < 0
43//
44// input parameter is a pointer to a "TGraph" with at least 6 points
45//
46// returns a "TVectorD" ("empty" in case any problems encountered):
47// ellipse[0] = "X0"
48// ellipse[1] = "Y0"
49// ellipse[2] = "A"
50// ellipse[3] = "B"
51// ellipse[4] = "C"
52// ellipse[5] = "D"
53// ellipse[6] = "E"
54// ellipse[7] = "F"
55//
57{
58 TVectorD ellipse;
59
60 if (!g) return ellipse; // just a precaution
61 if (g->GetN() < 6) return ellipse; // just a precaution
62
63 int i;
64 double tmp;
65
66 int N = g->GetN();
67 double xmin, xmax, ymin, ymax, X0, Y0;
68 g->ComputeRange(xmin, ymin, xmax, ymax);
69#if 1 /* 0 or 1 */
70 X0 = (xmax + xmin) / 2.0;
71 Y0 = (ymax + ymin) / 2.0;
72#else /* 0 or 1 */
73 X0 = Y0 = 0.0;
74#endif /* 0 or 1 */
75
76 TMatrixD D1(N, 3); // quadratic part of the design matrix
77 TMatrixD D2(N, 3); // linear part of the design matrix
78
79 for (i = 0; i < N; i++) {
80 double x = (g->GetX())[i] - X0;
81 double y = (g->GetY())[i] - Y0;
82 D1[i][0] = x * x;
83 D1[i][1] = x * y;
84 D1[i][2] = y * y;
85 D2[i][0] = x;
86 D2[i][1] = y;
87 D2[i][2] = 1.0;
88 }
89
90 // quadratic part of the scatter matrix
92 // combined part of the scatter matrix
94 // linear part of the scatter matrix
96 S3.Invert(&tmp); S3 *= -1.0;
97 if (tmp == 0.0) {
98 std::cout << "fit_ellipse : linear part of the scatter matrix is singular!" << std::endl;
99 return ellipse;
100 }
101 // for getting a2 from a1
103 // reduced scatter matrix
104 TMatrixD M(S2, TMatrixD::kMult, T); M += S1;
105 // premultiply by inv(C1)
106 for (i = 0; i < 3; i++) {
107 tmp = M[0][i] / 2.0;
108 M[0][i] = M[2][i] / 2.0;
109 M[2][i] = tmp;
110 M[1][i] *= -1.0;
111 }
112 // solve eigensystem
113 TMatrixDEigen eig(M); // note: eigenvectors are not normalized
114 const TMatrixD &evec = eig.GetEigenVectors();
115 // const TVectorD &eval = eig.GetEigenValuesRe();
116 if ((eig.GetEigenValuesIm()).Norm2Sqr() != 0.0) {
117 std::cout << "fit_ellipse : eigenvalues have nonzero imaginary parts!" << std::endl;
118 return ellipse;
119 }
120 // evaluate a’Ca (in order to find the eigenvector for min. pos. eigenvalue)
121 for (i = 0; i < 3; i++) {
122 tmp = 4.0 * evec[0][i] * evec[2][i] - evec[1][i] * evec[1][i];
123 if (tmp > 0.0) break;
124 }
125 if (i > 2) {
126 std::cout << "fit_ellipse : no min. pos. eigenvalue found!" << std::endl;
127 // i = 2;
128 return ellipse;
129 }
130 // eigenvector for min. pos. eigenvalue
131 TVectorD a1(TMatrixDColumn_const(evec, i));
132 tmp = a1.Norm2Sqr();
133 if (tmp > 0.0) {
134 a1 *= 1.0 / std::sqrt(tmp); // normalize this eigenvector
135 } else {
136 std::cout << "fit_ellipse : eigenvector for min. pos. eigenvalue is NULL!" << std::endl;
137 return ellipse;
138 }
139 TVectorD a2(T * a1);
140
141 // ellipse coefficients
142 ellipse.ResizeTo(8);
143 ellipse[0] = X0; // "X0"
144 ellipse[1] = Y0; // "Y0"
145 ellipse[2] = a1[0]; // "A"
146 ellipse[3] = a1[1]; // "B"
147 ellipse[4] = a1[2]; // "C"
148 ellipse[5] = a2[0]; // "D"
149 ellipse[6] = a2[1]; // "E"
150 ellipse[7] = a2[2]; // "F"
151
152 return ellipse;
153}
154
155//
156// http://mathworld.wolfram.com/Ellipse.html
157// http://mathworld.wolfram.com/QuadraticCurve.html
158// http://mathworld.wolfram.com/ConicSection.html
159//
160// "Using the Ellipse to Fit and Enclose Data Points"
161// Charles F. Van Loan
162// http://www.cs.cornell.edu/cv/OtherPdf/Ellipse.pdf
163//
164// input parameter is a reference to a "TVectorD" which describes
165// an ellipse according to the equation:
166// 0 = A * (x - X0)^2 + B * (x - X0) * (y - Y0) + C * (y - Y0)^2
167// + D * (x - X0) + E * (y - Y0) + F
168// conic[0] = "X0"
169// conic[1] = "Y0"
170// conic[2] = "A"
171// conic[3] = "B"
172// conic[4] = "C"
173// conic[5] = "D"
174// conic[6] = "E"
175// conic[7] = "F"
176//
177// returns a "TVectorD" ("empty" in case any problems encountered):
178// ellipse[0] = ellipse's "x" center ("x0")
179// ellipse[1] = ellipse's "y" center ("y0")
180// ellipse[2] = ellipse's "semimajor" axis along "x" ("a" > 0)
181// ellipse[3] = ellipse's "semiminor" axis along "y" ("b" > 0)
182// ellipse[4] = ellipse's axes rotation angle ("theta" = -45 ... 135 degrees)
183//
185{
186 TVectorD ellipse;
187
188 if (conic.GetNrows() != 8) {
189 std::cout << "ConicToParametric : improper input vector length!" << std::endl;
190 return ellipse;
191 }
192
193 double a, b, theta;
194 double x0 = conic[0]; // = X0
195 double y0 = conic[1]; // = Y0
196
197 // http://mathworld.wolfram.com/Ellipse.html
198 double A = conic[2];
199 double B = conic[3] / 2.0;
200 double C = conic[4];
201 double D = conic[5] / 2.0;
202 double F = conic[6] / 2.0;
203 double G = conic[7];
204
205 double J = B * B - A * C;
206 double Delta = A * F * F + C * D * D + J * G - 2.0 * B * D * F;
207 double I = - (A + C);
208
209 // http://mathworld.wolfram.com/QuadraticCurve.html
210 if (!( (Delta != 0.0) && (J < 0.0) && (I != 0.0) && (Delta / I < 0.0) )) {
211 std::cout << "ConicToParametric : ellipse (real) specific constraints not met!" << std::endl;
212 return ellipse;
213 }
214
215 x0 += (C * D - B * F) / J;
216 y0 += (A * F - B * D) / J;
217
218 double tmp = std::sqrt((A - C) * (A - C) + 4.0 * B * B);
219 a = std::sqrt(2.0 * Delta / J / (I + tmp));
220 b = std::sqrt(2.0 * Delta / J / (I - tmp));
221
222 theta = 0.0;
223 if (B != 0.0) {
224 tmp = (A - C) / 2.0 / B;
225 theta = -45.0 * (std::atan(tmp) / TMath::PiOver2());
226 if (tmp < 0.0) { theta -= 45.0; } else { theta += 45.0; }
227 if (A > C) theta += 90.0;
228 } else if (A > C) theta = 90.0;
229
230 // try to keep "a" > "b"
231 if (a < b) { tmp = a; a = b; b = tmp; theta -= 90.0; }
232 // try to keep "theta" = -45 ... 135 degrees
233 if (theta < -45.0) theta += 180.0;
234 if (theta > 135.0) theta -= 180.0;
235
236 // ellipse coefficients
237 ellipse.ResizeTo(5);
238 ellipse[0] = x0; // ellipse's "x" center
239 ellipse[1] = y0; // ellipse's "y" center
240 ellipse[2] = a; // ellipse's "semimajor" axis along "x"
241 ellipse[3] = b; // ellipse's "semiminor" axis along "y"
242 ellipse[4] = theta; // ellipse's axes rotation angle (in degrees)
243
244 return ellipse;
245}
246
247//
248// creates a test TGraph with an ellipse
249//
250TGraph *TestGraphDLSF(bool randomize = false)
251{
252 // define the test ellipse
253 double x0 = 4; // ellipse's "x" center
254 double y0 = 3; // ellipse's "y" center
255 double a = 2; // ellipse's "semimajor" axis along "x" (> 0)
256 double b = 1; // ellipse's "semiminor" axis along "y" (> 0)
257 double theta = 100; // ellipse's axes rotation angle (-45 ... 135 degrees)
258
259 // gRandom->SetSeed(0);
260 if (randomize) {
261 x0 = 10.0 - 20.0 * gRandom->Rndm();
262 y0 = 10.0 - 20.0 * gRandom->Rndm();
263 a = 0.5 + 4.5 * gRandom->Rndm();
264 b = 0.5 + 4.5 * gRandom->Rndm();
265 theta = 180.0 - 360.0 * gRandom->Rndm();
266 }
267
268 const int n = 100; // number of points
269 double x[n], y[n];
270 double dt = TMath::TwoPi() / double(n);
271 double tmp;
272 theta *= TMath::PiOver2() / 90.0; // degrees -> radians
273 for (int i = 0; i < n; i++) {
274 x[i] = a * (std::cos(dt * double(i)) + 0.1 * gRandom->Rndm() - 0.05);
275 y[i] = b * (std::sin(dt * double(i)) + 0.1 * gRandom->Rndm() - 0.05);
276 // rotate the axes
277 tmp = x[i];
278 x[i] = x[i] * std::cos(theta) - y[i] * std::sin(theta);
279 y[i] = y[i] * std::cos(theta) + tmp * std::sin(theta);
280 // shift the center
281 x[i] += x0;
282 y[i] += y0;
283 }
284
285 // create the test TGraph
286 TGraph *g = static_cast<TGraph *>(gROOT->FindObject("g"));
287 if (g) delete g;
288 g = new TGraph(n, x, y);
289 g->SetNameTitle("g", "test ellipse");
290
291 return g;
292}
293
294//
295// "ROOT Script" entry point (the same name as the "filename's base")
296//
298{
299 if (!g) g = TestGraphDLSF(true); // create a "random" ellipse
300
301 // fit the TGraph
302 TVectorD conic = fit_ellipse(g);
303 TVectorD ellipse = ConicToParametric(conic);
304
305 if ( ellipse.GetNrows() == 5 ) {
306 std::cout << std::endl;
307 std::cout << "x0 = " << ellipse[0] << std::endl;
308 std::cout << "y0 = " << ellipse[1] << std::endl;
309 std::cout << "a = " << ellipse[2] << std::endl;
310 std::cout << "b = " << ellipse[3] << std::endl;
311 std::cout << "theta = " << ellipse[4] << std::endl;
312 std::cout << std::endl;
313 }
314
315 // draw everything
316 auto c1 = new TCanvas("c1","c1", 1000, 800);
317 c1->SetGrid(1, 1);
318 g->Draw("A*");
319 if ( ellipse.GetNrows() == 5 ) {
320 TEllipse *e = new TEllipse(ellipse[0], ellipse[1], // "x0", "y0"
321 ellipse[2], ellipse[3], // "a", "b"
322 0, 360,
323 ellipse[4]); // "theta" (in degrees)
324 e->SetFillStyle(0); // hollow
325 e->Draw();
326 }
327 c1->Modified();
328 c1->Update(); // make sure it's really drawn
329}
330
331// end of file fitEllipseTGraphDLSF.cxx by Silesius Anonymus
#define b(i)
Definition RSha256.hxx:100
#define S1(x)
Definition RSha256.hxx:89
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
#define N
float xmin
float ymin
float xmax
float ymax
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
#define gROOT
Definition TROOT.h:406
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
The Canvas class.
Definition TCanvas.h:23
Draw Ellipses.
Definition TEllipse.h:23
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
TMatrixDEigen.
const TMatrixD & GetEigenVectors() const
const TVectorD & GetEigenValuesIm() const
TMatrixT< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant.
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:559
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition TVectorT.cxx:294
Int_t GetNrows() const
Definition TVectorT.h:73
Element Norm2Sqr() const
Compute the square of the 2-norm SUM{ v[i]^2 }.
Definition TVectorT.cxx:584
TVectorD fit_ellipse(TGraph *g)
TGraph * TestGraphDLSF(bool randomize=false)
TVectorD ConicToParametric(const TVectorD &conic)
void fitEllipseTGraphDLSF(TGraph *g=nullptr)
Double_t y[n]
Definition legend1.C:17
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
#define F(x, y, z)
#define I(x, y, z)
#define G(x, y, z)
constexpr Double_t PiOver2()
Definition TMath.h:51
constexpr Double_t TwoPi()
Definition TMath.h:44