0023178: Intersection of cylinders fails to produce results
[occt.git] / src / IntPatch / IntPatch_ImpImpIntersection_2.gxx
1 // Created on: 1992-05-07
2 // Created by: Jacques GOUSSARD
3 // Copyright (c) 1992-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <IntPatch_WLine.hxx>
18
19 static 
20   Standard_Integer SetQuad(const Handle(Adaptor3d_HSurface)& theS,
21                            GeomAbs_SurfaceType& theTS,
22                            IntSurf_Quadric& theQuad);
23
24 //=======================================================================
25 //function : IntPatch_ImpImpIntersection
26 //purpose  : 
27 //=======================================================================
28 IntPatch_ImpImpIntersection::IntPatch_ImpImpIntersection ():
29         done(Standard_False)
30 {
31 }
32 //=======================================================================
33 //function : IntPatch_ImpImpIntersection
34 //purpose  : 
35 //=======================================================================
36 IntPatch_ImpImpIntersection::IntPatch_ImpImpIntersection
37        (const Handle(Adaptor3d_HSurface)&  S1,
38         const Handle(Adaptor3d_TopolTool)& D1,
39         const Handle(Adaptor3d_HSurface)&  S2,
40         const Handle(Adaptor3d_TopolTool)& D2,
41         const Standard_Real TolArc,
42         const Standard_Real TolTang,
43         const Standard_Boolean theIsReqToKeepRLine)
44 {
45   Perform(S1,D1,S2,D2,TolArc,TolTang, theIsReqToKeepRLine);
46 }
47 //=======================================================================
48 //function : Perform
49 //purpose  : 
50 //=======================================================================
51 void IntPatch_ImpImpIntersection::Perform(const Handle(Adaptor3d_HSurface)&  S1,
52                                           const Handle(Adaptor3d_TopolTool)& D1,
53                                           const Handle(Adaptor3d_HSurface)&  S2,
54                                           const Handle(Adaptor3d_TopolTool)& D2,
55                                           const Standard_Real TolArc,
56                                           const Standard_Real TolTang,
57                                           const Standard_Boolean theIsReqToKeepRLine)
58 {
59   done = Standard_False;
60   spnt.Clear();
61   slin.Clear();
62
63   Standard_Boolean isPostProcessingRequired = Standard_True;
64
65   empt = Standard_True;
66   tgte = Standard_False;
67   oppo = Standard_False;
68
69   Standard_Boolean all1 = Standard_False;
70   Standard_Boolean all2 = Standard_False;
71   Standard_Boolean SameSurf = Standard_False;
72   Standard_Boolean multpoint = Standard_False;
73
74   Standard_Boolean nosolonS1 = Standard_False;
75   // indique s il y a des points sur restriction du carreau 1
76   Standard_Boolean nosolonS2 = Standard_False;
77   // indique s il y a des points sur restriction du carreau 2
78   Standard_Integer i, nbpt, nbseg;
79   IntPatch_SequenceOfSegmentOfTheSOnBounds edg1,edg2;
80   IntPatch_SequenceOfPathPointOfTheSOnBounds pnt1,pnt2;
81   //
82   // On commence par intersecter les supports des surfaces
83   IntSurf_Quadric quad1, quad2;
84   IntPatch_ArcFunction AFunc;
85   const Standard_Real Tolang = 1.e-8;
86   GeomAbs_SurfaceType typs1, typs2;
87   Standard_Boolean bEmpty = Standard_False;
88   //
89   const Standard_Integer iT1 = SetQuad(S1, typs1, quad1);
90   const Standard_Integer iT2 = SetQuad(S2, typs2, quad2);
91   //
92   if (!iT1 || !iT2) {
93     Standard_ConstructionError::Raise();
94     return;
95   }
96   //
97   const Standard_Boolean bReverse = iT1 > iT2;
98   const Standard_Integer iTT = iT1*10 + iT2;
99   //
100   switch (iTT) {
101     case 11: { // Plane/Plane
102       if (!IntPP(quad1, quad2, Tolang, TolTang, SameSurf, slin)) {
103         return;
104       }
105       break;
106     }
107     //
108     case 12:
109     case 21: { // Plane/Cylinder
110       Standard_Real VMin, VMax, H;
111       //
112       const Handle(Adaptor3d_HSurface)& aSCyl = bReverse ? S2 : S1;
113       VMin = aSCyl->FirstVParameter();
114       VMax = aSCyl->LastVParameter();
115       H = (Precision::IsNegativeInfinite(VMin) || 
116            Precision::IsPositiveInfinite(VMax)) ? 0 : (VMax - VMin);
117       //
118       if (!IntPCy(quad1, quad2, Tolang, TolTang, bReverse, empt, slin, H)) {
119         return;
120       }
121       bEmpty = empt;
122       break;
123     }
124     //
125     case 13:
126     case 31: { // Plane/Cone
127       if (!IntPCo(quad1, quad2, Tolang, TolTang, bReverse, empt, multpoint, slin, spnt)) {
128         return;
129       }
130       bEmpty = empt;
131       break;
132     }
133     //
134     case 14:
135     case 41: { // Plane/Sphere
136       if (!IntPSp(quad1, quad2, Tolang, TolTang, bReverse, empt, slin, spnt)) {
137         return;
138       }
139       bEmpty = empt;
140       break;
141     }
142     //
143     case 15:
144     case 51: { // Plane/Torus
145       if (!IntPTo(quad1, quad2, TolTang, bReverse, empt, slin)) {
146         return;
147       }
148       bEmpty = empt;
149       break;
150     }
151     //
152     case 22:
153       { // Cylinder/Cylinder
154         Standard_Boolean isDONE = Standard_False;
155         Bnd_Box2d aBox1, aBox2;
156
157         const Standard_Real aU1f = S1->FirstUParameter();
158         Standard_Real aU1l = S1->LastUParameter();
159         const Standard_Real aU2f = S2->FirstUParameter();
160         Standard_Real aU2l = S2->LastUParameter();
161
162         const Standard_Real anUperiod = 2.0*M_PI;
163
164         if(aU1l - aU1f > anUperiod)
165           aU1l = aU1f + anUperiod;
166
167         if(aU2l - aU2f > anUperiod)
168           aU2l = aU2f + anUperiod;
169
170         aBox1.Add(gp_Pnt2d(aU1f, S1->FirstVParameter()));
171         aBox1.Add(gp_Pnt2d(aU1l, S1->LastVParameter()));
172         aBox2.Add(gp_Pnt2d(aU2f, S2->FirstVParameter()));
173         aBox2.Add(gp_Pnt2d(aU2l, S2->LastVParameter()));
174
175         // Resolution is too big if the cylinder radius is
176         // too small. Therefore, we shall bounde its value above. 
177         // Here, we use simple constant.
178         const Standard_Real a2DTol = Min(1.0e-4, Min( S1->UResolution(TolTang),
179                                             S2->UResolution(TolTang)));
180
181         //The bigger range the bigger nunber of points in Walking-line (WLine)
182         //we will be able to add and, consequently, we will obtain more
183         //precise intersection line.
184         //Every point of WLine is determined as function from U1-parameter,
185         //where U1 is U-parameter on 1st quadric.
186         //Therefore, we should use quadric with bigger range as 1st parameter
187         //in IntCyCy() function.
188         //On the other hand, there is no point in reversing in case of
189         //analytical intersection (when result is line, ellipse, point...).
190         //This result is independent of the arguments order.
191         Standard_Boolean isReversed = ((aU2l - aU2f) < (aU1l - aU1f));
192
193         if(isReversed)
194         {
195           isDONE = IntCyCy(quad2, quad1, TolTang, a2DTol, aBox2, aBox1,
196                                     Standard_True, empt, SameSurf, multpoint, slin, spnt);
197         }
198         else
199         {
200           isDONE = IntCyCy(quad1, quad2, TolTang, a2DTol, aBox1, aBox2,
201                                     Standard_False, empt, SameSurf, multpoint, slin, spnt);
202         }
203
204         if (!isDONE)
205         {
206           return;
207         }
208
209         bEmpty = empt;
210         if(!slin.IsEmpty())
211         {
212           const Handle(IntPatch_WLine)& aWLine = 
213                                     Handle(IntPatch_WLine)::DownCast(slin.Value(1));
214
215           if(!aWLine.IsNull())
216           {//No geometric solution
217             isPostProcessingRequired = Standard_False;
218           }
219         }
220
221         break;
222       }
223     //
224     case 23:
225     case 32: { // Cylinder/Cone
226       if (!IntCyCo(quad1, quad2, TolTang, bReverse, empt, multpoint, slin, spnt)) {
227         return;
228       }
229       bEmpty = empt;
230       break;
231     }
232     //
233     case 24:
234     case 42: { // Cylinder/Sphere
235       if (!IntCySp(quad1, quad2, TolTang, bReverse, empt, multpoint, slin, spnt)) {
236         return;
237       }
238       bEmpty = empt;
239       break;
240     }
241     //
242     case 25:
243     case 52: { // Cylinder/Torus
244       if (!IntCyTo(quad1, quad2, TolTang, bReverse, empt, slin)) {
245         return;
246       }
247       bEmpty = empt;
248       break;
249     }
250     //
251     case 33: { // Cone/Cone
252       if (!IntCoCo(quad1, quad2, TolTang, empt, SameSurf, multpoint, slin, spnt)) {
253         return;
254       }
255       bEmpty = empt;
256       break;
257     }
258     //
259     case 34:
260     case 43: { // Cone/Sphere
261       if (!IntCoSp(quad1, quad2, TolTang, bReverse, empt, multpoint, slin, spnt)) {
262         return;
263       }
264       bEmpty = empt;
265       break;
266     }
267     //
268     case 35:
269     case 53: { // Cone/Torus
270       if (!IntCoTo(quad1, quad2, TolTang, bReverse, empt, slin)) {
271         return;
272       }
273       break;
274     }
275     //
276     case 44: { // Sphere/Sphere
277       if (!IntSpSp(quad1, quad2, TolTang, empt, SameSurf, slin, spnt)) {
278         return;
279       }
280       bEmpty = empt;
281       break;
282     }
283     //
284     case 45:
285     case 54: { // Sphere/Torus
286       if (!IntSpTo(quad1, quad2, TolTang, bReverse, empt, slin)) {
287         return;
288       }
289       bEmpty = empt;
290       break;
291     }
292     //
293     case 55: { // Torus/Torus
294       if (!IntToTo(quad1, quad2, TolTang, SameSurf, empt, slin)) {
295         return;
296       }
297       bEmpty = empt;
298       break;
299     }
300     //
301     default: {
302       Standard_ConstructionError::Raise();
303       break;
304     }
305   }
306   //
307   if (bEmpty) {
308     done = Standard_True;
309     return;
310   }
311   //
312
313   if(isPostProcessingRequired)
314   {
315     if (!SameSurf) {
316       AFunc.SetQuadric(quad2);
317       AFunc.Set(S1);
318     
319       solrst.Perform(AFunc, D1, TolArc, TolTang);
320       if (!solrst.IsDone()) {
321         return;
322       }
323
324       if (solrst.AllArcSolution() && typs1 == typs2) {
325         all1 = Standard_True;
326       }
327       nbpt = solrst.NbPoints();
328       nbseg= solrst.NbSegments();
329       for (i=1; i<= nbpt; i++) {
330         pnt1.Append(solrst.Point(i));
331       }
332       for (i=1; i<= nbseg; i++) {
333         edg1.Append(solrst.Segment(i));
334       }
335       nosolonS1 = (nbpt == 0) && (nbseg == 0);
336
337       if (nosolonS1 && all1) {  // cas de face sans restrictions
338         all1 = Standard_False;
339       }
340     }//if (!SameSurf) {
341     else {
342       nosolonS1 = Standard_True;
343     }
344
345     if (!SameSurf) {
346       AFunc.SetQuadric(quad1);
347       AFunc.Set(S2);
348
349       solrst.Perform(AFunc, D2, TolArc, TolTang);
350       if (!solrst.IsDone()) {
351         return;
352       }
353     
354       if (solrst.AllArcSolution() && typs1 == typs2) {
355         all2 = Standard_True;
356       }
357       nbpt = solrst.NbPoints();
358       nbseg= solrst.NbSegments();
359       for (i=1; i<= nbpt; i++) {
360         pnt2.Append(solrst.Point(i));
361       }
362     
363       for (i=1; i<= nbseg; i++) {
364         edg2.Append(solrst.Segment(i));
365       }
366       nosolonS2 = (nbpt == 0) && (nbseg == 0);
367
368       if (nosolonS2 && all2) {  // cas de face sans restrictions
369         all2 = Standard_False;
370       }
371     }// if (!SameSurf) {
372     else {
373       nosolonS2 = Standard_True;
374     }
375     //
376     if (SameSurf || (all1 && all2)) {
377       // faces "paralleles" parfaites
378       empt = Standard_False;
379       tgte = Standard_True;
380       slin.Clear();
381       spnt.Clear();
382
383       gp_Pnt Ptreference;
384
385       switch (typs1) {
386       case GeomAbs_Plane:      {
387         Ptreference = (S1->Plane()).Location();
388       }
389         break;
390       case GeomAbs_Cylinder: {
391         Ptreference = ElSLib::Value(0.,0.,S1->Cylinder());
392       }
393         break;
394       case GeomAbs_Sphere:      {
395         Ptreference = ElSLib::Value(M_PI/4.,M_PI/4.,S1->Sphere());
396       }
397         break;
398       case GeomAbs_Cone:      {
399         Ptreference = ElSLib::Value(0.,10.,S1->Cone());
400       }
401         break;
402       case GeomAbs_Torus:      {
403         Ptreference = ElSLib::Value(0.,0.,S1->Torus());
404       }
405         break;
406       default: 
407         break;
408       }
409       //
410       oppo = quad1.Normale(Ptreference).Dot(quad2.Normale(Ptreference)) < 0.0;
411       done = Standard_True;
412       return;
413     }// if (SameSurf || (all1 && all2)) {
414
415     if (!nosolonS1 || !nosolonS2) {
416       empt = Standard_False;
417       // C est la qu il faut commencer a bosser...
418       PutPointsOnLine(S1,S2,pnt1, slin, Standard_True, D1, quad1,quad2,
419                       multpoint,TolArc);
420     
421       PutPointsOnLine(S1,S2,pnt2, slin, Standard_False,D2, quad2,quad1,
422                       multpoint,TolArc);
423
424       if (edg1.Length() != 0) {
425         ProcessSegments(edg1,slin,quad1,quad2,Standard_True,TolArc);
426       }
427
428       if (edg2.Length() != 0) {
429         ProcessSegments(edg2,slin,quad1,quad2,Standard_False,TolArc);
430       }
431
432       if (edg1.Length() !=0 || edg2.Length() !=0) {
433         //      ProcessRLine(slin,S1,S2,TolArc);
434         ProcessRLine(slin,quad1,quad2,TolArc, theIsReqToKeepRLine);
435       }
436     }//if (!nosolonS1 || !nosolonS2) {
437     else {
438       empt = ((slin.Length()==0) && (spnt.Length()==0));
439     }
440   }
441   
442   Standard_Integer  nblin = slin.Length(),
443                     aNbPnt = spnt.Length();
444   //
445   //modified by NIZNHY-PKV Tue Sep 06 10:03:35 2011f
446   if (aNbPnt) {
447     IntPatch_SequenceOfPoint aSIP;
448     //
449     for(i=1; i<=aNbPnt; ++i) {
450       Standard_Real aU1, aV1, aU2, aV2;
451       gp_Pnt2d aP2D;
452       TopAbs_State aState1, aState2;
453       //
454       const IntPatch_Point& aIP=spnt(i);
455       aIP.Parameters(aU1, aV1, aU2, aV2);
456       //
457       aP2D.SetCoord(aU1, aV1);
458       aState1=D1->Classify(aP2D, TolArc);
459       //
460       aP2D.SetCoord(aU2, aV2);
461       aState2=D2->Classify(aP2D, TolArc);
462       //
463       if(aState1!=TopAbs_OUT && aState2!=TopAbs_OUT) {
464         aSIP.Append(aIP);
465       }
466     }
467     //
468     spnt.Clear();
469     //
470     aNbPnt=aSIP.Length();
471     for(i=1; i<=aNbPnt; ++i) {
472       const IntPatch_Point& aIP=aSIP(i);
473       spnt.Append(aIP);
474     }
475     //
476   }//  if (aNbPnt) {
477   //modified by NIZNHY-PKV Tue Sep 06 10:18:20 2011t
478   //
479   for(i=1; i<=nblin; i++) {
480     IntPatch_IType thetype = slin.Value(i)->ArcType();
481     if(  (thetype ==  IntPatch_Ellipse)
482        ||(thetype ==  IntPatch_Circle)
483        ||(thetype ==  IntPatch_Lin)
484        ||(thetype ==  IntPatch_Parabola)
485        ||(thetype ==  IntPatch_Hyperbola)) { 
486       Handle(IntPatch_GLine)& glin = *((Handle(IntPatch_GLine)*)&slin.Value(i));
487     glin->ComputeVertexParameters(TolArc); 
488     }
489     else if(thetype == IntPatch_Analytic) { 
490       Handle(IntPatch_ALine)& aligold = *((Handle(IntPatch_ALine)*)&slin.Value(i));
491       aligold->ComputeVertexParameters(TolArc);
492     }
493     else if(thetype == IntPatch_Restriction) {
494       Handle(IntPatch_RLine)& rlig = *((Handle(IntPatch_RLine)*)&slin.Value(i));
495       rlig->ComputeVertexParameters(TolArc); 
496     }
497   }
498   //
499   //----------------------------------------------------------------
500   //-- On place 2 vertex sur les courbes de GLine qui n en 
501   //-- contiennent pas. 
502   for(i=1; i<=nblin; i++) {
503     gp_Pnt P;
504     IntPatch_Point point;
505     Standard_Real u1,v1,u2,v2; 
506     if(slin.Value(i)->ArcType() == IntPatch_Circle) { 
507       const Handle(IntPatch_GLine)& glin = *((Handle(IntPatch_GLine)*)&slin.Value(i));
508       if(glin->NbVertex() == 0) { 
509         gp_Circ Circ = glin->Circle();
510         P=ElCLib::Value(0.0,Circ);
511         quad1.Parameters(P,u1,v1);
512         quad2.Parameters(P,u2,v2);
513         point.SetValue(P,TolArc,Standard_False);
514         point.SetParameters(u1,v1,u2,v2);
515         point.SetParameter(0.0);
516         glin->AddVertex(point);
517
518         P=ElCLib::Value(0.0,Circ);
519         quad1.Parameters(P,u1,v1);
520         quad2.Parameters(P,u2,v2);
521         point.SetValue(P,TolArc,Standard_False);
522         point.SetParameters(u1,v1,u2,v2);
523         point.SetParameter(M_PI+M_PI);
524         glin->AddVertex(point);
525       }
526     }
527     
528     else if(slin.Value(i)->ArcType() == IntPatch_Ellipse) { 
529       const Handle(IntPatch_GLine)& glin = *((Handle(IntPatch_GLine)*)&slin.Value(i));
530       if(glin->NbVertex() == 0) { 
531         gp_Elips Elips = glin->Ellipse();
532         P=ElCLib::Value(0.0,Elips);
533         quad1.Parameters(P,u1,v1);
534         quad2.Parameters(P,u2,v2);
535         point.SetValue(P,TolArc,Standard_False);
536         point.SetParameters(u1,v1,u2,v2);
537         point.SetParameter(0.0);
538         glin->AddVertex(point);
539
540         P=ElCLib::Value(0.0,Elips);
541         quad1.Parameters(P,u1,v1);
542         quad2.Parameters(P,u2,v2);
543         point.SetValue(P,TolArc,Standard_False);
544         point.SetParameters(u1,v1,u2,v2);
545         point.SetParameter(M_PI+M_PI);
546         glin->AddVertex(point);
547       }
548     }
549   }
550   done = Standard_True;
551 }
552
553 //=======================================================================
554 //function : SetQuad
555 //purpose  : 
556 //=======================================================================
557 Standard_Integer SetQuad(const Handle(Adaptor3d_HSurface)& theS,
558                          GeomAbs_SurfaceType& theTS,
559                          IntSurf_Quadric& theQuad)
560 {
561   theTS = theS->GetType();
562   Standard_Integer iRet = 0;
563   switch (theTS) {
564   case GeomAbs_Plane:
565     theQuad.SetValue(theS->Plane());
566     iRet = 1;
567     break;
568   case GeomAbs_Cylinder:
569     theQuad.SetValue(theS->Cylinder());
570     iRet = 2;
571     break;
572   case GeomAbs_Cone:
573     theQuad.SetValue(theS->Cone());
574     iRet = 3;
575     break;
576   case GeomAbs_Sphere:
577     theQuad.SetValue(theS->Sphere());
578     iRet = 4;
579     break;
580   case GeomAbs_Torus:
581     theQuad.SetValue(theS->Torus());
582     iRet = 5;
583     break;
584   default:
585     break;
586   }
587   //
588   return iRet;
589 }
590