// Copyright (c) 1996-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 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.

// Jeannine PANCIATICI le 06/06/96
// Igor FEOKTISTOV 14/12/98 - correction of Approximate() and Init().
// Approximation d une MultiLine de points decrite par le tool MLineTool.
// avec criteres variationnels


#define No_Standard_RangeError
#define No_Standard_OutOfRange
#define No_Standard_DimensionError
#define No_Standard_ConstructionError

#include <Standard_Stream.hxx>

#include <AppParCurves.hxx>
#include <AppParCurves_Constraint.hxx>
#include <AppParCurves_HArray1OfConstraintCouple.hxx>
#include <AppParCurves_Array1OfMultiPoint.hxx>
#include <AppParCurves_MultiPoint.hxx>
#include <AppParCurves_MultiBSpCurve.hxx>
#include <Convert_CompPolynomialToPoles.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 <TColgp_Array1OfVec.hxx>
#include <TColgp_Array1OfVec2d.hxx>
#include <TColStd_Array1OfInteger.hxx>
#include <TColStd_HArray1OfInteger.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <TColStd_HArray2OfReal.hxx>
#include <StdFail_NotDone.hxx>
#include <Standard_SStream.hxx>
#include <Standard_NoSuchObject.hxx>
#include <Precision.hxx>
//#include <Smoothing.h>

#if defined(WNT)
# include <stdio.h>
# include <stdarg.h>
#endif  /* WNT */


//
//=======================================================================
//function : AppParCurves_Variational
//purpose  : Initialization of   the   fields.
//=======================================================================
//
 AppParCurves_Variational::AppParCurves_Variational(const MultiLine& SSP,
                                                    const Standard_Integer FirstPoint, 
                                                    const Standard_Integer LastPoint,
                                                    const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
                                                    const Standard_Integer MaxDegree,
                                                    const Standard_Integer MaxSegment,
                                                    const GeomAbs_Shape Continuity, 
                                                    const Standard_Boolean WithMinMax, 
                                                    const Standard_Boolean WithCutting, 
                                                    const Standard_Real Tolerance, 
                                                    const Standard_Integer NbIterations):
                                        mySSP(SSP),
                                        myFirstPoint(FirstPoint),
                                        myLastPoint(LastPoint),
                                        myConstraints(TheConstraints),
                                        myMaxDegree(MaxDegree),
                                        myMaxSegment(MaxSegment),
                                        myNbIterations(NbIterations),
                                        myTolerance(Tolerance),
                                        myContinuity(Continuity),
                                        myWithMinMax(WithMinMax),
                                        myWithCutting(WithCutting)


