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