SphinxExamples/src/Core/Common/FilterAndParallelizeImageRegion/Code.cxx

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
 
#include "itkLogImageFilter.h"
#include "itkRandomImageSource.h"
#include "itkImageDuplicator.h"
#include "itkMultiThreaderBase.h"
#include "itkImageRegionIterator.h"
 
constexpr unsigned int Dimension = 2;
using PixelType = unsigned int;
using ImageType = itk::Image<PixelType, Dimension>;
 
// calculate log(1+x), where x is pixel value, using LogImageFilter
void
log1xViaLogImageFilter(ImageType::Pointer & image)
{
  // LogImageFilter calculates log(x), so we have to modify the image first
  // by increase its every pixel value by 1, and then apply log filter to it
  itk::ImageRegionIterator<ImageType> it(image, image->GetBufferedRegion());
  for (; !it.IsAtEnd(); ++it)
  {
    it.Set(1 + it.Get());
  }
 
  // classic filter declaration and invocation
  using LogType = itk::LogImageFilter<ImageType, ImageType>;
  auto logF = LogType::New();
  logF->SetInput(image);
  logF->SetInPlace(true);
  logF->Update();
  image = logF->GetOutput();
  image->DisconnectPipeline();
}
 
// calculate log(1+x), where x is pixel value, using ParallelizeImageRegion
void
log1xViaParallelizeImageRegion(ImageType::Pointer & image)
{
  itk::MultiThreaderBase::Pointer mt = itk::MultiThreaderBase::New();
  // ParallelizeImageRegion invokes the provided lambda function in parallel
  // each invocation will contain a piece of the overall region
  mt->ParallelizeImageRegion<Dimension>(
    image->GetBufferedRegion(),
    // here we creat an ad-hoc lambda function to process the region pieces
    // the lambda will have access to variable 'image' from the outer function
    // it will have parameter 'region', which needs to be processed
    [image](const ImageType::RegionType & region) {
      itk::ImageRegionIterator<ImageType> it(image, region);
      for (; !it.IsAtEnd(); ++it)
      {
        it.Set(std::log(1 + it.Get()));
      }
    },
    nullptr); // we don't have a filter whose progress needs to be updated
}
 
int
main()
{
  int result = EXIT_SUCCESS;
 
  // create an image
  ImageType::RegionType region = { { 0, 0 }, { 50, 20 } }; // indices zero, size 50x20
  using RandomSourceType = itk::RandomImageSource<ImageType>;
  auto randomImageSource = RandomSourceType::New();
  randomImageSource->SetSize(region.GetSize());
  // we don't want overflow on 1+x operation, so set max pixel value
  randomImageSource->SetMax(itk::NumericTraits<PixelType>::max() - 1);
  randomImageSource->SetNumberOfWorkUnits(1); // to produce deterministic results
  randomImageSource->Update();
 
  ImageType::Pointer image = randomImageSource->GetOutput();
  image->DisconnectPipeline();
 
  // create another image, to be passed to the alternative method
  using DuplicatorType = itk::ImageDuplicator<ImageType>;
  auto duplicator = DuplicatorType::New();
  duplicator->SetInputImage(image);
  duplicator->Update();
 
  ImageType::Pointer clonedImage = duplicator->GetOutput();
  clonedImage->DisconnectPipeline();
 
  // invoke the two functions
  log1xViaLogImageFilter(image);
  log1xViaParallelizeImageRegion(clonedImage);
 
  // compare to make sure the results are the same
  unsigned                                 diffCount = 0;
  itk::ImageRegionConstIterator<ImageType> it1(image, region);
  itk::ImageRegionConstIterator<ImageType> it2(clonedImage, region);
  for (; !it1.IsAtEnd(); ++it1, ++it2)
  {
    if (it1.Get() != it2.Get())
    {
      std::cerr << "Pixel values are different at index " << it1.GetIndex() << it1.Get() << " vs. " << it2.Get()
                << std::endl;
      //<< "\n\tlog1xViaLogImageFilter's value: " << it1.Get()
      //<< "\n\tlog1xViaParallelizeImageRegion: " << it2.Get() << std::endl;
      diffCount++;
      result = EXIT_FAILURE;
    }
  }
 
  if (diffCount == 0)
  {
    std::cout << "LogImageFilter and ParallelizeImageRegion generate the same result." << std::endl;
  }
  else
  {
    std::cout << "Discrepancy! " << diffCount << " pixels out of " << region.GetNumberOfPixels() << " are different."
              << std::endl;
  }
  return result;
}