0031258: Mesh - OCCT 7.4.0 VIS get wrong render data
[occt.git] / src / BRepMesh / BRepMesh_NURBSRangeSplitter.cxx
1 // Created on: 2016-04-19
2 // Copyright (c) 2016 OPEN CASCADE SAS
3 // Created by: Oleg AGASHIN
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 <BRepMesh_NURBSRangeSplitter.hxx>
17
18 #include <algorithm>
19 #include <BRepMesh_GeomTool.hxx>
20 #include <GeomAdaptor_Curve.hxx>
21 #include <GeomLib.hxx>
22 #include <IMeshData_Edge.hxx>
23 #include <IMeshData_Wire.hxx>
24 #include <NCollection_Handle.hxx>
25
26 namespace
27 {
28   class AnalyticalFilter
29   {
30   public:
31     //! Constructor.
32     AnalyticalFilter(
33       const IMeshData::IFaceHandle&             theDFace,
34       const GeomAbs_IsoType                     theIsoType,
35       const Handle(IMeshData::SequenceOfReal)&  theParams,
36       const Handle(IMeshData::SequenceOfReal)&  theControlParams,
37       const Handle(IMeshData::MapOfReal)&       theParamsForbiddenToRemove,
38       const Handle(IMeshData::MapOfReal)&       theControlParamsForbiddenToRemove)
39       : myDFace(theDFace),
40         mySurface(myDFace->GetSurface()->ChangeSurface().Surface().Surface()),
41         myIsoU(theIsoType == GeomAbs_IsoU),
42         myParams(theParams),
43         myControlParams(theControlParams),
44         myParamsForbiddenToRemove(theParamsForbiddenToRemove),
45         myControlParamsForbiddenToRemove(theControlParamsForbiddenToRemove),
46         myAllocator(new NCollection_IncAllocator(IMeshData::MEMORY_BLOCK_SIZE_HUGE)),
47         myControlParamsToRemove(new IMeshData::MapOfReal(1, myAllocator))
48     {
49     }
50
51     //! Returns map of parameters supposed to be removed.
52     const Handle(IMeshData::MapOfReal)& GetControlParametersToRemove(
53       const IMeshTools_Parameters& theParameters)
54     {
55       myParameters = theParameters;
56
57       Standard_Integer aStartIndex, aEndIndex;
58       if (myIsoU)
59       {
60         aStartIndex = 1;
61         aEndIndex = myParams->Length();
62       }
63       else
64       {
65         aStartIndex = 2;
66         aEndIndex = myParams->Length() - 1;
67       }
68
69       for (Standard_Integer i = aStartIndex; i <= aEndIndex; ++i)
70       {
71         myCurrParam = myParams->Value(i);
72         myIso = new GeomAdaptor_Curve(myIsoU ? mySurface->UIso(myCurrParam) : mySurface->VIso(myCurrParam));
73
74         myPrevControlParam = myControlParams->Value(1);
75         myIso->D1(myPrevControlParam, myPrevControlPnt, myPrevControlVec);
76         for (Standard_Integer j = 2; j <= myControlParams->Length();)
77         {
78           j += checkControlPointAndMoveOn(j);
79         }
80       }
81
82       return myControlParamsToRemove;
83     }
84
85   private:
86
87     //! Checks the given control point for deviation.
88     //! Returns number of steps to be used to move point iterator.
89     Standard_Integer checkControlPointAndMoveOn(const Standard_Integer theIndex)
90     {
91       Standard_Integer aMoveSteps = 0;
92       myCurrControlParam = myControlParams->Value(theIndex);
93       myIso->D1(myCurrControlParam, myCurrControlPnt, myCurrControlVec);
94
95       const Standard_Real aMidParam = 0.5 * (myPrevControlParam + myCurrControlParam);
96       const gp_Pnt aMidPnt = myIso->Value(aMidParam);
97
98       const Standard_Real aSqDist = BRepMesh_GeomTool::SquareDeflectionOfSegment(
99         myPrevControlPnt, myCurrControlPnt, aMidPnt);
100
101       Standard_Real anAngle = 0.0;
102       
103       if ((myPrevControlVec.SquareMagnitude() > Precision::SquareConfusion()) &&
104           (myCurrControlVec.SquareMagnitude() > Precision::SquareConfusion()))
105       {
106         anAngle = myPrevControlVec.Angle(myCurrControlVec);
107       }
108
109       const Standard_Real aSqMaxDeflection = myDFace->GetDeflection() *
110         myDFace->GetDeflection();
111
112       if (((aSqDist > aSqMaxDeflection) || (anAngle > myParameters.AngleInterior)) &&
113           aSqDist > myParameters.MinSize * myParameters.MinSize)
114       {
115         // insertion 
116         myControlParams->InsertBefore(theIndex, aMidParam);
117       }
118       else
119       {
120         // Here we should leave at least 3 parameters as far as
121         // we must have at least one parameter related to surface
122         // internals in order to prevent movement of triangle body
123         // outside the surface in case of highly curved ones, e.g.
124         // BSpline springs.
125         if (((aSqDist < aSqMaxDeflection) || (anAngle < myParameters.AngleInterior)) &&
126             myControlParams->Length() > 3 && theIndex < myControlParams->Length())
127         {
128           // Remove too dense points
129           const Standard_Real aTmpParam = myControlParams->Value(theIndex + 1);
130           if (checkParameterForDeflectionAndUpdateCache(aTmpParam))
131           {
132             ++aMoveSteps;
133           }
134         }
135
136         myPrevControlParam = myCurrControlParam;
137         myPrevControlPnt   = myCurrControlPnt;
138         myPrevControlVec   = myCurrControlVec;
139
140         ++aMoveSteps;
141       }
142
143       return aMoveSteps;
144     }
145
146     //! Checks whether the given param suits specified deflection. Updates cache.
147     Standard_Boolean checkParameterForDeflectionAndUpdateCache(const Standard_Real theParam)
148     {
149       gp_Pnt aTmpPnt;
150       gp_Vec aTmpVec;
151       myIso->D1(theParam, aTmpPnt, aTmpVec);
152
153       const Standard_Real aTmpMidParam = 0.5 * (myPrevControlParam + theParam);
154       const gp_Pnt        aTmpMidPnt = myIso->Value(aTmpMidParam);
155
156       // Lets check next parameter.
157       // If it also fits deflection, we can remove previous parameter.
158       const Standard_Real aSqDist = BRepMesh_GeomTool::SquareDeflectionOfSegment(
159         myPrevControlPnt, aTmpPnt, aTmpMidPnt);
160
161       if (aSqDist < myDFace->GetDeflection() * myDFace->GetDeflection())
162       {
163         // Lets check parameters for angular deflection.
164         if (myPrevControlVec.SquareMagnitude() < gp::Resolution() ||
165             aTmpVec.SquareMagnitude()          < gp::Resolution() ||
166             myPrevControlVec.Angle(aTmpVec)    < myParameters.AngleInterior)
167         {
168           // For current Iso line we can remove this parameter.
169           myControlParamsToRemove->Add(myCurrControlParam);
170           myCurrControlParam = theParam;
171           myCurrControlPnt   = aTmpPnt;
172           myCurrControlVec   = aTmpVec;
173           return Standard_True;
174         }
175         else
176         {
177           // We have found a place on the surface refusing 
178           // removement of this parameter.
179           myParamsForbiddenToRemove       ->Add(myCurrParam);
180           myControlParamsForbiddenToRemove->Add(myCurrControlParam);
181         }
182       }
183
184       return Standard_False;
185     }
186
187   private:
188
189     IMeshData::IFaceHandle                myDFace;
190     Handle(Geom_Surface)                  mySurface;
191     Standard_Boolean                      myIsoU;
192     Handle(IMeshData::SequenceOfReal)     myParams;
193     Handle(IMeshData::SequenceOfReal)     myControlParams;
194
195     Handle(IMeshData::MapOfReal)          myParamsForbiddenToRemove;
196     Handle(IMeshData::MapOfReal)          myControlParamsForbiddenToRemove;
197
198     Handle(NCollection_IncAllocator)      myAllocator;
199     Handle(IMeshData::MapOfReal)          myControlParamsToRemove;
200
201
202     IMeshTools_Parameters                 myParameters;
203     NCollection_Handle<GeomAdaptor_Curve> myIso;
204
205     Standard_Real                         myCurrParam;
206
207     Standard_Real                         myCurrControlParam;
208     gp_Pnt                                myCurrControlPnt;
209     gp_Vec                                myCurrControlVec;
210
211     Standard_Real                         myPrevControlParam;
212     gp_Pnt                                myPrevControlPnt;
213     gp_Vec                                myPrevControlVec;
214   };
215
216   //! Adds param to map if it fits specified range.
217   inline Standard_Boolean addParam(
218     const Standard_Real&                           theParam,
219     const std::pair<Standard_Real, Standard_Real>& theRange,
220     IMeshData::IMapOfReal&                         theParams)
221   {
222     if (theParam < theRange.first ||
223         theParam > theRange.second)
224     {
225       return Standard_False;
226     }
227
228     theParams.Add(theParam);
229     return Standard_True;
230   }
231
232   //! Initializes parameters map using CN intervals.
233   inline Standard_Boolean initParamsFromIntervals(
234     const TColStd_Array1OfReal&                    theIntervals,
235     const std::pair<Standard_Real, Standard_Real>& theRange,
236     const Standard_Boolean                         isSplitIntervals,
237     IMeshData::IMapOfReal&                         theParams)
238   {
239     Standard_Boolean isAdded = Standard_False;
240
241     for (Standard_Integer i = theIntervals.Lower(); i <= theIntervals.Upper(); ++i)
242     {
243       const Standard_Real aStartParam = theIntervals.Value(i);
244       if (addParam(aStartParam, theRange, theParams))
245       {
246         isAdded = Standard_True;
247       }
248
249       if (isSplitIntervals && i < theIntervals.Upper())
250       {
251         const Standard_Real aMidParam = (aStartParam + theIntervals.Value(i + 1)) / 2.;
252         if (addParam(aMidParam, theRange, theParams))
253         {
254           isAdded = Standard_True;
255         }
256       }
257     }
258
259     return isAdded;
260   }
261
262   //! Checks whether intervals should be split.
263   //! Returns true in case if it is impossible to compute normal 
264   //! directly on intervals, false is returned elsewhere.
265   Standard_Boolean toSplitIntervals (const Handle (Geom_Surface)&  theSurf,
266                                      const TColStd_Array1OfReal  (&theIntervals)[2])
267   {
268     Standard_Integer aIntervalU = theIntervals[0].Lower ();
269     for (; aIntervalU <= theIntervals[0].Upper (); ++aIntervalU)
270     {
271       const Standard_Real aParamU = theIntervals[0].Value(aIntervalU);
272       if (Precision::IsInfinite (aParamU))
273         continue;
274
275       Standard_Integer aIntervalV = theIntervals[1].Lower ();
276       for (; aIntervalV <= theIntervals[1].Upper (); ++aIntervalV)
277       {
278         gp_Dir aNorm;
279         const Standard_Real aParamV = theIntervals[1].Value(aIntervalV);
280         if (Precision::IsInfinite (aParamV))
281           continue;
282
283         if (GeomLib::NormEstim (theSurf, gp_Pnt2d (aParamU, aParamV), Precision::Confusion (), aNorm) != 0)
284         {
285           return Standard_True;
286         }
287         // TODO: do not split intervals if there is no normal in the middle of interval.
288       }
289     }
290
291     return Standard_False;
292   }
293 }
294
295 //=======================================================================
296 // Function: AdjustRange
297 // Purpose : 
298 //=======================================================================
299 void BRepMesh_NURBSRangeSplitter::AdjustRange()
300 {
301   BRepMesh_DefaultRangeSplitter::AdjustRange();
302   mySurfaceType = GetSurface()->GetType();
303
304   if (mySurfaceType == GeomAbs_BezierSurface)
305   {
306     const std::pair<Standard_Real, Standard_Real>& aRangeU = GetRangeU();
307     const std::pair<Standard_Real, Standard_Real>& aRangeV = GetRangeV();
308
309     myIsValid = !(aRangeU.first  < -0.5 ||
310                   aRangeU.second >  1.5 ||
311                   aRangeV.first  < -0.5 ||
312                   aRangeV.second >  1.5);
313   }
314 }
315
316 //=======================================================================
317 // Function: GenerateSurfaceNodes
318 // Purpose : 
319 //=======================================================================
320 Handle(IMeshData::ListOfPnt2d) BRepMesh_NURBSRangeSplitter::GenerateSurfaceNodes(
321   const IMeshTools_Parameters& theParameters) const
322 {
323   if (!initParameters())
324   {
325     return Handle(IMeshData::ListOfPnt2d)();
326   }
327
328   const std::pair<Standard_Real, Standard_Real>& aRangeU = GetRangeU();
329   const std::pair<Standard_Real, Standard_Real>& aRangeV = GetRangeV();
330   const std::pair<Standard_Real, Standard_Real>& aDelta  = GetDelta ();
331
332   const Standard_Real                 aDefFace = GetDFace()->GetDeflection();
333   const Handle(BRepAdaptor_HSurface)& gFace    = GetSurface();
334   Handle(Geom_Surface)                aSurface = gFace->ChangeSurface().Surface().Surface();
335
336   const Handle(NCollection_IncAllocator) aTmpAlloc =
337     new NCollection_IncAllocator(IMeshData::MEMORY_BLOCK_SIZE_HUGE);
338
339   const Handle(IMeshData::SequenceOfReal) aParams[2] = {
340     computeGrainAndFilterParameters(GetParametersU(), gFace->UResolution(aDefFace),
341       (aRangeU.second - aRangeU.first), aDelta.first,  theParameters, aTmpAlloc),
342
343     computeGrainAndFilterParameters(GetParametersV(), gFace->VResolution(aDefFace),
344       (aRangeV.second - aRangeV.first), aDelta.second, theParameters, aTmpAlloc)
345   };
346
347   // check intermediate isolines
348   Handle(IMeshData::MapOfReal) aFixedParams[2] = {
349     new IMeshData::MapOfReal(1, aTmpAlloc),
350     new IMeshData::MapOfReal(1, aTmpAlloc)
351   };
352
353   const Handle(IMeshData::MapOfReal) aParamsToRemove[2] = {
354     AnalyticalFilter(GetDFace(), GeomAbs_IsoV, aParams[1], aParams[0],
355       aFixedParams[1], aFixedParams[0]).GetControlParametersToRemove(theParameters),
356
357     AnalyticalFilter(GetDFace(), GeomAbs_IsoU, aParams[0], aParams[1],
358       aFixedParams[0], aFixedParams[1]).GetControlParametersToRemove(theParameters),
359   };
360
361   aParamsToRemove[0]->Subtract(*aFixedParams[0]);
362   aParamsToRemove[1]->Subtract(*aFixedParams[1]);
363
364   // insert nodes of the regular grid
365   Handle(IMeshData::ListOfPnt2d) aNodes = new IMeshData::ListOfPnt2d(
366     new NCollection_IncAllocator(IMeshData::MEMORY_BLOCK_SIZE_HUGE));
367
368   // insert nodes of the regular grid
369   for (Standard_Integer i = 1; i <= aParams[0]->Length(); ++i)
370   {
371     const Standard_Real aParam1 = aParams[0]->Value(i);
372     if (aParamsToRemove[0]->Contains(aParam1))
373     {
374       continue;
375     }
376
377     for (Standard_Integer j = 1; j <= aParams[1]->Length(); ++j)
378     {
379       const Standard_Real aParam2 = aParams[1]->Value(j);
380       if (aParamsToRemove[1]->Contains(aParam2))
381       {
382         continue;
383       }
384
385       aNodes->Append(gp_Pnt2d(aParam1, aParam2));
386     }
387   }
388
389   return aNodes;
390 }
391
392 //=======================================================================
393 // Function: initParameters
394 // Purpose : 
395 //=======================================================================
396 Standard_Boolean BRepMesh_NURBSRangeSplitter::initParameters() const
397 {
398   const Handle(BRepAdaptor_HSurface)& aSurface = GetSurface();
399
400   const GeomAbs_Shape aContinuity = GeomAbs_CN;
401   const std::pair<Standard_Integer, Standard_Integer> aIntervalsNb(
402     aSurface->NbUIntervals(aContinuity),
403     aSurface->NbVIntervals(aContinuity)
404   );
405
406   TColStd_Array1OfReal aIntervals[2] = {
407     TColStd_Array1OfReal(1, aIntervalsNb.first  + 1),
408     TColStd_Array1OfReal(1, aIntervalsNb.second + 1)
409   };
410
411   aSurface->UIntervals(aIntervals[0], aContinuity);
412   aSurface->VIntervals(aIntervals[1], aContinuity);
413
414   const Standard_Boolean isSplitIntervals = toSplitIntervals (
415     aSurface->ChangeSurface().Surface().Surface(), aIntervals);
416
417   if (!initParamsFromIntervals(aIntervals[0], GetRangeU(), isSplitIntervals,
418                                const_cast<IMeshData::IMapOfReal&>(GetParametersU())))
419   {
420     //if (!grabParamsOfEdges (Edge_Frontier, Param_U))
421     {
422       return Standard_False;
423     }
424   }
425
426   if (!initParamsFromIntervals(aIntervals[1], GetRangeV(), isSplitIntervals,
427                                const_cast<IMeshData::IMapOfReal&>(GetParametersV())))
428   {
429     //if (!grabParamsOfEdges (Edge_Frontier, Param_V))
430     {
431       return Standard_False;
432     }
433   }
434
435   return grabParamsOfEdges(Edge_Internal, Param_U | Param_V);
436 }
437
438 //=======================================================================
439 //function : grabParamsOfInternalEdges
440 //purpose  : 
441 //=======================================================================
442 Standard_Boolean BRepMesh_NURBSRangeSplitter::grabParamsOfEdges (
443   const EdgeType         theEdgeType,
444   const Standard_Integer theParamDimensionFlag) const
445 {
446   if ((theParamDimensionFlag & (Param_U | Param_V)) == 0)
447   {
448     return Standard_False;
449   }
450
451   const IMeshData::IFaceHandle& aDFace = GetDFace ();
452   for (Standard_Integer aWireIt = 0; aWireIt < aDFace->WiresNb (); ++aWireIt)
453   {
454     const IMeshData::IWireHandle& aDWire = aDFace->GetWire (aWireIt);
455     for (Standard_Integer aEdgeIt = 0; aEdgeIt < aDWire->EdgesNb (); ++aEdgeIt)
456     {
457       const IMeshData::IEdgePtr& aDEdge = aDWire->GetEdge (aEdgeIt);
458       for (Standard_Integer aPCurveIt = 0; aPCurveIt < aDEdge->PCurvesNb (); ++aPCurveIt)
459       {
460         const IMeshData::IPCurveHandle& aDPCurve = aDEdge->GetPCurve (aPCurveIt);
461         if (aDPCurve->GetFace () == aDFace)
462         {
463           if (theEdgeType == Edge_Internal && !aDPCurve->IsInternal ())
464           {
465             continue;
466           }
467
468           for (Standard_Integer aPointIt = 0; aPointIt < aDPCurve->ParametersNb (); ++aPointIt)
469           {
470             const gp_Pnt2d& aPnt2d = aDPCurve->GetPoint (aPointIt);
471             if (theParamDimensionFlag & Param_U)
472             {
473               const_cast<IMeshData::IMapOfReal&>(GetParametersU ()).Add (aPnt2d.X ());
474             }
475
476             if (theParamDimensionFlag & Param_V)
477             {
478               const_cast<IMeshData::IMapOfReal&>(GetParametersV ()).Add (aPnt2d.Y ());
479             }
480           }
481         }
482       }
483     }
484   }
485
486   return Standard_True;
487 }
488
489 //=======================================================================
490 //function : computeGrainAndFilterParameters
491 //purpose  : 
492 //=======================================================================
493 Handle(IMeshData::SequenceOfReal) BRepMesh_NURBSRangeSplitter::computeGrainAndFilterParameters(
494   const IMeshData::IMapOfReal&            theSourceParams,
495   const Standard_Real                     theTol2d,
496   const Standard_Real                     theRangeDiff,
497   const Standard_Real                     theDelta,
498   const IMeshTools_Parameters&            theParameters,
499   const Handle(NCollection_IncAllocator)& theAllocator) const
500 {
501   // Sort and filter sequence of parameters
502   Standard_Real aMinDiff = Precision::PConfusion();
503   if (theDelta < 1.)
504   {
505     aMinDiff /= theDelta;
506   }
507
508   aMinDiff = Max(theParameters.MinSize, aMinDiff);
509
510   const Standard_Real aDiffMaxLim = 0.1 * theRangeDiff;
511   const Standard_Real aDiffMinLim = Max(0.005 * theRangeDiff,
512                                         2. * theTol2d);
513   const Standard_Real aDiff = Max(theParameters.MinSize,
514                                   Min(aDiffMaxLim, aDiffMinLim));
515   return filterParameters(theSourceParams, aMinDiff, aDiff, theAllocator);
516 }
517
518 //=======================================================================
519 //function : filterParameters
520 //purpose  : 
521 //=======================================================================
522 Handle(IMeshData::SequenceOfReal) BRepMesh_NURBSRangeSplitter::filterParameters(
523   const IMeshData::IMapOfReal&            theParams,
524   const Standard_Real                     theMinDist,
525   const Standard_Real                     theFilterDist,
526   const Handle(NCollection_IncAllocator)& theAllocator) const
527 {
528   Handle(IMeshData::SequenceOfReal) aResult = new IMeshData::SequenceOfReal(theAllocator);
529
530   // Sort sequence of parameters
531   const Standard_Integer anInitLen = theParams.Extent();
532
533   if (anInitLen < 1)
534   {
535     return aResult;
536   }
537
538   TColStd_Array1OfReal aParamArray(1, anInitLen);
539   Standard_Integer j;
540   for (j = 1; j <= anInitLen; j++)
541     aParamArray(j) = theParams(j);
542
543   std::sort(aParamArray.begin(), aParamArray.end());
544
545   // mandatory pre-filtering using the first (minimal) filter value
546   Standard_Integer aParamLength = 1;
547   for (j = 2; j <= anInitLen; j++)
548   {
549     if ((aParamArray(j) - aParamArray(aParamLength)) > theMinDist)
550     {
551       if (++aParamLength < j)
552         aParamArray(aParamLength) = aParamArray(j);
553     }
554   }
555
556   //perform filtering on series
557   Standard_Real aLastAdded, aLastCandidate;
558   Standard_Boolean isCandidateDefined = Standard_False;
559   aLastAdded = aParamArray(1);
560   aLastCandidate = aLastAdded;
561   aResult->Append(aLastAdded);
562
563   for (j = 2; j < aParamLength; j++)
564   {
565     Standard_Real aVal = aParamArray(j);
566     if (aVal - aLastAdded > theFilterDist)
567     {
568       //adds the parameter
569       if (isCandidateDefined)
570       {
571         aLastAdded = aLastCandidate;
572         isCandidateDefined = Standard_False;
573         j--;
574       }
575       else
576       {
577         aLastAdded = aVal;
578       }
579       aResult->Append(aLastAdded);
580       continue;
581     }
582
583     aLastCandidate = aVal;
584     isCandidateDefined = Standard_True;
585   }
586   aResult->Append(aParamArray(aParamLength));
587
588   return aResult;
589 }