[C] Ellipse-fitting algorithm / Algorithme d'alignement d'ellipse

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 :oops:

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" :D).. 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:

---

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:

---

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:

---

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 :


Pièce jointe 197265

Pièce jointe 197266



[EDIT]

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

• Résultats gros fichier :
  
  Citation:

  ---

  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 :
  
  Citation:

  ---

  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:

  ---

  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]