{
// Verifications:
    if (myMaxDegree < 1) Standard_DomainError::Raise();
    myMaxDegree = Min (30, myMaxDegree);
//
    if (myMaxSegment < 1) Standard_DomainError::Raise();
//
    if (myWithMinMax != 0 && myWithMinMax !=1 ) Standard_DomainError::Raise();
    if (myWithCutting != 0 && myWithCutting !=1 ) Standard_DomainError::Raise();
//
    myIsOverConstr = Standard_False;
    myIsCreated    = Standard_False;
    myIsDone       = Standard_False;  
    switch (myContinuity) {
    case GeomAbs_C0:
       myNivCont=0;
       break ;
    case GeomAbs_C1: 
       myNivCont=1;
       break ;
    case GeomAbs_C2:  
       myNivCont=2;
       break ;
    default:
    Standard_ConstructionError::Raise();
  }
//
    myNbP2d = ToolLine::NbP2d(SSP); 
    myNbP3d = ToolLine::NbP3d(SSP);
    myDimension = 2 * myNbP2d + 3* myNbP3d ;
//
    myPercent[0]=0.4;
    myPercent[1]=0.2;
    myPercent[2]=0.4;
    myKnots= new TColStd_HArray1OfReal(1,2);
    myKnots->SetValue(1,0.);
    myKnots->SetValue(2,1.);
    
//  Declaration
//  
    mySmoothCriterion = new AppParCurves_MyCriterion(mySSP, myFirstPoint, myLastPoint);
    myParameters = new TColStd_HArray1OfReal(myFirstPoint, myLastPoint);
    myNbPoints=myLastPoint-myFirstPoint+1;
    if (myNbPoints <= 0) Standard_ConstructionError::Raise();
// 
    myTabPoints= new TColStd_HArray1OfReal(1,myDimension*myNbPoints);
//
//  Table of Points initialization
//
    Standard_Integer ipoint,jp2d,jp3d,index;
    TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
    TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));    
    gp_Pnt2d P2d;
    gp_Pnt P3d;
    index=1;
    
    for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
      {

	if(myNbP2d !=0 && myNbP3d ==0 ) {
                          ToolLine::Value(mySSP,ipoint,TabP2d);

                          for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
			      {  P2d = TabP2d.Value(jp2d);
                                 
                                 myTabPoints->SetValue(index++,P2d.X());
                                 myTabPoints->SetValue(index++,P2d.Y());
			      }
                         } 
        if(myNbP3d !=0 && myNbP2d == 0 ) {
                          ToolLine::Value(mySSP,ipoint,TabP3d);
                          
                          for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)

			      {  P3d=TabP3d.Value(jp3d);
                                 
                                 myTabPoints->SetValue(index++,P3d.X());
                                 myTabPoints->SetValue(index++,P3d.Y()); 
                                 myTabPoints->SetValue(index++,P3d.Z());

			      }
                         }
        if(myNbP3d !=0 && myNbP2d != 0 ) {
                          ToolLine::Value(mySSP,ipoint,TabP3d,TabP2d);
                          
                          for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)

			      {  P3d=TabP3d.Value(jp3d);
                                 
                                 myTabPoints->SetValue(index++,P3d.X());
                                 myTabPoints->SetValue(index++,P3d.Y()); 
                                 myTabPoints->SetValue(index++,P3d.Z());

			      }
                          for ( jp2d = 1 ; jp2d <= myNbP2d ; jp2d++)

			      {  P2d=TabP2d.Value(jp2d);
                                 
                                 myTabPoints->SetValue(index++,P2d.X());
                                 myTabPoints->SetValue(index++,P2d.Y()); 
			      }
                         }
      }
    Init();
}
//
//=======================================================================
//function : Init
//purpose  : Initializes the tables of constraints
//           and verifies if the problem is not over-constrained.
//           This method is used in the Create and the method SetConstraint. 
//=======================================================================
//
void AppParCurves_Variational::Init()
{
			      
    Standard_Integer ipoint,jp2d,jp3d,index,jndex;
    Standard_Integer CurMultyPoint;
    TColgp_Array1OfVec TabV3d(1, Max(1,myNbP3d));
    TColgp_Array1OfVec2d TabV2d(1, Max(1,myNbP2d));
    TColgp_Array1OfVec TabV3dcurv(1, Max(1,myNbP3d));
    TColgp_Array1OfVec2d TabV2dcurv(1, Max(1,myNbP2d));

    gp_Vec Vt3d, Vc3d;
    gp_Vec2d Vt2d, Vc2d;

    myNbConstraints=myConstraints->Length();
    if (myNbConstraints < 0) Standard_ConstructionError::Raise(); 

    myTypConstraints = new TColStd_HArray1OfInteger(1,Max(1,2*myNbConstraints));
    myTabConstraints = new TColStd_HArray1OfReal(1,Max(1,2*myDimension*myNbConstraints));
    myTtheta = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));
    myTfthet = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));

	
