[Insight-users] How to deform segmentation to compare registration result??
MINWOO PARK
mw770824 at yahoo.com
Fri Nov 13 12:39:43 EST 2009
I tried and it doesn't have salt and pepper like noise any more.
However NN interpolate assigns spatially close pixel value to the unknown
pixel, segmentation label can jump around from labels to labels.
So for my purpose, I think I need to divide all the labels, make binary
images (if I have 10 different labels then 10 binary images), transform each
binary image, and measure overlapping. But I think that is really time
consuming procedure.
Do you have any idea???
Thanks!!!
On 11/13/09 11:16 AM, "Bill Lorensen" <bill.lorensen at gmail.com> wrote:
> 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",qPrint
> able(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",qPrint
> able(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",qPrint
> able(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