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