0023024: Update headers of OCCT files
[occt.git] / src / AppParCurves / AppParCurves_Function.gxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
3 //
4 // The content of this file is subject to the Open CASCADE Technology Public
5 // License Version 6.5 (the "License"). You may not use the content of this file
6 // except in compliance with the License. Please obtain a copy of the License
7 // at http://www.opencascade.org and read it completely before using this file.
8 //
9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11 //
12 // The Original Code and all software distributed under the License is
13 // distributed on an "AS IS" basis, without warranty of any kind, and the
14 // Initial Developer hereby disclaims all such warranties, including without
15 // limitation, any warranties of merchantability, fitness for a particular
16 // purpose or non-infringement. Please see the License for the specific terms
17 // and conditions governing the rights and limitations under the License.
18
19 // Lpa, le 20/09/91
20
21
22 // Calcul de la valeur de F et grad_F, connaissant le parametrage.
23 // Cette fonction, appelee par le gradient conjugue, calcul F et 
24 // DF(ui, Poles(ui)) ce qui implique un calcul des nouveaux poles 
25 //  a chaque appel.
26
27 #define No_Standard_RangeError
28 #define No_Standard_OutOfRange
29
30
31
32 #include <AppParCurves_MultiCurve.hxx>
33 #include <AppParCurves_MultiPoint.hxx>
34 #include <TColStd_HArray1OfInteger.hxx>
35 #include <gp_Pnt.hxx>
36 #include <gp_Pnt2d.hxx>
37 #include <gp_Vec.hxx>
38 #include <gp_Vec2d.hxx>
39 #include <TColgp_Array1OfPnt.hxx>
40 #include <TColgp_Array1OfPnt2d.hxx>
41 #include <AppParCurves_ConstraintCouple.hxx>
42
43 AppParCurves_Function::
44   AppParCurves_Function(const MultiLine& SSP,
45          const Standard_Integer FirstPoint,
46          const Standard_Integer LastPoint,
47          const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
48          const math_Vector& Parameters,
49          const Standard_Integer Deg) :
50          MyMultiLine(SSP),
51          MyMultiCurve(Deg+1),          
52          myParameters(Parameters.Lower(), Parameters.Upper()),
53          ValGrad_F(FirstPoint, LastPoint),
54          MyF(FirstPoint, LastPoint, 
55              1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
56          PTLX(FirstPoint, LastPoint, 
57              1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
58          PTLY(FirstPoint, LastPoint, 
59              1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
60          PTLZ(FirstPoint, LastPoint, 
61              1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
62          A(FirstPoint, LastPoint, 1, Deg+1),
63          DA(FirstPoint, LastPoint, 1, Deg+1),
64          MyLeastSquare(SSP, FirstPoint, LastPoint, 
65                        FirstConstraint(TheConstraints, FirstPoint),
66                        LastConstraint(TheConstraints, LastPoint), Deg+1)
67 {
68   Standard_Integer i;
69   for (i=Parameters.Lower(); i<=Parameters.Upper();i++)
70     myParameters(i)=Parameters(i);
71   FirstP = FirstPoint;
72   LastP = LastPoint;
73   myConstraints = TheConstraints;
74   NbP = LastP-FirstP+1;
75   Adeb = FirstP;
76   Afin = LastP;
77   Degre = Deg;
78   Contraintes = Standard_False;
79   Standard_Integer myindex;
80   Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
81   AppParCurves_ConstraintCouple mycouple;
82   AppParCurves_Constraint Cons;
83
84   for (i = low; i <= upp; i++) {
85     mycouple = TheConstraints->Value(i);
86     Cons = mycouple.Constraint();
87     myindex = mycouple.Index();
88     if (myindex == FirstP) {
89       if (Cons >= 1) Adeb = Adeb+1;
90     }
91     else if (myindex == LastP) {
92       if (Cons >= 1) Afin = Afin-1;
93     }
94     else {
95       if (Cons >= 1) Contraintes = Standard_True;
96     }
97   }
98
99   Standard_Integer nb3d = ToolLine::NbP3d(SSP);
100   Standard_Integer nb2d = ToolLine::NbP2d(SSP);
101   Standard_Integer mynb3d= nb3d, mynb2d=nb2d;
102   if (nb3d == 0) mynb3d = 1;
103   if (nb2d == 0) mynb2d = 1;
104
105   NbCu = nb3d+nb2d;
106   tabdim = new TColStd_HArray1OfInteger(0, NbCu-1);
107
108   if (Contraintes) {
109     for (i = 1; i <= NbCu; i++) {
110       if (i <= nb3d) tabdim->SetValue(i-1, 3);
111       else tabdim->SetValue(i-1, 2);
112     }
113
114     TColgp_Array1OfPnt TabP(1, mynb3d);
115     TColgp_Array1OfPnt2d TabP2d(1, mynb2d);
116     
117     for ( i = FirstP; i <= LastP; i++) {
118       if (nb3d != 0 && nb2d != 0) ToolLine::Value(SSP, i, TabP, TabP2d);
119       else if (nb3d != 0)         ToolLine::Value(SSP, i, TabP);
120       else                        ToolLine::Value(SSP, i, TabP2d);
121       for (Standard_Integer j = 1; j <= NbCu; j++) {
122         if (tabdim->Value(j-1) == 3) {
123           TabP(j).Coord(PTLX(i, j), PTLY(i, j),PTLZ(i, j));
124         }
125         else {
126           TabP2d(j).Coord(PTLX(i, j), PTLY(i, j));
127         }
128       }
129     }
130   }
131 }
132
133
134 AppParCurves_Constraint AppParCurves_Function::FirstConstraint
135   (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
136    const Standard_Integer FirstPoint) const
137 {
138   Standard_Integer i, myindex;
139   Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
140   AppParCurves_ConstraintCouple mycouple;
141   AppParCurves_Constraint Cons = AppParCurves_NoConstraint;
142
143   for (i = low; i <= upp; i++) {
144     mycouple = TheConstraints->Value(i);
145     Cons = mycouple.Constraint();
146     myindex = mycouple.Index();
147     if (myindex == FirstPoint) {
148       break;
149     }
150   }
151   return Cons;
152 }
153
154
155 AppParCurves_Constraint AppParCurves_Function::LastConstraint
156   (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
157    const Standard_Integer LastPoint) const
158 {
159   Standard_Integer i, myindex;
160   Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
161   AppParCurves_ConstraintCouple mycouple;
162   AppParCurves_Constraint Cons = AppParCurves_NoConstraint;
163
164   for (i = low; i <= upp; i++) {
165     mycouple = TheConstraints->Value(i);
166     Cons = mycouple.Constraint();
167     myindex = mycouple.Index();
168     if (myindex == LastPoint) {
169       break;
170     }
171   }
172   return Cons;
173 }
174
175
176
177
178 Standard_Boolean AppParCurves_Function::Value (const math_Vector& X, 
179                                                Standard_Real& F) {
180
181   myParameters = X;
182
183   // Resolution moindres carres:
184   // ===========================
185   MyLeastSquare.Perform(myParameters);
186   if (!(MyLeastSquare.IsDone())) { 
187     Done = Standard_False;
188     return Standard_False;
189   }
190   if (!Contraintes) {
191     MyLeastSquare.Error(FVal, ERR3d, ERR2d);
192     F = FVal;
193   }
194
195   // Resolution avec contraintes:
196   // ============================
197   else { 
198     Standard_Integer Npol = Degre+1;
199 //    Standard_Boolean Ext = Standard_True;
200     Standard_Integer Ci, i, j, dimen;
201     Standard_Real AA, BB, CC, AIJ, FX, FY, FZ, Fi;
202     math_Vector PTCXCI(1, Npol), PTCYCI(1, Npol), PTCZCI(1, Npol);
203     ERR3d = ERR2d = 0.0;
204     
205     MyMultiCurve = MyLeastSquare.BezierValue();
206
207     A = MyLeastSquare.FunctionMatrix();
208     ResolCons Resol(MyMultiLine, MyMultiCurve, FirstP, LastP, myConstraints,
209                     A, MyLeastSquare.DerivativeFunctionMatrix());
210     if (!Resol.IsDone()) {
211       Done = Standard_False;
212       return Standard_False;
213     }
214
215     // Calcul de F = Sum||C(ui)-Ptli||2  sur toutes les courbes :
216     // ========================================================================
217     FVal = 0.0;
218     
219     for (Ci = 1; Ci <= NbCu; Ci++) {
220       dimen = tabdim->Value(Ci-1);
221       for (j = 1; j <= Npol; j++) {
222         if (dimen == 3){ 
223           MyMultiCurve.Value(j).Point(Ci).Coord(PTCXCI(j),PTCYCI(j),PTCZCI(j));
224         }
225         else{ 
226           MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCXCI(j), PTCYCI(j));
227         }
228       }
229       
230       // Calcul de F:
231       // ============
232       for (i = Adeb; i <= Afin; i++) {
233         AA = 0.0; BB = 0.0; CC = 0.0;
234         for (j = 1; j <= Npol; j++) {
235           AIJ = A(i, j);
236           AA += AIJ*PTCXCI(j);
237           BB += AIJ*PTCYCI(j);
238           if (dimen == 3) { 
239             CC += AIJ*PTCZCI(j);
240           }
241         }
242         FX = AA-PTLX(i, Ci);
243         FY = BB-PTLY(i, Ci);
244         MyF(i,Ci) = FX*FX + FY*FY;
245         if (dimen == 3) {
246           FZ = CC-PTLZ(i,Ci);
247           MyF(i, Ci) += FZ*FZ;
248           Fi = MyF(i, Ci);
249           if (Sqrt(Fi) > ERR3d) ERR3d = Sqrt(Fi);
250         }
251         else {
252           Fi = MyF(i, Ci);
253           if (Sqrt(Fi) > ERR2d) ERR2d = Sqrt(Fi);
254         }
255         FVal += Fi;
256       }
257     }
258     F = FVal;
259   }  
260   return Standard_True;
261 }
262
263
264
265
266 void AppParCurves_Function::Perform(const math_Vector& X) {
267   Standard_Integer j;
268
269   myParameters = X;
270   // Resolution moindres carres:
271   // ===========================
272   MyLeastSquare.Perform(myParameters);
273
274   if (!(MyLeastSquare.IsDone())) { 
275     Done = Standard_False;
276     return;
277   }
278
279   for(j = myParameters.Lower(); j <= myParameters.Upper(); j++) {
280     ValGrad_F(j) = 0.0;
281   }
282
283   if (!Contraintes) {
284     MyLeastSquare.ErrorGradient(ValGrad_F, FVal, ERR3d, ERR2d);
285   }
286   else {
287     Standard_Integer Pi, Ci, i, k, dimen;
288     Standard_Integer  Npol = Degre+1;
289     Standard_Real Scal, AA, BB, CC, DAA, DBB, DCC;
290     Standard_Real FX, FY, FZ, AIJ, DAIJ, px, py, pz, Fi;
291     AppParCurves_Constraint Cons=AppParCurves_NoConstraint;
292     math_Matrix Grad_F(FirstP, LastP, 1, NbCu, 0.0);
293     math_Vector PTCXCI(1, Npol), PTCYCI(1, Npol), PTCZCI(1, Npol);
294     math_Vector PTCOXCI(1, Npol), PTCOYCI(1, Npol), PTCOZCI(1, Npol);
295 //    Standard_Boolean Ext = Standard_True;
296     ERR3d = ERR2d = 0.0;
297
298     math_Matrix PTCOX(1, Npol, 1, NbCu), PTCOY(1, Npol, 1, NbCu), 
299                 PTCOZ(1, Npol,1, NbCu);
300     math_Matrix PTCX(1, Npol, 1, NbCu), PTCY(1, Npol, 1, NbCu), 
301                 PTCZ(1, Npol,1, NbCu);
302     Standard_Integer Inc;
303
304     MyMultiCurve = MyLeastSquare.BezierValue();
305
306     for (Ci =1; Ci <= NbCu; Ci++) {
307       dimen = tabdim->Value(Ci-1);
308       for (j = 1; j <= Npol; j++) {
309         if (dimen == 3){ 
310           MyMultiCurve.Value(j).Point(Ci).Coord(PTCOX(j, Ci),
311                                                 PTCOY(j, Ci),
312                                                 PTCOZ(j, Ci));
313         }
314         else{ 
315           MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCOX(j, Ci), PTCOY(j, Ci));
316           PTCOZ(j, Ci) = 0.0;
317         }
318       }
319     }
320
321     A = MyLeastSquare.FunctionMatrix();
322     DA = MyLeastSquare.DerivativeFunctionMatrix();
323     
324     // Resolution avec contraintes:
325     // ============================
326     
327     ResolCons Resol(MyMultiLine, MyMultiCurve, FirstP, LastP, 
328                     myConstraints, A, DA);
329     if (!Resol.IsDone()) {
330       Done = Standard_False;
331       return;
332     }
333     
334     
335     // Calcul de F = Sum||C(ui)-Ptli||2 et du gradient non contraint de F pour
336     // chaque point PointIndex.
337     // ========================================================================
338     FVal = 0.0;
339     for(j = FirstP; j <= LastP; j++) {
340       ValGrad_F(j) = 0.0;
341     }
342
343     math_Matrix TrA(A.LowerCol(), A.UpperCol(), A.LowerRow(), A.UpperRow());
344     math_Matrix TrDA(DA.LowerCol(), DA.UpperCol(), DA.LowerRow(), DA.UpperRow());
345     math_Matrix RESTM(A.LowerCol(), A.UpperCol(), A.LowerCol(), A.UpperCol());
346
347     const math_Matrix& K = Resol.ConstraintMatrix();
348     const math_Matrix& DK = Resol.ConstraintDerivative(MyMultiLine, X, Degre, DA);
349     math_Matrix TK(K.LowerCol(), K.UpperCol(), K.LowerRow(), K.UpperRow());
350     TK = K.Transposed();
351     const math_Vector& Vardua = Resol.Duale();
352     math_Matrix KK(K.LowerCol(), K.UpperCol(), Vardua.Lower(), Vardua.Upper());
353     KK = (K.Transposed())*(Resol.InverseMatrix());
354     math_Matrix DTK(DK.LowerCol(), DK.UpperCol(), DK.LowerRow(), DK.UpperRow());
355     DTK = DK.Transposed();
356     TrA = A.Transposed();
357     TrDA = DA.Transposed();
358     RESTM = ((A.Transposed()*A).Inverse());
359
360     math_Vector DPTCO(1, K.ColNumber());
361     math_Matrix DPTCO1(FirstP, LastP, 1, K.ColNumber());
362     math_Vector DKPTC(1, K.RowNumber());
363
364
365
366
367     FVal = 0.0;
368     for (Ci = 1; Ci <= NbCu; Ci++) {
369       dimen = tabdim->Value(Ci-1);
370       for (j = 1; j <= Npol; j++) {
371         if (dimen == 3){ 
372           MyMultiCurve.Value(j).Point(Ci).Coord(PTCX(j, Ci), 
373                                                 PTCY(j, Ci), 
374                                                 PTCZ(j, Ci));
375         }
376         else{ 
377           MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCX(j, Ci), PTCY(j,Ci));
378           PTCZ(j, Ci) = 0.0;
379         }
380       }
381     }
382
383     
384     // Calcul du gradient sans contraintes:
385     // ====================================
386
387     for (Ci = 1; Ci <= NbCu; Ci++) {
388       dimen = tabdim->Value(Ci-1);
389       for (i = Adeb; i <= Afin; i++) {
390         AA = 0.0; BB = 0.0; CC = 0.0; DAA = 0.0; DBB = 0.0; DCC = 0.0;
391         for (j = 1; j <= Npol; j++) {
392           AIJ = A(i, j); DAIJ = DA(i, j);
393           px = PTCX(j, Ci); py = PTCY(j, Ci);
394           AA += AIJ*px;  BB += AIJ*py;
395           DAA += DAIJ*px;  DBB += DAIJ*py;
396           if (dimen == 3) { 
397             pz = PTCZ(j, Ci);
398             CC += AIJ*pz;  DCC += DAIJ*pz;
399           }
400         }
401         FX = AA-PTLX(i, Ci);
402         FY = BB-PTLY(i, Ci);
403         MyF(i,Ci) = FX*FX + FY*FY;
404         Grad_F(i, Ci) = 2.0*(DAA*FX + DBB*FY);
405         if (dimen == 3) {
406           FZ = CC-PTLZ(i,Ci);
407           MyF(i, Ci) += FZ*FZ;
408           Grad_F(i, Ci) += 2.0*DCC*FZ;
409           Fi = MyF(i, Ci);
410           if (Sqrt(Fi) > ERR3d) ERR3d = Sqrt(Fi);
411         }
412         else {
413           Fi = MyF(i, Ci);
414           if (Sqrt(Fi) > ERR2d) ERR2d = Sqrt(Fi);
415         }
416         FVal += Fi;
417         ValGrad_F(i) += Grad_F(i, Ci);
418       }
419     }
420
421
422     // Calcul de DK*PTC:
423     // =================
424     for (i = 1; i <= K.RowNumber(); i++) {
425       Inc = 0;
426       for (Ci = 1; Ci <= NbCu; Ci++) {
427         dimen = tabdim->Value(Ci-1);
428         DKPTC(i) = 0.0;
429         for (j = 1; j <= Npol; j++) {
430           DKPTC(i) += DK(i, j+Inc)*PTCX(j, Ci)+ DK(i, j+Inc+Npol)*PTCY(j, Ci);
431           if (dimen == 3) {
432             DKPTC(i) += DK(i, j+Inc+2*Npol)*PTCZ(j, Ci);
433           }
434         }
435         if (dimen == 3) Inc = Inc +3*Npol;
436         else Inc = Inc +2*Npol;
437       }
438     }
439     
440     math_Vector DERR(DTK.LowerRow(), DTK.UpperRow());
441     DERR = (DTK)*Vardua-KK* ((DKPTC) + K*(DTK)*Vardua);
442
443     // rajout du gradient avec contraintes:
444     // ====================================
445     // dPTCO1/duk = [d(TA)/duk*[A*PTCO-PTL] + TA*dA/duk*PTCO]
446
447
448     Inc = 0;
449
450     math_Vector Errx(A.LowerRow(), A.UpperRow());
451     math_Vector Erry(A.LowerRow(), A.UpperRow());
452     math_Vector Errz(A.LowerRow(), A.UpperRow());
453     math_Vector Scalx(DA.LowerRow(), DA.UpperRow());
454     math_Vector Scaly(DA.LowerRow(), DA.UpperRow());
455     math_Vector Scalz(DA.LowerRow(), DA.UpperRow());
456     math_Vector Erruzax(PTCXCI.Lower(), PTCXCI.Upper());
457     math_Vector Erruzay(PTCYCI.Lower(), PTCYCI.Upper());
458     math_Vector Erruzaz(PTCZCI.Lower(), PTCZCI.Upper());
459     math_Vector TrDAPI(TrDA.LowerRow(), TrDA.UpperRow());
460     math_Vector TrAPI(TrA.LowerRow(), TrA.UpperRow());
461
462     for (Ci = 1; Ci <= NbCu; Ci++) {
463       dimen = tabdim->Value(Ci-1);
464       PTCOXCI = PTCOX.Col(Ci);
465       PTCOYCI = PTCOY.Col(Ci);
466       PTCOZCI = PTCOZ.Col(Ci);
467       PTCXCI = PTCX.Col(Ci);
468       PTCYCI = PTCY.Col(Ci);
469       PTCZCI = PTCZ.Col(Ci);
470
471
472       Errx = (A*PTCOXCI - PTLX.Col(Ci));
473       Erry = (A*PTCOYCI - PTLY.Col(Ci));
474       Errz = (A*PTCOZCI - PTLZ.Col(Ci));
475       Scalx = (DA*PTCOXCI);   // Scal = DA * PTCO
476       Scaly = (DA*PTCOYCI);
477       Scalz = (DA*PTCOZCI);
478       Erruzax = (PTCXCI - PTCOXCI);
479       Erruzay = (PTCYCI - PTCOYCI);
480       Erruzaz = (PTCZCI - PTCOZCI);
481       
482       for (Pi = FirstP; Pi <= LastP; Pi++) {
483         TrDAPI = (TrDA.Col(Pi));
484         TrAPI = (TrA.Col(Pi));
485         Standard_Real Taa = TrAPI*A.Row(Pi);
486         Scal = 0.0;
487         for (j = 1; j <= Npol; j++) {
488           DPTCO1(Pi, j + Inc) = (TrDAPI*Errx(Pi)+TrAPI*Scalx(Pi))(j);
489           DPTCO1(Pi, j + Inc+ Npol) = (TrDAPI*Erry(Pi)+TrAPI*Scaly(Pi))(j);
490           Scal += DPTCO1(Pi, j+Inc)* Taa*Erruzax(j) + DPTCO1(Pi, j+Inc+Npol)*Taa*Erruzay(j);
491           if (dimen == 3) {
492             DPTCO1(Pi, j + Inc+ 2*Npol) = (TrDAPI*Errz(Pi)+TrAPI*Scalz(Pi))(j);
493             Scal += DPTCO1(Pi, j+Inc+2*Npol)*Taa*Erruzaz(j);
494           }
495         }
496         ValGrad_F(Pi) = ValGrad_F(Pi) - 2*Scal;
497       }
498       if (dimen == 3) Inc = Inc + 3*Npol;
499       else Inc = Inc +2*Npol;
500     }
501
502
503     // on calcule DPTCO = - RESTM * DPTCO1:
504     
505     // Calcul de DPTCO/duk:
506     // dPTCO/duk = -Inv(T(A)*A)*[d(TA)/duk*[A*PTCO-PTL] + TA*dA/duk*PTCO]
507
508     Standard_Integer low=myConstraints->Lower(), upp=myConstraints->Upper();
509     Inc = 0;
510     for (Pi = FirstP; Pi <= LastP; Pi++) {
511       for (i = low; i <= upp; i++) {
512         if (myConstraints->Value(i).Index() == Pi) {
513           Cons = myConstraints->Value(i).Constraint();
514           break;
515         }
516       }
517       if (Cons >= 1) {
518         Inc = 0;
519         for (Ci = 1; Ci <= NbCu; Ci++) {
520           dimen = tabdim->Value(Ci-1);
521           for (j = 1; j <= Npol; j++) {
522             DPTCO(j+Inc) = 0.0;
523             DPTCO(j+Inc+Npol) = 0.0;
524             if (dimen == 3) DPTCO(j+Inc+2*Npol) = 0.0;
525             for (k = 1; k <= Npol; k++) {
526               DPTCO(j+Inc) = DPTCO(j+Inc) -RESTM(j, k) * DPTCO1(Pi, j+Inc);
527               DPTCO(j+Inc+Npol)=DPTCO(j+Inc+Npol)-RESTM(j, k)*DPTCO1(Pi,j+Inc+Npol);
528               if (dimen == 3) {
529                 DPTCO(j+Inc+2*Npol) = DPTCO(j+Inc+2*Npol) 
530                   -RESTM(j, k) * DPTCO1(Pi, j+Inc+2*Npol);
531               }
532             }
533           }
534           if (dimen == 3) Inc += 3*Npol;
535           else Inc += 2*Npol;
536         }
537         
538         DERR = DERR-KK*K*DPTCO;
539         
540         Inc = 0;
541         for (Ci = 1; Ci <= NbCu; Ci++) {
542           dimen = tabdim->Value(Ci-1);
543           PTCOXCI = PTCOX.Col(Ci);
544           PTCOYCI = PTCOY.Col(Ci);
545           PTCOZCI = PTCOZ.Col(Ci);
546           PTCXCI = PTCX.Col(Ci);
547           PTCYCI = PTCY.Col(Ci);
548           PTCZCI = PTCZ.Col(Ci);
549           Erruzax = (PTCXCI - PTCOXCI);
550           Erruzay = (PTCYCI - PTCOYCI);
551           Erruzaz = (PTCZCI - PTCOZCI);
552           Scal = 0.0;
553           
554           for (j = 1; j <= Npol ; j++) {
555             Scal = (A(Pi, j)*Erruzax(j)) * (A(Pi, j)*DERR(j+Inc)) + 
556               (A(Pi, j)*Erruzay(j)) * (A(Pi, j)*DERR(j+Inc+Npol));
557             if (dimen == 3) {
558               Scal += (A(Pi, j)*Erruzax(j)) * (A(Pi, j)*DERR(j+Inc+2*Npol)); 
559             }
560           }
561           
562           ValGrad_F(Pi) = ValGrad_F(Pi) + 2*Scal;
563           if (dimen == 3) Inc = Inc +3*Npol;
564           else Inc = Inc + 2*Npol;
565         }
566       }
567     }
568     
569   }
570 }
571
572
573 Standard_Integer AppParCurves_Function::NbVariables() const{ 
574   return NbP;
575 }
576
577
578 Standard_Boolean AppParCurves_Function::Gradient (const math_Vector& X,
579                                                   math_Vector& G) {
580
581   Perform(X);
582   G = ValGrad_F;
583
584   return Standard_True;
585 }
586
587
588 Standard_Boolean AppParCurves_Function::Values (const math_Vector& X, 
589                                                 Standard_Real& F, 
590                                                 math_Vector& G) {
591
592
593   Perform(X);
594   F = FVal;
595   G = ValGrad_F;
596   return Standard_True;
597 }
598
599
600 const AppParCurves_MultiCurve& AppParCurves_Function::CurveValue() {
601   if (!Contraintes)  MyMultiCurve = MyLeastSquare.BezierValue();
602   return MyMultiCurve;
603 }
604
605
606 Standard_Real AppParCurves_Function::Error(const Standard_Integer IPoint,
607                                      const Standard_Integer CurveIndex) const {
608   return Sqrt(MyF(IPoint, CurveIndex));
609 }
610
611 Standard_Real AppParCurves_Function::MaxError3d() const
612 {
613   return ERR3d;
614 }
615
616 Standard_Real AppParCurves_Function::MaxError2d() const
617 {
618   return ERR2d;
619 }
620
621
622
623 const math_Vector& AppParCurves_Function::NewParameters() const
624 {
625   return myParameters;
626 }