Skip to content

Commit

Permalink
ENH: Parallelize the codes over the z slice being raster.
Browse files Browse the repository at this point in the history 
THIS NEEDS TO BE FULLY TESTED
It will probably break the unit test for the Sample SurfaceMesh.

Signed-off-by: Michael Jackson <[email protected]>
  • Loading branch information
@imikejackson
imikejackson committed Nov 28, 2024
1 parent 6363ea9 commit 5dd0c12
Show file tree
Hide file tree
Showing 3 changed files with 126 additions and 57 deletions.
163 changes: 119 additions & 44 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,7 @@ 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 All @@ -134,26 +222,25 @@ Result<> RegularGridSampleSurfaceMesh::operator()()
return result;
}

Check failure on line 223 in src/Plugins/SimplnxCore/src/SimplnxCore/Filters/Algorithms/RegularGridSampleSurfaceMesh.cpp

View workflow job for this annotation

GitHub Actions / clang_format_pr 

code should be clang-formatted [-Wclang-format-violations]


/////////////////////////////////////////////////////////////////////////////
// RASTER THE PIXELS BASED ON POINT IN POLYGON
DataPath edgeAmPath = edgeDataPath.createChildPath(inputValues.EdgeAttributeMatrixName);
DataPath sliceIdDataPath = edgeAmPath.createChildPath(inputValues.SliceIdArrayName);
auto& edgeGeom = m_DataStructure.getDataRefAs<EdgeGeom>(edgeDataPath);
auto& sliceId = m_DataStructure.getDataRefAs<Int32Array>(sliceIdDataPath);
INodeGeometry0D::SharedVertexList& verts = edgeGeom.getVerticesRef();
INodeGeometry1D::SharedEdgeList& edges = edgeGeom.getEdgesRef();
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];

// 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.
ParallelTaskAlgorithm taskRunner;
taskRunner.setParallelizationEnabled(true);

int32 currentSliceId = 0;
int32 totalSlices = static_cast<int32>((inputValues.Zend - inputValues.Zstart) / inputValues.SliceResolution);
Expand All @@ -176,36 +263,24 @@ 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++)
{
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++)
{
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++;
}
}
// 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 5dd0c12

Please sign in to comment.