1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195 |
- #pragma once
- #include <Mathematics/ConvexHull3.h>
- #include <Mathematics/MinimumAreaBox2.h>
- #include <thread>
- namespace WwiseGTE
- {
- template <typename InputType, typename ComputeType>
- class MinimumVolumeBox3
- {
- public:
-
-
-
-
-
-
-
- MinimumVolumeBox3(unsigned int numThreads = 1, bool threadProcessEdges = false)
- :
- mNumThreads(numThreads),
- mThreadProcessEdges(threadProcessEdges),
- mNumPoints(0),
- mPoints(nullptr),
- mComputePoints(nullptr),
- mUseRotatingCalipers(true),
- mVolume((InputType)0),
- mZero(0),
- mOne(1),
- mNegOne(-1),
- mHalf((InputType)0.5)
- {
- }
-
-
-
- OrientedBox3<InputType> operator()(int numPoints, Vector3<InputType> const* points,
- InputType epsilon = (InputType)0,
- bool useRotatingCalipers = !std::is_floating_point<ComputeType>::value)
- {
- mNumPoints = numPoints;
- mPoints = points;
- mUseRotatingCalipers = useRotatingCalipers;
- mHull.clear();
- mUniqueIndices.clear();
-
- ConvexHull3<InputType, ComputeType> ch3;
- ch3(mNumPoints, mPoints, epsilon);
- int dimension = ch3.GetDimension();
- OrientedBox3<InputType> itMinBox;
- if (dimension == 0)
- {
-
-
- itMinBox.center = mPoints[0];
- itMinBox.axis[0] = Vector3<InputType>::Unit(0);
- itMinBox.axis[1] = Vector3<InputType>::Unit(1);
- itMinBox.axis[2] = Vector3<InputType>::Unit(2);
- itMinBox.extent[0] = (InputType)0;
- itMinBox.extent[1] = (InputType)0;
- itMinBox.extent[2] = (InputType)0;
- mHull.resize(1);
- mHull[0] = 0;
- return itMinBox;
- }
- if (dimension == 1)
- {
-
-
-
-
- Line3<InputType> const& line = ch3.GetLine();
- InputType tmin = (InputType)0, tmax = (InputType)0;
- int imin = 0, imax = 0;
- for (int i = 0; i < mNumPoints; ++i)
- {
- Vector3<InputType> diff = mPoints[i] - line.origin;
- InputType t = Dot(diff, line.direction);
- if (t > tmax)
- {
- tmax = t;
- imax = i;
- }
- else if (t < tmin)
- {
- tmin = t;
- imin = i;
- }
- }
- itMinBox.center = line.origin + ((InputType)0.5) * (tmin + tmax) * line.direction;
- itMinBox.extent[0] = ((InputType)0.5) * (tmax - tmin);
- itMinBox.extent[1] = (InputType)0;
- itMinBox.extent[2] = (InputType)0;
- itMinBox.axis[0] = line.direction;
- ComputeOrthogonalComplement(1, &itMinBox.axis[0]);
- mHull.resize(2);
- mHull[0] = imin;
- mHull[1] = imax;
- return itMinBox;
- }
- if (dimension == 2)
- {
-
-
-
- Plane3<InputType> const& plane = ch3.GetPlane();
-
-
- Vector3<InputType> origin = mPoints[0];
- Vector3<InputType> basis[3];
- basis[0] = plane.normal;
- ComputeOrthogonalComplement(1, basis);
-
- std::vector<Vector2<InputType>> projection(mNumPoints);
- for (int i = 0; i < mNumPoints; ++i)
- {
- Vector3<InputType> diff = mPoints[i] - origin;
- projection[i][0] = Dot(basis[1], diff);
- projection[i][1] = Dot(basis[2], diff);
- }
-
- MinimumAreaBox2<InputType, ComputeType> mab2;
- OrientedBox2<InputType> rectangle = mab2(mNumPoints, &projection[0]);
-
- itMinBox.center = origin + rectangle.center[0] * basis[1] + rectangle.center[1] * basis[2];
- itMinBox.axis[0] = rectangle.axis[0][0] * basis[1] + rectangle.axis[0][1] * basis[2];
- itMinBox.axis[1] = rectangle.axis[1][0] * basis[1] + rectangle.axis[1][1] * basis[2];
- itMinBox.axis[2] = basis[0];
- itMinBox.extent[0] = rectangle.extent[0];
- itMinBox.extent[1] = rectangle.extent[1];
- itMinBox.extent[2] = (InputType)0;
- mHull = mab2.GetHull();
- return itMinBox;
- }
-
-
- ETManifoldMesh const& mesh = ch3.GetHullMesh();
- mHull.resize(3 * mesh.GetTriangles().size());
- int h = 0;
- for (auto const& element : mesh.GetTriangles())
- {
- for (int i = 0; i < 3; ++i, ++h)
- {
- int index = element.first.V[i];
- mHull[h] = index;
- mUniqueIndices.insert(index);
- }
- }
- mComputePoints = ch3.GetQuery().GetVertices();
- Box minBox, minBoxEdges;
- minBox.volume = mNegOne;
- minBoxEdges.volume = mNegOne;
- if (mThreadProcessEdges)
- {
- std::thread doEdges(
- [this, &mesh, &minBoxEdges]()
- {
- ProcessEdges(mesh, minBoxEdges);
- });
- ProcessFaces(mesh, minBox);
- doEdges.join();
- }
- else
- {
- ProcessEdges(mesh, minBoxEdges);
- ProcessFaces(mesh, minBox);
- }
- if (minBoxEdges.volume != mNegOne
- && minBoxEdges.volume < minBox.volume)
- {
- minBox = minBoxEdges;
- }
- ConvertTo(minBox, itMinBox);
- mComputePoints = nullptr;
- return itMinBox;
- }
-
-
- OrientedBox3<InputType> operator()(int numPoints, Vector3<InputType> const* points,
- int numIndices, int const* indices,
- bool useRotatingCalipers = !std::is_floating_point<ComputeType>::value)
- {
- mNumPoints = numPoints;
- mPoints = points;
- mUseRotatingCalipers = useRotatingCalipers;
- mUniqueIndices.clear();
-
-
- ETManifoldMesh mesh;
- int numTriangles = numIndices / 3;
- for (int t = 0; t < numTriangles; ++t)
- {
- int v0 = *indices++;
- int v1 = *indices++;
- int v2 = *indices++;
- mesh.Insert(v0, v1, v2);
- }
-
-
- mHull.resize(3 * mesh.GetTriangles().size());
- int h = 0;
- for (auto const& element : mesh.GetTriangles())
- {
- for (int i = 0; i < 3; ++i, ++h)
- {
- int index = element.first.V[i];
- mHull[h] = index;
- mUniqueIndices.insert(index);
- }
- }
-
- std::vector<Vector3<ComputeType>> computePoints(mNumPoints);
- for (auto i : mUniqueIndices)
- {
- for (int j = 0; j < 3; ++j)
- {
- computePoints[i][j] = (ComputeType)mPoints[i][j];
- }
- }
- OrientedBox3<InputType> itMinBox;
- mComputePoints = &computePoints[0];
- Box minBox, minBoxEdges;
- minBox.volume = mNegOne;
- minBoxEdges.volume = mNegOne;
- if (mThreadProcessEdges)
- {
- std::thread doEdges(
- [this, &mesh, &minBoxEdges]()
- {
- ProcessEdges(mesh, minBoxEdges);
- });
- ProcessFaces(mesh, minBox);
- doEdges.join();
- }
- else
- {
- ProcessEdges(mesh, minBoxEdges);
- ProcessFaces(mesh, minBox);
- }
- if (minBoxEdges.volume != mNegOne && minBoxEdges.volume < minBox.volume)
- {
- minBox = minBoxEdges;
- }
- ConvertTo(minBox, itMinBox);
- mComputePoints = nullptr;
- return itMinBox;
- }
-
- inline int GetNumPoints() const
- {
- return mNumPoints;
- }
- inline Vector3<InputType> const* GetPoints() const
- {
- return mPoints;
- }
- inline std::vector<int> const& GetHull() const
- {
- return mHull;
- }
- inline InputType GetVolume() const
- {
- return mVolume;
- }
- private:
- struct Box
- {
- Vector3<ComputeType> P, U[3];
- ComputeType sqrLenU[3], range[3][2], volume;
- };
- struct ExtrudeRectangle
- {
- Vector3<ComputeType> U[2];
- std::array<int, 4> index;
- ComputeType sqrLenU[2], area;
- };
-
- void ProcessFaces(ETManifoldMesh const& mesh, Box& minBox)
- {
-
- auto const& tmap = mesh.GetTriangles();
- auto const& emap = mesh.GetEdges();
-
-
-
-
-
- std::vector<Vector3<ComputeType>> normal(tmap.size());
- std::map<std::shared_ptr<Triangle>, int> triNormalMap;
- int index = 0;
- for (auto const& element : tmap)
- {
- auto tri = element.second;
- Vector3<ComputeType> const& v0 = mComputePoints[tri->V[0]];
- Vector3<ComputeType> const& v1 = mComputePoints[tri->V[1]];
- Vector3<ComputeType> const& v2 = mComputePoints[tri->V[2]];
- Vector3<ComputeType> edge1 = v1 - v0;
- Vector3<ComputeType> edge2 = v2 - v0;
- normal[index] = Cross(edge2, edge1);
- if (normal[index] == Vector3<ComputeType>(0))
- continue;
- triNormalMap[tri] = index++;
- }
-
-
-
-
-
-
- unsigned int numFaces = static_cast<unsigned int>(tmap.size());
- if (mNumThreads > 1 && numFaces >= mNumThreads)
- {
-
-
- std::vector<std::shared_ptr<Triangle>> triangles;
- triangles.reserve(numFaces);
- for (auto const& element : tmap)
- {
- triangles.push_back(element.second);
- }
-
- unsigned int numFacesPerThread = numFaces / mNumThreads;
- std::vector<unsigned int> imin(mNumThreads), imax(mNumThreads);
- std::vector<Box> localMinBox(mNumThreads);
- for (unsigned int t = 0; t < mNumThreads; ++t)
- {
- imin[t] = t * numFacesPerThread;
- imax[t] = imin[t] + numFacesPerThread - 1;
- localMinBox[t].volume = mNegOne;
- }
- imax[mNumThreads - 1] = numFaces - 1;
-
- std::vector<std::thread> process(mNumThreads);
- for (unsigned int t = 0; t < mNumThreads; ++t)
- {
- process[t] = std::thread([this, t, &imin, &imax, &triangles,
- &normal, &triNormalMap, &emap, &localMinBox]()
- {
- for (unsigned int i = imin[t]; i <= imax[t]; ++i)
- {
- auto const& supportTri = triangles[i];
- ProcessFace(supportTri, normal, triNormalMap, emap, localMinBox[t]);
- }
- });
- }
-
- for (unsigned int t = 0; t < mNumThreads; ++t)
- {
- process[t].join();
-
- if (minBox.volume == mNegOne || localMinBox[t].volume < minBox.volume)
- {
- minBox = localMinBox[t];
- }
- }
- }
- else
- {
- for (auto const& element : tmap)
- {
- auto const& supportTri = element.second;
- ProcessFace(supportTri, normal, triNormalMap, emap, minBox);
- }
- }
- }
-
-
- void ProcessEdges(ETManifoldMesh const& mesh, Box& minBox)
- {
-
-
-
-
-
- int index = mesh.GetTriangles().begin()->first.V[0];
- Vector3<ComputeType> const& origin = mComputePoints[index];
- Vector3<ComputeType> U[3];
- std::array<ComputeType, 3> sqrLenU;
- auto const& emap = mesh.GetEdges();
- auto e2 = emap.begin(), end = emap.end();
- for (; e2 != end; ++e2)
- {
- U[2] = mComputePoints[e2->first.V[1]] - mComputePoints[e2->first.V[0]];
- auto e1 = e2;
- for (++e1; e1 != end; ++e1)
- {
- U[1] = mComputePoints[e1->first.V[1]] - mComputePoints[e1->first.V[0]];
- if (Dot(U[1], U[2]) != mZero)
- {
- continue;
- }
- sqrLenU[1] = Dot(U[1], U[1]);
- auto e0 = e1;
- for (++e0; e0 != end; ++e0)
- {
- U[0] = mComputePoints[e0->first.V[1]] - mComputePoints[e0->first.V[0]];
- sqrLenU[0] = Dot(U[0], U[0]);
- if (Dot(U[0], U[1]) != mZero || Dot(U[0], U[2]) != mZero)
- {
- continue;
- }
-
-
-
- U[2] = Cross(U[0], U[1]);
- sqrLenU[2] = sqrLenU[0] * sqrLenU[1];
-
-
- std::array<ComputeType, 3> umin, umax;
- for (int j = 0; j < 3; ++j)
- {
- umin[j] = mZero;
- umax[j] = mZero;
- }
- for (auto i : mUniqueIndices)
- {
- Vector3<ComputeType> diff = mComputePoints[i] - origin;
- for (int j = 0; j < 3; ++j)
- {
- ComputeType dot = Dot(diff, U[j]);
- if (dot < umin[j])
- {
- umin[j] = dot;
- }
- else if (dot > umax[j])
- {
- umax[j] = dot;
- }
- }
- }
- ComputeType volume = (umax[0] - umin[0]) / sqrLenU[0];
- volume *= (umax[1] - umin[1]) / sqrLenU[1];
- volume *= (umax[2] - umin[2]);
-
- if (minBox.volume == mOne || volume < minBox.volume)
- {
-
-
-
- if (DotCross(U[0], U[1], U[2]) < mZero)
- {
- U[2] = -U[2];
- }
- minBox.P = origin;
- for (int j = 0; j < 3; ++j)
- {
- minBox.U[j] = U[j];
- minBox.sqrLenU[j] = sqrLenU[j];
- for (int k = 0; k < 3; ++k)
- {
- minBox.range[k][0] = umin[k];
- minBox.range[k][1] = umax[k];
- }
- }
- minBox.volume = volume;
- }
- }
- }
- }
- }
-
- typedef ETManifoldMesh::Triangle Triangle;
- void ProcessFace(std::shared_ptr<Triangle> const& supportTri,
- std::vector<Vector3<ComputeType>> const& normal,
- std::map<std::shared_ptr<Triangle>, int> const& triNormalMap,
- ETManifoldMesh::EMap const& emap, Box& localMinBox)
- {
-
- Vector3<ComputeType> const& supportNormal = normal[triNormalMap.find(supportTri)->second];
- static const Vector3<ComputeType> sZero = { 0.0, 0.0, 0.0 };
- if (supportNormal == sZero)
- return;
-
-
-
- std::vector<int> polyline(mNumPoints);
- int polylineStart = -1;
- for (auto const& edgeElement : emap)
- {
- auto const& edge = *edgeElement.second;
- auto const& tri0 = edge.T[0].lock();
- auto const& tri1 = edge.T[1].lock();
- auto const& normal0 = normal[triNormalMap.find(tri0)->second];
- auto const& normal1 = normal[triNormalMap.find(tri1)->second];
- ComputeType dot0 = Dot(supportNormal, normal0);
- ComputeType dot1 = Dot(supportNormal, normal1);
- std::shared_ptr<Triangle> tri;
- if (dot0 < mZero && dot1 >= mZero)
- {
- tri = tri0;
- }
- else if (dot1 < mZero && dot0 >= mZero)
- {
- tri = tri1;
- }
- if (tri)
- {
-
-
-
-
-
- for (int j0 = 2, j1 = 0; j1 < 3; j0 = j1++)
- {
- if (tri->V[j1] == edge.V[0])
- {
- if (tri->V[j0] == edge.V[1])
- {
- polyline[edge.V[1]] = edge.V[0];
- }
- else
- {
- polyline[edge.V[0]] = edge.V[1];
- }
- polylineStart = edge.V[0];
- break;
- }
- }
- }
- }
- if (polylineStart == -1)
- return;
-
-
-
- std::vector<int> closedPolyline(mNumPoints);
- int numClosedPolyline = 0;
- int v = polylineStart;
- for (auto& cp : closedPolyline)
- {
- cp = v;
- ++numClosedPolyline;
- v = polyline[v];
- if (v == polylineStart)
- {
- break;
- }
- }
- closedPolyline.resize(numClosedPolyline);
-
-
- RemoveCollinearPoints(supportNormal, closedPolyline);
-
-
- Box faceBox;
- if (mUseRotatingCalipers)
- {
- ComputeBoxForFaceOrderN(supportNormal, closedPolyline, faceBox);
- }
- else
- {
- ComputeBoxForFaceOrderNSqr(supportNormal, closedPolyline, faceBox);
- }
-
- if (localMinBox.volume == mNegOne || faceBox.volume < localMinBox.volume)
- {
- localMinBox = faceBox;
- }
- }
-
-
-
-
- void RemoveCollinearPoints(Vector3<ComputeType> const& N, std::vector<int>& polyline)
- {
- std::vector<int> tmpPolyline = polyline;
- int const numPolyline = static_cast<int>(polyline.size());
- int numNoncollinear = 0;
- Vector3<ComputeType> ePrev =
- mComputePoints[tmpPolyline[0]] - mComputePoints[tmpPolyline.back()];
- for (int i0 = 0, i1 = 1; i0 < numPolyline; ++i0)
- {
- Vector3<ComputeType> eNext =
- mComputePoints[tmpPolyline[i1]] - mComputePoints[tmpPolyline[i0]];
- ComputeType tsp = DotCross(ePrev, eNext, N);
- if (tsp != mZero)
- {
- polyline[numNoncollinear++] = tmpPolyline[i0];
- }
- ePrev = eNext;
- if (++i1 == numPolyline)
- {
- i1 = 0;
- }
- }
- polyline.resize(numNoncollinear);
- }
-
- void ComputeBoxForFaceOrderNSqr(Vector3<ComputeType> const& N, std::vector<int> const& polyline, Box& box)
- {
-
-
-
-
-
-
-
- box.P = mComputePoints[polyline[0]];
- box.U[2] = N;
- box.sqrLenU[2] = Dot(N, N);
- box.range[1][0] = mZero;
- box.volume = mNegOne;
- int const numPolyline = static_cast<int>(polyline.size());
- for (int i0 = numPolyline - 1, i1 = 0; i1 < numPolyline; i0 = i1++)
- {
-
-
- Vector3<ComputeType> const& P = mComputePoints[polyline[i0]];
- Vector3<ComputeType> E = mComputePoints[polyline[i1]] - mComputePoints[polyline[i0]];
- Vector3<ComputeType> U1 = Cross(N, E);
- Vector3<ComputeType> U0 = Cross(U1, N);
-
-
- ComputeType min0 = mZero, max0 = mZero, max1 = mZero;
- for (int j = 0; j < numPolyline; ++j)
- {
- Vector3<ComputeType> diff = mComputePoints[polyline[j]] - P;
- ComputeType dot = Dot(U0, diff);
- if (dot < min0)
- {
- min0 = dot;
- }
- else if (dot > max0)
- {
- max0 = dot;
- }
- dot = Dot(U1, diff);
- if (dot > max1)
- {
- max1 = dot;
- }
- }
-
-
-
-
- ComputeType sqrLenU1 = Dot(U1, U1);
- ComputeType volume = (max0 - min0) * max1 / sqrLenU1;
- if (box.volume == mNegOne || volume < box.volume)
- {
- box.P = P;
- box.U[0] = U0;
- box.U[1] = U1;
- box.sqrLenU[0] = sqrLenU1 * box.sqrLenU[2];
- box.sqrLenU[1] = sqrLenU1;
- box.range[0][0] = min0;
- box.range[0][1] = max0;
- box.range[1][1] = max1;
- box.volume = volume;
- }
- }
-
- box.range[2][0] = mZero;
- box.range[2][1] = mZero;
- for (auto i : mUniqueIndices)
- {
- Vector3<ComputeType> diff = mComputePoints[i] - box.P;
- ComputeType height = Dot(box.U[2], diff);
- if (height < box.range[2][0])
- {
- box.range[2][0] = height;
- }
- else if (height > box.range[2][1])
- {
- box.range[2][1] = height;
- }
- }
-
- box.volume *= (box.range[2][1] - box.range[2][0]) / box.sqrLenU[2];
- }
-
- void ComputeBoxForFaceOrderN(Vector3<ComputeType> const& N, std::vector<int> const& polyline, Box& box)
- {
-
-
-
-
-
-
-
-
-
-
- std::vector<bool> visited(polyline.size());
- std::fill(visited.begin(), visited.end(), false);
-
-
-
-
-
-
-
-
- ExtrudeRectangle minRct =
- SmallestRectangle((int)polyline.size() - 1, 0, N, polyline);
- visited[minRct.index[0]] = true;
-
- ExtrudeRectangle rct = minRct;
- for (size_t i = 0; i < polyline.size(); ++i)
- {
- std::array<std::pair<ComputeType, int>, 4> A;
- int numA;
- if (!ComputeAngles(N, polyline, rct, A, numA))
- {
-
-
- break;
- }
-
- std::array<int, 4> sort = SortAngles(A, numA);
-
-
- if (!UpdateSupport(A, numA, sort, N, polyline, visited, rct))
- {
-
-
- break;
- }
- if (rct.area < minRct.area)
- {
- minRct = rct;
- }
- }
-
-
- box.P = mComputePoints[polyline[minRct.index[0]]];
- box.U[0] = minRct.U[0];
- box.U[1] = minRct.U[1];
- box.U[2] = N;
- box.sqrLenU[0] = minRct.sqrLenU[0];
- box.sqrLenU[1] = minRct.sqrLenU[1];
- box.sqrLenU[2] = Dot(box.U[2], box.U[2]);
-
-
- box.range[0][0] = Dot(box.U[0], mComputePoints[polyline[minRct.index[3]]] - box.P);
- box.range[0][1] = Dot(box.U[0], mComputePoints[polyline[minRct.index[1]]] - box.P);
- box.range[1][0] = mZero;
- box.range[1][1] = Dot(box.U[1], mComputePoints[polyline[minRct.index[2]]] - box.P);
-
- box.range[2][0] = mZero;
- box.range[2][1] = mZero;
- for (auto i : mUniqueIndices)
- {
- Vector3<ComputeType> diff = mComputePoints[i] - box.P;
- ComputeType height = Dot(box.U[2], diff);
- if (height < box.range[2][0])
- {
- box.range[2][0] = height;
- }
- else if (height > box.range[2][1])
- {
- box.range[2][1] = height;
- }
- }
-
- box.volume =
- (box.range[0][1] - box.range[0][0]) *
- ((box.range[1][1] - box.range[1][0]) / box.sqrLenU[1]) *
- ((box.range[2][1] - box.range[2][0]) / box.sqrLenU[2]);
- }
-
- ExtrudeRectangle SmallestRectangle(int i0, int i1, Vector3<ComputeType> const& N, std::vector<int> const& polyline)
- {
- ExtrudeRectangle rct;
- Vector3<ComputeType> E = mComputePoints[polyline[i1]] - mComputePoints[polyline[i0]];
- rct.U[1] = Cross(N, E);
- rct.U[0] = Cross(rct.U[1], N);
- rct.index = { i1, i1, i1, i1 };
- rct.sqrLenU[0] = Dot(rct.U[0], rct.U[0]);
- rct.sqrLenU[1] = Dot(rct.U[1], rct.U[1]);
- Vector3<ComputeType> const& origin = mComputePoints[polyline[i1]];
- Vector2<ComputeType> support[4];
- for (int j = 0; j < 4; ++j)
- {
- support[j] = { mZero, mZero };
- }
- int i = 0;
- for (auto p : polyline)
- {
- Vector3<ComputeType> diff = mComputePoints[p] - origin;
- Vector2<ComputeType> v = { Dot(rct.U[0], diff), Dot(rct.U[1], diff) };
-
-
-
-
-
-
- if (v[0] > support[1][0] ||
- (v[0] == support[1][0] && v[1] > support[1][1]))
- {
-
-
- rct.index[1] = i;
- support[1] = v;
- }
- if (v[1] > support[2][1] ||
- (v[1] == support[2][1] && v[0] < support[2][0]))
- {
-
-
- rct.index[2] = i;
- support[2] = v;
- }
- if (v[0] < support[3][0] ||
- (v[0] == support[3][0] && v[1] < support[3][1]))
- {
-
-
- rct.index[3] = i;
- support[3] = v;
- }
- ++i;
- }
-
-
-
- ComputeType scaledWidth = support[1][0] - support[3][0];
- ComputeType scaledHeight = support[2][1];
- rct.area = scaledWidth * scaledHeight / rct.sqrLenU[1];
- return rct;
- }
-
-
-
-
- bool ComputeAngles(Vector3<ComputeType> const& N,
- std::vector<int> const& polyline, ExtrudeRectangle const& rct,
- std::array<std::pair<ComputeType, int>, 4>& A, int& numA) const
- {
- int const numPolyline = static_cast<int>(polyline.size());
- numA = 0;
- for (int k0 = 3, k1 = 0; k1 < 4; k0 = k1++)
- {
- if (rct.index[k0] != rct.index[k1])
- {
-
-
- int lookup = (k0 & 1);
- Vector3<ComputeType> D = ((k0 & 2) ? -rct.U[lookup] : rct.U[lookup]);
- int j0 = rct.index[k0], j1 = j0 + 1;
- if (j1 == numPolyline)
- {
- j1 = 0;
- }
- Vector3<ComputeType> E = mComputePoints[polyline[j1]] - mComputePoints[polyline[j0]];
- Vector3<ComputeType> Eperp = Cross(N, E);
- E = Cross(Eperp, N);
- Vector3<ComputeType> DxE = Cross(D, E);
- ComputeType csqrlen = Dot(DxE, DxE);
- ComputeType dsqrlen = rct.sqrLenU[lookup];
- ComputeType esqrlen = Dot(E, E);
- ComputeType sinThetaSqr = csqrlen / (dsqrlen * esqrlen);
- A[numA++] = std::make_pair(sinThetaSqr, k0);
- }
- }
- return numA > 0;
- }
-
-
-
- std::array<int, 4> SortAngles(std::array<std::pair<ComputeType, int>, 4> const& A, int numA) const
- {
- std::array<int, 4> sort = { 0, 1, 2, 3 };
- if (numA > 1)
- {
- if (numA == 2)
- {
- if (A[sort[0]].first > A[sort[1]].first)
- {
- std::swap(sort[0], sort[1]);
- }
- }
- else if (numA == 3)
- {
- if (A[sort[0]].first > A[sort[1]].first)
- {
- std::swap(sort[0], sort[1]);
- }
- if (A[sort[0]].first > A[sort[2]].first)
- {
- std::swap(sort[0], sort[2]);
- }
- if (A[sort[1]].first > A[sort[2]].first)
- {
- std::swap(sort[1], sort[2]);
- }
- }
- else
- {
- if (A[sort[0]].first > A[sort[1]].first)
- {
- std::swap(sort[0], sort[1]);
- }
- if (A[sort[2]].first > A[sort[3]].first)
- {
- std::swap(sort[2], sort[3]);
- }
- if (A[sort[0]].first > A[sort[2]].first)
- {
- std::swap(sort[0], sort[2]);
- }
- if (A[sort[1]].first > A[sort[3]].first)
- {
- std::swap(sort[1], sort[3]);
- }
- if (A[sort[1]].first > A[sort[2]].first)
- {
- std::swap(sort[1], sort[2]);
- }
- }
- }
- return sort;
- }
- bool UpdateSupport(std::array<std::pair<ComputeType, int>, 4> const& A, int numA,
- std::array<int, 4> const& sort, Vector3<ComputeType> const& N, std::vector<int> const& polyline,
- std::vector<bool>& visited, ExtrudeRectangle& rct)
- {
-
-
- int const numPolyline = static_cast<int>(polyline.size());
- auto const& amin = A[sort[0]];
- for (int k = 0; k < numA; ++k)
- {
- auto const& a = A[sort[k]];
- if (a.first == amin.first)
- {
- if (++rct.index[a.second] == numPolyline)
- {
- rct.index[a.second] = 0;
- }
- }
- else
- {
- break;
- }
- }
- int bottom = rct.index[amin.second];
- if (visited[bottom])
- {
-
- return false;
- }
- visited[bottom] = true;
-
- std::array<int, 4> nextIndex;
- for (int k = 0; k < 4; ++k)
- {
- nextIndex[k] = rct.index[(amin.second + k) % 4];
- }
- rct.index = nextIndex;
-
- int j1 = rct.index[0], j0 = j1 - 1;
- if (j1 < 0)
- {
- j1 = numPolyline - 1;
- }
- Vector3<ComputeType> E =
- mComputePoints[polyline[j1]] - mComputePoints[polyline[j0]];
- rct.U[1] = Cross(N, E);
- rct.U[0] = Cross(rct.U[1], N);
- rct.sqrLenU[0] = Dot(rct.U[0], rct.U[0]);
- rct.sqrLenU[1] = Dot(rct.U[1], rct.U[1]);
-
- Vector3<ComputeType> diff[2] =
- {
- mComputePoints[polyline[rct.index[1]]]
- - mComputePoints[polyline[rct.index[3]]],
- mComputePoints[polyline[rct.index[2]]]
- - mComputePoints[polyline[rct.index[0]]]
- };
- ComputeType scaledWidth = Dot(rct.U[0], diff[0]);
- ComputeType scaledHeight = Dot(rct.U[1], diff[1]);
- rct.area = scaledWidth * scaledHeight / rct.sqrLenU[1];
- return true;
- }
-
-
-
- void ConvertTo(Box const& minBox, OrientedBox3<InputType>& itMinBox)
- {
- Vector3<ComputeType> center = minBox.P;
- for (int i = 0; i < 3; ++i)
- {
- ComputeType average = mHalf * (minBox.range[i][0] + minBox.range[i][1]);
- center += (average / minBox.sqrLenU[i]) * minBox.U[i];
- }
-
-
- Vector3<ComputeType> sqrExtent;
- for (int i = 0; i < 3; ++i)
- {
- sqrExtent[i] = mHalf * (minBox.range[i][1] - minBox.range[i][0]);
- sqrExtent[i] *= sqrExtent[i];
- sqrExtent[i] /= minBox.sqrLenU[i];
- }
- for (int i = 0; i < 3; ++i)
- {
- itMinBox.center[i] = (InputType)center[i];
- itMinBox.extent[i] = std::sqrt((InputType)sqrExtent[i]);
-
-
-
-
- Vector3<ComputeType> const& axis = minBox.U[i];
- ComputeType cmax = std::max<ComputeType>(std::abs((double)axis[0]), std::abs((double)axis[1]));
- cmax = std::max<ComputeType>(cmax, std::abs((double)axis[2]));
- ComputeType invCMax = mOne / cmax;
- for (int j = 0; j < 3; ++j)
- {
- itMinBox.axis[i][j] = (InputType)(axis[j] * invCMax);
- }
- Normalize(itMinBox.axis[i]);
- }
- mVolume = (InputType)minBox.volume;
- }
-
-
-
-
- unsigned int mNumThreads;
- bool mThreadProcessEdges;
-
- int mNumPoints;
- Vector3<InputType> const* mPoints;
-
-
-
- Vector3<ComputeType> const* mComputePoints;
-
-
- std::vector<int> mHull;
-
-
-
-
- std::set<int> mUniqueIndices;
-
-
- bool mUseRotatingCalipers;
-
-
-
-
- InputType mVolume;
-
-
- ComputeType mZero, mOne, mNegOne, mHalf;
- };
- }
|