0025368: BREPExtrma DistShapeShape gives wrong result for Sphere and Line
[occt.git] / src / Extrema / Extrema_ExtElCS.cxx
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 //  Modified by skv - Thu Jul  7 14:37:05 2005 OCC9134
16
17 #include <Extrema_ExtElCS.ixx>
18 #include <Extrema_ExtPElS.hxx>
19 #include <Extrema_ExtPElC.hxx>
20 #include <Extrema_ExtElC.hxx>
21 #include <Extrema_POnCurv.hxx>
22 #include <Standard_NotImplemented.hxx>
23 #include <StdFail_InfiniteSolutions.hxx>
24 #include <Precision.hxx>
25 #include <ElSLib.hxx>
26 #include <ElCLib.hxx>
27 #include <gp_Vec.hxx>
28 #include <IntAna_Quadric.hxx>
29 #include <IntAna_IntConicQuad.hxx>
30
31
32 Extrema_ExtElCS::Extrema_ExtElCS() 
33 {
34   myDone = Standard_False;
35   myIsPar = Standard_False;
36 }
37
38
39 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
40                                  const gp_Pln& S)
41 {
42   Perform(C, S);
43 }
44
45
46
47 void Extrema_ExtElCS::Perform(const gp_Lin& C,
48                               const gp_Pln& S)
49 {
50   myDone = Standard_True;
51   myIsPar = Standard_False;
52
53   if (C.Direction().IsNormal(S.Axis().Direction(), 
54                              Precision::Angular())) {
55     mySqDist = new TColStd_HArray1OfReal(1, 1);
56     mySqDist->SetValue(1, S.SquareDistance(C));
57     myIsPar = Standard_True;
58   }
59   else {
60     myNbExt = 0;
61   }
62   
63 }
64
65
66 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
67                                  const gp_Cylinder& S)
68 {
69   Perform(C, S);
70 }
71
72
73
74 void Extrema_ExtElCS::Perform(const gp_Lin& C,
75                               const gp_Cylinder& S)
76 {
77   myDone = Standard_False;
78   myNbExt = 0;
79   myIsPar = Standard_False;
80
81   gp_Ax3 Pos = S.Position();
82   gp_Pnt Origin = Pos.Location();
83   gp_Pnt LineOrig = C.Location();
84
85   Standard_Real radius = S.Radius();
86   Extrema_ExtElC Extrem(gp_Lin(Pos.Axis()), C, Precision::Angular());
87   if (Extrem.IsParallel()) {
88     // Line direction is similar to cylinder axis of rotation.
89     mySqDist = new TColStd_HArray1OfReal(1, 1);
90     myPoint1 = new Extrema_HArray1OfPOnCurv(1, 1);
91     myPoint2 = new Extrema_HArray1OfPOnSurf(1, 1);
92     Standard_Real aDist = sqrt(Extrem.SquareDistance(1)) - radius;
93     mySqDist->SetValue(1, aDist * aDist);
94     Standard_Real u, v, w;
95     gp_Vec aVec(LineOrig, Origin);
96     gp_Vec aDirVec(C.Direction());
97     w = aVec*aDirVec;
98     gp_Pnt LinPoint = LineOrig.Translated(w * aDirVec);
99     Extrema_POnCurv PonC(w, LinPoint);
100     myPoint1->SetValue(1, PonC);
101     gp_Pnt CylPoint;
102     gp_Vec OrigToLine(Origin, LinPoint);
103     if (OrigToLine.Magnitude() <= gp::Resolution())
104     {
105       u = 0.;
106       v = 0.;
107       CylPoint = ElSLib::Value(u, v, S);
108     }
109     else
110     {
111       OrigToLine.Normalize();
112       CylPoint = Origin.Translated(radius * OrigToLine);
113       ElSLib::CylinderParameters(Pos, radius, CylPoint, u, v);
114     }
115     Extrema_POnSurf PonS(u, v, CylPoint);
116     myPoint2->SetValue(1, PonS);
117     myDone = Standard_True;
118     myIsPar = Standard_True;
119   }
120   else {
121     Standard_Integer i, aStartIdx = 0;
122
123     Extrema_POnCurv myPOnC1, myPOnC2;
124     Extrem.Points(1, myPOnC1, myPOnC2);
125     gp_Pnt PonAxis = myPOnC1.Value();
126     gp_Pnt PC = myPOnC2.Value();
127
128     // line intersects the cylinder
129     if (radius - PonAxis.Distance(PC) > Precision::PConfusion())
130     {
131       IntAna_Quadric theQuadric(S);
132       IntAna_IntConicQuad Inters(C, theQuadric);
133       if (Inters.IsDone())
134       {
135         myNbExt = Inters.NbPoints();
136         aStartIdx = myNbExt;
137         if (myNbExt > 0)
138         {
139           // Not more than 2 additional points from perpendiculars.
140           mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
141           myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
142           myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);
143           Standard_Real u, v, w;
144           for (i = 1; i <= myNbExt; i++)
145           {
146             mySqDist->SetValue(i, 0.);
147             gp_Pnt P_int = Inters.Point(i);
148             w = Inters.ParamOnConic(i);
149             Extrema_POnCurv PonC(w, P_int);
150             myPoint1->SetValue(i, PonC);
151             ElSLib::CylinderParameters(Pos, radius, P_int, u, v);
152             Extrema_POnSurf PonS(u, v, P_int);
153             myPoint2->SetValue(i, PonS);
154           }
155         }
156       }
157     }
158
159     // line is tangent or outside of the cylinder
160     Extrema_ExtPElS ExPS(PC, S, Precision::Confusion());
161       if (ExPS.IsDone())
162       {
163         if (aStartIdx == 0)
164         {
165           myNbExt = ExPS.NbExt();
166           mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
167           myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
168           myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
169         }
170         else
171           myNbExt += ExPS.NbExt();
172
173         for (i = aStartIdx + 1; i <= myNbExt; i++) {
174           myPoint1->SetValue(i, myPOnC2);
175           myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
176           mySqDist->SetValue(i,(myPOnC2.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
177         }
178       }
179     myDone = Standard_True;
180   }
181
182 }
183
184
185
186 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
187                                  const gp_Cone& S)
188 {
189   Perform(C, S);
190 }
191
192
193
194 //void Extrema_ExtElCS::Perform(const gp_Lin& C,
195 //                            const gp_Cone& S)
196 void Extrema_ExtElCS::Perform(const gp_Lin& ,
197                               const gp_Cone& )
198 {
199   Standard_NotImplemented::Raise();
200
201 }
202
203
204
205 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
206                                  const gp_Sphere& S)
207 {
208   Perform(C, S);
209 }
210
211
212
213 void Extrema_ExtElCS::Perform(const gp_Lin& C,
214                               const gp_Sphere& S)
215 {
216   // In case of intersection - return four points:
217   // 2 intersection points and 2 perpendicular.
218   // No intersection - only min and max.
219
220   myDone = Standard_False;
221   myNbExt = 0;
222   myIsPar = Standard_False;
223   Standard_Integer aStartIdx = 0;
224
225   gp_Pnt aCenter = S.Location();
226
227   Extrema_ExtPElC Extrem(aCenter, C, Precision::Angular(), RealFirst(), RealLast());
228
229   Standard_Integer i;
230   if (Extrem.IsDone() &&
231       Extrem.NbExt() > 0)
232   {
233     Extrema_POnCurv myPOnC1 =  Extrem.Point(1);
234     if (myPOnC1.Value().Distance(aCenter) <= S.Radius())
235     {
236       IntAna_IntConicQuad aLinSphere(C, S);
237       if (aLinSphere.IsDone())
238       {
239         myNbExt = aLinSphere.NbPoints();
240         aStartIdx = myNbExt;
241         // Not more than 2 additional points from perpendiculars.
242         mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
243         myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
244         myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);
245
246         for (i = 1; i <= myNbExt; i++)
247         {
248           Extrema_POnCurv aCPnt(aLinSphere.ParamOnConic(i), aLinSphere.Point(i));
249
250           Standard_Real u,v;
251           ElSLib::Parameters(S, aLinSphere.Point(i), u, v);
252           Extrema_POnSurf aSPnt(u, v, aLinSphere.Point(i));
253
254           myPoint1->SetValue(i, aCPnt);
255           myPoint2->SetValue(i, aSPnt);
256           mySqDist->SetValue(i,(aCPnt.Value()).SquareDistance(aSPnt.Value()));
257         }
258       }
259     }
260
261     Extrema_ExtPElS ExPS(myPOnC1.Value(), S, Precision::Confusion());
262     if (ExPS.IsDone())
263     {
264       if (aStartIdx == 0)
265       {
266         myNbExt = ExPS.NbExt();
267         mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
268         myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
269         myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
270       }
271       else
272         myNbExt += ExPS.NbExt();
273
274       for (i = aStartIdx + 1; i <= myNbExt; i++)
275       {
276         myPoint1->SetValue(i, myPOnC1);
277         myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
278         mySqDist->SetValue(i,(myPOnC1.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
279       }
280     }
281   }
282   myDone = Standard_True;
283 }
284
285
286 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
287                                  const gp_Torus& S)
288 {
289   Perform(C, S);
290 }
291
292
293
294 //void Extrema_ExtElCS::Perform(const gp_Lin& C,
295 //                            const gp_Torus& S)
296 void Extrema_ExtElCS::Perform(const gp_Lin& ,
297                               const gp_Torus& )
298 {
299   Standard_NotImplemented::Raise();
300
301 }
302
303
304 //        Circle-?
305
306 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
307                                  const gp_Pln& S)
308 {
309   Perform(C, S);
310 }
311
312
313
314 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
315 //                            const gp_Pln& S)
316 void Extrema_ExtElCS::Perform(const gp_Circ& ,
317                               const gp_Pln& )
318 {
319   Standard_NotImplemented::Raise();
320
321 }
322
323
324
325 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
326                                  const gp_Cylinder& S)
327 {
328   Perform(C, S);
329 }
330
331
332
333 //  Modified by skv - Thu Jul  7 14:37:05 2005 OCC9134 Begin
334 // Implementation of the method.
335 void Extrema_ExtElCS::Perform(const gp_Circ& C,
336                               const gp_Cylinder& S)
337 {
338   myDone  = Standard_False;
339   myIsPar = Standard_False;
340   myNbExt = 0;
341
342   // Get an axis line of the cylinder.
343   gp_Lin anAxis(S.Axis());
344
345   // Compute extrema between the circle and the line.
346   Extrema_ExtElC anExtC(anAxis, C, 0.);
347
348   if (anExtC.IsDone()) {
349     if (anExtC.IsParallel()) {
350       myIsPar =     Standard_True;
351       mySqDist = new TColStd_HArray1OfReal(1, 1);
352       Standard_Real aDist = sqrt (anExtC.SquareDistance(1)) - S.Radius();
353       mySqDist->SetValue(1, aDist * aDist);
354     } else {
355       Standard_Integer aNbExt   = anExtC.NbExt();
356       Standard_Integer i;
357       Standard_Integer aCurI    = 1;
358       Standard_Real    aTolConf = Precision::Confusion();
359       Standard_Real    aCylRad  = S.Radius();
360
361       // Check whether two objects have intersection points
362       IntAna_Quadric aCylQuad(S);
363       IntAna_IntConicQuad aCircCylInter(C, aCylQuad);
364       Standard_Integer aNbInter = aCircCylInter.NbPoints();
365       if (!aCircCylInter.IsDone())
366         aNbInter = 0;
367
368       // Compute the extremas.
369       myNbExt  =     2*aNbExt + aNbInter;
370       mySqDist  = new TColStd_HArray1OfReal(1, myNbExt);
371       myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
372       myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
373
374       for (i = 1; i <= aNbExt; i++) {
375         Extrema_POnCurv aPOnAxis;
376         Extrema_POnCurv aPOnCirc;
377         Standard_Real   aSqDist = anExtC.SquareDistance(i);
378         Standard_Real   aDist = sqrt (aSqDist);
379
380         anExtC.Points(i, aPOnAxis, aPOnCirc);
381
382         if (aSqDist <= (aTolConf * aTolConf)) {
383           myNbExt -= 2;
384           continue;
385         }
386
387         gp_Dir aDir(aPOnAxis.Value().XYZ().Subtracted(aPOnCirc.Value().XYZ()));
388         Standard_Real aShift[2] = { aDist + aCylRad, aDist - aCylRad };
389         Standard_Integer j;
390
391         for (j = 0; j < 2; j++) {
392           gp_Vec aVec(aDir);
393           gp_Pnt aPntOnCyl;
394
395           aVec.Multiply(aShift[j]);
396           aPntOnCyl = aPOnCirc.Value().Translated(aVec);
397
398           Standard_Real aU;
399           Standard_Real aV;
400
401           ElSLib::Parameters(S, aPntOnCyl, aU, aV);
402
403           Extrema_POnSurf aPOnSurf(aU, aV, aPntOnCyl);
404
405           myPoint1->SetValue(aCurI, aPOnCirc);
406           myPoint2->SetValue(aCurI, aPOnSurf);
407           mySqDist->SetValue(aCurI++, aShift[j] * aShift[j]);
408         }
409       }
410
411       // Adding intersection points to the list of extremas
412       for (i=1; i<=aNbInter; i++)
413       {
414         Standard_Real aU;
415         Standard_Real aV;
416
417         gp_Pnt aInterPnt = aCircCylInter.Point(i);
418
419         aU = ElCLib::Parameter(C, aInterPnt);
420         Extrema_POnCurv aPOnCirc(aU, aInterPnt);
421
422         ElSLib::Parameters(S, aInterPnt, aU, aV);
423         Extrema_POnSurf aPOnCyl(aU, aV, aInterPnt);
424         myPoint1->SetValue(aCurI, aPOnCirc);
425         myPoint2->SetValue(aCurI, aPOnCyl);
426         mySqDist->SetValue(aCurI++, 0.0);
427       }
428     }
429
430     myDone = Standard_True;
431   }
432 }
433 //  Modified by skv - Thu Jul  7 14:37:05 2005 OCC9134 End
434
435
436
437 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
438                                  const gp_Cone& S)
439 {
440   Perform(C, S);
441 }
442
443
444
445 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
446 //                       const gp_Cone& S)
447 void Extrema_ExtElCS::Perform(const gp_Circ& ,
448                          const gp_Cone& )
449 {
450   Standard_NotImplemented::Raise();
451
452 }
453
454
455
456 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
457                                  const gp_Sphere& S)
458 {
459   Perform(C, S);
460 }
461
462
463
464 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
465 //                            const gp_Sphere& S)
466 void Extrema_ExtElCS::Perform(const gp_Circ& ,
467                               const gp_Sphere& )
468 {
469   Standard_NotImplemented::Raise();
470
471 }
472
473 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
474                                  const gp_Torus& S)
475 {
476   Perform(C, S);
477 }
478
479
480
481 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
482 //                            const gp_Torus& S)
483 void Extrema_ExtElCS::Perform(const gp_Circ& ,
484                               const gp_Torus& )
485 {
486   Standard_NotImplemented::Raise();
487
488 }
489
490 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Hypr& C,
491                                  const gp_Pln& S)
492 {
493   Perform(C, S);
494 }
495
496
497
498 void Extrema_ExtElCS::Perform(const gp_Hypr& C,
499                               const gp_Pln& S)
500 {
501   myDone = Standard_True;
502   myIsPar = Standard_False;
503
504   gp_Ax2 Pos = C.Position();
505   gp_Dir NHypr = Pos.Direction();
506   gp_Dir NPln = S.Axis().Direction();
507
508   if (NHypr.IsParallel(NPln, Precision::Angular())) {
509
510     mySqDist = new TColStd_HArray1OfReal(1, 1);
511     mySqDist->SetValue(1, S.SquareDistance(C.Location()));
512     myIsPar = Standard_True;
513
514   }
515   else {
516
517     gp_Dir XDir = Pos.XDirection();
518     gp_Dir YDir = Pos.YDirection();
519
520     Standard_Real A = C.MinorRadius()*(NPln.Dot(YDir)); 
521     Standard_Real B = C.MajorRadius()*(NPln.Dot(XDir)); 
522
523     if(Abs(B) > Abs(A)) {
524       Standard_Real T = -0.5 * Log((A+B)/(B-A));
525       gp_Pnt Ph = ElCLib::HyperbolaValue(T, Pos, C.MajorRadius(), C.MinorRadius());
526       Extrema_POnCurv PC(T, Ph);
527       myPoint1 = new Extrema_HArray1OfPOnCurv(1,1);
528       myPoint1->SetValue(1, PC);
529
530       mySqDist = new TColStd_HArray1OfReal(1, 1);
531       mySqDist->SetValue(1, S.SquareDistance(Ph));
532
533       Standard_Real U, V;
534       ElSLib::PlaneParameters(S.Position(), Ph, U, V);
535       gp_Pnt Pp = ElSLib::PlaneValue(U, V, S.Position());
536       Extrema_POnSurf PS(U, V, Pp);
537       myPoint2 = new Extrema_HArray1OfPOnSurf(1,1);
538       myPoint2->SetValue(1, PS);
539
540       myNbExt = 1;
541     }
542     else {
543       myNbExt = 0;
544     }
545
546   }
547   
548 }
549
550
551 Standard_Boolean Extrema_ExtElCS::IsDone() const
552 {
553   return myDone;
554 }
555
556
557 Standard_Integer Extrema_ExtElCS::NbExt() const
558 {
559   if (myIsPar) StdFail_InfiniteSolutions::Raise();
560   return myNbExt;
561 }
562
563 Standard_Real Extrema_ExtElCS::SquareDistance(const Standard_Integer N) const
564 {
565   if (myIsPar && N != 1) StdFail_InfiniteSolutions::Raise();
566   return mySqDist->Value(N);
567 }
568
569
570 void Extrema_ExtElCS::Points(const Standard_Integer N,
571                              Extrema_POnCurv& P1,
572                              Extrema_POnSurf& P2) const
573 {
574   if (myIsPar) StdFail_InfiniteSolutions::Raise();
575   P1 = myPoint1->Value(N);
576   P2 = myPoint2->Value(N);
577 }
578
579
580 Standard_Boolean Extrema_ExtElCS::IsParallel() const
581 {
582   return myIsPar;
583 }