Skip to content

Parallelize ContourExtractor2DImageFilter #2371

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
dzenanz wants to merge 2 commits into from
Closed

Parallelize ContourExtractor2DImageFilter #2371

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
53c5423
PERF: using constant boundary condition instead of filling the edge
dzenanz Mar 2, 2021
07e880a
ENH: parallelize per label contour extraction
dzenanz Mar 3, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 6 additions & 8 deletions Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
#include <unordered_map>
#include <deque>
#include <list>
#include <mutex>

namespace itk
{
Expand Down Expand Up @@ -310,6 +311,7 @@ class ITK_TEMPLATE_EXPORT ContourExtractor2DImageFilter
VertexToContourMap m_ContourStarts;
VertexToContourMap m_ContourEnds;
SizeValueType m_NumberOfContoursCreated = 0;
InputRealType m_Isovalue = 0;
};

void
Expand All @@ -318,16 +320,12 @@ class ITK_TEMPLATE_EXPORT ContourExtractor2DImageFilter
void
FillOutputs(ContourData & contourData);

// The number of outputs we have allocated capacity for
unsigned int m_NumberOutputsAllocated;

// The number of outputs we have written out so far
unsigned int m_NumberOutputsWritten;
bool m_Interpolate = false; // whether contour positions will be interpolated (yes for single, no for LabelContours)

// The number of labels we have yet to write outputs for
unsigned int m_NumberLabelsRemaining;
// to allow outputting paths in sorted order after parallel computation
std::map<InputRealType, std::vector<OutputPathPointer>> m_Paths;

bool m_Interpolate = false; // whether contour positions will be interpolated (yes for single, no for LabelContours)
std::mutex m_PathsMutex; // to prevent simultanous write access to path map
};
} // end namespace itk

Expand Down
140 changes: 56 additions & 84 deletions Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@

#include <algorithm>

#include "itkConstantBoundaryCondition.h"
#include "itkConstShapedNeighborhoodIterator.h"
#include "itkTotalProgressReporter.h"
#include "itkContourExtractor2DImageFilter.h"
Expand All @@ -43,16 +44,14 @@ template <typename TInputImage>
void
ContourExtractor2DImageFilter<TInputImage>::GenerateData()
{
m_NumberOutputsAllocated = 0;
m_NumberOutputsWritten = 0;
m_Paths.clear();

if (m_LabelContours) // each label has one or more contours
{
this->GenerateDataForLabels();
}
else // simple case of a single iso-value
{
m_NumberLabelsRemaining = 1;
m_Interpolate = true;
const InputRegionType region{ this->GetInput()->GetRequestedRegion() };

Expand All @@ -68,9 +67,21 @@ ContourExtractor2DImageFilter<TInputImage>::GenerateData()
shrunkRegion.GetNumberOfPixels());
}

if (m_NumberOutputsWritten != m_NumberOutputsAllocated)

SizeValueType outputCount = 0;
for (auto const & pathIt : m_Paths)
{
this->SetNumberOfIndexedOutputs(m_NumberOutputsWritten);
outputCount += pathIt.second.size();
}
this->SetNumberOfIndexedOutputs(outputCount);

outputCount = 0;
for (auto const & pathIt : m_Paths) // hashmap
{
for (auto const & path : pathIt.second) // vector
{
this->SetNthOutput(outputCount++, path);
}
}
}

