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