Warning, /include/opencascade/BVH_DistanceField.lxx is written in an unsupported language. File is not indexed.
0001 // Created on: 2014-09-06
0002 // Created by: Denis BOGOLEPOV
0003 // Copyright (c) 2013-2014 OPEN CASCADE SAS
0004 //
0005 // This file is part of Open CASCADE Technology software library.
0006 //
0007 // This library is free software; you can redistribute it and/or modify it under
0008 // the terms of the GNU Lesser General Public License version 2.1 as published
0009 // by the Free Software Foundation, with special exception defined in the file
0010 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0011 // distribution for complete text of the license and disclaimer of any warranty.
0012 //
0013 // Alternatively, this file may be used under the terms of Open CASCADE
0014 // commercial license or contractual agreement.
0015
0016 #include <BVH_Triangulation.hxx>
0017 #include <OSD_Parallel.hxx>
0018 #include <BVH_Distance.hxx>
0019
0020 // =======================================================================
0021 // function : BVH_DistanceField
0022 // purpose :
0023 // =======================================================================
0024 template<class T, int N>
0025 BVH_DistanceField<T, N>::BVH_DistanceField (const Standard_Integer theMaximumSize,
0026 const Standard_Boolean theComputeSign)
0027 : myDimensionX(0),
0028 myDimensionY(0),
0029 myDimensionZ(0),
0030 myMaximumSize (theMaximumSize),
0031 myComputeSign (theComputeSign),
0032 myIsParallel (Standard_False)
0033 {
0034 Standard_STATIC_ASSERT (N == 3 || N == 4);
0035
0036 myVoxelData = new T[myMaximumSize * myMaximumSize * myMaximumSize];
0037 }
0038
0039 // =======================================================================
0040 // function : ~BVH_DistanceField
0041 // purpose :
0042 // =======================================================================
0043 template<class T, int N>
0044 BVH_DistanceField<T, N>::~BVH_DistanceField()
0045 {
0046 delete [] myVoxelData;
0047 }
0048
0049 #if defined (_WIN32) && defined (max)
0050 #undef max
0051 #endif
0052
0053 #include <limits>
0054
0055 #define BVH_DOT3(A, B) (A.x() * B.x() + A.y() * B.y() + A.z() * B.z())
0056
0057 namespace BVH
0058 {
0059 //=======================================================================
0060 //function : DistanceToBox
0061 //purpose : Computes squared distance from point to box
0062 //=======================================================================
0063 template<class T, int N>
0064 T DistanceToBox (const typename VectorType<T, N>::Type& thePnt,
0065 const typename VectorType<T, N>::Type& theMin,
0066 const typename VectorType<T, N>::Type& theMax)
0067 {
0068 Standard_STATIC_ASSERT (N == 3 || N == 4);
0069
0070 T aNearestX = Min (Max (thePnt.x(), theMin.x()), theMax.x());
0071 T aNearestY = Min (Max (thePnt.y(), theMin.y()), theMax.y());
0072 T aNearestZ = Min (Max (thePnt.z(), theMin.z()), theMax.z());
0073
0074 if (aNearestX == thePnt.x()
0075 && aNearestY == thePnt.y()
0076 && aNearestZ == thePnt.z())
0077 {
0078 return static_cast<T> (0);
0079 }
0080
0081 aNearestX -= thePnt.x();
0082 aNearestY -= thePnt.y();
0083 aNearestZ -= thePnt.z();
0084
0085 return aNearestX * aNearestX +
0086 aNearestY * aNearestY +
0087 aNearestZ * aNearestZ;
0088 }
0089
0090 //=======================================================================
0091 //function : DirectionToNearestPoint
0092 //purpose : Computes squared distance from point to triangle
0093 // ======================================================================
0094 template<class T, int N>
0095 typename VectorType<T, N>::Type DirectionToNearestPoint (
0096 const typename VectorType<T, N>::Type& thePoint,
0097 const typename VectorType<T, N>::Type& theVertA,
0098 const typename VectorType<T, N>::Type& theVertB,
0099 const typename VectorType<T, N>::Type& theVertC)
0100 {
0101 Standard_STATIC_ASSERT (N == 3 || N == 4);
0102
0103 const typename VectorType<T, N>::Type aAB = theVertB - theVertA;
0104 const typename VectorType<T, N>::Type aAC = theVertC - theVertA;
0105 const typename VectorType<T, N>::Type aAP = thePoint - theVertA;
0106
0107 const T aABdotAP = BVH_DOT3 (aAB, aAP);
0108 const T aACdotAP = BVH_DOT3 (aAC, aAP);
0109
0110 if (aABdotAP <= static_cast<T> (0) && aACdotAP <= static_cast<T> (0))
0111 {
0112 return aAP;
0113 }
0114
0115 const typename VectorType<T, N>::Type aBC = theVertC - theVertB;
0116 const typename VectorType<T, N>::Type aBP = thePoint - theVertB;
0117
0118 const T aBAdotBP = -BVH_DOT3 (aAB, aBP);
0119 const T aBCdotBP = BVH_DOT3 (aBC, aBP);
0120
0121 if (aBAdotBP <= static_cast<T> (0) && aBCdotBP <= static_cast<T> (0))
0122 {
0123 return aBP;
0124 }
0125
0126 const typename VectorType<T, N>::Type aCP = thePoint - theVertC;
0127
0128 const T aCBdotCP = -BVH_DOT3 (aBC, aCP);
0129 const T aCAdotCP = -BVH_DOT3 (aAC, aCP);
0130
0131 if (aCAdotCP <= static_cast<T> (0) && aCBdotCP <= static_cast<T> (0))
0132 {
0133 return aCP;
0134 }
0135
0136 const T aACdotBP = BVH_DOT3 (aAC, aBP);
0137
0138 const T aVC = aABdotAP * aACdotBP + aBAdotBP * aACdotAP;
0139
0140 if (aVC <= static_cast<T> (0) && aABdotAP >= static_cast<T> (0) && aBAdotBP >= static_cast<T> (0))
0141 {
0142 return aAP - aAB * (aABdotAP / (aABdotAP + aBAdotBP));
0143 }
0144
0145 const T aABdotCP = BVH_DOT3 (aAB, aCP);
0146
0147 const T aVA = aBAdotBP * aCAdotCP - aABdotCP * aACdotBP;
0148
0149 if (aVA <= static_cast<T> (0) && aBCdotBP >= static_cast<T> (0) && aCBdotCP >= static_cast<T> (0))
0150 {
0151 return aBP - aBC * (aBCdotBP / (aBCdotBP + aCBdotCP));
0152 }
0153
0154 const T aVB = aABdotCP * aACdotAP + aABdotAP * aCAdotCP;
0155
0156 if (aVB <= static_cast<T> (0) && aACdotAP >= static_cast<T> (0) && aCAdotCP >= static_cast<T> (0))
0157 {
0158 return aAP - aAC * (aACdotAP / (aACdotAP + aCAdotCP));
0159 }
0160
0161 const T aNorm = static_cast<T> (1.0) / (aVA + aVB + aVC);
0162
0163 const T aU = aVA * aNorm;
0164 const T aV = aVB * aNorm;
0165
0166 return thePoint - (theVertA * aU + theVertB * aV + theVertC * (static_cast<T> (1.0) - aU - aV));
0167 }
0168
0169 //=======================================================================
0170 //function : SquareDistanceToPoint
0171 //purpose : Abstract class to compute squared distance from point to BVH tree
0172 //=======================================================================
0173 template<class T, int N, class BVHSetType>
0174 class SquareDistanceToPoint : public BVH_Distance<T,N,typename VectorType<T, N>::Type, BVHSetType>
0175 {
0176 public:
0177 typedef typename VectorType<T, N>::Type BVH_VecNt;
0178
0179 public:
0180 SquareDistanceToPoint()
0181 : BVH_Distance<T, N, BVH_VecNt, BVHSetType>(),
0182 myIsOutside (Standard_True)
0183 {}
0184
0185 public:
0186
0187 //! IsOutside
0188 Standard_Boolean IsOutside() const { return myIsOutside; }
0189
0190 public:
0191
0192 //! Defines the rules for node rejection
0193 virtual Standard_Boolean RejectNode (const BVH_VecNt& theCMin,
0194 const BVH_VecNt& theCMax,
0195 T& theMetric) const Standard_OVERRIDE
0196 {
0197 theMetric = DistanceToBox<T, N> (this->myObject, theCMin, theCMax);
0198 return theMetric > this->myDistance;
0199 }
0200
0201 public:
0202
0203 //! Redefine the Stop to never stop the selection
0204 virtual Standard_Boolean Stop() const Standard_OVERRIDE
0205 {
0206 return Standard_False;
0207 }
0208
0209 protected:
0210
0211 Standard_Boolean myIsOutside;
0212 };
0213
0214 //=======================================================================
0215 //function : PointTriangulationSquareDistance
0216 //purpose : Computes squared distance from point to BVH triangulation
0217 //=======================================================================
0218 template<class T, int N>
0219 class PointTriangulationSquareDistance: public SquareDistanceToPoint <T,N, BVH_Triangulation<T, N>>
0220 {
0221 public:
0222 typedef typename VectorType<T, N>::Type BVH_VecNt;
0223
0224 public:
0225 //! Constructor
0226 PointTriangulationSquareDistance()
0227 : SquareDistanceToPoint<T, N, BVH_Triangulation<T, N>>()
0228 {
0229 }
0230
0231 public:
0232 // Accepting the element
0233 virtual Standard_Boolean Accept (const Standard_Integer theIndex,
0234 const T&) Standard_OVERRIDE
0235 {
0236 const BVH_Vec4i aTriangle = this->myBVHSet->Elements[theIndex];
0237
0238 const BVH_VecNt aVertex0 = this->myBVHSet->Vertices[aTriangle.x()];
0239 const BVH_VecNt aVertex1 = this->myBVHSet->Vertices[aTriangle.y()];
0240 const BVH_VecNt aVertex2 = this->myBVHSet->Vertices[aTriangle.z()];
0241
0242 const BVH_VecNt aDirection =
0243 DirectionToNearestPoint<T, N> (this->myObject, aVertex0, aVertex1, aVertex2);
0244
0245 const T aDistance = BVH_DOT3 (aDirection, aDirection);
0246
0247 if (aDistance < this->myDistance)
0248 {
0249 this->myDistance = aDistance;
0250
0251 BVH_VecNt aTrgEdges[] = { aVertex1 - aVertex0, aVertex2 - aVertex0 };
0252
0253 BVH_VecNt aTrgNormal;
0254
0255 aTrgNormal.x () = aTrgEdges[0].y () * aTrgEdges[1].z () - aTrgEdges[0].z () * aTrgEdges[1].y ();
0256 aTrgNormal.y () = aTrgEdges[0].z () * aTrgEdges[1].x () - aTrgEdges[0].x () * aTrgEdges[1].z ();
0257 aTrgNormal.z () = aTrgEdges[0].x () * aTrgEdges[1].y () - aTrgEdges[0].y () * aTrgEdges[1].x ();
0258
0259 this->myIsOutside = BVH_DOT3 (aTrgNormal, aDirection) > 0;
0260
0261 return Standard_True;
0262 }
0263
0264 return Standard_False;
0265 }
0266 };
0267
0268 //=======================================================================
0269 //function : SquareDistanceToObject
0270 //purpose : Computes squared distance from point to BVH triangulation
0271 //=======================================================================
0272 template<class T, int N>
0273 T SquareDistanceToObject (BVH_Object<T, N>* theObject,
0274 const typename VectorType<T, N>::Type& thePnt, Standard_Boolean& theIsOutside)
0275 {
0276 Standard_STATIC_ASSERT (N == 3 || N == 4);
0277
0278 T aMinDistance = std::numeric_limits<T>::max();
0279
0280 BVH_Triangulation<T, N>* aTriangulation =
0281 dynamic_cast<BVH_Triangulation<T, N>*> (theObject);
0282
0283 Standard_ASSERT_RETURN (aTriangulation != NULL,
0284 "Error: Unsupported BVH object (non triangulation)", aMinDistance);
0285
0286 const opencascade::handle<BVH_Tree<T, N> >& aBVH = aTriangulation->BVH();
0287 if (aBVH.IsNull())
0288 {
0289 return Standard_False;
0290 }
0291
0292 PointTriangulationSquareDistance<T,N> aDistTool;
0293 aDistTool.SetObject (thePnt);
0294 aDistTool.SetBVHSet (aTriangulation);
0295 aDistTool.ComputeDistance();
0296 theIsOutside = aDistTool.IsOutside();
0297 return aDistTool.Distance();
0298 }
0299
0300 //=======================================================================
0301 //function : PointGeometrySquareDistance
0302 //purpose : Computes squared distance from point to BVH geometry
0303 //=======================================================================
0304 template<class T, int N>
0305 class PointGeometrySquareDistance: public SquareDistanceToPoint <T,N,BVH_Geometry<T, N>>
0306 {
0307 public:
0308 typedef typename VectorType<T, N>::Type BVH_VecNt;
0309
0310 public:
0311 //! Constructor
0312 PointGeometrySquareDistance()
0313 : SquareDistanceToPoint<T, N, BVH_Geometry<T, N>>()
0314 {
0315 }
0316
0317 public:
0318 // Accepting the element
0319 virtual Standard_Boolean Accept (const Standard_Integer theIndex,
0320 const T&) Standard_OVERRIDE
0321 {
0322 Standard_Boolean isOutside = Standard_True;
0323 const T aDistance = SquareDistanceToObject (
0324 this->myBVHSet->Objects ()(theIndex).operator->(), this->myObject, isOutside);
0325
0326 if (aDistance < this->myDistance)
0327 {
0328 this->myDistance = aDistance;
0329 this->myIsOutside = isOutside;
0330
0331 return Standard_True;
0332 }
0333 return Standard_False;
0334 }
0335 };
0336
0337 //=======================================================================
0338 //function : SquareDistanceToGeomerty
0339 //purpose : Computes squared distance from point to BVH geometry
0340 //=======================================================================
0341 template<class T, int N>
0342 T SquareDistanceToGeomerty (BVH_Geometry<T, N>& theGeometry,
0343 const typename VectorType<T, N>::Type& thePnt, Standard_Boolean& theIsOutside)
0344 {
0345 Standard_STATIC_ASSERT (N == 3 || N == 4);
0346
0347 const BVH_Tree<T, N, BVH_BinaryTree>* aBVH = theGeometry.BVH().get();
0348
0349 if (aBVH == NULL)
0350 {
0351 return Standard_False;
0352 }
0353
0354 PointGeometrySquareDistance<T,N> aDistTool;
0355 aDistTool.SetObject (thePnt);
0356 aDistTool.SetBVHSet (&theGeometry);
0357 aDistTool.ComputeDistance();
0358 theIsOutside = aDistTool.IsOutside();
0359 return aDistTool.Distance();
0360 }
0361 }
0362
0363 #undef BVH_DOT3
0364
0365 //! Tool object for parallel construction of distance field (uses Intel TBB).
0366 template<class T, int N>
0367 class BVH_ParallelDistanceFieldBuilder
0368 {
0369 private:
0370
0371 //! Input BVH geometry.
0372 BVH_Geometry<T, N>* myGeometry;
0373
0374 //! Output distance field.
0375 BVH_DistanceField<T, N>* myOutField;
0376
0377 public:
0378
0379 BVH_ParallelDistanceFieldBuilder (BVH_DistanceField<T, N>* theOutField, BVH_Geometry<T, N>* theGeometry)
0380 : myGeometry (theGeometry),
0381 myOutField (theOutField)
0382 {
0383 //
0384 }
0385
0386 void operator() (const Standard_Integer theIndex) const
0387 {
0388 myOutField->BuildSlices (*myGeometry, theIndex, theIndex + 1);
0389 }
0390 };
0391
0392 // =======================================================================
0393 // function : BuildSlices
0394 // purpose : Performs building of distance field for the given Z slices
0395 // =======================================================================
0396 template<class T, int N>
0397 void BVH_DistanceField<T, N>::BuildSlices (BVH_Geometry<T, N>& theGeometry,
0398 const Standard_Integer theStartSlice, const Standard_Integer theFinalSlice)
0399 {
0400 for (Standard_Integer aZ = theStartSlice; aZ < theFinalSlice; ++aZ)
0401 {
0402 for (Standard_Integer aY = 0; aY < myDimensionY; ++aY)
0403 {
0404 for (Standard_Integer aX = 0; aX < myDimensionX; ++aX)
0405 {
0406 BVH_VecNt aCenter;
0407
0408 aCenter.x() = myCornerMin.x() + myVoxelSize.x() * (aX + static_cast<T> (0.5));
0409 aCenter.y() = myCornerMin.y() + myVoxelSize.y() * (aY + static_cast<T> (0.5));
0410 aCenter.z() = myCornerMin.z() + myVoxelSize.z() * (aZ + static_cast<T> (0.5));
0411
0412 Standard_Boolean isOutside = Standard_True;
0413
0414 const T aDistance = sqrt (
0415 BVH::SquareDistanceToGeomerty<T, N> (theGeometry, aCenter, isOutside));
0416
0417 Voxel (aX, aY, aZ) = (!myComputeSign || isOutside) ? aDistance : -aDistance;
0418 }
0419 }
0420 }
0421 }
0422
0423 // =======================================================================
0424 // function : Build
0425 // purpose : Builds 3D distance field from BVH geometry
0426 // =======================================================================
0427 template<class T, int N>
0428 Standard_Boolean BVH_DistanceField<T, N>::Build (BVH_Geometry<T, N>& theGeometry)
0429 {
0430 if (theGeometry.Size() == 0)
0431 {
0432 return Standard_False;
0433 }
0434
0435 const BVH_VecNt aGlobalBoxSize = theGeometry.Box().Size();
0436
0437 const T aMaxBoxSide = Max (Max (aGlobalBoxSize.x(), aGlobalBoxSize.y()), aGlobalBoxSize.z());
0438
0439 myDimensionX = static_cast<Standard_Integer> (myMaximumSize * aGlobalBoxSize.x() / aMaxBoxSide);
0440 myDimensionY = static_cast<Standard_Integer> (myMaximumSize * aGlobalBoxSize.y() / aMaxBoxSide);
0441 myDimensionZ = static_cast<Standard_Integer> (myMaximumSize * aGlobalBoxSize.z() / aMaxBoxSide);
0442
0443 myDimensionX = Min (myMaximumSize, Max (myDimensionX, 16));
0444 myDimensionY = Min (myMaximumSize, Max (myDimensionY, 16));
0445 myDimensionZ = Min (myMaximumSize, Max (myDimensionZ, 16));
0446
0447 const BVH_VecNt aGlobalBoxMin = theGeometry.Box().CornerMin();
0448 const BVH_VecNt aGlobalBoxMax = theGeometry.Box().CornerMax();
0449
0450 const Standard_Integer aVoxelOffset = 2;
0451
0452 myCornerMin.x() = aGlobalBoxMin.x() - aVoxelOffset * aGlobalBoxSize.x() / (myDimensionX - 2 * aVoxelOffset);
0453 myCornerMin.y() = aGlobalBoxMin.y() - aVoxelOffset * aGlobalBoxSize.y() / (myDimensionY - 2 * aVoxelOffset);
0454 myCornerMin.z() = aGlobalBoxMin.z() - aVoxelOffset * aGlobalBoxSize.z() / (myDimensionZ - 2 * aVoxelOffset);
0455
0456 myCornerMax.x() = aGlobalBoxMax.x() + aVoxelOffset * aGlobalBoxSize.x() / (myDimensionX - 2 * aVoxelOffset);
0457 myCornerMax.y() = aGlobalBoxMax.y() + aVoxelOffset * aGlobalBoxSize.y() / (myDimensionY - 2 * aVoxelOffset);
0458 myCornerMax.z() = aGlobalBoxMax.z() + aVoxelOffset * aGlobalBoxSize.z() / (myDimensionZ - 2 * aVoxelOffset);
0459
0460 myVoxelSize.x() = (myCornerMax.x() - myCornerMin.x()) / myDimensionX;
0461 myVoxelSize.y() = (myCornerMax.y() - myCornerMin.y()) / myDimensionY;
0462 myVoxelSize.z() = (myCornerMax.z() - myCornerMin.z()) / myDimensionZ;
0463
0464 OSD_Parallel::For (0, myDimensionZ, BVH_ParallelDistanceFieldBuilder<T, N> (this, &theGeometry), !IsParallel());
0465
0466 return Standard_True;
0467 }
