How to use ITK multi threading primitives to do parallel execution of loop iteration

br_cpvc · November 30, 2022, 9:21am

I am currently working on a patch which accelerates the execution speed of the ANTs Atropos application. I have until now been using openmp to implement the acceleration. The suggested code changes and initial related discussion can be found in PR https://github.com/ANTsX/ANTs/pull/1452

In essence the change accelerates a loop using the omp parallel c++ pragma macro.

In the PR it was suggested by @gdevenyi (thank you) that I reimplement the code using ITK’s multithreading features as described in https://itk.org/ItkSoftwareGuide.pdf

I have been trying to wrap my head around the ITK multi-threading framework, but cannot seem to find a way to convert my current openmp accelerated code https://github.com/ANTsX/ANTs/commit/4c7f75a3f4a56e22387bf09097b983515ad854a5 into code that utilizes to ITK building blocks described in the documentation.

If anyone have any ideas, suggestion, or resources such as examples or tutorials, it will be of greatly help.

dzenanz · November 30, 2022, 3:52pm

Here are a few examples of itk::MultiThreaderBase::ParallelizeArray in action:

github.com

InsightSoftwareConsortium/ITK/blob/777e66d45981f287a88a7f63e645b637588288f5/Modules/Filtering/Path/include/itkContourExtractor2DImageFilter.hxx#L382-L390


      
          itk::MultiThreaderBase::Pointer mt = this->GetMultiThreader();
          mt->ParallelizeArray(
            0,
            allLabels.size(),
            [this, &allLabels, &input, &labelsRegions, totalPixelCount, &labelsContoursOutput](SizeValueType i) -> void {
              const InputPixelType label{ allLabels[i] };
              this->CreateSingleContour(label, input, labelsRegions[label], totalPixelCount, labelsContoursOutput[label]);
            },
            nullptr);

github.com

InsightSoftwareConsortium/ITK/blob/23da3dfa79152648f3d38050ee147df47f41ce0a/Modules/Core/Common/test/itkMultiThreaderParallelizeArrayTest.cxx#L73-L74


      
          mt->ParallelizeArray(
            1, size, [&vec](int i) { vec[i] = i; }, progressPO);

github.com

InsightSoftwareConsortium/ITK/blob/23da3dfa79152648f3d38050ee147df47f41ce0a/Modules/Segmentation/SuperPixel/include/itkSLICImageFilter.hxx#L694-L695


      
          this->GetMultiThreader()->ParallelizeArray(
            0, numberOfClusters, [this](SizeValueType idx) { this->ThreadedConnectivity(idx); }, this);

github.com

InsightSoftwareConsortium/ITK/blob/b9005b83b315cdf61a345037c49dfd84f7bde386/Modules/Registration/Metricsv4/include/itkManifoldParzenWindowsPointSetFunction.hxx#L53-L67


      
          m_MultiThreader->ParallelizeArray(
            0,
            points->Size(),
            [&](SizeValueType index) {
              CovarianceMatrixType covariance(PointDimension, PointDimension);
          
          
    covariance.SetIdentity();
              covariance *= this->m_KernelSigma;
          
          
    inputGaussians[index] = GaussianType::New();
              inputGaussians[index]->SetMeasurementVectorSize(PointDimension);
              inputGaussians[index]->SetMean(points->ElementAt(index));
              inputGaussians[index]->SetCovariance(covariance);
            },
            nullptr);

br_cpvc · December 5, 2022, 10:06am

Hi @dzenanz. Thank you for your reply, I really appreciate your help and the pointers to the code examples. I have tried my best to think of ways to apply the examples to my problem, but I still fail to see how. As fare as I can see, the example that you have linked do parallel acceleration of an integer index iterationer, and when I try to apply this it seems hard to make a generel solution that works for multiple dimensionalties.

Maybe a parallelized version of the Iterator’s, or some sort of wrapper function exists and can be applied?

I’m trying to replace/accelerate the loop iteration:

ConstNeighborhoodIterator<ClassifiedImageType> ItO(
  radius, this->GetOutput(), this->GetOutput()->GetRequestedRegion());
ImageRegionIterator<RealImageType> ItS(this->m_SumPosteriorProbabilityImage,
                                       this->m_SumPosteriorProbabilityImage->GetRequestedRegion());
for (ItO.GoToBegin(), ItS.GoToBegin(); !ItO.IsAtEnd(); ++ItO, ++ItS)

in https://github.com/ANTsX/ANTs/blob/master/ImageSegmentation/antsAtroposSegmentationImageFilter.hxx#L1991-L1995

I have been looking at the ParallelizeImageRegion construct as used in useparallelizeimageregion maybe this construct can be used? I however also for this construct fail to see how the code should be put together (how ParallelizeImageRegion should be used together with the ConstNeighborhoodIterator and ImageRegionIterator).

