correction of gantry tilted ct image

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:

const unsigned int dim= 3;
typedef itk::AffineTransform< float, dim> TransformType;

TransformType::Pointer transform = TransformType::New(); 
transform->Shear(int axis1, int axis2, TParametersValueType coef); // ?

itk::AffineTransform<float, dim>::OutputVectorType translation;
translation[0] = ;
translation[1] = ;
translation[2] = ;
transform->Translate(translation); //?

How can I determine the direction axis and shear factor for the Shear() function and in which direction should I translate the dataset?
Thank you!

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?

if (!m_ForceOrthogonalDirection)
{
  for (j = 0; j < TOutputImage::ImageDimension; ++j)
  {
    direction[j][this->m_NumberOfDimensionsInImage] = dirN[j] / dirNnorm;
  }
}
output->SetOrigin(origin);
output->SetSpacing(spacing);
output->SetDirection(direction);
output->SetLargestPossibleRegion(largestRegion);

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();

Do you get a non-black image if you replace exporter->SetInput(itkIm); by exporter->SetInput(seriesReader->GetOutput());?

What if you replace it by exporter->SetInput(resampleFilter->GetOutput());?

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:

ItkImageType::Pointer itkIm = seriesReader->GetOutput();
itkIm ->SetDirection(direction);
itkIm ->SetSpacing(spacing);
itkIm ->SetOrigin(origin);

ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
resampleFilter->SetInput(itkIm);
resampleFilter->SetOutputOrigin(origin);
resampleFilter->SetOutputSpacing(spacing);
resampleFilter->SetOutputDirection(direction);
resampleFilter->SetSize(outputSize);
resampleFilter->Update();

//itkIm = resampleFilter->GetOutput();
//itkIm->Update();

itk::VTKImageExport<ItkImageType>::Pointer exporter = itk::VTKImageExport<ItkImageType>::New();
exporter->SetInput(resampleFilter->GetOutput());
exporter->ReleaseDataFlagOn();
exporter->Update();

If I set following pipeline, I get only black backgrounds

ItkImageType::Pointer itkIm = seriesReader->GetOutput();
itkIm ->SetDirection(direction);
itkIm ->SetSpacing(spacing);
itkIm ->SetOrigin(origin);

ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
resampleFilter->SetInput(seriesReader->GetOutput());
resampleFilter->SetOutputOrigin(origin);
resampleFilter->SetOutputSpacing(spacing);
resampleFilter->SetOutputDirection(direction);
resampleFilter->SetSize(outputSize);
resampleFilter->Update();

itkIm = resampleFilter->GetOutput();
itkIm->Update();

itk::VTKImageExport<ItkImageType>::Pointer exporter = itk::VTKImageExport<ItkImageType>::New();
exporter->SetInput(resampleFilter->GetOputput());
exporter->ReleaseDataFlagOn();
exporter->Update();

If I set following pipeline, I can not read the image dataset:

ItkImageType::Pointer itkIm = seriesReader->GetOutput();
itkIm ->SetDirection(direction);
itkIm ->SetSpacing(spacing);
itkIm ->SetOrigin(origin);

ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
resampleFilter->SetInput(seriesReader->GetOutput());
resampleFilter->SetOutputOrigin(origin);
resampleFilter->SetOutputSpacing(spacing);
resampleFilter->SetOutputDirection(direction);
resampleFilter->SetSize(outputSize);
resampleFilter->Update();

itkIm = resampleFilter->GetOutput();
itkIm->Update();

itk::VTKImageExport<ItkImageType>::Pointer exporter = itk::VTKImageExport<ItkImageType>::New();
exporter->SetInput(itkIm);
exporter->ReleaseDataFlagOn();
exporter->Update();

EDIT
What if I use the itkIm, which contains new direction as a reference image in resampleFilter and give still the seriesReader as input?

ItkImageType::Pointer itkIm = seriesReader->GetOutput();
itkIm ->SetDirection(direction);
itkIm ->SetSpacing(spacing);
itkIm ->SetOrigin(origin);

ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
resampleFilter->SetInput(seriesReader->GetOutput());
resampleFilter->UseReferenceImageOn();
resampleFilter->SetReferenceImage(itkIm);
resampleFilter->SetOutputOrigin(origin);
resampleFilter->SetOutputSpacing(spacing);
resampleFilter->SetOutputDirection(direction);
resampleFilter->SetSize(outputSize);
resampleFilter->Update();

itk::VTKImageExport<ItkImageType>::Pointer exporter = itk::VTKImageExport<ItkImageType>::New();
exporter->SetInput(resampleFilter->GetOputput());
exporter->ReleaseDataFlagOn();
exporter->Update();

So I can read and visualize images and the image direction in resamleFilter is still correct:

DIRECTION in resampleFilter : 
1 0 0
0 0.984808 0.173648
0 -0.173648 0.984808

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.

itk::WriteImage does not exist in version 3.20 therefore I had used itk::ImageFileWriter

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.

Does not matter if I write the nrrd image before or after tilt correction, the images are represented in slicer correctly. How it is possible?

Maybe the code was modified to force orthogonal orientation with ITKv4, so 3.20 doesn’t have it?

The version 3.20 does also the same mathematical operations as version 4 to correct the tilt in series reader :slight_smile: 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.