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

&gt;&gt; &gt;<br>
&gt;&gt; &gt; ------------------------------<br>
&gt;&gt; &gt;<br>
&gt;&gt; &gt; _______________________________________________<br>
&gt;&gt; &gt; Insight-users mailing list<br>
&gt;&gt; &gt; <a href="mailto:Insight-users@itk.org">Insight-users@itk.org</a><br>
&gt;&gt; &gt; <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
&gt;&gt; &gt;<br>
&gt;&gt; &gt;<br>
&gt;&gt; &gt; End of Insight-users Digest, Vol 65, Issue 56<br>
&gt;&gt; &gt; *********************************************<br>
&gt;&gt; &gt;<br>
&gt;<br>
&gt;<br>
&gt;<br>
&gt; --<br>
&gt; Patrik Brynolfsson<br>
&gt;<br>
&gt;<br>
&gt;<br>
&gt; _____________________________________<br>
&gt; Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
&gt;<br>
&gt; Visit other Kitware open-source projects at<br>
&gt; <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
&gt;<br>
&gt; Please keep messages on-topic and check the ITK FAQ at:<br>
&gt; <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
&gt;<br>
&gt; Follow this link to subscribe/unsubscribe:<br>
&gt; <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
&gt;<br>
&gt;<br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br>Patrik Brynolfsson<br><br><br>