5ae4555e9ac70a13799229f929df35e6f94fb529
[occt.git] / src / AppBlend / AppBlend_AppSurf.gxx
1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
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.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
15 #include <AppDef_MultiLine.hxx>
16 #include <AppDef_MultiPointConstraint.hxx>
17 #include <AppParCurves_MultiBSpCurve.hxx>
18 #include <AppParCurves_MultiCurve.hxx>
19 #include <AppDef_BSplineCompute.hxx>
20 #include <AppDef_Compute.hxx>
21 #include <AppParCurves_Constraint.hxx>
22 #include <Approx_MCurvesToBSpCurve.hxx>
23 #include <TColgp_Array1OfPnt.hxx>
24 #include <TColgp_Array1OfPnt2d.hxx>
25 #include <TColgp_Array1OfVec.hxx>
26 #include <TColgp_Array1OfVec2d.hxx>
27 #include <gp_Vec.hxx>
28 #include <gp_Vec2d.hxx>
29 #include <gp_Pnt.hxx>
30 #include <gp_Pnt2d.hxx>
31 #include <math_Vector.hxx>
32 #include <BSplCLib.hxx>
33
34 #include <StdFail_NotDone.hxx>
35 #include <AppParCurves_HArray1OfConstraintCouple.hxx>
36 #include <AppDef_TheVariational.hxx>
37
38 static   Standard_Boolean scal = 1;
39 #ifdef DEB
40 extern Standard_Boolean AppBlend_GetContextSplineApprox(); 
41 extern Standard_Boolean AppBlend_GetContextApproxWithNoTgt(); 
42 #endif
43
44 //  modified by EAP (Edward AGAPOV) Fri Jan 4 2002, bug OCC9
45 //  --- keep pipe parametrized like path
46
47
48 //=======================================================================
49 //function : AppBlend_AppSurf
50 //purpose  : 
51 //=======================================================================
52
53 AppBlend_AppSurf::AppBlend_AppSurf ():done(Standard_False) {}
54
55
56 //=======================================================================
57 //function : AppBlend_AppSurf
58 //purpose  : 
59 //=======================================================================
60
61 AppBlend_AppSurf::AppBlend_AppSurf (const Standard_Integer Degmin,
62                                     const Standard_Integer Degmax,
63                                     const Standard_Real Tol3d,
64                                     const Standard_Real Tol2d,
65                                     const Standard_Integer NbIt,
66                                     const Standard_Boolean KnownParameters):
67        done(Standard_False),dmin(Degmin),dmax(Degmax),
68        tol3d(Tol3d),tol2d(Tol2d),nbit(NbIt),knownp(KnownParameters)
69 {
70   continuity = GeomAbs_C2;
71   paramtype = Approx_ChordLength;
72   critweights[0]=0.4;
73   critweights[1]=0.2;
74   critweights[2]=0.4;
75 }
76
77 //=======================================================================
78 //function : Init
79 //purpose  : 
80 //=======================================================================
81
82 void AppBlend_AppSurf::Init (const Standard_Integer Degmin,
83                              const Standard_Integer Degmax,
84                              const Standard_Real Tol3d,
85                              const Standard_Real Tol2d,
86                              const Standard_Integer NbIt,
87                              const Standard_Boolean KnownParameters)
88 {
89   done  = Standard_False;
90   dmin  = Degmin;
91   dmax  = Degmax;
92   tol3d = Tol3d;
93   tol2d = Tol2d;
94   nbit  = NbIt;
95   knownp = KnownParameters;
96   continuity = GeomAbs_C2;
97   paramtype = Approx_ChordLength;
98   critweights[0]=0.4;
99   critweights[1]=0.2;
100   critweights[2]=0.4;
101 }
102
103 //=======================================================================
104 //function : CriteriumWeight
105 //purpose  : returns the Weights associed  to the criterium used in
106 //           the  optimization.
107 //=======================================================================
108 //
109 void AppBlend_AppSurf::CriteriumWeight(Standard_Real& W1, Standard_Real& W2, Standard_Real& W3) const 
110 {
111   W1 = critweights[0];
112   W2 = critweights[1];
113   W3 = critweights[2] ;
114 }
115 //=======================================================================
116 //function : SetCriteriumWeight
117 //purpose  : 
118 //=======================================================================
119
120 void AppBlend_AppSurf::SetCriteriumWeight(const Standard_Real W1, const Standard_Real W2, const Standard_Real W3)
121 {
122   if (W1 < 0 || W2 < 0 || W3 < 0 ) Standard_DomainError::Raise();
123   critweights[0] = W1;
124   critweights[1] = W2;
125   critweights[2] = W3;
126 }
127 //=======================================================================
128 //function : SetContinuity
129 //purpose  : 
130 //=======================================================================
131
132 void AppBlend_AppSurf::SetContinuity (const GeomAbs_Shape TheCont)
133 {
134   continuity = TheCont;
135 }
136
137 //=======================================================================
138 //function : Continuity
139 //purpose  : 
140 //=======================================================================
141
142 GeomAbs_Shape AppBlend_AppSurf::Continuity () const
143 {
144   return continuity;
145 }
146
147 //=======================================================================
148 //function : SetParType
149 //purpose  : 
150 //=======================================================================
151
152 void AppBlend_AppSurf::SetParType (const Approx_ParametrizationType ParType)
153 {
154   paramtype = ParType;
155 }
156
157 //=======================================================================
158 //function : ParType
159 //purpose  : 
160 //=======================================================================
161
162 Approx_ParametrizationType AppBlend_AppSurf::ParType () const
163 {
164   return paramtype;
165 }
166
167
168 //=======================================================================
169 //function : Perform
170 //purpose  : 
171 //=======================================================================
172
173 void AppBlend_AppSurf::Perform(const Handle(TheLine)& Lin,
174                                TheSectionGenerator& F,
175                                const Standard_Boolean SpApprox)
176
177 {
178   InternalPerform(Lin, F, SpApprox, Standard_False);
179 }
180
181 //=======================================================================
182 //function : PerformSmoothing
183 //purpose  : 
184 //=======================================================================
185
186 void AppBlend_AppSurf::PerformSmoothing(const Handle(TheLine)& Lin,
187                                           TheSectionGenerator& F)
188
189 {
190   InternalPerform(Lin, F, Standard_True, Standard_True);
191 }
192
193 //=======================================================================
194 //function : InternalPerform
195 //purpose  : 
196 //=======================================================================
197
198 void AppBlend_AppSurf::InternalPerform(const Handle(TheLine)& Lin,
199                                        TheSectionGenerator& F,
200                                        const Standard_Boolean SpApprox,
201                                        const Standard_Boolean UseSmoothing)
202
203 {
204   done = Standard_False;
205   if (Lin.IsNull()) {return;}
206   Standard_Integer i,j,k,NbPoint;
207   Standard_Integer NbUPoles,NbUKnots,NbPoles2d,NbVPoles;
208   Standard_Boolean withderiv;
209   AppParCurves_Constraint Cfirst,Clast;
210
211   Standard_Real mytol3d,mytol2d;
212   gp_XYZ newDv;
213
214   seqPoles2d.Clear();
215
216   NbPoint=Lin->NbPoints();
217   AppDef_MultiPointConstraint multP;
218   AppDef_MultiLine multL(NbPoint);
219
220   F.GetShape(NbUPoles,NbUKnots,udeg,NbPoles2d);
221
222   tabUKnots  = new TColStd_HArray1OfReal (1,NbUKnots);
223   tabUMults  = new TColStd_HArray1OfInteger (1,NbUKnots);
224
225   F.Knots(tabUKnots->ChangeArray1());
226   F.Mults(tabUMults->ChangeArray1());
227
228   TColgp_Array1OfPnt tabAppP(1,NbUPoles);
229   TColgp_Array1OfVec tabAppV(1,NbUPoles);
230
231   TColgp_Array1OfPnt2d tabP2d(1,Max(1,NbPoles2d));
232   TColgp_Array1OfVec2d tabV2d(1,Max(1,NbPoles2d));
233
234   TColStd_Array1OfReal tabW(1,NbUPoles),tabDW(1,NbUPoles);
235
236   TColgp_Array1OfPnt2d tabAppP2d(1,NbPoles2d+NbUPoles); // points2d + poids
237   TColgp_Array1OfVec2d tabAppV2d(1,NbPoles2d+NbUPoles); 
238
239
240   AppParCurves_MultiBSpCurve multC;
241
242 //  Standard_Boolean SpApprox = Standard_False;
243
244   withderiv = F.Section(Lin->Point(1),tabAppP,tabAppV,tabP2d,tabV2d,
245                         tabW,tabDW);
246
247 #if DEB
248 //  SpApprox = AppBlend_GetContextSplineApprox();
249   if(AppBlend_GetContextApproxWithNoTgt()) withderiv = Standard_False;
250 #endif
251
252
253   for (j=1; j<=NbPoles2d; j++) {
254     tabAppP2d(j) = tabP2d(j);
255     if (withderiv) {
256       tabAppV2d(j) = tabV2d(j);
257     }
258   }
259   for (j=1; j<=NbUPoles; j++) {
260     // pour les courbes rationnelles il faut multiplier les poles par
261     // leurs poids respectifs
262     if (withderiv) {
263       tabAppV2d(NbPoles2d+j).SetCoord(tabDW(j),0.);
264       newDv.SetLinearForm(tabDW(j),tabAppP(j).XYZ(),tabW(j),tabAppV(j).XYZ());
265       tabAppV(j).SetXYZ(newDv);
266     }
267     tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
268     tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
269   }
270     
271   if (withderiv) {
272     multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d,tabAppV,tabAppV2d);
273     Cfirst = AppParCurves_TangencyPoint;
274   }
275   else {
276     multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
277     Cfirst = AppParCurves_PassPoint;
278   }
279   multL.SetValue(1,multP);
280
281   for (i=2; i<=NbPoint-1; i++) {
282     if (SpApprox) {
283       F.Section(Lin->Point(i),tabAppP,tabP2d,tabW);
284       for (j=1; j<=NbPoles2d; j++) {
285         tabAppP2d(j) = tabP2d(j);
286       }
287       for (j=1; j<=NbUPoles; j++) {
288         // pour les courbes rationnelles il faut multiplier les poles par
289         // leurs poids respectifs
290         tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
291         tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
292       }
293       multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
294       multL.SetValue(i,multP);
295     }
296 // ***********************
297     else {
298       withderiv = F.Section(Lin->Point(i),tabAppP,tabAppV,tabP2d,tabV2d,
299                             tabW,tabDW);
300 #if DEB
301       if(AppBlend_GetContextApproxWithNoTgt()) withderiv = Standard_False;
302 #endif
303       
304       for (j=1; j<=NbPoles2d; j++) {
305         tabAppP2d(j) = tabP2d(j);
306         if (withderiv) {
307           tabAppV2d(j) = tabV2d(j);
308         }
309       }
310       for (j=1; j<=NbUPoles; j++) {
311         // pour les courbes rationnelles il faut multiplier les poles par
312         // leurs poids respectifs
313         if (withderiv) {
314           tabAppV2d(NbPoles2d+j).SetCoord(tabDW(j),0.);
315           newDv.SetLinearForm(tabDW(j),tabAppP(j).XYZ(),tabW(j),tabAppV(j).XYZ());
316           tabAppV(j).SetXYZ(newDv);
317         }
318         tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
319         tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
320       }
321       if (withderiv) {
322         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d,tabAppV,tabAppV2d);
323       }
324       else {
325         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
326       }
327       multL.SetValue(i,multP);
328     }
329 // ******************************
330   }
331   
332   withderiv = F.Section(Lin->Point(NbPoint),tabAppP,tabAppV,tabP2d,tabV2d,
333                         tabW,tabDW);
334 #if DEB
335       if(AppBlend_GetContextApproxWithNoTgt()) withderiv = Standard_False;
336 #endif
337
338   for (j=1; j<=NbPoles2d; j++) {
339     tabAppP2d(j) = tabP2d(j);
340     if (withderiv) {
341       tabAppV2d(j) = tabV2d(j);
342     }
343   }
344   for (j=1; j<=NbUPoles; j++) {
345     // pour les courbes rationnelles il faut multiplier les poles par
346     // leurs poids respectifs
347     if (withderiv) {
348       tabAppV2d(NbPoles2d+j).SetCoord(tabDW(j),0.);
349       newDv.SetLinearForm(tabDW(j),tabAppP(j).XYZ(),tabW(j),tabAppV(j).XYZ());
350       tabAppV(j).SetXYZ(newDv);
351     }
352     tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
353     tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
354   }
355
356   if (withderiv) {
357     multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d,tabAppV,tabAppV2d);
358     Clast = AppParCurves_TangencyPoint;
359   }
360   else {
361     multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
362     Clast = AppParCurves_PassPoint;
363   }
364   multL.SetValue(NbPoint,multP);
365
366   //IFV 04.06.07 occ13904
367   if(NbPoint == 2) {
368     dmin = 1;
369     if(Cfirst == AppParCurves_PassPoint && Clast == AppParCurves_PassPoint) {
370       dmax = 1;
371     }
372   }
373
374
375   if (!SpApprox) {
376     AppDef_Compute theapprox (dmin,dmax,tol3d,tol2d,nbit, Standard_True, paramtype);
377     if (knownp) {
378       math_Vector theParams(1,NbPoint);
379
380       // On recale les parametres entre 0 et 1.
381       theParams(1) = 0.;
382       theParams(NbPoint) = 1.;
383       Standard_Real Uf = F.Parameter(Lin->Point(1));
384       Standard_Real Ul = F.Parameter(Lin->Point(NbPoint))-Uf;
385       for (i=2; i<NbPoint; i++) {
386         theParams(i) = (F.Parameter(Lin->Point(i))-Uf)/Ul;;
387       }
388       AppDef_Compute theAppDef(theParams,dmin,dmax,tol3d,tol2d,nbit,
389                                  Standard_True, Standard_True);
390       theapprox = theAppDef;
391     }
392     theapprox.SetConstraints(Cfirst,Clast);
393     theapprox.Perform(multL);
394
395     Standard_Real TheTol3d, TheTol2d;
396     mytol3d = mytol2d = 0.0;
397     for (Standard_Integer Index=1; Index<=theapprox.NbMultiCurves(); Index++) {
398       theapprox.Error(Index, TheTol3d, TheTol2d);
399       mytol3d = Max(TheTol3d, mytol3d);
400       mytol2d = Max(TheTol2d, mytol2d);
401     }
402 #ifdef DEB
403     cout << " Tolerances obtenues  --> 3d : "<< mytol3d << endl;
404     cout << "                      --> 2d : "<< mytol2d << endl;
405 #endif
406     multC = theapprox.SplineValue();
407   }  
408
409   else {
410     if(!UseSmoothing) {
411       Standard_Boolean UseSquares = Standard_False;
412       if(nbit == 0) UseSquares = Standard_True;
413       AppDef_BSplineCompute theapprox (dmin,dmax,tol3d,tol2d,nbit,Standard_True, paramtype,
414                                        UseSquares);
415       if(continuity == GeomAbs_C0) {
416         theapprox.SetContinuity(0);
417       }
418       if(continuity == GeomAbs_C1) {
419         theapprox.SetContinuity(1);
420       }
421       else if(continuity == GeomAbs_C2) {
422         theapprox.SetContinuity(2);
423       }
424       else {
425         theapprox.SetContinuity(3);
426       }
427
428       theapprox.SetConstraints(Cfirst,Clast);
429
430       if (knownp) {
431         math_Vector theParams(1,NbPoint);
432         // On recale les parametres entre 0 et 1.
433         theParams(1) = 0.;
434         theParams(NbPoint) = 1.;
435         Standard_Real Uf = F.Parameter(Lin->Point(1));
436         Standard_Real Ul = F.Parameter(Lin->Point(NbPoint))-Uf;
437         for (i=2; i<NbPoint; i++) {
438           theParams(i) = (F.Parameter(Lin->Point(i))-Uf)/Ul;;
439         }
440
441         theapprox.Init(dmin,dmax,tol3d,tol2d,nbit,Standard_True,
442                        Approx_IsoParametric,Standard_True);
443         theapprox.SetParameters(theParams);
444       }
445       theapprox.Perform(multL);
446       theapprox.Error(mytol3d,mytol2d);
447 #ifdef DEB
448       cout << " Tolerances obtenues  --> 3d : "<< mytol3d << endl;
449       cout << "                      --> 2d : "<< mytol2d << endl;
450 #endif    
451       tol3dreached = mytol3d;
452       tol2dreached = mytol2d;
453       multC = theapprox.Value();
454     }
455     else {
456       //Variational algo
457       Handle(AppParCurves_HArray1OfConstraintCouple) TABofCC = 
458         new AppParCurves_HArray1OfConstraintCouple(1, NbPoint);
459       AppParCurves_Constraint  Constraint=AppParCurves_NoConstraint;
460
461       for(i = 1; i <= NbPoint; ++i) {
462         AppParCurves_ConstraintCouple ACC(i,Constraint);
463         TABofCC->SetValue(i,ACC);
464       }
465       
466       TABofCC->ChangeValue(1).SetConstraint(Cfirst);
467       TABofCC->ChangeValue(NbPoint).SetConstraint(Clast);
468
469       AppDef_TheVariational Variation(multL, 1, NbPoint, TABofCC);
470
471 //===================================
472       Standard_Integer theMaxSegments = 1000;
473       Standard_Boolean theWithMinMax = Standard_False;
474       Standard_Boolean theWithCutting = Standard_True;
475 //===================================      
476
477       Variation.SetMaxDegree(dmax);
478       Variation.SetContinuity(continuity);
479       Variation.SetMaxSegment(theMaxSegments);
480
481       Variation.SetTolerance(tol3d);
482       Variation.SetWithMinMax(theWithMinMax);
483       Variation.SetWithCutting(theWithCutting);
484       Variation.SetNbIterations(nbit);
485
486       Variation.SetCriteriumWeight(critweights[0], critweights[1], critweights[2]);
487
488       if(!Variation.IsCreated()) {
489         return;
490       }
491   
492       if(Variation.IsOverConstrained()) {
493         return;
494       }
495
496       try {
497         Variation.Approximate();
498       }
499       catch (Standard_Failure) {
500         return;
501       }
502
503       if(!Variation.IsDone()) {
504         return;
505       }
506
507       mytol3d = Variation.MaxError();
508       mytol2d = 0.;
509 #ifdef DEB
510       cout << " Tolerances obtenues  --> 3d : "<< mytol3d << endl;
511       cout << "                      --> 2d : "<< mytol2d << endl;
512 #endif    
513       tol3dreached = mytol3d;
514       tol2dreached = mytol2d;
515       multC = Variation.Value();
516     }
517   }
518
519   vdeg = multC.Degree();
520   NbVPoles = multC.NbPoles();
521   
522   tabPoles   = new TColgp_HArray2OfPnt (1,NbUPoles,1,NbVPoles);
523   tabWeights = new TColStd_HArray2OfReal (1,NbUPoles,1,NbVPoles);
524   tabVKnots  = new TColStd_HArray1OfReal (multC.Knots().Lower(),
525                                           multC.Knots().Upper());
526   tabVKnots->ChangeArray1() = multC.Knots();
527
528   if (knownp && !UseSmoothing) {
529     BSplCLib::Reparametrize(F.Parameter(Lin->Point(1)),
530                             F.Parameter(Lin->Point(NbPoint)),
531                             tabVKnots->ChangeArray1());
532   }
533
534   tabVMults  = new TColStd_HArray1OfInteger (multC.Multiplicities().Lower(),
535                                              multC.Multiplicities().Upper());
536   tabVMults->ChangeArray1() = multC.Multiplicities();
537
538   
539   TColgp_Array1OfPnt newtabP(1,NbVPoles);
540   Handle(TColgp_HArray1OfPnt2d) newtabP2d = 
541     new TColgp_HArray1OfPnt2d(1,NbVPoles);
542   for (j=1; j <=NbUPoles; j++) {
543     multC.Curve(j,newtabP);
544     multC.Curve(j+NbUPoles+NbPoles2d,newtabP2d->ChangeArray1());
545     for (k=1; k<=NbVPoles; k++) {
546       // pour les courbes rationnelles il faut maintenant diviser
547       // les poles par leurs poids respectifs
548       tabPoles->ChangeValue(j,k).
549         SetXYZ(newtabP(k).XYZ()/newtabP2d->Value(k).X());
550       tabWeights->SetValue(j,k,newtabP2d->Value(k).X());
551     }
552   }
553
554   for (j=1; j<=NbPoles2d; j++) {
555     newtabP2d = new TColgp_HArray1OfPnt2d(1,NbVPoles);
556     multC.Curve(NbUPoles+j,newtabP2d->ChangeArray1());
557     seqPoles2d.Append(newtabP2d);
558   }
559   
560   done = Standard_True;
561 }
562
563
564 //=======================================================================
565 //function : Perform
566 //purpose  : 
567 //=======================================================================
568
569 void AppBlend_AppSurf::Perform(const Handle(TheLine)& Lin,
570                                TheSectionGenerator& F,
571                                const Standard_Integer NbMaxP)
572 {
573   done = Standard_False;
574   if (Lin.IsNull()) {return;}
575   Standard_Integer i,j,k;
576   Standard_Integer NbUPoles,NbUKnots,NbPoles2d,NbVPoles;
577   Standard_Boolean withderiv;
578   AppParCurves_Constraint Cfirst=AppParCurves_NoConstraint,Clast=AppParCurves_NoConstraint;
579
580   Standard_Real mytol3d = 0.0, mytol2d = 0.0;
581   gp_XYZ newDv;
582
583   seqPoles2d.Clear();
584
585   Standard_Integer NbPointTot = Lin->NbPoints();
586
587   F.GetShape(NbUPoles,NbUKnots,udeg,NbPoles2d);
588
589   tabUKnots  = new TColStd_HArray1OfReal (1,NbUKnots);
590   tabUMults  = new TColStd_HArray1OfInteger (1,NbUKnots);
591
592   F.Knots(tabUKnots->ChangeArray1());
593   F.Mults(tabUMults->ChangeArray1());
594
595   TColgp_Array1OfPnt tabAppP(1,NbUPoles);
596   TColgp_Array1OfVec tabAppV(1,NbUPoles);
597   Standard_Real X,Y,Z,DX,DY,DZ;
598   X = Y = Z = RealLast();
599   DX = DY = DZ = RealFirst();
600
601   TColgp_Array1OfPnt2d tabP2d(1,Max(1,NbPoles2d));
602   TColgp_Array1OfVec2d tabV2d(1,Max(1,NbPoles2d));
603   TColStd_Array1OfReal X2d(1,Max(1,NbPoles2d));X2d.Init(RealLast());
604   TColStd_Array1OfReal Y2d(1,Max(1,NbPoles2d));Y2d.Init(RealLast());
605   TColStd_Array1OfReal DX2d(1,Max(1,NbPoles2d));DX2d.Init(RealFirst());
606   TColStd_Array1OfReal DY2d(1,Max(1,NbPoles2d));DY2d.Init(RealFirst());
607
608   TColStd_Array1OfReal tabW(1,NbUPoles),tabDW(1,NbUPoles);
609
610   TColgp_Array1OfPnt2d tabAppP2d(1,NbPoles2d+NbUPoles); // points2d + poids
611   TColgp_Array1OfVec2d tabAppV2d(1,NbPoles2d+NbUPoles); 
612
613   // On calcule les boites de chaque ligne (box for all lines)
614   for(i = 1; i <= NbPointTot; i++){
615     F.Section(Lin->Point(i),tabAppP,tabAppV,tabP2d,tabV2d,tabW,tabDW);
616     Standard_Real x,y,z;
617     for(j = 1; j <= NbUPoles; j++){
618       tabAppP(j).Coord(x,y,z);
619       if(x < X) X = x; if(x > DX) DX = x;
620       if(y < Y) Y = y; if(y > DY) DY = y;
621       if(z < Z) Z = z; if(z > DZ) DZ = z;
622     }
623     for(j = 1; j <= NbPoles2d; j++){
624       tabP2d(j).Coord(x,y);
625       if(x < X2d(j)) X2d(j) = x; if(x > DX2d(j)) DX2d(j) = x;
626       if(y < Y2d(j)) Y2d(j) = y; if(y > DY2d(j)) DY2d(j) = y;
627     }
628   }
629   // On calcule pour chaque ligne la transformation vers 0 1.
630   Standard_Real seuil = 1000.*tol3d;
631   Standard_Real seuil2d = 1000.*tol2d;
632   if((DX - X) < seuil ){ DX = 1.; X = 0.; }
633   else{ DX = 1./(DX - X); X *= -DX; }
634   if((DY - Y) < seuil){ DY = 1.; Y = 0.; }
635   else{ DY = 1./(DY - Y); Y *= -DY; }
636   if((DZ - Z) < seuil){ DZ = 1.; Z = 0.; }
637   else{ DZ = 1./(DZ - Z); Z *= -DZ; }
638   for(j = 1; j <= NbPoles2d; j++){
639     if((DX2d(j) - X2d(j)) < seuil2d){ DX2d(j) = 1.; X2d(j) = 0.; }
640     else{ DX2d(j) = 1./(DX2d(j) - X2d(j)); X2d(j) *= -DX2d(j); }
641     if((DY2d(j) - Y2d(j)) < seuil2d){ DY2d(j) = 1.; Y2d(j) = 0.; }
642     else{ DY2d(j) = 1./(DY2d(j) - Y2d(j)); Y2d(j) *= -DY2d(j); }
643   }
644   if(!scal){
645     DX = 1.; X = 0.;
646     DY = 1.; Y = 0.;
647     DZ = 1.; Z = 0.;
648     for(j = 1; j <= NbPoles2d; j++){
649       DX2d(j) = 1.; X2d(j) = 0.;
650       DY2d(j) = 1.; Y2d(j) = 0.;
651     }
652   }
653 //  modified by eap Thu Jan  3 14:45:22 2002 ___BEGIN___
654   // Keep "inter-troncons" parameters, not only first and last
655 //  Standard_Real Ufirst=0,Ulast=0;
656   TColStd_SequenceOfReal aParamSeq;
657    if (knownp) {
658 //     Ufirst = F.Parameter(Lin->Point(1));
659 //     Ulast = F.Parameter(Lin->Point(NbPointTot));
660      aParamSeq.Append( F.Parameter (Lin->Point(1)) );
661   }    
662 //  modified by EAP Thu Jan  3 14:45:41 2002 ___END___
663
664   Approx_MCurvesToBSpCurve concat;
665
666   //On calcule le nombre de troncons.
667   Standard_Integer nbtronc = NbPointTot/NbMaxP;
668   Standard_Integer reste = NbPointTot - (nbtronc * NbMaxP);
669   // On regarde si il faut prendre un troncon de plus.
670   Standard_Integer nmax = NbMaxP;
671   if(nbtronc > 0 && reste > 0){
672     nmax = NbPointTot/(nbtronc + 1);
673     if(nmax > (2*NbMaxP)/3) {
674       nbtronc++;
675       reste = NbPointTot - (nbtronc * nmax);
676     }
677     else nmax = NbMaxP;
678   }
679   else if(nbtronc == 0){
680     nbtronc = 1;
681     nmax = reste;
682     reste = 0;
683   }
684
685   // Approximate each "troncon" with nb of Bezier's using AppDef_Compute
686   // and concat them into BSpline with Approx_MCurvesToBSpCurve 
687
688   TColStd_Array1OfInteger troncsize(1,nbtronc);
689   TColStd_Array1OfInteger troncstart(1,nbtronc);
690
691   Standard_Integer rab = reste/nbtronc + 1;
692   Standard_Integer start = 1;
693   Standard_Integer itronc ;
694   for( itronc = 1; itronc <= nbtronc; itronc++){
695     troncstart(itronc) = start;
696     Standard_Integer rabrab = Min(rab,reste);
697     if(reste > 0){ reste -= rabrab; }
698     troncsize(itronc) = nmax + rabrab + 1;
699     start += (nmax + rabrab);
700   }
701   troncsize(nbtronc) = troncsize(nbtronc) - 1;
702   for(itronc = 1; itronc <= nbtronc; itronc++){
703     Standard_Integer NbPoint = troncsize(itronc); 
704     Standard_Integer StPoint = troncstart(itronc);
705     AppDef_MultiPointConstraint multP;
706     AppDef_MultiLine multL(NbPoint);
707     
708     for (i=1; i<=NbPoint; i++) {
709       Standard_Integer iLin = StPoint + i - 1;
710       Standard_Real x,y,z;
711       withderiv = F.Section(Lin->Point(iLin),tabAppP,tabAppV,tabP2d,tabV2d,
712                             tabW,tabDW);
713 #if DEB
714       if(AppBlend_GetContextApproxWithNoTgt()) withderiv = Standard_False;
715 #endif
716       
717       for (j=1; j<=NbPoles2d; j++) {
718         tabP2d(j).Coord(x,y);
719         tabAppP2d(j).SetCoord(DX2d(j)*x+X2d(j),DY2d(j)*y+Y2d(j));
720         if (withderiv) {
721           tabV2d(j).Coord(x,y);
722           tabAppV2d(j).SetCoord(DX2d(j)*x,DY2d(j)*y);
723         }
724       }
725       for (j=1; j<=NbUPoles; j++) {
726         // pour les courbes rationnelles il faut multiplier les poles par
727         // leurs poids respectifs
728         if (withderiv) {
729           tabAppV2d(NbPoles2d+j).SetCoord(tabDW(j),0.);
730           newDv.SetLinearForm(tabDW(j),tabAppP(j).XYZ(),tabW(j),tabAppV(j).XYZ());
731           tabAppV(j).SetCoord(DX*newDv.X(),DY*newDv.Y(),DZ*newDv.Z());
732         }
733         tabAppP(j).SetXYZ(tabAppP(j).XYZ() * tabW(j));
734         tabAppP2d(NbPoles2d+j).SetCoord(tabW(j),0.);
735         tabAppP(j).Coord(x,y,z);
736         tabAppP(j).SetCoord(DX*x+X,DY*y+Y,DZ*z+Z);
737       }
738       if (withderiv) {
739         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d,tabAppV,tabAppV2d);
740         if(i == 1) Cfirst = AppParCurves_TangencyPoint;
741         else if(i == NbPoint) Clast = AppParCurves_TangencyPoint;
742       }
743       else {
744         multP = AppDef_MultiPointConstraint(tabAppP,tabAppP2d);
745         if(i == 1) Cfirst = AppParCurves_PassPoint;
746         else if(i == NbPoint) Clast = AppParCurves_PassPoint;
747       }
748       multL.SetValue(i,multP);
749     }
750     
751
752   //IFV 04.06.07 occ13904
753     if(NbPoint == 2) {
754       dmin = 1;
755       if(Cfirst == AppParCurves_PassPoint && Clast == AppParCurves_PassPoint) {
756         dmax = 1;
757       }
758     }
759
760 //  modified by EAP Thu Jan  3 15:44:13 2002 ___BEGIN___
761     Standard_Real Ufloc=0., Ulloc=0.;
762     AppDef_Compute theapprox (dmin,dmax,tol3d,tol2d,nbit);
763     if (knownp) {
764       math_Vector theParams(1,NbPoint);
765       // On recale les parametres entre 0 et 1.
766       /*Standard_Real*/ Ufloc = F.Parameter(Lin->Point(StPoint));
767       /*Standard_Real*/ Ulloc = F.Parameter(Lin->Point(StPoint+NbPoint-1));
768 //  modified by EAP Thu Jan  3 15:45:17 2002 ___END___
769       for (i=1; i <= NbPoint; i++) {
770         Standard_Integer iLin = StPoint + i - 1;
771         theParams(i) = (F.Parameter(Lin->Point(iLin))-Ufloc)/(Ulloc - Ufloc);
772       }
773       AppDef_Compute theAppDef1(theParams,dmin,dmax,tol3d,tol2d,nbit, Standard_True,Standard_True);
774       theapprox = theAppDef1;
775     }
776     theapprox.SetConstraints(Cfirst,Clast);
777     theapprox.Perform(multL);
778
779 //  modified by EAP Thu Jan  3 16:00:43 2002 ___BEGIN___
780     // To know internal parameters if multicurve is approximated by several Bezier's
781     TColStd_SequenceOfReal aPoleDistSeq;
782     Standard_Real aWholeDist=0;
783 //  modified by EAP Thu Jan  3 16:45:48 2002 ___END___
784     Standard_Real TheTol3d, TheTol2d;
785     for (Standard_Integer Index=1; Index<=theapprox.NbMultiCurves(); Index++) {
786       AppParCurves_MultiCurve& mucu = theapprox.ChangeValue(Index);
787       theapprox.Error(Index, TheTol3d, TheTol2d);
788       mytol3d = Max(TheTol3d/DX, mytol3d);
789       mytol3d = Max(TheTol3d/DY, mytol3d);
790       mytol3d = Max(TheTol3d/DZ, mytol3d);
791       for(j = 1; j <= NbUPoles; j++){
792         mucu.Transform(j,
793                        -X/DX,1./DX,
794                        -Y/DY,1./DY,
795                        -Z/DZ,1./DZ);
796       }
797       for(j = 1; j <= NbPoles2d; j++){
798         mucu.Transform2d(j + NbUPoles,
799                          -X2d(j)/DX2d(j),1./DX2d(j),
800                          -Y2d(j)/DY2d(j),1./DY2d(j));
801         mytol2d = Max(TheTol2d/DX2d(j), mytol2d);
802         mytol2d = Max(TheTol2d/DY2d(j), mytol2d);
803       }
804       concat.Append(mucu);
805       
806 //  modified by EAP Thu Jan  3 15:45:23 2002 ___BEGIN___
807       if (knownp && theapprox.NbMultiCurves() > 1) 
808       {
809         gp_Pnt aFirstPole = mucu.Pole(Index, 1);
810         gp_Pnt aLastPole  = mucu.Pole(Index, mucu.NbPoles());
811         aPoleDistSeq.Append (aFirstPole.Distance(aLastPole));
812         aWholeDist += aPoleDistSeq.Last();
813       }
814     }
815     if (knownp)
816     {
817       Standard_Integer iDist;
818       Standard_Real iU = Ufloc;
819       for (iDist=1; iDist<aPoleDistSeq.Length(); iDist++)
820       {
821         iU += aPoleDistSeq(iDist) / aWholeDist * (Ulloc - Ufloc);
822         //cout << "Internal: " << iU << endl;
823         aParamSeq.Append(iU);
824       }
825       aParamSeq.Append(Ulloc);
826     }
827 //  modified by EAP Thu Jan  3 15:45:27 2002 ___END___
828   }
829 #ifdef DEB
830   cout << "   Tolerances obtenues  --> 3d : "<< mytol3d << endl;
831   cout << "                        --> 2d : "<< mytol2d << endl;
832 #endif
833   tol3dreached = mytol3d;
834   tol2dreached = mytol2d;
835   concat.Perform();
836   const AppParCurves_MultiBSpCurve& multC = concat.Value();
837   vdeg = multC.Degree();
838   NbVPoles = multC.NbPoles();
839   
840   tabPoles   = new TColgp_HArray2OfPnt (1,NbUPoles,1,NbVPoles);
841   tabWeights = new TColStd_HArray2OfReal (1,NbUPoles,1,NbVPoles);
842   tabVKnots  = new TColStd_HArray1OfReal (multC.Knots().Lower(),
843                                           multC.Knots().Upper());
844   tabVKnots->ChangeArray1() = multC.Knots();
845   
846   if (knownp) {
847 //  modified by EAP Fri Jan  4 12:07:30 2002 ___BEGIN___
848     if (aParamSeq.Length() != tabVKnots->Length())
849     {
850       BSplCLib::Reparametrize(F.Parameter(Lin->Point(1)),
851                               F.Parameter(Lin->Point(Lin->NbPoints())),
852                               tabVKnots->ChangeArray1()
853                               );
854 #ifdef DEB
855       cout << "Warning: AppBlend_AppSurf::Perform(), bad length of aParamSeq: " <<
856         aParamSeq.Length() << " instead of " << tabVKnots->Length() << endl;
857 #endif
858     }
859     else
860     {
861       Standard_Integer iKnot, iTabKnot = tabVKnots->Lower();
862       for (iKnot=1; iKnot<=aParamSeq.Length(); iKnot++, iTabKnot++)
863       {
864         //cout << "Replace " << tabVKnots->Value(iTabKnot) << " with " << aParamSeq(iKnot) << endl;
865         tabVKnots->SetValue(iTabKnot, aParamSeq(iKnot));
866       }
867     }
868 //  modified by EAP Fri Jan  4 12:07:35 2002 ___END___
869   }
870   
871   tabVMults  = new TColStd_HArray1OfInteger (multC.Multiplicities().Lower(),
872                                              multC.Multiplicities().Upper());
873   tabVMults->ChangeArray1() = multC.Multiplicities();
874   
875   
876   TColgp_Array1OfPnt newtabP(1,NbVPoles);
877   Handle(TColgp_HArray1OfPnt2d) newtabP2d = 
878     new TColgp_HArray1OfPnt2d(1,NbVPoles);
879   for (j=1; j <=NbUPoles; j++) {
880     multC.Curve(j,newtabP);
881     multC.Curve(j+NbUPoles+NbPoles2d,newtabP2d->ChangeArray1());
882     for (k=1; k<=NbVPoles; k++) {
883       // pour les courbes rationnelles il faut maintenant diviser
884       // les poles par leurs poids respectifs
885       tabPoles->ChangeValue(j,k).
886         SetXYZ(newtabP(k).XYZ()/newtabP2d->Value(k).X());
887       tabWeights->SetValue(j,k,newtabP2d->Value(k).X());
888     }
889   }
890   
891   for (j=1; j<=NbPoles2d; j++) {
892     newtabP2d = new TColgp_HArray1OfPnt2d(1,NbVPoles);
893     multC.Curve(NbUPoles+j,newtabP2d->ChangeArray1());
894     seqPoles2d.Append(newtabP2d);
895   }
896   
897   done = Standard_True;
898 }
899
900
901 //=======================================================================
902 //function : SurfShape
903 //purpose  : 
904 //=======================================================================
905
906 void AppBlend_AppSurf::SurfShape (Standard_Integer& UDegree,
907                                   Standard_Integer& VDegree,
908                                   Standard_Integer& NbUPoles,
909                                   Standard_Integer& NbVPoles,
910                                   Standard_Integer& NbUKnots,
911                                   Standard_Integer& NbVKnots) const
912 {
913   if (!done) {StdFail_NotDone::Raise();}
914   UDegree  = udeg;
915   VDegree  = vdeg;
916   NbUPoles = tabPoles->ColLength();
917   NbVPoles = tabPoles->RowLength();
918   NbUKnots = tabUKnots->Length();
919   NbVKnots = tabVKnots->Length();
920 }
921
922
923 void AppBlend_AppSurf::Surface(TColgp_Array2OfPnt& TPoles,
924                                TColStd_Array2OfReal& TWeights,
925                                TColStd_Array1OfReal& TUKnots,
926                                TColStd_Array1OfReal& TVKnots,
927                                TColStd_Array1OfInteger& TUMults,
928                                TColStd_Array1OfInteger& TVMults) const
929
930 {
931   if (!done) {StdFail_NotDone::Raise();}
932   TPoles   = tabPoles->Array2();
933   TWeights = tabWeights->Array2();
934   TUKnots  = tabUKnots->Array1();
935   TUMults  = tabUMults->Array1();
936   TVKnots  = tabVKnots->Array1();
937   TVMults  = tabVMults->Array1();
938 }
939
940 //=======================================================================
941 //function : Curves2dShape
942 //purpose  : 
943 //=======================================================================
944
945 void AppBlend_AppSurf::Curves2dShape(Standard_Integer& Degree,
946                                      Standard_Integer& NbPoles,
947                                      Standard_Integer& NbKnots) const
948 {
949   if (!done) {StdFail_NotDone::Raise();}
950   if (seqPoles2d.Length() == 0) {Standard_DomainError::Raise();}
951   Degree = vdeg;
952   NbPoles = tabPoles->ColLength();
953   NbKnots = tabVKnots->Length();
954 }
955
956 //=======================================================================
957 //function : Curve2d
958 //purpose  : 
959 //=======================================================================
960
961 void AppBlend_AppSurf::Curve2d(const Standard_Integer Index,
962                                TColgp_Array1OfPnt2d& TPoles,
963                                TColStd_Array1OfReal& TKnots,
964                                TColStd_Array1OfInteger& TMults) const
965 {
966   if (!done) {StdFail_NotDone::Raise();}
967   if (seqPoles2d.Length() == 0) {Standard_DomainError::Raise();}
968   TPoles = seqPoles2d(Index)->Array1();
969   TKnots  = tabVKnots->Array1();
970   TMults  = tabVMults->Array1();
971 }
972
973 //=======================================================================
974 //function : TolCurveOnSurf
975 //purpose  : 
976 //=======================================================================
977
978 Standard_Real AppBlend_AppSurf::TolCurveOnSurf(const Standard_Integer) const
979 {
980   return tol3dreached; //On ne s'embete pas !!
981 }
982                                         
983
984