Do you have any idears or thoughts about how and/or if this sort of acceleration of the ConstNeighborhoodIterator and ImageRegionIterator can be achieved using the ITK primitives?

dzenanz · December 5, 2022, 2:43pm

I thought you were trying to parallelize for (unsigned int k = 0; k < zSize; k++).

For parallelizing region iteration, ParallelizeImageRegion is the right choice. In essence, the old code goes into the lambda function of ParallelizeImageRegion, and it iterates over the lambda’s parameter (partial region) instead of the whole region. Some other changes might be needed to enable parallel computation. Here is a simpler example:

github.com

InsightSoftwareConsortium/ITK/blob/23da3dfa79152648f3d38050ee147df47f41ce0a/Modules/Segmentation/SuperPixel/include/itkSLICImageFilter.hxx#L636-L641


      
          this->GetMultiThreader()->template ParallelizeImageRegion<ImageDimension>(
            outputImage->GetRequestedRegion(),
            [this](const OutputImageRegionType & outputRegionForThread) {
              this->ThreadedUpdateClusters(outputRegionForThread);
            },
            this);

Here is a non-trivial example:

github.com

InsightSoftwareConsortium/ITK/blob/23da3dfa79152648f3d38050ee147df47f41ce0a/Modules/Segmentation/ConnectedComponents/include/itkRelabelComponentImageFilter.hxx#L231-L267


      
          // In parallel apply the relabling map
          this->GetMultiThreader()->template ParallelizeImageRegion<ImageDimension>(
            output->GetRequestedRegion(),
            [this, &relabelMap](const RegionType & outputRegionForThread) {
              auto                  outputRequestedRegion = this->GetOutput()->GetRequestedRegion();
              TotalProgressReporter report(this, outputRequestedRegion.GetNumberOfPixels(), 100, 0.5f);
          
          
    ImageScanlineIterator<OutputImageType>     oit(this->GetOutput(), outputRegionForThread);
              ImageScanlineConstIterator<InputImageType> it(this->GetInput(), outputRegionForThread);
          
          
    auto mapIt = relabelMap.cbegin();
          
          
    while (!oit.IsAtEnd())
              {
                while (!oit.IsAtEndOfLine())
                {
                  const auto && inputValue = it.Get();
          
          
        if (mapIt->first != inputValue)
                  {

This file has been truncated. show original

and the refactoring change from the old to the new, parallel implementation:

github.com/InsightSoftwareConsortium/ITK

PERF: Refactor RelabelComponentImageFilter

committed 02:20AM - 25 Jan 20 UTC

blowekamp

+216 -151

Parallelized the computation of the label counting with parallel reduction of th…e map. Parallelized applying the output production step. Modernizes the use of STL data structures and loops. Stream-lined cases. Save iterators used for searching to exploit coherence in images. Updated some type to avoid casts. The filter clearly maps input pixels to output pixels, so those type are used. And the number of objects type has been changed to SizeValueType.

Also take a look at ITK5MigrationGuide.

br_cpvc · December 6, 2022, 1:32pm

Hi @dzenanz Thank your so much, this was exactly what I needed to get going, with offset in the https://github.com/InsightSoftwareConsortium/ITK/blob/23da3dfa79152648f3d38050ee147df47f41ce0a/Modules/Segmentation/ConnectedComponents/include/itkRelabelComponentImageFilter.hxx#L231-L267 example, I have managed to rework the code into something that gives the same both output and execution acceleration as in my original openmp implementation, the new code however is a lot more generic and should also work for all dimensions.

I have made the following branch with my changes: https://github.com/ANTsX/ANTs/compare/master…br-cpvc:ANTs:upstream_20221206_1415

As this is my first go at using the ITK multithreading primitives I’m not entirely sure that my implementation is ITK best practice. If you could glans over these changes and comment on the implementation that I have made, I would be of great value. Any suggestions for improvements are very welcome.

dzenanz · December 6, 2022, 2:20pm

From a glance, you got it right. It is supposed to be that easy.

br_cpvc · December 7, 2022, 9:34am

@dzenanz Thank for so much for taking the time to do a final review and comment on my changes.

I’m glad to inform you that I now have submitted this change as a new PR for the ANTs project: PERF: accelerating antsAtroposSegmentationImageFilter.hxx function GetPosteriorProbabilityImage using ITK ParallelizeImageRegion by br-cpvc · Pull Request #1462 · ANTsX/ANTs · GitHub My hope is that this can be merged, so that every one can benefit from this.

Thank you again for helping me wrap my head around the ITK multi-threading primitives and realize how this should be implemented.

Also special thanks to @gdevenyi for sending me here with my questions.

For everyone in the projects: keep up the great work that you all are doing, thanks a bunch, I’m very grateful