<div dir="ltr"><div>Hi Richard,</div><div><br></div><div>I'm only going to touch on one part of your email now... hopefully I'll be able to</div><div>address some of the other points later.</div><div><br></div><div>You cannot generate the correct 3D surface from the image data alone. The</div><div>DICOM images are oriented images, which means that, if you want to know</div><div>where things are, you need not only the image data, but also "orientation"</div><div>and "position" information for the images. And of course you are already</div><div>familiar with this... so let's jump ahead to what you actually need to _do_ in</div><div>order to generate the correct surface.</div><div><br></div><div>For the vtkDICOMReader, it should be easy:</div><div>1) read the image (obviously!)<br></div><div>2) call the GetPatientMatrix() method to get position/orientation info</div><div>3) run marching cubes on the image to get a surface</div><div>4) correct the surface by using vtkTransformPolyDataFilter to apply the</div><div> aforementioned PatientMatrix to the surface</div><div><br></div><div>These same basic steps can be used with any of the DICOM readers,</div><div>but of course each reader is slightly different and the devil is in the</div><div>details.</div><div><br></div><div>For example, here are some tricky details for vtkDICOMReader (and it is</div><div>the least "tricky" of the three readers, IMHO):</div><div><br></div><div>1) If you look at the orientation returned by vtkDICOMReader::GetPatientMatrix()</div><div>and compare it to the DICOM file metadata, you will notice that it is different!</div><div>This is because the vtkDICOMReader (like other VTK readers) rasterizes its</div><div>output from bottom to top, to historically match the glDrawPixels function.</div><div>Since native DICOM is rasterized top-to-bottom, this requires a vertical flip of</div><div>the orientation matrix.</div><div><br></div><div>For this reason, I strongly recommend that anyone using vtkDICOMReader in</div><div>a "serious" application call the SetMemoryRowOrderToFileNative() method,</div><div>which will cause the reader to rasterize the data from top-to-bottom and will</div><div>cause the matrix returned by GetPatientMatrix() to exactly match the original</div><div>DICOM meta data.</div><div><br></div><div>2) The vtkDICOMReader always uses IPP sorting for a volumetric series,</div><div>unless you specifically call SortingOff(), in which case it uses the order of</div><div>the files. I do not recommend calling SortingOff(). Relying on IPP is best.</div><div><br></div><div>3) For isosurfaces, note that RescaleIntercept/RescaleSlope can vary from</div><div>slice-to-slice in CT and PET. For such images, I wrote a filter called</div><div>vtkDICOMApplyRescale that can be applied to the data immediately after<br></div><div>the reader (but please read the doxygen docs for this filter before using it).</div><div><br></div><div>4) An even more pernicious problem for CTs is that, if the scan is done with</div><div>a tilted gantry, the resulting series does not form a rectangular "volume".</div><div>For such images, the vtkDICOMCTRectifier can resample the data so that</div><div>it can be properly used within VTK.</div><div><br></div><div>Hopefully this info will be of some help. Let me know if you want</div><div>clarification on any of the above points.</div><div><br></div><div>Cheers,</div><div> - David</div><div><br></div><div><br></div><br><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Apr 1, 2016 at 6:52 AM, mbcx9rb9 <span dir="ltr"><<a href="mailto:richard.j.brown@live.co.uk" target="_blank">richard.j.brown@live.co.uk</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Hi all,<br>
<br>
I'm having some troubles with orientation and location of my DICOMs when<br>
loaded with vtkDICOMReader and vtkGDCMImageReader.<br>
<br>
To illustrate my problems, I opened the same files with both the readers and<br>
performed marching cubes on them to display 3D surfaces. I also displayed an<br>
annotated orientation cube to allow me to compare these surfaces with<br>
Slicer3D (which I'm assuming correctly calculates orientation etc.).<br>
<br>
In the ideal world, the 3D surfaces should superpose perfectly, but this is<br>
not the case.<br>
<br>
My pipeline for reading DICOMs with vtkGDCMImageReader is:<br>
gdcm::IPPSorter <-- Calculates the proper<br>
z-spacing<br>
vtkGDCMImageReader<br>
vtkImageChangeInfromation <-- Applies the corrected z-spacing that<br>
was calculated by IPPSorter<br>
<br>
My pipeline for reading DICOMs with vtkDICOMReader is simply:<br>
vtkDICOMReader<br>
and then actor->SetUserMatrix(reader->GetPatientMatrix())<br>
<br>
My problem is that the 3D structures are not in the same position nor the<br>
same orientation. Furthermore, one might agree with Slicer3D for one image<br>
but not for the next...<br>
<br>
I guess my main question is - What is the correct way to systematically open<br>
DICOM files (using either vtkDICOMReader or vtkGDCMImageReader, I'm not<br>
fussy)?<br>
<br>
Below are two images, one of a human skull and one of a pig skull (chosen at<br>
random). Each time, the red is rendered from vtkGDCMImageReader, the black<br>
is rendered from vtkDICOMReader (turns black when I use SetUserMatrix, not<br>
sure why...) and the last is from Slicer3D.<br>
<br>
With the human skull, we see that vtkGDCMImageReader agrees with slicer,<br>
whereas vtkDICOMReader is upside-down (Superior-Inferior reflection).<br>
However for the pig, vtkGDCMImageReader is upside-down (Superior-Inferior<br>
reflection) (PLUS something funny has happened to the spacing), whereas<br>
vtkDICOMReader is facing the anterior direction as opposed to the superior<br>
direction.<br>
<br>
<<a href="http://vtk.1045678.n5.nabble.com/file/n5737481/Screen_Shot_2016-04-01_at_12.png" rel="noreferrer" target="_blank">http://vtk.1045678.n5.nabble.com/file/n5737481/Screen_Shot_2016-04-01_at_12.png</a>><br>
<<a href="http://vtk.1045678.n5.nabble.com/file/n5737481/Screen_Shot_2016-04-01_at_14.png" rel="noreferrer" target="_blank">http://vtk.1045678.n5.nabble.com/file/n5737481/Screen_Shot_2016-04-01_at_14.png</a>><br>
<br>
Below is the code I have used for both. Any pointers on the correct way to<br>
consistently read DICOM images correctly would be appreciated.<br>
<br>
Regards,<br>
Richard<br>
<br>
P.S. A very important question - is my vtkAnnotatedCubeActor correctly<br>
annotated? If not, the comparison between my 3D structures and Slicer3D is<br>
useless until this is corrected...<br>
<br>
#include <vtkSmartPointer.h><br>
#include <vtkDICOMReader.h><br>
#include <vtkDICOMFileSorter.h><br>
#include <vtkStringArray.h><br>
#include <vtkRenderWindow.h><br>
#include <vtkRenderer.h><br>
#include <vtkRenderWindowInteractor.h><br>
#include <vtkInteractorStyleTrackballCamera.h><br>
#include <vtkMarchingCubes.h><br>
#include <vtkPolyDataNormals.h><br>
#include <vtkStripper.h><br>
#include <vtkPolyDataMapper.h><br>
#include <vtkProperty.h><br>
#include <vtkImageChangeInformation.h><br>
#include <vtkAnnotatedCubeActor.h><br>
#include <vtkOrientationMarkerWidget.h><br>
<br>
#include <gdcmIPPSorter.h><br>
#include <vtkGDCMImageReader.h><br>
<br>
#include <QStringList><br>
#include <QDirIterator><br>
<br>
int main()<br>
{<br>
<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
//<br>
// Get filenames<br>
//<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
<br>
// Get all DCM files in directory<br>
vtkSmartPointer<vtkStringArray> filenames =<br>
vtkSmartPointer<vtkStringArray>::New();<br>
QDirIterator iterator("/Scan/Directory/",<br>
QStringList("*.dcm"), QDir::Files,<br>
QDirIterator::Subdirectories);<br>
while (iterator.hasNext()) {<br>
filenames->InsertNextValue(iterator.next().toUtf8().data());<br>
}<br>
<br>
// Sort<br>
vtkSmartPointer <vtkDICOMFileSorter> sorter = vtkSmartPointer<br>
<vtkDICOMFileSorter>::New();<br>
sorter->SetInputFileNames(filenames);<br>
sorter->Update();<br>
<br>
// Get most numerous series<br>
int numSeries = sorter->GetNumberOfSeries();<br>
int mostNumerousSeriesIndex = 0;<br>
<br>
for (int i=0; i<numSeries; i++) {<br>
if (sorter->GetFileNamesForSeries(i)->GetNumberOfValues() ><br>
sorter->GetFileNamesForSeries(mostNumerousSeriesIndex)->GetNumberOfValues())<br>
{<br>
mostNumerousSeriesIndex = i;<br>
}<br>
}<br>
<br>
vtkStringArray *seriesFilenames =<br>
sorter->GetFileNamesForSeries(mostNumerousSeriesIndex);<br>
<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
//<br>
// Setup renderer etc.<br>
//<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
<br>
vtkSmartPointer<vtkRenderer> renderer =<br>
vtkSmartPointer<vtkRenderer>::New();<br>
renderer->SetGradientBackground(1);<br>
renderer->SetBackground(.5,.5,.5);<br>
renderer->SetBackground2(0,0,0);<br>
<br>
vtkSmartPointer<vtkRenderWindow> renWin =<br>
vtkSmartPointer<vtkRenderWindow>::New();<br>
renWin->AddRenderer(renderer);<br>
renWin->SetSize(1000,1000);<br>
<br>
vtkSmartPointer<vtkInteractorStyleTrackballCamera> interactorStyle =<br>
vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New();<br>
<br>
vtkSmartPointer<vtkRenderWindowInteractor> renWinInt =<br>
vtkSmartPointer<vtkRenderWindowInteractor>::New();<br>
renWinInt->SetRenderWindow(renWin);<br>
renWinInt->SetInteractorStyle(interactorStyle);<br>
<br>
<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
//<br>
// vtkGDCMImageReader<br>
//<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
<br>
// need to convert from vtkStringArray to std::vector<std::string> to be<br>
able to sort...<br>
std::vector<std::string> seriesFilenames_stdVector;<br>
for (int i=0; i<seriesFilenames->GetNumberOfValues(); i++) {<br>
<br>
seriesFilenames_stdVector.push_back(seriesFilenames->GetValue(i).c_str());<br>
}<br>
<br>
// sort and get z-spacing<br>
gdcm::IPPSorter IPPs;<br>
IPPs.SetComputeZSpacing(true);<br>
IPPs.SetZSpacingTolerance(1e-3);<br>
IPPs.Sort(seriesFilenames_stdVector);<br>
<br>
// Convert the sorted std::vector back to vtkStringArray to be able to<br>
use with vtkGDCMImageReader...<br>
vtkSmartPointer<vtkStringArray> seriesFilenames_vtkStringArray_sorted =<br>
vtkSmartPointer<vtkStringArray>::New();<br>
for (size_t i=0; i<IPPs.GetFilenames().size(); i++) {<br>
<br>
seriesFilenames_vtkStringArray_sorted->InsertNextValue(IPPs.GetFilenames().at(i));<br>
}<br>
<br>
// GDCM reader<br>
vtkSmartPointer<vtkGDCMImageReader> gdcmReader =<br>
vtkSmartPointer<vtkGDCMImageReader>::New();<br>
gdcmReader->SetFileNames(seriesFilenames_vtkStringArray_sorted);<br>
gdcmReader->FileLowerLeftOn();<br>
<br>
// correct z-spacing<br>
vtkSmartPointer<vtkImageChangeInformation> changeInfo =<br>
vtkSmartPointer<vtkImageChangeInformation>::New();<br>
changeInfo->SetInputConnection(gdcmReader->GetOutputPort());<br>
changeInfo->SetOutputSpacing(gdcmReader->GetOutput()->GetSpacing()[0],<br>
gdcmReader->GetOutput()->GetSpacing()[1], IPPs.GetZSpacing());<br>
changeInfo->Update();<br>
<br>
// Marching cubes<br>
vtkSmartPointer<vtkMarchingCubes> gdcmMC =<br>
vtkSmartPointer<vtkMarchingCubes>::New();<br>
gdcmMC->SetInputConnection(changeInfo->GetOutputPort());<br>
gdcmMC->SetNumberOfContours(1);<br>
gdcmMC->GenerateValues(1,700.,3000.);<br>
<br>
vtkSmartPointer<vtkPolyDataNormals> gdcmNormals =<br>
vtkSmartPointer<vtkPolyDataNormals>::New();<br>
gdcmNormals->SetInputConnection(gdcmMC->GetOutputPort());<br>
gdcmNormals->SetFeatureAngle(60.0);<br>
<br>
vtkSmartPointer<vtkStripper> gdcmStripper =<br>
vtkSmartPointer<vtkStripper>::New();<br>
gdcmStripper->SetInputConnection(gdcmNormals->GetOutputPort());<br>
<br>
vtkSmartPointer<vtkPolyDataMapper> gdcmMapper =<br>
vtkSmartPointer<vtkPolyDataMapper>::New();<br>
gdcmMapper->SetInputConnection(gdcmStripper->GetOutputPort());<br>
gdcmMapper->ScalarVisibilityOff();<br>
gdcmMapper->Update();<br>
<br>
vtkSmartPointer<vtkActor> gdcmActor = vtkSmartPointer<vtkActor>::New();<br>
gdcmActor->SetMapper(gdcmMapper);<br>
gdcmActor->GetProperty()->SetColor(1.0,0.0,0.0);<br>
gdcmActor->GetProperty()->SetOpacity(1.);<br>
<br>
renderer->AddActor(gdcmActor);<br>
<br>
<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
//<br>
// vtkDICOMReader<br>
//<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
<br>
// DCM reader<br>
vtkSmartPointer<vtkDICOMReader> dcmReader =<br>
vtkSmartPointer<vtkDICOMReader>::New();<br>
dcmReader->SetFileNames(seriesFilenames);<br>
dcmReader->Update();<br>
<br>
// Marching cubes<br>
vtkSmartPointer<vtkMarchingCubes> dcmMC =<br>
vtkSmartPointer<vtkMarchingCubes>::New();<br>
dcmMC->SetInputConnection(dcmReader->GetOutputPort());<br>
dcmMC->SetNumberOfContours(1);<br>
dcmMC->GenerateValues(1,700.,3000.);<br>
<br>
vtkSmartPointer<vtkPolyDataNormals> dcmNormals =<br>
vtkSmartPointer<vtkPolyDataNormals>::New();<br>
dcmNormals->SetInputConnection(dcmMC->GetOutputPort());<br>
dcmNormals->SetFeatureAngle(60.0);<br>
<br>
vtkSmartPointer<vtkStripper> dcmStripper =<br>
vtkSmartPointer<vtkStripper>::New();<br>
dcmStripper->SetInputConnection(dcmNormals->GetOutputPort());<br>
<br>
vtkSmartPointer<vtkPolyDataMapper> dcmMapper =<br>
vtkSmartPointer<vtkPolyDataMapper>::New();<br>
dcmMapper->SetInputConnection(dcmStripper->GetOutputPort());<br>
dcmMapper->ScalarVisibilityOff();<br>
dcmMapper->Update();<br>
<br>
vtkSmartPointer<vtkActor> dcmActor = vtkSmartPointer<vtkActor>::New();<br>
dcmActor->SetUserMatrix(dcmReader->GetPatientMatrix());<br>
dcmActor->SetMapper(dcmMapper);<br>
dcmActor->GetProperty()->SetColor(0.0,1.0,0.0);<br>
dcmActor->GetProperty()->SetOpacity(1.);<br>
<br>
renderer->AddActor(dcmActor);<br>
<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
//<br>
// Orientation cube<br>
//<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
<br>
vtkSmartPointer<vtkAnnotatedCubeActor> annotatedCubeActor =<br>
vtkSmartPointer<vtkAnnotatedCubeActor>::New();<br>
annotatedCubeActor->SetXPlusFaceText ("L (+X)");<br>
annotatedCubeActor->SetXMinusFaceText("R (-X)");<br>
annotatedCubeActor->SetYPlusFaceText ("P (+Y)");<br>
annotatedCubeActor->SetYMinusFaceText("A (-Y)");<br>
annotatedCubeActor->SetZPlusFaceText ("S (+Z)");<br>
annotatedCubeActor->SetZMinusFaceText("I (-Z)");<br>
annotatedCubeActor->SetZFaceTextRotation(90);<br>
annotatedCubeActor->GetCubeProperty()->SetColor(0.5, 1.0, 1.0);<br>
annotatedCubeActor->GetTextEdgesProperty()->SetLineWidth(1.0);<br>
annotatedCubeActor->GetTextEdgesProperty()->SetDiffuse(0.0);<br>
annotatedCubeActor->GetTextEdgesProperty()->SetAmbient(1.0);<br>
annotatedCubeActor->GetTextEdgesProperty()->SetColor(0.18, 0.28, 0.23);<br>
annotatedCubeActor->SetFaceTextScale(0.2);<br>
<br>
annotatedCubeActor->GetXPlusFaceProperty()->SetColor(0.0,0.0,1.0);<br>
annotatedCubeActor->GetXPlusFaceProperty()->SetInterpolationToFlat();<br>
annotatedCubeActor->GetXMinusFaceProperty()->SetColor(0.0,0.0,1.0);<br>
annotatedCubeActor->GetXMinusFaceProperty()->SetInterpolationToFlat();<br>
<br>
annotatedCubeActor->GetYPlusFaceProperty()->SetColor(0.0,1.0,0.0);<br>
annotatedCubeActor->GetYPlusFaceProperty()->SetInterpolationToFlat();<br>
annotatedCubeActor->GetYMinusFaceProperty()->SetColor(0.0,1.0,0.0);<br>
annotatedCubeActor->GetYMinusFaceProperty()->SetInterpolationToFlat();<br>
<br>
annotatedCubeActor->GetZPlusFaceProperty()->SetColor(1.0,0.0,0.0);<br>
annotatedCubeActor->GetZPlusFaceProperty()->SetInterpolationToFlat();<br>
annotatedCubeActor->GetZMinusFaceProperty()->SetColor(1.0,0.0,0.0);<br>
annotatedCubeActor->GetZMinusFaceProperty()->SetInterpolationToFlat();<br>
<br>
vtkSmartPointer<vtkOrientationMarkerWidget> widget =<br>
vtkSmartPointer<vtkOrientationMarkerWidget>::New();<br>
widget->SetOutlineColor( 0.9300, 0.5700, 0.1300 );<br>
widget->SetOrientationMarker( annotatedCubeActor );<br>
widget->SetInteractor( renWinInt );<br>
widget->SetViewport( 0.0, 0.0, 0.2, 0.2 );<br>
widget->SetEnabled( 1 );<br>
widget->InteractiveOn();<br>
<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
//<br>
// Render<br>
//<br>
//<br>
--------------------------------------------------------------------------------------------<br>
//<br>
<br>
renderer->ResetCamera();<br>
renWin->Render();<br>
renWinInt->Start();<br>
return 0;<br>
}<br></blockquote></div></div></div>