[Insight-users] Vital DICOM information/tags when registering images?
Bill Lorensen
bill.lorensen at gmail.com
Thu Sep 17 13:19:40 EDT 2009
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
>
>
>
More information about the Insight-users
mailing list