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