<html><body><div style="color:#000; background-color:#fff; font-family:HelveticaNeue, Helvetica Neue, Helvetica, Arial, Lucida Grande, sans-serif;font-size:10pt"><div>Hi Miro,</div><div><br></div><div>Thank you for the guide. I'm new to vtk so while messing around, I managed to apply the followings to a set of dicom images (in a very crude way) and save it as a .vtk file.</div><div>1. Marching Cubes - get mesh</div><div>2. Decimation - reduce the number of triangles</div><div>3. Smoothing - smooth the mesh</div><div>4. Cropping - select the largest region using connectivity filter</div><div><br></div><div>I guess my next step is to find out how to do pre-processing and refining of the defected mesh, and also making statistical models (if there are enough training images...).</div><div><br></div><div>By the way, here's my crude code based on the the &nbsp;MarchingCubes Example at http://www.vtk.org/Wiki/VTK/Examples/Cxx/Modelling/MarchingCubes". Hopefully
 it can be a small help to other novice like me out there.</div><div><br></div><div>/////////////////////////////////</div><div>int main(int argc, char *argv[])</div><div>{</div><div>&nbsp; vtkSmartPointer&lt;vtkImageData&gt; volume =</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkImageData&gt;::New();</div><div>&nbsp; double isoValue;</div><div><br></div><div>&nbsp; if (argc &lt; 3) &nbsp;&nbsp;</div><div>&nbsp; &nbsp; {</div><div><span class="Apple-tab-span" style="white-space:pre">        </span>// Create a sphere if there's no dicom input</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkSphereSource&gt; sphereSource =&nbsp;</div><div>&nbsp; &nbsp; &nbsp; vtkSmartPointer&lt;vtkSphereSource&gt;::New();</div><div>&nbsp; &nbsp; sphereSource-&gt;SetPhiResolution(20);</div><div>&nbsp; &nbsp; sphereSource-&gt;SetThetaResolution(20);</div><div>&nbsp; &nbsp; sphereSource-&gt;Update();</div><div>&nbsp; &nbsp;&nbsp;</div><div>&nbsp; &nbsp; double
 bounds[6];</div><div>&nbsp; &nbsp; sphereSource-&gt;GetOutput()-&gt;GetBounds(bounds);</div><div>&nbsp; &nbsp; for (unsigned int i = 0; i &lt; 6; i += 2)</div><div>&nbsp; &nbsp; &nbsp; {</div><div>&nbsp; &nbsp; &nbsp; double range = bounds[i+1] - bounds[i];</div><div>&nbsp; &nbsp; &nbsp; bounds[i] &nbsp; = bounds[i] - .1 * range;</div><div>&nbsp; &nbsp; &nbsp; bounds[i+1] = bounds[i+1] + .1 * range;</div><div>&nbsp; &nbsp; &nbsp; }</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkVoxelModeller&gt; voxelModeller =&nbsp;</div><div>&nbsp; &nbsp; &nbsp; vtkSmartPointer&lt;vtkVoxelModeller&gt;::New();</div><div>&nbsp; &nbsp; voxelModeller-&gt;SetSampleDimensions(50,50,50);</div><div>&nbsp; &nbsp; voxelModeller-&gt;SetModelBounds(bounds);</div><div>&nbsp; &nbsp; voxelModeller-&gt;SetScalarTypeToFloat();</div><div>&nbsp; &nbsp; voxelModeller-&gt;SetMaximumDistance(.1);</div><div>&nbsp; &nbsp;&nbsp;</div><div>&nbsp; &nbsp;
 voxelModeller-&gt;SetInputConnection(sphereSource-&gt;GetOutputPort());</div><div>&nbsp; &nbsp; voxelModeller-&gt;Update();</div><div>&nbsp; &nbsp; isoValue = 0.5;</div><div>&nbsp; &nbsp; volume-&gt;DeepCopy(voxelModeller-&gt;GetOutput());</div><div>&nbsp; &nbsp; }</div><div>&nbsp; else</div><div>&nbsp; &nbsp; {</div><div>&nbsp; &nbsp; // Read dicom image when dicom image folder path is specified</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkDICOMImageReader&gt; reader =</div><div>&nbsp; &nbsp; &nbsp; vtkSmartPointer&lt;vtkDICOMImageReader&gt;::New();</div><div>&nbsp; &nbsp; reader-&gt;SetDirectoryName(argv[1]);</div><div>&nbsp; &nbsp; reader-&gt;Update();</div><div>&nbsp; &nbsp; volume-&gt;DeepCopy(reader-&gt;GetOutput());</div><div>&nbsp; &nbsp; isoValue = atof(argv[2]);</div><div>&nbsp; &nbsp; }</div><div><br></div><div>&nbsp; // Use MarchingCubes to generate iso-surface</div><div>&nbsp; vtkSmartPointer&lt;vtkMarchingCubes&gt; surface
 =&nbsp;</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkMarchingCubes&gt;::New();</div><div><br></div><div>&nbsp; #if VTK_MAJOR_VERSION &lt;= 5</div><div>&nbsp; &nbsp; surface-&gt;SetInput(volume);</div><div>&nbsp; #else</div><div>&nbsp; &nbsp; surface-&gt;SetInputData(volume);</div><div>&nbsp; #endif</div><div>&nbsp; surface-&gt;ComputeNormalsOn();</div><div>&nbsp; surface-&gt;ComputeScalarsOn();</div><div>&nbsp; surface-&gt;SetValue(0, isoValue);</div><div>&nbsp; std::cout &lt;&lt; "Marching Cube finished..." &lt;&lt; std::endl;</div><div><br></div><div>&nbsp; // Create polydata from iso-surface</div><div>&nbsp; vtkSmartPointer&lt;vtkPolyData&gt; marched =</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkPolyData&gt;::New();</div><div>&nbsp; surface-&gt;Update();</div><div>&nbsp; marched-&gt;DeepCopy(surface-&gt;GetOutput());</div><div>&nbsp; std::cout&lt;&lt;"Number of points: " &lt;&lt; marched-&gt;GetNumberOfPoints() &lt;&lt;
 std::endl;</div><div>&nbsp;&nbsp;</div><div>&nbsp; // Decimation to reduce the number of triangles</div><div>&nbsp; vtkSmartPointer&lt;vtkDecimatePro&gt; decimator =&nbsp;</div><div>&nbsp; &nbsp; vtkDecimatePro::New();</div><div>&nbsp; decimator-&gt;SetInputData(marched);</div><div>&nbsp; decimator-&gt;SetTargetReduction(0.5);</div><div>&nbsp; decimator-&gt;SetPreserveTopology(1);</div><div>&nbsp; decimator-&gt;Update();</div><div>&nbsp; std::cout &lt;&lt; "Decimation finished...." &lt;&lt; std::endl;</div><div><span class="Apple-tab-span" style="white-space:pre">                </span></div><div>&nbsp; // Smoothing</div><div>&nbsp; vtkSmartPointer&lt;vtkSmoothPolyDataFilter&gt; smoother = vtkSmartPointer&lt;vtkSmoothPolyDataFilter&gt;::New();</div><div>&nbsp; smoother-&gt;SetInputData(decimator-&gt;GetOutput());</div><div>&nbsp; smoother-&gt;SetNumberOfIterations(5);</div><div>&nbsp; smoother-&gt;SetFeatureAngle(60);</div><div>&nbsp;
 smoother-&gt;SetRelaxationFactor(0.05);</div><div>&nbsp; smoother-&gt;FeatureEdgeSmoothingOff();</div><div>&nbsp; std::cout &lt;&lt; "Smoothing mesh finished...." &lt;&lt; std::endl;</div><div><br></div><div>&nbsp;// Select the largest region</div><div>&nbsp; vtkSmartPointer&lt;vtkPolyDataConnectivityFilter&gt; connectivityFilter =</div><div><span class="Apple-tab-span" style="white-space:pre">        </span>vtkSmartPointer&lt;vtkPolyDataConnectivityFilter&gt;::New();</div><div>&nbsp; connectivityFilter-&gt;SetInputConnection(smoother-&gt;GetOutputPort());</div><div>&nbsp; connectivityFilter-&gt;ScalarConnectivityOff();</div><div>&nbsp; connectivityFilter-&gt;SetExtractionModeToLargestRegion();</div><div>&nbsp; connectivityFilter-&gt;Update();</div><div>&nbsp;&nbsp;</div><div>&nbsp;// Create final polygon mesh</div><div>&nbsp; vtkSmartPointer&lt;vtkPolyData&gt; mesh =</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkPolyData&gt;::New();</div><div>&nbsp;
 mesh-&gt;ShallowCopy(connectivityFilter-&gt;GetOutput());</div><div>&nbsp; std::cout&lt;&lt;"Number of points in the final polygon: " &lt;&lt; mesh-&gt;GetNumberOfPoints() &lt;&lt; std::endl;</div><div>&nbsp;</div><div>&nbsp; // Visualization</div><div>&nbsp; vtkSmartPointer&lt;vtkRenderer&gt; renderer =&nbsp;</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkRenderer&gt;::New();</div><div>&nbsp; renderer-&gt;SetBackground(.1, .2, .3);</div><div><br></div><div>&nbsp; vtkSmartPointer&lt;vtkRenderWindow&gt; renderWindow =&nbsp;</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkRenderWindow&gt;::New();</div><div>&nbsp; renderWindow-&gt;AddRenderer(renderer);</div><div>&nbsp; vtkSmartPointer&lt;vtkRenderWindowInteractor&gt; interactor =&nbsp;</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkRenderWindowInteractor&gt;::New();</div><div>&nbsp; interactor-&gt;SetRenderWindow(renderWindow);</div><div><br></div><div>&nbsp; vtkSmartPointer&lt;vtkPolyDataMapper&gt; mapper
 =&nbsp;</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkPolyDataMapper&gt;::New();</div><div>&nbsp; mapper-&gt;SetInputData(mesh);</div><div>&nbsp; mapper-&gt;ScalarVisibilityOff();</div><div><br></div><div>&nbsp; vtkSmartPointer&lt;vtkActor&gt; actor =&nbsp;</div><div>&nbsp; &nbsp; vtkSmartPointer&lt;vtkActor&gt;::New();</div><div>&nbsp; actor-&gt;SetMapper(mapper);</div><div><br></div><div>&nbsp; renderer-&gt;AddActor(actor);</div><div><br></div><div>&nbsp; renderWindow-&gt;Render();</div><div>&nbsp; interactor-&gt;Start();</div><div><br></div><div>&nbsp; // Save the mesh as a .vtk file</div><div>&nbsp; vtkSmartPointer&lt;vtkPolyDataWriter&gt; writer = vtkPolyDataWriter::New();</div><div>&nbsp; writer-&gt;SetFileName("mesh.vtk");</div><div>&nbsp; writer-&gt;SetInputData(mesh);</div><div>&nbsp; writer-&gt;SetFileTypeToASCII();</div><div>&nbsp; writer-&gt;Write();</div><div>&nbsp;&nbsp;</div><div>&nbsp; return
 EXIT_SUCCESS;</div><div>}</div><div><br></div><div><br></div><div>Best Regards,</div><div><span></span></div><div>Harold</div><div class="yahoo_quoted" style="display: block;"> <br> <br> <div style="font-family: HelveticaNeue, 'Helvetica Neue', Helvetica, Arial, 'Lucida Grande', sans-serif; font-size: 10pt;"> <div style="font-family: HelveticaNeue, 'Helvetica Neue', Helvetica, Arial, 'Lucida Grande', sans-serif; font-size: 12pt;"> <div dir="ltr"> <font size="2" face="Arial"> On Wednesday, November 27, 2013 12:06 PM, Miro Drahos &lt;mdrahos@robodoc.com&gt; wrote:<br> </font> </div>  <div class="y_msg_container">This is basically a segmentation problem, which in general is a <br clear="none">difficult problem.<br clear="none">One such pipeline could be as follows:<br clear="none">1.) image preprocessing (cropping, smoothing, etc).<br clear="none">2.) thresholding to get a filled binary mask of your object (vertebra)<br clear="none">3.) running Marching
 Cubes to get a mesh.<br clear="none"><br clear="none">This is a very crude and basic approach, but something to get you <br clear="none">started. Note that marching cubes might give you a mesh with topological <br clear="none">defects (holes, handles), from the noise in your data.<br clear="none"><br clear="none">There are a variety of other sophisticated algorithms, e.g. using <br clear="none">statistical methods if you have some already segmented data (like <br clear="none">creating a mean shape from the previously segmented data, and using PCA <br clear="none">to do shape analysis, or some machine learning approach).<br clear="none"><br clear="none">Another approach would be to use deformable models, i.e. define forces <br clear="none">that would (iteratively) drive an initial mesh to fit the image. The <br clear="none">forces would be based on image intensity, gradient and/or laplacian, and <br clear="none">internal forces to keep the mesh regular
 and relatively smooth.<br clear="none"><br clear="none">There is a lot that has been done in this domain, so look at some <br clear="none">articles on segmentation.<br clear="none"><br clear="none">HTH,<br clear="none">Miro<br clear="none"><br clear="none"><div class="yqt5310183674" id="yqtfd95729"><br clear="none">On 11/25/2013 06:20 AM, Boogiedoll wrote:<br clear="none">&gt; Hi,<br clear="none">&gt; I would like to ask the experts about common way in vtk to create a spine bone (e.g L1) mesh from a CT DICOM stack, including the pre-processing.<br clear="none">&gt;<br clear="none">&gt; Thanks &amp; Regards,<br clear="none">&gt; Harold<br clear="none">&gt; _______________________________________________<br clear="none">&gt; Powered by www.kitware.com<br clear="none">&gt;<br clear="none">&gt; Visit other Kitware open-source projects at <a shape="rect" href="http://www.kitware.com/opensource/opensource.html"
 target="_blank">http://www.kitware.com/opensource/opensource.html</a><br clear="none">&gt;<br clear="none">&gt; Please keep messages on-topic and check the VTK FAQ at: <a shape="rect" href="http://www.vtk.org/Wiki/VTK_FAQ" target="_blank">http://www.vtk.org/Wiki/VTK_FAQ</a><br clear="none">&gt;<br clear="none">&gt; Follow this link to subscribe/unsubscribe:<br clear="none">&gt; <a shape="rect" href="http://www.vtk.org/mailman/listinfo/vtkusers" target="_blank">http://www.vtk.org/mailman/listinfo/vtkusers</a><br clear="none"><br clear="none"></div><br><br></div>  </div> </div>  </div> </div></body></html>