Hello everyone.
I use ITK (5.4.4) and VTK (9.5.2) in my code. I use some VTK classes (like vtkImageStencil) which use some mesh (vtkPolyData) and i want the result in an ITK image.
To convert result of the vtkImageStencil in an ITK image, i use VTKImageToImageFilter.
When i get the output of the VTKImageToImageFilter in my itk image, it seems that the data buffer is still managed by the vtkImageStencil. Indeed when i delete the vtkImageStencil and after try to access the data buffer of my itk image, the application crash.
I try to duplicate the output of the VTKImageToImageFilter, before filling my ITK Image, it seems to fix my issue.
I see nowhere in the documentation that the user must copy the buffer/the Stencil keep the control of the output buffer.
Am i using the right classes in the right way ?
Is this a lack of detail in the documentation, or a bug in the expected behavior of the VTK classes?
Thank
Here a minimal sample:
#include “itkImage.h”
#include “itkVTKImageToImageFilter.h”
#include <vtkNew.h>
#include <vtkImageStencil.h>
#include <vtkPolyDataToImageStencil.h>
#include <vtkPointData.h>
#include <vtkPolyData.h>
#include <vtkDataArray.h>
#include <vtkImageData.h>
#include <vtkSphereSource.h>
#include // juste use for qDebug
#include “itkMinimumMaximumImageCalculator.h”
class Voxelisation
{
public:
Voxelisation() = default;
~Voxelisation() = default;
void setInputData(vtkPolyData* p_inputData);
void execute(itk::Image<unsigned char, 3>::Pointer& p_outputImage);
std::array<double, 3> m_spacing {0.5, 0.5, 0.5};
unsigned char m_backgroundValue {0};
unsigned char m_foregroundValue {255};
vtkNew<vtkImageStencil> m_vtkImageStencil;
};
void Voxelisation::setInputData(vtkPolyData* p_inputData)
{
vtkNew whiteImage;
std::array<double, 3> origin;
{
std::array<double, 6> bounds;
p_inputData->GetBounds(bounds.data());
whiteImage->SetSpacing(m_spacing.data());
std::array<int, 3> dim;
for (size_t i = 0; i < 3; i++)
{
dim[i] = static_cast(ceil((bounds[i * 2 + 1] - bounds[i * 2]) / m_spacing[i]));
}
whiteImage->SetDimensions(dim.data());
whiteImage->SetExtent(0, dim[0] - 1, 0, dim[1] - 1, 0, dim[2] - 1);
origin[0] = bounds[0] + m_spacing[0] / 2;
origin[1] = bounds[2] + m_spacing[1] / 2;
origin[2] = bounds[4] + m_spacing[2] / 2;
whiteImage->SetOrigin(origin.data());
whiteImage->AllocateScalars(VTK_UNSIGNED_CHAR, 1);
const vtkIdType count = whiteImage->GetNumberOfPoints();
for (vtkIdType i = 0; i < count; ++i)
{
whiteImage->GetPointData()->GetScalars()->SetTuple1(i, m_foregroundValue);
}
}
vtkNew<vtkPolyDataToImageStencil> vtkPolyDataToImageStencil;
{
vtkPolyDataToImageStencil->SetInputData(p_inputData);
vtkPolyDataToImageStencil->SetOutputOrigin(origin.data());
vtkPolyDataToImageStencil->SetOutputSpacing(m_spacing.data());
vtkPolyDataToImageStencil->SetOutputWholeExtent(whiteImage->GetExtent());
vtkPolyDataToImageStencil->Update();
}
m_vtkImageStencil->SetInputData(whiteImage);
m_vtkImageStencil->SetStencilConnection(vtkPolyDataToImageStencil->GetOutputPort());
m_vtkImageStencil->ReverseStencilOff();
m_vtkImageStencil->SetBackgroundValue(m_backgroundValue);
}
void Voxelisation::execute(itk::Image<unsigned char, 3>::Pointer& p_outputImage)
{
m_vtkImageStencil->Update();
auto VTKImageToImageFilter = itk::VTKImageToImageFilter<itk::Image<unsigned char, 3>>::New();
VTKImageToImageFilter->SetInput(m_vtkImageStencil->GetOutput());
VTKImageToImageFilter->UpdateLargestPossibleRegion();
qDebug() << "refCount inside Voxelisation" << VTKImageToImageFilter->GetOutput()->GetReferenceCount(); //1
p_outputImage = VTKImageToImageFilter->GetOutput();
qDebug() << "refCount p_outputImage inside Voxelisation" << p_outputImage->GetReferenceCount(); // 2
if( p_outputImage == nullptr)
{
qDebug() << "outputImage is nullptr";
}
else {
qDebug() << "outputImage is not null"; // not null
}
}
int main(int argc, char *argv)
{
itk::Image<unsigned char, 3>::Pointer outputImage = nullptr;
vtkNew<vtkPolyData> mesh;
//Fill the Mesh with a Sphre inside
vtkNew<vtkSphereSource> sphereSource;
sphereSource->SetPhiResolution(21);
sphereSource->SetThetaResolution(21);
sphereSource->SetRadius(20);
sphereSource->SetCenter(0,0,0);
sphereSource->SetOutput(mesh);
sphereSource->Update();
Voxelisation voxelAlgo;
voxelAlgo.setInputData(mesh);
voxelAlgo.execute(outputImage);
if(outputImage == nullptr)
{
qDebug() << "outputImage == nullptr";
return 1;
}
//Check refcount of the output image
qDebug() << "ref count of the output image" << outputImage->GetReferenceCount(); //1
//check some value inside the output image
{
auto imageCalculatorFilter = itk::MinimumMaximumImageCalculator<itk::Image<unsigned char, 3>>::New();
imageCalculatorFilter->SetImage(outputImage);
imageCalculatorFilter->SetRegion(outputImage->GetLargestPossibleRegion());
imageCalculatorFilter->ComputeMaximum();
const unsigned char maximum = imageCalculatorFilter->GetMaximum();
qDebug() << "maximum before delete" << maximum; // value is 255
}
// now delete the stencil object
voxelAlgo.m_vtkImageStencil->Delete();
//Check refcount of the output image again
// must be the same
qDebug() << "ref count of the output image again" << outputImage->GetReferenceCount(); //1
//check some value inside the output image again
{
auto imageCalculatorFilter = itk::MinimumMaximumImageCalculator<itk::Image<unsigned char, 3>>::New();
imageCalculatorFilter->SetImage(outputImage);
imageCalculatorFilter->SetRegion(outputImage->GetLargestPossibleRegion());
qDebug() << "before recalcul the maximum after delete";
imageCalculatorFilter->ComputeMaximum(); //CRASH here
qDebug() << "Crash just before";
const unsigned char maximum = imageCalculatorFilter->GetMaximum();
qDebug() << "maximum after delete" << maximum;
}
}
Result of this code is:
refCount inside Voxelisation 1
refCount p_outputImage inside Voxelisation 2
outputImage is not null
ref count of the output image 1
maximum before delete 255
ref count of the output image again 1
before recalcul the maximum after delete