[Insight-users] Vital DICOM information/tags when registering images?
Patrik Brynolfsson
patrik.brynolfsson at gmail.com
Thu Sep 17 16:48:33 EDT 2009
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090917/638515a7/attachment-0001.htm>
More information about the Insight-users
mailing list