[Insight-users] How to deform segmentation to compare registration result??
Bill Lorensen
bill.lorensen at gmail.com
Fri Nov 13 12:43:31 EST 2009
What will you do when they overlap?
On Fri, Nov 13, 2009 at 12:39 PM, MINWOO PARK <mw770824 at yahoo.com> wrote:
>
> 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