[FORTRAN] Ellipse fitting algorithm

[EDIT]
OK je vais présenter un peu mieux..

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/... . 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'horzontale) 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.

[/EDIT]

OK.. J'ai retrouvé le code.. Mais c'est du Fortran...

J'essayerais de le traduire en C un de ces jours... Et de donner un exemple d'utilisation

(j'ai gardé la bibliothèque, mais perdu le prog de départ.. Je crois que j'ai une sortie papier quelque part, mais j'aurais pas accès avant un bon mois..)

D'abord les 2 routines de mahs nécessaires (moindres carrés généralisés à 10 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
 
C***********************************************************************C
C         GENERALIZED LEAST-SQUARE FOR UP TO 10 PARAMETERS              C
C                     (Jean SOUVIRON - DAO Victoria - 1985)             C
C        (METHOD OF R. SEDGEWICK in ALGORITHMS                          C
C                                   ADDISON-WESLEY PUBLISHING CO.       C
C                                   1984                             )  C
C-----------------------------------------------------------------------C
C                                                                       C
C     INPUT :                                                           C
C           NARAY = NUMBER OF DIFFERENT PARAMETERS+1                    C
C            NDIM = REAL MAX DIMENSION IN THE MAIN PROGRAM              C
C          NTOTAL = TOTAL NUMBER OF POINTS                              C
C          TABOBS = ARRAY OF PARAMETERS FOR EACH POINT                  C
C                   (LAST VECTOR HAS TO BE THE DATA TO BE FITTED)       C
C     OUTPUT :                                                          C
C          FACTOR = ARRAY OF COEFFICIENTS                               C
C                       (DIMENSION MUST BE BUILD IN IN MAIN PROGRAM)    C
C          SIGMA                                                        C
C          IERR   ERROR FLAG                                            C
C                                                                       C
C***********************************************************************C
        SUBROUTINE LEAST(TABOBS,NDIM,NARAY,NTOTAL,FACTOR,SIGMA,IERR)
        DIMENSION TABOBS(NDIM,NARAY),FACTOR(1),COEFF(10,11)
        DIMENSION SOLU(10)
        DOUBLE PRECISION COEFF,SOLU
        DO 757 J=1,10
        SOLU(J)=0.D0
        DO 756 I=1,11
        COEFF(J,I)=0.D0
 756    CONTINUE
 757    CONTINUE
        SIGMA=0.
        NREAL=NARAY-1
        DO 1 I=1,NREAL
        DO 2 J=1,NARAY
        T=0.0
        DO 3 K=1,NTOTAL
        Z1=TABOBS(K,I)
        Z2=TABOBS(K,J)
        T=T+Z1*Z2
3       CONTINUE
        COEFF(I,J)=DBLE(T)
2       CONTINUE
1       CONTINUE
C-----------------
C   SOLVING BY GAUSSIAN ELIMINATION
C-----------------
        CALL LNSYS(COEFF,SOLU,NREAL,IERR)
C-----------------
C  SET UP COEFFICIENTS
C-----------------
        DO 4 I=1,NREAL
        FACTOR(I)=SOLU(I)
4       CONTINUE
        DO 7 I=1,NTOTAL
        ZZ=TABOBS(I,NARAY)
        DO 8 J=1,NREAL
        ZZ=ZZ-FACTOR(J)*TABOBS(I,J)
8       CONTINUE
        IF(ZZ.LE.0.)ZZ=-ZZ
        SIGMA=SIGMA+ZZ
7       CONTINUE
        SIGMA=(SIGMA/FLOATJ(NTOTAL))*SQRT(3.14159/2.)
        RETURN
        END
 
C***********************************************************C
C       SOLVE LINEAR SYSTEM VIA GAUSSIAN ELIMINATION        C
C          ORDER UP TO 10                                   C
C             ( 10.01.1979 K. BANSE )                       C
C                              (SACLAY - HPCCD - 1982/84)   C
C***********************************************************C
        SUBROUTINE LNSYS(ARRAY,X,NORDER,IERR)
        DIMENSION ARRAY(10,11), X(10)
        DOUBLE PRECISION ARRAY, SAVE,RMAX,X
C
C
        IERR=0
        DO 100 K=1,NORDER-1
C--------------------------------- GET LARGEST PIVOTAL ELEMENT
        SAVE=0.D0
        DO 30 J=K,NORDER
        RMAX=DABS(ARRAY(J,K))
        IF(SAVE.GE.RMAX)GO TO 30
        SAVE=RMAX
        MXROW=J
30      CONTINUE
C--------------------------------- EXCHANGE ROWS
        DO 40 J=K,NORDER+1
        SAVE=ARRAY(K,J)
        ARRAY(K,J)=ARRAY(MXROW,J)
        ARRAY(MXROW,J)=SAVE
40      CONTINUE
C--------------------------------- TEST IF ARRAY SINGULAR
        IF(DABS(ARRAY(K,K)).LT.10.D-12) GO TO 300
C--------------------------------- SUBSTRACT ROW FROM LOWER ROWS
        DO 60 I=K+1,NORDER
        SAVE=ARRAY(I,K)/ARRAY(K,K)
        ARRAY(I,K)=0.D0
        DO 61 J=K+1,NORDER+1
        ARRAY(I,J)=ARRAY(I,J)-SAVE*ARRAY(K,J)
61      CONTINUE
60      CONTINUE
100     CONTINUE
C
C       CALCULATE VECTOR X       
C
        X(NORDER)=ARRAY(NORDER,NORDER+1)/ARRAY(NORDER,NORDER)
        DO 200 NN=2,NORDER
        I=NORDER+1-NN
        SAVE=ARRAY(I,NORDER+1)
        DO 150 J=I+1,NORDER
        SAVE=SAVE-ARRAY(I,J)*X(J)
150     CONTINUE
        X(I)=SAVE/ARRAY(I,I)
200     CONTINUE
        GO TO 333
C
C    
300     DO 310 N=1,NORDER
        X(N)=0.0
310     CONTINUE
        IERR=-1
C
333     RETURN 
        END

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
 
C*******************************************************************C
C               SUBROUTINE OF ELLIPSES FITTING                      C
C                                                                   C
C            IT TAKES THE GENERAL EQUATION THAT AN ELLIPSE          C
C        HAS TO SATISFY AND FIT BY LEAST-SQUARE THE COEFFICIENTS    C
C        IF THE AXIS WERE ALIGNED WITH X AND Y , THE EQUATION       C
C        WOULD BE :   (X-XC)**2/A2 + (Y-YC)**2/B2 = 1               C
C        WITH A CERTAIN ANGLE ALPH THE COORDINATES TRANSFORMS       C
C        INTO :   X' = COS (APLH) * (X-XC) - SIN (ALPH) * (Y-YC)    C
C                 Y' = SIN (ALPH) * (X-XC) + COS (ALPH) * (Y-YC)    C
C        REPLACING IN THE FORMER EQUATION GIVES AN EQUATION WHOSE   C
C        FORM IS :   A (X-XC)**2 + B (Y-YC)**2 + C X*Y = 1          C
C        A VERY SIMPLE CALCULUS GIVES THE ANSWER FOR THE VALUE      C
C        OF THE DIFFERENT PARAMETERS WHEN THE COEFFICIENTS A,B,C    C
C        ARE DETERMINED BY LEAST-SQUARE.                            C
C                                                                   C
C        THIS ROUTINE ALLOWS THE USER TO FIT ELLIPSES WHERE :       C
C             - THE CENTER IS EITHER KNOWN OR UNKNOWN               C
C             - THE POSITION ANGLE IS EITHER KNOWN OR UNKNOWN       C
C             - THE ELLIPTICITY IS EITHER KNONW OR UNKNOWN          C
C                                                                   C
C        IT IS LIMITED TO 10000 POINTS .                            C
C                                                                   C
C                                                                   C
C    INPUT :                                                        C
C        XX : TABLE OF X COORDINATES                                C
C        YY : TABLE OF Y COORDINATES                                C
C        XCENTR,YCENTR : COORDINATES OF CENTER                      C
C        NPOIN : NUMBER OF POINTS THAT HAVE TO BE FITTED            C
C        INDIC : FLAG FOR DETERMINATION OF CENTER (0) OR NO (1)     C
C        INDIA : FLAG FOR DETERMINATION OF ANGLE (0) OR NO (1)      C
C                (in case of no determination, the value of the     C
C                angle has to be passed by ANG)                     C
C        INDIE : FLAG FOR DETERMINATION OF ELLIPTICITY (0) OR NO(1) C
C                (in case of no determination, the value of the     C
C                ellipticity has to be passed by ELL)               C
C                                                                   C
C    OUTPUT :                                                       C
C        RMAJ,RAVE,ELL,ANG : PARAMETERS OF THE ELLIPSE              C
C        XCENTR,YCENTR : EVENTUALLY NEW COORDINATES CENTER          C
C        SIGM,IERR : SIGMA AND ERROR FLAG                           C
C                                                                   C
C        ( Jean SOUVIRON - DAO Victoria 1985/86)                    C
C                                                                   C
C*******************************************************************C
      SUBROUTINE FITELL(XX,YY,XCENTR,YCENTR,NPOIN,INDIC,INDIA,INDIE,
     +RMAJ,RAVE,ELL,ANG,SIGM,IERR)
      DOUBLE PRECISION F1,F2,F3,G1,G2
      DIMENSION TABFIT(10000,4),XX(1),YY(1),FACTO(3) 
      DIMENSION TABFI1(10000,6),FACTO1(5)
      IERR=0
      IF(INDIC.EQ.1)GO TO 1000
C-----------------------------C
C  CENTER HAS TO BE FOUND     C
C-----------------------------C
200   CORRX=0.
      CORRY=0.
      DO 2500 KKK=1,NPOIN
      TABFI1(KKK,1)=(XX(KKK)-XCENTR)*(XX(KKK)-XCENTR)
      TABFI1(KKK,2)=(YY(KKK)-YCENTR)*(YY(KKK)-YCENTR)
      TABFI1(KKK,3)=(XX(KKK)-XCENTR)*(YY(KKK)-YCENTR)
      TABFI1(KKK,4)=XX(KKK)-XCENTR
      TABFI1(KKK,5)=YY(KKK)-YCENTR
      TABFI1(KKK,6)=1.
2500  CONTINUE
      CALL LEAST(TABFI1,10000,6,NPOIN,FACTO1,ASIGM,IER)
      IF(IER.LT.0)GO TO 5
      DENUM=4.*FACTO1(1)*FACTO1(2)-FACTO1(3)*FACTO1(3)
      IF(DENUM.EQ.0.)GO TO 3
      CORRX=(-2.*FACTO1(2)*FACTO1(4)+FACTO1(3)*FACTO1(5))/DENUM
      XCENTR=XCENTR+CORRX
      CORRY=(FACTO1(3)*FACTO1(4)-2.*FACTO1(1)*FACTO1(5))/DENUM
      YCENTR=YCENTR+CORRY
C
C     CORRECTS THE FACTORS TO TAKE INTO ACCOUNT THE BAD CENTER
C
       CORREC=FACTO1(1)*CORRX*CORRX+FACTO1(2)*CORRY*CORRY
       CORREC=CORREC+FACTO1(3)*CORRX*CORRY
       IF((1.+CORREC).EQ.0.)GO TO 4
       COR=1.-(CORREC/(1.+CORREC))
       FACTO(1)=FACTO1(1)*COR
       FACTO(2)=FACTO1(2)*COR
       FACTO(3)=FACTO1(3)*COR
       GO TO 2000
C-----------------------------C
C  CENTER ALREADY DETERMINED  C
C-----------------------------C
1000  CONTINUE
      DO 4440 KKK=1,NPOIN
      TABFIT(KKK,1)=(XX(KKK)-XCENTR)*(XX(KKK)-XCENTR)
      TABFIT(KKK,2)=(YY(KKK)-YCENTR)*(YY(KKK)-YCENTR)
      TABFIT(KKK,3)=(XX(KKK)-XCENTR)*(YY(KKK)-YCENTR)
      TABFIT(KKK,4)=1.
4440  CONTINUE
      CALL LEAST(TABFIT,10000,4,NPOIN,FACTO,ASIGM,IER)
      IF(IER.LT.0)GO TO 5
C-----------------------------------------C
C     DETERMINATION OF THE PARAMETERS     C
C-----------------------------------------C
2000  CONTINUE
      F1=DBLE(FACTO(1))
      F2=DBLE(FACTO(2))
      F3=DBLE(FACTO(3))
      FACT=180./3.14159
      SIG=ASIGM
      IF(INDIA.EQ.1)GO TO 5000
      IF(INDIE.EQ.1)GO TO 6000
C************************************************************  NORMAL CASE
C
C--------------   COMPONENT ON x-AXIS
C
      G1=F3*F3+(F1-F2)*(F1-F2)
      IF(G1.LE.0.D0)GO TO 1
      G2=(2./(F1+F2+DBLE(SQRT(SNGL(G1)))))
      ZZ1=SNGL(G2)
      IF(ZZ1.LE.0.)GO TO 1
      RMAJ=SQRT(ZZ1)
      A=RMAJ
C
C---------------   COMPONENT ON y-AXIS
C
      ZZ=SNGL(F1+F2-DBLE(1./(RMAJ*RMAJ)))
      IF(ZZ.LE.0.)GO TO 2
      RMIN=1./SQRT(ZZ)
      B=RMIN
      IREVER=0
      IF(RMAJ.GE.RMIN)GO TO 4419
      IREVER=1
      TT=RMAJ
      RMAJ=RMIN
      RMIN=TT
C
C---------------     ELLIPTICITY , AND MEAN AXIS
C
4419  ELL=1.-RMIN/RMAJ
      RAVE=SQRT(RMIN*RMAJ)
C
C---------------    ANGLE
C
C  1)    x/HORIZONTAL
C    SIGN - BECAUSE IT IS THE REVERSE THAT WE WANT (FROM NORMAL TO ELLIPSE)
C
      G1=-F3/(F1-F2+DBLE(1./(A*A))-DBLE(1./(B*B)))
      ARGUM=SNGL(G1)
      ANG=ATAN(ARGUM)*FACT
