TriangleMeshToBinaryImageFilter not giving correct result after padding image and cropping output

[nicktasios]

Description

Due to a bug with TriangleMeshToBinaryImageFilter where it can get corrupt near the edges, I tried padding the image and then cropping it by 1 pixel border. Unfortunately, that shifts the final results as you can see here:

Green is with padding/cropping, and white is without.

Steps to Reproduce

Run the following code

t_mesh_filter.py                                                                        
import itk                                                                              

if __name__ == '__main__':                                                              
    should_pad_and_crop = True                                                          

    TCoordinate = itk.D                                                                 
    Dimension = 3                                                                       
    TMesh = itk.Mesh[TCoordinate, Dimension].New()                                      
    sphere = itk.RegularSphereMeshSource[TMesh].New()                                   
    sphere.SetResolution(4)                                                             
    sphere_mesh = sphere.GetOutput()                                                    

    mesh_writer = itk.MeshFileWriter[TMesh].New()                                       
    mesh_writer.SetFileName( "sphere.vtk" )                                             
    mesh_writer.SetInput(sphere_mesh)                                                   
    mesh_writer.Update()                                                                

    TMesh = itk.Mesh[itk.SS, Dimension].New()                                           
    mesh_reader = itk.MeshFileReader[TMesh].New()                                       
    mesh_reader.SetFileName("sphere.vtk")                                               
    mesh_reader.Update()                                                                
    sphere_mesh = mesh_reader.GetOutput()                                               

    TPixel = itk.SS                                                                     
    TImage = itk.Image[TPixel, Dimension]                                               

    image = itk.Image[TPixel, Dimension].New()                                          
    region = itk.ImageRegion[Dimension]()                                               
    region.SetSize([128, 128, 128])                                                     
    region.SetIndex([0, 0, 0])                                                          
    image.SetRegions(region)                                                            
    image.Allocate()                                                                    
    image.SetOrigin([-1.28, -1.28, -1.28])                                              
    image.SetSpacing([0.02, 0.02, 0.02])                                                

    if should_pad_and_crop:                                                             
        pad_filter = itk.ConstantPadImageFilter[TImage, TImage].New()                   
        pad_filter.SetInput(image)                                                      
        pad_filter.SetPadLowerBound((1,1,1))                                            
        pad_filter.SetPadUpperBound((1,1,1))                                            
        pad_filter.Update()                                                             

        mesh_to_image_filter = itk.TriangleMeshToBinaryImageFilter[TMesh, TImage].New() 
        mesh_to_image_filter.SetInput(sphere_mesh)                                      
        mesh_to_image_filter.SetInfoImage(pad_filter.GetOutput())                       
        mesh_to_image_filter.Update()                                                   

        crop_filter = itk.CropImageFilter[TImage, TImage].New()                         
        crop_filter.SetInput(mesh_to_image_filter.GetOutput())                          
        crop_filter.SetUpperBoundaryCropSize((1,1,1))                                   
        crop_filter.SetLowerBoundaryCropSize((1,1,1))                                   
        crop_filter.Update()                                                            

        itk.imwrite(crop_filter.GetOutput(), "sphere.nii.gz")                           
    else:                                                                               
        mesh_to_image_filter = itk.TriangleMeshToBinaryImageFilter[TMesh, TImage].New() 
        mesh_to_image_filter.SetInput(sphere_mesh)                                      
        mesh_to_image_filter.SetInfoImage(image)                                        
        mesh_to_image_filter.Update()                                                   

        itk.imwrite(mesh_to_image_filter.GetOutput(), "sphere.nii.gz")                  

Expected behavior

Cropping and padding should not affect the output.

Versions

Python ITK '5.2.1'

Environment

Linux, Python 3.8.10

[thewtex]

@dzenanz any thoughts on how to approach this?

[dzenanz]

I found a workaround (below). I guess the problem is in TriangleMeshToBinaryImageFilter. I will circle back to it later.

import itk

