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, <b>which they did</b> so I can't think of any reaon why there is a difference! <div>
<br></div><div>To check that the information is the same in the image loaded from file and the image passed by matlab I used this code:</div><div><br></div><div><div> typedef itk::ImageRegionConstIterator<FixedImageType> ConstIterator1;</div>
<div> typedef itk::ImageRegionConstIterator<FixedImageType> ConstIterator2;</div><div><br></div><div> ConstIterator1 dicomIterator(fixedDicomReader->GetOutput(),</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> fixedDicomReader->GetOutput()->GetRequestedRegion());</div>
<div> ConstIterator2 matlabIterator(fixedMatlabReader->GetOutput(),</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> fixedMatlabReader->GetOutput()->GetRequestedRegion());</div><div><br></div>
<div> dicomIterator.GoToBegin();</div><div> matlabIterator.GoToBegin();</div><div><br></div><div> for(; !dicomIterator.IsAtEnd() ; ++dicomIterator, ++matlabIterator )</div><div> {</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> if( dicomIterator.Get() == matlabIterator.Get() )</div>
<div><span class="Apple-tab-span" style="white-space:pre"> </span> {}</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> else</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> {</div>
<div><span class="Apple-tab-span" style="white-space:pre"> </span> mexPrintf("Index: %i differs\n",dicomIterator.GetIndex());</div><div><span class="Apple-tab-span" style="white-space:pre"> </span> }</div><div>
}</div><div>
<br></div><div>I used the image->Print() method to compare them as well, and the only difference there were the time stamps and the line:<div><br></div><div>Matlab image:</div><div>Container manages memory: false</div>
<div><br></div><div>image from file:</div><div>Container manages memory: true</div><div><br></div><div>Is there anything else I can check that might be wrong? I am confused.</div><div><br></div><div><br><br><div class="gmail_quote">
2009/9/16 Patrik Brynolfsson <span dir="ltr"><<a href="mailto:patrik.brynolfsson@gmail.com" target="_blank">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">
Exelent, I will check that out!<div><div></div><div><br><br><div class="gmail_quote">2009/9/16 Bill Lorensen <span dir="ltr"><<a href="mailto:bill.lorensen@gmail.com" target="_blank">bill.lorensen@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
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>
<div><div></div><div><<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 other.<br>
> I'm checking if the two different methods of reading the images result in<br>
> the same pixel values right now, so I will report back on that when I have<br>
> the results. The difference between the methods are huge for some 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 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> ImportFilterType;<br>
>> > ImportFilterType::Pointer fixedImageReader = ImportFilterType::New();<br>
>> > ImportFilterType::Pointer movingImageReader = 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>
>> > fixedImageReader->SetImportPointer((PixelType*)mxGetData(prhs[0]),mxGetNumberOfElements(prhs[0]),false);<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 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 and<br>
>> > they<br>
>> > are identical to when reading from DICOM files. I read the DICOM 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 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 registration.<br>
>> >> > This<br>
>> >> > means that instead of loading the images from files I will load them<br>
>> >> > from<br>
>> >> > memory as 3D matrices. All the DICOM information needs to be 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 from<br>
>> >> files.<br>
>> >> > However, the result of the registration is different from that when<br>
>> >> loading<br>
>> >> > the images from files, with all other components set the same, (of<br>
>> >> > course<br>
>> >> I<br>
>> >> > reinitialize the seed). Not by much, but they should be exactly the<br>
>> >> > same!<br>
>> >> I<br>
>> >> > use Mattes mutual information metric, linear interpolation, 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 the<br>
>> >> > data.<br>
>> >> If<br>
>> >> > I change the input from memory to reading the dicom files I get the<br>
>> >> > same<br>
>> >> > result as with a stand-alone .exe-application, so this is where the<br>
>> >> > error<br>
>> >> > is.<br>
>> >> > Does anyone have any idea as to what other information I need to 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>
>> > <<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>
</div></div></blockquote></div><br><br clear="all"><br></div></div>-- <br><font color="#888888">Patrik Brynolfsson<br><br><br>
</font></blockquote></div><br><br clear="all"><br>-- <br>Patrik Brynolfsson<br><br><br>
</div></div>
</div>