diff --git a/Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.h b/Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.h index c5c7926610f..ab274096095 100644 --- a/Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.h +++ b/Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.h @@ -24,6 +24,7 @@ #include #include #include +#include namespace itk { @@ -310,6 +311,7 @@ class ITK_TEMPLATE_EXPORT ContourExtractor2DImageFilter VertexToContourMap m_ContourStarts; VertexToContourMap m_ContourEnds; SizeValueType m_NumberOfContoursCreated = 0; + InputRealType m_Isovalue = 0; }; void @@ -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> 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 diff --git a/Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.hxx b/Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.hxx index c4c4e570d50..07d42e56eb8 100644 --- a/Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.hxx +++ b/Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.hxx @@ -20,6 +20,7 @@ #include +#include "itkConstantBoundaryCondition.h" #include "itkConstShapedNeighborhoodIterator.h" #include "itkTotalProgressReporter.h" #include "itkContourExtractor2DImageFilter.h" @@ -43,8 +44,7 @@ template void ContourExtractor2DImageFilter::GenerateData() { - m_NumberOutputsAllocated = 0; - m_NumberOutputsWritten = 0; + m_Paths.clear(); if (m_LabelContours) // each label has one or more contours { @@ -52,7 +52,6 @@ ContourExtractor2DImageFilter::GenerateData() } else // simple case of a single iso-value { - m_NumberLabelsRemaining = 1; m_Interpolate = true; const InputRegionType region{ this->GetInput()->GetRequestedRegion() }; @@ -68,9 +67,21 @@ ContourExtractor2DImageFilter::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); + } } } @@ -84,6 +95,7 @@ ContourExtractor2DImageFilter::CreateSingleContour(const InputImage SizeValueType totalNumberOfPixels) { ContourData contourData; + contourData.m_Isovalue = lowerIsovalue * 0.5 + upperIsovalue * 0.5; TotalProgressReporter progress(this, totalNumberOfPixels); @@ -103,13 +115,23 @@ ContourExtractor2DImageFilter::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; + + using BoundaryConditionType = ConstantBoundaryCondition; + BoundaryConditionType boundaryCondition; + boundaryCondition.SetConstant(lowerIsovalue); // a value distinct from the current label + + using SquareIterator = ConstShapedNeighborhoodIterator; 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 *>((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); @@ -290,7 +312,6 @@ ContourExtractor2DImageFilter::GenerateDataForLabels() const typename std::vector::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). @@ -350,61 +371,25 @@ ContourExtractor2DImageFilter::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(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(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(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(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); } @@ -585,25 +570,10 @@ template void ContourExtractor2DImageFilter::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(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(this->MakeOutput(m_NumberOutputsWritten).GetPointer()); - this->SetNthOutput(m_NumberOutputsWritten, output.GetPointer()); - } typename VertexListType::Pointer path{ const_cast(output->GetVertexList()) }; path->Initialize(); // use std::vector version of 'reserve()' instead of @@ -632,7 +602,9 @@ ContourExtractor2DImageFilter::FillOutputs(ContourData & contourDat ++itC; } } - output->Modified(); + + std::lock_guard pathProtector(m_PathsMutex); + m_Paths[contourData.m_Isovalue].push_back(output); } }