Skip to content

Commit

Permalink
ENH: Parallelize the codes over the z slice being rastered.
Browse files Browse the repository at this point in the history 
Signed-off-by: Michael Jackson <[email protected]>
  • Loading branch information
@imikejackson
imikejackson committed Nov 28, 2024
1 parent 6363ea9 commit 7711bfc
Show file tree
Hide file tree
Showing 3 changed files with 154 additions and 50 deletions.
184 changes: 147 additions & 37 deletions src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<nx::core::Point3Df, 2> GetEdgeCoordinates(usize edgeId, INodeGeometry0D::SharedVertexList& verts, INodeGeometry1D::SharedEdgeList& edges)
inline std::array<nx::core::Point3Df, 2> GetEdgeCoordinates(usize edgeId, const INodeGeometry0D::SharedVertexList& verts, const INodeGeometry1D::SharedEdgeList& edges)
{
usize v0Idx = edges[edgeId * 2];
usize v1Idx = edges[edgeId * 2 + 1];
Expand All @@ -23,7 +25,8 @@ inline std::array<nx::core::Point3Df, 2> 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<usize>& edgeIndices, const Point3Df& point, INodeGeometry0D::SharedVertexList& verts, INodeGeometry1D::SharedEdgeList& edges)
bool pointInPolygon(const EdgeGeom& edgeGeom, const std::vector<usize>& edgeIndices, const Point3Df& point, const INodeGeometry0D::SharedVertexList& verts,
const INodeGeometry1D::SharedEdgeList& edges)
{
size_t intersections = 0;
size_t numEdges = edgeIndices.size();
Expand Down Expand Up @@ -57,6 +60,97 @@ bool pointInPolygon(const EdgeGeom& edgeGeom, const std::vector<usize>& 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<usize> 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<int32> 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<std::chrono::milliseconds>(now - start).count() > 1000)
// {
// const float progress = static_cast<float>(i) / static_cast<float>(m_TotalPoints) * 100.0f;
// m_FilterAlg->updateProgress(fmt::format("Processing {}: {}% completed", m_DataArrayPtr->getName(), static_cast<int>(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

// -----------------------------------------------------------------------------
Expand Down Expand Up @@ -106,13 +200,6 @@ void RegularGridSampleSurfaceMesh::generatePoints(std::vector<Point3Df>& 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;
Expand Down Expand Up @@ -140,21 +227,25 @@ Result<> RegularGridSampleSurfaceMesh::operator()()
auto& sliceId = m_DataStructure.getDataRefAs<Int32Array>(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<ImageGeom>(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<Int32Array>(m_InputValues->FeatureIdsArrayPath)->getDataStoreRef();
auto featureIds = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->FeatureIdsArrayPath); //->getDataStoreRef();

std::vector<usize> edgeIndices;
edgeIndices.reserve(1024); // Reserve some space in the vector. This is just a guess.

auto nthreads = static_cast<int32>(std::thread::hardware_concurrency()); // Returns ZERO if not defined on this platform
ParallelTaskAlgorithm taskRunner;
taskRunner.setParallelizationEnabled(true);

int32 currentSliceId = 0;
int32 totalSlices = static_cast<int32>((inputValues.Zend - inputValues.Zstart) / inputValues.SliceResolution);
// Loop over each slice that generated a polygon for the outline of the mesh
Expand All @@ -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<int32>& rasterBuffer, usize offset)
{
// We lock access to the DataArray since I don't think DataArray is thread safe.
std::lock_guard<std::mutex> lock(m_ProgressMessage_Mutex);
auto& dataStore = featureIds.getDataStoreRef();
for(usize idx = 0; idx < rasterBuffer.size(); idx++)
{
dataStore[offset + idx] = rasterBuffer[idx];
}
}
7 changes: 7 additions & 0 deletions src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<int32>& rasterBuffer, usize offset);

protected:
void generatePoints(std::vector<Point3Df>& points) override;
Expand All @@ -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
13 changes: 0 additions & 13 deletions src/simplnx/Utilities/GeometryUtilities.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -647,8 +647,6 @@ GeometryUtilities::SliceTriangleReturnType GeometryUtilities::SliceTriangleGeome

// Create the plane
slice_helper::Plane plane(planeNormal, pointOnPlane);
int32 localEdgeCount = 0;
std::set<slice_helper::Edge> uniqueEdges;
// Loop over each Triangle and get edges/vertices of any intersection
for(usize triIdx = 0; triIdx < numTris; triIdx++)
{
Expand All @@ -658,20 +656,10 @@ GeometryUtilities::SliceTriangleReturnType GeometryUtilities::SliceTriangleGeome
{
regionId = triRegionIdPtr->operator[](triIdx);
}

// std::array<nx::core::Point3Df, 3> faceVertices;
// triangle.getFaceCoordinates(triIdx, faceVertices);

std::array<nx::core::Point3Df, 3> 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]);
Expand All @@ -688,7 +676,6 @@ GeometryUtilities::SliceTriangleReturnType GeometryUtilities::SliceTriangleGeome
regionIds.push_back(regionId);
}
edgeCounter++;
localEdgeCount++;
}
}
sliceIndex++;
Expand Down

0 comments on commit 7711bfc

Please sign in to comment.