[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