<br clear="all">Hi All,<br>I&#39;m trying to run the EM algorithm in ITK on my image.<br>My image is a simple scalar 3D volume.<br>I wrote the following code, based on the itkSoftwareGuide example, but it seems that the estimator does not run at all.<br>
The output estimated values are exactly the same as the initialization.<br>Thanks,<br>Moti<br><br><br>*********************************** Begining of code ****************************************************<br>&nbsp;&nbsp;&nbsp; unsigned int numberOfClasses = 2;<br>
&nbsp;&nbsp;&nbsp; typedef itk::Vector&lt; PixelType, 1 &gt; MeasurementVectorType;<br><br>&nbsp;&nbsp;&nbsp; typedef itk::Statistics::ScalarImageToListAdaptor &lt;ImageType&gt; itkScalarImageToListAdaptorType;<br>&nbsp;&nbsp;&nbsp; itkScalarImageToListAdaptorType::Pointer im2list = itkScalarImageToListAdaptorType::New();<br>
&nbsp;&nbsp;&nbsp; im2list-&gt;SetImage (sample_volume);<br><br><br><br>&nbsp;&nbsp;&nbsp; typedef itkScalarImageToListAdaptorType SampleType;<br>&nbsp;&nbsp;&nbsp; <br>&nbsp;&nbsp;&nbsp; <br>&nbsp;&nbsp;&nbsp; typedef itk::Array&lt; double &gt; ParametersType;<br>&nbsp;&nbsp;&nbsp; ParametersType params( 2 );<br>
<br>&nbsp;&nbsp;&nbsp; std::vector&lt; ParametersType &gt; initialParameters( numberOfClasses );<br>&nbsp;&nbsp;&nbsp; params[0] = 60;<br>&nbsp;&nbsp;&nbsp; params[1] = 50;<br>&nbsp;&nbsp;&nbsp; initialParameters[0] = params;<br><br>&nbsp;&nbsp;&nbsp; params[0] = 250;<br>&nbsp;&nbsp;&nbsp; params[1] = 50;<br>&nbsp;&nbsp;&nbsp; initialParameters[1] = params;<br>
<br>&nbsp;&nbsp;&nbsp; typedef itk::Statistics::GaussianMixtureModelComponent&lt; SampleType &gt;&nbsp; ComponentType;<br><br>&nbsp;&nbsp;&nbsp; std::vector&lt; ComponentType::Pointer &gt; components;<br>&nbsp;&nbsp;&nbsp; for ( unsigned int i = 0 ; i &lt; numberOfClasses ; i++ )<br>
&nbsp;&nbsp;&nbsp; {<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; components.push_back( ComponentType::New() );<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; (components[i])-&gt;SetSample( im2list );<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; (components[i])-&gt;SetParameters( initialParameters[i] );<br>&nbsp;&nbsp;&nbsp; }<br><br><br><br>&nbsp;&nbsp;&nbsp; typedef itk::Statistics::ExpectationMaximizationMixtureModelEstimator&lt; SampleType &gt; EstimatorType;<br>
&nbsp;&nbsp;&nbsp; EstimatorType::Pointer estimator = EstimatorType::New();<br>&nbsp;&nbsp;&nbsp; estimator-&gt;SetSample( im2list );<br>&nbsp;&nbsp;&nbsp; estimator-&gt;SetMaximumIteration( 200 );<br><br>&nbsp;&nbsp;&nbsp; itk::Array&lt; double &gt; initialProportions(numberOfClasses);<br>
&nbsp;&nbsp;&nbsp; initialProportions[0] = 0.5;<br>&nbsp;&nbsp;&nbsp; initialProportions[1] = 0.5;<br><br>&nbsp;&nbsp;&nbsp; estimator-&gt;SetInitialProportions( initialProportions );<br>&nbsp;&nbsp;&nbsp; &nbsp; <br>&nbsp;&nbsp;&nbsp; for ( unsigned int i = 0 ; i &lt; numberOfClasses ; i++)<br>&nbsp;&nbsp;&nbsp; {<br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; estimator-&gt;AddComponent( (ComponentType::Superclass*)<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; (components[i]).GetPointer() );<br>&nbsp;&nbsp;&nbsp; }<br>&nbsp;&nbsp;&nbsp; &nbsp; <br>&nbsp;&nbsp;&nbsp; estimator-&gt;Update();<br>&nbsp;&nbsp;&nbsp; for ( unsigned int i = 0 ; i &lt; numberOfClasses ; i++ )<br>
&nbsp;&nbsp;&nbsp; {<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; std::cout &lt;&lt; &quot;Cluster[&quot; &lt;&lt; i &lt;&lt; &quot;]&quot; &lt;&lt; std::endl;<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; std::cout &lt;&lt; &quot;&nbsp;&nbsp;&nbsp; Parameters:&quot; &lt;&lt; std::endl;<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; std::cout &lt;&lt; &quot;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &quot; &lt;&lt; (components[i])-&gt;GetFullParameters() <br>
&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; &lt;&lt; std::endl;<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; std::cout &lt;&lt; &quot;&nbsp;&nbsp;&nbsp; Proportion: &quot;;<br>&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp; std::cout &lt;&lt; &quot;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &quot; &lt;&lt; (*estimator-&gt;GetProportions())[i] &lt;&lt; std::endl;<br>&nbsp;&nbsp;&nbsp; }<br>
<br>************************************ end of code **********************************************************<br>-- <br>__<br>Moti Freiman, Ph.D Student.<br>Medical Image Processing and Computer-Assisted Surgery Laboratory.<br>
School of Computer Science and Engineering.<br>The Hebrew University of Jerusalem Givat Ram, Jerusalem 91904, Israel<br>Phone: +(972)-2-658-5371 (laboratory)<br>WWW site: <a href="http://www.cs.huji.ac.il/~freiman">http://www.cs.huji.ac.il/~freiman</a>