[Insight-users] Vital DICOM information/tags when registering images?
Bill Lorensen
bill.lorensen at gmail.com
Wed Sep 16 10:33:23 EDT 2009
Matlab stores images in Fortran (column) order. C/C++ in row order.
This may account for the 90 degree difference.
On Wed, Sep 16, 2009 at 9:27 AM, Patrik Brynolfsson
<patrik.brynolfsson at gmail.com> wrote:
> Hello Andriy,
> thanks for you tips. I reinitialize the seed every time, so the two
> different methods are consistent between runs but not between each other.
> I'm checking if the two different methods of reading the images result in
> the same pixel values right now, so I will report back on that when I have
> the results. The difference between the methods are huge for some images,
> and reading the images from files constantly generate correct results
> whereas the resuts when the images are passed from matlab are sometimes
> really bad (rotated 90 degrees or something like that).
> I was thinking of comparing the images using an iterator and comparing pixel
> by pixel. Is there another way of comparing the images that is easier?
> Thanx for your input.
> //Patrik
>
> 2009/9/16 Andriy Fedorov <fedorov at bwh.harvard.edu>
>>
>> Patrick,
>>
>> Registration has the metric sampling component, which is initialized
>> with a random seed by default. If you want to have consistent results
>> you might need to initialize registration with the same seed.
>>
>> Another guess is to check that there is no type cast involved when you
>> read the image in Matlab. This might result in loss of precision.
>>
>> Not sure if any of these is the source of the problem for you, just
>> some ideas to investigate.
>>
>> Andriy Fedorov
>>
>>
>>
>> > Date: Wed, 16 Sep 2009 09:02:59 +0200
>> > From: Patrik Brynolfsson <patrik.brynolfsson at gmail.com>
>> > Subject: Re: [Insight-users] Vital DICOM information/tags when
>> > registering images?
>> > To: Bill Lorensen <bill.lorensen at gmail.com>, ITK
>> > <insight-users at itk.org>
>> > Message-ID:
>> > <8d0068770909160002r1e0aa86el32f78b0d027c9368 at mail.gmail.com>
>> > Content-Type: text/plain; charset="iso-8859-1"
>> >
>> > Hello Bill, thanx for answering!
>> >
>> > The filter I use are:
>> > Registration part (same for both cases):
>> >
>> > itkImageRegistrationMethod
>> > itkMattesMutualInformationImageToImageMetric
>> > itkLinearInterpolateImageFunction
>> > itkVersorRigid3DTransform
>> > itkVersorRigid3DTransformOptimizer
>> > itkCenteredTransformInitializer
>> >
>> > Reading images from memory:
>> >
>> > itkImportImageFilter
>> > itkChangeInformationImageFilter
>> >
>> > Reading images from DICOM:
>> >
>> > itkImageSeriesReader
>> > itkGDCMImageIO
>> > itkGDCMSeriesFileNames
>> >
>> > *The code to read images from memory is:*
>> >
>> > typedef itk::ImportImageFilter<PixelType,Dimension> ImportFilterType;
>> > ImportFilterType::Pointer fixedImageReader = ImportFilterType::New();
>> > ImportFilterType::Pointer movingImageReader = ImportFilterType::New();
>> >
>> > const mwSize* dimsFix = mxGetDimensions(prhs[0]);
>> > const mwSize* dimsMov = mxGetDimensions(prhs[1]);
>> >
>> > ImportFilterType::SizeType sizeFix;
>> > ImportFilterType::SizeType sizeMov;
>> > sizeFix[0] = dimsFix[0];
>> > sizeFix[1] = dimsFix[1];
>> > sizeFix[2] = dimsFix[2];
>> >
>> > sizeMov[0] = dimsMov[0];
>> > sizeMov[1] = dimsMov[1];
>> > sizeMov[2] = dimsMov[2];
>> >
>> > ImportFilterType::IndexType start;
>> > start.Fill( 0 );
>> > ImportFilterType::RegionType regionFix;
>> > regionFix.SetIndex( start );
>> > regionFix.SetSize( sizeFix );
>> > ImportFilterType::RegionType regionMov;
>> > regionMov.SetIndex( start );
>> > regionMov.SetSize( sizeMov );
>> > fixedImageReader->SetRegion( regionFix );
>> > movingImageReader->SetRegion(regionMov );
>> >
>> > //Import image from memory, passed from MATLAB
>> >
>> > fixedImageReader->SetImportPointer((PixelType*)mxGetData(prhs[0]),mxGetNumberOfElements(prhs[0]),false);
>> >
>> > movingImageReader->SetImportPointer((PixelType*)mxGetData(prhs[1]),mxGetNumberOfElements(prhs[1]),false);
>> > fixedImageReader->Update();
>> > movingImageReader->Update();
>> >
>> > *...and the code to set the relevant info is:*
>> >
>> > double *fixedSpacingInput = mxGetPr(prhs[2]);
>> > double *movingSpacingInput = mxGetPr(prhs[3]);
>> >
>> > double *fixedOrientationInput = mxGetPr(prhs[4]);
>> > double *movingOrientationInput = mxGetPr(prhs[5]);
>> >
>> > double *fixedOriginInput = mxGetPr(prhs[6]);
>> > double *movingOriginInput = mxGetPr(prhs[7]);
>> >
>> > **Code for casting input to correct types are skipped, not important
>> > for
>> > this question**
>> >
>> > typedef itk::ChangeInformationImageFilter< FixedImageType >
>> > ChangeInfoType;
>> > ChangeInfoType::Pointer changeFixedInfo = ChangeInfoType::New();
>> > ChangeInfoType::Pointer changeMovingInfo = ChangeInfoType::New();
>> >
>> > changeFixedInfo->SetInput(fixedImageReader->GetOutput());
>> > changeMovingInfo->SetInput(movingImageReader->GetOutput());
>> >
>> > changeFixedInfo->SetOutputSpacing(fixedSpacing);
>> > changeMovingInfo->SetOutputSpacing(movingSpacing);
>> > changeFixedInfo->SetChangeSpacing(true);
>> > changeMovingInfo->SetChangeSpacing(true);
>> >
>> > changeFixedInfo->SetOutputDirection(fixedOrientation);
>> > changeMovingInfo->SetOutputDirection(movingOrientation);
>> > changeFixedInfo->SetChangeDirection(true);
>> > changeMovingInfo->SetChangeDirection(true);
>> >
>> > changeFixedInfo->SetOutputOrigin(fixedOrigin);
>> > changeMovingInfo->SetOutputOrigin(movingOrigin);
>> > changeFixedInfo->SetChangeOrigin(true);
>> > changeMovingInfo->SetChangeOrigin(true);
>> >
>> > changeFixedInfo->Update();
>> > changeMovingInfo->Update();
>> >
>> > registration->SetFixedImage( changeFixedInfo->GetOutput() );
>> > registration->SetMovingImage( changeMovingInfo->GetOutput() );
>> > registration->SetFixedImageRegion(
>> > changeFixedInfo->GetOutput()->GetBufferedRegion()
>> > );
>> >
>> > As a check I print the origin, spacing and direction for the images and
>> > they
>> > are identical to when reading from DICOM files. I read the DICOM series
>> > just
>> > as described in chapter 7.12 in the ITK manual.
>> >
>> > Any ideas to what might be wrong here?
>> >
>> > Thanks in advance!
>> >
>> > //Patrik Brynolfsson
>> >
>> >
>> >
>> > 2009/9/16 Bill Lorensen <bill.lorensen at gmail.com>
>> >
>> >> You should get the same answers. Exactly what filters are you using in
>> >> each of the cases. I can't tell from your description.
>> >>
>> >> Bill
>> >>
>> >> On Tue, Sep 15, 2009 at 4:38 PM, Patrik Brynolfsson
>> >> <patrik.brynolfsson at gmail.com> wrote:
>> >> > Hello,
>> >> > I'm working on a Matlab mex-implementation of an ITK registration.
>> >> > This
>> >> > means that instead of loading the images from files I will load them
>> >> > from
>> >> > memory as 3D matrices. All the DICOM information needs to be passed
>> >> > to
>> >> ITK
>> >> > from Matlab as well, and this is where my problems begin. The
>> >> > information
>> >> > that I thought was essential for a correct registration is:
>> >> > * Origin
>> >> > * Spacing
>> >> > * Direction
>> >> > I set these using the "ChangeInformationImageFilter" and it is
>> >> > working;
>> >> the
>> >> > images have the origin, spacing and directional cosines that are
>> >> > assigned
>> >> to
>> >> > them, and they are the same as if the dicom images were loaded from
>> >> files.
>> >> > However, the result of the registration is different from that when
>> >> loading
>> >> > the images from files, with all other components set the same, (of
>> >> > course
>> >> I
>> >> > reinitialize the seed). Not by much, but they should be exactly the
>> >> > same!
>> >> I
>> >> > use Mattes mutual information metric, linear interpolation, Versor
>> >> > rigid
>> >> 3D
>> >> > transform and the accompanying Versor rigid 3D optimizer together
>> >> > with
>> >> the
>> >> > "CenteredTransformInitializer".
>> >> > I know that the difference in results are due to the way I read the
>> >> > data.
>> >> If
>> >> > I change the input from memory to reading the dicom files I get the
>> >> > same
>> >> > result as with a stand-alone .exe-application, so this is where the
>> >> > error
>> >> > is.
>> >> > Does anyone have any idea as to what other information I need to pass
>> >> from
>> >> > matlab to get identical results?
>> >> > Thanks in advance
>> >> > --
>> >> > Patrik Brynolfsson
>> >> >
>> >> >
>> >> >
>> >> > _____________________________________
>> >> > Powered by www.kitware.com
>> >> >
>> >> > Visit other Kitware open-source projects at
>> >> > http://www.kitware.com/opensource/opensource.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
>> >> >
>> >> >
>> >>
>> >
>> >
>> >
>> > --
>> > Patrik Brynolfsson
>> > -------------- next part --------------
>> > An HTML attachment was scrubbed...
>> > URL:
>> > <http://www.itk.org/pipermail/insight-users/attachments/20090916/f0f7f731/attachment.htm>
>> >
>> > ------------------------------
>> >
>> > _______________________________________________
>> > Insight-users mailing list
>> > Insight-users at itk.org
>> > http://www.itk.org/mailman/listinfo/insight-users
>> >
>> >
>> > End of Insight-users Digest, Vol 65, Issue 56
>> > *********************************************
>> >
>
>
>
> --
> Patrik Brynolfsson
>
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.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