Logo ROOT  
Reference Guide
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_t 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 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_t i;
64 Double_t tmp;
65
66 Int_t N = g->GetN();
67 Double_t 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_t x = (g->GetX())[i] - X0;
81 Double_t 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_t a, b, theta;
194 Double_t x0 = conic[0]; // = X0
195 Double_t y0 = conic[1]; // = Y0
196
197 // http://mathworld.wolfram.com/Ellipse.html
198 Double_t A = conic[2];
199 Double_t B = conic[3] / 2.0;
200 Double_t C = conic[4];
201 Double_t D = conic[5] / 2.0;
202 Double_t F = conic[6] / 2.0;
203 Double_t G = conic[7];
204
205 Double_t J = B * B - A * C;
206 Double_t Delta = A * F * F + C * D * D + J * G - 2.0 * B * D * F;
207 Double_t 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_t 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//
251 Int_t i;
252
253 // define the test ellipse
254 Double_t x0 = 4; // ellipse's "x" center
255 Double_t y0 = 3; // ellipse's "y" center
256 Double_t a = 2; // ellipse's "semimajor" axis along "x" (> 0)
257 Double_t b = 1; // ellipse's "semiminor" axis along "y" (> 0)
258 Double_t theta = 100; // ellipse's axes rotation angle (-45 ... 135 degrees)
259
260 // gRandom->SetSeed(0);
261 if (randomize) {
262 x0 = 10.0 - 20.0 * gRandom->Rndm();
263 y0 = 10.0 - 20.0 * gRandom->Rndm();
264 a = 0.5 + 4.5 * gRandom->Rndm();
265 b = 0.5 + 4.5 * gRandom->Rndm();
266 theta = 180.0 - 360.0 * gRandom->Rndm();
267 }
268
269 const Int_t n = 100; // number of points
270 Double_t x[n], y[n];
272 Double_t tmp;
273 theta *= TMath::PiOver2() / 90.0; // degrees -> radians
274 for (i = 0; i < n; i++) {
275 x[i] = a * (std::cos(dt * Double_t(i)) + 0.1 * gRandom->Rndm() - 0.05);
276 y[i] = b * (std::sin(dt * Double_t(i)) + 0.1 * gRandom->Rndm() - 0.05);
277 // rotate the axes
278 tmp = x[i];
279 x[i] = x[i] * std::cos(theta) - y[i] * std::sin(theta);
280 y[i] = y[i] * std::cos(theta) + tmp * std::sin(theta);
281 // shift the center
282 x[i] += x0;
283 y[i] += y0;
284 }
285
286 // create the test TGraph
287 TGraph *g = ((TGraph *)(gROOT->FindObject("g")));
288 if (g) delete g;
289 g = new TGraph(n, x, y);
290 g->SetNameTitle("g", "test ellipse");
291
292 return g;
293}
294
295//
296// "ROOT Script" entry point (the same name as the "filename's base")
297//
299{
300 if (!g) g = TestGraphDLSF(kTRUE); // create a "random" ellipse
301
302 // fit the TGraph
303 TVectorD conic = fit_ellipse(g);
304 TVectorD ellipse = ConicToParametric(conic);
305
306#if 1 /* 0 or 1 */
307 if ( ellipse.GetNrows() == 5 ) {
308 std::cout << std::endl;
309 std::cout << "x0 = " << ellipse[0] << std::endl;
310 std::cout << "y0 = " << ellipse[1] << std::endl;
311 std::cout << "a = " << ellipse[2] << std::endl;
312 std::cout << "b = " << ellipse[3] << std::endl;
313 std::cout << "theta = " << ellipse[4] << std::endl;
314 std::cout << std::endl;
315 }
316#endif /* 0 or 1 */
317
318#if 1 /* 0 or 1 */
319 // draw everything
320 TCanvas *c = ((TCanvas *)(gROOT->GetListOfCanvases()->FindObject("c")));
321 if (c) { c->Clear(); } else { c = new TCanvas("c", "c"); }
322 c->SetGrid(1, 1);
323 g->Draw("A*");
324 if ( ellipse.GetNrows() == 5 ) {
325 TEllipse *e = new TEllipse(ellipse[0], ellipse[1], // "x0", "y0"
326 ellipse[2], ellipse[3], // "a", "b"
327 0, 360,
328 ellipse[4]); // "theta" (in degrees)
329 e->SetFillStyle(0); // hollow
330 e->Draw();
331 }
332 c->Modified(); c->Update(); // make sure it's really drawn
333#endif /* 0 or 1 */
334
335 return;
336}
337
338// end of file fitEllipseTGraphDLSF.cxx by Silesius Anonymus
#define b(i)
Definition: RSha256.hxx:100
#define S1(x)
Definition: RSha256.hxx:89
#define c(i)
Definition: RSha256.hxx:101
#define g(i)
Definition: RSha256.hxx:105
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define N
float xmin
Definition: THbookFile.cxx:93
float ymin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
float ymax
Definition: THbookFile.cxx:93
double cos(double)
double atan(double)
double sqrt(double)
double sin(double)
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:27
Draw Ellipses.
Definition: TEllipse.h:24
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
TMatrixDEigen.
Definition: TMatrixDEigen.h:27
const TMatrixD & GetEigenVectors() const
Definition: TMatrixDEigen.h:55
const TVectorD & GetEigenValuesIm() const
Definition: TMatrixDEigen.h:57
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
Definition: TMatrixT.cxx:1399
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:541
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:295
Int_t GetNrows() const
Definition: TVectorT.h:75
Element Norm2Sqr() const
Compute the square of the 2-norm SUM{ v[i]^2 }.
Definition: TVectorT.cxx:585
TGraph * TestGraphDLSF(Bool_t randomize=kFALSE)
TVectorD fit_ellipse(TGraph *g)
TVectorD ConicToParametric(const TVectorD &conic)
void fitEllipseTGraphDLSF(TGraph *g=((TGraph *) 0))
Double_t y[n]
Definition: legend1.C:17
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)
static double B[]
static double A[]
static double C[]
double T(double x)
Definition: ChebyshevPol.h:34
constexpr Double_t PiOver2()
Definition: TMath.h:52
constexpr Double_t TwoPi()
Definition: TMath.h:45
auto * a
Definition: textangle.C:12