# 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;
typedef itk::ResampleImageFilter<OutputImageType, OutputImageType> ResampleFilterType;

ResampleFilterType::Pointer resampleFilter = ResampleFilterType::New();
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();

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->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->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->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)
{

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;

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]);
}

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