Pipeline confusion?

[Theoretialperson]

I’m having facing a problem that possibly stems for a misunderstanding of how pipelines work and get updated. I have the following piece of code,

labelfile = labelfilelist[0]                                                         
labelpath = os.path.join(label_path, labelfile)                                      
MeshType = itk.Mesh[itk.SS,3]                                                        
meshReader = itk.MeshFileReader[MeshType].New()                                      
meshReader.SetFileName(labelpath)                                                    
                                                                                     
pad_filter = itk.ConstantPadImageFilter[ImageType, ImageType].New()                  
pad_filter.SetInput(reader.GetOutput())                                              
pad_filter.SetPadLowerBound((1,1,1))                                                 
pad_filter.SetPadUpperBound((1,1,1))                                                 
                                                                                     
mesh_to_image_filter = itk.TriangleMeshToBinaryImageFilter[MeshType, ImageType].New()
mesh_to_image_filter.SetInput(meshReader.GetOutput())                                
mesh_to_image_filter.SetInfoImage(pad_filter.GetOutput())                            
                                                                                     
crop_filter = itk.CropImageFilter[ImageType, ImageType].New()                        
crop_filter.SetInput(mesh_to_image_filter.GetOutput())                               
crop_filter.SetUpperBoundaryCropSize((1,1,1))                                        
crop_filter.SetLowerBoundaryCropSize((1,1,1))                                        
crop_filter.Update()                                                                 
                                                                                     
cropped_image = crop_filter.GetOutput()                                              
image = itk.GetArrayViewFromImage(cropped_image).astype(np.short)                    
itk_image = itk.GetImageFromArray(image)                                             
itk_image.SetMetaDataDictionary(reader.GetMetaDataDictionary())                      
                                                                                     
header = itk.ChangeInformationImageFilter[ImageType].New()                           
header.SetInput(itk_image)                                                           
header.SetOutputSpacing(reader.GetOutput().GetSpacing())                             
header.ChangeSpacingOn()                                                             
header.SetOutputOrigin(reader.GetOutput().GetOrigin())                               
header.ChangeOriginOn()                                                              
header.SetOutputDirection(reader.GetOutput().GetDirection())                         
header.ChangeDirectionOn()                                                           
header.UpdateOutputInformation()                                                     
                                                                                     
writer = itk.ImageFileWriter[ImageType].New()                                        
writer.SetInput(header.GetOutput())                                                  
segmentation_name = os.path.join(output_label, series_name + '_' +ann+ '.nii.gz')    
writer.SetFileName(segmentation_name)                                                
writer.Update()                                                                      

which gives me an error:

Traceback (most recent call last):                                                                                                
  File "dicom_to_nifti.py", line 276, in <module>                                                                                 
    separate_labels(spec_data_dar, spec_annotation_directory, output, case_number, phase)                                         
  File "dicom_to_nifti.py", line 127, in separate_labels                                                                          
    crop_filter.Update()                                                                                                          
RuntimeError: /work/ITK-source/ITK/Modules/Filtering/ImageGrid/include/itkCropImageFilter.hxx:76:                                 
itk::ERROR: itk::ERROR: CropImageFilter(0x55577c544990): The input image's size [0, 0, 0] is less than the total of the crop size!

If I Update() the mesh_to_image_filter before passing its output to the crop_filter, it does not complain anymore, but I don’t understand why I have to in the first place.

[dzenanz (Dženan Zukić)]

pad_filter needs to be updated, in order to have the output image’s metadata (size, start index etc) properly computed. As mesh_to_image_filter depends on it, it calls the update in its update method.

I guess that crop_filter tries to calculate image metadata before ensuring its input is up to date, thus causing the problem.

[Theoretialperson]

I guess this is where the error originates from:

