cc3211ca8ee436324ba6e74df15aa8fe2b4e4de5
[occt.git] / src / StdPrs / StdPrs_Isolines.cxx
1 // Created on: 2014-10-14
2 // Created by: Anton POLETAEV
3 // Copyright (c) 2013-2014 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16 #include <StdPrs_Isolines.hxx>
17
18 #include <Adaptor3d_IsoCurve.hxx>
19 #include <BRepTools.hxx>
20 #include <BRep_Tool.hxx>
21 #include <ElCLib.hxx>
22 #include <ElSLib.hxx>
23 #include <GCE2d_MakeLine.hxx>
24 #include <gce_MakeLin2d.hxx>
25 #include <gce_MakePln.hxx>
26 #include <gce_MakeLin.hxx>
27 #include <GCPnts_AbscissaPoint.hxx>
28 #include <GCPnts_QuasiUniformDeflection.hxx>
29 #include <Geom_BezierSurface.hxx>
30 #include <Geom_BSplineSurface.hxx>
31 #include <Geom_Plane.hxx>
32 #include <Geom_Line.hxx>
33 #include <Geom2d_Line.hxx>
34 #include <Geom2dAdaptor_Curve.hxx>
35 #include <Geom2dInt_GInter.hxx>
36 #include <GeomAdaptor_Curve.hxx>
37 #include <GeomAPI_IntCS.hxx>
38 #include <GeomLib.hxx>
39 #include <GeomLib_Tool.hxx>
40 #include <gp_Lin2d.hxx>
41 #include <Graphic3d_ArrayOfSegments.hxx>
42 #include <Graphic3d_ArrayOfPolylines.hxx>
43 #include <Hatch_Hatcher.hxx>
44 #include <IntRes2d_IntersectionPoint.hxx>
45 #include <NCollection_List.hxx>
46 #include <NCollection_QuickSort.hxx>
47 #include <ProjLib.hxx>
48 #include <Prs3d_IsoAspect.hxx>
49 #include <Prs3d_NListOfSequenceOfPnt.hxx>
50 #include <Prs3d_NListIteratorOfListOfSequenceOfPnt.hxx>
51 #include <Poly_Array1OfTriangle.hxx>
52 #include <Poly_Triangulation.hxx>
53 #include <StdPrs_DeflectionCurve.hxx>
54 #include <StdPrs_ToolRFace.hxx>
55 #include <TColgp_Array1OfPnt.hxx>
56 #include <TColgp_SequenceOfPnt2d.hxx>
57 #include <Standard_ErrorHandler.hxx>
58 #include <Geom_Surface.hxx>
59
60 typedef NCollection_Sequence<Handle(TColgp_HSequenceOfPnt)> Prs3d_NSequenceOfSequenceOfPnt;
61
62 namespace
63 {
64   const gp_Lin2d isoU (const Standard_Real theU) { return gp_Lin2d (gp_Pnt2d (theU, 0.0), gp::DY2d()); }
65   const gp_Lin2d isoV (const Standard_Real theV) { return gp_Lin2d (gp_Pnt2d (0.0, theV), gp::DX2d()); }
66
67   //! Assembles array of primitives for sequence of polyine points.
68   //! @param thePoints [in] the polyline points.
69   //! @return array of primitives.
70   template <typename T>
71   inline Handle(T) primitivesForPolyline (const Prs3d_NSequenceOfSequenceOfPnt& thePoints)
72   {
73     if (thePoints.IsEmpty())
74     {
75       return Handle(T)();
76     }
77
78     Standard_Integer aNbBounds   = thePoints.Size();
79     Standard_Integer aNbVertices = 0;
80     for (Prs3d_NSequenceOfSequenceOfPnt::Iterator anIt (thePoints); anIt.More(); anIt.Next())
81     {
82       aNbVertices += anIt.Value()->Length();
83     }
84     Handle(T) aPrimitives = new T (aNbVertices, aNbBounds);
85     for (NCollection_Sequence<Handle(TColgp_HSequenceOfPnt)>::Iterator anIt (thePoints); anIt.More(); anIt.Next())
86     {
87       const Handle(TColgp_HSequenceOfPnt)& aPoints = anIt.Value();
88       aPrimitives->AddBound (aPoints->Length());
89       for (Standard_Integer anI = 1; anI <= aPoints->Length(); ++anI)
90       {
91         aPrimitives->AddVertex (aPoints->Value (anI));
92       }
93     }
94
95     return aPrimitives;
96   }
97
98   //! Reoder and adjust to the limit a curve's parameter values.
99   //! @param theCurve [in] the curve.
100   //! @param theLimit [in] the parameter limit value.
101   //! @param theFirst [in/out] the first parameter value.
102   //! @param theLast  [in/out] the last parameter value.
103   static void findLimits (const Adaptor3d_Curve& theCurve,
104                           const Standard_Real    theLimit,
105                           Standard_Real&         theFirst,
106                           Standard_Real&         theLast)
107   {
108     theFirst = Max (theCurve.FirstParameter(), theFirst);
109     theLast  = Min (theCurve.LastParameter(), theLast);
110
111     Standard_Boolean isFirstInf = Precision::IsNegativeInfinite (theFirst);
112     Standard_Boolean isLastInf  = Precision::IsPositiveInfinite (theLast);
113
114     if (!isFirstInf && !isLastInf)
115     {
116       return;
117     }
118
119     gp_Pnt aP1, aP2;
120     Standard_Real aDelta = 1.0;
121     if (isFirstInf && isLastInf)
122     {
123       do
124       {
125         aDelta *= 2.0;
126         theFirst = -aDelta;
127         theLast  =  aDelta;
128         theCurve.D0 (theFirst, aP1);
129         theCurve.D0 (theLast,  aP2);
130       }
131       while (aP1.Distance (aP2) < theLimit);
132     }
133     else if (isFirstInf)
134     {
135       theCurve.D0 (theLast, aP2);
136       do
137       {
138         aDelta *= 2.0;
139         theFirst = theLast - aDelta;
140         theCurve.D0 (theFirst, aP1);
141       }
142       while (aP1.Distance (aP2) < theLimit);
143     }
144     else if (isLastInf)
145     {
146       theCurve.D0 (theFirst, aP1);
147       do
148       {
149         aDelta *= 2.0;
150         theLast = theFirst + aDelta;
151         theCurve.D0 (theLast, aP2);
152       }
153       while (aP1.Distance (aP2) < theLimit);
154     }
155   }
156
157 }
158
159 //==================================================================
160 // function : AddOnTriangulation
161 // purpose  :
162 //==================================================================
163 void StdPrs_Isolines::AddOnTriangulation (const Handle(Prs3d_Presentation)& thePresentation,
164                                           const TopoDS_Face&                theFace,
165                                           const Handle(Prs3d_Drawer)&       theDrawer)
166 {
167   const Standard_Integer aNbIsoU = theDrawer->UIsoAspect()->Number();
168   const Standard_Integer aNbIsoV = theDrawer->VIsoAspect()->Number();
169   if (aNbIsoU < 1 && aNbIsoV < 1)
170   {
171     return;
172   }
173
174   // Evalute parameters for uv isolines.
175   TColStd_SequenceOfReal aUIsoParams;
176   TColStd_SequenceOfReal aVIsoParams;
177   UVIsoParameters (theFace, aNbIsoU, aNbIsoV, theDrawer->MaximalParameterValue(), aUIsoParams, aVIsoParams);
178
179   // Access surface definition.
180   TopLoc_Location aLocSurface;
181   Handle(Geom_Surface) aSurface = BRep_Tool::Surface (theFace, aLocSurface);
182
183   // Access triangulation.
184   TopLoc_Location aLocTriangulation;
185   const Handle(Poly_Triangulation)& aTriangulation = BRep_Tool::Triangulation (theFace, aLocTriangulation);
186   if (aTriangulation.IsNull())
187   {
188     return;
189   }
190
191   // Setup equal location for surface and triangulation.
192   if (!aLocTriangulation.IsEqual (aLocSurface))
193   {
194     aSurface = Handle (Geom_Surface)::DownCast (
195       aSurface->Transformed ((aLocSurface / aLocTriangulation).Transformation()));
196   }
197
198   AddOnTriangulation (thePresentation,
199                       aTriangulation,
200                       aSurface,
201                       aLocTriangulation,
202                       theDrawer,
203                       aUIsoParams,
204                       aVIsoParams);
205 }
206
207 //==================================================================
208 // function : AddOnTriangulation
209 // purpose  :
210 //==================================================================
211 void StdPrs_Isolines::AddOnTriangulation (const Handle(Prs3d_Presentation)& thePresentation,
212                                           const Handle(Poly_Triangulation)& theTriangulation,
213                                           const Handle(Geom_Surface)&       theSurface,
214                                           const TopLoc_Location&            theLocation,
215                                           const Handle(Prs3d_Drawer)&       theDrawer,
216                                           const TColStd_SequenceOfReal&     theUIsoParams,
217                                           const TColStd_SequenceOfReal&     theVIsoParams)
218 {
219   const Standard_Integer aNbIsoU = theUIsoParams.Length();
220   const Standard_Integer aNbIsoV = theVIsoParams.Length();
221
222   Prs3d_NSequenceOfSequenceOfPnt aUPolylines;
223   Prs3d_NSequenceOfSequenceOfPnt aVPolylines;
224
225   const Poly_Array1OfTriangle& aTriangles = theTriangulation->Triangles();
226   const TColgp_Array1OfPnt&    aNodes     = theTriangulation->Nodes();
227   const TColgp_Array1OfPnt2d&  aUVNodes   = theTriangulation->UVNodes();
228
229   TColStd_Array1OfInteger aUIsoIndexes (1, aNbIsoU);
230   TColStd_Array1OfInteger aVIsoIndexes (1, aNbIsoV);
231   aUIsoIndexes.Init (-1);
232   aVIsoIndexes.Init (-1);
233
234   for (Standard_Integer anI = aTriangles.Lower(); anI <= aTriangles.Upper(); ++anI)
235   {
236     Standard_Integer aNodeIdxs[3];
237     aTriangles.Value (anI).Get (aNodeIdxs[0], aNodeIdxs[1],aNodeIdxs[2]);
238     const gp_Pnt aNodesXYZ[3] = { aNodes.Value (aNodeIdxs[0]),
239                                   aNodes.Value (aNodeIdxs[1]),
240                                   aNodes.Value (aNodeIdxs[2]) };
241     const gp_Pnt2d aNodesUV[3] = { aUVNodes.Value (aNodeIdxs[0]),
242                                    aUVNodes.Value (aNodeIdxs[1]),
243                                    aUVNodes.Value (aNodeIdxs[2]) };
244
245     // Evaluate polyline points for u isolines.
246     for (Standard_Integer anIsoIdx = 1; anIsoIdx <= aNbIsoU; ++anIsoIdx)
247     {
248       gp_Pnt aSegment[2];
249       const gp_Lin2d anIsolineUV = isoU (theUIsoParams.Value (anIsoIdx));
250
251       // Find intersections with triangle in uv space and its projection on triangulation.
252       if (!findSegmentOnTriangulation (theSurface, anIsolineUV, aNodesXYZ, aNodesUV, aSegment))
253       {
254         continue;
255       }
256
257       if (aUIsoIndexes.Value (anIsoIdx) == -1)
258       {
259         aUPolylines.Append (new TColgp_HSequenceOfPnt());
260         aUIsoIndexes.SetValue (anIsoIdx, aUPolylines.Size());
261       }
262
263       Handle(TColgp_HSequenceOfPnt) anIsoPnts = aUPolylines.ChangeValue (aUIsoIndexes.Value (anIsoIdx));
264       anIsoPnts->Append (theLocation.IsIdentity() ? aSegment[0] : aSegment[0].Transformed (theLocation));
265       anIsoPnts->Append (theLocation.IsIdentity() ? aSegment[1] : aSegment[1].Transformed (theLocation));
266     }
267
268     // Evaluate polyline points for v isolines.
269     for (Standard_Integer anIsoIdx = 1; anIsoIdx <= aNbIsoV; ++anIsoIdx)
270     {
271       gp_Pnt aSegment[2];
272       const gp_Lin2d anIsolineUV = isoV (theVIsoParams.Value (anIsoIdx));
273
274       if (!findSegmentOnTriangulation (theSurface, anIsolineUV, aNodesXYZ, aNodesUV, aSegment))
275       {
276         continue;
277       }
278
279       if (aVIsoIndexes.Value (anIsoIdx) == -1)
280       {
281         aVPolylines.Append (new TColgp_HSequenceOfPnt());
282         aVIsoIndexes.SetValue (anIsoIdx, aVPolylines.Size());
283       }
284
285       Handle(TColgp_HSequenceOfPnt) anIsoPnts = aVPolylines.ChangeValue (aVIsoIndexes.Value (anIsoIdx));
286       anIsoPnts->Append (theLocation.IsIdentity() ? aSegment[0] : aSegment[0].Transformed (theLocation));
287       anIsoPnts->Append (theLocation.IsIdentity() ? aSegment[1] : aSegment[1].Transformed (theLocation));
288     }
289   }
290
291   // Add primitive arrays for isoline segments.
292   Handle(Graphic3d_ArrayOfSegments) aUPrimitives = primitivesForPolyline<Graphic3d_ArrayOfSegments> (aUPolylines);
293   Handle(Graphic3d_ArrayOfSegments) aVPrimitives = primitivesForPolyline<Graphic3d_ArrayOfSegments> (aVPolylines);
294
295   if (!aUPrimitives.IsNull())
296   {
297     Handle(Graphic3d_Group) aGroup = Prs3d_Root::NewGroup (thePresentation);
298     aGroup->SetPrimitivesAspect (theDrawer->UIsoAspect()->Aspect());
299     aGroup->AddPrimitiveArray (aUPrimitives);
300   }
301
302   if (!aVPrimitives.IsNull())
303   {
304     Handle(Graphic3d_Group) aGroup = Prs3d_Root::NewGroup (thePresentation);
305     aGroup->SetPrimitivesAspect (theDrawer->VIsoAspect()->Aspect());
306     aGroup->AddPrimitiveArray (aVPrimitives);
307   }
308 }
309
310 //==================================================================
311 // function : AddOnSurface
312 // purpose  :
313 //==================================================================
314 void StdPrs_Isolines::AddOnSurface (const Handle(Prs3d_Presentation)& thePresentation,
315                                     const TopoDS_Face&                theFace,
316                                     const Handle(Prs3d_Drawer)&       theDrawer,
317                                     const Standard_Real               theDeflection)
318 {
319   const Standard_Integer aNbIsoU = theDrawer->UIsoAspect()->Number();
320   const Standard_Integer aNbIsoV = theDrawer->VIsoAspect()->Number();
321   if (aNbIsoU < 1 && aNbIsoV < 1)
322   {
323     return;
324   }
325
326   // Evalute parameters for uv isolines.
327   TColStd_SequenceOfReal aUIsoParams;
328   TColStd_SequenceOfReal aVIsoParams;
329   UVIsoParameters (theFace, aNbIsoU, aNbIsoV, theDrawer->MaximalParameterValue(), aUIsoParams, aVIsoParams);
330
331   BRepAdaptor_Surface aSurface (theFace);
332   AddOnSurface (thePresentation,
333                 new BRepAdaptor_HSurface (aSurface),
334                 theDrawer,
335                 theDeflection,
336                 aUIsoParams,
337                 aVIsoParams);
338 }
339
340 //==================================================================
341 // function : AddOnSurface
342 // purpose  :
343 //==================================================================
344 void StdPrs_Isolines::AddOnSurface (const Handle(Prs3d_Presentation)&   thePresentation,
345                                     const Handle(BRepAdaptor_HSurface)& theSurface,
346                                     const Handle(Prs3d_Drawer)&         theDrawer,
347                                     const Standard_Real                 theDeflection,
348                                     const TColStd_SequenceOfReal&       theUIsoParams,
349                                     const TColStd_SequenceOfReal&       theVIsoParams)
350 {
351   // Choose a deflection for sampling edge uv curves.
352   Standard_Real aUVLimit = theDrawer->MaximalParameterValue();
353   Standard_Real aUmin  = Max (theSurface->FirstUParameter(), -aUVLimit);
354   Standard_Real aUmax  = Min (theSurface->LastUParameter(),   aUVLimit);
355   Standard_Real aVmin  = Max (theSurface->FirstVParameter(), -aUVLimit);
356   Standard_Real aVmax  = Min (theSurface->LastVParameter(),   aUVLimit);
357   Standard_Real aSamplerDeflection = Max (aUmax - aUmin, aVmax - aVmin) * theDrawer->DeviationCoefficient();
358   Standard_Real aHatchingTolerance = RealLast();
359   Prs3d_NSequenceOfSequenceOfPnt aUPolylines;
360   Prs3d_NSequenceOfSequenceOfPnt aVPolylines;
361
362   try
363   {
364     OCC_CATCH_SIGNALS
365     // Determine edge points for trimming uv hatch region.
366     TColgp_SequenceOfPnt2d aTrimPoints;
367     StdPrs_ToolRFace anEdgeTool (theSurface);
368     for (anEdgeTool.Init(); anEdgeTool.More(); anEdgeTool.Next())
369     {
370       TopAbs_Orientation anOrientation = anEdgeTool.Orientation();
371       if (anOrientation != TopAbs_FORWARD && anOrientation != TopAbs_REVERSED)
372       {
373         continue;
374       }
375
376       Adaptor2d_Curve2dPtr anEdgeCurve = anEdgeTool.Value();
377       if (anEdgeCurve->GetType() != GeomAbs_Line)
378       {
379         GCPnts_QuasiUniformDeflection aSampler (*anEdgeCurve, aSamplerDeflection);
380         if (!aSampler.IsDone())
381         {
382 #ifdef OCCT_DEBUG
383           std::cout << "Cannot evaluate curve on surface" << std::endl;
384 #endif
385           continue;
386         }
387
388         Standard_Integer aNumberOfPoints = aSampler.NbPoints();
389         if (aNumberOfPoints < 2)
390         {
391           continue;
392         }
393
394         for (Standard_Integer anI = 1; anI < aNumberOfPoints; ++anI)
395         {
396           gp_Pnt2d aP1 (aSampler.Value (anI    ).X(), aSampler.Value (anI    ).Y());
397           gp_Pnt2d aP2 (aSampler.Value (anI + 1).X(), aSampler.Value (anI + 1).Y());
398
399           aHatchingTolerance = Min (aP1.SquareDistance (aP2), aHatchingTolerance);
400
401           aTrimPoints.Append (anOrientation == TopAbs_FORWARD ? aP1 : aP2);
402           aTrimPoints.Append (anOrientation == TopAbs_FORWARD ? aP2 : aP1);
403         }
404       }
405       else
406       {
407         Standard_Real aU1 = anEdgeCurve->FirstParameter();
408         Standard_Real aU2 = anEdgeCurve->LastParameter();
409
410         // MSV 17.08.06 OCC13144: U2 occured less than U1, to overcome it
411         // ensure that distance U2-U1 is not greater than aLimit*2,
412         // if greater then choose an origin and use aLimit to define
413         // U1 and U2 anew.
414         Standard_Real anOrigin = 0.0;
415
416         if (!Precision::IsNegativeInfinite (aU1) || !Precision::IsPositiveInfinite(aU2))
417         {
418           if (Precision::IsNegativeInfinite (aU1))
419           {
420             anOrigin = aU2 - aUVLimit;
421           }
422           else if (Precision::IsPositiveInfinite (aU2))
423           {
424             anOrigin = aU1 + aUVLimit;
425           }
426           else
427           {
428             anOrigin = (aU1 + aU2) * 0.5;
429           }
430         }
431
432         aU1 = Max (anOrigin - aUVLimit, aU1);
433         aU2 = Min (anOrigin + aUVLimit, aU2);
434
435         gp_Pnt2d aP1 = anEdgeCurve->Value (aU1);
436         gp_Pnt2d aP2 = anEdgeCurve->Value (aU2);
437
438         aHatchingTolerance = Min (aP1.SquareDistance(aP2), aHatchingTolerance);
439
440         aTrimPoints.Append (anOrientation == TopAbs_FORWARD ? aP1 : aP2);
441         aTrimPoints.Append (anOrientation == TopAbs_FORWARD ? aP2 : aP1);
442       }
443     }
444
445     // Compute a hatching tolerance.
446     aHatchingTolerance *= 0.1;
447     aHatchingTolerance = Max (Precision::Confusion(), aHatchingTolerance);
448     aHatchingTolerance = Min (1.0E-5, aHatchingTolerance);
449
450     // Load isolines into hatcher.
451     Hatch_Hatcher aHatcher (aHatchingTolerance, anEdgeTool.IsOriented());
452
453     for (Standard_Integer anIso = 1; anIso <= theUIsoParams.Length(); ++anIso)
454     {
455       aHatcher.AddXLine (theUIsoParams.Value (anIso));
456     }
457     for (Standard_Integer anIso = 1; anIso <= theVIsoParams.Length(); ++anIso)
458     {
459       aHatcher.AddYLine (theVIsoParams.Value (anIso));
460     }
461
462     // Trim hatching region.
463     for (Standard_Integer anI = 1; anI <= aTrimPoints.Length(); anI += 2)
464     {
465       aHatcher.Trim (aTrimPoints (anI), aTrimPoints (anI + 1));
466     }
467
468     // Use surface definition for evaluation of Bezier, B-spline surface.
469     // Use isoline adapter for other types of surfaces.
470     GeomAbs_SurfaceType  aSurfType = theSurface->GetType();
471     Handle(Geom_Surface) aBSurface;
472     GeomAdaptor_Curve    aBSurfaceCurve;
473     Adaptor3d_IsoCurve   aCanonicalCurve;
474     if (aSurfType == GeomAbs_BezierSurface)
475     {
476       aBSurface = theSurface->Bezier();
477     }
478     else if (aSurfType == GeomAbs_BSplineSurface)
479     {
480       aBSurface = theSurface->BSpline();
481     }
482     else
483     {
484       aCanonicalCurve.Load (theSurface);
485     }
486
487     // For each isoline: compute its segments.
488     for (Standard_Integer anI = 1; anI <= aHatcher.NbLines(); anI++)
489     {
490       Standard_Real anIsoParam = aHatcher.Coordinate (anI);
491       Standard_Boolean isIsoU  = aHatcher.IsXLine (anI);
492
493       // For each isoline's segment: evaluate its points.
494       for (Standard_Integer aJ = 1; aJ <= aHatcher.NbIntervals (anI); aJ++)
495       {
496         Standard_Real aSegmentP1 = aHatcher.Start (anI, aJ);
497         Standard_Real aSegmentP2 = aHatcher.End (anI, aJ);
498
499         if (!aBSurface.IsNull())
500         {
501           aBSurfaceCurve.Load (isIsoU ? aBSurface->UIso (anIsoParam) : aBSurface->VIso (anIsoParam));
502
503           findLimits (aBSurfaceCurve, aUVLimit, aSegmentP1, aSegmentP2);
504
505           if (aSegmentP2 - aSegmentP1 <= Precision::Confusion())
506           {
507             continue;
508           }
509         }
510         else
511         {
512           aCanonicalCurve.Load (isIsoU ? GeomAbs_IsoU : GeomAbs_IsoV, anIsoParam, aSegmentP1, aSegmentP2);
513
514           findLimits (aCanonicalCurve, aUVLimit, aSegmentP1, aSegmentP2);
515
516           if (aSegmentP2 - aSegmentP1 <= Precision::Confusion())
517           {
518             continue;
519           }
520         }
521         Adaptor3d_Curve* aCurve = aBSurface.IsNull() ? (Adaptor3d_Curve*) &aCanonicalCurve
522           : (Adaptor3d_Curve*) &aBSurfaceCurve;
523
524         Handle(TColgp_HSequenceOfPnt) aPoints = new TColgp_HSequenceOfPnt();
525         StdPrs_DeflectionCurve::Add (thePresentation,
526                                      *aCurve,
527                                      aSegmentP1,
528                                      aSegmentP2,
529                                      theDeflection,
530                                      aPoints->ChangeSequence(),
531                                      theDrawer->DeviationAngle(),
532                                      Standard_False);
533         if (aPoints->IsEmpty())
534         {
535           continue;
536         }
537
538         if (isIsoU)
539         {
540           aUPolylines.Append (aPoints);
541         }
542         else
543         {
544           aVPolylines.Append (aPoints);
545         }
546       }
547     }
548   }
549   catch (Standard_Failure)
550   {
551     // ...
552   }
553
554   // Add primitive arrays for isoline segments.
555   Handle(Graphic3d_ArrayOfPolylines) aUPrimitives = primitivesForPolyline<Graphic3d_ArrayOfPolylines> (aUPolylines);
556   Handle(Graphic3d_ArrayOfPolylines) aVPrimitives = primitivesForPolyline<Graphic3d_ArrayOfPolylines> (aVPolylines);
557
558   if (!aUPrimitives.IsNull())
559   {
560     Handle(Graphic3d_Group) aGroup = Prs3d_Root::NewGroup (thePresentation);
561     aGroup->SetPrimitivesAspect (theDrawer->UIsoAspect()->Aspect());
562     aGroup->AddPrimitiveArray (aUPrimitives);
563   }
564
565   if (!aVPrimitives.IsNull())
566   {
567     Handle(Graphic3d_Group) aGroup = Prs3d_Root::NewGroup (thePresentation);
568     aGroup->SetPrimitivesAspect (theDrawer->VIsoAspect()->Aspect());
569     aGroup->AddPrimitiveArray (aVPrimitives);
570   }
571 }
572
573 //==================================================================
574 // function : UVIsoParameters
575 // purpose  :
576 //==================================================================
577 void StdPrs_Isolines::UVIsoParameters (const TopoDS_Face&      theFace,
578                                        const Standard_Integer  theNbIsoU,
579                                        const Standard_Integer  theNbIsoV,
580                                        const Standard_Real     theUVLimit,
581                                        TColStd_SequenceOfReal& theUIsoParams,
582                                        TColStd_SequenceOfReal& theVIsoParams)
583 {
584   Standard_Real aUmin = 0.0;
585   Standard_Real aUmax = 0.0;
586   Standard_Real aVmin = 0.0;
587   Standard_Real aVmax = 0.0;
588
589   BRepTools::UVBounds (theFace, aUmin, aUmax, aVmin, aVmax);
590
591   aUmin = Max (aUmin, -theUVLimit);
592   aUmax = Min (aUmax,  theUVLimit);
593   aVmin = Max (aVmin, -theUVLimit);
594   aVmax = Min (aVmax,  theUVLimit);
595
596   TopLoc_Location aLocation;
597   Handle(Geom_Surface) aSurface = BRep_Tool::Surface (theFace, aLocation);
598
599   const Standard_Boolean isUClosed = aSurface->IsUClosed();
600   const Standard_Boolean isVClosed = aSurface->IsVClosed();
601
602   if (!isUClosed)
603   {
604     aUmin = aUmin + (aUmax - aUmin) / 1000.0;
605     aUmax = aUmax - (aUmax - aUmin) / 1000.0;
606   }
607
608   if (!isVClosed)
609   {
610     aVmin = aVmin + (aVmax - aVmin) / 1000.0;
611     aVmax = aVmax - (aVmax - aVmin) / 1000.0;
612   }
613
614   Standard_Real aUstep = (aUmax - aUmin) / (1 + theNbIsoU);
615   Standard_Real aVstep = (aVmax - aVmin) / (1 + theNbIsoV);
616
617   for (Standard_Integer anIso = 1; anIso <= theNbIsoU; ++anIso)
618   {
619     theUIsoParams.Append (aUmin + aUstep * anIso);
620   }
621
622   for (Standard_Integer anIso = 1; anIso <= theNbIsoV; ++anIso)
623   {
624     theVIsoParams.Append (aVmin + aVstep * anIso);
625   }
626 }
627
628 //==================================================================
629 // function : FindSegmentOnTriangulation
630 // purpose  :
631 //==================================================================
632 Standard_Boolean StdPrs_Isolines::findSegmentOnTriangulation (const Handle(Geom_Surface)& theSurface,
633                                                               const gp_Lin2d&             theIsoline,
634                                                               const gp_Pnt*               theNodesXYZ,
635                                                               const gp_Pnt2d*             theNodesUV,
636                                                               gp_Pnt*                     theSegment)
637 {
638   Standard_Integer aNPoints = 0;
639
640   for (Standard_Integer aLinkIter = 0; aLinkIter < 3 && aNPoints < 2; ++aLinkIter)
641   {
642     // ...
643     // Check that uv isoline crosses the triangulation link in parametric space
644     // ...
645
646     const gp_Pnt2d& aNodeUV1 = theNodesUV[aLinkIter];
647     const gp_Pnt2d& aNodeUV2 = theNodesUV[(aLinkIter + 1) % 3];
648     const gp_Pnt& aNode1 = theNodesXYZ[aLinkIter];
649     const gp_Pnt& aNode2 = theNodesXYZ[(aLinkIter + 1) % 3];
650
651     // Compute distance of uv points to isoline taking into consideration their relative
652     // location against the isoline (left or right). Null value for a node means that the
653     // isoline crosses the node. Both positive or negative means that the isoline does not
654     // cross the segment.
655     Standard_Boolean isLeftUV1 = (theIsoline.Direction().XY() ^ gp_Vec2d (theIsoline.Location(), aNodeUV1).XY()) > 0.0;
656     Standard_Boolean isLeftUV2 = (theIsoline.Direction().XY() ^ gp_Vec2d (theIsoline.Location(), aNodeUV2).XY()) > 0.0;
657     Standard_Real aDistanceUV1 = isLeftUV1 ? theIsoline.Distance (aNodeUV1) : -theIsoline.Distance (aNodeUV1);
658     Standard_Real aDistanceUV2 = isLeftUV2 ? theIsoline.Distance (aNodeUV2) : -theIsoline.Distance (aNodeUV2);
659
660     // Isoline crosses first point of an edge.
661     if (Abs (aDistanceUV1) < Precision::PConfusion())
662     {
663       theSegment[aNPoints++] = aNode1;
664       continue;
665     }
666
667     // Isoline crosses second point of an edge.
668     if (Abs (aDistanceUV2) < Precision::PConfusion())
669     {
670       theSegment[aNPoints++] = aNode2;
671       aLinkIter++;
672       continue;
673     }
674
675     // Isoline does not cross the triangle link.
676     if (aDistanceUV1 * aDistanceUV2 > 0.0)
677     {
678       continue;
679     }
680
681     // Isoline crosses degenerated link.
682     if (aNode1.SquareDistance (aNode2) < Precision::PConfusion())
683     {
684       theSegment[aNPoints++] = aNode1;
685       continue;
686     }
687
688     // ...
689     // Derive cross-point from parametric coordinates
690     // ...
691
692     Standard_Real anAlpha = Abs (aDistanceUV1) / (Abs (aDistanceUV1) + Abs (aDistanceUV2));
693
694     gp_Pnt aCross (0.0, 0.0, 0.0);
695
696     // Is surface definition available?
697     if (theSurface.IsNull())
698     {
699       // Do linear interpolation of point coordinates using
700       // triangulation nodes.
701       aCross.SetX (aNode1.X() + anAlpha * (aNode2.X() - aNode1.X()));
702       aCross.SetY (aNode1.Y() + anAlpha * (aNode2.Y() - aNode1.Y()));
703       aCross.SetZ (aNode1.Z() + anAlpha * (aNode2.Z() - aNode1.Z()));
704     }
705     else
706     {
707       // Do linear interpolation of point coordinates by triangulation nodes.
708       Standard_Real aCrossU = aNodeUV1.X() + anAlpha * (aNodeUV2.X() - aNodeUV1.X());
709       Standard_Real aCrossV = aNodeUV1.Y() + anAlpha * (aNodeUV2.Y() - aNodeUV1.Y());
710
711       // Get 3d point on surface.
712       Handle(Geom_Curve) anIso1, anIso2, anIso3;
713       Standard_Real aPntOnNode1Iso = 0.0;
714       Standard_Real aPntOnNode2Iso = 0.0;
715       Standard_Real aPntOnNode3Iso = 0.0;
716
717       if (theIsoline.Direction().X() == 0.0)
718       {
719         aPntOnNode1Iso = aNodeUV1.X();
720         aPntOnNode2Iso = aNodeUV2.X();
721         aPntOnNode3Iso = aCrossU;
722         anIso1 = theSurface->VIso (aNodeUV1.Y());
723         anIso2 = theSurface->VIso (aNodeUV2.Y());
724         anIso3 = theSurface->VIso (aCrossV);
725       }
726       else if (theIsoline.Direction().Y() == 0.0)
727       {
728         aPntOnNode1Iso = aNodeUV1.Y();
729         aPntOnNode2Iso = aNodeUV2.Y();
730         aPntOnNode3Iso = aCrossV;
731         anIso1 = theSurface->UIso (aNodeUV1.X());
732         anIso2 = theSurface->UIso (aNodeUV2.X());
733         anIso3 = theSurface->UIso (aCrossU);
734       }
735
736       GeomAdaptor_Curve aCurveAdaptor1 (anIso1);
737       GeomAdaptor_Curve aCurveAdaptor2 (anIso2);
738       GeomAdaptor_Curve aCurveAdaptor3 (anIso3);
739       Standard_Real aLength1 = GCPnts_AbscissaPoint::Length (aCurveAdaptor1, aPntOnNode1Iso, aPntOnNode3Iso, 1e-2);
740       Standard_Real aLength2 = GCPnts_AbscissaPoint::Length (aCurveAdaptor2, aPntOnNode2Iso, aPntOnNode3Iso, 1e-2);
741       if (Abs (aLength1) < Precision::Confusion() || Abs (aLength2) < Precision::Confusion())
742       {
743         theSegment[aNPoints++] = (aNode2.XYZ() - aNode1.XYZ()) * anAlpha + aNode1.XYZ();
744         continue;
745       }
746
747       aCross = (aNode2.XYZ() - aNode1.XYZ()) * (aLength1 / (aLength1 + aLength2)) + aNode1.XYZ();
748     }
749
750     theSegment[aNPoints++] = aCross;
751   }
752
753   return aNPoints == 2;
754 }