Discussion :

[souviron34]

Après 7 ans (!!!!) d'intervalle depuis ma promesse de traduire le code Fortran mis dans cette rubrique ici-même ([FORTRAN] Ellipse fitting algorithm), je finis par la tenir

L'"ellipse-fitting", ou en .. franglais... le "fit d'ellipse", est l'action de trouver les paramètres d'une ellipse définie par un ensemble de points 2D isolés et bruités, que l'on pense figurer sur une ellipse.

En exemple, voir l'image donnée sur le post http://www.developpez.net/forums/m2213840-1/... . Ou les pièces jointes ci-dessous...

En fait, on a au départ une image (ou un graphique) brut. Certains pixels (positions) sont plus ou moins répartis suivant ce qu'à l'oeil on pense être une ellipse. Cet algorithme trouve les paramètres (1/2 petit et grand axe, ellipticité, centre, angle par rapport à l'horizontale) en prenant l'ensemble des coordonnées des points concernés et en appliquant une méthode des moindres carrés sur la développée de l'équation de l'ellipse.

En fait, le principe est extrêmement simple :

• Une ellipse alignée en X,Y centrée à l'origine a pour équation :
  
  X²/a² + Y²/b² = 1
  
  Où a et b sont respectivement appelés demi-grand axe et demi-petit axe.
  
  
• Si maintenant le centre n'est pas à l'origine, l'équation devient :
  
  (X-Xc)²/a² + (Y-Yc)²/b² = 1
  
  
• Si maintenant l'ellipse n'est pas alignée avec X,Y mais est inclinée d'un angle alpha, l'équation devient :
  
  Xa²/a² + Ya²/b² = 1
  
  Où
  Xa = (X-Xc) cos alpha + (Y-Yc) sin alpha
  Ya = -(X-Xc) sin alpha + (Y-Yc) cos alpha

La solution se fait donc simplement en développant l'équation (on a donc des termes en x², y², xy, x, y si on ne connaît pas le centre exactement, ou en x², y², xy si on connaît le centre).

On construit une matrice correspondant pour chaque point de mesure aux x², y², xy, x, y (respectivement x², y², xy), et en ajoutant une colonne qui contient simplement la valeur "1".

Un algorithme de moindre-carrés sur l'ensemble de cette matrice donne les coefficients de chacune des colonnes, et donc un simple rétro-calcul donne et les 1/2 grand et petit axe a et b, l'ellipticité 1-b/a, l'angle d'inclinaison avec l'horizontale (de l'ellipse VERS X et non l'inverse) , et le sigma total de l'ensemble des mesures par rapport à la courbe obtenue.

Des exemples sont fournis en pièces jointes. La complexité est linéaire en temps et en espace (on va fabriquer une matrice N*6 doubles au maximum).


Le recours aux moindre-carrés sous cette forme permet de trouver des paramètres cohérents même avec très peu de points (voir exemple ci-dessous)(en astrophysique, domaine d'origine de ce travail, cela permet de descendre au moins 2 magnitudes plus bas que des méthodes plus sophistiquées comme les coefficients de Fourier)



Par rapport à la version Fortran, tant qu'à faire j'ai profité du passage en C pour faire un moindre-carrés généralisé à n'importe quel nombre de paramètres en entrée.. (plutôt que tenter de trouver le bug (sans doute des inversions colonnes/lignes entre Fortran et C) dans la traduction de la routine d'élimination gaussienne originale, j'ai trouvé la correspondante déjà traduite dans "Numerical Recipes" ).. Voici donc les routines de maths initiales (moindres carrés généralisés à N paramètres, et résolution d'un système linéaire par élimination de gaussienne) :

Code C :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
 
#ifndef PI
#define PI 3.14159265359
#endif
#ifndef Boolean
#define Boolean char
#endif
#ifndef SUCCESS
#define SUCCESS 1
#endif
#ifndef ERROR
#define ERROR 0
#endif
#ifndef True
#define True 1
#define False 0
#endif
 
 
#define ACCURACY 1.0e-12
 
 
 
/*
***********************************************************
       SOLVE LINEAR SYSTEM VIA GAUSSIAN ELIMINATION
 
 Source code taken from "Numerical Recipes" in Java at
    http://introcs.cs.princeton.edu/java/95linear/GaussianElimination.java.html
 
 and ported to C by Jean Souviron 2016/01/06
 
 
***********************************************************
*/
int LinearSystemByGaussian (double *A, double *X, int Order )
{
  int    N=Order, Length=Order+1, max, i, p, j ;
  double t, sum, alpha ;
 
 
  for (p = 0; p < N; p++) {
 
    /* find pivot row and swap */
    max = p;
    for (i = p + 1; i < N; i++) {
      if (fabs(A[i*Length+p]) > fabs(A[max*Length+p])) {
	max = i;
      }
    }
 
    /* Exchange rows */
    for ( i = 0 ; i < Length ; i++ )
      {
	t = A[p*Length+i] ;
	A[p*Length+i] = A[max*Length+i] ;
	A[max*Length+i] = t ;
      }
 
 
    /* singular or nearly singular */
    if (fabs(A[p*Length+p]) <= ACCURACY) {
      for ( j = 0 ; j < Order ; j++ )
	X[j] = 0.0 ;
      return 0 ;
    }
 
    /* pivot within A and b */
    for (i = p + 1; i < N; i++) {
      alpha = A[i*Length+p] / A[p*Length+p];
      A[i*Length+N] -= alpha * A[p*Length+N];
      for (j = p; j < N; j++) {
	A[i*Length+j] -= alpha * A[p*Length+j];
      }
    }
  }
 
  /* back substitution */
  for (i = N - 1; i >= 0; i--) {
    sum = 0.0;
    for (j = i + 1; j < N; j++) {
      sum += A[i*Length+j] * X[j];
    }
    X[i] = (A[i*Length+N] - sum) / A[i*Length+i];
  }
 
  return 1;
}
 
 
/*
***********************************************************************
      GENERALIZED LEAST-SQUARE FOR ANY NUMBER OF PARAMETERS
 
                      (Jean SOUVIRON - DAO Victoria - 1985)
 
         (Method of R. SEDGEWICK in "Algorithms"
                                    addison-wesley publishing co. 1984)
 
      Modified 2015/12/12 Jean Souviron
         Porting to C
	 Upgrading to take any number of parameters
 
-----------------------------------------------------------------------
 
      INPUT :
                 NDims = number of different parameters+1 (for the data to be fitted to)
                  NPts = total number of points
           InputMatrix = array of parameters for each point
                         (last vector has to be the data to be fitted)
      OUTPUT :
          OutputCoeffs = array of coefficients (dimension NDims-1 )
                 Sigma
 
      RETURN VALUE :
           1 if everything is OK
           0 if an error occurred (memory or singular array)
 
 
  Example : if we want to fit an ellipse with a group of points according to
            equation x2/a2 + y2/b2 = 1
 
            the input matrix will be N pts long, and 3 columns wide, with for
            each point the x2, the y2, and 1.
            the output coeffs will be 2 items, 1/a2 and 1/b2
            the sigma will be computed by applying these coeffs to the first 2 columns
	    for each point, and comparing with the result (1).
            
 
***********************************************************************
*/
int LeastSquares ( double *InputMatrix, int NDims, int NPts, 
		   double *OutputCoeffs, double *Sigma )
{
  double *Coeffs=NULL ;
  double  T, Z1, Z2 ;
  int     i,j, k, io = 0, ko, NParams ;
 
 
  NParams = NDims-1 ;
 
/*
  Allocations
*/
  if ( (Coeffs = (double *)calloc(NParams*NDims, sizeof(double))) == NULL )
    {
      return 0 ;
    }
 
/*
  Initializations
*/
  io = 0 ;
  for ( i = 0 ; i < NParams ; i++ )
    {
      OutputCoeffs[i] = 0.0 ;
 
      for ( j = 0 ; j < NDims ; j++ )
	{
	  T = 0.0 ;
 
	  for ( k = 0 ; k < NPts ; k++ )
	    {
	      ko = k * NDims ;
 
	      Z1 = InputMatrix[ko+i] ;
	      Z2 = InputMatrix[ko+j] ;
 
	      T = T + (Z1*Z2) ;
	    }
 
	  Coeffs[io+j]= T ;
	}
      io = io + NDims ;
    }
 
/*
-----------------
   Solving by gaussian elimination
-----------------
*/
  i = LinearSystemByGaussian ( Coeffs, OutputCoeffs, NParams ) ;
 
  free ( Coeffs );
  if ( i != 1 )
    return 0 ;
 
/*
-----------------
  Computes sigma
-----------------
*/
  *Sigma = 0.0 ;
  for ( i = 0 ; i < NPts ; i++ )
    { 
      Z1 = InputMatrix[i*NDims+NParams]; /* The value to be compared with */
 
      for ( j = 0 ; j < NParams ; j++ )
	{
	  Z1= Z1 - (OutputCoeffs[j]*InputMatrix[i*NDims+j]) ;
	}
 
      if( Z1 < 0.0 ) 
	Z1 = -Z1;
 
      *Sigma = *Sigma + Z1 ;
    }
 
  *Sigma = (*Sigma/(double)NPts)*sqrt(PI/2.0) ;
 
  return 1 ;
}

Et ensuite la vraie routine de fit d'ellipse :

Code C :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
 
/*
*******************************************************************
               SUBROUTINE OF ELLIPSES FITTING
 
        It takes the general equation that an ellipse has to 
        satisfy and fit by least-square the coefficients.
 
        If the axis were aligned with x and y , the equation
        would be :   
	           (x-xc)**2/a2 + (y-yc)**2/b2 = 1
 
        With a certain angle alpha the coordinates transforms
        into :   
                 x' = cos (aplha) * (x-xc) - sin (alpha) * (y-yc)
                 y' = sin (alpha) * (x-xc) + cos (alpha) * (y-yc)
 
        Replacing in the former equation gives an equation whose
        form is :   
                 a (x-xc)**2 + b (y-yc)**2 + c x*y = 1
 
        A very simple calculus gives the answer for the value
        of the different parameters when the coefficients a,b,c
        are determined by least-square.
 
        This routine allows the user to fit ellipses where :
             - the center is either known or unknown
             - the position angle is either known or unknown
             - the ellipticity is either knonw or unknown
 
 
    INPUT :
        XX              : table of X coordinates
        YY              : table of Y coordinates
        NPts            : number of points that have to be fitted
        XCenter,YCenter : coordinates of center (input/output)
        FixedCenter     : flag to compute the center (0) or take it from input (1)
                          (in case of no determination the values of the coordinates
                           have to be passed through XCenter and YCenter)
        FixedAngle      : flag to compute the axis angle (0) or take it from input (1)
	                  (in case of no determination, the value of the
                           angle has to be passed through ANG)
        FixedEll        : flag to compute the ellipticity (0) or take it from the input (1)
                          (in case of no determination, the value of the
                           ellipticity has to be passed through ELL)
 
    OUTPUT :
        RMajor, RMinor, RAverage,
        Ellipticity, Angle : parameters of the ellipse
 
        XCenter,YCenter    : eventually new coordinates center
 
        Sigma              : sigma on the whole set
 
 
    RETURN VALUE :
 
        1 on success
        0 if memory error
       -1 to -10 if another error occured
 
 
        ( Jean SOUVIRON - DAO Victoria 1985/86)
 
   Modified 2015/12/13 by Jean Souviron
     Porting to C
 
*******************************************************************
*/
int FitEllipse ( double *XX, double *YY, int NPts,
		 double *XCenter, double *YCenter,
		 Boolean FixedCenter, Boolean FixedAngle, Boolean FixedEll,
		 double *RMajor, double *RMinor, double *RAverage, 
		 double *Ellipticity, double *Angle,
		 double *Sigma )
{
  double  F1, F2, F3, G1, CorrX, CorrY, Correc, CorrecFactor, Denum ;
  double  Radian2Degree, A, B, Sinus, Alpha, Z, BoverA ;
  double *Table=NULL, Facto[5] ;
  int     i, NDims, Base, Status=SUCCESS ;
  Boolean IsReverse = False ; 
 
/*
--- Adjusting of parameters and memory allocation depending on the case
*/
  if ( FixedCenter ) 
      NDims = 4 ;
  else
      NDims = 6 ;
 
  if ( (Table =  (double *) calloc ( NDims*NPts, sizeof(double) )) == NULL )
    {
      return 0 ;
    }
 
/*
--- Determination of the ellipse parameters through a least-square
*/ 
  switch ( NDims )
    {
/*
-----------------------------
  CENTER HAS TO BE FOUND
-----------------------------
*/
      case 6 :
 
	  CorrX = 0.0 ;
	  CorrY = 0.0 ;
 
	  /*
	    Builds the matrix
	  */
	  for ( i = 0 ; i < NPts ; i++ )
	    {
	      Base = i * NDims ;
 
	      A = XX[i] - (*XCenter) ;
	      B = YY[i] - (*YCenter) ;
 
	      Table[Base] = A*A ;
	      Table[Base+1] = B*B ;
	      Table[Base+2] = A*B ;
	      Table[Base+3] = A ;
	      Table[Base+4] = B ;
	      Table[Base+5] = 1.0 ;
	    }
 
	  /*
	    Computes least-squares
	  */
	  if ( LeastSquares (Table, NDims, NPts, Facto, Sigma) != 1 )
	    {
	      Status = -1 ;
	      break ;
	    }
 
	  /*
	    Computes the denominator in the formula
	  */
	  Denum = (4.0 * Facto[0] * Facto[1]) - (Facto[2] * Facto[2]) ;
	  if ( Denum == 0.0)
	    {
	      Status = -2 ;
	      break ;
	    }
 
	  /* 
	     Computes and corrects center coordinates
	  */
	  CorrX = ( (-2.0*Facto[1]*Facto[3]) + (Facto[2]*Facto[4]) ) / Denum ;
	  *XCenter = *XCenter + CorrX ;
 
	  CorrY = ((Facto[2]*Facto[3]) - (2.0*Facto[0]*Facto[4])) / Denum ;
	  *YCenter = *YCenter + CorrY ;
 
	  /* 
	     Corrects the factors to take into account 
	     the "approximated" center
	  */
	  Correc = (Facto[0]*CorrX*CorrX) + (Facto[1]*CorrY*CorrY) + 
	           (Facto[2]*CorrX*CorrY) ;
 
	  if ( (1.0+Correc) == 0.0 )
	    {
	      Status = -3 ;
	      break ;
	    }
 
	  CorrecFactor = 1.0 - (Correc/(1.0+Correc)) ;
 
	  Facto[0] = Facto[0] * CorrecFactor ;
	  Facto[1] = Facto[1] * CorrecFactor ;
	  Facto[2] = Facto[2] * CorrecFactor ;
	  break ;
 
/*
-----------------------------
  CENTER ALREADY DETERMINED
-----------------------------
*/
     case 4 :
     default :
 
       /*
	 Builds the matrix
       */
       for ( i = 0 ; i < NPts ; i++ )
	 {
	   Base = i * NDims ;
 
	   A = XX[i] - (*XCenter) ;
	   B = YY[i] - (*YCenter) ;
 
	   Table[Base] = A*A ;
	   Table[Base+1] = B*B ;
	   Table[Base+2] = A*B ;
	   Table[Base+3] = 1.0 ;
	 }
 
       /*
	 Computes the least-squares
       */
       Status = LeastSquares (Table, NDims, NPts, Facto, Sigma);
       break ;
    }
 
  free ( Table );
  if ( Status != 1 )
    return Status ;
 
/*
-----------------------------------------
     DETERMINATION OF THE PARAMETERS
-----------------------------------------
*/
  F1 = Facto[0] ;
  F2 = Facto[1] ;
  F3 = Facto[2] ;
 
  Radian2Degree = 180.0 / PI ;
 
 
/*
***************
   CASE OF KNOWN ANGLE
**************
*/
  if ( FixedAngle )
    {
      Alpha = -(*Angle-90.0)/Radian2Degree ;
      Sinus = 2.0 * sin(Alpha) * cos(Alpha) ;
      if ( Sinus == 0.0)
	Sinus = 0.000001 ;
 
      G1 = F1 + F2 ;
      Z = (F3 / Sinus) + G1 ;
 
      if ( Z > 0.0 )
	{
	  *RMajor = sqrt(2.0/Z) ;
 
	  Z = G1 - (F3/Sinus) ;
	  if (Z > 0.0)
	    {
	      *RMinor = sqrt(2.0/Z) ;
	      if (*RMajor < *RMinor )
		{
		  Z = *RMajor ;
		  *RMajor = *RMinor ;
		  *RMinor = Z ;
		}
 
	      *Ellipticity = 1.0 - (*RMinor/ (*RMajor)) ;
	      *RAverage = sqrt(*RMinor*(*RMajor)) ;
	    }
	  else
	      Status = -4 ;
	}
      else
	  Status = -5 ;
    }
  else
/*
***************
   CASE OF KNOWN ELLIPTICITY
**************
*/
  if ( FixedEll )
    {
      BoverA = 1.0 - *Ellipticity ;
      Z = F1 + F2 ;
      if ( Z > 0.0 )
	{
	  *RMinor = sqrt( (1.0+(BoverA*BoverA)) / Z ) ;
	  *RMajor = *RMinor / BoverA ;
	  *RAverage = sqrt(*RMajor * (*RMinor)) ;
 
	  Z = (F1-F2) + (((BoverA*BoverA)-1.0)/(*RMinor *(*RMinor))) ;
	  if ( Z != 0.0 )
	      *Angle = (atan(- F3/Z)*Radian2Degree) + 90.0 ;
	  else
	      Status = -7 ;
	}
      else
	  Status = -6 ;
    }
/*
***************
   NORMAL CASE
**************
*/
  else
    {
      /*
	-------------- component on X-axis
      */
      Z = (F3*F3) + ((F1-F2)*(F1-F2)) ;
      if( Z > 0.0 )
	{
	  G1 = 2.0 / (F1+F2+sqrt(Z)) ;
	  if ( G1 > 0.0 )
	    {
	      *RMajor = sqrt(G1) ;
 
	      A = *RMajor ;
	      /*
		--------------- component on Y-axis
	      */
	      Z = F1 + F2 - (1.0/(A * A)) ;
	      if ( Z > 0.0 )
		{
		  *RMinor = 1.0 / sqrt(Z) ;
 
		  B = *RMinor ;
 
		  if ( *RMajor < *RMinor )
		    {
		      IsReverse = 1 ;
		      *RMajor = B ;
		      *RMinor = A ;
		    }
		}
	      else
		  Status = -10 ;
	    }
	  else
	      Status = -9 ;
	}
      else
	  Status = -8 ;
 
      /*
	If everything is OK
      */
      if ( Status > 0 )
	{
	  A = *RMajor ;
	  B = *RMinor ;
 
	  /*--- Ellipticity and mean axis ---*/
	  *Ellipticity = 1.0 - (*RMinor / (*RMajor)) ;
	  *RAverage = sqrt(*RMinor * (*RMajor)) ;
 
	  /*--- Angle ---*/
	  /*
	     1) x/HORIZONTAL
	        Sign - because it is the reverse that what we want 
		(from normal to ellipse)
	  */
	  Z = - F3 / (F1 - F2 + (1.0/(A*A)) - (1.0/(B*B))) ;
	  *Angle = atan(Z) * Radian2Degree ;
 
	  /*
	     2) Checks if X is major axis
	  */
	  if ( IsReverse )
	    {
	      if ( *Angle >= 0.0 )
		*Angle = *Angle - 90.0 ;
	      else
		*Angle = *Angle + 90.0 ;
	    }
 
	  /*
	     3) Angle major axis/Vertical, unclockwise
	  */
	  *Angle = *Angle + 90.0 ;
	}
    }
 
  return Status ;
}

Il est à noter que l'on doit donner une estimation du centre initalement (cela peut par exemple être le barycentre des points). Elle sera corrigée si on demande à déterminer le centre.

On peut par drapeau à l'appel préciser si on veut travailler avec un centre fixe ou le déterminer, avec une ellipticité fixe ou la déterminer, avec un angle fixe ou le déterminer.

L'angle est en degrés décimaux.

Le tracé de l'ellipse correspondant aux paramètres trouvés est

Code C :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 
	   newangle = -Angle * PI / 180.0 ;   /* le "-" est necessaire car là on veut l'angle de X VERS l'ellipse */
	   newcos = cos (newangle) ;
	   newsin = sin (newangle) ;
 
           delta = PI/180.0 ;
	   for ( i = 0 ; i < 360 ; i++ )
	     {
	       cosinus = cos ((double)i*delta) ;
	       sinus = sin ( (double)i*delta) ;
 
	       x1 = XCen + (RMaj * cosinus * newcos) - (RMin * sinus * newsin) ;
	       y1 = YCen + (RMaj * cosinus * newsin) + (RMin * sinus * newcos) ;
 
              /* trace (x,y) */
           }

Exemples :








[EDIT]

Données des graphes ci-dessus et petit prog d'exemple :

• Résultats gros fichier :
  

  Starts with a (380.159 370.736) center
  Processed 5170 points.
  Found RMaj 360.175 RMin 269.758 RAve 311.705
  Found ellipticty 0.251036 and angle 162.756 at center (x=395.952,y=372.379)
  with a sigma of 0.130287

• Résultats petit fichier :
  

  Starts with a (428.628 323.899) center
  Processed 744 points.
  Found RMaj 617.274 RMin 399.203 RAve 496.405
  Found ellipticty 0.353281 and angle 174.513 at center (x=388.07,y=366.353)
  with a sigma of 0.122929

• Code exemple (C) :
  

  Code :

  Sélectionner tout - Visualiser dans une fenêtre à part

  1
  2
  3
  4
  5
  6
  7
  8
  9
  10
  11
  12
  13
  14
  15
  16
  17
  18
  19
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
  100
  101
  102
  103
  104
  105
  106
  107
  108
  109
  110
  111
  112
  113
  114
  115
  116
  117
  118
  119
  120
  121
  122
  123
  124
  125
  126
  127
  128
  129
  130
  131
  132
  133
  134
  135
  136
  137
  138
  139
  140
  141
   
  /*
  =======================
   
  =======================
  */
  int main ( int argc, char **argv )
  {
    FILE     *in=(FILE *)NULL ;
    char      line[100] ;
    int       i, j, N = 0 ;
    double   *X=NULL, *Y=NULL, XCen, YCen, RMaj, RMin, RAve, Ellip, Angle, Sig  ;
    double    x1, y1, delta, newangle, sinus, cosinus, newcos, newsin ;
   
  /*
  --- Opens the file
  */
     if ( argc == 1 )
       {
         fprintf(stderr, "\nProgram should be called as : %s Filename",argv[0]);
   
         fprintf(stderr, "\nFilename containing the points as floating points, in 2 columns x y !!!!%c", 7 ); 
         return ERROR ;
       }
   
     in = fopen(argv[1], "r");   
     if (in == (FILE *)NULL)
       {
  	fprintf(stderr, "\nError when trying to open file %s ...", argv[1]);
  	return ERROR ;
       }
   
     /* Now load the data */	
     N = 0 ;
     while ( fgets(line, 100, in) != (char *)NULL )
       {
         if ( (line[0] != '#') && (line[0] != '/') )
  	 N = N + 1 ;
       }
   
     if ( N < 3 )
       {
         fclose(in);
         fprintf(stderr, "\nError.. Not enough points in file %s ...", argv[1]);
         return ERROR ;
       }
   
     if ( (X = (double *) calloc ( N, sizeof(double) )) != NULL )
       {
         if ( (Y = (double *) calloc ( N, sizeof(double) )) != NULL )
  	 {
  	   rewind(in);
   
  	   i = 0 ;
  	   while ( fgets(line, 100, in) != (char *)NULL )
  	     {
  	       if ( (line[0] != '#') && (line[0] != '/') )
  		 {
  		   sscanf ( line, "%lf %lf", &X[i], &Y[i] );
  		   i++ ;
  		 }
  	     }
  	 }
         else
  	 {
  	   free ( X);
  	   N = 0 ;
  	 }
       }
     else
       N = 0 ;
   
     fclose(in);
   
   
  /*
  ----------------------------------------------
   
    Appel de la routine de fit
   
  ----------------------------------------------
  */
     if ( N == 0 )
       {
         fprintf (stderr, "\nERROR !! A memory allocation error occured !! %c", 7);
         return ERROR ;
       }
   
     /* Computes the barycenter */
     XCen = 0.0 ;
     YCen = 0.0 ;
     for ( i = 0 ; i < N ; i++ )
       {
         XCen = XCen + X[i] ;
         YCen = YCen + Y[i] ;
       }
   
     XCen = XCen / (double)N ;
     YCen = YCen / (double)N ;
     fprintf ( stdout, "\n# Starts with a (%g %g) center", XCen,YCen);
     /* Calls the routine */
     if ( (j=FitEllipse ( X, Y, N, &XCen, &YCen, False, False, False,
  			&RMaj, &RMin, &RAve,
  			&Ellip, &Angle, &Sig )) == 1 )
       {
         fprintf ( stdout, "\n# Processed %d points.\n# Found RMaj %g RMin %g RAve %g",
  		 N,RMaj,RMin,RAve );
         fprintf ( stdout, "\n# Found ellipticty %g and angle %g at center (x=%g,y=%g)",
  		 Ellip,Angle,XCen,YCen);
         fprintf ( stdout, "\n# with a sigma of %g\n", Sig);
   
   
         /* Computing ellipse with computed parameters */
         delta = PI / 180.0 ;
         newangle = -Angle * PI / 180.0 ;
         newcos = cos (newangle) ;
         newsin = sin (newangle) ;
   
         for ( i = 0 ; i < 360 ; i++ )
  	 {
  	   cosinus = cos ((double)i*delta) ;
  	   sinus = sin ( (double)i*delta) ;
   
  	   x1 = XCen + (RMaj * cosinus * newcos) - (RMin * sinus * newsin) ;
  	   y1 = YCen + (RMaj * cosinus * newsin) + (RMin * sinus * newcos) ;
  	   fprintf ( stdout, "%g %g\n",x1,y1);
  	 }
   
         i = SUCCESS ;
       }
     else
       {
         fprintf ( stderr, "\nAN ERROR OCCURRED during the computation [error #%d]  !!!%c\n",j,7);
         i = ERROR ;
       }
   
     free ( Y );
     free ( X );
   
     return i ;
  }

[/EDIT]