C
C  2)  CHECK IF x IS MAJOR AXIS
C
      IF(IREVER.EQ.0)GO TO 10
      IF(ANG.GE.0.)GO TO 11
      ANG=ANG+90.
      GO TO 10
11    ANG=ANG-90.
C
C  3)  ANGLE MAJOR AXIS/VERTICAL, UNCLOCKWISE
C
10    ANG=ANG+90.
      RETURN
C********************************************************   CASE OF KNOWN ANGLE
C
C
5000  ANGL=-(ANG-90.)/FACT
      SINU=2.*SIN(ANGL)*COS(ANGL)
      IF(SINU.EQ.0.)SINU=0.000001
      G1=F1+F2
      G2=F1-F2
      ZZ=SNGL(F3/DBLE(SINU)+G1)
      IF(ZZ.LE.0.)GO TO 6
      AA=SQRT(2./ZZ)
      RMAJ=AA
      ZZ=SNGL(G1-F3/DBLE(SINU))
      IF(ZZ.LE.0.)GO TO 6
      BB=SQRT(2./ZZ)
      RMIN=BB
      IF(RMAJ.GE.RMIN)GO TO 5001
      TT=RMAJ
      RMAJ=RMIN
      RMIN=TT
5001  ELL=1.-RMIN/RMAJ
      RAVE=SQRT(RMIN*RMAJ)
      RETURN