if __name__ == '__main__':
  should_pad_and_crop = True

  TCoordinate = itk.F
  Dimension = 3
  TMesh = itk.Mesh[TCoordinate, Dimension].New()
  sphere = itk.RegularSphereMeshSource[TMesh].New()
  sphere.SetResolution(4)
  sphere_mesh = sphere.GetOutput()

  mesh_writer = itk.MeshFileWriter[TMesh].New()
  mesh_writer.SetFileName("sphere.vtk")
  mesh_writer.SetInput(sphere_mesh)
  mesh_writer.Update()

  TMesh = itk.Mesh[itk.SS, Dimension].New()
  mesh_reader = itk.MeshFileReader[TMesh].New()
  mesh_reader.SetFileName("sphere.vtk")
  mesh_reader.Update()
  sphere_mesh = mesh_reader.GetOutput()

  TPixel = itk.SS
  TImage = itk.Image[TPixel, Dimension]

  image = itk.Image[TPixel, Dimension].New()
  region = itk.ImageRegion[Dimension]()
  region.SetSize([128, 128, 128])
  region.SetIndex([0, 0, 0])
  image.SetRegions(region)
  image.Allocate()
  image.SetOrigin([-1.28, -1.28, -1.28])
  image.SetSpacing([0.02, 0.02, 0.02])

  crop_lower = 20
  crop_upper = 10

  if should_pad_and_crop:
    pad_filter = itk.ConstantPadImageFilter[TImage, TImage].New()
    pad_filter.SetInput(image)
    pad_filter.SetPadLowerBound((crop_lower, crop_lower, crop_lower))
    pad_filter.SetPadUpperBound((crop_upper, crop_upper, crop_upper))
    pad_filter.Update()
    itk.imwrite(pad_filter.GetOutput(), "spherePadded.nrrd")
    spherePadded = itk.imread("spherePadded.nrrd")

    mesh_to_image_filter = itk.TriangleMeshToBinaryImageFilter[TMesh, TImage].New()
    mesh_to_image_filter.SetInput(sphere_mesh)
    mesh_to_image_filter.SetInfoImage(spherePadded)
    mesh_to_image_filter.Update()

    crop_filter = itk.CropImageFilter[TImage, TImage].New()
    crop_filter.SetInput(mesh_to_image_filter.GetOutput())
    crop_filter.SetLowerBoundaryCropSize((crop_lower, crop_lower, crop_lower))
    crop_filter.SetUpperBoundaryCropSize((crop_upper, crop_upper, crop_upper))
    crop_filter.Update()

    itk.imwrite(crop_filter.GetOutput(), "sphere.nrrd")
  else:
    mesh_to_image_filter = itk.TriangleMeshToBinaryImageFilter[TMesh, TImage].New()
    mesh_to_image_filter.SetInput(sphere_mesh)
    mesh_to_image_filter.SetInfoImage(image)
    mesh_to_image_filter.Update()

    itk.imwrite(mesh_to_image_filter.GetOutput(), "sphere.nrrd")

[dzenanz]

C++ repro code, to be invoked as tester sphere.vtk sphere.nrrd 1 (the inputs are created by the Pytho n version above):

#include "itkMesh.h"
#include "itkMeshFileReader.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkTriangleMeshToBinaryImageFilter.h"
#include "itkConstantPadImageFilter.h"
#include "itkCropImageFilter.h"

