0025418: Debug output to be limited to OCC development environment
[occt.git] / src / ShapeCustom / ShapeCustom_Surface.cxx
1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
2 //
3 // This file is part of Open CASCADE Technology software library.
4 //
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
10 //
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
13
14 //abv 06.01.99 fix of misprint
15 //:p6 abv 26.02.99: make ConvertToPeriodic() return Null if nothing done
16 #include <ShapeCustom_Surface.ixx>
17
18 #include <gp_Ax3.hxx>
19 #include <gp_Pnt.hxx>
20 #include <gp_Vec.hxx>
21 #include <gp_Pln.hxx>
22 #include <gp_Cylinder.hxx>
23
24 #include <ElSLib.hxx>
25 #include <TColgp_Array1OfPnt.hxx>
26 #include <TColStd_Array1OfReal.hxx>
27 #include <TColgp_Array2OfPnt.hxx>
28 #include <TColStd_Array2OfReal.hxx>
29 #include <TColStd_Array1OfInteger.hxx>
30
31 #include <Geom_Curve.hxx>
32 #include <Geom_Plane.hxx>
33 #include <Geom_BSplineSurface.hxx>
34 #include <Geom_BezierSurface.hxx>
35 #include <Geom_SphericalSurface.hxx>
36 #include <Geom_CylindricalSurface.hxx>
37 #include <Geom_ConicalSurface.hxx>
38 #include <Geom_ToroidalSurface.hxx>
39 #include <GeomAdaptor_HSurface.hxx>
40 #include <GeomAdaptor_Surface.hxx>
41 #include <GeomAbs_SurfaceType.hxx>
42
43 #include <ShapeAnalysis_Geom.hxx>
44 #include <ShapeAnalysis_Surface.hxx>
45
46 //=======================================================================
47 //function : ShapeCustom_Surface
48 //purpose  : 
49 //=======================================================================
50
51 ShapeCustom_Surface::ShapeCustom_Surface() : myGap (0)
52 {
53 }
54
55 //=======================================================================
56 //function : ShapeCustom_Surface
57 //purpose  : 
58 //=======================================================================
59
60 ShapeCustom_Surface::ShapeCustom_Surface (const Handle(Geom_Surface)& S)
61      : myGap (0)
62 {
63   Init ( S );
64 }
65
66 //=======================================================================
67 //function : Init
68 //purpose  : 
69 //=======================================================================
70
71 void ShapeCustom_Surface::Init (const Handle(Geom_Surface)& S) 
72 {
73   mySurf = S;
74 }
75
76 //=======================================================================
77 //function : ConvertToAnalytical
78 //purpose  : 
79 //=======================================================================
80
81 Handle(Geom_Surface) ShapeCustom_Surface::ConvertToAnalytical (const Standard_Real tol,
82                                                                const Standard_Boolean substitute) 
83 {
84   Handle(Geom_Surface) newSurf;
85
86   Standard_Integer nUP, nVP, nCP, i, j , UDeg, VDeg;
87   Standard_Real U1, U2, V1, V2, C1, C2, DU, DV, U=0, V=0;
88   Handle(Geom_Curve) iso;
89   Standard_Boolean uClosed = Standard_True;
90
91   // seuls cas traites : BSpline et Bezier
92   Handle(Geom_BSplineSurface) theBSplneS = 
93     Handle(Geom_BSplineSurface)::DownCast(mySurf);
94   if (theBSplneS.IsNull()) {
95     Handle(Geom_BezierSurface) theBezierS = 
96       Handle(Geom_BezierSurface)::DownCast(mySurf);
97     if (!theBezierS.IsNull()) {    // Bezier :
98       nUP  = theBezierS->NbUPoles();   
99       nVP  = theBezierS->NbVPoles();   
100       UDeg = theBezierS->UDegree();
101       VDeg = theBezierS->VDegree();
102     }
103     else return newSurf;           // non reconnu : terminus
104   }
105   else {                           // BSpline :
106     nUP  = theBSplneS->NbUPoles();
107     nVP  = theBSplneS->NbVPoles();
108     UDeg = theBSplneS->UDegree();
109     VDeg = theBSplneS->VDegree();
110   }
111
112
113   mySurf->Bounds(U1, U2, V1, V2);
114 //  mySurf->Bounds(U1, U2, V1, V2);
115   TColgp_Array1OfPnt p1(1, 3), p2(1, 3), p3(1, 3);
116   TColStd_Array1OfReal R(1,3);
117   gp_Pnt origPnt, resPnt;  
118   gp_Vec origD1U, resD1U, resD1V;
119     
120   Standard_Boolean aCySpCo = Standard_False;
121   Standard_Boolean aToroid = Standard_False;
122   Standard_Boolean aPlanar = Standard_False;
123
124   if (nUP == 2 && nVP == 2) {
125     if (UDeg == 1 && VDeg == 1) aPlanar = Standard_True;
126   } else if (mySurf->IsUClosed()) {    // VRAI IsUClosed
127     if (mySurf->IsVClosed())   aToroid = Standard_True;
128     else                        aCySpCo = Standard_True;
129   } else {
130     if(mySurf->IsVClosed()) {    // VRAI IsVClosed
131       aCySpCo = Standard_True;
132       uClosed = Standard_False;
133     }
134   }
135
136   if (aPlanar) {
137 //    NearestPlane ...
138     TColgp_Array1OfPnt Pnts(1,4);
139     Pnts.SetValue(1,mySurf->Value(U1,V1));
140     Pnts.SetValue(2,mySurf->Value(U2,V1));
141     Pnts.SetValue(3,mySurf->Value(U1,V2));
142     Pnts.SetValue(4,mySurf->Value(U2,V2));
143     gp_Pln aPln;// Standard_Real Dmax;
144     Standard_Integer It = ShapeAnalysis_Geom::NearestPlane (Pnts,aPln,myGap/*Dmax*/);
145
146 //  ICI, on fabrique le plan, et zou
147     if (It == 0 || myGap/*Dmax*/ > tol) return newSurf;    //  pas un plan
148
149 //    IL RESTE a verifier l orientation ...
150 //    On regarde sur chaque surface les vecteurs P(U0->U1),P(V0->V1)
151 //    On prend la normale : les deux normales doivent etre dans le meme sens
152 //    Sinon, inverser la normale (pas le Pln entier !) et refaire la Plane
153     newSurf = new Geom_Plane (aPln);
154     gp_Vec uold (Pnts(1),Pnts(2));
155     gp_Vec vold (Pnts(1),Pnts(3));
156     gp_Vec nold = uold.Crossed (vold);
157     gp_Vec unew (newSurf->Value(U1,V1), newSurf->Value(U2,V1));
158     gp_Vec vnew (newSurf->Value(U1,V1), newSurf->Value(U1,V2));
159     gp_Vec nnew = unew.Crossed (vnew);
160     if (nold.Dot (nnew) < 0.0) {
161       gp_Ax3 ax3 = aPln.Position();
162       ax3.ZReverse();
163       ax3.XReverse();
164       aPln = gp_Pln (ax3);
165       newSurf = new Geom_Plane (aPln);
166     }
167
168     if (substitute) {
169       Init (newSurf);
170     }
171     return newSurf;
172
173
174   } else if (aCySpCo) {
175     if (!uClosed) {
176       C1 = U1; C2 = U2;
177       U1 = V1; U2 = V2;
178       V1 = C1; V2 = C2;
179       nCP = nUP; nUP = nVP; nVP = nCP;
180     }
181   
182     for (i=1; i<=3; i++) {
183       if      (i==1) V = V1;
184       else if (i==2) V = V2;
185       else if (i==3) V = 0.5*(V1+V2);
186
187       if(uClosed) iso = mySurf->VIso(V);
188       else        iso = mySurf->UIso(V);
189
190       iso->D0(U1, p1(i)); 
191       iso->D0(0.5*(U1+U2), p2(i));
192       p3(i).SetCoord(0.5*(p1(i).X()+p2(i).X()), 
193                      0.5*(p1(i).Y()+p2(i).Y()), 
194                      0.5*(p1(i).Z()+p2(i).Z()));
195       R(i) = p3(i).Distance(p1(i));
196 //      cout<<"sphere, i="<<i<<" V="<<V<<" R="<<R(i)<<" p1="<<p1(i).X()<<","<<p1(i).Y()<<","<<p1(i).Z()<<" p2="<<p2(i).X()<<","<<p2(i).Y()<<","<<p2(i).Z()<<" p3="<<p3(i).X()<<","<<p3(i).Y()<<","<<p3(i).Z()<<endl;
197     }
198       
199     iso->D1 (0.,origPnt,origD1U);
200     gp_Vec xVec(p3(3), p1(3));
201     gp_Vec aVec(p3(1), p3(2));
202 //      gp_Dir xDir(xVec);  ne sert pas. Null si R3 = 0
203     gp_Dir aDir(aVec);
204     gp_Ax3 aAx3 (p3(1),aDir,xVec);
205     //  CKY  3-FEV-1997 : verification du sens de description
206     //gp_Dir AXY = aAx3.YDirection(); // AXY not used (skl)
207     if (aAx3.YDirection().Dot (origD1U) < 0) {
208 #ifdef OCCT_DEBUG
209       cout<<" Surface Analytique : sens a inverser"<<endl;
210 #endif
211       aAx3.YReverse();    //  mais X reste !
212     }
213       
214     if (nVP > 2) {
215       if ((Abs(R(1)) < tol) && 
216           (Abs(R(2)) < tol) && 
217           (Abs(R(3)) > tol)) {
218 // deja fait      gp_Ax3 aAx3(p3(1), aDir, xVec);
219           //gp_Ax3 aAx3(p3(3), aDir);
220         Handle(Geom_SphericalSurface) anObject = 
221           new Geom_SphericalSurface(aAx3, R(3));
222         if (!uClosed) anObject->UReverse();
223         newSurf = anObject;         
224       }
225     }
226     else if (nVP == 2) {
227
228 // deja fait    gp_Ax3 aAx3(p3(1), aDir, xVec); 
229         //gp_Ax3 aAx3(p3(1), aDir);
230
231       if (Abs(R(2)-R(1)) < tol) {
232         Handle(Geom_CylindricalSurface) anObject = 
233           new Geom_CylindricalSurface(aAx3, R(1));
234         if (!uClosed) anObject->UReverse();
235         newSurf = anObject;
236       }
237       else {
238         gp_Vec aVec2(p1(1), p1(2));
239         Standard_Real angle = aVec.Angle(aVec2);
240         if (R(1) < R(2)) {
241           Handle(Geom_ConicalSurface) anObject = 
242             new Geom_ConicalSurface(aAx3, angle, R(1));
243           //if (!uClosed) anObject->UReverse();
244           anObject->UReverse();
245           newSurf = anObject;
246         }
247         else {
248           aDir.Reverse();
249           gp_Vec anotherXVec(p3(2), p1(2));
250           gp_Dir anotherXDir(anotherXVec);
251           gp_Ax3 anotherAx3(p3(2), aDir, anotherXDir);
252           Handle(Geom_ConicalSurface) anObject = 
253             new Geom_ConicalSurface(anotherAx3, angle, R(2));
254           //if (!uClosed) anObject->UReverse();
255           anObject->UReverse();
256           newSurf = anObject;
257         }
258       }
259     }
260   }
261   else if (aToroid) {
262     // test by iso U and isoV
263     Standard_Boolean  isFound = Standard_False;
264     for (j=1; (j<=2) && !isFound; j++) {
265       if (j==2) {
266         C1 = U1; C2 = U2;
267         U1 = V1; U2 = V2;
268         V1 = C1; V2 = C2;
269       }
270       for (i=1; i<=3; i++) {
271         if      (i==1) U = U1;
272         else if (i==2) U = 0.5*(U1+U2);
273         else if (i==3) U = 0.25*(U1+U2);
274
275         iso = mySurf->UIso(U);
276
277         iso->D0(V1, p1(i)); 
278         iso->D0(0.5*(V1+V2), p2(i));
279         p3(i).SetCoord(0.5*(p1(i).X()+p2(i).X()), 
280                        0.5*(p1(i).Y()+p2(i).Y()), 
281                        0.5*(p1(i).Z()+p2(i).Z()));
282         R(i) = p3(i).Distance(p1(i));
283       }
284       if ((Abs(R(1)-R(2))< tol) && 
285           (Abs(R(1)-R(3))< tol)) {
286         gp_Pnt p10(0.5*(p3(1).X()+p3(2).X()), 
287                    0.5*(p3(1).Y()+p3(2).Y()), 
288                    0.5*(p3(1).Z()+p3(2).Z()));
289         gp_Vec aVec(p10, p3(1));
290         gp_Vec aVec2(p10, p3(3));
291         Standard_Real RR1 = R(1), RR2 = R(2), RR3;
292         aVec ^= aVec2;
293
294         if (aVec.Magnitude() <= gp::Resolution()) aVec.SetCoord(0., 0., 1.);
295
296         gp_Dir aDir(aVec);
297
298         gp_Ax3 aAx3(p10, aDir);
299         RR1 = p10.Distance(p3(1));
300 //          modif empirique (pourtant NON DEMONTREE) : inverser roles RR1,RR2
301 //          CKY, 24-JAN-1997
302         if (RR1 < RR2)  {  RR3 = RR1; RR1 = RR2; RR2 = RR3;  }
303         Handle(Geom_ToroidalSurface) anObject = 
304           new Geom_ToroidalSurface(aAx3, RR1, RR2);
305         if (j==2) anObject->UReverse();
306         anObject->D1 (0.,0.,resPnt,resD1U,resD1V);
307 #ifdef OCCT_DEBUG
308         if (resD1U.Dot(origD1U) < 0 && j != 2)
309           cout<<" Tore a inverser !"<<endl;
310 #endif
311         newSurf = anObject;
312         isFound = Standard_True;
313       }
314     }
315   }
316   if (newSurf.IsNull()) return newSurf;
317   
318   //---------------------------------------------------------------------
319   //                 verification
320   //---------------------------------------------------------------------
321   
322   Handle(GeomAdaptor_HSurface) NHS = new GeomAdaptor_HSurface (newSurf);
323   GeomAdaptor_Surface& SurfAdapt = NHS->ChangeSurface();
324
325   const Standard_Integer NP = 21;
326   Standard_Real S = 0., T = 0.;  // U,V deja fait
327   gp_Pnt P3d, P3d2;
328   Standard_Boolean onSurface = Standard_True;
329
330   Standard_Real dis;  myGap = 0.;
331
332   DU = (U2-U1)/(NP-1);
333   DV = (V2-V1)/(NP-1);
334   for (j=1; (j<=NP) && onSurface; j++) {
335     V = V1 + DV*(j-1);
336
337     if(uClosed) iso = mySurf->VIso(V);
338     else        iso = mySurf->UIso(V);
339
340     for (i=1; i<=NP; i++) {
341       U = U1 + DU*(i-1);
342       iso->D0(U, P3d);
343       switch (SurfAdapt.GetType()){    
344           
345         case GeomAbs_Cylinder :
346           {
347             gp_Cylinder Cylinder = SurfAdapt.Cylinder();
348             ElSLib::Parameters( Cylinder, P3d, S, T);
349             break;
350           }
351         case GeomAbs_Cone :
352           {
353             gp_Cone Cone = SurfAdapt.Cone();
354             ElSLib::Parameters( Cone, P3d, S, T);
355             break;
356           }
357         case GeomAbs_Sphere :
358           {
359             gp_Sphere Sphere = SurfAdapt.Sphere();
360             ElSLib::Parameters( Sphere, P3d, S, T);
361             break;
362           }
363         case GeomAbs_Torus :
364           {
365             gp_Torus Torus = SurfAdapt.Torus();
366             ElSLib::Parameters( Torus, P3d, S, T);
367             break;
368           }
369         default:
370           break;
371       }
372
373       newSurf->D0(S, T, P3d2);
374
375       dis = P3d.Distance(P3d2);
376       if (dis > myGap) myGap = dis;
377
378       if (dis > tol) {
379         onSurface = Standard_False;
380         newSurf.Nullify();
381           // The presumption is rejected
382         break;
383       }
384     }
385   }
386   if (substitute && !NHS.IsNull()) {
387     Init (newSurf);
388   }
389   return newSurf;
390 }
391  
392 //%pdn 30 Nov 98: converting bspline surfaces with degree+1 at ends to periodic
393 // UKI60591, entity 48720
394 Handle(Geom_Surface) ShapeCustom_Surface::ConvertToPeriodic (const Standard_Boolean substitute,
395                                                              const Standard_Real preci)
396 {
397   Handle(Geom_Surface) newSurf;
398   Handle(Geom_BSplineSurface) BSpl = Handle(Geom_BSplineSurface)::DownCast(mySurf);
399   if (BSpl.IsNull()) return newSurf;
400   
401   ShapeAnalysis_Surface sas(mySurf);
402   Standard_Boolean uclosed = sas.IsUClosed(preci);
403   Standard_Boolean vclosed = sas.IsVClosed(preci);
404   
405   if ( ! uclosed && ! vclosed ) return newSurf;
406   
407   Standard_Boolean converted = Standard_False; //:p6
408
409   if ( uclosed && ! BSpl->IsUPeriodic() && BSpl->NbUPoles() >3 ) {
410     Standard_Boolean set = Standard_True;
411     // if degree+1 at ends, first change it to 1 by rearranging knots
412     if ( BSpl->UMultiplicity(1) == BSpl->UDegree() + 1 &&
413          BSpl->UMultiplicity(BSpl->NbUKnots()) == BSpl->UDegree() + 1 ) {
414       Standard_Integer nbUPoles = BSpl->NbUPoles();
415       Standard_Integer nbVPoles = BSpl->NbVPoles();
416       TColgp_Array2OfPnt oldPoles(1,nbUPoles,1,nbVPoles);
417       TColStd_Array2OfReal oldWeights(1,nbUPoles,1,nbVPoles);
418       Standard_Integer nbUKnots = BSpl->NbUKnots();
419       Standard_Integer nbVKnots = BSpl->NbVKnots();
420       TColStd_Array1OfReal oldUKnots(1,nbUKnots);
421       TColStd_Array1OfReal oldVKnots(1,nbVKnots);
422       TColStd_Array1OfInteger oldUMults(1,nbUKnots);
423       TColStd_Array1OfInteger oldVMults(1,nbVKnots);
424       
425       BSpl->Poles(oldPoles);
426       BSpl->Weights(oldWeights);
427       BSpl->UKnots(oldUKnots);
428       BSpl->VKnots(oldVKnots);
429       BSpl->UMultiplicities(oldUMults);
430       BSpl->VMultiplicities(oldVMults);
431       
432       TColStd_Array1OfReal newUKnots (1,nbUKnots+2);
433       TColStd_Array1OfInteger newUMults(1,nbUKnots+2);
434       Standard_Real a = 0.5 * ( BSpl->UKnot(2) - BSpl->UKnot(1) + 
435                                 BSpl->UKnot(nbUKnots) - BSpl->UKnot(nbUKnots-1) );
436       
437       newUKnots(1) = oldUKnots(1) - a;
438       newUKnots(nbUKnots+2) = oldUKnots(nbUKnots) + a;
439       newUMults(1) = newUMults(nbUKnots+2) = 1;
440       for (Standard_Integer i = 2; i<=nbUKnots+1; i++) {
441         newUKnots(i) = oldUKnots(i-1);
442         newUMults(i) = oldUMults(i-1);
443       }
444       newUMults(2) = newUMults(nbUKnots+1) = BSpl->UDegree();
445       Handle(Geom_BSplineSurface) res = new Geom_BSplineSurface(oldPoles,
446                                                                 oldWeights,
447                                                                 newUKnots,oldVKnots,
448                                                                 newUMults,oldVMults,
449                                                                 BSpl->UDegree(),BSpl->VDegree(),
450                                                                 BSpl->IsUPeriodic(),BSpl->IsVPeriodic());
451       BSpl = res;
452     }
453     else if ( BSpl->UMultiplicity(1) > BSpl->UDegree() ||
454               BSpl->UMultiplicity(BSpl->NbUKnots()) > BSpl->UDegree() + 1 ) set = Standard_False;
455     if ( set ) {
456       BSpl->SetUPeriodic(); // make periodic
457       converted = Standard_True;
458     }
459   }
460   
461   if ( vclosed && ! BSpl->IsVPeriodic() && BSpl->NbVPoles() >3 ) {      
462     Standard_Boolean set = Standard_True;
463     // if degree+1 at ends, first change it to 1 by rearranging knots
464     if ( BSpl->VMultiplicity(1) == BSpl->VDegree() + 1 &&
465          BSpl->VMultiplicity(BSpl->NbVKnots()) == BSpl->VDegree() + 1 ) {
466       Standard_Integer nbUPoles = BSpl->NbUPoles();
467       Standard_Integer nbVPoles = BSpl->NbVPoles();
468       TColgp_Array2OfPnt oldPoles(1,nbUPoles,1,nbVPoles);
469       TColStd_Array2OfReal oldWeights(1,nbUPoles,1,nbVPoles);
470       Standard_Integer nbUKnots = BSpl->NbUKnots();
471       Standard_Integer nbVKnots = BSpl->NbVKnots();
472       TColStd_Array1OfReal oldUKnots(1,nbUKnots);
473       TColStd_Array1OfReal oldVKnots(1,nbVKnots);
474       TColStd_Array1OfInteger oldUMults(1,nbUKnots);
475       TColStd_Array1OfInteger oldVMults(1,nbVKnots);
476       
477       BSpl->Poles(oldPoles);
478       BSpl->Weights(oldWeights);
479       BSpl->UKnots(oldUKnots);
480       BSpl->VKnots(oldVKnots);
481       BSpl->UMultiplicities(oldUMults);
482       BSpl->VMultiplicities(oldVMults);
483       
484       TColStd_Array1OfReal newVKnots (1,nbVKnots+2);
485       TColStd_Array1OfInteger newVMults(1,nbVKnots+2);
486       Standard_Real a = 0.5 * ( BSpl->VKnot(2) - BSpl->VKnot(1) + 
487                                 BSpl->VKnot(nbVKnots) - BSpl->VKnot(nbVKnots-1) );
488       
489       newVKnots(1) = oldVKnots(1) - a;
490       newVKnots(nbVKnots+2) = oldVKnots(nbVKnots) + a;
491       newVMults(1) = newVMults(nbVKnots+2) = 1;
492       for (Standard_Integer i = 2; i<=nbVKnots+1; i++) {
493         newVKnots(i) = oldVKnots(i-1);
494         newVMults(i) = oldVMults(i-1);
495       }
496       newVMults(2) = newVMults(nbVKnots+1) = BSpl->VDegree();
497       Handle(Geom_BSplineSurface) res = new Geom_BSplineSurface(oldPoles,
498                                                                 oldWeights,
499                                                                 oldUKnots,newVKnots,
500                                                                 oldUMults,newVMults,
501                                                                 BSpl->UDegree(),BSpl->VDegree(),
502                                                                 BSpl->IsUPeriodic(),BSpl->IsVPeriodic());
503       BSpl = res;
504     }
505     else if ( BSpl->VMultiplicity(1) > BSpl->VDegree() ||
506               BSpl->VMultiplicity(BSpl->NbVKnots()) > BSpl->VDegree() + 1 ) set = Standard_False;
507     if ( set ) {
508       BSpl->SetVPeriodic(); // make periodic
509       converted = Standard_True;
510     }
511   }
512   
513 #ifdef OCCT_DEBUG
514   cout << "Warning: ShapeCustom_Surface: Closed BSplineSurface is caused to be periodic" << endl;
515 #endif
516   if ( ! converted ) return newSurf;
517   newSurf = BSpl;
518   if ( substitute ) mySurf = newSurf;
519   return newSurf;
520 }