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