Logo ROOT   6.12/07
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
91  TMatrixD S1(TMatrixD::kAtA, D1);
92  // combined part of the scatter matrix
94  // linear part of the scatter matrix
95  TMatrixD S3(TMatrixD::kAtA, D2);
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];
271  Double_t dt = TMath::TwoPi() / Double_t(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
static double B[]
TVectorD ConicToParametric(const TVectorD &conic)
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:292
float xmin
Definition: THbookFile.cxx:93
TMatrixDEigen.
Definition: TMatrixDEigen.h:26
TVectorD fit_ellipse(TGraph *g)
constexpr Double_t TwoPi()
Definition: TMath.h:44
float ymin
Definition: THbookFile.cxx:93
const TVectorD & GetEigenValuesIm() const
Definition: TMatrixDEigen.h:57
double T(double x)
Definition: ChebyshevPol.h:34
TVectorT.
Definition: TMatrixTBase.h:77
#define gROOT
Definition: TROOT.h:402
#define N
Int_t GetNrows() const
Definition: TVectorT.h:75
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
virtual void SetFillStyle(Style_t fstyle)
Set the fill area style.
Definition: TAttFill.h:39
virtual void SetNameTitle(const char *name, const char *title)
Set all the TNamed parameters (name and title).
Definition: TNamed.cxx:154
static double A[]
#define G(x, y, z)
double cos(double)
TMatrixT.
Definition: TMatrixDfwd.h:22
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
virtual void SetGrid(Int_t valuex=1, Int_t valuey=1)
Definition: TPad.h:327
void Clear(Option_t *option="")
Remove all primitives from the canvas.
Definition: TCanvas.cxx:707
double sin(double)
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
Definition: TMatrixT.cxx:1396
virtual void Draw(Option_t *option="")
Draw this ellipse with its current attributes.
Definition: TEllipse.cxx:167
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:533
#define F(x, y, z)
float ymax
Definition: THbookFile.cxx:93
static double C[]
auto * a
Definition: textangle.C:12
Int_t GetN() const
Definition: TGraph.h:122
float xmax
Definition: THbookFile.cxx:93
Element Norm2Sqr() const
Compute the square of the 2-norm SUM{ v[i]^2 }.
Definition: TVectorT.cxx:582
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
const Bool_t kFALSE
Definition: RtypesCore.h:88
Double_t * GetX() const
Definition: TGraph.h:129
const TMatrixD & GetEigenVectors() const
Definition: TMatrixDEigen.h:55
The Canvas class.
Definition: TCanvas.h:31
double Double_t
Definition: RtypesCore.h:55
void fitEllipseTGraphDLSF(TGraph *g=((TGraph *) 0))
Double_t y[n]
Definition: legend1.C:17
double atan(double)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
TGraph * TestGraphDLSF(Bool_t randomize=kFALSE)
Draw Ellipses.
Definition: TEllipse.h:24
Double_t * GetY() const
Definition: TGraph.h:130
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual void ComputeRange(Double_t &xmin, Double_t &ymin, Double_t &xmax, Double_t &ymax) const
Compute the x/y range of the points in this graph.
Definition: TGraph.cxx:645
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2248
#define I(x, y, z)
const Bool_t kTRUE
Definition: RtypesCore.h:87
const Int_t n
Definition: legend1.C:16
constexpr Double_t PiOver2()
Definition: TMath.h:48
void Modified(Bool_t flag=1)
Definition: TPad.h:414
static constexpr double g