[Insight-users] Vital DICOM information/tags when registering images?
Patrik Brynolfsson
patrik.brynolfsson at gmail.com
Fri Sep 18 06:25:36 EDT 2009
I have found the source of the error now. There is a discrepancy in the
Origin, Spacing and Direction properties in the 6th decimal. If I set the
Origin, Spacing and Direction properties from the dicom images to the images
loaded from matlab I get identical results.
Is it possible that the dicom tags are read and handled in float and then
typecasted to double? I see no other reason why there should be a difference
in the 6th decimal place...
2009/9/17 Patrik Brynolfsson <patrik.brynolfsson at gmail.com>
> Yes, I set the seed, which make both methods consistent between runs but
> not between eachother. I'm thinking of changing the metric or using the
> whole image region instead of randomly assigning voxels, see if it gives
> better results. Tomorrow.
> Anyone has any idea what this might be? I only change the way the images
> are read into ITK, all other code is identical.
>
> All ideas appreciated!
>
> 2009/9/17 Bill Lorensen <bill.lorensen at gmail.com>
>
> How are you setting the random number seed? The MI metrics use a
>> random number generator.
>>
>> Bill
>>
>> On Thu, Sep 17, 2009 at 11:23 AM, Patrik Brynolfsson
>> <patrik.brynolfsson at gmail.com> wrote:
>> > Yes, the images were oriented differently, which I fixed. But the
>> problem is
>> > still there; there is still a discrepancy (although smaller) between the
>> > results. I used iterators to check if the two images contained the same
>> > information, which they did so I can't think of any reaon why there is a
>> > difference!
>> > To check that the information is the same in the image loaded from file
>> and
>> > the image passed by matlab I used this code:
>> > typedef itk::ImageRegionConstIterator<FixedImageType> ConstIterator1;
>> > typedef itk::ImageRegionConstIterator<FixedImageType> ConstIterator2;
>> > ConstIterator1 dicomIterator(fixedDicomReader->GetOutput(),
>> > fixedDicomReader->GetOutput()->GetRequestedRegion());
>> > ConstIterator2 matlabIterator(fixedMatlabReader->GetOutput(),
>> > fixedMatlabReader->GetOutput()->GetRequestedRegion());
>> > dicomIterator.GoToBegin();
>> > matlabIterator.GoToBegin();
>> > for(; !dicomIterator.IsAtEnd() ; ++dicomIterator, ++matlabIterator )
>> > {
>> > if( dicomIterator.Get() == matlabIterator.Get() )
>> > {}
>> > else
>> > {
>> > mexPrintf("Index: %i differs\n",dicomIterator.GetIndex());
>> > }
>> > }
>> > I used the image->Print() method to compare them as well, and the only
>> > difference there were the time stamps and the line:
>> > Matlab image:
>> > Container manages memory: false
>> > image from file:
>> > Container manages memory: true
>> > Is there anything else I can check that might be wrong? I am confused.
>> >
>> >
>> > 2009/9/16 Patrik Brynolfsson <patrik.brynolfsson at gmail.com>
>> >>
>> >> Exelent, I will check that out!
>> >>
>> >> 2009/9/16 Bill Lorensen <bill.lorensen at gmail.com>
>> >>>
>> >>> 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
>> >>> >
>> >>> >
>> >>
>> >>
>> >>
>> >> --
>> >> Patrik Brynolfsson
>> >>
>> >>
>> >
>> >
>> >
>> > --
>> > Patrik Brynolfsson
>> >
>> >
>> >
>>
>
>
>
> --
> Patrik Brynolfsson
>
>
>
--
Patrik Brynolfsson
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090918/e00f0662/attachment-0001.htm>
More information about the Insight-users
mailing list