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
|
// File: AppParCurves_Variational_3.gxx
// Created: Tue Dec 8 09:08:54 1998
// Author: Igor FEOKTISTOV
// <ifv@paradox.nnov.matra-dtv.fr>
void AppParCurves_Variational::Project(const Handle(FEmTool_Curve)& C,
const TColStd_Array1OfReal& Ti,
TColStd_Array1OfReal& ProjTi,
TColStd_Array1OfReal& Distance,
Standard_Integer& NumPoints,
Standard_Real& MaxErr,
Standard_Real& QuaErr,
Standard_Real& AveErr,
const Standard_Integer NbIterations) const
{
// Initialisation
const Standard_Real Seuil = 1.e-9, Eps = 1.e-12;
MaxErr = QuaErr = AveErr = 0.;
Standard_Integer Ipnt, NItCv, Iter, i, i0 = -myDimension, d0 = Distance.Lower() - 1;
Standard_Real TNew, Dist, T0, Dist0, F1, F2, Aux, DF, Ecart;
Standard_Boolean EnCour;
TColStd_Array1OfReal ValOfC(1, myDimension), FirstDerOfC(1, myDimension),
SecndDerOfC(1, myDimension);
for(Ipnt = 1; Ipnt <= ProjTi.Length(); Ipnt++) {
i0 += myDimension;
TNew = Ti(Ipnt);
EnCour = Standard_True;
NItCv = 0;
Iter = 0;
C->D0(TNew, ValOfC);
Dist = 0;
for(i = 1; i <= myDimension; i++) {
Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
Dist += Aux * Aux;
}
Dist = Sqrt(Dist);
// ------- Newton's method for solving (C'(t),C(t) - P) = 0
while( EnCour ) {
Iter++;
T0 = TNew;
Dist0 = Dist;
C->D2(TNew, SecndDerOfC);
C->D1(TNew, FirstDerOfC);
F1 = F2 = 0.;
for(i = 1; i <= myDimension; i++) {
Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
DF = FirstDerOfC(i);
F1 += Aux*DF; // (C'(t),C(t) - P)
F2 += DF*DF + Aux * SecndDerOfC(i); // ((C'(t),C(t) - P))'
}
if(Abs(F2) < Eps)
EnCour = Standard_False;
else {
// Formula of Newton x(k+1) = x(k) - F(x(k))/F'(x(k))
TNew -= F1 / F2;
if(TNew < 0.) TNew = 0.;
if(TNew > 1.) TNew = 1.;
// Analysis of result
C->D0(TNew, ValOfC);
Dist = 0;
for(i = 1; i <= myDimension; i++) {
Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
Dist += Aux * Aux;
}
Dist = Sqrt(Dist);
Ecart = Dist0 - Dist;
if(Ecart <= -Seuil) {
// Pas d'amelioration on s'arrete
EnCour = Standard_False;
TNew = T0;
Dist = Dist0;
}
else if(Ecart <= Seuil)
// Convergence
NItCv++;
else
NItCv = 0;
if((NItCv >= 2) || (Iter >= NbIterations)) EnCour = Standard_False;
}
}
ProjTi(Ipnt) = TNew;
Distance(d0 + Ipnt) = Dist;
if(Dist > MaxErr) {
MaxErr = Dist;
NumPoints = Ipnt;
}
QuaErr += Dist * Dist;
AveErr += Dist;
}
NumPoints = NumPoints + myFirstPoint - 1;// Setting NumPoints to interval [myFirstPoint, myLastPoint]
}
|
