<p>Hi,luis,<br> I have several questions about ITK deformable registration.<br> First, I use the Demons for 3D deformable registration, when I set the parameter "SetStandardDeviations" smaller(for example 1.0), <br>
the image whose size is 256*176*16 will have black regions in each slice, but when I set this paremeter bigger (for example 10.0),<br> only the first and the last slice have black regions, so I do not know why this phenomena will happen? <br>
second, since the phenomena happened, so I use the "MirrorPadImageFilter" to pad the image in the slice direction in order to <br> make its size 256*176*18, and set the "SetStandardDeviations" parameter 10.0, then after registration, <br>
I use the WarpImageFilter to warp the movingImage according to the deformable field, at last I use the "RegionOfInterestImageFilter" extract the image from second slice to the seventeenth slice. Thus I can get better result, but I don't know this method is whether correct.<br>
In order to explain these problem, I copy the part programs in the email, and the result images in the inclosure.<br> Thank you for your answers as soon as possible.<br> <br> With best regards!<br> <br> HuaJiang,<br>
</p>
<p>The Program:(this is the part of the my program, the read and write parts arenot attached)<br>FixedImageCasterType::Pointer fixedImageCaster = FixedImageCasterType::New();<br> MovingImageCasterType::Pointer movingImageCaster = MovingImageCasterType::New();<br>
fixedImageCaster->SetInput( fixedImage );<br> movingImageCaster->SetInput( movingImage ); <br>/*<br> HistogramMatchFilterType::Pointer matcher = HistogramMatchFilterType::New();<br> matcher->SetInput( movingImageCaster->GetOutput() );<br>
matcher->SetReferenceImage( fixedImageCaster->GetOutput() );<br> <br> matcher->SetNumberOfHistogramLevels( 1024 );<br> matcher->SetNumberOfMatchPoints( 7 );<br> matcher->ThresholdAtMeanIntensityOn();<br>
*/<br> <br> <br> Demons3RegistrationFilterType::Pointer filter = Demons3RegistrationFilterType::New();<br> <br> filter->SetStandardDeviations( 10.0 );<br> MultiResRegistrationFilterType:: Pointer multires = MultiResRegistrationFilterType::New();<br>
<br> multires->SetRegistrationFilter( filter );<br> multires->SetNumberOfLevels( 3 );<br> <br> unsigned int nIterations[4] = {32, 16, 4 };<br> multires->SetNumberOfIterations( nIterations );<br> DemonsCommandIterationUpdate::Pointer observer = DemonsCommandIterationUpdate::New();<br>
filter->AddObserver(itk::IterationEvent(), observer);<br> <br> CommandResolutionLevelUpdate::Pointer levelobserver = CommandResolutionLevelUpdate::New();<br> multires->AddObserver( itk::IterationEvent(), levelobserver );<br>
<br> // Pad the image.<br> unsigned long upFactor[3] = {0,0,1};<br> unsigned long lowFactor[3] = {0,0,1};<br> InternalPixelType defaultPixel = 0;<br> // define warp filter.<br> WarperType::Pointer warper = WarperType::New();<br>
typedef itk::LinearInterpolateImageFunction< InternalImageType, double > LinearInterpolatorType;<br> LinearInterpolatorType::Pointer linearInterpolator = LinearInterpolatorType::New();<br> <br> warper->SetInterpolator( linearInterpolator );<br>
ConstantPadImageFilterType::Pointer constantPadFixedFilter; <br> ConstantPadImageFilterType::Pointer constantPadMovingFilter; <br> MirrorPadImageFilterType::Pointer mirrorPadFixedFilter;<br> MirrorPadImageFilterType::Pointer mirrorPadMovingFilter; <br>
if (m_padMode == ConstantPad)<br> {<br> constantPadFixedFilter = ConstantPadImageFilterType::New();<br> constantPadMovingFilter = ConstantPadImageFilterType::New();<br> <br> constantPadFixedFilter->SetInput(fixedImageCaster->GetOutput());<br>
constantPadMovingFilter->SetInput(movingImageCaster->GetOutput());<br> <br> constantPadFixedFilter->SetPadLowerBound(lowFactor);<br> constantPadFixedFilter->SetPadUpperBound(upFactor);<br> constantPadFixedFilter->SetConstant(defaultPixel);<br>
constantPadMovingFilter->SetPadLowerBound(lowFactor);<br> constantPadMovingFilter->SetPadUpperBound(upFactor);<br> constantPadMovingFilter->SetConstant(defaultPixel);<br> multires->SetFixedImage( constantPadFixedFilter->GetOutput() );<br>
multires->SetMovingImage( constantPadMovingFilter->GetOutput() );<br> <br> try<br> {<br> multires->Update();<br> }<br> catch( itk::ExceptionObject & excp )<br> {<br> std::cerr << excp << std::endl;<br>
return ;<br> }<br> warper->SetInput(constantPadMovingFilter->GetOutput());<br> warper->SetOutputSpacing( constantPadFixedFilter->GetOutput()->GetSpacing() );<br> warper->SetOutputOrigin( constantPadFixedFilter->GetOutput()->GetOrigin() );<br>
<br> warper->SetDeformationField( multires->GetOutput() );<br> warper->Update();<br> }<br> else if (m_padMode == MirrorPad)<br> {<br> mirrorPadFixedFilter = MirrorPadImageFilterType::New();<br> mirrorPadMovingFilter = MirrorPadImageFilterType::New();</p>
<p> mirrorPadFixedFilter->SetInput(fixedImageCaster->GetOutput());<br> mirrorPadMovingFilter->SetInput(movingImageCaster->GetOutput());<br> <br> mirrorPadFixedFilter->SetPadLowerBound(lowFactor);<br> mirrorPadFixedFilter->SetPadUpperBound(upFactor);<br>
<br> mirrorPadMovingFilter->SetPadLowerBound(lowFactor);<br> mirrorPadMovingFilter->SetPadUpperBound(upFactor);<br> <br> multires->SetFixedImage( mirrorPadFixedFilter->GetOutput() );<br> multires->SetMovingImage( mirrorPadMovingFilter->GetOutput() );<br>
<br> try<br> {<br> multires->Update();<br> }<br> catch( itk::ExceptionObject & excp )<br> {<br> std::cerr << excp << std::endl;<br> return ;<br> }<br> <br> warper->SetInput(mirrorPadMovingFilter->GetOutput());<br>
warper->SetOutputSpacing( mirrorPadFixedFilter->GetOutput()->GetSpacing() );<br> warper->SetOutputOrigin( mirrorPadFixedFilter->GetOutput()->GetOrigin() );<br> <br> warper->SetDeformationField( multires->GetOutput() );<br>
warper->Update();<br> }<br> <br> // extract the image.<br> RegionOfInterestImageFilterType::Pointer regionExtractFitler<br> = RegionOfInterestImageFilterType::New();<br> <br> regionExtractFitler->SetInput(warper->GetOutput());<br>
InternalImageType::RegionType desiredRegion ;<br> InternalImageType::SizeType size = fixedImage->GetLargestPossibleRegion().GetSize();<br> InternalImageType::IndexType index;<br> index[0] = 0;<br> index[1] = 0;<br> index[2] = 0;<br>
desiredRegion.SetSize(size);<br> desiredRegion.SetIndex(index);<br> regionExtractFitler->SetRegionOfInterest(desiredRegion);<br> InternalImageCasterType::Pointer internalCaster = InternalImageCasterType::New();<br>
internalCaster->SetInput(regionExtractFitler->GetOutput());<br> internalCaster->Update();</p>