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