Expand All @@ -84,6 +95,7 @@ ContourExtractor2DImageFilter<TInputImage>::CreateSingleContour(const InputImage
SizeValueType totalNumberOfPixels)
{
ContourData contourData;
contourData.m_Isovalue = lowerIsovalue * 0.5 + upperIsovalue * 0.5;

TotalProgressReporter progress(this, totalNumberOfPixels);

Expand All @@ -103,13 +115,23 @@ ContourExtractor2DImageFilter<TInputImage>::CreateSingleContour(const InputImage
// with the center pixel at the top-left. So we only activate the
// coresponding offsets, and only query pixels 4, 5, 7, and 8 with the
// iterator's GetPixel method.
using SquareIterator = ConstShapedNeighborhoodIterator<InputImageType>;

using BoundaryConditionType = ConstantBoundaryCondition<InputImageType, InputImageType>;
BoundaryConditionType boundaryCondition;
boundaryCondition.SetConstant(lowerIsovalue); // a value distinct from the current label

using SquareIterator = ConstShapedNeighborhoodIterator<InputImageType, BoundaryConditionType>;
const typename SquareIterator::RadiusType radius{ { 1, 1 } };
SquareIterator it(radius, image, shrunkRegion);
const InputOffsetType none{ { 0, 0 } };
const InputOffsetType right{ { 1, 0 } };
const InputOffsetType down{ { 0, 1 } };
const InputOffsetType diag{ { 1, 1 } };

// it.SetBoundaryCondition(boundaryCondition); // doesn't compile
auto nit = static_cast<NeighborhoodIterator<InputImageType, BoundaryConditionType> *>((void *)&it);
nit->SetBoundaryCondition(boundaryCondition);

const InputOffsetType none{ { 0, 0 } };
const InputOffsetType right{ { 1, 0 } };
const InputOffsetType down{ { 0, 1 } };
const InputOffsetType diag{ { 1, 1 } };
it.ActivateOffset(none);
it.ActivateOffset(right);
it.ActivateOffset(down);
Expand Down Expand Up @@ -290,7 +312,6 @@ ContourExtractor2DImageFilter<TInputImage>::GenerateDataForLabels()
const typename std::vector<InputPixelType>::const_iterator last{ std::unique(allLabels.begin(), allLabels.end()) };
allLabels.erase(last, allLabels.end());
}
m_NumberLabelsRemaining = allLabels.size(); // We haven't processed any yet

// Compute bounding box for each label. These will be [inclusive, inclusive] ranges in each coordinate, not
// [inclusive, exclusive).
Expand Down Expand Up @@ -350,61 +371,25 @@ ContourExtractor2DImageFilter<TInputImage>::GenerateDataForLabels()

m_Interpolate = false;

for (SizeValueType i = 0; i < allLabels.size(); i++)
{
const InputPixelType label{ allLabels[i] };
const InputRealType previousLabel = i > 0 ? allLabels[i - 1] : label - 1;
const InputRealType followingLabel = i < allLabels.size() - 1 ? allLabels[i + 1] : label + 1;

// Use the bounding box for this label
const BoundingBox & bbox{ bboxes[label] };
const IndexType min{ bbox.first };
const IndexType max{ bbox.second };
const InputPixelType differentLabel = previousLabel;

// Set boundary values in largerRegion to be distinct from label if they will be looked at.
if (min[0] == left_top[0])
{
const IndexType stripeIndex{ min[0] - 1, min[1] - 1 };
const SizeType stripeSize{ 1, static_cast<SizeValueType>(max[1] - min[1] + 3) };
const InputRegionType stripeRegion{ stripeIndex, stripeSize };
const ImageRegionRange stripeRange{ *largerImage, stripeRegion };
std::fill(stripeRange.begin(), stripeRange.end(), differentLabel);
}
if (min[1] == left_top[1])
{
const IndexType stripeIndex{ min[0] - 1, min[1] - 1 };
const SizeType stripeSize{ static_cast<SizeValueType>(max[0] - min[0] + 3), 1 };
const InputRegionType stripeRegion{ stripeIndex, stripeSize };
const ImageRegionRange stripeRange{ *largerImage, stripeRegion };
std::fill(stripeRange.begin(), stripeRange.end(), differentLabel);
}
if (max[0] == right_bot[0])
{
const IndexType stripeIndex{ max[0] + 1, min[1] - 1 };
const SizeType stripeSize{ 1, static_cast<SizeValueType>(max[1] - min[1] + 3) };
const InputRegionType stripeRegion{ stripeIndex, stripeSize };
const ImageRegionRange stripeRange{ *largerImage, stripeRegion };
std::fill(stripeRange.begin(), stripeRange.end(), differentLabel);
}
if (max[1] == right_bot[1])
{
const IndexType stripeIndex{ min[0] - 1, max[1] + 1 };
const SizeType stripeSize{ static_cast<SizeValueType>(max[0] - min[0] + 3), 1 };
const InputRegionType stripeRegion{ stripeIndex, stripeSize };
const ImageRegionRange stripeRange{ *largerImage, stripeRegion };
std::fill(stripeRange.begin(), stripeRange.end(), differentLabel);
}

// this does not work if labels are floats such as 0.1, 0.23, 0.31, 0.7, etc.
// this->CreateSingleContour(largerImage, labelsRegions[label], label - 0.5, label + 0.5, totalPixelCount);

this->CreateSingleContour(largerImage,
labelsRegions[label],
0.5 * previousLabel + 0.5 * label,
0.5 * label + 0.5 * followingLabel,
totalPixelCount);
}
itk::MultiThreaderBase::Pointer mt = this->GetMultiThreader();
mt->ParallelizeArray(
0,
allLabels.size(),
[this, &allLabels, largerImage, &labelsRegions, totalPixelCount](SizeValueType i) {
const InputPixelType label{ allLabels[i] };
const InputRealType previousLabel = i > 0 ? allLabels[i - 1] : label - 1;
const InputRealType followingLabel = i < allLabels.size() - 1 ? allLabels[i + 1] : label + 1;

// this does not work if labels are floats such as 0.1, 0.23, 0.31, 0.7, etc.
// this->CreateSingleContour(largerImage, labelsRegions[label], label - 0.5, label + 0.5, totalPixelCount);

this->CreateSingleContour(largerImage,
labelsRegions[label],
0.5 * previousLabel + 0.5 * label,
0.5 * label + 0.5 * followingLabel,
totalPixelCount);
},
nullptr);
}


