Bonjour à tous,

j'ai un nuage de points que je veux modéliser par une surface. J'ai développé un algo qui fait une régression d'ordre 1(plan) et qui fonctionne. J'ai ensuite trouvé sur internet les formules pour la modélisation d'une surface d'ordre 2

http://support.articque.com/samples/doc/surf_t2.htm

Voici l'algo mais il ne marche pas

où est la faute?
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
function [a,b,c,d,e,f,g,h,i] = reg_ordre2(N,X,Y,Z)
 
Xi = 0;
Yi = 0;
Zi = 0;
Xi2 = 0;
Yi2 = 0;
Xi3 = 0;
Yi3 = 0;
Xi4 = 0;
Yi4 = 0;
 
XiYi = 0;
XiYi2 = 0;
XiYi3 = 0;
XiYi4 = 0;
 
 
Xi2Yi = 0;
Xi2Yi2 = 0;
Xi2Yi3 = 0;
Xi2Yi4 = 0;
 
Xi3Yi = 0;
Xi3Yi2 = 0;
Xi3Yi3 = 0;
Xi3Yi4 = 0;
 
Xi4Yi = 0;
Xi4Yi2 = 0;
Xi4Yi3 = 0;
Xi4Yi4 = 0;
 
XiZi = 0;
YiZi = 0;
XiYiZi = 0;
Xi2Zi = 0;
Yi2Zi = 0;
Xi2YiZi = 0;
XiYi2Zi = 0;
Xi2Yi2Zi = 0;
 
 
for i = 1:N
    Xi = Xi + N*X(i);
    Yi = Yi + N*Y(i);
    Xi2 = Xi2 + N*X(i)^2;
    Yi2 = Yi2 + N*Y(i)^2;
    Xi3 = Xi3 + N*X(i)^3;
    Yi3 = Yi3 + N*Y(i)^3;
    Xi4 = Xi4 + N*X(i)^4;
    Yi4 = Yi4 + N*Y(i)^4;
end
 
for i = 1:N
    for t = 1:N
        XiYi = XiYi + X(i)*Y(t);
        XiYi2 = XiYi2 + X(i)*Y(t)^2;
        XiYi3 = XiYi3 + X(i)*Y(t)^3;
        XiYi4 = XiYi4 + X(i)*Y(t)^4;
 
        Xi2Yi = Xi2Yi + X(i)^2*Y(t);
        Xi2Yi2 = Xi2Yi2 + X(i)^2*Y(t)^2;
        Xi2Yi3 = Xi2Yi3 + X(i)^2*Y(t)^3;
        Xi2Yi4 = Xi2Yi4 + X(i)^2*Y(t)^4;
 
        Xi3Yi = Xi3Yi + X(i)^3*Y(t);
        Xi3Yi2 = Xi3Yi2 + X(i)^3*Y(t)^2;
        Xi3Yi3 = Xi3Yi3 + X(i)^3*Y(t)^3;
        Xi3Yi4 = Xi3Yi4 + X(i)^3*Y(t)^4;
 
        Xi4Yi = Xi4Yi + X(i)^4*Y(t);
        Xi4Yi2 = Xi4Yi2 + X(i)^4*Y(t)^2;
        Xi4Yi3 = Xi4Yi3 + X(i)^4*Y(t)^3;
        Xi4Yi4 = Xi4Yi4 + X(i)^4*Y(t)^4;
 
        XiZi = XiZi + X(i)*Z(i,t);
        YiZi = YiZi + Y(t)*Z(i,t);
        XiYiZi = XiYiZi + Y(t)*X(i)*Z(i,t);
        Xi2YiZi = Xi2YiZi + Y(t)*X(i)^2*Z(i,t);
        XiYi2Zi = XiYi2Zi + Y(t)^2*X(i)*Z(i,t);
        Xi2Yi2Zi = Xi2Yi2Zi + Y(t)^2*X(i)^2*Z(i,t);
        Xi2Zi = Xi2Zi + X(i)^2*Z(i,t);
        Yi2Zi = Yi2Zi + Y(t)^2*Z(i,t);
        Zi = Zi + Z(i,t);
    end
end
 
P = [Xi2    XiYi   Xi2Yi  Xi3    XiYi2  Xi3Yi  Xi2Yi2 Xi3Yi2 Xi;
     XiYi   Yi2    XiYi2  Xi2Yi  Yi3    Xi2Yi2 XiYi3  Xi2Yi3 Yi;
     XiYi   XiYi2  Xi2Yi2 Xi3Yi  XiYi3  Xi3Yi2 Xi2Yi3 Xi3Yi3 XiYi
     Xi3    Xi2Yi  Xi3Yi  Xi4    Xi2Yi2 Xi4Yi  Xi3Yi2 Xi4Yi2 Xi2;
     XiYi2  Yi3    XiYi3  Xi2Yi2 Yi4    Xi2Yi3 XiYi4  Xi2Yi4 Yi2;
     Xi3Yi  Xi2Yi2 Xi3Yi2 Xi4    Xi2Yi3 Xi4Yi2 Xi3Yi3 Xi4Yi3 Xi2Yi;
     Xi2Yi2 XiYi3  Xi2Yi3 Xi3Yi2 XiYi4  Xi3Yi3 Xi2Yi4 Xi3Yi4 XiYi2;
     Xi3Yi2 Xi2Yi3 Xi3Yi3 Xi4Yi2 Xi2Yi4 Xi4Yi3 Xi3Yi4 Xi4Yi4 Xi2Yi2;
     Xi     Yi     XiYi   Xi2    Yi2    Xi2Yi  XiYi2  Xi2Yi2 N^2];
 
V = [XiZi;
     YiZi;
     XiYiZi;
     Xi2Zi;
     Yi2Zi;
     Xi2YiZi;
     XiYi2Zi;
     Xi2Yi2Zi;
     Zi];
 
R = inv(P)*V;
 
a = R(1);
b = R(2);
c = R(3);
d = R(4);
e = R(5);
f = R(6);
g = R(7);
h = R(8);
i = R(9);
 
return