src/AppParCurves/AppParCurves_Variational_3.gxx

   1 // File:        AppParCurves_Variational_3.gxx
   2 // Created:     Tue Dec  8 09:08:54 1998
   3 // Author:      Igor FEOKTISTOV
   4 //              <ifv@paradox.nnov.matra-dtv.fr>
   5
   6
   7
   8 void AppParCurves_Variational::Project(const Handle(FEmTool_Curve)& C,
   9                                        const TColStd_Array1OfReal& Ti,
  10                                        TColStd_Array1OfReal& ProjTi,
  11                                        TColStd_Array1OfReal& Distance,
  12                                        Standard_Integer& NumPoints,
  13                                        Standard_Real& MaxErr,
  14                                        Standard_Real& QuaErr,
  15                                        Standard_Real& AveErr,
  16                                        const Standard_Integer NbIterations) const
  17
  18 {
  19   // Initialisation
  20
  21   const Standard_Real Seuil = 1.e-9, Eps = 1.e-12;
  22
  23   MaxErr = QuaErr = AveErr = 0.;
  24
  25   Standard_Integer Ipnt, NItCv, Iter, i, i0 = -myDimension, d0 = Distance.Lower() - 1;
  26
  27   Standard_Real TNew, Dist, T0, Dist0, F1, F2, Aux, DF, Ecart;
  28
  29   Standard_Boolean EnCour;
  30
  31   TColStd_Array1OfReal ValOfC(1, myDimension), FirstDerOfC(1, myDimension),
  32                        SecndDerOfC(1, myDimension);
  33
  34   for(Ipnt = 1; Ipnt <= ProjTi.Length(); Ipnt++) {
  35
  36     i0 += myDimension;
  37
  38     TNew = Ti(Ipnt);
  39
  40     EnCour = Standard_True;
  41     NItCv = 0;
  42     Iter = 0;
  43     C->D0(TNew, ValOfC);
  44
  45     Dist = 0;
  46     for(i = 1; i <= myDimension; i++) {
  47       Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
  48       Dist += Aux * Aux;
  49     }
  50     Dist = Sqrt(Dist);
  51
  52     // ------- Newton's method for solving (C'(t),C(t) - P) = 0
  53
  54     while( EnCour ) {
  55
  56       Iter++;
  57       T0 = TNew;
  58       Dist0 = Dist;
  59
  60       C->D2(TNew, SecndDerOfC);
  61       C->D1(TNew, FirstDerOfC);
  62
  63       F1 = F2 = 0.;
  64       for(i = 1; i <= myDimension; i++) {
  65         Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
  66         DF = FirstDerOfC(i);
  67         F1 += Aux*DF;          // (C'(t),C(t) - P)
  68         F2 += DF*DF + Aux * SecndDerOfC(i); // ((C'(t),C(t) - P))'
  69       }
  70
  71       if(Abs(F2) < Eps)
  72         EnCour = Standard_False;
  73       else {
  74         // Formula of Newton x(k+1) = x(k) - F(x(k))/F'(x(k))
  75         TNew -= F1 / F2;
  76         if(TNew < 0.) TNew = 0.;
  77         if(TNew > 1.) TNew = 1.;
  78
  79
  80         // Analysis of result
  81
  82         C->D0(TNew, ValOfC);
  83
  84         Dist = 0;
  85         for(i = 1; i <= myDimension; i++) {
  86           Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
  87           Dist += Aux * Aux;
  88         }
  89         Dist = Sqrt(Dist);
  90
  91         Ecart = Dist0 - Dist;
  92
  93         if(Ecart <= -Seuil) {
  94           // Pas d'amelioration on s'arrete
  95           EnCour = Standard_False;
  96           TNew = T0;
  97           Dist = Dist0;
  98         }
  99         else if(Ecart <= Seuil)
 100           // Convergence
 101           NItCv++;
 102         else
 103           NItCv = 0;
 104
 105         if((NItCv >= 2) || (Iter >= NbIterations)) EnCour = Standard_False;
 106
 107       }
 108     }
 109
 110
 111     ProjTi(Ipnt) = TNew;
 112     Distance(d0 + Ipnt) = Dist;
 113     if(Dist > MaxErr) {
 114       MaxErr = Dist;
 115       NumPoints = Ipnt;
 116     }
 117     QuaErr += Dist * Dist;
 118     AveErr += Dist;
 119   }
 120
 121   NumPoints = NumPoints + myFirstPoint - 1;// Setting NumPoints to interval [myFirstPoint, myLastPoint]
 122
 123 }