Hello All,
I have some 3D CT brain images in nii.gz format. I would like to use one image as fixed image and another as moving image.
I also have some blood vessels in vtk format in Image Coordinate System and in World coordinates System.
This is my input 1.fixed image 2. moving image 3. blood vessels in moving image(Image or world)
I want to get 1.the output image 2.corresponding blood vessel in the output image in Image Coordinate System.
To put it simply, I want to get the corresponding coordinates in Image Coordinate System in the output image through the input coordinates in moving image.
I modified ImageRegistration20.cxx to get this output, the code can run, but there is a problem with the output vessel coordinates.
I tried finaltransform and inverseFinalTransform to transform the coordinates.
I also tried to use Image coordinates as input to directly transform and use World coordinates as input to transform and then convert to Image coordinates.
But I still can’t get the correct corresponding coordinates in the Image Coordinate System. So i want to find the right way.
Any help, guidance or tutorials will greatly be appreciated. Thank you.
here’s the code:
/*=========================================================================
*
* Copyright Insight Software Consortium
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
// Software Guide : BeginCommandLineArgs
// INPUTS: {brainweb1e1a10f20.mha}
// INPUTS: {brainweb1e1a10f20Rot10Tx15.mha}
// ARGUMENTS: ImageRegistration20Output.mhd
// Software Guide : EndCommandLineArgs
// Software Guide : BeginLatex
//
// This example illustrates the use of the \doxygen{AffineTransform}
// for performing registration in $3D$.
//
// \index{itk::AffineTransform}
//
// Software Guide : EndLatex
#include "itkImageRegistrationMethod.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkCenteredTransformInitializer.h"
// Software Guide : BeginLatex
//
// Let's start by including the header file of the AffineTransform.
//
// \index{itk::AffineTransform!header}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
#include "itkAffineTransform.h"
// Software Guide : EndCodeSnippet
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkResampleImageFilter.h"
#include "itkCastImageFilter.h"
#include "itkSubtractImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
//
// The following piece of code implements an observer
// that will monitor the evolution of the registration process.
//
#include "itkCommand.h"
#include "vtkSmartPointer.h"
#include "vtkUnstructuredGrid.h"
#include "vtkUnstructuredGridReader.h"
#include "vtkUnstructuredGridWriter.h"
#include "vtkPoints.h"
#include "vtkPolyLine.h"
#include "itkPointSet.h"
#include <iostream>
#include <fstream>
class CommandIterationUpdate : public itk::Command
{
public:
typedef CommandIterationUpdate Self;
typedef itk::Command Superclass;
typedef itk::SmartPointer<Self> Pointer;
itkNewMacro( Self );
protected:
CommandIterationUpdate() {};
public:
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef const OptimizerType * OptimizerPointer;
void Execute(itk::Object *caller, const itk::EventObject & event) ITK_OVERRIDE
{
Execute( (const itk::Object *)caller, event);
}
void Execute(const itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE
{
OptimizerPointer optimizer = static_cast< OptimizerPointer >( object );
if( ! itk::IterationEvent().CheckEvent( &event ) )
{
return;
}
std::cout << optimizer->GetCurrentIteration() << " ";
std::cout << optimizer->GetValue() << " ";
std::cout << optimizer->GetCurrentPosition() << std::endl;
}
};
int main( int argc, char *argv[] )
{
if( argc < 4 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " fixedImageFile movingImageFile " << std::endl;
std::cerr << " outputImagefile [differenceBeforeRegistration] " << std::endl;
std::cerr << " [differenceAfterRegistration] " << std::endl;
std::cerr << " [physicalinputpoints] [physicaloutputpoints] " << std::endl;
return EXIT_FAILURE;
}
// Software Guide : BeginLatex
//
// We define then the types of the images to be registered.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
const unsigned int Dimension = 3;
typedef unsigned short PixelType;
typedef itk::Image< PixelType, Dimension > FixedImageType;
typedef itk::Image< PixelType, Dimension > MovingImageType;
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// The transform type is instantiated using the code below. The template
// parameters of this class are the representation type of the space
// coordinates and the space dimension.
//
// \index{itk::AffineTransform!Instantiation}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::AffineTransform<
double,
Dimension > TransformType;
// Software Guide : EndCodeSnippet
typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
typedef itk::MeanSquaresImageToImageMetric<
FixedImageType,
MovingImageType > MetricType;
typedef itk:: LinearInterpolateImageFunction<
MovingImageType,
double > InterpolatorType;
typedef itk::ImageRegistrationMethod<
FixedImageType,
MovingImageType > RegistrationType;
MetricType::Pointer metric = MetricType::New();
OptimizerType::Pointer optimizer = OptimizerType::New();
InterpolatorType::Pointer interpolator = InterpolatorType::New();
RegistrationType::Pointer registration = RegistrationType::New();
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
// Software Guide : BeginLatex
//
// The transform object is constructed below and passed to the registration
// method.
//
// \index{itk::AffineTransform!New()}
// \index{itk::AffineTransform!Pointer}
// \index{itk::RegistrationMethod!SetTransform()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
TransformType::Pointer transform = TransformType::New();
registration->SetTransform( transform );
// Software Guide : EndCodeSnippet
typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
fixedImageReader->SetFileName( argv[1] );
movingImageReader->SetFileName( argv[2] );
registration->SetFixedImage( fixedImageReader->GetOutput() );
registration->SetMovingImage( movingImageReader->GetOutput() );
fixedImageReader->Update();
registration->SetFixedImageRegion(
fixedImageReader->GetOutput()->GetBufferedRegion() );
// Software Guide : BeginLatex
//
// In this example, we again use the
// \doxygen{CenteredTransformInitializer} helper class in order to compute
// a reasonable value for the initial center of rotation and the
// translation. The initializer is set to use the center of mass of each
// image as the initial correspondence correction.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
typedef itk::CenteredTransformInitializer<
TransformType, FixedImageType,
MovingImageType > TransformInitializerType;
TransformInitializerType::Pointer initializer
= TransformInitializerType::New();
initializer->SetTransform( transform );
initializer->SetFixedImage( fixedImageReader->GetOutput() );
initializer->SetMovingImage( movingImageReader->GetOutput() );
initializer->MomentsOn();
initializer->InitializeTransform();
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Now we pass the parameters of the current transform as the initial
// parameters to be used when the registration process starts.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
registration->SetInitialTransformParameters(
transform->GetParameters() );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Keeping in mind that the scale of units in scaling, rotation and
// translation are quite different, we take advantage of the scaling
// functionality provided by the optimizers. We know that the first $N
// \times N$ elements of the parameters array correspond to the rotation
// matrix factor, and the last $N$ are the components of the translation to
// be applied after multiplication with the matrix is performed.
//
// Software Guide : EndLatex
double translationScale = 1.0 / 1000.0;
/* if( argc > 8 )
{
translationScale = atof( argv[8] );
} */
// Software Guide : BeginCodeSnippet
typedef OptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = 1.0;
optimizerScales[4] = 1.0;
optimizerScales[5] = 1.0;
optimizerScales[6] = 1.0;
optimizerScales[7] = 1.0;
optimizerScales[8] = 1.0;
optimizerScales[9] = translationScale;
optimizerScales[10] = translationScale;
optimizerScales[11] = translationScale;
optimizer->SetScales( optimizerScales );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We also set the usual parameters of the optimization method. In this
// case we are using an
// \doxygen{RegularStepGradientDescentOptimizer}. Below, we define the
// optimization parameters like initial step length, minimal step length
// and number of iterations. These last two act as stopping criteria for
// the optimization.
//
// Software Guide : EndLatex
double steplength = 0.1;
/* if( argc > 6 )
{
steplength = atof( argv[6] );
} */
unsigned int maxNumberOfIterations = 30;
/* if( argc > 7 )
{
maxNumberOfIterations = atoi( argv[7] );
} */
// Software Guide : BeginCodeSnippet
optimizer->SetMaximumStepLength( steplength );
optimizer->SetMinimumStepLength( 0.0001 );
optimizer->SetNumberOfIterations( maxNumberOfIterations );
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// We also set the optimizer to do minimization by calling the
// \code{MinimizeOn()} method.
//
// \index{itk::Regular\-Step\-Gradient\-Descent\-Optimizer!MinimizeOn()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
optimizer->MinimizeOn();
// Software Guide : EndCodeSnippet
// Create the Command observer and register it with the optimizer.
//
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
optimizer->AddObserver( itk::IterationEvent(), observer );
// Software Guide : BeginLatex
//
// Finally we trigger the execution of the registration method by calling
// the \code{Update()} method. The call is placed in a \code{try/catch}
// block in case any exceptions are thrown.
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
try
{
registration->Update();
std::cout << "Optimizer stop condition: "
<< registration->GetOptimizer()->GetStopConditionDescription()
<< std::endl;
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
return EXIT_FAILURE;
}
// Software Guide : EndCodeSnippet
// Software Guide : BeginLatex
//
// Once the optimization converges, we recover the parameters from the
// registration method. This is done with the
// \code{GetLastTransformParameters()} method. We can also recover the
// final value of the metric with the \code{GetValue()} method and the
// final number of iterations with the \code{GetCurrentIteration()}
// method.
//
// \index{itk::RegistrationMethod!GetValue()}
// \index{itk::RegistrationMethod!GetCurrentIteration()}
// \index{itk::RegistrationMethod!GetLastTransformParameters()}
//
// Software Guide : EndLatex
// Software Guide : BeginCodeSnippet
OptimizerType::ParametersType finalParameters =
registration->GetLastTransformParameters();
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();
const double bestValue = optimizer->GetValue();
// Software Guide : EndCodeSnippet
// Print out results
//
std::cout << "Result = " << std::endl;
std::cout << " Iterations = " << numberOfIterations << std::endl;
std::cout << " Metric value = " << bestValue << std::endl;
// The following code is used to dump output images to files.
// They illustrate the final results of the registration.
// We will resample the moving image and write out the difference image
// before and after registration. We will also rescale the intensities of the
// difference images, so that they look better!
typedef itk::ResampleImageFilter<
MovingImageType,
FixedImageType > ResampleFilterType;
TransformType::Pointer finalTransform = TransformType::New();
finalTransform->SetParameters( finalParameters );
finalTransform->SetFixedParameters( transform->GetFixedParameters() );
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetTransform( finalTransform );
resampler->SetInput( movingImageReader->GetOutput() );
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resampler->SetOutputOrigin( fixedImage->GetOrigin() );
resampler->SetOutputSpacing( fixedImage->GetSpacing() );
resampler->SetOutputDirection( fixedImage->GetDirection() );
resampler->SetDefaultPixelValue( 1 );
typedef signed short OutputPixelType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
typedef itk::CastImageFilter<
FixedImageType,
OutputImageType > CastFilterType;
typedef itk::ImageFileWriter< OutputImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
CastFilterType::Pointer caster = CastFilterType::New();
writer->SetFileName( argv[3] );
caster->SetInput( resampler->GetOutput() );
writer->SetInput( caster->GetOutput() );
writer->Update();
typedef itk::SubtractImageFilter<
FixedImageType,
FixedImageType,
FixedImageType > DifferenceFilterType;
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
difference->SetInput1( fixedImageReader->GetOutput() );
difference->SetInput2( resampler->GetOutput() );
WriterType::Pointer writer2 = WriterType::New();
typedef itk::RescaleIntensityImageFilter<
FixedImageType,
OutputImageType > RescalerType;
RescalerType::Pointer intensityRescaler = RescalerType::New();
intensityRescaler->SetInput( difference->GetOutput() );
intensityRescaler->SetOutputMinimum( 0 );
intensityRescaler->SetOutputMaximum( 255 );
writer2->SetInput( intensityRescaler->GetOutput() );
resampler->SetDefaultPixelValue( 1 );
// Compute the difference image between the
// fixed and resampled moving image.
if( argc > 5 )
{
writer2->SetFileName( argv[5] );
writer2->Update();
}
typedef itk::IdentityTransform< double, Dimension > IdentityTransformType;
IdentityTransformType::Pointer identity = IdentityTransformType::New();
// Compute the difference image between the
// fixed and moving image before registration.
if( argc > 4 )
{
resampler->SetTransform( identity );
writer2->SetFileName( argv[4] );
writer2->Update();
}
TransformType::Pointer inverseFinalTransform = TransformType::New();
finalTransform->GetInverse(inverseFinalTransform);
if (argc > 7)
{
typedef itk::Vector<float, Dimension> PixelType;
typedef itk::PointSet<PixelType, Dimension> PointSetType;
typedef PointSetType::PointType PointType;
PointType OutputPoint;
fstream DetectFile(argv[6], ios::in);
if (!DetectFile)
{
std::cout << "Cannot find the VTK file!" << endl;
return -1;
}
DetectFile.close();
vtkSmartPointer<vtkPoints> iPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkPolyLine> iLine = vtkSmartPointer<vtkPolyLine>::New();
vtkSmartPointer<vtkUnstructuredGridReader> iVtkReader = vtkSmartPointer<vtkUnstructuredGridReader>::New();
vtkSmartPointer<vtkUnstructuredGridWriter> iVtkWriter = vtkSmartPointer<vtkUnstructuredGridWriter>::New();
vtkSmartPointer<vtkUnstructuredGrid> iGridRead = vtkSmartPointer<vtkUnstructuredGrid>::New();
vtkSmartPointer<vtkUnstructuredGrid> iGridWrite = vtkSmartPointer<vtkUnstructuredGrid>::New();
iVtkReader->SetFileName(argv[6]);
iVtkReader->Update();
iGridRead = iVtkReader->GetOutput();
int nPointNum = iGridRead->GetMaxCellSize();
double dCord[3];
for (int i = 0; i < nPointNum; i++)
{
iGridRead->GetPoint(i, dCord);
OutputPoint = inverseFinalTransform->TransformPoint(dCord);
iPoints->InsertNextPoint(OutputPoint[0], OutputPoint[1], OutputPoint[2]);
}
iLine->GetPointIds()->SetNumberOfIds(nPointNum);
for (int i = 0; i < nPointNum; i++)
{
iLine->GetPointIds()->SetId(i, i);
}
iGridWrite->Allocate(1, 1);
iGridWrite->InsertNextCell(iLine->GetCellType(), iLine->GetPointIds());
iGridWrite->SetPoints(iPoints);
iVtkWriter->SetInputData(iGridWrite);
iVtkWriter->SetFileName(argv[7]);
iVtkWriter->Write();
}
return EXIT_SUCCESS;
}