<div dir="ltr">Massinissa,<div><br></div><div>The fact that you are compiling for Debug, explains why the code is taking a long time to run.</div><div><br></div><div>Please give it a try at compiling for Release, that can easily make a difference of 10X or 20X</div>
<div>in the execution time. Most of the C++ template optimizations that ITK relies on, only happen</div><div>when you compile for Release.</div><div><br></div><div>Once you run it in Release mode, please share the output of the Command observer class.</div>
<div>This should help illuminate the progression of the transform and the metric values.</div><div><br></div><div> Thanks</div><div><br></div><div> Luis</div><div><br></div><div><br></div></div><div class="gmail_extra">
<br><br><div class="gmail_quote">On Fri, Dec 13, 2013 at 2:59 PM, Massinissa Bandou <span dir="ltr"><<a href="mailto:Massinissa.Bandou@usherbrooke.ca" target="_blank">Massinissa.Bandou@usherbrooke.ca</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Thank you for replying Mr Luis!!!!<br>
<br>
I'm using win7 x64 and I'm compiling on debug mode.<br>
<br>
The *source* image is from a series of dicom<br>
format: float<br>
image dimension: 120 x 120 x 137<br>
voxel dimension: 0.5 x 0.5 x 0.596756<br>
<br>
The *target* image is a raw file.<br>
format: short<br>
image dimension: 256 x 256 x 256<br>
voxel dimension: 0.308 x 0.308 x 0.308<br>
<br>
I decided to register 2 same images and use the same code as<br>
*Registration/ImageRegistration8.cxx*<br>
but I got the following error:<br>
<br>
<<a href="http://itk-users.7.n7.nabble.com/file/n32985/error.png" target="_blank">http://itk-users.7.n7.nabble.com/file/n32985/error.png</a>><br>
<br>
and my code is:<br>
<br>
int main(int argc, char *argv[])<br>
{<br>
<br>
vtkSmartPointer<vtkImageReader> imagetarget =<br>
vtkSmartPointer<vtkImageReader>::New();<br>
<br>
imagetarget->SetFileName("C:/Users/Massi/Desktop/Engineering/IRM_CT_FIDUCIALS/CT_(130802_m_FDGF18_004ct)/image/CT_volume.raw");<br>
imagetarget->SetFileDimensionality(3);<br>
imagetarget->SetDataExtent(0,511,0,511,0,511);<br>
imagetarget->SetDataSpacing(0.154,0.154,0.154);<br>
imagetarget->SetDataByteOrderToLittleEndian();<br>
imagetarget->SetDataScalarTypeToShort();<br>
imagetarget->Update();<br>
<br>
vtkSmartPointer<vtkImageReslice> resliceMyTargetImage =<br>
vtkSmartPointer<vtkImageReslice>::New();<br>
resliceMyTargetImage ->SetInputConnection(imagetarget->GetOutputPort());<br>
resliceMyTargetImage->SetOutputExtent(0,255, 0,255, 0,255);<br>
double p[3];<br>
for(unsigned int i=0;i<3;i++){<br>
p[i] = (512 * 0.154) / 256;<br>
}<br>
resliceMyTargetImage->SetOutputSpacing(p[0],p[1],p[2]);<br>
resliceMyTargetImage->SetOutputDimensionality(3);<br>
resliceMyTargetImage->SetInterpolationModeToCubic();<br>
resliceMyTargetImage->AutoCropOutputOn();<br>
resliceMyTargetImage->SetResliceAxesOrigin(0,0,0);<br>
resliceMyTargetImage->Update();<br>
<br>
vtkSmartPointer<vtkDICOMImageReader> imagesource =<br>
vtkSmartPointer<vtkDICOMImageReader>::New();<br>
<br>
imagesource->SetDirectoryName("C:/Users/Massi/Downloads/M016_Labpet8_20130425/T_6");<br>
imagesource->Update();<br>
<br>
<br>
typedef itk::Image<float, 3> FixedImageType;<br>
typedef itk::Image<short, 3> MovingImageType;<br>
<br>
/*transform vtk image to itk image*/<br>
typedef itk::VTKImageToImageFilter<FixedImageType> FixedImageReaderType;<br>
typedef itk::VTKImageToImageFilter<MovingImageType> MovingImageReaderType;<br>
<br>
FixedImageReaderType::Pointer fixedImageReader =<br>
FixedImageReaderType::New();<br>
MovingImageReaderType::Pointer movingImageReader =<br>
MovingImageReaderType::New();<br>
fixedImageReader->SetInput(imagesource->GetOutput()); // source<br>
movingImageReader->SetInput(resliceMyTargetImage->GetOutput()); //target<br>
<br>
typedef itk::VersorRigid3DTransform< double > TransformType;<br>
typedef itk::VersorRigid3DTransformOptimizer<br>
OptimizerType;<br>
typedef itk::MeanSquaresImageToImageMetric< FixedImageType, MovingImageType<br>
> MetricType;<br>
typedef itk:: LinearInterpolateImageFunction< MovingImageType, double ><br>
InterpolatorType;<br>
typedef itk::ImageRegistrationMethod< FixedImageType, MovingImageType ><br>
RegistrationType;<br>
<div class="im"><br>
MetricType::Pointer metric = MetricType::New();<br>
</div><div class="im"> OptimizerType::Pointer optimizer = OptimizerType::New();<br>
InterpolatorType::Pointer interpolator = InterpolatorType::New();<br>
RegistrationType::Pointer registration = RegistrationType::New();<br>
</div> registration->SetMetric( metric );<br>
registration->SetOptimizer( optimizer );<br>
registration->SetInterpolator( interpolator );<br>
<div class="im"><br>
TransformType::Pointer transform = TransformType::New();<br>
</div> registration->SetTransform( transform );<br>
<br>
fixedImageReader->Update();<br>
<br>
registration->SetFixedImageRegion(fixedImageReader->GetOutput()->GetBufferedRegion());<br>
<br>
typedef<br>
itk::CenteredTransformInitializer<TransformType,FixedImageType,MovingImageType><br>
TransformInitializerType;<br>
TransformInitializerType::Pointer initializer =<br>
TransformInitializerType::New();<br>
<br>
initializer->SetTransform(transform);<br>
initializer->SetFixedImage(fixedImageReader->GetOutput());<br>
initializer->SetMovingImage(movingImageReader->GetOutput());<br>
initializer->MomentsOn();<br>
initializer->InitializeTransform();<br>
<br>
typedef TransformType::VersorType VersorType;<br>
typedef VersorType::VectorType VectorType;<br>
VersorType rotation;<br>
VectorType axis;<br>
axis[0] = 0.0;<br>
axis[1] = 0.0;<br>
axis[2] = 1.0;<br>
const double angle = 0;<br>
rotation.Set( axis, angle );<br>
transform->SetRotation( rotation );<br>
<br>
registration->SetInitialTransformParameters( transform->GetParameters() );<br>
<br>
typedef OptimizerType::ScalesType OptimizerScalesType;<br>
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters() );<br>
const double translationScale = 1.0 / 1000.0;<br>
optimizerScales[0] = 1.0;<br>
optimizerScales[1] = 1.0;<br>
optimizerScales[2] = 1.0;<br>
optimizerScales[3] = translationScale;<br>
optimizerScales[4] = translationScale;<br>
optimizerScales[5] = translationScale;<br>
optimizer->SetScales( optimizerScales );<br>
optimizer->SetMaximumStepLength( 0.2000 );<br>
optimizer->SetMinimumStepLength( 0.0001 );<br>
optimizer->SetNumberOfIterations( 200 );<br>
<br>
CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();<br>
optimizer->AddObserver( itk::IterationEvent(), observer );<br>
<br>
try{<br>
<div class="im"> registration->Update();<br>
std::cout << "Optimizer stop condition: "<<<br>
registration->GetOptimizer()->GetStopConditionDescription()<< std::endl;<br>
}<br>
</div> catch( itk::ExceptionObject & err ){<br>
std::cerr << "ExceptionObject caught !" << std::endl;<br>
std::cerr << err << std::endl;<br>
}<br>
<br>
OptimizerType::ParametersType finalParameters =<br>
registration->GetLastTransformParameters();<br>
const double versorX = finalParameters[0];<br>
const double versorY = finalParameters[1];<br>
const double versorZ = finalParameters[2];<br>
const double finalTranslationX = finalParameters[3];<br>
const double finalTranslationY = finalParameters[4];<br>
const double finalTranslationZ = finalParameters[5];<br>
const unsigned int numberOfIterations = optimizer->GetCurrentIteration();<br>
const double bestValue = optimizer->GetValue();<br>
<br>
std::cout << std::endl << std::endl;<br>
std::cout << " Result = " << std::endl;<br>
std::cout << " versor X = " << versorX << std::endl;<br>
std::cout << " versor Y = " << versorY << std::endl;<br>
std::cout << " versor Z = " << versorZ << std::endl;<br>
std::cout << " Translation X = " << finalTranslationX << std::endl;<br>
std::cout << " Translation Y = " << finalTranslationY << std::endl;<br>
std::cout << " Translation Z = " << finalTranslationZ << std::endl;<br>
std::cout << " Iterations = " << numberOfIterations << std::endl;<br>
std::cout << " Metric value = " << bestValue << std::endl;<br>
<br>
transform->SetParameters( finalParameters );<br>
TransformType::MatrixType matrix = transform->GetMatrix();<br>
TransformType::OffsetType offset = transform->GetOffset();<br>
std::cout << "Matrix = " << std::endl << matrix << std::endl;<br>
std::cout << "Offset = " << std::endl << offset << std::endl;<br>
<br>
typedef itk::ResampleImageFilter< MovingImageType,FixedImageType ><br>
<div class="im">ResampleFilterType;<br>
TransformType::Pointer finalTransform = TransformType::New();<br>
</div> finalTransform->SetCenter( transform->GetCenter() );<br>
finalTransform->SetParameters( finalParameters );<br>
finalTransform->SetFixedParameters( transform->GetFixedParameters() );<br>
ResampleFilterType::Pointer resampler = ResampleFilterType::New();<br>
resampler->SetTransform( finalTransform );<br>
resampler->SetInput( movingImageReader->GetOutput() );<br>
FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();<br>
resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );<br>
resampler->SetOutputOrigin( fixedImage->GetOrigin() );<br>
resampler->SetOutputSpacing( fixedImage->GetSpacing() );<br>
resampler->SetOutputDirection( fixedImage->GetDirection() );<br>
resampler->SetDefaultPixelValue( 100 );<br>
<br>
typedef unsigned char OutputPixelType;<br>
typedef itk::Image< OutputPixelType, 3 > OutputImageType;<br>
typedef itk::CastImageFilter< FixedImageType, OutputImageType ><br>
CastFilterType;<br>
typedef itk::ImageFileWriter< OutputImageType > WriterType;<br>
//WriterType::Pointer writer = WriterType::New();<br>
CastFilterType::Pointer caster = CastFilterType::New();<br>
//writer->SetFileName( "result.raw" );<br>
caster->SetInput( resampler->GetOutput() );<br>
//writer->SetInput( caster->GetOutput() );<br>
//writer->Update();<br>
<br>
typedef<br>
itk::SubtractImageFilter<FixedImageType,FixedImageType,FixedImageType ><br>
DifferenceFilterType;<br>
DifferenceFilterType::Pointer difference = DifferenceFilterType::New();<br>
typedef itk::RescaleIntensityImageFilter<FixedImageType,OutputImageType ><br>
RescalerType;<br>
RescalerType::Pointer intensityRescaler = RescalerType::New();<br>
intensityRescaler->SetInput( difference->GetOutput() );<br>
intensityRescaler->SetOutputMinimum( 0 );<br>
intensityRescaler->SetOutputMaximum( 255 );<br>
difference->SetInput1( fixedImageReader->GetOutput() );<br>
difference->SetInput2( resampler->GetOutput() );<br>
resampler->SetDefaultPixelValue( 1 );<br>
<br>
<br>
/*connect to VTK*/<br>
typedef itk::ImageToVTKImageFilter<OutputImageType> ConnectorType;<br>
<div class="im"> ConnectorType::Pointer connector = ConnectorType::New();<br>
</div> connector->SetInput(caster->GetOutput());<br>
connector->Update();<br>
<br>
vtkSmartPointer<vtkGPUVolumeRayCastMapper> volumeMapper =<br>
vtkSmartPointer<vtkGPUVolumeRayCastMapper>::New();<br>
volumeMapper->SetInput(connector->GetOutput());<br>
volumeMapper->AutoAdjustSampleDistancesOff();<br>
volumeMapper->SetBlendModeToComposite();<br>
<br>
/********************set color & opacity*******************/<br>
vtkSmartPointer<vtkVolumeProperty> volumeProperty =<br>
vtkSmartPointer<vtkVolumeProperty>::New();<br>
vtkSmartPointer<vtkPiecewiseFunction> compositeOpacity =<br>
vtkSmartPointer<vtkPiecewiseFunction>::New();<br>
vtkSmartPointer<vtkColorTransferFunction> color =<br>
vtkSmartPointer<vtkColorTransferFunction>::New();<br>
<br>
<br>
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[0],0.0);<br>
<br>
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*3/255,0.0);<br>
<br>
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*15/255,0.1);<br>
<br>
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*70/255,0.2);<br>
<br>
compositeOpacity->AddPoint(connector->GetOutput()->GetScalarRange()[1]*75/255,0.3);<br>
<br>
<br>
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*0/255,0.0,0.0,0.0);<br>
<br>
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*3/255,1.0,0.3,0.3);//<br>
<br>
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*20/255,0.0,0.0,1.0);<br>
<br>
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*50/255,0.9,0.37,0.27);<br>
<br>
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*80.0/255,0.15,0.15,0.5);<br>
<br>
color->AddRGBPoint(connector->GetOutput()->GetScalarRange()[1]*120.0/255,1.0,1.0,1.0);<br>
<br>
volumeProperty->SetAmbient(0.5);<br>
volumeProperty->SetDiffuse(0.8);<br>
volumeProperty->SetSpecular(8.0);<br>
volumeProperty->SetInterpolationTypeToLinear();<br>
volumeProperty->SetColor(color);<br>
volumeProperty->SetScalarOpacity(compositeOpacity);<br>
volumeProperty->ShadeOn();<br>
<br>
vtkSmartPointer<vtkVolume> volume = vtkSmartPointer<vtkVolume>::New();<br>
volume->SetMapper(volumeMapper);<br>
volume->SetProperty(volumeProperty);<br>
<br>
vtkSmartPointer<vtkRenderWindow> renderWindow =<br>
vtkSmartPointer<vtkRenderWindow>::New();<br>
vtkSmartPointer<vtkRenderer> renderer =<br>
vtkSmartPointer<vtkRenderer>::New();<br>
renderWindow->AddRenderer(renderer);<br>
renderer->AddVolume(volume);<br>
renderer->SetBackground(0.9,0.9,0.9);<br>
renderer->GradientBackgroundOn();<br>
vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =<br>
vtkSmartPointer<vtkRenderWindowInteractor>::New();<br>
renderWindowInteractor->SetRenderWindow(renderWindow);<br>
<br>
/*Set camera position*/<br>
renderer->GetActiveCamera()->Azimuth(30);<br>
renderer->GetActiveCamera()->Elevation(30);<br>
renderer->ResetCamera();<br>
renderWindow->Render();<br>
renderWindowInteractor->Start();<br>
<br>
return EXIT_SUCCESS;<br>
}<br>
<br>
<br>
<br>
<br>
--<br>
View this message in context: <a href="http://itk-users.7.n7.nabble.com/registration-using-Mutual-Information-tp32972p32985.html" target="_blank">http://itk-users.7.n7.nabble.com/registration-using-Mutual-Information-tp32972p32985.html</a><br>
<div class="im HOEnZb">Sent from the ITK - Users mailing list archive at Nabble.com.<br>
_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Kitware offers ITK Training Courses, for more information visit:<br>
<a href="http://www.kitware.com/products/protraining.php" target="_blank">http://www.kitware.com/products/protraining.php</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at:<br>
<a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
</div><div class="HOEnZb"><div class="h5">_______________________________________________<br>
Community mailing list<br>
<a href="mailto:Community@itk.org">Community@itk.org</a><br>
<a href="http://public.kitware.com/cgi-bin/mailman/listinfo/community" target="_blank">http://public.kitware.com/cgi-bin/mailman/listinfo/community</a><br>
</div></div></blockquote></div><br></div>