[occt.git] / src / AppParCurves / AppParCurves_Function.gxx

// Copyright (c) 1995-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

// Lpa, le 20/09/91


// Calcul de la valeur de F et grad_F, connaissant le parametrage.
// Cette fonction, appelee par le gradient conjugue, calcul F et 
// DF(ui, Poles(ui)) ce qui implique un calcul des nouveaux poles 
//  a chaque appel.

#define No_Standard_RangeError
#define No_Standard_OutOfRange


#include <AppParCurves_MultiCurve.hxx>
#include <AppParCurves_MultiPoint.hxx>
#include <TColStd_HArray1OfInteger.hxx>
#include <gp_Pnt.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Vec.hxx>
#include <gp_Vec2d.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array1OfPnt2d.hxx>
#include <AppParCurves_ConstraintCouple.hxx>

AppParCurves_Function::
  AppParCurves_Function(const MultiLine& SSP,
	 const Standard_Integer FirstPoint,
	 const Standard_Integer LastPoint,
	 const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
	 const math_Vector& Parameters,
	 const Standard_Integer Deg) :
         MyMultiLine(SSP),
         MyMultiCurve(Deg+1),          
         myParameters(Parameters.Lower(), Parameters.Upper()),
         ValGrad_F(FirstPoint, LastPoint),
         MyF(FirstPoint, LastPoint, 
	     1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
         PTLX(FirstPoint, LastPoint, 
	     1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
         PTLY(FirstPoint, LastPoint, 
	     1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
         PTLZ(FirstPoint, LastPoint, 
	     1, ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP), 0.0),
         A(FirstPoint, LastPoint, 1, Deg+1),
         DA(FirstPoint, LastPoint, 1, Deg+1),
         MyLeastSquare(SSP, FirstPoint, LastPoint, 
		       FirstConstraint(TheConstraints, FirstPoint),
		       LastConstraint(TheConstraints, LastPoint), Deg+1)
{
  Standard_Integer i;
  for (i=Parameters.Lower(); i<=Parameters.Upper();i++)
    myParameters(i)=Parameters(i);
  FirstP = FirstPoint;
  LastP = LastPoint;
  myConstraints = TheConstraints;
  NbP = LastP-FirstP+1;
  Adeb = FirstP;
  Afin = LastP;
  Degre = Deg;
  Contraintes = Standard_False;
  Standard_Integer myindex;
  Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
  AppParCurves_ConstraintCouple mycouple;
  AppParCurves_Constraint Cons;

  for (i = low; i <= upp; i++) {
    mycouple = TheConstraints->Value(i);
    Cons = mycouple.Constraint();
    myindex = mycouple.Index();
    if (myindex == FirstP) {
      if (Cons >= 1) Adeb = Adeb+1;
    }
    else if (myindex == LastP) {
      if (Cons >= 1) Afin = Afin-1;
    }
    else {
      if (Cons >= 1) Contraintes = Standard_True;
    }
  }

  Standard_Integer nb3d = ToolLine::NbP3d(SSP);
  Standard_Integer nb2d = ToolLine::NbP2d(SSP);
  Standard_Integer mynb3d= nb3d, mynb2d=nb2d;
  if (nb3d == 0) mynb3d = 1;
  if (nb2d == 0) mynb2d = 1;

  NbCu = nb3d+nb2d;
  tabdim = new TColStd_HArray1OfInteger(0, NbCu-1);

  if (Contraintes) {
    for (i = 1; i <= NbCu; i++) {
      if (i <= nb3d) tabdim->SetValue(i-1, 3);
      else tabdim->SetValue(i-1, 2);
    }

    TColgp_Array1OfPnt TabP(1, mynb3d);
    TColgp_Array1OfPnt2d TabP2d(1, mynb2d);
    
    for ( i = FirstP; i <= LastP; i++) {
      if (nb3d != 0 && nb2d != 0) ToolLine::Value(SSP, i, TabP, TabP2d);
      else if (nb3d != 0)         ToolLine::Value(SSP, i, TabP);
      else                        ToolLine::Value(SSP, i, TabP2d);
      for (Standard_Integer j = 1; j <= NbCu; j++) {
	if (tabdim->Value(j-1) == 3) {
	  TabP(j).Coord(PTLX(i, j), PTLY(i, j),PTLZ(i, j));
	}
	else {
	  TabP2d(j).Coord(PTLX(i, j), PTLY(i, j));
	}
      }
    }
  }
}


AppParCurves_Constraint AppParCurves_Function::FirstConstraint
  (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
   const Standard_Integer FirstPoint) const
{
  Standard_Integer i, myindex;
  Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
  AppParCurves_ConstraintCouple mycouple;
  AppParCurves_Constraint Cons = AppParCurves_NoConstraint;

  for (i = low; i <= upp; i++) {
    mycouple = TheConstraints->Value(i);
    Cons = mycouple.Constraint();
    myindex = mycouple.Index();
    if (myindex == FirstPoint) {
      break;
    }
  }
  return Cons;
}


AppParCurves_Constraint AppParCurves_Function::LastConstraint
  (const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
   const Standard_Integer LastPoint) const
{
  Standard_Integer i, myindex;
  Standard_Integer low = TheConstraints->Lower(), upp= TheConstraints->Upper();
  AppParCurves_ConstraintCouple mycouple;
  AppParCurves_Constraint Cons = AppParCurves_NoConstraint;

  for (i = low; i <= upp; i++) {
    mycouple = TheConstraints->Value(i);
    Cons = mycouple.Constraint();
    myindex = mycouple.Index();
    if (myindex == LastPoint) {
      break;
    }
  }
  return Cons;
}


Standard_Boolean AppParCurves_Function::Value (const math_Vector& X, 
					       Standard_Real& F) {

  myParameters = X;

  // Resolution moindres carres:
  // ===========================
  MyLeastSquare.Perform(myParameters);
  if (!(MyLeastSquare.IsDone())) { 
    Done = Standard_False;
    return Standard_False;
  }
  if (!Contraintes) {
    MyLeastSquare.Error(FVal, ERR3d, ERR2d);
    F = FVal;
  }

  // Resolution avec contraintes:
  // ============================
  else { 
    Standard_Integer Npol = Degre+1;
//    Standard_Boolean Ext = Standard_True;
    Standard_Integer Ci, i, j, dimen;
    Standard_Real AA, BB, CC, AIJ, FX, FY, FZ, Fi;
    math_Vector PTCXCI(1, Npol), PTCYCI(1, Npol), PTCZCI(1, Npol);
    ERR3d = ERR2d = 0.0;
    
    MyMultiCurve = MyLeastSquare.BezierValue();

    A = MyLeastSquare.FunctionMatrix();
    ResolCons Resol(MyMultiLine, MyMultiCurve, FirstP, LastP, myConstraints,
		    A, MyLeastSquare.DerivativeFunctionMatrix());
    if (!Resol.IsDone()) {
      Done = Standard_False;
      return Standard_False;
    }

    // Calcul de F = Sum||C(ui)-Ptli||2  sur toutes les courbes :
    // ========================================================================
    FVal = 0.0;
    
    for (Ci = 1; Ci <= NbCu; Ci++) {
      dimen = tabdim->Value(Ci-1);
      for (j = 1; j <= Npol; j++) {
	if (dimen == 3){ 
	  MyMultiCurve.Value(j).Point(Ci).Coord(PTCXCI(j),PTCYCI(j),PTCZCI(j));
	}
	else{ 
	  MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCXCI(j), PTCYCI(j));
	}
      }
      
      // Calcul de F:
      // ============
      for (i = Adeb; i <= Afin; i++) {
	AA = 0.0; BB = 0.0; CC = 0.0;
	for (j = 1; j <= Npol; j++) {
	  AIJ = A(i, j);
	  AA += AIJ*PTCXCI(j);
	  BB += AIJ*PTCYCI(j);
	  if (dimen == 3) { 
	    CC += AIJ*PTCZCI(j);
	  }
	}
	FX = AA-PTLX(i, Ci);
	FY = BB-PTLY(i, Ci);
	MyF(i,Ci) = FX*FX + FY*FY;
	if (dimen == 3) {
	  FZ = CC-PTLZ(i,Ci);
	  MyF(i, Ci) += FZ*FZ;
	  Fi = MyF(i, Ci);
	  if (Sqrt(Fi) > ERR3d) ERR3d = Sqrt(Fi);
	}
	else {
	  Fi = MyF(i, Ci);
	  if (Sqrt(Fi) > ERR2d) ERR2d = Sqrt(Fi);
	}
	FVal += Fi;
      }
    }
    F = FVal;
  }  
  return Standard_True;
}


void AppParCurves_Function::Perform(const math_Vector& X) {
  Standard_Integer j;

  myParameters = X;
  // Resolution moindres carres:
  // ===========================
  MyLeastSquare.Perform(myParameters);

  if (!(MyLeastSquare.IsDone())) { 
    Done = Standard_False;
    return;
  }

  for(j = myParameters.Lower(); j <= myParameters.Upper(); j++) {
    ValGrad_F(j) = 0.0;
  }

  if (!Contraintes) {
    MyLeastSquare.ErrorGradient(ValGrad_F, FVal, ERR3d, ERR2d);
  }
  else {
    Standard_Integer Pi, Ci, i, k, dimen;
    Standard_Integer  Npol = Degre+1;
    Standard_Real Scal, AA, BB, CC, DAA, DBB, DCC;
    Standard_Real FX, FY, FZ, AIJ, DAIJ, px, py, pz, Fi;
    AppParCurves_Constraint Cons=AppParCurves_NoConstraint;
    math_Matrix Grad_F(FirstP, LastP, 1, NbCu, 0.0);
    math_Vector PTCXCI(1, Npol), PTCYCI(1, Npol), PTCZCI(1, Npol);
    math_Vector PTCOXCI(1, Npol), PTCOYCI(1, Npol), PTCOZCI(1, Npol);
//    Standard_Boolean Ext = Standard_True;
    ERR3d = ERR2d = 0.0;

    math_Matrix PTCOX(1, Npol, 1, NbCu), PTCOY(1, Npol, 1, NbCu), 
                PTCOZ(1, Npol,1, NbCu);
    math_Matrix PTCX(1, Npol, 1, NbCu), PTCY(1, Npol, 1, NbCu), 
                PTCZ(1, Npol,1, NbCu);
    Standard_Integer Inc;

    MyMultiCurve = MyLeastSquare.BezierValue();

    for (Ci =1; Ci <= NbCu; Ci++) {
      dimen = tabdim->Value(Ci-1);
      for (j = 1; j <= Npol; j++) {
	if (dimen == 3){ 
	  MyMultiCurve.Value(j).Point(Ci).Coord(PTCOX(j, Ci),
						PTCOY(j, Ci),
						PTCOZ(j, Ci));
	}
	else{ 
	  MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCOX(j, Ci), PTCOY(j, Ci));
	  PTCOZ(j, Ci) = 0.0;
	}
      }
    }

    A = MyLeastSquare.FunctionMatrix();
    DA = MyLeastSquare.DerivativeFunctionMatrix();
    
    // Resolution avec contraintes:
    // ============================
    
    ResolCons Resol(MyMultiLine, MyMultiCurve, FirstP, LastP, 
		    myConstraints, A, DA);
    if (!Resol.IsDone()) {
      Done = Standard_False;
      return;
    }
    
    
    // Calcul de F = Sum||C(ui)-Ptli||2 et du gradient non contraint de F pour
    // chaque point PointIndex.
    // ========================================================================
    FVal = 0.0;
    for(j = FirstP; j <= LastP; j++) {
      ValGrad_F(j) = 0.0;
    }

    math_Matrix TrA(A.LowerCol(), A.UpperCol(), A.LowerRow(), A.UpperRow());
    math_Matrix TrDA(DA.LowerCol(), DA.UpperCol(), DA.LowerRow(), DA.UpperRow());
    math_Matrix RESTM(A.LowerCol(), A.UpperCol(), A.LowerCol(), A.UpperCol());

    const math_Matrix& K = Resol.ConstraintMatrix();
    const math_Matrix& DK = Resol.ConstraintDerivative(MyMultiLine, X, Degre, DA);
    math_Matrix TK(K.LowerCol(), K.UpperCol(), K.LowerRow(), K.UpperRow());
    TK = K.Transposed();
    const math_Vector& Vardua = Resol.Duale();
    math_Matrix KK(K.LowerCol(), K.UpperCol(), Vardua.Lower(), Vardua.Upper());
    KK = (K.Transposed())*(Resol.InverseMatrix());
    math_Matrix DTK(DK.LowerCol(), DK.UpperCol(), DK.LowerRow(), DK.UpperRow());
    DTK = DK.Transposed();
    TrA = A.Transposed();
    TrDA = DA.Transposed();
    RESTM = ((A.Transposed()*A).Inverse());

    math_Vector DPTCO(1, K.ColNumber());
    math_Matrix DPTCO1(FirstP, LastP, 1, K.ColNumber());
    math_Vector DKPTC(1, K.RowNumber());


    FVal = 0.0;
    for (Ci = 1; Ci <= NbCu; Ci++) {
      dimen = tabdim->Value(Ci-1);
      for (j = 1; j <= Npol; j++) {
	if (dimen == 3){ 
	  MyMultiCurve.Value(j).Point(Ci).Coord(PTCX(j, Ci), 
						PTCY(j, Ci), 
						PTCZ(j, Ci));
	}
	else{ 
	  MyMultiCurve.Value(j).Point2d(Ci).Coord(PTCX(j, Ci), PTCY(j,Ci));
	  PTCZ(j, Ci) = 0.0;
	}
      }
    }

    
    // Calcul du gradient sans contraintes:
    // ====================================

    for (Ci = 1; Ci <= NbCu; Ci++) {
      dimen = tabdim->Value(Ci-1);
      for (i = Adeb; i <= Afin; i++) {
	AA = 0.0; BB = 0.0; CC = 0.0; DAA = 0.0; DBB = 0.0; DCC = 0.0;
	for (j = 1; j <= Npol; j++) {
	  AIJ = A(i, j); DAIJ = DA(i, j);
	  px = PTCX(j, Ci); py = PTCY(j, Ci);
	  AA += AIJ*px;  BB += AIJ*py;
	  DAA += DAIJ*px;  DBB += DAIJ*py;
	  if (dimen == 3) { 
	    pz = PTCZ(j, Ci);
	    CC += AIJ*pz;  DCC += DAIJ*pz;
	  }
	}
	FX = AA-PTLX(i, Ci);
	FY = BB-PTLY(i, Ci);
	MyF(i,Ci) = FX*FX + FY*FY;
	Grad_F(i, Ci) = 2.0*(DAA*FX + DBB*FY);
	if (dimen == 3) {
	  FZ = CC-PTLZ(i,Ci);
	  MyF(i, Ci) += FZ*FZ;
	  Grad_F(i, Ci) += 2.0*DCC*FZ;
	  Fi = MyF(i, Ci);
	  if (Sqrt(Fi) > ERR3d) ERR3d = Sqrt(Fi);
	}
	else {
	  Fi = MyF(i, Ci);
	  if (Sqrt(Fi) > ERR2d) ERR2d = Sqrt(Fi);
	}
	FVal += Fi;
	ValGrad_F(i) += Grad_F(i, Ci);
      }
    }


    // Calcul de DK*PTC:
    // =================
    for (i = 1; i <= K.RowNumber(); i++) {
      Inc = 0;
      for (Ci = 1; Ci <= NbCu; Ci++) {
	dimen = tabdim->Value(Ci-1);
	DKPTC(i) = 0.0;
	for (j = 1; j <= Npol; j++) {
	  DKPTC(i) += DK(i, j+Inc)*PTCX(j, Ci)+ DK(i, j+Inc+Npol)*PTCY(j, Ci);
	  if (dimen == 3) {
	    DKPTC(i) += DK(i, j+Inc+2*Npol)*PTCZ(j, Ci);
	  }
	}
	if (dimen == 3) Inc = Inc +3*Npol;
	else Inc = Inc +2*Npol;
      }
    }
    
    math_Vector DERR(DTK.LowerRow(), DTK.UpperRow());
    DERR = (DTK)*Vardua-KK* ((DKPTC) + K*(DTK)*Vardua);

    // rajout du gradient avec contraintes:
    // ====================================
    // dPTCO1/duk = [d(TA)/duk*[A*PTCO-PTL] + TA*dA/duk*PTCO]


    Inc = 0;

    math_Vector Errx(A.LowerRow(), A.UpperRow());
    math_Vector Erry(A.LowerRow(), A.UpperRow());
    math_Vector Errz(A.LowerRow(), A.UpperRow());
    math_Vector Scalx(DA.LowerRow(), DA.UpperRow());
    math_Vector Scaly(DA.LowerRow(), DA.UpperRow());
    math_Vector Scalz(DA.LowerRow(), DA.UpperRow());
    math_Vector Erruzax(PTCXCI.Lower(), PTCXCI.Upper());
    math_Vector Erruzay(PTCYCI.Lower(), PTCYCI.Upper());
    math_Vector Erruzaz(PTCZCI.Lower(), PTCZCI.Upper());
    math_Vector TrDAPI(TrDA.LowerRow(), TrDA.UpperRow());
    math_Vector TrAPI(TrA.LowerRow(), TrA.UpperRow());

    for (Ci = 1; Ci <= NbCu; Ci++) {
      dimen = tabdim->Value(Ci-1);
      PTCOXCI = PTCOX.Col(Ci);
      PTCOYCI = PTCOY.Col(Ci);
      PTCOZCI = PTCOZ.Col(Ci);
      PTCXCI = PTCX.Col(Ci);
      PTCYCI = PTCY.Col(Ci);
      PTCZCI = PTCZ.Col(Ci);


      Errx = (A*PTCOXCI - PTLX.Col(Ci));
      Erry = (A*PTCOYCI - PTLY.Col(Ci));
      Errz = (A*PTCOZCI - PTLZ.Col(Ci));
      Scalx = (DA*PTCOXCI);   // Scal = DA * PTCO
      Scaly = (DA*PTCOYCI);
      Scalz = (DA*PTCOZCI);
      Erruzax = (PTCXCI - PTCOXCI);
      Erruzay = (PTCYCI - PTCOYCI);
      Erruzaz = (PTCZCI - PTCOZCI);
      
      for (Pi = FirstP; Pi <= LastP; Pi++) {
	TrDAPI = (TrDA.Col(Pi));
	TrAPI = (TrA.Col(Pi));
	Standard_Real Taa = TrAPI*A.Row(Pi);
	Scal = 0.0;
	for (j = 1; j <= Npol; j++) {
	  DPTCO1(Pi, j + Inc) = (TrDAPI*Errx(Pi)+TrAPI*Scalx(Pi))(j);
	  DPTCO1(Pi, j + Inc+ Npol) = (TrDAPI*Erry(Pi)+TrAPI*Scaly(Pi))(j);
	  Scal += DPTCO1(Pi, j+Inc)* Taa*Erruzax(j) + DPTCO1(Pi, j+Inc+Npol)*Taa*Erruzay(j);
	  if (dimen == 3) {
	    DPTCO1(Pi, j + Inc+ 2*Npol) = (TrDAPI*Errz(Pi)+TrAPI*Scalz(Pi))(j);
	    Scal += DPTCO1(Pi, j+Inc+2*Npol)*Taa*Erruzaz(j);
	  }
	}
	ValGrad_F(Pi) = ValGrad_F(Pi) - 2*Scal;
      }
      if (dimen == 3) Inc = Inc + 3*Npol;
      else Inc = Inc +2*Npol;
    }


    // on calcule DPTCO = - RESTM * DPTCO1:
    
    // Calcul de DPTCO/duk:
    // dPTCO/duk = -Inv(T(A)*A)*[d(TA)/duk*[A*PTCO-PTL] + TA*dA/duk*PTCO]

    Standard_Integer low=myConstraints->Lower(), upp=myConstraints->Upper();
    Inc = 0;
    for (Pi = FirstP; Pi <= LastP; Pi++) {
      for (i = low; i <= upp; i++) {
	if (myConstraints->Value(i).Index() == Pi) {
	  Cons = myConstraints->Value(i).Constraint();
	  break;
	}
      }
      if (Cons >= 1) {
	Inc = 0;
	for (Ci = 1; Ci <= NbCu; Ci++) {
	  dimen = tabdim->Value(Ci-1);
	  for (j = 1; j <= Npol; j++) {
	    DPTCO(j+Inc) = 0.0;
	    DPTCO(j+Inc+Npol) = 0.0;
	    if (dimen == 3) DPTCO(j+Inc+2*Npol) = 0.0;
	    for (k = 1; k <= Npol; k++) {
	      DPTCO(j+Inc) = DPTCO(j+Inc) -RESTM(j, k) * DPTCO1(Pi, j+Inc);
	      DPTCO(j+Inc+Npol)=DPTCO(j+Inc+Npol)-RESTM(j, k)*DPTCO1(Pi,j+Inc+Npol);
	      if (dimen == 3) {
		DPTCO(j+Inc+2*Npol) = DPTCO(j+Inc+2*Npol) 
		  -RESTM(j, k) * DPTCO1(Pi, j+Inc+2*Npol);
	      }
	    }
	  }
	  if (dimen == 3) Inc += 3*Npol;
	  else Inc += 2*Npol;
	}
	
	DERR = DERR-KK*K*DPTCO;
	
	Inc = 0;
	for (Ci = 1; Ci <= NbCu; Ci++) {
	  dimen = tabdim->Value(Ci-1);
	  PTCOXCI = PTCOX.Col(Ci);
	  PTCOYCI = PTCOY.Col(Ci);
	  PTCOZCI = PTCOZ.Col(Ci);
	  PTCXCI = PTCX.Col(Ci);
	  PTCYCI = PTCY.Col(Ci);
	  PTCZCI = PTCZ.Col(Ci);
	  Erruzax = (PTCXCI - PTCOXCI);
	  Erruzay = (PTCYCI - PTCOYCI);
	  Erruzaz = (PTCZCI - PTCOZCI);
	  Scal = 0.0;
	  
	  for (j = 1; j <= Npol ; j++) {
	    Scal = (A(Pi, j)*Erruzax(j)) * (A(Pi, j)*DERR(j+Inc)) + 
	      (A(Pi, j)*Erruzay(j)) * (A(Pi, j)*DERR(j+Inc+Npol));
	    if (dimen == 3) {
	      Scal += (A(Pi, j)*Erruzax(j)) * (A(Pi, j)*DERR(j+Inc+2*Npol)); 
	    }
	  }
	  
	  ValGrad_F(Pi) = ValGrad_F(Pi) + 2*Scal;
	  if (dimen == 3) Inc = Inc +3*Npol;
	  else Inc = Inc + 2*Npol;
	}
      }
    }
    
  }
}


Standard_Integer AppParCurves_Function::NbVariables() const{ 
  return NbP;
}


Standard_Boolean AppParCurves_Function::Gradient (const math_Vector& X,
						  math_Vector& G) {

  Perform(X);
  G = ValGrad_F;

  return Standard_True;
}


Standard_Boolean AppParCurves_Function::Values (const math_Vector& X, 
						Standard_Real& F, 
						math_Vector& G) {


  Perform(X);
  F = FVal;
  G = ValGrad_F;
  return Standard_True;
}


const AppParCurves_MultiCurve& AppParCurves_Function::CurveValue() {
  if (!Contraintes)  MyMultiCurve = MyLeastSquare.BezierValue();
  return MyMultiCurve;
}


Standard_Real AppParCurves_Function::Error(const Standard_Integer IPoint,
				     const Standard_Integer CurveIndex) const {
  return Sqrt(MyF(IPoint, CurveIndex));
}

Standard_Real AppParCurves_Function::MaxError3d() const
{
  return ERR3d;
}

Standard_Real AppParCurves_Function::MaxError2d() const
{
  return ERR2d;
}


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