//
// Table of types initialization
      Standard_Integer iconstr;
      index=1;
      jndex=1;
      CurMultyPoint = 1;
      myNbPassPoints=0;
      myNbTangPoints=0;
      myNbCurvPoints=0;
      AppParCurves_Constraint valcontr;

      for ( iconstr = myConstraints->Lower() ; iconstr <= myConstraints->Upper() ; iconstr++)
	{
        ipoint=(myConstraints->Value(iconstr)).Index();

          valcontr=(myConstraints->Value(iconstr)).Constraint();
          switch (valcontr) {
                             case AppParCurves_NoConstraint:
	                          CurMultyPoint -= myNbP3d * 6 + myNbP2d * 2;
                                  break ;
                             case AppParCurves_PassPoint: 
                                  myTypConstraints->SetValue(index++,ipoint);                                    
                                  myTypConstraints->SetValue(index++,0);
                                  myNbPassPoints++;
                                  if(myNbP2d !=0 ) jndex=jndex+4*myNbP2d;
                                  if(myNbP3d !=0 ) jndex=jndex+6*myNbP3d;
                                  break ;
			     case AppParCurves_TangencyPoint:
                                  myTypConstraints->SetValue(index++,ipoint);
                                  myTypConstraints->SetValue(index++,1); 
                                  myNbTangPoints++;
	                          if(myNbP2d !=0 && myNbP3d == 0 ) 
                                    {
                                     if (ToolLine::Tangency(mySSP,ipoint,TabV2d) == Standard_False)
                                          Standard_ConstructionError::Raise();
                                     for (jp2d=1;jp2d<=myNbP2d;jp2d++)
                  			 {  
                                          Vt2d=TabV2d.Value(jp2d);
                                          Vt2d.Normalize();
                                          myTabConstraints->SetValue(jndex++,Vt2d.X());
                                          myTabConstraints->SetValue(jndex++,Vt2d.Y());
					  jndex=jndex+2;
					  InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
					 }
                                     } 
                                   if(myNbP3d !=0 && myNbP2d == 0) 
                                     {
                                      if (ToolLine::Tangency(mySSP,ipoint,TabV3d) == Standard_False)
                                          Standard_ConstructionError::Raise();
                                      for (jp3d=1;jp3d<=myNbP3d;jp3d++)
		                	  {  
                                           Vt3d=TabV3d.Value(jp3d);
                                           Vt3d.Normalize();    
                                           myTabConstraints->SetValue(jndex++,Vt3d.X());
                                                 
                                           myTabConstraints->SetValue(jndex++,Vt3d.Y()); 
                                                 
                                           myTabConstraints->SetValue(jndex++,Vt3d.Z());
					   jndex=jndex+3;
					   InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
                 			  }
                                      }
                                   if(myNbP3d !=0 && myNbP2d != 0) 
                                     {
                                      if (ToolLine::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False)
                                          Standard_ConstructionError::Raise();
                                      for (jp3d=1;jp3d<=myNbP3d;jp3d++)
		                	  {  
                                           Vt3d=TabV3d.Value(jp3d);
                                           Vt3d.Normalize();    
                                           myTabConstraints->SetValue(jndex++,Vt3d.X());
                                           myTabConstraints->SetValue(jndex++,Vt3d.Y()); 
                                           myTabConstraints->SetValue(jndex++,Vt3d.Z());
					   jndex=jndex+3;
					   InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
                 			  }

                                     for (jp2d=1;jp2d<=myNbP2d;jp2d++)
                  			 {  
                                          Vt2d=TabV2d.Value(jp2d);
                                          Vt2d.Normalize();
                                          myTabConstraints->SetValue(jndex++,Vt2d.X());
                                          myTabConstraints->SetValue(jndex++,Vt2d.Y());
					  jndex=jndex+2;
					  InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
			                 }
                                     }

                                       
                                       break ;

                             case AppParCurves_CurvaturePoint:
                                   myTypConstraints->SetValue(index++,ipoint);
                                   myTypConstraints->SetValue(index++,2);
                                   myNbCurvPoints++;
	                           if(myNbP2d !=0 && myNbP3d == 0) 
                                       {
                                       if (ToolLine::Tangency(mySSP,ipoint,TabV2d) == Standard_False )
                                           Standard_ConstructionError::Raise();
                                       if (ToolLine::Curvature(mySSP,ipoint,TabV2dcurv) == Standard_False)
                                           Standard_ConstructionError::Raise();
                                       for (jp2d=1;jp2d<=myNbP2d;jp2d++)
                  			   {  
                                            Vt2d=TabV2d.Value(jp2d);
                                            Vt2d.Normalize();
                                            Vc2d=TabV2dcurv.Value(jp2d);
                                            if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
                                                 Standard_ConstructionError::Raise();
                                            myTabConstraints->SetValue(jndex++,Vt2d.X());
                                            myTabConstraints->SetValue(jndex++,Vt2d.Y());
                                            myTabConstraints->SetValue(jndex++,Vc2d.X());
                                            myTabConstraints->SetValue(jndex++,Vc2d.Y());
					    InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
			                   }
                                       } 

                                    if(myNbP3d !=0 && myNbP2d == 0 ) 
                                       {
                                       if (ToolLine::Tangency(mySSP,ipoint,TabV3d) == Standard_False )
                                           Standard_ConstructionError::Raise();
				       if (ToolLine::Curvature(mySSP,ipoint,TabV3dcurv) == Standard_False)
					   Standard_ConstructionError::Raise();
                                       for (jp3d=1;jp3d<=myNbP3d;jp3d++)
		                	   {  
                                            Vt3d=TabV3d.Value(jp3d);
                                            Vt3d.Normalize();
                                            Vc3d=TabV3dcurv.Value(jp3d);
                                            if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False) 
                                                 Standard_ConstructionError::Raise();
                                            myTabConstraints->SetValue(jndex++,Vt3d.X());
                                            myTabConstraints->SetValue(jndex++,Vt3d.Y()); 
                                            myTabConstraints->SetValue(jndex++,Vt3d.Z());
                                            myTabConstraints->SetValue(jndex++,Vc3d.X());
                                            myTabConstraints->SetValue(jndex++,Vc3d.Y()); 
                                            myTabConstraints->SetValue(jndex++,Vc3d.Z());
					    InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
                 	        	   }
                                       }
                                    if(myNbP3d !=0 && myNbP2d != 0 ) 
                                       {
                                       if (ToolLine::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False )
                                           Standard_ConstructionError::Raise();
				       if (ToolLine::Curvature(mySSP,ipoint,TabV3dcurv,TabV2dcurv) == Standard_False)
                                           Standard_ConstructionError::Raise();
                                       for (jp3d=1;jp3d<=myNbP3d;jp3d++)
		                	   {  
                                            Vt3d=TabV3d.Value(jp3d);
                                            Vt3d.Normalize();
                                            Vc3d=TabV3dcurv.Value(jp3d);
                                            if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False) 
                                                 Standard_ConstructionError::Raise();
                                            myTabConstraints->SetValue(jndex++,Vt3d.X());
                                            myTabConstraints->SetValue(jndex++,Vt3d.Y()); 
                                            myTabConstraints->SetValue(jndex++,Vt3d.Z());
                                            myTabConstraints->SetValue(jndex++,Vc3d.X());
                                            myTabConstraints->SetValue(jndex++,Vc3d.Y()); 
                                            myTabConstraints->SetValue(jndex++,Vc3d.Z());
					    InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
                 	        	   }
                                       for (jp2d=1;jp2d<=myNbP2d;jp2d++)
                  			   {  
                                            Vt2d=TabV2d.Value(jp2d);
                                            Vt2d.Normalize();
                                            Vc2d=TabV2dcurv.Value(jp2d);
                                            if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
                                                 Standard_ConstructionError::Raise();
                                            myTabConstraints->SetValue(jndex++,Vt2d.X());
                                            myTabConstraints->SetValue(jndex++,Vt2d.Y());
                                            myTabConstraints->SetValue(jndex++,Vc2d.X());
                                            myTabConstraints->SetValue(jndex++,Vc2d.Y());
					    InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
			                   }

				      }
                                      break ; 
                                 default:
                                      Standard_ConstructionError::Raise();
                                }
	  CurMultyPoint += myNbP3d * 6 + myNbP2d * 2;
           }