C*************************************************   CASE OF KNOWN ELLIPTICITY
C
C
6000  BOVERA=1.-ELL
      G1=F1+F2
      IF(G1.LE.0.D0)GO TO 7
      RMIN=SQRT((1.+BOVERA*BOVERA)/SNGL(G1))
      RMAJ=RMIN/BOVERA
      RAVE=SQRT(RMAJ*RMIN)
      ZZ=SNGL((F1-F2)+DBLE((BOVERA*BOVERA-1.)/(RMIN*RMIN)))
      IF(ZZ.EQ.0.)GO TO 7
      ARGUM=-(SNGL(F3))/ZZ
      ANG=ATAN(ARGUM)*FACT
      ANG=ANG+90.
      RETURN
C
C----------- ERRORS
C
1     IERR=-1
      RETURN
2     IERR=-2
      RETURN
3     IERR=-3
      RETURN
4     IERR=-4
      RETURN
5     IERR=-5
      RETURN
6     IERR=-6
      RETURN
7     IERR=-7
      RETURN
      END

Bon désolé c'est du Fortran.... On verra plus tard..

Citation:

Et de donner un exemple d'utilisation

Avant de donner un exemple d'utilisation, il aurait été bon de dire exactement ce qu'est l'ellipse fitting, parce que mis à part ceux qui connaissent ce que c'est, je ne suis pas convaincu que les autres le savent :D . Juste dire qu'à partir d'un groupe de points, ça trouve la meilleure ellipse englobante.

