Salut,
Pour ma première contribution (c++), je vous propose une régression circulaire sur un ensemble de points (x; y).
Voilà la déclaration de la classe Matrix (matrix.h) :
et voilà l'implémentation de la classe Matrix (matrix.cpp) :
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 #ifndef __MATRIX__ #define __MATRIX__ #include <vector> struct Matrix { /// @brief Solve the system AX=Y where A and Y are known and X is unknown. /// @param [in] A A n*n matrix where n = A.size() / 2 = Y.size(). /// @param [in] X The matrix to find. X is resized if needed. /// @param [in] Y A 1*n matrix where n = A.size() / 2 = Y.size(). /// @throw std::runtime_error If an error occured. static void gaussjordan(std::vector<double> const & A, std::vector<double> & X, std::vector<double> const & Y); /// @brief Process to the regression of the points (xlist; ylist) by a circle. /// @param [in] xlist List of the x coordinates of the points that should fit a circle. /// @param [in] ylist List of the y coordinates of the points that should fit a circle. /// @param [in] xcenter The x coordinate of the center of the circle computed from the points (xlist; ylist). /// @param [in] ycenter The y coordinate of the center of the circle computed from the points (xlist; ylist). /// @param [in] radius The radius of the circle computed from the points (xlist; ylist). /// @param [in] rmse The root mean square error of the regression computed as the root squared sum of the difference between the computed radius and the distances between the computed center and the points (xlist; ylist) /// @throw std::runtime_error If an error occured. static void circularRegression(std::vector<double> const & xlist, std::vector<double> const & ylist, double & xcenter, double & ycenter, double & radius, double & rmse); }; // Matrix #endif
Le tout s'utilise très facilement (en c++ je rappelle) via le code :
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 #include <limits> #include <stdexcept> #include <utility> #include <map> #include <cmath> #include "matrix.h" void Matrix::gaussjordan(std::vector<double> const & A, std::vector<double> & X, std::vector<double> const & Y) { // get the size of the matrix and check system input unsigned int const n = Y.size(); unsigned int const m = n + 1; if(A.size() != n * n) { throw std::runtime_error("Bad size with A and Y arguments"); } // resize the ouptut if(X.size() != n) { X.resize(n); } std::map<unsigned int, unsigned int> swapped; for(unsigned int row = 0; row < n; ++row) { swapped[row] = row; } // here is the system to solve std::vector<double> system(n * n + n); // array of pointers on each row of the system std::vector<double *> rows(n); // setup the system std::vector<double>::const_iterator a = A.begin(); std::vector<double>::const_iterator y = Y.begin(); std::vector<double>::iterator s = system.begin(); for(unsigned int row = 0; row < n; ++row) { // keep the address of the first item of the current row rows[row] = &(*s); // set A for(unsigned int col = 0; col < n; ++col, ++a, ++s) { *s = *a; } // set Y *s++ = *y++; } // find X to solve the system AX=Y for(unsigned int ipivot = 0; ipivot < n; ++ipivot) { double pivot = rows[ipivot][ipivot]; // if the current pivot is null, search for another line with a valid pivot (for the column ipivot) if(fabs(pivot) <= std::numeric_limits<double>::epsilon()) { // browse next row, searching for a valid pivot at column ipivot for(unsigned int row = ipivot + 1; row < n; ++row) { // if found, swap the rows pointer (not the memory) pivot = rows[row][ipivot]; if(fabs(pivot) > std::numeric_limits<double>::epsilon()) { std::swap(rows[ipivot], rows[row]); std::swap(swapped[ipivot], swapped[row]); break; } } } // if a valid pivot was not found, then the matrix can not be inversed if(fabs(pivot) <= std::numeric_limits<double>::epsilon()) { throw std::runtime_error("The system can not be solved"); } // normalize the row of the pivot so as the pivot is set to 1.0 for(unsigned int col = 0; col < m; ++col) { rows[ipivot][col] /= pivot; } // for all row, except the one of the pivot, set the item of the column "ipivot" ot 0 for(unsigned int row = 0; row < n; ++row) { if(row == ipivot) { continue; } double const factor = rows[row][ipivot]; for(unsigned int col = 0; col < m; ++col) { rows[row][col] -= factor * rows[ipivot][col]; } } } // now the system is solved, X is the last column of the system for(unsigned int row = 0; row < n; ++row) { X[row] = system[swapped[row] * m + m - 1]; } } void Matrix::circularRegression(std::vector<double> const & xlist, std::vector<double> const & ylist, double & xcenter, double & ycenter, double & radius, double & rmse) { xcenter = 0.0; ycenter = 0.0; radius = 0.0; rmse = 0.0; // compute the A (3x3 matrix) and B (3x1 matrix) std::vector<double> A(9); std::vector<double> B(3); std::fill(A.begin(), A.end(), 0); std::fill(B.begin(), B.end(), 0); for(unsigned int index = 0; index < xlist.size(); ++index) { double const & x = xlist[index]; double const & y = ylist[index]; double const xx = x * x; double const yy = y * y; double const xy = x * y; double const xxx = xx * x; double const xxy = xx * y; double const yyy = yy * y; double const xyy = x * yy; A[0] += 2 * xx; A[1] += 2 * xy; A[2] += x; A[3] += 2 * xy; A[4] += 2 * yy; A[5] += y; A[6] += 2 * x; A[7] += 2 * y; A[8] += 1; B[0] += xxx + xyy; B[1] += xxy + yyy; B[2] += xx + yy; } // solve the system AX=B using a the gauss-jordan method std::vector<double> X(3); gaussjordan(A, X, B); // circle details extraction xcenter = X[0]; ycenter = X[1]; radius = sqrt(X[2] + xcenter * xcenter + ycenter * ycenter); // compute the root mean square for(unsigned int index = 0; index < xlist.size(); ++index) { double const & x = xlist[index]; double const & y = ylist[index]; double const distance = sqrt((x - xcenter) * (x - xcenter) + (y - ycenter) * (y - ycenter)); rmse += (radius - distance) * (radius - distance); } rmse /= (double)xlist.size(); rmse = sqrt(rmse); }
Et un main (main.cpp) pour tester tout ça :
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 // variables std::vector<double> xlist; // ... The list of the x coordinates of the points that match a circle. std::vector<double> ylist; // ... The list of the y coordinates of the points that match a circle. double xcenter = 0.0; // To store the x coordinate of the center of the circle computed by the function. double ycenter = 0.0; // To store the y coordinate of the center of the circle computed by the function. double radius = 0.0; // To store the radius of the circle computed by the function. double rmse = 0.0; // To store the root mean square error resulted from the regression. try { // regression of the points by a circle Matrix::circularRegression(xlist, ylist, xcenter, ycenter, radius, rmse); } catch(std::exception const & error) { std::cout << error.what() << std::endl; } catch(...) { std::cout << "An undefined error has occured" << std::endl; }
Le code ne requiert aucune autre dépendance que la Standart Template Library (native sur tous les compilateurs), pas de boost, opencv ou autres. Vous compilez sous gcc par la commande suivante :
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 #include <iostream> #include <iomanip> #include <cmath> #include <stdio.h> #include <stdlib.h> #include "matrix.h" struct Circle { double xcenter; double ycenter; double radius; }; // Circle int main() { // for std::stw unsigned int const setwv = 5; // definition of the circle Circle const expected = {20.5, -3.7, 8.0}; Circle computed = {0.0, 0.0, 0.0}; // lists of the points on the circle used to make the regression std::vector<double> xlist; std::vector<double> ylist; // display the informations about the expected circle std::cout << " ---------------------------------------------------------------------"; std::cout << std::endl; std::cout << "||" << std::right << std::setw(2 * setwv) << "Expected circle : " << std::endl; std::cout << "||" << std::right << std::setw(4 * setwv) << "x center = " << expected.xcenter << std::endl; std::cout << "||" << std::right << std::setw(4 * setwv) << "y center = " << expected.ycenter << std::endl; std::cout << "||" << std::right << std::setw(4 * setwv) << "radius = " << expected.radius << std::endl; // display the titles of the informations about the points on the circle std::cout << " ---------------------------------------------------------------------"; std::cout << std::endl; std::cout << "||" << std::right << std::setw(2 * setwv) << "x original" << "|" << std::right << std::setw(2 * setwv) << "x noised" << "|" << std::right << std::setw(2 * setwv) << "noise" << std::right << std::setw(5) << "||" << std::right << std::setw(2 * setwv) << "y original" << "|" << std::right << std::setw(2 * setwv) << "y noised" << "|" << std::right << std::setw(2 * setwv) << "noise" << "||"; std::cout << std::endl; std::cout << " ---------------------------------------------------------------------"; std::cout << std::endl; // compute the points used to make the regression for(unsigned int angle = 0; angle < 360; angle += 12) { double const radian = 2 * 3.1415 * (double)angle / 360.0; double const x = expected.xcenter + expected.radius * cos(radian); double const y = expected.ycenter + expected.radius * sin(radian); double const xnoise = 2.0 * (double)(rand() % 200 - 100) / 100.0; double const ynoise = 2.0 * (double)(rand() % 200 - 100) / 100.0; xlist.push_back(x + xnoise); ylist.push_back(y + ynoise); // values std::cout << "||" << std::right << std::setw(2 * setwv) << x<< "|" << std::right << std::setw(2 * setwv) << xlist.back() << "|" << std::right << std::setw(2 * setwv) << xnoise << std::right << std::setw(5) << "||" << std::right << std::setw(2 * setwv) << y << "|" << std::right << std::setw(2 * setwv) << ylist.back() << "|" << std::right << std::setw(2 * setwv) << ynoise << "||"; std::cout << std::endl; } double rmse = 0.0; try { // regression of the points by a circle Matrix::circularRegression(xlist, ylist, computed.xcenter, computed.ycenter, computed.radius, rmse); } catch(std::exception const & error) { std::cout << error.what() << std::endl; return 0; } catch(...) { std::cout << "An undefined error has occured" << std::endl; return 0; } // display the informations about the computed circle std::cout << " ---------------------------------------------------------------------"; std::cout << std::endl; std::cout << "||" << std::right << std::setw(2 * setwv) << "Computed circle : " << std::endl; std::cout << "||" << std::right << std::setw(4 * setwv) << "x center = " << computed.xcenter << std::endl; std::cout << "||" << std::right << std::setw(4 * setwv) << "y center = " << computed.ycenter << std::endl; std::cout << "||" << std::right << std::setw(4 * setwv) << "radius = " << computed.radius << std::endl; std::cout << "||" << std::right << std::setw(4 * setwv) << "rmse = " << rmse << std::endl; std::cout << " ---------------------------------------------------------------------"; std::cout << std::endl; return 0; }
g++ main.cpp matrix.cpp -o test
La sortie est du style :
Bon je veux pas faire de la pub pour mon blog, mais je veux pas non plus ré-expliquer le fonctionnement du code. Le lien vers des explications détaillées est le suivant :
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 --------------------------------------------------------------------- ||Expected circle : || x center = 20.5 || y center = -3.7 || radius = 8 --------------------------------------------------------------------- ||x original| x noised| noise ||y original| y noised| noise|| --------------------------------------------------------------------- || 28.5| 27.32| -1.18 || -3.7| -4.36| -0.66|| || 28.3252| 29.0052| 0.68 || -2.03675| -2.03675| 0|| || 27.8084| 29.1884| 1.38 || -0.446197| 0.0338029| 0.48|| || 26.9722| 26.5322| -0.44 || 1.00216| 2.16216| 1.16|| || 25.8532| 27.0932| 1.24 || 2.24503| 1.52503| -0.72|| || 24.5002| 24.6002| 0.1 || 3.22808| 4.12808| 0.9|| || 22.9724| 22.5924| -0.38 || 3.90836| 2.44836| -1.46|| || 21.3366| 22.5566| 1.22 || 4.25614| 4.07614| -0.18|| || 19.6642| 21.5642| 1.9 || 4.25622| 5.09622| 0.84|| || 18.0283| 16.5683| -1.46 || 3.90859| 2.62859| -1.28|| || 16.5004| 18.3204| 1.82 || 3.22845| 1.30845| -1.92|| || 15.1474| 15.1874| 0.04 || 2.24552| 3.30552| 1.06|| || 14.0282| 13.8682| -0.16 || 1.00276| 2.64276| 1.64|| || 13.1919| 11.6119| -1.58 || -0.44552| -0.12552| 0.32|| || 12.675| 13.035| 0.36 || -2.03603| -2.13603| -0.1|| || 12.5| 11.44| -1.06 || -3.69926| -3.17926| 0.52|| || 12.6747| 14.0947| 1.42 || -5.36252| -4.60252| 0.76|| || 13.1913| 12.5713| -0.62 || -6.95313| -6.71313| 0.24|| || 14.0273| 13.3673| -0.66 || -8.40156| -8.42156| -0.02|| || 15.1463| 13.8463| -1.3 || -9.64453| -9.76453| -0.12|| || 16.4991| 16.5591| 0.06 || -10.6277| -12.4077| -1.78|| || 18.0269| 18.4669| 0.44 || -11.3081| -10.6481| 0.66|| || 19.6627| 19.1227| -0.54 || -11.6561| -12.3761| -0.72|| || 21.3351| 22.1551| 0.82 || -11.6563| -11.4363| 0.22|| || 22.971| 22.031| -0.94 || -11.3088| -11.9488| -0.64|| || 24.4989| 25.4389| 0.94 || -10.6288| -11.7488| -1.12|| || 25.8521| 25.0921| -0.76 || -9.64602| -8.50602| 1.14|| || 26.9714| 25.7114| -1.26 || -8.40336| -9.22336| -0.82|| || 27.8078| 28.2678| 0.46 || -6.95516| -6.13516| 0.82|| || 28.3249| 28.9049| 0.58 || -5.3647| -3.8047| 1.56|| --------------------------------------------------------------------- ||Computed circle : || x center = 20.3688 || y center = -3.76401 || radius = 8.25942 || rmse = 1.01834 ---------------------------------------------------------------------
http://floriansella.free.fr/index.ph...y=ti&billid=92
Histoire de pas faire un truc creux, la régression d'un ensemble de points par un cercle (xc; yc; R) fonctionne comme la régression d'un ensemble de points par une droite. C'est une régression au sens des moindres carrés et on parvient à un système matriciel qu'il suffit de résoudre par la méthode basique de Gauss-Jordan (pivot de Gauss).
Flo.
Partager