diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.cpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.cpp index 505097db3a..2f870fe205 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.cpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.cpp @@ -7,14 +7,16 @@ #include "simplnx/DataStructure/Geometry/EdgeGeom.hpp" #include "simplnx/DataStructure/Geometry/ImageGeom.hpp" #include "simplnx/DataStructure/Geometry/RectGridGeom.hpp" +#include "simplnx/Utilities/ParallelTaskAlgorithm.hpp" using namespace nx::core; namespace { + // ---------------------------------------------------------------------------- // -inline std::array GetEdgeCoordinates(usize edgeId, INodeGeometry0D::SharedVertexList& verts, INodeGeometry1D::SharedEdgeList& edges) +inline std::array GetEdgeCoordinates(usize edgeId, const INodeGeometry0D::SharedVertexList& verts, const INodeGeometry1D::SharedEdgeList& edges) { usize v0Idx = edges[edgeId * 2]; usize v1Idx = edges[edgeId * 2 + 1]; @@ -23,7 +25,8 @@ inline std::array GetEdgeCoordinates(usize edgeId, INodeG // ---------------------------------------------------------------------------- // Helper function to check if a point lies inside a polygon using ray-casting -bool pointInPolygon(const EdgeGeom& edgeGeom, const std::vector& edgeIndices, const Point3Df& point, INodeGeometry0D::SharedVertexList& verts, INodeGeometry1D::SharedEdgeList& edges) +bool pointInPolygon(const EdgeGeom& edgeGeom, const std::vector& edgeIndices, const Point3Df& point, const INodeGeometry0D::SharedVertexList& verts, + const INodeGeometry1D::SharedEdgeList& edges) { size_t intersections = 0; size_t numEdges = edgeIndices.size(); @@ -57,6 +60,97 @@ bool pointInPolygon(const EdgeGeom& edgeGeom, const std::vector& edgeIndi } return (intersections % 2) == 1; } + +// ---------------------------------------------------------------------------- +// +class SampleSurfaceMeshSliceImpl +{ +public: + SampleSurfaceMeshSliceImpl() = delete; + SampleSurfaceMeshSliceImpl(const SampleSurfaceMeshSliceImpl&) = default; + + SampleSurfaceMeshSliceImpl(RegularGridSampleSurfaceMesh* filterAlg, const EdgeGeom& edgeGeom, int32 currentSliceId, usize imageGeomIdx, const ImageGeom& imageGeom, const Int32Array& sliceIds, + Int32Array& featureIds, const std::atomic_bool& shouldCancel) + : m_FilterAlg(filterAlg) + , m_EdgeGeom(edgeGeom) + , m_CurrentSliceId(currentSliceId) + , m_ImageGeomIdx(imageGeomIdx) + , m_ImageGeom(imageGeom) + , m_SliceIds(sliceIds) + , m_FeatureIds(featureIds) + , m_ShouldCancel(shouldCancel) + { + } + SampleSurfaceMeshSliceImpl(SampleSurfaceMeshSliceImpl&&) = default; // Move Constructor Not Implemented + SampleSurfaceMeshSliceImpl& operator=(const SampleSurfaceMeshSliceImpl&) = delete; // Copy Assignment Not Implemented + SampleSurfaceMeshSliceImpl& operator=(SampleSurfaceMeshSliceImpl&&) = delete; // Move Assignment Not Implemented + + ~SampleSurfaceMeshSliceImpl() = default; + + void operator()() const + { + auto start = std::chrono::steady_clock::now(); + usize numEdges = m_EdgeGeom.getNumberOfEdges(); + std::vector edgeIndices; + edgeIndices.reserve(1024); // Reserve some space in the vector. This is just a guess. + SizeVec3 dimensions = m_ImageGeom.getDimensions(); + size_t cellsPerSlice = dimensions[0] * dimensions[1]; + const INodeGeometry0D::SharedVertexList& verts = m_EdgeGeom.getVerticesRef(); + const INodeGeometry1D::SharedEdgeList& edges = m_EdgeGeom.getEdgesRef(); + + // Loop over all edges and find the edges that are just for the current Slice Id + for(usize edgeIdx = 0; edgeIdx < numEdges; edgeIdx++) + { + int32 sliceIndex = m_SliceIds[edgeIdx]; + if(m_CurrentSliceId == sliceIndex) + { + edgeIndices.push_back(edgeIdx); + } + } + + std::vector featureIds(cellsPerSlice, 0); + + // Now that we have the edges that are on this slice, iterate over all + // voxels on this slice + for(size_t planeIdx = 0; planeIdx < cellsPerSlice; planeIdx++) + { + Point3Df imagePoint = m_ImageGeom.getCoordsf(m_ImageGeomIdx + planeIdx); + imagePoint[2] = 0.0f; // Force this down to the zero plane. + + if(pointInPolygon(m_EdgeGeom, edgeIndices, imagePoint, verts, edges)) + { + // featureIds[m_ImageGeomIdx + planeIdx] = 1; + featureIds[planeIdx] = 1; // Parallel version + } + } + + m_FilterAlg->sendThreadSafeUpdate(m_FeatureIds, featureIds, m_ImageGeomIdx); + + // for(usize i = 0; i < m_TotalPoints; i++) + // { + // auto now = std::chrono::steady_clock::now(); + // //// Only send updates every 1 second + // if(std::chrono::duration_cast(now - start).count() > 1000) + // { + // const float progress = static_cast(i) / static_cast(m_TotalPoints) * 100.0f; + // m_FilterAlg->updateProgress(fmt::format("Processing {}: {}% completed", m_DataArrayPtr->getName(), static_cast(progress))); + // start = std::chrono::steady_clock::now(); + // } + // } + } + +private: + RegularGridSampleSurfaceMesh* m_FilterAlg = nullptr; + ; + const EdgeGeom& m_EdgeGeom; + int32 m_CurrentSliceId; + usize m_ImageGeomIdx; + const ImageGeom m_ImageGeom; + const Int32Array& m_SliceIds; + Int32Array& m_FeatureIds; + const std::atomic_bool& m_ShouldCancel; +}; + } // namespace // ----------------------------------------------------------------------------- @@ -106,13 +200,6 @@ void RegularGridSampleSurfaceMesh::generatePoints(std::vector& points) // ----------------------------------------------------------------------------- Result<> RegularGridSampleSurfaceMesh::operator()() { - - // SampleSurfaceMeshInputValues inputs; - // inputs.TriangleGeometryPath = m_InputValues->TriangleGeometryPath; - // inputs.SurfaceMeshFaceLabelsArrayPath = m_InputValues->SurfaceMeshFaceLabelsArrayPath; - // inputs.FeatureIdsArrayPath = m_InputValues->FeatureIdsArrayPath; - // return execute(inputs); - // Slice the Triangle Geometry SliceTriangleGeometryInputValues inputValues; inputValues.SliceRange = 1; @@ -140,21 +227,25 @@ Result<> RegularGridSampleSurfaceMesh::operator()() auto& sliceId = m_DataStructure.getDataRefAs(sliceIdDataPath); INodeGeometry0D::SharedVertexList& verts = edgeGeom.getVerticesRef(); INodeGeometry1D::SharedEdgeList& edges = edgeGeom.getEdgesRef(); - usize numEdges = edgeGeom.getNumberOfEdges(); + // usize numEdges = edgeGeom.getNumberOfEdges(); // Get the Image Geometry that is the sampling Grid auto& imageGeom = m_DataStructure.getDataRefAs(m_InputValues->ImageGeometryOutputPath); FloatVec3 origin = imageGeom.getOrigin(); FloatVec3 spacing = imageGeom.getSpacing(); SizeVec3 dimensions = imageGeom.getDimensions(); - size_t cellsPerSlice = dimensions[0] * dimensions[1]; + // size_t cellsPerSlice = dimensions[0] * dimensions[1]; // Get the Feature Ids array - auto& featureIds = m_DataStructure.getDataAs(m_InputValues->FeatureIdsArrayPath)->getDataStoreRef(); + auto featureIds = m_DataStructure.getDataRefAs(m_InputValues->FeatureIdsArrayPath); //->getDataStoreRef(); std::vector edgeIndices; edgeIndices.reserve(1024); // Reserve some space in the vector. This is just a guess. + auto nthreads = static_cast(std::thread::hardware_concurrency()); // Returns ZERO if not defined on this platform + ParallelTaskAlgorithm taskRunner; + taskRunner.setParallelizationEnabled(true); + int32 currentSliceId = 0; int32 totalSlices = static_cast((inputValues.Zend - inputValues.Zstart) / inputValues.SliceResolution); // Loop over each slice that generated a polygon for the outline of the mesh @@ -176,36 +267,55 @@ Result<> RegularGridSampleSurfaceMesh::operator()() continue; } - // We should probably parallelize over each slice - // Loop over all edges and find the edges that are just for the current Slice Id - for(usize edgeIdx = 0; edgeIdx < numEdges; edgeIdx++) + if(true) { - int32 sliceIndex = sliceId[edgeIdx]; - if(currentSliceId == sliceIndex) - { - edgeIndices.push_back(edgeIdx); - } + taskRunner.execute(SampleSurfaceMeshSliceImpl(this, edgeGeom, currentSliceId, possibleIndex.value(), imageGeom, sliceId, featureIds, m_ShouldCancel)); } - - // Now that we have the edges that are on this slice, iterate over all - // voxels on this slice - size_t imageGeomIdx = possibleIndex.value(); - int32 hitCount = 0; - for(size_t planeIdx = 0; planeIdx < cellsPerSlice; planeIdx++) + else { - Point3Df imagePoint = imageGeom.getCoordsf(imageGeomIdx + planeIdx); - imagePoint[2] = 0.0f; // Force this down to the zero plane. - - if(pointInPolygon(edgeGeom, edgeIndices, imagePoint, verts, edges)) - { - featureIds[imageGeomIdx + planeIdx] = 1; - hitCount++; - } + // We should probably parallelize over each slice + // Loop over all edges and find the edges that are just for the current Slice Id + // for(usize edgeIdx = 0; edgeIdx < numEdges; edgeIdx++) + // { + // int32 sliceIndex = sliceId[edgeIdx]; + // if(currentSliceId == sliceIndex) + // { + // edgeIndices.push_back(edgeIdx); + // } + // } + // + // // Now that we have the edges that are on this slice, iterate over all + // // voxels on this slice + // size_t imageGeomIdx = possibleIndex.value(); + // for(size_t planeIdx = 0; planeIdx < cellsPerSlice; planeIdx++) + // { + // Point3Df imagePoint = imageGeom.getCoordsf(imageGeomIdx + planeIdx); + // imagePoint[2] = 0.0f; // Force this down to the zero plane. + // + // if(pointInPolygon(edgeGeom, edgeIndices, imagePoint, verts, edges)) + // { + // featureIds[imageGeomIdx + planeIdx] = 1; + // } + // } + // edgeIndices.clear(); + // edgeIndices.reserve(1024); } - // fmt::print("[{}] edgeIndices.size(): {} Voxels in Polygon: {}\n", currentSliceId, edgeIndices.size(), hitCount); - edgeIndices.clear(); - edgeIndices.reserve(1024); currentSliceId++; } + + taskRunner.wait(); // This will spill over if the number of DataArrays to process does not divide evenly by the number of threads. + return {}; } + +// ----------------------------------------------------------------------------- +void RegularGridSampleSurfaceMesh::sendThreadSafeUpdate(Int32Array& featureIds, const std::vector& rasterBuffer, usize offset) +{ + // We lock access to the DataArray since I don't think DataArray is thread safe. + std::lock_guard lock(m_ProgressMessage_Mutex); + auto& dataStore = featureIds.getDataStoreRef(); + for(usize idx = 0; idx < rasterBuffer.size(); idx++) + { + dataStore[offset + idx] = rasterBuffer[idx]; + } +} diff --git a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.hpp b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.hpp index 3b050054d6..5bfcf62fb4 100644 --- a/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.hpp +++ b/src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.hpp @@ -43,6 +43,7 @@ class SIMPLNXCORE_EXPORT RegularGridSampleSurfaceMesh : public SampleSurfaceMesh Result<> operator()(); const std::atomic_bool& getCancel(); + void sendThreadSafeUpdate(Int32Array& m_FeatureIds, const std::vector& rasterBuffer, usize offset); protected: void generatePoints(std::vector& points) override; @@ -52,5 +53,11 @@ class SIMPLNXCORE_EXPORT RegularGridSampleSurfaceMesh : public SampleSurfaceMesh const RegularGridSampleSurfaceMeshInputValues* m_InputValues = nullptr; const std::atomic_bool& m_ShouldCancel; const IFilter::MessageHandler& m_MessageHandler; + + // Thread safe Progress Message + mutable std::mutex m_ProgressMessage_Mutex; + usize m_ProgressCounter = 0; + usize m_LastProgressInt = 0; + std::chrono::steady_clock::time_point m_InitialTime = std::chrono::steady_clock::now(); }; } // namespace nx::core diff --git a/src/simplnx/Utilities/GeometryUtilities.cpp b/src/simplnx/Utilities/GeometryUtilities.cpp index 5132f113c4..9369971728 100644 --- a/src/simplnx/Utilities/GeometryUtilities.cpp +++ b/src/simplnx/Utilities/GeometryUtilities.cpp @@ -647,8 +647,6 @@ GeometryUtilities::SliceTriangleReturnType GeometryUtilities::SliceTriangleGeome // Create the plane slice_helper::Plane plane(planeNormal, pointOnPlane); - int32 localEdgeCount = 0; - std::set uniqueEdges; // Loop over each Triangle and get edges/vertices of any intersection for(usize triIdx = 0; triIdx < numTris; triIdx++) { @@ -658,20 +656,10 @@ GeometryUtilities::SliceTriangleReturnType GeometryUtilities::SliceTriangleGeome { regionId = triRegionIdPtr->operator[](triIdx); } - - // std::array faceVertices; - // triangle.getFaceCoordinates(triIdx, faceVertices); - std::array faceVertices = GetFaceCoordinates(triIdx, triVertStore, triEdgeStore); // Compute the intersection slice_helper::Edge intersectionEdge = IntersectTriangleWithPlane(faceVertices[0], faceVertices[1], faceVertices[2], plane); - - if(intersectionEdge.valid) - { - uniqueEdges.insert(intersectionEdge); - } - if(intersectionEdge.valid) { slicedVerts.push_back(intersectionEdge.start[0]); @@ -688,7 +676,6 @@ GeometryUtilities::SliceTriangleReturnType GeometryUtilities::SliceTriangleGeome regionIds.push_back(regionId); } edgeCounter++; - localEdgeCount++; } } sliceIndex++;