// OverConstraint Detection
        Standard_Integer MaxSeg;
        if(myWithCutting == Standard_True) MaxSeg = myMaxSegment ;
        else MaxSeg = 1;
        if (((myMaxDegree-myNivCont)*MaxSeg-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
             {
               myIsOverConstr =Standard_True;
               myIsCreated = Standard_False;
             }
        else {
	  InitSmoothCriterion();
	  myIsCreated = Standard_True;
	}
              
          
}
//
//=======================================================================
//function : Approximate
//purpose  : Makes the approximation with the current fields.
//=======================================================================
//
void AppParCurves_Variational::Approximate()

{
  if (myIsCreated == Standard_False )  StdFail_NotDone:: Raise();


  Standard_Real WQuadratic, WQuality;

  TColStd_Array1OfReal Ecarts(myFirstPoint, myLastPoint); 

  mySmoothCriterion->GetWeight(WQuadratic, WQuality);
  
  Handle(FEmTool_Curve) TheCurve;  

  mySmoothCriterion->GetCurve(TheCurve);

//---------------------------------------------------------------------

  TheMotor(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);


  if(myWithMinMax && myTolerance < myMaxError)
    Adjusting(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
    
//---------------------------------------------------------------------

  Standard_Integer jp2d,jp3d,index,ipole, 
                   NbElem = TheCurve->NbElements();

  TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
  TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
  Standard_Real debfin[2] = {-1., 1};
  
  gp_Pnt2d P2d;
  gp_Pnt P3d;
  
  index=0;

  {
    Handle(TColStd_HArray2OfReal) PolynomialIntervalsPtr =
      new TColStd_HArray2OfReal(1,NbElem,1,2) ;

    Handle(TColStd_HArray1OfInteger) NbCoeffPtr = 
      new TColStd_HArray1OfInteger(1, myMaxSegment);

    Standard_Integer size = myMaxSegment * (myMaxDegree + 1) * myDimension ;
    Handle(TColStd_HArray1OfReal) CoeffPtr = new TColStd_HArray1OfReal(1,size);

    CoeffPtr->Init(0.);

    Handle(TColStd_HArray1OfReal) IntervallesPtr = 
      new TColStd_HArray1OfReal(1, NbElem + 1);

    IntervallesPtr->ChangeArray1() = TheCurve->Knots();

    TheCurve->GetPolynom(CoeffPtr->ChangeArray1());

    Standard_Integer ii;
    
    for(ii = 1; ii <= NbElem; ii++)  
      NbCoeffPtr->SetValue(ii, TheCurve->Degree(ii)+1);
      

    for (ii = PolynomialIntervalsPtr->LowerRow() ;
	 ii <= PolynomialIntervalsPtr->UpperRow() ;ii++) 
      {
	PolynomialIntervalsPtr->SetValue(ii,1,debfin[0]) ;
	PolynomialIntervalsPtr->SetValue(ii,2,debfin[1]) ;
      }
/*   
	printf("\n =========== Parameters for converting\n");
	printf("nb_courbes : %d \n", NbElem);
	printf("tab_option[4] : %d \n", myNivCont);
	printf("myDimension : %d \n", myDimension);
	printf("myMaxDegree : %d \n", myMaxDegree);
	printf("\n NbcoeffPtr :\n");
	for(ii = 1; ii <= NbElem; ii++)  printf("NbCoeffPtr(%d) = %d \n", ii, NbCoeffPtr->Value(ii));
	size = NbElem*(myMaxDegree + 1) * myDimension;
	printf("\n CoeffPtr :\n");
	for(ii = 1; ii <= size; ii++)  printf("CoeffPtr(%d) = %.8e \n", ii, CoeffPtr->Value(ii));
 	printf("\n PolinimialIntervalsPtr :\n");
        for (ii = PolynomialIntervalsPtr->LowerRow() ;
	     ii <= PolynomialIntervalsPtr->UpperRow() ;ii++) 
            {
             printf(" %d : [%f, %f] \n", ii, PolynomialIntervalsPtr->Value(ii,1),
		    PolynomialIntervalsPtr->Value(ii,2)) ;
            }
 	printf("\n IntervallesPtr :\n");
        for (ii = IntervallesPtr->Lower();
	     ii <= IntervallesPtr->Upper() - 1; ii++) 
            {
             printf(" %d : [%f, %f] \n", ii, IntervallesPtr->Value(ii),
		    IntervallesPtr->Value(ii+1)) ;
            }
*/  
    Convert_CompPolynomialToPoles AConverter(NbElem, 
					     myNivCont,
					     myDimension,
					     myMaxDegree,
					     NbCoeffPtr,
					     CoeffPtr,
					     PolynomialIntervalsPtr,
					     IntervallesPtr) ; 
    if (AConverter.IsDone()) 
      {
	Handle(TColStd_HArray2OfReal) PolesPtr ;
	Handle(TColStd_HArray1OfInteger) Mults;
	Standard_Integer NbPoles=AConverter.NbPoles();
//	Standard_Integer Deg=AConverter.Degree();
	AppParCurves_Array1OfMultiPoint TabMU(1,NbPoles);
	AConverter.Poles(PolesPtr) ;
	AConverter.Knots(myKnots) ;
	AConverter.Multiplicities(Mults) ;
	
	for (ipole=PolesPtr->LowerRow();ipole<=PolesPtr->UpperRow();ipole++)
	  {
	    index=PolesPtr->LowerCol();
/*	    if(myNbP2d !=0 ) 
	      {
		for (jp2d=1;jp2d<=myNbP2d;jp2d++)
		  {
		    P2d.SetX(PolesPtr->Value(ipole,index++));
		    P2d.SetY(PolesPtr->Value(ipole,index++));  
		    TabP2d.SetValue(jp2d,P2d);
		  }
	      }*/
	    if(myNbP3d !=0 ) 
	      {
		for (jp3d=1;jp3d<=myNbP3d;jp3d++)
		  {
		    //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
		    P3d.SetX(PolesPtr->Value(ipole,index++));
		    //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
		    P3d.SetY(PolesPtr->Value(ipole,index++));
		    //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
		    P3d.SetZ(PolesPtr->Value(ipole,index++)); 
		    TabP3d.SetValue(jp3d,P3d);
		  } 
	      }
	    if(myNbP2d !=0 ) 
	      {
		for (jp2d=1;jp2d<=myNbP2d;jp2d++)
		  {
		    P2d.SetX(PolesPtr->Value(ipole,index++));
		    P2d.SetY(PolesPtr->Value(ipole,index++));  
		    TabP2d.SetValue(jp2d,P2d);
		  }
	      }
	    if(myNbP2d !=0 && myNbP3d !=0) 
	      {
		AppParCurves_MultiPoint aMultiPoint(TabP3d,TabP2d);
		TabMU.SetValue(ipole,aMultiPoint);                     
	      }
	    else if (myNbP2d !=0)
	      {
		AppParCurves_MultiPoint aMultiPoint(TabP2d);
		TabMU.SetValue(ipole,aMultiPoint);                     
	      }
	    else 
	      {
		
		AppParCurves_MultiPoint aMultiPoint(TabP3d);
		TabMU.SetValue(ipole,aMultiPoint);                     
	      }
	    
	  }
	AppParCurves_MultiBSpCurve aCurve(TabMU,myKnots->Array1(),Mults->Array1());
	myMBSpCurve=aCurve; 
	myIsDone = Standard_True;
      }
  }

}
//
//=======================================================================
//function : IsCreated
//purpose  : returns True if the creation is done 
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::IsCreated() const 
{
  return myIsCreated;
}
//
//=======================================================================
//function : IsDone
//purpose  : returns True if the  approximation is ok
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::IsDone() const 
{
  return myIsDone;
}
//
//=======================================================================
//function : IsOverConstrained
//purpose  : returns True if the problem is overconstrained
//           in this case, approximation cannot be done.
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::IsOverConstrained() const 
{
  return myIsOverConstr;
}
//
//=======================================================================
//function : Value
//purpose  : returns all the BSpline curves approximating the
//    	     MultiLine SSP after minimization of the parameter.

//=======================================================================
//
AppParCurves_MultiBSpCurve AppParCurves_Variational::Value() const 
{ 
  if (myIsDone == Standard_False)  StdFail_NotDone::Raise(); 
  return myMBSpCurve;

}
//
//=======================================================================
//function : MaxError
//purpose  : returns the maximum of the distances between 
//           the points of the multiline and the approximation 
//           curves.
//=======================================================================
//
Standard_Real AppParCurves_Variational::MaxError() const 
{
  if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
  return myMaxError;
}
//
//=======================================================================
//function : MaxErrorIndex
//purpose  : returns the index of the MultiPoint of ErrorMax
//=======================================================================
//
Standard_Integer AppParCurves_Variational::MaxErrorIndex() const 
{
  if (myIsDone == Standard_False)  StdFail_NotDone::Raise(); 
  return myMaxErrorIndex;
}
//
//=======================================================================
//function : QuadraticError
//purpose  :  returns the quadratic average of the distances between 
//            the points of the multiline and the approximation 
//            curves.
//=======================================================================
//
Standard_Real AppParCurves_Variational::QuadraticError() const 
{
  if (myIsDone == Standard_False)  StdFail_NotDone::Raise(); 
  return myCriterium[0];
}
//
//=======================================================================
//function : Distance
//purpose  : returns the distances between the points of the 
//           multiline and the approximation curves.
//=======================================================================
//
void AppParCurves_Variational::Distance(math_Matrix& mat)

{
      if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
      Standard_Integer ipoint,jp2d,jp3d,index;
      TColgp_Array1OfPnt TabP3d(1,Max(1,myNbP3d));
      TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
      Standard_Integer j0 = mat.LowerCol() - myFirstPoint;

      gp_Pnt2d P2d;
      gp_Pnt P3d;


      gp_Pnt Pt3d;
      gp_Pnt2d Pt2d;

      for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
        {
          index=1;
          if(myNbP3d !=0 ) {
                            ToolLine::Value(mySSP,ipoint,TabP3d);

                            for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)

		        	      {  P3d=TabP3d.Value(jp3d);
                                         myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt3d);
                                         mat(index++, j0 + ipoint)=P3d.Distance(Pt3d);                                 

			              }
                            }
	  if(myNbP2d !=0 ) {
                            if(myNbP3d == 0 ) ToolLine::Value(mySSP,ipoint,TabP2d);
			    else ToolLine::Value(mySSP,ipoint,TabP3d,TabP2d);
                            for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)

	        		      {  P2d = TabP2d.Value(jp2d);
                                         myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt2d);
                                         mat(index++, j0 + ipoint)=P2d.Distance(Pt2d);
			              }
                           } 
      }

}
//
//=======================================================================
//function : AverageError
//purpose  : returns the average error between          
//           the MultiLine and the approximation.
//=======================================================================
//
Standard_Real AppParCurves_Variational::AverageError() const 
{
  if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
  return myAverageError;
}
//
//=======================================================================
//function : Parameters
//purpose  : returns the parameters uses to the approximations
//=======================================================================
//
const Handle(TColStd_HArray1OfReal)& AppParCurves_Variational::Parameters() const 
{
  if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
  return myParameters;
}
//
//=======================================================================
//function : Knots
//purpose  : returns the knots uses to the approximations
//=======================================================================
//
const Handle(TColStd_HArray1OfReal)& AppParCurves_Variational::Knots() const 
{
  if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
  return myKnots;
}
//
//=======================================================================
//function : Criterium
//purpose  : returns the values of the quality criterium.
//=======================================================================
//
void AppParCurves_Variational::Criterium(Standard_Real& VFirstOrder, Standard_Real& VSecondOrder, Standard_Real& VThirdOrder) const 
{     
  if (myIsDone == Standard_False)  StdFail_NotDone::Raise();
  VFirstOrder=myCriterium[1] ;
  VSecondOrder=myCriterium[2];
  VThirdOrder=myCriterium[3];
}
//
//=======================================================================
//function : CriteriumWeight
//purpose  : returns the Weights (as percent) associed  to the criterium used in
//           the  optimization.
//=======================================================================
//
void AppParCurves_Variational::CriteriumWeight(Standard_Real& Percent1, Standard_Real& Percent2, Standard_Real& Percent3) const 
{
  Percent1 = myPercent[0];
  Percent2 = myPercent[1];
  Percent3 = myPercent[2] ;
}
//
//=======================================================================
//function : MaxDegree
//purpose  : returns the Maximum Degree used in the approximation 
//=======================================================================
//
Standard_Integer AppParCurves_Variational::MaxDegree() const 
{
   return myMaxDegree;
}
//
//=======================================================================
//function : MaxSegment
//purpose  : returns the Maximum of segment used in the approximation
//=======================================================================
//
Standard_Integer AppParCurves_Variational::MaxSegment() const 
{
    return myMaxSegment;
 }
