0024919: ShapeFix can crash time-to-time due to inproper interator handling inside...
[occt.git] / src / Approx / Approx_SameParameter.cxx
CommitLineData
b311480e 1// Created on: 1995-06-06
2// Created by: Xavier BENVENISTE
3// Copyright (c) 1995-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
17// Modified by skv - Wed Jun 2 11:49:59 2004 OCC5898
18
19#include <Approx_SameParameter.ixx>
20#include <TColStd_Array1OfReal.hxx>
21#include <BSplCLib.hxx>
22#include <Adaptor3d_CurveOnSurface.hxx>
23#include <Geom2dAdaptor_Curve.hxx>
24#include <Geom2dAdaptor_HCurve.hxx>
25#include <GeomAdaptor_Curve.hxx>
26#include <GeomAdaptor_HCurve.hxx>
27#include <GeomAdaptor_Surface.hxx>
28#include <GeomAdaptor_HSurface.hxx>
7fd59977 29#include <GCPnts_QuasiUniformDeflection.hxx>
30#include <Extrema_LocateExtPC.hxx>
31#include <AdvApprox_ApproxAFunction.hxx>
32#include <GeomLib_MakeCurvefromApprox.hxx>
33#include <Precision.hxx>
a86d3ec0 34#include <Extrema_ExtPC.hxx>
7fd59977 35
36#ifdef DEB
37#ifdef DRAW
38#include <DrawTrSurf.hxx>
39#endif
40#include <Geom2d_BSplineCurve.hxx>
41#include <stdio.h>
42static Standard_Boolean Voir = Standard_False;
43static Standard_Boolean AffichFw = Standard_False;
44static Standard_Integer NbCurve = 0;
45#endif
46//
0d969553 47// allows testing if Extrema produces correct results/
7fd59977 48
49
50static void ProjectPointOnCurve(const Standard_Real InitValue,
a86d3ec0 51 const gp_Pnt APoint,
52 const Standard_Real Tolerance,
53 const Standard_Integer NumIteration,
54 const Adaptor3d_Curve& Curve,
55 Standard_Boolean& Status,
56 Standard_Real& Result)
7fd59977 57{
58 Standard_Integer num_iter = 0,
a86d3ec0 59 not_done = 1,
60 ii ;
61
7fd59977 62 gp_Pnt a_point ;
63 gp_Vec vector,
a86d3ec0 64 d1,
65 d2 ;
7fd59977 66 Standard_Real func,
a86d3ec0 67 func_derivative,
68 param = InitValue ;
7fd59977 69 Status = Standard_False ;
70 Standard_Real Toler = 1.0e-12;
71 do {
72 num_iter += 1 ;
73 Curve.D2(param,
a86d3ec0 74 a_point,
75 d1,
76 d2) ;
7fd59977 77 for (ii = 1 ; ii <= 3 ; ii++) {
78 vector.SetCoord(ii, APoint.Coord(ii) - a_point.Coord(ii)) ;
79 }
80 func = vector.Dot(d1) ;
81 func_derivative = vector.Dot(d2) ;
82 func_derivative -= d1.Dot(d1) ;
83 if ( Abs(func) < Tolerance * d1.Magnitude()) {
84 not_done = 0 ;
85 Status = Standard_True ;
86 }
87 else
a86d3ec0 88 { // fixing a bug PRO18577 : avoid divizion by zero
89 if( Abs(func_derivative) > Toler ) {
90 param -= func / func_derivative ;
7fd59977 91 }
a86d3ec0 92 param = Max(param,Curve.FirstParameter()) ;
93 param = Min(param,Curve.LastParameter()) ;
94 //Status = Standard_True ;
95 }
7fd59977 96 }
97 while (not_done && num_iter <= NumIteration) ;
98 Result = param ;
99}
a86d3ec0 100
7fd59977 101
102
103//=======================================================================
104//class : Approx_SameParameter_Evaluator
105//purpose :
106//=======================================================================
107
108class Approx_SameParameter_Evaluator : public AdvApprox_EvaluatorFunction
109{
a86d3ec0 110public:
7fd59977 111 Approx_SameParameter_Evaluator (const TColStd_Array1OfReal& theFlatKnots,
a86d3ec0 112 const TColStd_Array1OfReal& thePoles,
113 const Handle(Adaptor2d_HCurve2d)& theHCurve2d)
7fd59977 114 : FlatKnots(theFlatKnots), Poles(thePoles), HCurve2d(theHCurve2d) {}
115
116 virtual void Evaluate (Standard_Integer *Dimension,
a86d3ec0 117 Standard_Real StartEnd[2],
118 Standard_Real *Parameter,
119 Standard_Integer *DerivativeRequest,
120 Standard_Real *Result, // [Dimension]
121 Standard_Integer *ErrorCode);
122
123private:
7fd59977 124 const TColStd_Array1OfReal& FlatKnots;
125 const TColStd_Array1OfReal& Poles;
126 Handle(Adaptor2d_HCurve2d) HCurve2d;
127};
128
129void Approx_SameParameter_Evaluator::Evaluate (Standard_Integer *,/*Dimension*/
a86d3ec0 130 Standard_Real /*StartEnd*/[2],
131 Standard_Real *Parameter,
132 Standard_Integer *DerivativeRequest,
133 Standard_Real *Result,
134 Standard_Integer *ReturnCode)
7fd59977 135{
136 gp_Pnt2d Point ;
137 gp_Vec2d Vector ;
138 Standard_Integer extrap_mode[2] ;
139 extrap_mode[0] = extrap_mode[1] = 3;
140 Standard_Real eval_result[2] ;
141 Standard_Real *PolesArray =
142 (Standard_Real *) &Poles(Poles.Lower()) ;
143 //
144 // evaluate the 1D bspline that represents the change in parameterization
145 //
146 BSplCLib::Eval(*Parameter,
a86d3ec0 147 Standard_False,
148 *DerivativeRequest,
149 extrap_mode[0],
150 3,
151 FlatKnots,
152 1,
153 PolesArray[0],
154 eval_result[0]) ;
155
156
7fd59977 157 if (*DerivativeRequest == 0){
158 HCurve2d->D0(eval_result[0],Point);
159 Point.Coord(Result[0],Result[1]);
160 }
161 else if (*DerivativeRequest == 1){
162 HCurve2d->D1(eval_result[0], Point, Vector);
163 Vector.Multiply(eval_result[1]);
164 Vector.Coord(Result[0],Result[1]);
165 }
166 ReturnCode[0] = 0 ;
167}
168
169static Standard_Real ComputeTolReached(const Handle(Adaptor3d_HCurve)& c3d,
a86d3ec0 170 const Adaptor3d_CurveOnSurface& cons,
171 const Standard_Integer nbp)
7fd59977 172{
173 Standard_Real d2 = 0.;
a86d3ec0 174 const Standard_Real first = c3d->FirstParameter();
175 const Standard_Real last = c3d->LastParameter();
176 for(Standard_Integer i = 0; i <= nbp; i++){
177 Standard_Real t = IntToReal(i)/IntToReal(nbp);
7fd59977 178 Standard_Real u = first*(1.-t) + last*t;
179 gp_Pnt Pc3d = c3d->Value(u);
180 gp_Pnt Pcons = cons.Value(u);
181 if (Precision::IsInfinite(Pcons.X()) ||
a86d3ec0 182 Precision::IsInfinite(Pcons.Y()) ||
183 Precision::IsInfinite(Pcons.Z())) {
184 d2=Precision::Infinite();
185 break;
7fd59977 186 }
187 Standard_Real temp = Pc3d.SquareDistance(Pcons);
188 if(temp > d2) d2 = temp;
189 }
190 d2 = 1.5*sqrt(d2);
191 if(d2<1.e-7) d2 = 1.e-7;
192 return d2;
193}
194
195static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
a86d3ec0 196 const TColStd_Array1OfReal& Poles,
197 const Standard_Integer nbp,
198 const TColStd_Array1OfReal& pc3d,
199 // const TColStd_Array1OfReal& pcons,
200 const TColStd_Array1OfReal& ,
201 const Handle(Adaptor3d_HCurve)& c3d,
202 const Adaptor3d_CurveOnSurface& cons,
203 Standard_Real& tol,
204 const Standard_Real oldtol)
7fd59977 205{
206 Standard_Real d = tol;
207 Standard_Integer extrap_mode[2] ;
208 extrap_mode[0] = extrap_mode[1] = 3;
209 Standard_Integer i;
210#ifdef DEB
211 if (Voir) {
212 cout<<endl;
0d969553
Y
213 cout<<"Control the change of variable : "<<endl;
214 cout<<"yawn mesured by projection : "<<d<<endl;
215 cout<<"Number of points : "<<nbp<<endl;
7fd59977 216 }
217#endif
218#if 0
219 Standard_Real glis = 0., dglis = 0.;
220 for(i = 1; i <= nbp; i++){
221 Standard_Real tc3d = pc3d(i);
222 gp_Pnt Pc3d = c3d->Value(tc3d);
223 Standard_Real tcons;
224 BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
a86d3ec0 225 3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
7fd59977 226 gp_Pnt Pcons = cons.Value(tcons);
227 Standard_Real temp = Pc3d.SquareDistance(Pcons);
228 if(temp >= dglis) dglis = temp;
229 temp = Abs(tcons-pcons(i));
230 if(temp >= glis) glis = temp;
231 }
232 dglis = sqrt(dglis);
233#ifdef DEB
234 if ( Voir) {
0d969553
Y
235 cout<<"shift of parameter to the imposed points : "<<glis<<endl;
236 cout<<"shift distance at the imposed points : "<<dglis<<endl;
7fd59977 237 }
238#endif
239 dglis = 0.;
240 for(i = 1; i < nbp; i++){
241 Standard_Real tc3d = 0.5*(pc3d(i)+pc3d(i+1));
242 gp_Pnt Pc3d = c3d->Value(tc3d);
243 Standard_Real tcons;
244 BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
a86d3ec0 245 3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
7fd59977 246 gp_Pnt Pcons = cons.Value(tcons);
247 Standard_Real temp = Pc3d.SquareDistance(Pcons);
248 if(temp >= dglis) dglis = temp;
249 }
250 dglis = sqrt(dglis);
251#ifdef DEB
252 if (Voir)
253 cout<<"distance de glissement en milieu d intervals : "<<dglis<<endl;
254#endif
255#endif
256
257 Standard_Real d2 = 0.;
258 Standard_Integer nn = 2*nbp;
259 Standard_Real unsurnn = 1./nn;
a86d3ec0 260 // Modified by skv - Wed Jun 2 11:49:59 2004 OCC5898 Begin
261 // Correction of the interval of valid values. This condition has no sensible
262 // grounds. But it is better then the old one (which is commented out) because
263 // it fixes the bug OCC5898. To develop more or less sensible criterion it is
264 // necessary to deeply investigate this problem which is not possible in frames
265 // of debugging.
266
267 // Standard_Real firstborne= 2*pc3d(1)-pc3d(nbp);
268 // Standard_Real lastborne= 2*pc3d(nbp)-pc3d(1);
7fd59977 269 Standard_Real firstborne= 3.*pc3d(1) - 2.*pc3d(nbp);
270 Standard_Real lastborne = 3.*pc3d(nbp) - 2.*pc3d(1);
a86d3ec0 271 // Modified by skv - Wed Jun 2 11:50:03 2004 OCC5898 End
7fd59977 272 for(i = 0; i <= nn; i++){
273 Standard_Real t = unsurnn*i;
274 Standard_Real tc3d = pc3d(1)*(1.-t) + pc3d(nbp)*t;
275 gp_Pnt Pc3d = c3d->Value(tc3d);
276 Standard_Real tcons;
277 BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
a86d3ec0 278 3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
7fd59977 279 if (tcons < firstborne || tcons > lastborne) {
280 tol=Precision::Infinite();
281 return Standard_False;
282 }
283 gp_Pnt Pcons = cons.Value(tcons);
284 Standard_Real temp = Pc3d.SquareDistance(Pcons);
285 if(temp > d2) d2 = temp;
286 }
287 tol = sqrt(d2);
288#ifdef DEB
289 if (Voir)
0d969553 290 cout<<"distance max on "<<nn<<" points : "<<tol<<endl<<endl;
7fd59977 291#endif
292 return ((tol <= d) || (tol > 0.8 * oldtol));
293}
294
295
296//=======================================================================
297//function : Approx_SameParameter
298//purpose :
299//=======================================================================
300
301Approx_SameParameter::Approx_SameParameter(const Handle(Geom_Curve)& C3D,
a86d3ec0 302 const Handle(Geom2d_Curve)& C2D,
303 const Handle(Geom_Surface)& S,
304 const Standard_Real Tol):
305mySameParameter(Standard_True), myDone(Standard_False)
7fd59977 306{
307 myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
308 myC3d = new GeomAdaptor_HCurve(C3D);
309 mySurf = new GeomAdaptor_HSurface(S);
310 Build(Tol);
311}
312
313
314//=======================================================================
315//function : Approx_SameParameter
316//purpose :
317//=======================================================================
318
319Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)& C3D,
a86d3ec0 320 const Handle(Geom2d_Curve)& C2D,
321 const Handle(Adaptor3d_HSurface)& S,
322 const Standard_Real Tol):
323mySameParameter(Standard_True), myDone(Standard_False)
7fd59977 324{
325 myC3d = C3D;
326 mySurf = S;
327 myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
328 Build(Tol);
329}
330
331
332//=======================================================================
333//function : Approx_SameParameter
334//purpose :
335//=======================================================================
336
337Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)& C3D,
a86d3ec0 338 const Handle(Adaptor2d_HCurve2d)& C2D,
339 const Handle(Adaptor3d_HSurface)& S,
340 const Standard_Real Tol):
341mySameParameter(Standard_True), myDone(Standard_False)
7fd59977 342{
343 myC3d = C3D;
344 mySurf = S;
345 myHCurve2d = C2D;
346 Build(Tol);
347}
348
349
350//=======================================================================
351//function : Build
352//purpose :
353//=======================================================================
7fd59977 354void Approx_SameParameter::Build(const Standard_Real Tolerance)
355{
a86d3ec0 356 const Standard_Real anErrorMAX = 1.0e15;
357 const Standard_Integer aMaxArraySize = 1000;
358 const Standard_Integer NCONTROL = 22;
359
7fd59977 360 Standard_Integer ii ;
361 Adaptor3d_CurveOnSurface CurveOnSurface(myHCurve2d,mySurf);
362 Standard_Real fcons = CurveOnSurface.FirstParameter();
363 Standard_Real lcons = CurveOnSurface.LastParameter();
364 Standard_Real fc3d = myC3d->FirstParameter();
365 Standard_Real lc3d = myC3d->LastParameter();
366
367 GeomAbs_Shape Continuity = myHCurve2d->Continuity();
368
369 if(Continuity > GeomAbs_C1) Continuity = GeomAbs_C1;
370
0d969553
Y
371 //Control tangents at the extremities to know if the
372 //reparametring is possible and calculate the tangents
373 //at the extremities of the function of change of variable.
d20d815b 374 Standard_Real tangent[2] = { 0.0, 0.0 };
7fd59977 375 gp_Pnt Pcons,Pc3d;
376 gp_Vec Vcons,Vc3d;
377
378 Standard_Real Tol = Tolerance;
379 Standard_Real Tol2 = Tol * Tol;
380 Standard_Real Tolp = myC3d->Resolution(Tol), deltamin = 50*Tolp;
381
382 Standard_Real besttol2 = Tol2;
383 Standard_Boolean extrok = 0;
384
385 extrok = 1;
386 CurveOnSurface.D1(fcons,Pcons,Vcons);
387 myC3d->D1(fc3d,Pc3d,Vc3d);
388 Standard_Real dist2 = Pcons.SquareDistance(Pc3d);
389 Standard_Real dmax2 = dist2;
390
391 Standard_Real magVcons = Vcons.Magnitude();
392 if (magVcons > 1.e-12){
393 tangent[0] = Vc3d.Magnitude() / magVcons;
394 }
395 else extrok = 0;
396
397 CurveOnSurface.D1(lcons,Pcons,Vcons);
398 myC3d->D1(lc3d,Pc3d,Vc3d);
399 dist2 = Pcons.SquareDistance(Pc3d);
400
401 if(dist2 > dmax2) dmax2 = dist2;
402 magVcons = Vcons.Magnitude();
403 if (magVcons > 1.e-12){
404 tangent[1] = Vc3d.Magnitude() / magVcons;
405 }
406 else extrok = 0;
407
408
409 if(dmax2 > besttol2) besttol2 = dmax2;
410
0d969553
Y
411 //Take a multiple of the sample pof CheckShape,
412 //at least the control points will be correct. No comment!!!
7fd59977 413
7fd59977 414#ifdef DEB
415 Standard_Integer nbcoups = 0;
416#endif
a86d3ec0 417
7fd59977 418 Standard_Boolean interpolok = 0;
419 Standard_Real tolsov = 1.e200;
0d969553
Y
420 //Take parameters with constant step on the curve on surface
421 //and on curve 3d.
7fd59977 422 Standard_Real deltacons = lcons - fcons;
423 deltacons /= (NCONTROL);
424 Standard_Real deltac3d = lc3d - fc3d;
425 deltac3d /= (NCONTROL);
426
427 Standard_Real wcons = fcons;
428 Standard_Real wc3d = fc3d;
a86d3ec0 429
430 Standard_Real qpcons[aMaxArraySize], qnewpcons[aMaxArraySize],
431 qpc3d[aMaxArraySize], qnewpc3d[aMaxArraySize];
7fd59977 432 Standard_Real * pcons = qpcons; Standard_Real * newpcons = qnewpcons;
433 Standard_Real * pc3d = qpc3d; Standard_Real * newpc3d = qnewpc3d;
434
435 for ( ii = 0 ; ii < NCONTROL; ii++) {
436 pcons[ii] = wcons;
437 pc3d[ii] = wc3d;
438 wcons += deltacons;
439 wc3d += deltac3d;
440 }
441 pcons[NCONTROL] = lcons;
442 pc3d[NCONTROL] = lc3d;
443
444 Standard_Integer New_NCONTROL = NCONTROL;
445 if(Continuity < GeomAbs_C1) {
a86d3ec0 446 Standard_Integer NbInt = myHCurve2d->NbIntervals(GeomAbs_C1) + 1;
447 TColStd_Array1OfReal Param_de_decoupeC1 (1, NbInt);
448 myHCurve2d->Intervals(Param_de_decoupeC1, GeomAbs_C1);
449 TColStd_SequenceOfReal new_par;
450 Standard_Integer inter = 1;
451 ii =1;
452 new_par.Append(fcons);
453
3f5bebe8 454 while(inter <= NbInt && Param_de_decoupeC1(inter) <= fcons + deltamin) inter++;
455 while(NbInt > 0 && Param_de_decoupeC1(NbInt) >= lcons - deltamin) NbInt--;
a86d3ec0 456
3f5bebe8 457 while(inter <= NbInt || (ii < NCONTROL && inter <= Param_de_decoupeC1.Length()) ) {
a86d3ec0 458 if(Param_de_decoupeC1(inter) < pcons[ii]) {
459 new_par.Append(Param_de_decoupeC1(inter));
460 if((pcons[ii] - Param_de_decoupeC1(inter)) <= deltamin) {
461 ii++;
462 if(ii > NCONTROL) {ii = NCONTROL;}
463 }
464 inter++;
465 }
466 else {
467 if((Param_de_decoupeC1(inter) - pcons[ii]) > deltamin) {
468 new_par.Append(pcons[ii]);
469 }
470 ii++;
471 }
472 }
473
474 new_par.Append(lcons);
475 New_NCONTROL = new_par.Length() - 1;
476 //simple protection if New_NCONTROL > allocated elements in array
477 if (New_NCONTROL > aMaxArraySize) {
478 mySameParameter = Standard_False;
479 return;
480 }
481 for(ii = 1; ii <= New_NCONTROL; ii++){
482 pcons[ii] = pc3d[ii] = new_par.Value(ii + 1);
483 }
484 pc3d[New_NCONTROL] = lc3d;
485 }
486
487
7fd59977 488 Extrema_LocateExtPC Projector;
489 Projector.Initialize(myC3d->Curve(),fc3d,lc3d,Tol);
a86d3ec0 490
7fd59977 491 Standard_Integer count = 1;
492 Standard_Real previousp = fc3d, initp=0, curp;//, deltamin = 50*Tolp;
493 Standard_Real bornesup = lc3d - deltamin;
494 Standard_Boolean projok = 0,
495 use_parameter ;
496 for (ii = 1; ii < New_NCONTROL; ii++){
497 CurveOnSurface.D0(pcons[ii],Pcons);
498 myC3d->D0(pc3d[ii],Pc3d);
499 dist2 = Pcons.SquareDistance(Pc3d);
500 use_parameter = (dist2 <= Tol2 && (pc3d[ii] > pc3d[count-1] + deltamin)) ;
501 if(use_parameter) {
a86d3ec0 502
7fd59977 503 if(dist2 > dmax2) dmax2 = dist2;
504 initp = previousp = pc3d[count] = pc3d[ii];
505 pcons[count] = pcons[ii];
506 count++;
507 }
508 else {
509 if(!projok) initp = pc3d[ii];
510 projok = mySameParameter = Standard_False;
511 Projector.Perform(Pcons, initp);
512 if (Projector.IsDone()) {
a86d3ec0 513 curp = Projector.Point().Parameter();
514 Standard_Real dist_2 = Projector.SquareDistance();
515 if(dist_2 > besttol2) besttol2 = dist_2;
516 projok = 1;
7fd59977 517 }
a86d3ec0 518 else
519 {
520 ProjectPointOnCurve(initp,Pcons,Tol,30,myC3d->Curve(),projok,curp);
7fd59977 521 }
a86d3ec0 522
523 if(projok)
524 {
525 if(curp > previousp + deltamin && curp < bornesup){
526 initp = previousp = pc3d[count] = curp;
527 pcons[count] = pcons[ii];
528 count++;
529 }
7fd59977 530 }
a86d3ec0 531 else
532 {
533 Extrema_ExtPC PR(Pcons,myC3d->Curve(),fc3d,lc3d,Tol);
534 if(PR.IsDone())
535 {
536 const Standard_Integer aNbExt = PR.NbExt();
537 if(aNbExt > 0)
538 {
539 Standard_Integer anIndMin = 0;
540 Standard_Real aDistMin = RealLast();
541 for(Standard_Integer i = 1; i <= aNbExt; i++)
542 {
543 const gp_Pnt &aP = PR.Point(i).Value();
544 Standard_Real aDist2 = aP.SquareDistance(Pcons);
545 if(aDist2 < aDistMin)
546 {
547 aDistMin = aDist2;
548 anIndMin = i;
549 }
550 }
551 curp = PR.Point(anIndMin).Parameter();
552 if(curp > previousp + deltamin && curp < bornesup)
553 {
554 initp = previousp = pc3d[count] = curp;
555 pcons[count] = pcons[ii];
556 count++;
557 projok = Standard_True;
558 }
559 }
560 }
561 }
562
563 if(!projok)
564 {
565 //Projector
566#ifdef DEB
567 // JAG
568 cout << "Projection not done" << endl;
7fd59977 569#endif
570 }
571 }
572 }
a86d3ec0 573
7fd59977 574 if(mySameParameter){
575 myTolReached = 1.5*sqrt(dmax2);
576 return;
577 }
a86d3ec0 578
0d969553 579 if(!extrok) { // If not already SameP and tangent to mill, abandon.
7fd59977 580 mySameParameter = Standard_False;
581#ifdef DEB
0d969553 582 cout<<"SameParameter problem : zero tangent to extremities"<<endl;
7fd59977 583#endif
584 return;
585 }
586
587 pcons[count] = lcons;
588 pc3d[count] = lc3d;
589
590#ifdef DEB
591 if (AffichFw) {
592 char Name[17];
593 Name[0]='\0';
594 TColgp_Array1OfPnt2d DEBP2d (0,count);
595 TColStd_Array1OfInteger DEBMults(0,count);
596 DEBMults.Init(1); DEBMults(0) = 2; DEBMults(count) = 2;
597 TColStd_Array1OfReal DEBKnots(0,count);
598 for (Standard_Integer DEBi = 0; DEBi <= count; DEBi++) {
599 DEBKnots(DEBi) = DEBi;
600 DEBP2d (DEBi) = gp_Pnt2d(pc3d[DEBi],pcons[DEBi]);
601 }
602 Handle(Geom2d_BSplineCurve) DEBBS =
603 new Geom2d_BSplineCurve(DEBP2d,DEBKnots,DEBMults,1);
91322f44 604 Sprintf(Name,"DEBC2d_%d",++NbCurve);
7fd59977 605#ifdef DRAW
606 DrawTrSurf::Set(Name,DEBBS);
607#endif
608 }
609#endif
a86d3ec0 610
611 Standard_Boolean hasCountChanged = Standard_False;
612
613 while(!interpolok)
614 {
0d969553 615 // The tables and their limits for the interpolation.
7fd59977 616 Standard_Integer num_knots = count + 7;
617 Standard_Integer num_poles = count + 3;
618 TColStd_Array1OfReal Paramc3d(*pc3d,1,count+1);
619 TColStd_Array1OfReal Paramcons(*pcons,1,count+1);
620 TColStd_Array1OfInteger ContactOrder(1,num_poles) ;
621 TColStd_Array1OfReal Poles(1,num_poles) ;
622 TColStd_Array1OfReal InterpolationParameters(1,num_poles) ;
623 TColStd_Array1OfReal FlatKnots(1,num_knots) ;
a86d3ec0 624
0d969553 625 // Fill tables taking attention to end values.
7fd59977 626 ContactOrder.Init(0);
627 ContactOrder(2) = ContactOrder(num_poles - 1) = 1;
a86d3ec0 628
7fd59977 629 FlatKnots(1) = FlatKnots(2) = FlatKnots(3) = FlatKnots(4) = fc3d;
630 FlatKnots(num_poles + 1) = FlatKnots(num_poles + 2) =
631 FlatKnots(num_poles + 3) = FlatKnots(num_poles + 4) = lc3d;
a86d3ec0 632
7fd59977 633 Poles(1) = fcons; Poles(num_poles) = lcons;
634 Poles(2) = tangent[0]; Poles(num_poles - 1) = tangent[1];
a86d3ec0 635
7fd59977 636 InterpolationParameters(1) = InterpolationParameters(2) = fc3d;
637 InterpolationParameters(num_poles - 1) = InterpolationParameters(num_poles) = lc3d;
a86d3ec0 638
7fd59977 639 for (ii = 3; ii <= num_poles - 2; ii++) {
640 Poles(ii) = Paramcons(ii - 1);
641 InterpolationParameters(ii) = FlatKnots(ii+2) = Paramc3d(ii - 1);
642 }
643 Standard_Integer inversion_problem;
644 BSplCLib::Interpolate(3,FlatKnots,InterpolationParameters,ContactOrder,
a86d3ec0 645 1,Poles(1),inversion_problem);
7fd59977 646 if(inversion_problem) {
647 Standard_ConstructionError::Raise();
648 }
649
650 //-------------------------------------------
651#ifdef DEB
a86d3ec0 652 if (AffichFw) {
653 nbcoups ++;
654 char Name[17];
655 Name[0] = '\0';
656 Standard_Integer nnn = 100;
657 TColgp_Array1OfPnt2d DEBP2d (0,nnn);
658 TColStd_Array1OfInteger DEBMults(0,nnn);
659 DEBMults.Init(1); DEBMults(0) = 2; DEBMults(nnn) = 2;
660 TColStd_Array1OfReal DEBKnots(0,nnn);
661 Standard_Real du = (lc3d - fc3d) / nnn;
662 Standard_Real u3d = fc3d;
663 Standard_Integer extrap_mode[2] ;
664 extrap_mode[0] = extrap_mode[1] = 3;
665 Standard_Real eval_result[2] ;
666 Standard_Integer DerivativeRequest = 0;
667 Standard_Real *PolesArray =
668 (Standard_Real *) &Poles(Poles.Lower()) ;
669
670 for (Standard_Integer DEBi = 0; DEBi <= nnn; DEBi++) {
671 DEBKnots(DEBi) = DEBi;
672 BSplCLib::Eval(u3d,
673 Standard_False,
674 DerivativeRequest,
675 extrap_mode[0],
676 3,
677 FlatKnots,
678 1,
679 PolesArray[0],
680 eval_result[0]) ;
681
682 DEBP2d (DEBi) = gp_Pnt2d(u3d,eval_result[0]);
683 u3d += du;
684 }
7fd59977 685
a86d3ec0 686 Handle(Geom2d_BSplineCurve) DEBBS =
687 new Geom2d_BSplineCurve(DEBP2d,DEBKnots,DEBMults,1);
688 Sprintf(Name,"DEBC2d_%d_%d",NbCurve,nbcoups );
7fd59977 689#ifdef DRAW
a86d3ec0 690 DrawTrSurf::Set(Name,DEBBS);
7fd59977 691#endif
a86d3ec0 692 }
7fd59977 693#endif
a86d3ec0 694 //-------------------------------------------
7fd59977 695
a86d3ec0 696 //-------------------------------------------
7fd59977 697 // Test if par2d(par3d) is monotonous function or not ----- IFV, Jan 2000
698 // and try to insert new point to improve BSpline interpolation
699
700 Standard_Integer extrap_mode[2] ;
701 extrap_mode[0] = extrap_mode[1] = 3;
702 Standard_Real eval_result[2] ;
703 Standard_Integer DerivativeRequest = 0;
704 Standard_Real *PolesArray =
705 (Standard_Real *) &Poles(Poles.Lower()) ;
706
707 Standard_Integer newcount = 0;
708 for (ii = 0; ii < count; ii++) {
a86d3ec0 709
7fd59977 710 newpcons[newcount] = pcons[ii];
711 newpc3d[newcount] = pc3d[ii];
712 newcount++;
713
a86d3ec0 714 if(count - ii + newcount == aMaxArraySize) continue;
7fd59977 715
716 BSplCLib::Eval(.5*(pc3d[ii]+pc3d[ii+1]), Standard_False, DerivativeRequest,
a86d3ec0 717 extrap_mode[0], 3, FlatKnots, 1, PolesArray[0], eval_result[0]);
718
7fd59977 719 if(eval_result[0] < pcons[ii] || eval_result[0] > pcons[ii+1]) {
a86d3ec0 720 Standard_Real ucons = 0.5*(pcons[ii]+pcons[ii+1]);
721 Standard_Real uc3d = 0.5*(pc3d[ii]+pc3d[ii+1]);
722
723 CurveOnSurface.D0(ucons,Pcons);
724 Projector.Perform(Pcons, uc3d);
725 if (Projector.IsDone()) {
726 curp = Projector.Point().Parameter();
727 Standard_Real dist_2 = Projector.SquareDistance();
728 if(dist_2 > besttol2) besttol2 = dist_2;
729 projok = 1;
730 }
731 else {
732 ProjectPointOnCurve(uc3d,Pcons,Tol,30,myC3d->Curve(),projok,curp);
733 }
734 if(projok){
735 if(curp > pc3d[ii] + deltamin && curp < pc3d[ii+1] - deltamin){
736 newpc3d[newcount] = curp;
737 newpcons[newcount] = ucons;
738 newcount ++;
739 }
740 }
741 else {
7fd59977 742#ifdef DEB
a86d3ec0 743 // JAG
744 cout << "Projection not done" << endl;
7fd59977 745#endif
a86d3ec0 746 }
7fd59977 747 }
a86d3ec0 748
7fd59977 749 }
750
751 newpc3d[newcount] = pc3d[count];
752 newpcons[newcount] = pcons[count];
753 Standard_Real * temp;
754 temp = pc3d;
755 pc3d = newpc3d;
756 newpc3d = temp;
757 temp = pcons;
758 pcons = newpcons;
759 newpcons = temp;
760
a86d3ec0 761 if((count != newcount) && newcount < aMaxArraySize) { count = newcount; continue;}
7fd59977 762
763 count = newcount;
764
765 Standard_Real algtol = sqrt(besttol2);
766
767 interpolok = Check (FlatKnots, Poles, count+1, Paramc3d, Paramcons,
a86d3ec0 768 myC3d, CurveOnSurface, algtol, tolsov);
7fd59977 769
770 if (Precision::IsInfinite(algtol)) {
771 mySameParameter = Standard_False;
772#ifdef DEB
0d969553 773 cout<<"SameParameter problem : function of interpolation of parametration at mills !!"<<endl;
7fd59977 774#endif
775 return;
776 }
777
778 tolsov = algtol;
779
a86d3ec0 780 interpolok = (interpolok || count >= aMaxArraySize);
7fd59977 781
782 if(interpolok) {
a86d3ec0 783 Standard_Real besttol = sqrt(besttol2);
7fd59977 784#ifdef DEB
785 if (Voir) {
a86d3ec0 786 if(algtol > besttol){
787 cout<<"SameParameter : Tol can't be reached before approx"<<endl;
788 }
7fd59977 789 }
790#endif
791 Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
792 tol1d = new TColStd_HArray1OfReal(1,2) ;
793 tol1d->SetValue(1, mySurf->UResolution(besttol));
794 tol1d->SetValue(2, mySurf->VResolution(besttol));
795
796 Approx_SameParameter_Evaluator ev (FlatKnots, Poles, myHCurve2d);
797 AdvApprox_ApproxAFunction anApproximator(2,0,0,tol1d,tol2d,tol3d,fc3d,lc3d,
a86d3ec0 798 Continuity,11,40,ev);
7fd59977 799
800 if (anApproximator.IsDone() || anApproximator.HasResult()) {
a86d3ec0 801 Adaptor3d_CurveOnSurface ACS = CurveOnSurface;
802 GeomLib_MakeCurvefromApprox aCurveBuilder(anApproximator) ;
803 Handle(Geom2d_BSplineCurve) aC2d = aCurveBuilder.Curve2dFromTwo1d(1,2) ;
804 Handle(Adaptor2d_HCurve2d) aHCurve2d = new Geom2dAdaptor_HCurve(aC2d);
805 CurveOnSurface.Load(aHCurve2d);
806
807 myTolReached = ComputeTolReached(myC3d,CurveOnSurface,NCONTROL);
808
809 if(myTolReached > anErrorMAX)
810 {
811 //This tolerance is too big. Probably, we will not can get
812 //edge with sameparameter in this case.
813
814 myDone = Standard_False;
815 return;
816 }
817
818 if( (myTolReached < 250.0*besttol) ||
819 (count >= aMaxArraySize-2) ||
820 !hasCountChanged) //if count does not change after adding new point
821 //(else we can have circularity)
822 {
823 myCurve2d = aC2d;
824 myHCurve2d = new Geom2dAdaptor_HCurve(myCurve2d);
825 myDone = Standard_True;
826 }
827 else
828 {
829 interpolok = Standard_False;
830 CurveOnSurface = ACS;
831 }
7fd59977 832 }
833 }
a86d3ec0 834
835 if(!interpolok)
836 {
7fd59977 837#ifdef DEB
838 if (Voir)
a86d3ec0 839 cout<<"SameParameter : Not enough points, enrich"<<endl;
7fd59977 840#endif
841
a86d3ec0 842 newcount = 0;
7fd59977 843 for(Standard_Integer n = 0; n < count; n++){
a86d3ec0 844 newpc3d[newcount] = pc3d[n];
845 newpcons[newcount] = pcons[n];
846 newcount ++;
847
848 if(count - n + newcount == aMaxArraySize) continue;
849
850 Standard_Real ucons = 0.5*(pcons[n]+pcons[n+1]);
851 Standard_Real uc3d = 0.5*(pc3d[n]+pc3d[n+1]);
852
853 CurveOnSurface.D0(ucons,Pcons);
854 Projector.Perform(Pcons, uc3d);
855 if (Projector.IsDone()) {
856 curp = Projector.Point().Parameter();
857 Standard_Real dist_2 = Projector.SquareDistance();
858 if(dist_2 > besttol2) besttol2 = dist_2;
859 projok = 1;
860 }
861 else {
862 ProjectPointOnCurve(uc3d,Pcons,Tol,30,myC3d->Curve(),projok,curp);
863 }
864 if(projok){
865 if(curp > pc3d[n] + deltamin && curp < pc3d[n+1] - deltamin){
866 newpc3d[newcount] = curp;
867 newpcons[newcount] = ucons;
868 newcount ++;
869 }
870 }
871 else {
7fd59977 872#ifdef DEB
a86d3ec0 873 // JAG
874 cout << "Projection not done" << endl;
7fd59977 875#endif
a86d3ec0 876 }
7fd59977 877 }
878 newpc3d[newcount] = pc3d[count];
879 newpcons[newcount] = pcons[count];
880 Standard_Real * tempx;
881 tempx = pc3d;
882 pc3d = newpc3d;
883 newpc3d = tempx;
884 tempx = pcons;
885 pcons = newpcons;
886 newpcons = tempx;
a86d3ec0 887
888 if(count != newcount)
889 {
890 count = newcount;
891 hasCountChanged = Standard_True;
892 }
893 else
894 {
895 hasCountChanged = Standard_False;
896 }
7fd59977 897 }
898 }
899}