Citation:

Je crois que j'ai une sortie papier quelque part, mais j'aurais pas accès avant un bon mois..

Un mois ! A quand même :aie:

Tu as une idée de la complexité ?

Citation:

Envoyé par PRomu@ld

Avant de donner un exemple d'utilisation, il aurait été bon de dire exactement ce qu'est l'ellipse fitting, parce que mis à part ceux qui connaissent ce que c'est, je ne suis pas convaincu que les autres le savent :D . Juste dire qu'à partir d'un groupe de points, ça trouve la meilleure ellipse englobante.

OK Je vais ajouter ça..

Citation:

Envoyé par PRomu@ld

Tu as une idée de la complexité ?

Non, du tout.. Mais ce doit être en O(N) + chouia...

En fait, il y a la résolution du système linéaire, mais sinon c'est assez linéaire...

Salut,

pour la traduction de "ellipse fitting", le terme "ajustement d'ellipse" est peut-etre plus adapte, non ?

Est ce que se serait possible de detailler plus les grandes lignes de l'algo ? J'imagine qu'on minimise une certaine fonction, mais laquelle ?
J'imagine que la solution n'st pas unique, comment est ce que c'est gere par la methode ?

Citation:

Envoyé par Kangourou

J'imagine qu'on minimise une certaine fonction, mais laquelle ?

Tout est dans l'en-tête. Comme j'ai dit je ne pourrais pas donner plus de détails avant la semaine prochaine ou la suivante.

Mais c'est TRES SIMPLE :

Prendre un pixel "considéré" comme appartenant à une ellipse, de demi grand axe a et de demi petit axe b, tournée d'un angle alpha par rapport à l'horizontale, et dont le centre est (Xc,Yc).

Sans la rotation on a
Code:

1 2 3 4 (x-xc)^2 + (y-yc)^2 --------- -------- = 1 a^2 b^2
Développé, ça donne x2/a2 + y2/b2 + ....

Donc le principe est de calculer et stocker pour tout pixel les x2, y2, xy, etc..
et d'appliquer les moindres carrés sur cela.

Les coefficients donnent le résultat.

Je scannerais mon papier à ce sujet comme je le disais soit semaine prochaine ou suivante..

Et normalement la résolution devrait être unique , car système linéaire avec autant d'équations que d'inconnues.

[EDIT]
Ai remplacé "moins d'équations" par "autant d'équations" ...
[/EDIT]