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