template <typename TInputImage, typename TOutputImage>
void
CropImageFilter<TInputImage, TOutputImage>::VerifyInputInformation() ITKv5_CONST
{
  Superclass::VerifyInputInformation();

  const TInputImage * inputPtr = this->GetInput();

  InputImageSizeType input_sz = inputPtr->GetLargestPossibleRegion().GetSize();

  for (unsigned int i = 0; i < InputImageDimension; ++i)
  {
    if (input_sz[i] < (m_UpperBoundaryCropSize[i] + m_LowerBoundaryCropSize[i]))
    {
      itkExceptionMacro("The input image's size " << input_sz << " is less than the total of the crop size!");
    }
  }
}

Wouldn’t it be a bug if I have to explicitly call Update?

Also, as you see, I’m padding and then cropping by 1 pixel. This is so that I can bypass a bug in TriangleMeshToBinaryImageFilter. Again, I would expect that in reality the padding and cropping would not affect the result of TriangleMeshToBinaryImageFilter, but as you can see in the picture, the output seems shifted (red being the one where I cropped and padded, and green without cropping and padding).

[dzenanz (Dženan Zukić)]

VerifyInputInformation isn’t supposed to call Update. It might/should call inputPtr->UpdateOutputInformation();. Does this fix it?

TriangleMeshToBinaryImageFilter might assume that LargestPossibleRegion starts with index 0. In this case it is wrong. Can you check that this is the problem?

[Theoretialperson]

Sorry, I’m not compiling ITK, I’m just using the Python version. I was looking at the code on Github.

It does not seem to ignore the starting index:

/** Update */
template <typename TInputMesh, typename TOutputImage>
void
TriangleMeshToBinaryImageFilter<TInputMesh, TOutputImage>::GenerateData()
{
  itkDebugMacro(<< "TriangleMeshToBinaryImageFilter::Update() called");

  // Get the input and output pointers
  OutputImagePointer OutputImage = this->GetOutput();
  if (m_InfoImage == nullptr)
  {
    if (m_Size[0] == 0 || m_Size[1] == 0 || m_Size[2] == 0)
    {
      itkExceptionMacro(<< "Must Set Image Size");
    }

    typename OutputImageType::RegionType region;

    region.SetSize(m_Size);
    region.SetIndex(m_Index);

    OutputImage->SetLargestPossibleRegion(region); //
    OutputImage->SetBufferedRegion(region);        // set the region
    OutputImage->SetRequestedRegion(region);       //
    OutputImage->SetSpacing(m_Spacing);            // set spacing
    OutputImage->SetOrigin(m_Origin);              //   and origin
    OutputImage->SetDirection(m_Direction);        // direction cosines
  }
  else
  {
    itkDebugMacro(<< "Using info image");
    m_InfoImage->Update();
    OutputImage->CopyInformation(m_InfoImage);
    OutputImage->SetRegions(m_InfoImage->GetLargestPossibleRegion());
    m_Size = m_InfoImage->GetLargestPossibleRegion().GetSize();
    m_Index = m_InfoImage->GetLargestPossibleRegion().GetIndex();
    m_Spacing = m_InfoImage->GetSpacing();
    m_Origin = m_InfoImage->GetOrigin();
    m_Direction = m_InfoImage->GetDirection();
  }

[dzenanz (Dženan Zukić)]

Can you submit these as two separate bugs on the issue tracker?

ITKCuberille is kind of a higher-quality TriangleMeshToBinaryImageFilter. That might alleviate your immediate issues.

[Theoretialperson]

Thanks I’ll have a look.

Here are the issues I created:

TriangleMeshToBinaryImageFilter not giving correct result after padding image and cropping output

CropImageFilter complaints about 0 image size

[Theoretialperson]

@dzenanz From my understanding, ITKCuberille does the inverse of what I want to do, i.e. it is a binary image / implicit surface to mesh filter.

Are there any alternatives or am is e.g. VTK my best alternative?

[dzenanz (Dženan Zukić)]

@dzenanz

TriangleMeshToBinaryImageFilter

This is the only filter in ITK. VTK has vtkPolyDataToImageStencil. There might be some other implementations, but I think the one in VTK is quite mature and heavily used.