[occt.git] / src / AppParCurves / AppParCurves_Projection.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 11/09/91


#define No_Standard_RangeError
#define No_Standard_OutOfRange

#include <AppParCurves_Constraint.hxx>
#include <StdFail_NotDone.hxx>
#include <AppParCurves_MultiPoint.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 <PLib.hxx>
#include <BSplCLib.hxx>


AppParCurves_Projection::
   AppParCurves_Projection(const MultiLine& SSP,
         const Standard_Integer FirstPoint,
         const Standard_Integer LastPoint,
	 const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
         math_Vector& Parameters,
         const Standard_Integer Deg,
	 const Standard_Real Tol3d,
	 const Standard_Real Tol2d,
	 const Standard_Integer NbIterations):
	 ParError(FirstPoint, LastPoint,0.0) {

  Standard_Boolean grad = Standard_True;
  Standard_Integer i, j, k, i2, l;
  Standard_Real UF, DU, Fval = 0.0, FU, DFU;
  Standard_Real MErr3d=0.0, MErr2d=0.0, 
                Gain3d = Tol3d, Gain2d=Tol2d;
  Standard_Real EC, ECmax = 1.e10, Errsov3d, Errsov2d;
  Standard_Integer nbP3d = ToolLine::NbP3d(SSP);
  Standard_Integer nbP2d = ToolLine::NbP2d(SSP);
  Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
  Standard_Integer nbP = nbP3d + nbP2d;
  gp_Pnt Pt, P1, P2;
  gp_Pnt2d Pt2d, P12d, P22d;
  gp_Vec V1, V2, MyV;
  gp_Vec2d V12d, V22d, MyV2d;
  
  if (nbP3d == 0) mynbP3d = 1;
  if (nbP2d == 0) mynbP2d = 1;
  TColgp_Array1OfPnt TabP(1, mynbP3d);
  TColgp_Array1OfPnt2d TabP2d(1, mynbP2d);
  TColgp_Array1OfVec TabV(1, mynbP3d);
  TColgp_Array1OfVec2d TabV2d(1, mynbP2d);

  // Calcul de la fonction F= somme(||C(ui)-Ptli||2):
  // Appel a une fonction heritant de MultipleVarFunctionWithGradient
  // pour calculer F et grad_F.
  // ================================================================

  AppParCurves_ProFunction MyF(SSP, FirstPoint,LastPoint, TheConstraints, Parameters, Deg);


  ECmax = 0.0;
  // Projection:
  // ===========
  MyF.Value(Parameters, Fval);
  SCU = MyF.CurveValue();
  Standard_Integer deg = SCU.Degree();
  TColgp_Array1OfPnt   TabPole(1, deg+1),   TabCoef(1, deg+1);
  TColgp_Array1OfPnt2d TabPole2d(1, deg+1), TabCoef2d(1, deg+1);
  TColgp_Array1OfPnt    TheCoef(1, (deg+1)*mynbP3d);
  TColgp_Array1OfPnt2d  TheCoef2d(1, (deg+1)*mynbP2d);

  
  for (i = 1; i <= NbIterations; i++) {

    // Stockage des Poles des courbes:
    // ===============================
    i2 = 0;
    for (k = 1; k <= nbP3d; k++) {
      SCU.Curve(k, TabPole);
      BSplCLib::PolesCoefficients(TabPole, PLib::NoWeights(),
				  TabCoef, PLib::NoWeights());
      for (j=1; j<=deg+1; j++) TheCoef(j+i2) = TabCoef(j);
      i2 += deg+1;
    }
    i2 = 0;
    for (k = 1; k <= nbP2d; k++) {
      SCU.Curve(nbP3d+k, TabPole2d);
      BSplCLib::PolesCoefficients(TabPole2d, PLib::NoWeights(), 
				  TabCoef2d, PLib::NoWeights());
      for (j=1; j<=deg+1; j++) TheCoef2d(j+i2) = TabCoef2d(j);
      i2 += deg+1;
    }
    for (j = FirstPoint+1; j <= LastPoint-1; j++) {
      UF = Parameters(j);
      if (nbP != 0 && nbP2d != 0) ToolLine::Value(SSP, j, TabP, TabP2d);
      else if (nbP2d != 0)        ToolLine::Value(SSP, j, TabP2d);
      else                        ToolLine::Value(SSP, j, TabP);
      
      FU = 0.0;
      DFU = 0.0;
      i2 = 0;
      for (k = 1; k <= nbP3d; k++) {
	for (l=1; l<=deg+1; l++) TabCoef(l) = TheCoef(l+i2);
	i2 += deg+1;
	BSplCLib::CoefsD2(UF, TabCoef, PLib::NoWeights(), Pt, V1, V2);
	MyV = gp_Vec(Pt, TabP(k));
	FU += MyV*V1;
	DFU += V1.SquareMagnitude() + MyV*V2;
      }
      i2 = 0;
      for (k = 1; k <= nbP2d; k++) {
	for (l=1; l<=deg+1; l++) TabCoef2d(l) = TheCoef2d(l+i2);
	i2 += deg+1;
	BSplCLib::CoefsD2(UF, TabCoef2d, PLib::NoWeights(), Pt2d, V12d, V22d);
	MyV2d = gp_Vec2d(Pt2d, TabP2d(k));
	FU += MyV2d*V12d;
	DFU += V12d.SquareMagnitude() + MyV2d*V22d;
      }

      if (DFU <= RealEpsilon()) {
	// Verification du parametrage:
	EC = Abs(Parameters(j) - UF);
	if (EC > ECmax) ECmax = EC;
	break;
      }

      DU = -FU/DFU;
      DU = Sign(Min(5.e-02, Abs(DU)), DU);
      UF += DU;
      Parameters(j) = UF;
    }

    // Test de progression:
    // ====================
    Errsov3d = MErr3d;
    Errsov2d = MErr2d;

    MyF.Value(Parameters, Fval);
    SCU = MyF.CurveValue();
    MErr3d = MyF.MaxError3d();
    MErr2d = MyF.MaxError2d();

    if (MErr3d< Tol3d && MErr2d < Tol2d) {
      Done = Standard_True;
      break;
    }
    if (MErr3d > 100.0*Tol3d || MErr2d > 100.0*Tol2d) {
      Done = Standard_False;
      grad = Standard_False;
      break;
    }
    if (i >= 2) {
      Gain3d = Abs(Errsov3d-MErr3d);
      Gain2d = Abs(Errsov2d-MErr2d);
      if ((MErr3d-Gain3d*(NbIterations-i)*0.5) > Tol3d ||
	  (MErr2d-Gain2d*(NbIterations-i)*0.5) > Tol2d) {
	if (Gain3d < Tol3d/(2*NbIterations) &&
	    Gain2d < Tol2d/(2*NbIterations)) {
	  break;
	}
      }
    }

  }


  AvError = 0.;
  for (j = FirstPoint; j <= LastPoint; j++) {  
    // Recherche des erreurs maxi et moyenne a un index donne:
    for (k = 1; k <= nbP; k++) {
      ParError(j) = Max(ParError(j), MyF.Error(j, k));
    }
    AvError += ParError(j);
  }
  AvError = AvError/(LastPoint-FirstPoint+1);


  MError3d = MyF.MaxError3d();
  MError2d = MyF.MaxError2d();
  if (MError3d < Tol3d && MError2d < Tol2d) {
    Done = Standard_True;
  }

}