int
main(int argc, char * argv[])
{
  if (argc != 4)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0];
    std::cerr << " <InputMeshName> <OutputImageName> <cropAndPad>";
    std::cerr << std::endl;
    return EXIT_FAILURE;
  }

  const char * inputMeshName = argv[1];
  const char * outputImageName = argv[2];
  bool         cropAndPad = std::stoul(argv[3]);

  constexpr unsigned int Dimension = 3;
  using PixelType = signed short;

  using MeshType = itk::Mesh<PixelType, Dimension>;

  using MeshReaderType = itk::MeshFileReader<MeshType>;
  MeshReaderType::Pointer meshReader = MeshReaderType::New();
  meshReader->SetFileName(inputMeshName);
  meshReader->Update();
  auto sphereMesh = meshReader->GetOutput();

  using ImageType = itk::Image<PixelType, Dimension>;
  auto                  size = ImageType::SizeType::Filled(128);
  ImageType::RegionType region;
  region.SetSize(size);
  auto image = ImageType::New();
  image->SetRegions(region);
  image->Allocate();
  image->SetOrigin(ImageType::PointType::Filled(-1.28));
  image->SetSpacing(ImageType::SpacingType::Filled(0.02));

  using FilterType = itk::TriangleMeshToBinaryImageFilter<MeshType, ImageType>;
  FilterType::Pointer filter = FilterType::New();
  filter->SetInput(sphereMesh);

  constexpr int padRadius = 50;

  try
  {
    if (cropAndPad)
    {
      using PadType = itk::ConstantPadImageFilter<ImageType, ImageType>;
      auto padder = PadType::New();
      padder->SetInput(image);
      padder->SetPadBound(ImageType::SizeType::Filled(padRadius));
      padder->Update();
      itk::WriteImage(padder->GetOutput(), "spherePadded.nrrd", true);
      image = itk::ReadImage<ImageType>("spherePadded.nrrd");

      filter->SetInfoImage(image); // bug is not triggered
      filter->SetInfoImage(padder->GetOutput()); // bug is triggered
      filter->Update();

      using CropType = itk::CropImageFilter<ImageType, ImageType>;
      auto cropper = CropType::New();
      cropper->SetInput(filter->GetOutput());
      cropper->Update();
      itk::WriteImage(cropper->GetOutput(), outputImageName, true);
    }
    else
    {
      filter->SetInfoImage(image);
      filter->Update();
      itk::WriteImage(filter->GetOutput(), outputImageName, true);
    }
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

[dzenanz]

@nicktasios can you check whether #3009 fixes this for you? What about #2416?

[nicktasios]

Unfortunately I'm having trouble trying to build the python package from source. I followed the instructions and created a new python env, cloned ITKPythonPackage and also ITK, then started the build with:

./build-itk/bin/python setup.py bdist_wheel -- -DITK_SOURCE_DIR:PATH=/home/username/path/ITKPythonPackage/ITK

but I got the following error:

CMake Error at Modules/ThirdParty/TBB/itk-module-init.cmake:1 (find_package):
  Could not find a package configuration file provided by "TBB" with any of  
  the following names:                                                       

    TBBConfig.cmake                                                          
    tbb-config.cmake                                                         

  Add the installation prefix of "TBB" to CMAKE_PREFIX_PATH or set "TBB_DIR" 
  to a directory containing one of the above files.  If "TBB" provides a     
  separate development package or SDK, be sure it has been installed.        
Call Stack (most recent call first):                                         
  CMake/ITKModuleEnablement.cmake:338 (include)                              
  CMakeLists.txt:538 (include)                                                                                                

[dzenanz]

I guess we now require TBB for building binary packages. Can you try latest stable?

[nicktasios]

So I saw that the by default the build script uses the nightly-master branch so I instead tried building without specifying the ITK source. Unfortunately, after hours of building I got the following, not very informative error:

[ 84%] Generating itkForward1DFFTImageFilterPython.cpp, ../../Generators/Python/itk/itkForward1DFFTImageFilterPython.py                                                               
/home/username/path/ITKPythonPackage/_skbuild/linux-x86_64-3.7/cmakebuild/ITKb/Wrapping/Typedefs/itkForward1DFFTImageFilter.i:82: Warning 401: Nothing known about base class 'itk::ImageToImageFilter< itk::Image< double,2 >,itk::Image< std::complex< double >,2 > >'. Ignored.                                                                           
/home/username/path/ITKPythonPackage/_skbuild/linux-x86_64-3.7/cmake-build/ITKb/Wrapping/Typedefs/itkForward1DFFTImageFilter.i:82: Warning 401: Maybe you forgot to instantiate 'itk::ImageToImageFilter< itk::Image< double,2 >,itk::Image< std::complex< double >,2 > >' using %template.                                                                   
gmake[5]: *** [Wrapping/Modules/ITKFFT/CMakeFiles/ITKFFTPython.dir/build.make:1738: Wrapping/Modules/ITKFFT/itkForward1DFFTImageFilterPython.cpp] Error 2                             
gmake[4]: *** [CMakeFiles/Makefile2:25792: Wrapping/Modules/ITKFFT/CMakeFiles/ITKFFTPython.dir/all] Error 2                                                                           
gmake[3]: *** [Makefile:172: all] Error 2                                                                                                                                             
gmake[2]: *** [CMakeFiles/ITK.dir/build.make:132: ITKp/src/ITK-stamp/ITK-build] Error 2                                                                                               
gmake[1]: *** [CMakeFiles/Makefile2:129: CMakeFiles/ITK.dir/all] Error 2                                                                                                              
gmake: *** [Makefile:150: all] Error 2                                                                                                                                                

[dzenanz]

I created an issue #3048 to track this. If @tbirdso does not get to it soon, I will.

[tbirdso]

Thanks for your work here @nicktasios and @dzenanz . I've added some initial discussion on this issue in #3048. @nicktasios, I'm wondering if you have a way of confirming whether the CMake variables ITK_WRAP_DOUBLE and ITK_WRAP_COMPLEX_DOUBLE are both dis/enabled? You could try building with both explicitly turned off (or on):

./build-itk/bin/python setup.py bdist_wheel -- -DITK_SOURCE_DIR:PATH=/home/username/path/ITKPythonPackage/ITK -DITK_WRAP_DOUBLE:BOOL=OFF -DITK_WRAP_COMPLEX_DOUBLE:BOOL=OFF

Several filters in the ITKFFT module require both real and complex pixel types to be enabled in order to be wrapped. Regardless of whether this is turns out to be the root cause, we should definitely be checking both values before we try to wrap affected classes -- I will put in a fix and update #3048.

[dzenanz]

ITK_WRAP_DOUBLE and ITK_WRAP_COMPLEX_DOUBLE

Both are disabled. I also commented in the issue.

[nicktasios]

@tbirdso I will check tomorrow. By the way, if I pass '-j N' to the build command, above, will that enable multithreaded building? I have to say it's quite a slow process at the moment.

[tbirdso]

@nicktasios With @dzenanz 's help we've submitted ITKPythonPackage PR 184 to fix the issue for now -- thank you for helping improve ITK 😄

By the way, if I pass '-j N' to the build command, above, will that enable multithreaded building?

I believe anything after the -- will be passed as a CMake argument (relevant line in the ITKPythonPackage script is here), so yes, please give it a try and enter an issue on the ITKPythonPackage repo if it doesn't work. Unfortunately building Python wrappings is a slow process, but multithreading will help a bit.

[dzenanz]

ninja uses parallelism by default, and is the preferred build system for most people.