Expand Down Expand Up @@ -585,25 +570,10 @@ template <typename TInputImage>
void
ContourExtractor2DImageFilter<TInputImage>::FillOutputs(ContourData & contourData)
{
--m_NumberLabelsRemaining;
if (m_NumberOutputsWritten + contourData.m_Contours.size() > m_NumberOutputsAllocated)
for (auto it = contourData.m_Contours.begin(); it != contourData.m_Contours.end(); ++it)
{
// We do not have enough capacity; increase capacity to what we need,
// plus a guess of one contour for each unprocessed label.
m_NumberOutputsAllocated = m_NumberOutputsWritten + contourData.m_Contours.size() + m_NumberLabelsRemaining;
this->SetNumberOfIndexedOutputs(m_NumberOutputsAllocated);
}
OutputPathPointer output = dynamic_cast<OutputPathType *>(this->MakeOutput(0).GetPointer());

for (auto it = contourData.m_Contours.begin(); it != contourData.m_Contours.end(); ++it, ++m_NumberOutputsWritten)
{
OutputPathPointer output{ this->GetOutput(m_NumberOutputsWritten) };
if (output.IsNull())
{
// Dynamic cast is OK because we know PathSource will make its templated
// class type
output = dynamic_cast<OutputPathType *>(this->MakeOutput(m_NumberOutputsWritten).GetPointer());
this->SetNthOutput(m_NumberOutputsWritten, output.GetPointer());
}
typename VertexListType::Pointer path{ const_cast<VertexListType *>(output->GetVertexList()) };
path->Initialize();
// use std::vector version of 'reserve()' instead of
Expand Down Expand Up @@ -632,7 +602,9 @@ ContourExtractor2DImageFilter<TInputImage>::FillOutputs(ContourData & contourDat
++itC;
}
}
output->Modified();

std::lock_guard<std::mutex> pathProtector(m_PathsMutex);
m_Paths[contourData.m_Isovalue].push_back(output);
}
}

Expand Down