//
//=======================================================================
//function : Continuity
//purpose  : returns the Continuity used in the approximation
//=======================================================================
//
GeomAbs_Shape AppParCurves_Variational::Continuity() const 
{
    return myContinuity;
}
//
//=======================================================================
//function : WithMinMax
//purpose  : returns if the  approximation  search to  minimize the
//           maximum Error or not.
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::WithMinMax() const 
{
    return myWithMinMax;
}
//
//=======================================================================
//function : WithCutting
//purpose  :  returns if the  approximation can insert new Knots or not.
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::WithCutting() const 
{
   return myWithCutting;
}
//
//=======================================================================
//function : Tolerance
//purpose  : returns the tolerance used in the approximation.
//=======================================================================
//
Standard_Real AppParCurves_Variational::Tolerance() const 
{
   return myTolerance;
}
//
//=======================================================================
//function : NbIterations
//purpose  : returns the number of iterations used in the approximation.
//=======================================================================
//
Standard_Integer AppParCurves_Variational::NbIterations() const 
{
   return myNbIterations;
}
//
//=======================================================================
//function : Dump
//purpose  : Prints on the stream o information on the current state 
//           of the object.
//=======================================================================
//
void AppParCurves_Variational::Dump(Standard_OStream& o) const 
{
  o << " \nVariational Smoothing " << endl;
  o << " Number of multipoints                   "  << myNbPoints << endl;
  o << " Number of 2d par multipoint "  << myNbP2d << endl;
  o << " Nombre of 3d par multipoint "  << myNbP3d << endl;
  o << " Number of PassagePoint      "  << myNbPassPoints << endl;
  o << " Number of TangencyPoints    "  << myNbTangPoints << endl;
  o << " Number of CurvaturePoints   "  << myNbCurvPoints << endl;
  o << " \nTolerance " << o.setf(ios::scientific) << setprecision(3) << setw(9) << myTolerance;
  if ( WithMinMax()) { o << "  as Max Error." << endl;}
  else { o << "  as size Error." << endl;}
  o << "CriteriumWeights : " << myPercent[0] << " , " 
    << myPercent[1] << " , " << myPercent[2] << endl;

  if (myIsDone ) {
    o << " MaxError             " << setprecision(3) << setw(9) << myMaxError << endl;
    o << " Index of  MaxError   " << myMaxErrorIndex << endl;
    o << " Average Error        " << setprecision(3) << setw(9) << myAverageError << endl;
    o << " Quadratic Error      " << setprecision(3) << setw(9) << myCriterium[0] << endl;
    o << " Tension Criterium    " << setprecision(3) << setw(9) << myCriterium[1] << endl;
    o << " Flexion  Criterium   " << setprecision(3) << setw(9) << myCriterium[2] << endl;
    o << " Jerk  Criterium      " << setprecision(3) << setw(9) << myCriterium[3] << endl;
    o << " NbSegments           "  << myKnots->Length()-1 << endl;
  }
  else
   { if (myIsOverConstr) o << "The probleme is overconstraint " << endl;
     else o << " Erreur dans l''approximation" << endl;
   }   
}
//
//=======================================================================
//function : SetConstraints
//purpose  : Define the constraints to approximate
//           If this value is incompatible with the others fields
//           this method modify nothing and returns false 
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::SetConstraints(const Handle(AppParCurves_HArray1OfConstraintCouple)& aConstraint)

