Warning, /include/opencascade/AppParCurves_Function.gxx is written in an unsupported language. File is not indexed.
0001 // Copyright (c) 1995-1999 Matra Datavision
0002 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0003 //
0004 // This file is part of Open CASCADE Technology software library.
0005 //
0006 // This library is free software; you can redistribute it and/or modify it under
0007 // the terms of the GNU Lesser General Public License version 2.1 as published
0008 // by the Free Software Foundation, with special exception defined in the file
0009 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0010 // distribution for complete text of the license and disclaimer of any warranty.
0011 //
0012 // Alternatively, this file may be used under the terms of Open CASCADE
0013 // commercial license or contractual agreement.
0014
0015 // Lpa, le 20/09/91
0016
0017
0018 // Calcul de la valeur de F et grad_F, connaissant le parametrage.
0019 // Cette fonction, appelee par le gradient conjugue, calcul F et
0020 // DF(ui, Poles(ui)) ce qui implique un calcul des nouveaux poles
0021 // a chaque appel.
0022
0023 #define No_Standard_RangeError
0024 #define No_Standard_OutOfRange
0025
0026
0027
0028 #include <AppParCurves_MultiCurve.hxx>
0029 #include <AppParCurves_MultiPoint.hxx>
0030 #include <TColStd_HArray1OfInteger.hxx>
0031 #include <gp_Pnt.hxx>
0032 #include <gp_Pnt2d.hxx>
0033 #include <gp_Vec.hxx>
0034 #include <gp_Vec2d.hxx>
0035 #include <TColgp_Array1OfPnt.hxx>
0036 #include <TColgp_Array1OfPnt2d.hxx>
0037 #include <AppParCurves_ConstraintCouple.hxx>
0038
0039 AppParCurves_Function::
0040 AppParCurves_Function(const MultiLine& SSP,
0041 const Standard_Integer FirstPoint,
0042 const Standard_Integer LastPoint,
0043 const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0044 const math_Vector& Parameters,
0045 const Standard_Integer Deg) :
0046 MyMultiLine(SSP),
0047 MyMultiCurve(Deg+1),
0048 myParameters(Parameters.Lower(), Parameters.Upper()),
0049 ValGrad_F(FirstPoint, LastPoint),
0050 MyF(FirstPoint, LastPoint,
0051 1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
0052 PTLX(FirstPoint, LastPoint,
0053 1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
0054 PTLY(FirstPoint, LastPoint,
0055 1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
0056 PTLZ(FirstPoint, LastPoint,
0057 1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
0058 A(FirstPoint, LastPoint, 1, Deg+1),
0059 DA(FirstPoint, LastPoint, 1, Deg+1),
0060 MyLeastSquare(SSP, FirstPoint, LastPoint,
0061 FirstConstraint(TheConstraints, FirstPoint),
0062 LastConstraint(TheConstraints, LastPoint), Deg+1)
0063 {
0064 Standard_Integer i;
0065 for (i=Parameters.Lower(); i<=Parameters.Upper();i++)
0066 myParameters(i)=Parameters(i);
0067 FirstP = FirstPoint;
0068 LastP = LastPoint;
0069 myConstraints = TheConstraints;
0070 NbP = LastP-FirstP+1;
0071 Adeb = FirstP;
0072 Afin = LastP;
0073 Degre = Deg;
0074 Contraintes = Standard_False;
0075 Standard_Integer myindex;
0076 Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
0077 AppParCurves_ConstraintCouple mycouple;
0078 AppParCurves_Constraint Cons;
0079
0080 for (i = low; i <= upp; i++) {
0081 mycouple = TheConstraints->Value(i);
0082 Cons = mycouple.Constraint();
0083 myindex = mycouple.Index();
0084 if (myindex == FirstP) {
0085 if (Cons >= 1) Adeb = Adeb+1;
0086 }
0087 else if (myindex == LastP) {
0088 if (Cons >= 1) Afin = Afin-1;
0089 }
0090 else {
0091 if (Cons >= 1) Contraintes = Standard_True;
0092 }
0093 }
0094
0095 Standard_Integer nb3d = ToolLine::NbP3d(SSP);
0096 Standard_Integer nb2d = ToolLine::NbP2d(SSP);
0097 Standard_Integer mynb3d= nb3d, mynb2d=nb2d;
0098 if (nb3d == 0) mynb3d = 1;
0099 if (nb2d == 0) mynb2d = 1;
0100
0101 NbCu = nb3d+nb2d;
0102 tabdim = new TColStd_HArray1OfInteger(0, NbCu-1);
0103
0104 if (Contraintes) {
0105 for (i = 1; i <= NbCu; i++) {
0106 if (i <= nb3d) tabdim->SetValue(i-1, 3);
0107 else tabdim->SetValue(i-1, 2);
0108 }
0109
0110 TColgp_Array1OfPnt TabP(1, mynb3d);
0111 TColgp_Array1OfPnt2d TabP2d(1, mynb2d);
0112
0113 for ( i = FirstP; i <= LastP; i++) {
0114 if (nb3d != 0 && nb2d != 0) ToolLine::Value(SSP, i, TabP, TabP2d);
0115 else if (nb3d != 0) ToolLine::Value(SSP, i, TabP);
0116 else ToolLine::Value(SSP, i, TabP2d);
0117 for (Standard_Integer j = 1; j <= NbCu; j++) {
0118 if (tabdim->Value(j-1) == 3) {
0119 TabP(j).Coord(PTLX(i, j), PTLY(i, j),PTLZ(i, j));
0120 }
0121 else {
0122 TabP2d(j).Coord(PTLX(i, j), PTLY(i, j));
0123 }
0124 }
0125 }
0126 }
0127 }
0128
0129
0130 AppParCurves_Constraint AppParCurves_Function::FirstConstraint
0131 (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0132 const Standard_Integer FirstPoint) const
0133 {
0134 Standard_Integer i, myindex;
0135 Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
0136 AppParCurves_ConstraintCouple mycouple;
0137 AppParCurves_Constraint Cons = AppParCurves_NoConstraint;
0138
0139 for (i = low; i <= upp; i++) {
0140 mycouple = TheConstraints->Value(i);
0141 Cons = mycouple.Constraint();
0142 myindex = mycouple.Index();
0143 if (myindex == FirstPoint) {
0144 break;
0145 }
0146 }
0147 return Cons;
0148 }
0149
0150
0151 AppParCurves_Constraint AppParCurves_Function::LastConstraint
0152 (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0153 const Standard_Integer LastPoint) const
0154 {
0155 Standard_Integer i, myindex;
0156 Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
0157 AppParCurves_ConstraintCouple mycouple;
0158 AppParCurves_Constraint Cons = AppParCurves_NoConstraint;
0159
0160 for (i = low; i <= upp; i++) {
0161 mycouple = TheConstraints->Value(i);
0162 Cons = mycouple.Constraint();
0163 myindex = mycouple.Index();
0164 if (myindex == LastPoint) {
0165 break;
0166 }
0167 }
0168 return Cons;
0169 }
0170
0171
0172
0173
0174 Standard_Boolean AppParCurves_Function::Value (const math_Vector& X,
0175 Standard_Real& F) {
0176
0177 myParameters = X;
0178
0179 // Resolution moindres carres:
0180 // ===========================
0181 MyLeastSquare.Perform(myParameters);
0182 if (!(MyLeastSquare.IsDone())) {
0183 Done = Standard_False;
0184 return Standard_False;
0185 }
0186 if (!Contraintes) {
0187 MyLeastSquare.Error(FVal, ERR3d, ERR2d);
0188 F = FVal;
0189 }
0190
0191 // Resolution avec contraintes:
0192 // ============================
0193 else {
0194 Standard_Integer Npol = Degre+1;
0195 // Standard_Boolean Ext = Standard_True;
0196 Standard_Integer Ci, i, j, dimen;
0197 Standard_Real AA, BB, CC, AIJ, FX, FY, FZ, Fi;
0198 math_Vector PTCXCI(1, Npol), PTCYCI(1, Npol), PTCZCI(1, Npol);
0199 ERR3d = ERR2d = 0.0;
0200
0201 MyMultiCurve = MyLeastSquare.BezierValue();
0202
0203 A = MyLeastSquare.FunctionMatrix();
0204 ResolCons Resol(MyMultiLine, MyMultiCurve, FirstP, LastP, myConstraints,
0205 A, MyLeastSquare.DerivativeFunctionMatrix());
0206 if (!Resol.IsDone()) {
0207 Done = Standard_False;
0208 return Standard_False;
0209 }
0210
0211 // Calcul de F = Sum||C(ui)-Ptli||2 sur toutes les courbes :
0212 // ========================================================================
0213 FVal = 0.0;
0214
0215 for (Ci = 1; Ci <= NbCu; Ci++) {
0216 dimen = tabdim->Value(Ci-1);
0217 for (j = 1; j <= Npol; j++) {
0218 if (dimen == 3){
0219 MyMultiCurve.Value(j).Point(Ci).Coord(PTCXCI(j),PTCYCI(j),PTCZCI(j));
0220 }
0221 else{
0222 MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCXCI(j), PTCYCI(j));
0223 }
0224 }
0225
0226 // Calcul de F:
0227 // ============
0228 for (i = Adeb; i <= Afin; i++) {
0229 AA = 0.0; BB = 0.0; CC = 0.0;
0230 for (j = 1; j <= Npol; j++) {
0231 AIJ = A(i, j);
0232 AA += AIJ*PTCXCI(j);
0233 BB += AIJ*PTCYCI(j);
0234 if (dimen == 3) {
0235 CC += AIJ*PTCZCI(j);
0236 }
0237 }
0238 FX = AA-PTLX(i, Ci);
0239 FY = BB-PTLY(i, Ci);
0240 MyF(i,Ci) = FX*FX + FY*FY;
0241 if (dimen == 3) {
0242 FZ = CC-PTLZ(i,Ci);
0243 MyF(i, Ci) += FZ*FZ;
0244 Fi = MyF(i, Ci);
0245 if (Sqrt(Fi) > ERR3d) ERR3d = Sqrt(Fi);
0246 }
0247 else {
0248 Fi = MyF(i, Ci);
0249 if (Sqrt(Fi) > ERR2d) ERR2d = Sqrt(Fi);
0250 }
0251 FVal += Fi;
0252 }
0253 }
0254 F = FVal;
0255 }
0256 return Standard_True;
0257 }
0258
0259
0260
0261
0262 void AppParCurves_Function::Perform(const math_Vector& X) {
0263 Standard_Integer j;
0264
0265 myParameters = X;
0266 // Resolution moindres carres:
0267 // ===========================
0268 MyLeastSquare.Perform(myParameters);
0269
0270 if (!(MyLeastSquare.IsDone())) {
0271 Done = Standard_False;
0272 return;
0273 }
0274
0275 for(j = myParameters.Lower(); j <= myParameters.Upper(); j++) {
0276 ValGrad_F(j) = 0.0;
0277 }
0278
0279 if (!Contraintes) {
0280 MyLeastSquare.ErrorGradient(ValGrad_F, FVal, ERR3d, ERR2d);
0281 }
0282 else {
0283 Standard_Integer Pi, Ci, i, k, dimen;
0284 Standard_Integer Npol = Degre+1;
0285 Standard_Real Scal, AA, BB, CC, DAA, DBB, DCC;
0286 Standard_Real FX, FY, FZ, AIJ, DAIJ, px, py, pz, Fi;
0287 AppParCurves_Constraint Cons=AppParCurves_NoConstraint;
0288 math_Matrix Grad_F(FirstP, LastP, 1, NbCu, 0.0);
0289 math_Vector PTCXCI(1, Npol), PTCYCI(1, Npol), PTCZCI(1, Npol);
0290 math_Vector PTCOXCI(1, Npol), PTCOYCI(1, Npol), PTCOZCI(1, Npol);
0291 // Standard_Boolean Ext = Standard_True;
0292 ERR3d = ERR2d = 0.0;
0293
0294 math_Matrix PTCOX(1, Npol, 1, NbCu), PTCOY(1, Npol, 1, NbCu),
0295 PTCOZ(1, Npol,1, NbCu);
0296 math_Matrix PTCX(1, Npol, 1, NbCu), PTCY(1, Npol, 1, NbCu),
0297 PTCZ(1, Npol,1, NbCu);
0298 Standard_Integer Inc;
0299
0300 MyMultiCurve = MyLeastSquare.BezierValue();
0301
0302 for (Ci =1; Ci <= NbCu; Ci++) {
0303 dimen = tabdim->Value(Ci-1);
0304 for (j = 1; j <= Npol; j++) {
0305 if (dimen == 3){
0306 MyMultiCurve.Value(j).Point(Ci).Coord(PTCOX(j, Ci),
0307 PTCOY(j, Ci),
0308 PTCOZ(j, Ci));
0309 }
0310 else{
0311 MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCOX(j, Ci), PTCOY(j, Ci));
0312 PTCOZ(j, Ci) = 0.0;
0313 }
0314 }
0315 }
0316
0317 A = MyLeastSquare.FunctionMatrix();
0318 DA = MyLeastSquare.DerivativeFunctionMatrix();
0319
0320 // Resolution avec contraintes:
0321 // ============================
0322
0323 ResolCons Resol(MyMultiLine, MyMultiCurve, FirstP, LastP,
0324 myConstraints, A, DA);
0325 if (!Resol.IsDone()) {
0326 Done = Standard_False;
0327 return;
0328 }
0329
0330
0331 // Calcul de F = Sum||C(ui)-Ptli||2 et du gradient non contraint de F pour
0332 // chaque point PointIndex.
0333 // ========================================================================
0334 FVal = 0.0;
0335 for(j = FirstP; j <= LastP; j++) {
0336 ValGrad_F(j) = 0.0;
0337 }
0338
0339 math_Matrix TrA(A.LowerCol(), A.UpperCol(), A.LowerRow(), A.UpperRow());
0340 math_Matrix TrDA(DA.LowerCol(), DA.UpperCol(), DA.LowerRow(), DA.UpperRow());
0341 math_Matrix RESTM(A.LowerCol(), A.UpperCol(), A.LowerCol(), A.UpperCol());
0342
0343 const math_Matrix& K = Resol.ConstraintMatrix();
0344 const math_Matrix& DK = Resol.ConstraintDerivative(MyMultiLine, X, Degre, DA);
0345 math_Matrix TK(K.LowerCol(), K.UpperCol(), K.LowerRow(), K.UpperRow());
0346 TK = K.Transposed();
0347 const math_Vector& Vardua = Resol.Duale();
0348 math_Matrix KK(K.LowerCol(), K.UpperCol(), Vardua.Lower(), Vardua.Upper());
0349 KK = (K.Transposed())*(Resol.InverseMatrix());
0350 math_Matrix DTK(DK.LowerCol(), DK.UpperCol(), DK.LowerRow(), DK.UpperRow());
0351 DTK = DK.Transposed();
0352 TrA = A.Transposed();
0353 TrDA = DA.Transposed();
0354 RESTM = ((A.Transposed()*A).Inverse());
0355
0356 math_Vector DPTCO(1, K.ColNumber());
0357 math_Matrix DPTCO1(FirstP, LastP, 1, K.ColNumber());
0358 math_Vector DKPTC(1, K.RowNumber());
0359
0360
0361
0362
0363 FVal = 0.0;
0364 for (Ci = 1; Ci <= NbCu; Ci++) {
0365 dimen = tabdim->Value(Ci-1);
0366 for (j = 1; j <= Npol; j++) {
0367 if (dimen == 3){
0368 MyMultiCurve.Value(j).Point(Ci).Coord(PTCX(j, Ci),
0369 PTCY(j, Ci),
0370 PTCZ(j, Ci));
0371 }
0372 else{
0373 MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCX(j, Ci), PTCY(j,Ci));
0374 PTCZ(j, Ci) = 0.0;
0375 }
0376 }
0377 }
0378
0379
0380 // Calcul du gradient sans contraintes:
0381 // ====================================
0382
0383 for (Ci = 1; Ci <= NbCu; Ci++) {
0384 dimen = tabdim->Value(Ci-1);
0385 for (i = Adeb; i <= Afin; i++) {
0386 AA = 0.0; BB = 0.0; CC = 0.0; DAA = 0.0; DBB = 0.0; DCC = 0.0;
0387 for (j = 1; j <= Npol; j++) {
0388 AIJ = A(i, j); DAIJ = DA(i, j);
0389 px = PTCX(j, Ci); py = PTCY(j, Ci);
0390 AA += AIJ*px; BB += AIJ*py;
0391 DAA += DAIJ*px; DBB += DAIJ*py;
0392 if (dimen == 3) {
0393 pz = PTCZ(j, Ci);
0394 CC += AIJ*pz; DCC += DAIJ*pz;
0395 }
0396 }
0397 FX = AA-PTLX(i, Ci);
0398 FY = BB-PTLY(i, Ci);
0399 MyF(i,Ci) = FX*FX + FY*FY;
0400 Grad_F(i, Ci) = 2.0*(DAA*FX + DBB*FY);
0401 if (dimen == 3) {
0402 FZ = CC-PTLZ(i,Ci);
0403 MyF(i, Ci) += FZ*FZ;
0404 Grad_F(i, Ci) += 2.0*DCC*FZ;
0405 Fi = MyF(i, Ci);
0406 if (Sqrt(Fi) > ERR3d) ERR3d = Sqrt(Fi);
0407 }
0408 else {
0409 Fi = MyF(i, Ci);
0410 if (Sqrt(Fi) > ERR2d) ERR2d = Sqrt(Fi);
0411 }
0412 FVal += Fi;
0413 ValGrad_F(i) += Grad_F(i, Ci);
0414 }
0415 }
0416
0417
0418 // Calcul de DK*PTC:
0419 // =================
0420 for (i = 1; i <= K.RowNumber(); i++) {
0421 Inc = 0;
0422 for (Ci = 1; Ci <= NbCu; Ci++) {
0423 dimen = tabdim->Value(Ci-1);
0424 DKPTC(i) = 0.0;
0425 for (j = 1; j <= Npol; j++) {
0426 DKPTC(i) += DK(i, j+Inc)*PTCX(j, Ci)+ DK(i, j+Inc+Npol)*PTCY(j, Ci);
0427 if (dimen == 3) {
0428 DKPTC(i) += DK(i, j+Inc+2*Npol)*PTCZ(j, Ci);
0429 }
0430 }
0431 if (dimen == 3) Inc = Inc +3*Npol;
0432 else Inc = Inc +2*Npol;
0433 }
0434 }
0435
0436 math_Vector DERR(DTK.LowerRow(), DTK.UpperRow());
0437 DERR = (DTK)*Vardua-KK* ((DKPTC) + K*(DTK)*Vardua);
0438
0439 // rajout du gradient avec contraintes:
0440 // ====================================
0441 // dPTCO1/duk = [d(TA)/duk*[A*PTCO-PTL] + TA*dA/duk*PTCO]
0442
0443
0444 Inc = 0;
0445
0446 math_Vector Errx(A.LowerRow(), A.UpperRow());
0447 math_Vector Erry(A.LowerRow(), A.UpperRow());
0448 math_Vector Errz(A.LowerRow(), A.UpperRow());
0449 math_Vector Scalx(DA.LowerRow(), DA.UpperRow());
0450 math_Vector Scaly(DA.LowerRow(), DA.UpperRow());
0451 math_Vector Scalz(DA.LowerRow(), DA.UpperRow());
0452 math_Vector Erruzax(PTCXCI.Lower(), PTCXCI.Upper());
0453 math_Vector Erruzay(PTCYCI.Lower(), PTCYCI.Upper());
0454 math_Vector Erruzaz(PTCZCI.Lower(), PTCZCI.Upper());
0455 math_Vector TrDAPI(TrDA.LowerRow(), TrDA.UpperRow());
0456 math_Vector TrAPI(TrA.LowerRow(), TrA.UpperRow());
0457
0458 for (Ci = 1; Ci <= NbCu; Ci++) {
0459 dimen = tabdim->Value(Ci-1);
0460 PTCOXCI = PTCOX.Col(Ci);
0461 PTCOYCI = PTCOY.Col(Ci);
0462 PTCOZCI = PTCOZ.Col(Ci);
0463 PTCXCI = PTCX.Col(Ci);
0464 PTCYCI = PTCY.Col(Ci);
0465 PTCZCI = PTCZ.Col(Ci);
0466
0467
0468 Errx = (A*PTCOXCI - PTLX.Col(Ci));
0469 Erry = (A*PTCOYCI - PTLY.Col(Ci));
0470 Errz = (A*PTCOZCI - PTLZ.Col(Ci));
0471 Scalx = (DA*PTCOXCI); // Scal = DA * PTCO
0472 Scaly = (DA*PTCOYCI);
0473 Scalz = (DA*PTCOZCI);
0474 Erruzax = (PTCXCI - PTCOXCI);
0475 Erruzay = (PTCYCI - PTCOYCI);
0476 Erruzaz = (PTCZCI - PTCOZCI);
0477
0478 for (Pi = FirstP; Pi <= LastP; Pi++) {
0479 TrDAPI = (TrDA.Col(Pi));
0480 TrAPI = (TrA.Col(Pi));
0481 Standard_Real Taa = TrAPI*A.Row(Pi);
0482 Scal = 0.0;
0483 for (j = 1; j <= Npol; j++) {
0484 DPTCO1(Pi, j + Inc) = (TrDAPI*Errx(Pi)+TrAPI*Scalx(Pi))(j);
0485 DPTCO1(Pi, j + Inc+ Npol) = (TrDAPI*Erry(Pi)+TrAPI*Scaly(Pi))(j);
0486 Scal += DPTCO1(Pi, j+Inc)* Taa*Erruzax(j) + DPTCO1(Pi, j+Inc+Npol)*Taa*Erruzay(j);
0487 if (dimen == 3) {
0488 DPTCO1(Pi, j + Inc+ 2*Npol) = (TrDAPI*Errz(Pi)+TrAPI*Scalz(Pi))(j);
0489 Scal += DPTCO1(Pi, j+Inc+2*Npol)*Taa*Erruzaz(j);
0490 }
0491 }
0492 ValGrad_F(Pi) = ValGrad_F(Pi) - 2*Scal;
0493 }
0494 if (dimen == 3) Inc = Inc + 3*Npol;
0495 else Inc = Inc +2*Npol;
0496 }
0497
0498
0499 // on calcule DPTCO = - RESTM * DPTCO1:
0500
0501 // Calcul de DPTCO/duk:
0502 // dPTCO/duk = -Inv(T(A)*A)*[d(TA)/duk*[A*PTCO-PTL] + TA*dA/duk*PTCO]
0503
0504 Standard_Integer low=myConstraints->Lower(), upp=myConstraints->Upper();
0505 Inc = 0;
0506 for (Pi = FirstP; Pi <= LastP; Pi++) {
0507 for (i = low; i <= upp; i++) {
0508 if (myConstraints->Value(i).Index() == Pi) {
0509 Cons = myConstraints->Value(i).Constraint();
0510 break;
0511 }
0512 }
0513 if (Cons >= 1) {
0514 Inc = 0;
0515 for (Ci = 1; Ci <= NbCu; Ci++) {
0516 dimen = tabdim->Value(Ci-1);
0517 for (j = 1; j <= Npol; j++) {
0518 DPTCO(j+Inc) = 0.0;
0519 DPTCO(j+Inc+Npol) = 0.0;
0520 if (dimen == 3) DPTCO(j+Inc+2*Npol) = 0.0;
0521 for (k = 1; k <= Npol; k++) {
0522 DPTCO(j+Inc) = DPTCO(j+Inc) -RESTM(j, k) * DPTCO1(Pi, j+Inc);
0523 DPTCO(j+Inc+Npol)=DPTCO(j+Inc+Npol)-RESTM(j, k)*DPTCO1(Pi,j+Inc+Npol);
0524 if (dimen == 3) {
0525 DPTCO(j+Inc+2*Npol) = DPTCO(j+Inc+2*Npol)
0526 -RESTM(j, k) * DPTCO1(Pi, j+Inc+2*Npol);
0527 }
0528 }
0529 }
0530 if (dimen == 3) Inc += 3*Npol;
0531 else Inc += 2*Npol;
0532 }
0533
0534 DERR = DERR-KK*K*DPTCO;
0535
0536 Inc = 0;
0537 for (Ci = 1; Ci <= NbCu; Ci++) {
0538 dimen = tabdim->Value(Ci-1);
0539 PTCOXCI = PTCOX.Col(Ci);
0540 PTCOYCI = PTCOY.Col(Ci);
0541 PTCOZCI = PTCOZ.Col(Ci);
0542 PTCXCI = PTCX.Col(Ci);
0543 PTCYCI = PTCY.Col(Ci);
0544 PTCZCI = PTCZ.Col(Ci);
0545 Erruzax = (PTCXCI - PTCOXCI);
0546 Erruzay = (PTCYCI - PTCOYCI);
0547 Erruzaz = (PTCZCI - PTCOZCI);
0548 Scal = 0.0;
0549
0550 for (j = 1; j <= Npol ; j++) {
0551 Scal = (A(Pi, j)*Erruzax(j)) * (A(Pi, j)*DERR(j+Inc)) +
0552 (A(Pi, j)*Erruzay(j)) * (A(Pi, j)*DERR(j+Inc+Npol));
0553 if (dimen == 3) {
0554 Scal += (A(Pi, j)*Erruzax(j)) * (A(Pi, j)*DERR(j+Inc+2*Npol));
0555 }
0556 }
0557
0558 ValGrad_F(Pi) = ValGrad_F(Pi) + 2*Scal;
0559 if (dimen == 3) Inc = Inc +3*Npol;
0560 else Inc = Inc + 2*Npol;
0561 }
0562 }
0563 }
0564
0565 }
0566 }
0567
0568
0569 Standard_Integer AppParCurves_Function::NbVariables() const{
0570 return NbP;
0571 }
0572
0573
0574 Standard_Boolean AppParCurves_Function::Gradient (const math_Vector& X,
0575 math_Vector& G) {
0576
0577 Perform(X);
0578 G = ValGrad_F;
0579
0580 return Standard_True;
0581 }
0582
0583
0584 Standard_Boolean AppParCurves_Function::Values (const math_Vector& X,
0585 Standard_Real& F,
0586 math_Vector& G) {
0587
0588
0589 Perform(X);
0590 F = FVal;
0591 G = ValGrad_F;
0592 return Standard_True;
0593 }
0594
0595
0596 const AppParCurves_MultiCurve& AppParCurves_Function::CurveValue() {
0597 if (!Contraintes) MyMultiCurve = MyLeastSquare.BezierValue();
0598 return MyMultiCurve;
0599 }
0600
0601
0602 Standard_Real AppParCurves_Function::Error(const Standard_Integer IPoint,
0603 const Standard_Integer CurveIndex) const {
0604 return Sqrt(MyF(IPoint, CurveIndex));
0605 }
0606
0607 Standard_Real AppParCurves_Function::MaxError3d() const
0608 {
0609 return ERR3d;
0610 }
0611
0612 Standard_Real AppParCurves_Function::MaxError2d() const
0613 {
0614 return ERR2d;
0615 }
0616
0617
0618
0619 const math_Vector& AppParCurves_Function::NewParameters() const
0620 {
0621 return myParameters;
0622 }