AppParCurves_MultiCurve AppParCurves_Projection::Value() const {
  return SCU;
}


Standard_Boolean AppParCurves_Projection::IsDone() const {
  return Done;
}


Standard_Real AppParCurves_Projection::Error(const Standard_Integer Index) const {
  return ParError(Index);
}

Standard_Real AppParCurves_Projection::AverageError() const {
  return AvError;
}

Standard_Real AppParCurves_Projection::MaxError3d() const {
  return MError3d;
}

Standard_Real AppParCurves_Projection::MaxError2d() const {
  return MError2d;
}
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 11/09/91
	16
	17
	18
	19	#define No_Standard_RangeError
	20	#define No_Standard_OutOfRange
	21
	22	#include <AppParCurves_Constraint.hxx>
	23	#include <StdFail_NotDone.hxx>
	24	#include <AppParCurves_MultiPoint.hxx>
	25	#include <gp_Pnt.hxx>
	26	#include <gp_Pnt2d.hxx>
	27	#include <gp_Vec.hxx>
	28	#include <gp_Vec2d.hxx>
	29	#include <TColgp_Array1OfPnt.hxx>
	30	#include <TColgp_Array1OfPnt2d.hxx>
	31	#include <TColgp_Array1OfVec.hxx>
	32	#include <TColgp_Array1OfVec2d.hxx>
	33	#include <PLib.hxx>
	34	#include <BSplCLib.hxx>
	35
	36
	37
	38
	39
	40	AppParCurves_Projection::
	41	AppParCurves_Projection(const MultiLine& SSP,
	42	const Standard_Integer FirstPoint,
	43	const Standard_Integer LastPoint,
	44	const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
	45	math_Vector& Parameters,
	46	const Standard_Integer Deg,
	47	const Standard_Real Tol3d,
	48	const Standard_Real Tol2d,
	49	const Standard_Integer NbIterations):
	50	ParError(FirstPoint, LastPoint,0.0) {
	51
	52	Standard_Boolean grad = Standard_True;
	53	Standard_Integer i, j, k, i2, l;
	54	Standard_Real UF, DU, Fval = 0.0, FU, DFU;
	55	Standard_Real MErr3d=0.0, MErr2d=0.0,
	56	Gain3d = Tol3d, Gain2d=Tol2d;
	57	Standard_Real EC, ECmax = 1.e10, Errsov3d, Errsov2d;
	58	Standard_Integer nbP3d = ToolLine::NbP3d(SSP);
	59	Standard_Integer nbP2d = ToolLine::NbP2d(SSP);
	60	Standard_Integer mynbP3d=nbP3d, mynbP2d=nbP2d;
	61	Standard_Integer nbP = nbP3d + nbP2d;
	62	gp_Pnt Pt, P1, P2;
	63	gp_Pnt2d Pt2d, P12d, P22d;
	64	gp_Vec V1, V2, MyV;
	65	gp_Vec2d V12d, V22d, MyV2d;
	66
	67	if (nbP3d == 0) mynbP3d = 1;
	68	if (nbP2d == 0) mynbP2d = 1;
	69	TColgp_Array1OfPnt TabP(1, mynbP3d);
	70	TColgp_Array1OfPnt2d TabP2d(1, mynbP2d);
	71	TColgp_Array1OfVec TabV(1, mynbP3d);
	72	TColgp_Array1OfVec2d TabV2d(1, mynbP2d);
	73
	74	// Calcul de la fonction F= somme(\|\|C(ui)-Ptli\|\|2):
	75	// Appel a une fonction heritant de MultipleVarFunctionWithGradient
	76	// pour calculer F et grad_F.
	77	// ================================================================
	78
79	AppParCurves_ProFunction MyF(SSP, FirstPoint,LastPoint, TheConstraints, Parameters, Deg);
80
81
82	ECmax = 0.0;
83	// Projection:
84	// ===========
85	MyF.Value(Parameters, Fval);
86	SCU = MyF.CurveValue();
87	Standard_Integer deg = SCU.Degree();
88	TColgp_Array1OfPnt TabPole(1, deg+1), TabCoef(1, deg+1);
89	TColgp_Array1OfPnt2d TabPole2d(1, deg+1), TabCoef2d(1, deg+1);
90	TColgp_Array1OfPnt TheCoef(1, (deg+1)*mynbP3d);
91	TColgp_Array1OfPnt2d TheCoef2d(1, (deg+1)*mynbP2d);
92
93
94	for (i = 1; i <= NbIterations; i++) {
95
96	// Stockage des Poles des courbes:
97	// ===============================
98	i2 = 0;
99	for (k = 1; k <= nbP3d; k++) {
100	SCU.Curve(k, TabPole);
101	BSplCLib::PolesCoefficients(TabPole, PLib::NoWeights(),
102	TabCoef, PLib::NoWeights());
103	for (j=1; j<=deg+1; j++) TheCoef(j+i2) = TabCoef(j);
104	i2 += deg+1;
105	}
106	i2 = 0;
107	for (k = 1; k <= nbP2d; k++) {
108	SCU.Curve(nbP3d+k, TabPole2d);
109	BSplCLib::PolesCoefficients(TabPole2d, PLib::NoWeights(),
110	TabCoef2d, PLib::NoWeights());
111	for (j=1; j<=deg+1; j++) TheCoef2d(j+i2) = TabCoef2d(j);
112	i2 += deg+1;
113	}
114	for (j = FirstPoint+1; j <= LastPoint-1; j++) {
115	UF = Parameters(j);
116	if (nbP != 0 && nbP2d != 0) ToolLine::Value(SSP, j, TabP, TabP2d);
117	else if (nbP2d != 0) ToolLine::Value(SSP, j, TabP2d);
118	else ToolLine::Value(SSP, j, TabP);
119
120	FU = 0.0;
121	DFU = 0.0;
122	i2 = 0;
123	for (k = 1; k <= nbP3d; k++) {
124	for (l=1; l<=deg+1; l++) TabCoef(l) = TheCoef(l+i2);
125	i2 += deg+1;
126	BSplCLib::CoefsD2(UF, TabCoef, PLib::NoWeights(), Pt, V1, V2);
127	MyV = gp_Vec(Pt, TabP(k));
128	FU += MyV*V1;
129	DFU += V1.SquareMagnitude() + MyV*V2;
130	}
131	i2 = 0;
132	for (k = 1; k <= nbP2d; k++) {
133	for (l=1; l<=deg+1; l++) TabCoef2d(l) = TheCoef2d(l+i2);
134	i2 += deg+1;
135	BSplCLib::CoefsD2(UF, TabCoef2d, PLib::NoWeights(), Pt2d, V12d, V22d);
136	MyV2d = gp_Vec2d(Pt2d, TabP2d(k));
137	FU += MyV2d*V12d;
138	DFU += V12d.SquareMagnitude() + MyV2d*V22d;
139	}
140
141	if (DFU <= RealEpsilon()) {
142	// Verification du parametrage:
143	EC = Abs(Parameters(j) - UF);
144	if (EC > ECmax) ECmax = EC;
145	break;
146	}
147
148	DU = -FU/DFU;
149	DU = Sign(Min(5.e-02, Abs(DU)), DU);
150	UF += DU;
151	Parameters(j) = UF;
152	}
153
154	// Test de progression:
155	// ====================
156	Errsov3d = MErr3d;
157	Errsov2d = MErr2d;
158
159	MyF.Value(Parameters, Fval);
160	SCU = MyF.CurveValue();
161	MErr3d = MyF.MaxError3d();
162	MErr2d = MyF.MaxError2d();
163
164	if (MErr3d< Tol3d && MErr2d < Tol2d) {
165	Done = Standard_True;
166	break;
167	}
168	if (MErr3d > 100.0Tol3d \|\| MErr2d > 100.0Tol2d) {
169	Done = Standard_False;
170	grad = Standard_False;
171	break;
172	}
173	if (i >= 2) {
174	Gain3d = Abs(Errsov3d-MErr3d);
175	Gain2d = Abs(Errsov2d-MErr2d);
176	if ((MErr3d-Gain3d(NbIterations-i)0.5) > Tol3d \|\|
177	(MErr2d-Gain2d(NbIterations-i)0.5) > Tol2d) {
178	if (Gain3d < Tol3d/(2*NbIterations) &&
179	Gain2d < Tol2d/(2*NbIterations)) {
180	break;
181	}
182	}
183	}
184
185	}
186
187
188
189	AvError = 0.;
190	for (j = FirstPoint; j <= LastPoint; j++) {
191	// Recherche des erreurs maxi et moyenne a un index donne:
192	for (k = 1; k <= nbP; k++) {
193	ParError(j) = Max(ParError(j), MyF.Error(j, k));
194	}
195	AvError += ParError(j);
196	}
197	AvError = AvError/(LastPoint-FirstPoint+1);
198
199
200	MError3d = MyF.MaxError3d();
201	MError2d = MyF.MaxError2d();
202	if (MError3d < Tol3d && MError2d < Tol2d) {
203	Done = Standard_True;
204	}
205
206	}
207
208
209
210	AppParCurves_MultiCurve AppParCurves_Projection::Value() const {
211	return SCU;
212	}
213
214
215	Standard_Boolean AppParCurves_Projection::IsDone() const {
216	return Done;
217	}
218
219
220	Standard_Real AppParCurves_Projection::Error(const Standard_Integer Index) const {
221	return ParError(Index);
222	}
223
224	Standard_Real AppParCurves_Projection::AverageError() const {
225	return AvError;
226	}
227
228	Standard_Real AppParCurves_Projection::MaxError3d() const {
229	return MError3d;
230	}
231
232	Standard_Real AppParCurves_Projection::MaxError2d() const {
233	return MError2d;
234	}
235
236