d69868fa6b96ae2356400043a20350f455acb091
[occt.git] / src / AppParCurves / AppParCurves_Variational_3.gxx
1 // Created on: 1998-12-08
2 // Created by: Igor FEOKTISTOV
3 // Copyright (c) 1998-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
5 //
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
10 //
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
13 //
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
20
21
22
23
24 void AppParCurves_Variational::Project(const Handle(FEmTool_Curve)& C,
25                                        const TColStd_Array1OfReal& Ti,
26                                        TColStd_Array1OfReal& ProjTi,
27                                        TColStd_Array1OfReal& Distance,
28                                        Standard_Integer& NumPoints,
29                                        Standard_Real& MaxErr,
30                                        Standard_Real& QuaErr,
31                                        Standard_Real& AveErr,
32                                        const Standard_Integer NbIterations) const
33
34 {
35   // Initialisation
36
37   const Standard_Real Seuil = 1.e-9, Eps = 1.e-12;
38
39   MaxErr = QuaErr = AveErr = 0.;
40
41   Standard_Integer Ipnt, NItCv, Iter, i, i0 = -myDimension, d0 = Distance.Lower() - 1;
42
43   Standard_Real TNew, Dist, T0, Dist0, F1, F2, Aux, DF, Ecart;
44
45   Standard_Boolean EnCour;
46
47   TColStd_Array1OfReal ValOfC(1, myDimension), FirstDerOfC(1, myDimension), 
48                        SecndDerOfC(1, myDimension);
49
50   for(Ipnt = 1; Ipnt <= ProjTi.Length(); Ipnt++) {
51
52     i0 += myDimension;
53
54     TNew = Ti(Ipnt);
55
56     EnCour = Standard_True;
57     NItCv = 0;
58     Iter = 0;
59     C->D0(TNew, ValOfC);
60
61     Dist = 0;
62     for(i = 1; i <= myDimension; i++) {
63       Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
64       Dist += Aux * Aux;
65     }
66     Dist = Sqrt(Dist);
67          
68     // ------- Newton's method for solving (C'(t),C(t) - P) = 0
69
70     while( EnCour ) {
71     
72       Iter++;
73       T0 = TNew;
74       Dist0 = Dist;
75
76       C->D2(TNew, SecndDerOfC);
77       C->D1(TNew, FirstDerOfC);
78
79       F1 = F2 = 0.;
80       for(i = 1; i <= myDimension; i++) {
81         Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
82         DF = FirstDerOfC(i);
83         F1 += Aux*DF;          // (C'(t),C(t) - P)
84         F2 += DF*DF + Aux * SecndDerOfC(i); // ((C'(t),C(t) - P))'
85       }
86
87       if(Abs(F2) < Eps) 
88         EnCour = Standard_False;
89       else {
90         // Formula of Newton x(k+1) = x(k) - F(x(k))/F'(x(k))
91         TNew -= F1 / F2;
92         if(TNew < 0.) TNew = 0.;
93         if(TNew > 1.) TNew = 1.;
94       
95
96         // Analysis of result
97
98         C->D0(TNew, ValOfC);
99
100         Dist = 0;
101         for(i = 1; i <= myDimension; i++) {
102           Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
103           Dist += Aux * Aux;
104         }
105         Dist = Sqrt(Dist);
106
107         Ecart = Dist0 - Dist;
108
109         if(Ecart <= -Seuil) {
110           // Pas d'amelioration on s'arrete
111           EnCour = Standard_False;
112           TNew = T0;
113           Dist = Dist0;
114         }
115         else if(Ecart <= Seuil) 
116           // Convergence
117           NItCv++;
118         else
119           NItCv = 0;
120
121         if((NItCv >= 2) || (Iter >= NbIterations)) EnCour = Standard_False;
122             
123       }
124     }
125
126     
127     ProjTi(Ipnt) = TNew;
128     Distance(d0 + Ipnt) = Dist;
129     if(Dist > MaxErr) {
130       MaxErr = Dist;
131       NumPoints = Ipnt;
132     }
133     QuaErr += Dist * Dist;
134     AveErr += Dist;
135   }
136
137   NumPoints = NumPoints + myFirstPoint - 1;// Setting NumPoints to interval [myFirstPoint, myLastPoint]
138
139 }