0030394: Modeling Algorithms - Empty result of offset operation (mode "Complete"...
[occt.git] / src / BVH / BVH_DistanceField.lxx
CommitLineData
df932fdf 1// Created on: 2014-09-06
2// Created by: Denis BOGOLEPOV
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
df932fdf 16#include <BVH_Triangulation.hxx>
17
18#ifdef HAVE_TBB
19 // On Windows, function TryEnterCriticalSection has appeared in Windows NT
20 // and is surrounded by #ifdef in MS VC++ 7.1 headers.
21 // Thus to use it we need to define appropriate macro saying that we will
22 // run on Windows NT 4.0 at least
23 #if defined(_WIN32) && !defined(_WIN32_WINNT)
24 #define _WIN32_WINNT 0x0501
25 #endif
26
27 #include <tbb/tbb.h>
28#endif
29
30// =======================================================================
31// function : BVH_DistanceField
32// purpose :
33// =======================================================================
34template<class T, int N>
35BVH_DistanceField<T, N>::BVH_DistanceField (const Standard_Integer theMaximumSize,
36 const Standard_Boolean theComputeSign)
37: myMaximumSize (theMaximumSize),
38 myComputeSign (theComputeSign)
39{
40 Standard_STATIC_ASSERT (N == 3 || N == 4);
41
42 myVoxelData = new T[myMaximumSize * myMaximumSize * myMaximumSize];
43}
44
45// =======================================================================
46// function : ~BVH_DistanceField
47// purpose :
48// =======================================================================
49template<class T, int N>
50BVH_DistanceField<T, N>::~BVH_DistanceField()
51{
52 delete [] myVoxelData;
53}
54
55#if defined (_WIN32) && defined (max)
56 #undef max
57#endif
58
59#include <limits>
60
61#define BVH_DOT3(A, B) (A.x() * B.x() + A.y() * B.y() + A.z() * B.z())
62
63namespace BVH
64{
65 //=======================================================================
66 //function : DistanceToBox
67 //purpose : Computes squared distance from point to box
68 //=======================================================================
69 template<class T, int N>
70 T DistanceToBox (const typename VectorType<T, N>::Type& thePnt,
71 const typename VectorType<T, N>::Type& theMin,
72 const typename VectorType<T, N>::Type& theMax)
73 {
74 Standard_STATIC_ASSERT (N == 3 || N == 4);
75
76 T aNearestX = Min (Max (thePnt.x(), theMin.x()), theMax.x());
77 T aNearestY = Min (Max (thePnt.y(), theMin.y()), theMax.y());
78 T aNearestZ = Min (Max (thePnt.z(), theMin.z()), theMax.z());
79
80 if (aNearestX == thePnt.x()
81 && aNearestY == thePnt.y()
82 && aNearestZ == thePnt.z())
83 {
84 return static_cast<T> (0);
85 }
86
87 aNearestX -= thePnt.x();
88 aNearestY -= thePnt.y();
89 aNearestZ -= thePnt.z();
90
91 return aNearestX * aNearestX +
92 aNearestY * aNearestY +
93 aNearestZ * aNearestZ;
94 }
95
96 //=======================================================================
97 //function : DirectionToNearestPoint
98 //purpose : Computes squared distance from point to triangle
99 // ======================================================================
100 template<class T, int N>
101 typename VectorType<T, N>::Type DirectionToNearestPoint (
102 const typename VectorType<T, N>::Type& thePoint,
103 const typename VectorType<T, N>::Type& theVertA,
104 const typename VectorType<T, N>::Type& theVertB,
105 const typename VectorType<T, N>::Type& theVertC)
106 {
107 Standard_STATIC_ASSERT (N == 3 || N == 4);
108
109 const typename VectorType<T, N>::Type aAB = theVertB - theVertA;
110 const typename VectorType<T, N>::Type aAC = theVertC - theVertA;
111 const typename VectorType<T, N>::Type aAP = thePoint - theVertA;
112
113 const T aABdotAP = BVH_DOT3 (aAB, aAP);
114 const T aACdotAP = BVH_DOT3 (aAC, aAP);
115
116 if (aABdotAP <= static_cast<T> (0) && aACdotAP <= static_cast<T> (0))
117 {
118 return aAP;
119 }
120
121 const typename VectorType<T, N>::Type aBC = theVertC - theVertB;
122 const typename VectorType<T, N>::Type aBP = thePoint - theVertB;
123
124 const T aBAdotBP = -BVH_DOT3 (aAB, aBP);
125 const T aBCdotBP = BVH_DOT3 (aBC, aBP);
126
127 if (aBAdotBP <= static_cast<T> (0) && aBCdotBP <= static_cast<T> (0))
128 {
129 return aBP;
130 }
131
132 const typename VectorType<T, N>::Type aCP = thePoint - theVertC;
133
134 const T aCBdotCP = -BVH_DOT3 (aBC, aCP);
135 const T aCAdotCP = -BVH_DOT3 (aAC, aCP);
136
137 if (aCAdotCP <= static_cast<T> (0) && aCBdotCP <= static_cast<T> (0))
138 {
139 return aCP;
140 }
141
142 const T aACdotBP = BVH_DOT3 (aAC, aBP);
143
144 const T aVC = aABdotAP * aACdotBP + aBAdotBP * aACdotAP;
145
146 if (aVC <= static_cast<T> (0) && aABdotAP >= static_cast<T> (0) && aBAdotBP >= static_cast<T> (0))
147 {
148 return aAP - aAB * (aABdotAP / (aABdotAP + aBAdotBP));
149 }
150
151 const T aABdotCP = BVH_DOT3 (aAB, aCP);
152
153 const T aVA = aBAdotBP * aCAdotCP - aABdotCP * aACdotBP;
154
155 if (aVA <= static_cast<T> (0) && aBCdotBP >= static_cast<T> (0) && aCBdotCP >= static_cast<T> (0))
156 {
157 return aBP - aBC * (aBCdotBP / (aBCdotBP + aCBdotCP));
158 }
159
160 const T aVB = aABdotCP * aACdotAP + aABdotAP * aCAdotCP;
161
162 if (aVB <= static_cast<T> (0) && aACdotAP >= static_cast<T> (0) && aCAdotCP >= static_cast<T> (0))
163 {
164 return aAP - aAC * (aACdotAP / (aACdotAP + aCAdotCP));
165 }
166
167 const T aNorm = static_cast<T> (1.0) / (aVA + aVB + aVC);
168
169 const T aU = aVA * aNorm;
170 const T aV = aVB * aNorm;
171
172 return thePoint - (theVertA * aU + theVertB * aV + theVertC * (static_cast<T> (1.0) - aU - aV));
173 }
174
175 //=======================================================================
176 //function : SquareDistanceToObject
177 //purpose : Computes squared distance from point to BVH triangulation
178 //=======================================================================
179 template<class T, int N>
180 T SquareDistanceToObject (BVH_Object<T, N>* theObject,
181 const typename VectorType<T, N>::Type& thePnt, Standard_Boolean& theIsOutside)
182 {
183 Standard_STATIC_ASSERT (N == 3 || N == 4);
184
185 T aMinDistance = std::numeric_limits<T>::max();
186
187 BVH_Triangulation<T, N>* aTriangulation =
188 dynamic_cast<BVH_Triangulation<T, N>*> (theObject);
189
190 Standard_ASSERT_RETURN (aTriangulation != NULL,
191 "Error: Unsupported BVH object (non triangulation)", aMinDistance);
192
f5b72419 193 const opencascade::handle<BVH_Tree<T, N> >& aBVH = aTriangulation->BVH();
df932fdf 194 if (aBVH.IsNull())
195 {
196 return Standard_False;
197 }
198
f5b72419 199 std::pair<Standard_Integer, T> aStack[BVH_Constants_MaxTreeDepth];
df932fdf 200
201 Standard_Integer aHead = -1;
202 Standard_Integer aNode = 0; // root node
203
204 for (;;)
205 {
206 BVH_Vec4i aData = aBVH->NodeInfoBuffer()[aNode];
207
208 if (aData.x() == 0) // if inner node
209 {
210 const T aDistToLft = DistanceToBox<T, N> (thePnt,
211 aBVH->MinPoint (aData.y()),
212 aBVH->MaxPoint (aData.y()));
213
214 const T aDistToRgh = DistanceToBox<T, N> (thePnt,
215 aBVH->MinPoint (aData.z()),
216 aBVH->MaxPoint (aData.z()));
217
218 const Standard_Boolean aHitLft = aDistToLft <= aMinDistance;
219 const Standard_Boolean aHitRgh = aDistToRgh <= aMinDistance;
220
221 if (aHitLft & aHitRgh)
222 {
223 aNode = (aDistToLft < aDistToRgh) ? aData.y() : aData.z();
224
225 aStack[++aHead] = std::pair<Standard_Integer, T> (
226 aDistToLft < aDistToRgh ? aData.z() : aData.y(), Max (aDistToLft, aDistToRgh));
227 }
228 else
229 {
230 if (aHitLft | aHitRgh)
231 {
232 aNode = aHitLft ? aData.y() : aData.z();
233 }
234 else
235 {
236 if (aHead < 0)
237 return aMinDistance;
238
239 std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
240
241 while (anInfo.second > aMinDistance)
242 {
243 if (aHead < 0)
244 return aMinDistance;
245
246 anInfo = aStack[aHead--];
247 }
248
249 aNode = anInfo.first;
250 }
251 }
252 }
253 else // if leaf node
254 {
255 for (Standard_Integer aTrgIdx = aData.y(); aTrgIdx <= aData.z(); ++aTrgIdx)
256 {
257 const BVH_Vec4i aTriangle = aTriangulation->Elements[aTrgIdx];
258
259 const typename VectorType<T, N>::Type aVertex0 = aTriangulation->Vertices[aTriangle.x()];
260 const typename VectorType<T, N>::Type aVertex1 = aTriangulation->Vertices[aTriangle.y()];
261 const typename VectorType<T, N>::Type aVertex2 = aTriangulation->Vertices[aTriangle.z()];
262
263 const typename VectorType<T, N>::Type aDirection =
264 DirectionToNearestPoint<T, N> (thePnt, aVertex0, aVertex1, aVertex2);
265
266 const T aDistance = BVH_DOT3 (aDirection, aDirection);
267
268 if (aDistance < aMinDistance)
269 {
270 aMinDistance = aDistance;
271
272 typename VectorType<T, N>::Type aTrgEdges[] = { aVertex1 - aVertex0,
273 aVertex2 - aVertex0 };
274
275 typename VectorType<T, N>::Type aTrgNormal;
276
277 aTrgNormal.x() = aTrgEdges[0].y() * aTrgEdges[1].z() - aTrgEdges[0].z() * aTrgEdges[1].y();
278 aTrgNormal.y() = aTrgEdges[0].z() * aTrgEdges[1].x() - aTrgEdges[0].x() * aTrgEdges[1].z();
279 aTrgNormal.z() = aTrgEdges[0].x() * aTrgEdges[1].y() - aTrgEdges[0].y() * aTrgEdges[1].x();
280
281 theIsOutside = BVH_DOT3 (aTrgNormal, aDirection) > 0;
282 }
283 }
284
285 if (aHead < 0)
286 return aMinDistance;
287
288 std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
289
290 while (anInfo.second > aMinDistance)
291 {
292 if (aHead < 0)
293 return aMinDistance;
294
295 anInfo = aStack[aHead--];
296 }
297
298 aNode = anInfo.first;
299 }
300 }
301 }
302
303 //=======================================================================
304 //function : SquareDistanceToGeomerty
305 //purpose : Computes squared distance from point to BVH geometry
306 //=======================================================================
307 template<class T, int N>
308 T SquareDistanceToGeomerty (BVH_Geometry<T, N>& theGeometry,
309 const typename VectorType<T, N>::Type& thePnt, Standard_Boolean& theIsOutside)
310 {
311 Standard_STATIC_ASSERT (N == 3 || N == 4);
312
f2474958 313 const BVH_Tree<T, N, BVH_BinaryTree>* aBVH = theGeometry.BVH().get();
df932fdf 314
f2474958 315 if (aBVH == NULL)
df932fdf 316 {
317 return Standard_False;
318 }
319
f5b72419 320 std::pair<Standard_Integer, T> aStack[BVH_Constants_MaxTreeDepth];
df932fdf 321
322 Standard_Integer aHead = -1;
323 Standard_Integer aNode = 0; // root node
324
325 T aMinDistance = std::numeric_limits<T>::max();
326
327 for (;;)
328 {
329 BVH_Vec4i aData = aBVH->NodeInfoBuffer()[aNode];
330
331 if (aData.x() == 0) // if inner node
332 {
333 const T aDistToLft = DistanceToBox<T, N> (thePnt,
334 aBVH->MinPoint (aData.y()),
335 aBVH->MaxPoint (aData.y()));
336
337 const T aDistToRgh = DistanceToBox<T, N> (thePnt,
338 aBVH->MinPoint (aData.z()),
339 aBVH->MaxPoint (aData.z()));
340
341 const Standard_Boolean aHitLft = aDistToLft <= aMinDistance;
342 const Standard_Boolean aHitRgh = aDistToRgh <= aMinDistance;
343
344 if (aHitLft & aHitRgh)
345 {
346 aNode = (aDistToLft < aDistToRgh) ? aData.y() : aData.z();
347
348 aStack[++aHead] = std::pair<Standard_Integer, T> (
349 aDistToLft < aDistToRgh ? aData.z() : aData.y(), Max (aDistToLft, aDistToRgh));
350 }
351 else
352 {
353 if (aHitLft | aHitRgh)
354 {
355 aNode = aHitLft ? aData.y() : aData.z();
356 }
357 else
358 {
359 if (aHead < 0)
360 return aMinDistance;
361
362 std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
363
364 while (anInfo.second > aMinDistance)
365 {
366 if (aHead < 0)
367 return aMinDistance;
368
369 anInfo = aStack[aHead--];
370 }
371
372 aNode = anInfo.first;
373 }
374 }
375 }
376 else // if leaf node
377 {
378 Standard_Boolean isOutside = Standard_True;
379
380 const T aDistance = SquareDistanceToObject (
381 theGeometry.Objects()(aNode).operator->(), thePnt, isOutside);
382
383 if (aDistance < aMinDistance)
384 {
385 aMinDistance = aDistance;
386 theIsOutside = isOutside;
387 }
388
389 if (aHead < 0)
390 return aMinDistance;
391
392 std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
393
394 while (anInfo.second > aMinDistance)
395 {
396 if (aHead < 0)
397 return aMinDistance;
398
399 anInfo = aStack[aHead--];
400 }
401
402 aNode = anInfo.first;
403 }
404 }
405 }
406}
407
408#undef BVH_DOT3
409
410#ifdef HAVE_TBB
411
412//! Tool object for parallel construction of distance field (uses Intel TBB).
413template<class T, int N>
414class BVH_ParallelDistanceFieldBuilder
415{
416private:
417
418 //! Input BVH geometry.
419 BVH_Geometry<T, N>* myGeometry;
420
421 //! Output distance field.
422 BVH_DistanceField<T, N>* myOutField;
423
424public:
425
426 BVH_ParallelDistanceFieldBuilder (BVH_DistanceField<T, N>* theOutField, BVH_Geometry<T, N>* theGeometry)
427 : myGeometry (theGeometry),
428 myOutField (theOutField)
429 {
430 //
431 }
432
433 void operator() (const tbb::blocked_range<Standard_Integer>& theRange) const
434 {
435 myOutField->BuildSlices (*myGeometry, theRange.begin(), theRange.end());
436 }
437};
438
439#endif
440
441// =======================================================================
442// function : BuildSlices
443// purpose : Performs building of distance field for the given Z slices
444// =======================================================================
445template<class T, int N>
446void BVH_DistanceField<T, N>::BuildSlices (BVH_Geometry<T, N>& theGeometry,
447 const Standard_Integer theStartSlice, const Standard_Integer theFinalSlice)
448{
449 for (Standard_Integer aZ = theStartSlice; aZ < theFinalSlice; ++aZ)
450 {
451 for (Standard_Integer aY = 0; aY < myDimensionY; ++aY)
452 {
453 for (Standard_Integer aX = 0; aX < myDimensionX; ++aX)
454 {
455 BVH_VecNt aCenter;
456
457 aCenter.x() = myCornerMin.x() + myVoxelSize.x() * (aX + static_cast<T> (0.5));
458 aCenter.y() = myCornerMin.y() + myVoxelSize.y() * (aY + static_cast<T> (0.5));
459 aCenter.z() = myCornerMin.z() + myVoxelSize.z() * (aZ + static_cast<T> (0.5));
460
461 Standard_Boolean isOutside = Standard_True;
462
463 const T aDistance = sqrt (
464 BVH::SquareDistanceToGeomerty<T, N> (theGeometry, aCenter, isOutside));
465
466 Voxel (aX, aY, aZ) = (!myComputeSign || isOutside) ? aDistance : -aDistance;
467 }
468 }
469 }
470}
471
472// =======================================================================
473// function : Build
474// purpose : Builds 3D distance field from BVH geometry
475// =======================================================================
476template<class T, int N>
477Standard_Boolean BVH_DistanceField<T, N>::Build (BVH_Geometry<T, N>& theGeometry)
478{
479 if (theGeometry.Size() == 0)
480 {
481 return Standard_False;
482 }
483
484 const BVH_VecNt aGlobalBoxSize = theGeometry.Box().Size();
485
486 const T aMaxBoxSide = Max (Max (aGlobalBoxSize.x(), aGlobalBoxSize.y()), aGlobalBoxSize.z());
487
488 myDimensionX = static_cast<Standard_Integer> (myMaximumSize * aGlobalBoxSize.x() / aMaxBoxSide);
489 myDimensionY = static_cast<Standard_Integer> (myMaximumSize * aGlobalBoxSize.y() / aMaxBoxSide);
490 myDimensionZ = static_cast<Standard_Integer> (myMaximumSize * aGlobalBoxSize.z() / aMaxBoxSide);
491
492 myDimensionX = Min (myMaximumSize, Max (myDimensionX, 16));
493 myDimensionY = Min (myMaximumSize, Max (myDimensionY, 16));
494 myDimensionZ = Min (myMaximumSize, Max (myDimensionZ, 16));
495
496 const BVH_VecNt aGlobalBoxMin = theGeometry.Box().CornerMin();
497 const BVH_VecNt aGlobalBoxMax = theGeometry.Box().CornerMax();
498
499 const Standard_Integer aVoxelOffset = 2;
500
501 myCornerMin.x() = aGlobalBoxMin.x() - aVoxelOffset * aGlobalBoxSize.x() / (myDimensionX - 2 * aVoxelOffset);
502 myCornerMin.y() = aGlobalBoxMin.y() - aVoxelOffset * aGlobalBoxSize.y() / (myDimensionY - 2 * aVoxelOffset);
503 myCornerMin.z() = aGlobalBoxMin.z() - aVoxelOffset * aGlobalBoxSize.z() / (myDimensionZ - 2 * aVoxelOffset);
504
505 myCornerMax.x() = aGlobalBoxMax.x() + aVoxelOffset * aGlobalBoxSize.x() / (myDimensionX - 2 * aVoxelOffset);
506 myCornerMax.y() = aGlobalBoxMax.y() + aVoxelOffset * aGlobalBoxSize.y() / (myDimensionY - 2 * aVoxelOffset);
507 myCornerMax.z() = aGlobalBoxMax.z() + aVoxelOffset * aGlobalBoxSize.z() / (myDimensionZ - 2 * aVoxelOffset);
508
509 myVoxelSize.x() = (myCornerMax.x() - myCornerMin.x()) / myDimensionX;
510 myVoxelSize.y() = (myCornerMax.y() - myCornerMin.y()) / myDimensionY;
511 myVoxelSize.z() = (myCornerMax.z() - myCornerMin.z()) / myDimensionZ;
512
513#ifdef HAVE_TBB
514
515 tbb::parallel_for (tbb::blocked_range<Standard_Integer> (0, myDimensionZ),
516 BVH_ParallelDistanceFieldBuilder<T, N> (this, &theGeometry));
517
518#else
519
520 BuildSlices (theGeometry, 0, myDimensionZ);
521
522#endif
523
524 return Standard_True;
525}