[Insight-users] How to deform segmentation to compare registration result??
Bill Lorensen
bill.lorensen at gmail.com
Fri Nov 13 11:16:00 EST 2009
If your "segmented" data is represented by discrete integer values and
not a continuous greyscale, then you should use the
NearestNeighborInterpolateImageFunction. The
LinearInterpolateImageFunction will try to interpolate the integer
segmentation labels and produce garbage.
Bill
On Fri, Nov 13, 2009 at 8:13 AM, MINWOO PARK <mw770824 at yahoo.com> wrote:
> Hi, I just started to use ITK and am comparing the performance of the
> registration algorithms.
> When I deform a brain label segmentation, I have many artifacts in the
> segmentation.
>
> The procedure is as follows.
> Let’s say I have 5, 3D brain images with ground truth segmentations.
> So I register all the brains to one specific brain then using the
> deformation field found, I transform all the segmentation to the same
> specific brain. Then I measure how well the each segmented region is
> overlapping each other.
> But problem I have now is that after I transform the segmentation, I observe
> lots of salt and pepper like noise here there inside each region.
>
>
> One more question is...how long does 3D BSplineTransform usually take??? I
> have 7*7*7 landmarks and it takes around 2 minutes in intel 2.4 mobile cpu
> 2GB memory.
>
> Thank you!!!
>
> Here is the code for Bspline deformation that I am using .......
>
> const unsigned int ImageDimension = 3;
>
>
> typedef itk::ImageFileWriter< ITKIMAGE3D > MovingWriterType;
> MovingWriterType::Pointer movingWriter = MovingWriterType::New();
>
> typedef itk::ResampleImageFilter< ITKIMAGE3D,
> ITKIMAGE3D > FilterType;
> ITKIMAGE3D::Pointer fixedImage=itkimage_src;
>
> FilterType::Pointer resampler = FilterType::New();
> typedef itk::LinearInterpolateImageFunction<
> ITKIMAGE3D, double > InterpolatorType;
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> resampler->SetInterpolator( interpolator );
> ITKIMAGE3D::SpacingType fixedSpacing = fixedImage->GetSpacing();
> ITKIMAGE3D::PointType fixedOrigin = fixedImage->GetOrigin();
> ITKIMAGE3D::DirectionType fixedDirection = fixedImage->GetDirection();
>
> resampler->SetOutputSpacing( fixedSpacing );
> resampler->SetOutputOrigin( fixedOrigin );
> resampler->SetOutputDirection( fixedDirection );
>
> ITKIMAGE3D::RegionType fixedRegion =
> fixedImage->GetLargestPossibleRegion();
> ITKIMAGE3D::SizeType fixedSize = fixedRegion.GetSize();
> resampler->SetSize( fixedSize );
> resampler->SetOutputStartIndex( fixedRegion.GetIndex() );
> resampler->SetInput(itkimage_dst );
> movingWriter->SetInput( resampler->GetOutput() );
> movingWriter->SetFileName(outputPath);
> /*
> We instantiate now the type of the BSplineDeformableTransform using as
> template param-
> eters the type for coordinates representation, the dimension of the
> space, and the order of the
> B-spline.*/
>
> TransformType3D::Pointer bsplineTransform = TransformType3D::New();
> /*Since we are using a B-spline of order 3, the coverage of the BSpling
> grid should exceed by one
> the spatial extent of the image on the lower region of image indices,
> and by two grid points on
> the upper region of image indices. We choose here to use a 8
> × 8 B-spline grid, from which only
> a 5 by 5 sub-grid will be covering the input image. If we use an input
> image of size 500 by 500
> pixels, and pixel spacing 2.0 × 2.0 then we need the 5 × 5 B-spline
> grid to cover a physical
> extent of 1000 × 1000 mm. This can be achieved by setting the pixel
> spacing of the B-spline
> grid to 250.0 × 250.0 mm. The origin of the B-spline grid must be set
> at one grid position away
> from the origin of the desired output image. All this is done with the
> following lines of code. */
>
> typedef TransformType3D::RegionType RegionType;
> RegionType bsplineRegion; //itk::BSplineDeformableTransform::RegionType
> RegionType::SizeType size;
> //itk::BSplineDeformableTransform::RegionType::SizeType
> const unsigned int numberOfGridNodesOutsideTheImageSupport = 3;
> const unsigned int numberOfGridNodesInsideTheImageSupportH = nH;
> const unsigned int numberOfGridNodesInsideTheImageSupportW = nW;
> const unsigned int numberOfGridNodesInsideTheImageSupportD = nD;
> const unsigned int numberOfGridNodesH =
> numberOfGridNodesInsideTheImageSupportH +
> numberOfGridNodesOutsideTheImageSupport;
> const unsigned int numberOfGridNodesW =
> numberOfGridNodesInsideTheImageSupportW +
> numberOfGridNodesOutsideTheImageSupport;
> const unsigned int numberOfGridNodesD =
> numberOfGridNodesInsideTheImageSupportD +
> numberOfGridNodesOutsideTheImageSupport;
>
>
>
> //We need to modify this size[0]=m_nW size[1]=m_nH
> //itk::BSplineDeformableTransform::RegionType::SizeType.Fill
> size[0]=numberOfGridNodesW;
> size[1]=numberOfGridNodesH;
> size[2]=numberOfGridNodesD;
> bsplineRegion.SetSize( size );
> //itk::BSplineDeformableTransform::RegionType::SetSize
> typedef TransformType3D::SpacingType SpacingType;
> //itk::BSplineDeformableTransform::SpacingType
> SpacingType spacing;
> spacing[0] = nHorSpace*fixedSpacing[0];//floor( fixedSpacing[0] *
> (fixedSize[0] - 1) / numberOfGridCellsW ); //convert pixel spacing to
> metric....
> spacing[1] = nVerSpace*fixedSpacing[1];//floor( fixedSpacing[1] *
> (fixedSize[1] - 1) / numberOfGridCellsH );
> spacing[2] = nSliSpace*fixedSpacing[2];//floor( fixedSpacing[1] *
> (fixedSize[1] - 1) / numberOfGridCellsH );
>
> typedef TransformType3D::OriginType OriginType;
> //itk::BSplineDeformableTransform::OriginType
> OriginType origin;
> origin = fixedOrigin - fixedDirection*spacing; //due to boundary...
>
> bsplineTransform->SetGridSpacing( spacing );
> //itk::BSplineDeformableTransform->SetGridSpacing
> bsplineTransform->SetGridOrigin( origin );
> //itk::BSplineDeformableTransform->SetGridOrigin
> bsplineTransform->SetGridRegion( bsplineRegion
> );//itk::BSplineDeformableTransform->SetGridRegion
> bsplineTransform->SetGridDirection( fixedDirection );
>
> typedef TransformType3D::ParametersType ParametersType;
> //itk::BSplineDeformableTransform::ParametersType
> const unsigned int numberOfParameters =
> bsplineTransform->GetNumberOfParameters();
> //itk::BSplineDeformableTransform->GetNumberOfParameters();
> const unsigned int numberOfNodes = numberOfParameters / ImageDimension;
> ParametersType parameters( numberOfParameters );
>
> /*
> The B-spline grid should now be fed with coeficients at each node. Since
> this is a two dimen-
> sional grid, each node should receive two coefficients. Each coefficient
> pair is representing a
> displacement vector at this node. The coefficients can be passed to the
> B-spline in the form of
> an array where the first set of elements are the first component of the
> displacements for all the
> nodes, and the second set of elemets is formed by the second component
> of the displacements
> for all the nodes.
> In this example we read such displacements from a file, but for
> convinience we have written this
> file using the pairs of (x, y) displacement for every node. The elements
> read from the file should
> therefore be reorganized when assigned to the elements of the array. We
> do this by storing all
> the odd elements from the file in the first block of the array, and all
> the even elements from the
> file in the second block of the array. Finally the array is passed to
> the B-spline transform using
> the SetParameters(). */
>
> int nPos=0;
> for(int z=-1;z<nD+2;z++)
> {
> for(int y=-1;y<nH+2;y++)
> {
> for(int x=-1;x<nW+2;x++)
> {
> if(z>=0 &&z< nD && y>=0 && y< nH && x>=0 && x<nW)
> {
> parameters[nPos]=disVector[z*nW*nH+y*nW+x].x;
> parameters[nPos+numberOfNodes]=disVector[z*nW*nH+y*nW+x].y;
> parameters[nPos+numberOfNodes*2]=disVector[z*nW*nH+y*nW+x].z;
> }
> else
> {
> parameters[nPos]=0;//disVector[y*nW+x].x;
> parameters[nPos+numberOfNodes]=0;//disVector[y*nW+x].y;
> parameters[nPos+numberOfNodes*2]=0;
> }
> nPos++;
> }
> }
> }
> bsplineTransform->SetParameters( parameters );
>
> resampler->SetDefaultPixelValue(0);
> resampler->SetTransform( bsplineTransform );
> try
> {
> movingWriter->Update();
> }
> catch( itk::ExceptionObject & excp )
> {
> std::cerr << "Exception thrown " << std::endl;
> std::cerr << excp << std::endl;
> return EXIT_FAILURE;
> }
>
> //We just transform segmentation image to evaluate.....
> if(m_strSegPath!="")
> {
> if(nLevel==0)
> {//input is fixed segmentation image...
> ITK_IMG_READER_TYPE3D::Pointer
> seg_reader=ITK_IMG_READER_TYPE3D::New();
> seg_reader->SetFileName(m_strSegPath);
>
> //we first compute maximum and minimum of the label...and
> remember it....
> typedef itk::MinimumMaximumImageCalculator<ITKIMAGE3D> calType;
> calType::Pointer cal=calType::New();
> cal->SetImage(seg_reader->GetOutput());
> seg_reader->Update();
> cal->ComputeMaximum();
> cal->ComputeMinimum();
> ITKIMAGE3D::PixelType maxval=cal->GetMaximum();
> ITKIMAGE3D::PixelType minval=cal->GetMinimum();
>
>
> typedef itk::RescaleIntensityImageFilter<
> ITKIMAGE3D, ITKIMAGE3D> RescaleFilterType;
> RescaleFilterType::Pointer enlarge_rescaler =
> RescaleFilterType::New();
> //This rescaler was used to avoid interpolation ambiguity
> between adjacent labels.
> // in hope for avoiding salt and pepper like noise....
> enlarge_rescaler->SetOutputMinimum( 0 );
> enlarge_rescaler->SetOutputMaximum( 30000 );
> enlarge_rescaler->SetInput(seg_reader->GetOutput());
>
> resampler->SetInput(enlarge_rescaler->GetOutput());
> ITK_IMG_WRITER_TYPE3D::Pointer
> seg_writer=ITK_IMG_WRITER_TYPE3D::New();
> QString strFileName;
> strFileName.sprintf("%s%s_seg_deformed%.2d.hdr",qPrintable(m_strOutPath),qPrintable(m_strFileName),nLevel);
> seg_writer->SetFileName(strFileName);
>
> RescaleFilterType::Pointer shorten_rescaler =
> RescaleFilterType::New();
>
> shorten_rescaler->SetOutputMinimum(minval);
> shorten_rescaler->SetOutputMaximum(maxval);
> shorten_rescaler->SetInput(resampler->GetOutput());
> seg_writer->SetInput(shorten_rescaler->GetOutput());
> seg_writer->Update();
> }
> else
> {//input is transformed segmentation image....
> QString strFileName;
> strFileName.sprintf("%s%s_seg_deformed%.2d.hdr",qPrintable(m_strOutPath),qPrintable(m_strFileName),nLevel-1);
>
> ITK_IMG_READER_TYPE3D::Pointer
> seg_reader=ITK_IMG_READER_TYPE3D::New();
>
> //we first compute maximum and minimum of the label...and
> remember it....
> typedef itk::MinimumMaximumImageCalculator<ITKIMAGE3D> calType;
> calType::Pointer cal=calType::New();
> cal->SetImage(seg_reader->GetOutput());
> seg_reader->Update();
> cal->ComputeMaximum();
> cal->ComputeMinimum();
> ITKIMAGE3D::PixelType maxval=cal->GetMaximum();
> ITKIMAGE3D::PixelType minval=cal->GetMinimum();
>
>
> typedef itk::RescaleIntensityImageFilter<
> ITKIMAGE3D, ITKIMAGE3D> RescaleFilterType;
> RescaleFilterType::Pointer enlarge_rescaler =
> RescaleFilterType::New();
>
> enlarge_rescaler->SetOutputMinimum( 0 );
> enlarge_rescaler->SetOutputMaximum( 30000 );
> enlarge_rescaler->SetInput(seg_reader->GetOutput());
>
> resampler->SetInput(enlarge_rescaler->GetOutput());
>
> seg_reader->SetFileName(strFileName);
>
> ITK_IMG_WRITER_TYPE3D::Pointer
> seg_writer=ITK_IMG_WRITER_TYPE3D::New();
> strFileName.sprintf("%s%s_seg_deformed%.2d.hdr",qPrintable(m_strOutPath),qPrintable(m_strFileName),nLevel);
> seg_writer->SetFileName(strFileName);
>
> RescaleFilterType::Pointer shorten_rescaler =
> RescaleFilterType::New();
>
> shorten_rescaler->SetOutputMinimum(minval);
> shorten_rescaler->SetOutputMaximum(maxval);
> shorten_rescaler->SetInput(resampler->GetOutput());
> seg_writer->SetInput(shorten_rescaler->GetOutput());
> seg_writer->Update();
> }
> }
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>
More information about the Insight-users
mailing list