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