Hello,
I need help by understanding the correct pipleline and determining the share factor of a gantry tilted ct dataset, so that I can apply a transform on it.
I know the exact angle of the gantry over read out DICOM Tag (-15°) and that I can use the Shear() function from itkAffineTransform class for the correction:
If you call reader->ForceOrthogonalDirectionOff();, the image direction should already be correct and you don’t need any manual processing. See the complete example.
Thank you for your response. Does it mean that this function corrects the tilt already during reading of images and visualize them series correctly (square)?
ForceOrthogonalDirectionOff() is implemented in version 5.2.1 and if I understand correctly it calculates the position difference of the first and last image in the dataset and shifts the images properly in the correct direction?
Yes. That was added specifically to correctly read gantry-tilted CT images. If other steps in the processing pipeline have trouble handling non-orthogonal images, you can resample the images into ortho-normal grid immediately after reading.
Since I am using an older version of the itk, I have implemented the ForceOrthogonalDirectionOff() function in my source code same as in itkImagesSeriesReader.hxx. As expected the image direction matrix is shifted:
Image direction before applying ForceOrthogonalDirectionOff on the image:
1 0 0
0 1 0
0 0 1
Image direction after applying ForceOrthogonalDirectionOff:
1 0 0
0 0.984808 0.173648
0 -0.173648 0.984808
The values in matrix -0.173648 and 0.173648 is exactly the sin of the tilt angle in radiance, if I calculate it separately: float sinOfTiltAngle = sin(vtkMath::RadiansFromDegrees(-10.0));
The image origin and spacing are remaining unchanged after applying the function. I guess it is ok? My goal is, doing the tilt correction on the fly without needing to save and load the ct dataset again. I will try this resampling (Resample), if you have no other recommendation.
Thank you!
When I was writing this, I had plenty of data to work with. Some of the series had 2D X-rays which allowed me to check whether the data looked correct for rendering and resampling into ortho-normal images.
I don’t remember whether the spacing changes, but it is computed by this line spacing[this->m_NumberOfDimensionsInImage] = dirNnorm / ( numberOfFiles - 1 );.
After resampling the image dataset as following and connecting it to the visualization pipeline, I get only a black the screen, which does not contains any oriented image (axial, sagittal,cornal). What I am doing wrong during resampling?
typedef short PixelType;
const unsigned int Dimension = 3;
typedef itk::OrientedImage<PixelType, Dimension> OutputImageType;
const OutputImageType* first = firstReader->GetOutput();
typedef itk::ResampleImageFilter<OutputImageType, OutputImageType> ResampleFilterType;
ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
resampleFilter->SetInput(seriesReader->GetOutput());
resampleFilter->SetOutputOrigin(imOrigin); // does not change
resampleFilter->SetOutputSpacing(imSpacing); // does not chaneg
resampleFilter->SetOutputDirection(imDirection); // shifted in correct direction
resampleFilter->SetSize(imOutputSize); // same as series reader size (512x512xnumberOfImages)
resampleFilter->Update();
ItkImageType::Pointer itkIm = seriesReader->GetOutput();
itkIm = resampleFilter->GetOutput();
itkIm ->Update();
itk::VTKImageExport<ItkImageType>::Pointer exporter = itk::VTKImageExport<ItkImageType>::New();
exporter->SetInput(itkIm);
exporter->ReleaseDataFlagOn();
exporter->Update();
If I use following pipeline as you mentioned, I get colored backgrounds without any image on them. I guess what I represented is only renderer from vtk:
Does resample filter work correctly? Try writing it via itk::WriteImage(resampleFilter->Update(), "temp.nrrd"); and visualize it using Slicer or ITK-SNAP.
If it does, this is more of a question for VTK than ITK.
I don’t think that resample filter works as expected. The written nrrd data is visualized as black (without any image) in slicer. Please see my whole code block:
if (gantry != 0.0)
{
ItkImageType::Pointer itkIm = seriesReader->GetOutput();
bool reverseOrder = true;
const auto numberOfFiles = static_cast<int>(m_Filenames.size());
const int firstFileName = (reverseOrder ? numberOfFiles - 1 : 0);
const int lastFileName = (reverseOrder ? 0 : numberOfFiles - 1);
typedef itk::OrientedImage<short, 3> ItkImageType;
using ReaderType = itk::ImageFileReader<ItkImageType>;
ReaderType::Pointer firstReader = ReaderType::New();
ReaderType::Pointer lastReader = ReaderType::New();
firstReader->SetFileName(m_Filenames[firstFileName].toStdString().c_str());
firstReader->Update();
firstReader->UpdateOutputInformation();
lastReader->SetFileName(m_Filenames[lastFileName].toStdString().c_str());
lastReader->Update();
lastReader->UpdateOutputInformation();
const ItkImageType *firstImage = firstReader->GetOutput();
typename ItkImageType::SpacingType spacing = firstImage->GetSpacing();
typename ItkImageType::PointType origin = firstImage->GetOrigin();
typename ItkImageType::DirectionType direction = firstImage->GetDirection();
typename ItkImageType::RegionType::ImageRegion largestRegion = firstImage->GetLargestPossibleRegion();
// Set size of the processed dataset same as the original dataset
ItkImageType::SizeType outputSize;
typedef ItkImageType::SizeType::SizeValueType SizeValueType;
const ItkImageType::SizeType &inputSize = largestRegion.GetSize();
outputSize[0] = inputSize[0];
outputSize[1] = inputSize[1];
outputSize[2] = numberOfFiles;
using SpacingScalarType = typename ItkImageType::SpacingValueType;
itk::Array<SpacingScalarType> position1(ItkImageType::ImageDimension);
position1.Fill(0.0f);
itk::Array<SpacingScalarType> positionN(ItkImageType::ImageDimension);
positionN.Fill(0.0f);
// Initialize the position to the origin returned by the reader
unsigned int j;
for (j = 0; j < ItkImageType::ImageDimension; j++)
{
position1[j] = static_cast<SpacingScalarType>(origin[j]);
}
const ItkImageType* last = lastReader->GetOutput();
// Initialize the position to the origin returned by the reader
for (j = 0; j < ItkImageType::ImageDimension; j++)
{
positionN[j] = static_cast<SpacingScalarType>(last->GetOrigin()[j]);
}
// Compute and set the inter slice spacing and last (usually third) axis of direction
itk::Vector<SpacingScalarType, ItkImageType::ImageDimension> dirN;
for (j = 0; j < ItkImageType::ImageDimension; ++j)
{
dirN[j] = positionN[j] - position1[j];
}
SpacingScalarType dirNnorm = dirN.GetNorm();
unsigned int numberOfDimensionsInImage = firstReader->GetImageIO()->GetNumberOfDimensions();
spacing[numberOfDimensionsInImage] = dirNnorm / (numberOfFiles - 1);
// Shift the image direction and set it to the output image
for (j = 0; j < ItkImageType::ImageDimension; ++j)
{
direction[j][numberOfDimensionsInImage] = dirN[j] / dirNnorm;
}
itkIm->SetDirection(direction); //Corrected image direction matrix
itkIm->SetSpacing(spacing);
resampleFilter = ResampleFilterType::New();
resampleFilter->SetInput(itkIm);
//resampleFilter->UseReferenceImageOn();
//resampleFilter->SetReferenceImage(itkIm);
resampleFilter->SetOutputOrigin(itkIm->GetOrigin());
resampleFilter->SetOutputSpacing(itkIm->GetSpacing());
resampleFilter->SetOutputDirection(itkIm->GetDirection());
resampleFilter->SetSize(outputSize);
resampleFilter->Update();
//itk::WriteImage(resampleFilter->Update(), "temp.nrrd");
using WriterType = itk::ImageFileWriter<ItkImageType>;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName("temp.nrrd");
writer->SetInput(resampleFilter->GetOutput());
writer->Update();
//itkIm = resampleFilter->GetOutput();
//itkIm->Update();
itk::VTKImageExport<ItkImageType>::Pointer exporter = itk::VTKImageExport<ItkImageType>::New();
exporter->SetInput(resampleFilter->GetOutput());
exporter->ReleaseDataFlagOn();
exporter->Update();
}
How about an even simpler test then: itk::WriteImage(itkIm, "temp.mha");. Then visualize this image in slicer. It should be correctly visualized with shear (gantry tilt). Current preview release handles such images correctly, and I think that was true for 5 years if not more.
Sure, you can use ImageFileWriter. Can you write the image just after reading the DICOM series, before you try to resample it? I suspect that your backporting of this feature to version 3.20 is not quite correct.
The version 3.20 does also the same mathematical operations as version 4 to correct the tilt in series reader Therefore the tilt is corrected automatically after reading it, without applying any additional operations. It means I had to be more careful during checking the functionality. Thank you for your support.
Therefore the tilt is corrected automatically after reading it, without applying any additional operations. It means I had to be more careful during checking the functionality. Thank you for your support.