{
   myConstraints=aConstraint;
   Init();
   if (myIsOverConstr ) return Standard_False;
   else return Standard_True;
}
//
//=======================================================================
//function : SetParameters
//purpose  : Defines the parameters used by the approximations.
//=======================================================================
//
void AppParCurves_Variational::SetParameters(const Handle(TColStd_HArray1OfReal)& param)
{
  myParameters->ChangeArray1() = param->Array1();
}
//
//=======================================================================
//function : SetKnots
//purpose  : Defines the knots used by the approximations
//    --          If this value is incompatible with the others fields
//    --          this method modify nothing and returns false
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::SetKnots(const Handle(TColStd_HArray1OfReal)& knots)
{
  myKnots->ChangeArray1() = knots->Array1();
  return Standard_True;
}
//
//=======================================================================
//function : SetMaxDegree
//purpose  : Define the Maximum Degree used in the approximation
//           If this value is incompatible with the others fields
//           this method modify nothing and returns false 
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::SetMaxDegree(const Standard_Integer Degree)
{
  if (((Degree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
       return Standard_False; 
  else
     {
          myMaxDegree=Degree;

	  InitSmoothCriterion();

          return Standard_True;
	}

}
//
//=======================================================================
//function : SetMaxSegment
//purpose  : Define the maximum number of segments used in the approximation
//           If this value is incompatible with the others fields
//           this method modify nothing and returns false 
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::SetMaxSegment(const Standard_Integer NbSegment)
{
  if ( myWithCutting == Standard_True && 
      ((myMaxDegree-myNivCont)*NbSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
       return Standard_False;
  else
    {
       myMaxSegment=NbSegment;
       return Standard_True;
    }   
}
//
//=======================================================================
//function : SetContinuity
//purpose  : Define the Continuity used in the approximation 
//           If this value is incompatible with the others fields
//           this method modify nothing and returns false 
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::SetContinuity(const GeomAbs_Shape C)
{
    Standard_Integer NivCont=0;
    switch (C) {
    case GeomAbs_C0:
       NivCont=0;
       break ;
    case GeomAbs_C1: 
       NivCont=1;
       break ;
    case GeomAbs_C2:  
       NivCont=2;
       break ;
    default:
    Standard_ConstructionError::Raise();
  }
  if (((myMaxDegree-NivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
       return Standard_False; 
  else
    {
         myContinuity= C;
         myNivCont=NivCont;

	 InitSmoothCriterion();
       return Standard_True;
    }
}
//
//=======================================================================
//function : SetWithMinMax
//purpose  : Define if the  approximation  search to  minimize the
//           maximum Error or not.
//=======================================================================
//
void AppParCurves_Variational::SetWithMinMax(const Standard_Boolean MinMax)
{
  myWithMinMax=MinMax;

  InitSmoothCriterion();
}
//
//=======================================================================
//function : SetWithCutting
//purpose  : Define if the  approximation can insert new Knots or not.
//           If this value is incompatible with the others fields
//           this method modify nothing and returns false 
//=======================================================================
//
Standard_Boolean AppParCurves_Variational::SetWithCutting(const Standard_Boolean Cutting)
{
  if (Cutting == Standard_False) 
      {
       if (((myMaxDegree-myNivCont)*myKnots->Length()-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
       return Standard_False;

       else
          {
        myWithCutting=Cutting;
	InitSmoothCriterion();
        return Standard_True;
      }
      }
  else
    {
      if (((myMaxDegree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
       return Standard_False;

       else
      {
        myWithCutting=Cutting;
	InitSmoothCriterion();
        return Standard_True;
      }
    }
}
//
//=======================================================================
//function : SetCriteriumWeight
//purpose  : define the Weights (as percent) associed to the criterium used in
//           the  optimization.
//=======================================================================
//
void AppParCurves_Variational::SetCriteriumWeight(const Standard_Real Percent1, const Standard_Real Percent2, const Standard_Real Percent3)
{
  if (Percent1 < 0 || Percent2 < 0 || Percent3 < 0 ) Standard_DomainError::Raise();
  Standard_Real Total = Percent1 + Percent2 + Percent3;
  myPercent[0] = Percent1/Total;
  myPercent[1] = Percent2/Total;
  myPercent[2] = Percent3/Total;

  InitSmoothCriterion();
}
//
//=======================================================================
//function : SetCriteriumWeight
//purpose  :  define the  Weight   (as  percent)  associed  to   the
//            criterium   Order used in   the optimization  : Others
//            weights are updated.
//=======================================================================
//
void AppParCurves_Variational::SetCriteriumWeight(const Standard_Integer Order, const Standard_Real Percent)
{
  if ( Percent < 0 ) Standard_DomainError::Raise();
  if ( Order < 1 || Order > 3 ) Standard_ConstructionError::Raise();
  myPercent[Order-1] = Percent;
  Standard_Real Total = myPercent[0] + myPercent[1] + myPercent[2];
  myPercent[0] = myPercent[0] / Total;
  myPercent[1] = myPercent[1] / Total;
  myPercent[2] = myPercent[2] / Total;

  InitSmoothCriterion();
  
}
//
//=======================================================================
//function : SetTolerance
//purpose  : define the tolerance used in the approximation.
//=======================================================================
//
void AppParCurves_Variational::SetTolerance(const Standard_Real Tol)
{
  myTolerance=Tol;
  InitSmoothCriterion();
}
//
//=======================================================================
//function : SetNbIterations
//purpose  : define the number of iterations used in the approximation.
//=======================================================================
//
void AppParCurves_Variational::SetNbIterations(const Standard_Integer Iter)
{
  myNbIterations=Iter;
}


// Private methods
#include <AppParCurves_Variational_1.gxx> // TheMotor
#include <AppParCurves_Variational_2.gxx> // Optimization
#include <AppParCurves_Variational_3.gxx> // Project
#include <AppParCurves_Variational_4.gxx> // ACR
#include <AppParCurves_Variational_5.gxx> // SplitCurve
#include <AppParCurves_Variational_6.gxx> // InitSmoothCriterion and other Init... methods
#include <AppParCurves_Variational_7.gxx> // Adjusting
#include <AppParCurves_Variational_8.gxx> // AssemblingConstraints
#include <AppParCurves_Variational_9.gxx> // InitTtheta and InitTfthet methods