[vtkusers] vtkDICOMReader and vtkGDCMImageReader orientation
mbcx9rb9
richard.j.brown at live.co.uk
Fri Apr 1 08:52:26 EDT 2016
Hi all,
I'm having some troubles with orientation and location of my DICOMs when
loaded with vtkDICOMReader and vtkGDCMImageReader.
To illustrate my problems, I opened the same files with both the readers and
performed marching cubes on them to display 3D surfaces. I also displayed an
annotated orientation cube to allow me to compare these surfaces with
Slicer3D (which I'm assuming correctly calculates orientation etc.).
In the ideal world, the 3D surfaces should superpose perfectly, but this is
not the case.
My pipeline for reading DICOMs with vtkGDCMImageReader is:
gdcm::IPPSorter <-- Calculates the proper
z-spacing
vtkGDCMImageReader
vtkImageChangeInfromation <-- Applies the corrected z-spacing that
was calculated by IPPSorter
My pipeline for reading DICOMs with vtkDICOMReader is simply:
vtkDICOMReader
and then actor->SetUserMatrix(reader->GetPatientMatrix())
My problem is that the 3D structures are not in the same position nor the
same orientation. Furthermore, one might agree with Slicer3D for one image
but not for the next...
I guess my main question is - What is the correct way to systematically open
DICOM files (using either vtkDICOMReader or vtkGDCMImageReader, I'm not
fussy)?
Below are two images, one of a human skull and one of a pig skull (chosen at
random). Each time, the red is rendered from vtkGDCMImageReader, the black
is rendered from vtkDICOMReader (turns black when I use SetUserMatrix, not
sure why...) and the last is from Slicer3D.
With the human skull, we see that vtkGDCMImageReader agrees with slicer,
whereas vtkDICOMReader is upside-down (Superior-Inferior reflection).
However for the pig, vtkGDCMImageReader is upside-down (Superior-Inferior
reflection) (PLUS something funny has happened to the spacing), whereas
vtkDICOMReader is facing the anterior direction as opposed to the superior
direction.
<http://vtk.1045678.n5.nabble.com/file/n5737481/Screen_Shot_2016-04-01_at_12.png>
<http://vtk.1045678.n5.nabble.com/file/n5737481/Screen_Shot_2016-04-01_at_14.png>
Below is the code I have used for both. Any pointers on the correct way to
consistently read DICOM images correctly would be appreciated.
Regards,
Richard
P.S. A very important question - is my vtkAnnotatedCubeActor correctly
annotated? If not, the comparison between my 3D structures and Slicer3D is
useless until this is corrected...
#include <vtkSmartPointer.h>
#include <vtkDICOMReader.h>
#include <vtkDICOMFileSorter.h>
#include <vtkStringArray.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkInteractorStyleTrackballCamera.h>
#include <vtkMarchingCubes.h>
#include <vtkPolyDataNormals.h>
#include <vtkStripper.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkImageChangeInformation.h>
#include <vtkAnnotatedCubeActor.h>
#include <vtkOrientationMarkerWidget.h>
#include <gdcmIPPSorter.h>
#include <vtkGDCMImageReader.h>
#include <QStringList>
#include <QDirIterator>
int main()
{
//
--------------------------------------------------------------------------------------------
//
//
// Get filenames
//
//
--------------------------------------------------------------------------------------------
//
// Get all DCM files in directory
vtkSmartPointer<vtkStringArray> filenames =
vtkSmartPointer<vtkStringArray>::New();
QDirIterator iterator("/Scan/Directory/",
QStringList("*.dcm"), QDir::Files,
QDirIterator::Subdirectories);
while (iterator.hasNext()) {
filenames->InsertNextValue(iterator.next().toUtf8().data());
}
// Sort
vtkSmartPointer <vtkDICOMFileSorter> sorter = vtkSmartPointer
<vtkDICOMFileSorter>::New();
sorter->SetInputFileNames(filenames);
sorter->Update();
// Get most numerous series
int numSeries = sorter->GetNumberOfSeries();
int mostNumerousSeriesIndex = 0;
for (int i=0; i<numSeries; i++) {
if (sorter->GetFileNamesForSeries(i)->GetNumberOfValues() >
sorter->GetFileNamesForSeries(mostNumerousSeriesIndex)->GetNumberOfValues())
{
mostNumerousSeriesIndex = i;
}
}
vtkStringArray *seriesFilenames =
sorter->GetFileNamesForSeries(mostNumerousSeriesIndex);
//
--------------------------------------------------------------------------------------------
//
//
// Setup renderer etc.
//
//
--------------------------------------------------------------------------------------------
//
vtkSmartPointer<vtkRenderer> renderer =
vtkSmartPointer<vtkRenderer>::New();
renderer->SetGradientBackground(1);
renderer->SetBackground(.5,.5,.5);
renderer->SetBackground2(0,0,0);
vtkSmartPointer<vtkRenderWindow> renWin =
vtkSmartPointer<vtkRenderWindow>::New();
renWin->AddRenderer(renderer);
renWin->SetSize(1000,1000);
vtkSmartPointer<vtkInteractorStyleTrackballCamera> interactorStyle =
vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New();
vtkSmartPointer<vtkRenderWindowInteractor> renWinInt =
vtkSmartPointer<vtkRenderWindowInteractor>::New();
renWinInt->SetRenderWindow(renWin);
renWinInt->SetInteractorStyle(interactorStyle);
//
--------------------------------------------------------------------------------------------
//
//
// vtkGDCMImageReader
//
//
--------------------------------------------------------------------------------------------
//
// need to convert from vtkStringArray to std::vector<std::string> to be
able to sort...
std::vector<std::string> seriesFilenames_stdVector;
for (int i=0; i<seriesFilenames->GetNumberOfValues(); i++) {
seriesFilenames_stdVector.push_back(seriesFilenames->GetValue(i).c_str());
}
// sort and get z-spacing
gdcm::IPPSorter IPPs;
IPPs.SetComputeZSpacing(true);
IPPs.SetZSpacingTolerance(1e-3);
IPPs.Sort(seriesFilenames_stdVector);
// Convert the sorted std::vector back to vtkStringArray to be able to
use with vtkGDCMImageReader...
vtkSmartPointer<vtkStringArray> seriesFilenames_vtkStringArray_sorted =
vtkSmartPointer<vtkStringArray>::New();
for (size_t i=0; i<IPPs.GetFilenames().size(); i++) {
seriesFilenames_vtkStringArray_sorted->InsertNextValue(IPPs.GetFilenames().at(i));
}
// GDCM reader
vtkSmartPointer<vtkGDCMImageReader> gdcmReader =
vtkSmartPointer<vtkGDCMImageReader>::New();
gdcmReader->SetFileNames(seriesFilenames_vtkStringArray_sorted);
gdcmReader->FileLowerLeftOn();
// correct z-spacing
vtkSmartPointer<vtkImageChangeInformation> changeInfo =
vtkSmartPointer<vtkImageChangeInformation>::New();
changeInfo->SetInputConnection(gdcmReader->GetOutputPort());
changeInfo->SetOutputSpacing(gdcmReader->GetOutput()->GetSpacing()[0],
gdcmReader->GetOutput()->GetSpacing()[1], IPPs.GetZSpacing());
changeInfo->Update();
// Marching cubes
vtkSmartPointer<vtkMarchingCubes> gdcmMC =
vtkSmartPointer<vtkMarchingCubes>::New();
gdcmMC->SetInputConnection(changeInfo->GetOutputPort());
gdcmMC->SetNumberOfContours(1);
gdcmMC->GenerateValues(1,700.,3000.);
vtkSmartPointer<vtkPolyDataNormals> gdcmNormals =
vtkSmartPointer<vtkPolyDataNormals>::New();
gdcmNormals->SetInputConnection(gdcmMC->GetOutputPort());
gdcmNormals->SetFeatureAngle(60.0);
vtkSmartPointer<vtkStripper> gdcmStripper =
vtkSmartPointer<vtkStripper>::New();
gdcmStripper->SetInputConnection(gdcmNormals->GetOutputPort());
vtkSmartPointer<vtkPolyDataMapper> gdcmMapper =
vtkSmartPointer<vtkPolyDataMapper>::New();
gdcmMapper->SetInputConnection(gdcmStripper->GetOutputPort());
gdcmMapper->ScalarVisibilityOff();
gdcmMapper->Update();
vtkSmartPointer<vtkActor> gdcmActor = vtkSmartPointer<vtkActor>::New();
gdcmActor->SetMapper(gdcmMapper);
gdcmActor->GetProperty()->SetColor(1.0,0.0,0.0);
gdcmActor->GetProperty()->SetOpacity(1.);
renderer->AddActor(gdcmActor);
//
--------------------------------------------------------------------------------------------
//
//
// vtkDICOMReader
//
//
--------------------------------------------------------------------------------------------
//
// DCM reader
vtkSmartPointer<vtkDICOMReader> dcmReader =
vtkSmartPointer<vtkDICOMReader>::New();
dcmReader->SetFileNames(seriesFilenames);
dcmReader->Update();
// Marching cubes
vtkSmartPointer<vtkMarchingCubes> dcmMC =
vtkSmartPointer<vtkMarchingCubes>::New();
dcmMC->SetInputConnection(dcmReader->GetOutputPort());
dcmMC->SetNumberOfContours(1);
dcmMC->GenerateValues(1,700.,3000.);
vtkSmartPointer<vtkPolyDataNormals> dcmNormals =
vtkSmartPointer<vtkPolyDataNormals>::New();
dcmNormals->SetInputConnection(dcmMC->GetOutputPort());
dcmNormals->SetFeatureAngle(60.0);
vtkSmartPointer<vtkStripper> dcmStripper =
vtkSmartPointer<vtkStripper>::New();
dcmStripper->SetInputConnection(dcmNormals->GetOutputPort());
vtkSmartPointer<vtkPolyDataMapper> dcmMapper =
vtkSmartPointer<vtkPolyDataMapper>::New();
dcmMapper->SetInputConnection(dcmStripper->GetOutputPort());
dcmMapper->ScalarVisibilityOff();
dcmMapper->Update();
vtkSmartPointer<vtkActor> dcmActor = vtkSmartPointer<vtkActor>::New();
dcmActor->SetUserMatrix(dcmReader->GetPatientMatrix());
dcmActor->SetMapper(dcmMapper);
dcmActor->GetProperty()->SetColor(0.0,1.0,0.0);
dcmActor->GetProperty()->SetOpacity(1.);
renderer->AddActor(dcmActor);
//
--------------------------------------------------------------------------------------------
//
//
// Orientation cube
//
//
--------------------------------------------------------------------------------------------
//
vtkSmartPointer<vtkAnnotatedCubeActor> annotatedCubeActor =
vtkSmartPointer<vtkAnnotatedCubeActor>::New();
annotatedCubeActor->SetXPlusFaceText ("L (+X)");
annotatedCubeActor->SetXMinusFaceText("R (-X)");
annotatedCubeActor->SetYPlusFaceText ("P (+Y)");
annotatedCubeActor->SetYMinusFaceText("A (-Y)");
annotatedCubeActor->SetZPlusFaceText ("S (+Z)");
annotatedCubeActor->SetZMinusFaceText("I (-Z)");
annotatedCubeActor->SetZFaceTextRotation(90);
annotatedCubeActor->GetCubeProperty()->SetColor(0.5, 1.0, 1.0);
annotatedCubeActor->GetTextEdgesProperty()->SetLineWidth(1.0);
annotatedCubeActor->GetTextEdgesProperty()->SetDiffuse(0.0);
annotatedCubeActor->GetTextEdgesProperty()->SetAmbient(1.0);
annotatedCubeActor->GetTextEdgesProperty()->SetColor(0.18, 0.28, 0.23);
annotatedCubeActor->SetFaceTextScale(0.2);
annotatedCubeActor->GetXPlusFaceProperty()->SetColor(0.0,0.0,1.0);
annotatedCubeActor->GetXPlusFaceProperty()->SetInterpolationToFlat();
annotatedCubeActor->GetXMinusFaceProperty()->SetColor(0.0,0.0,1.0);
annotatedCubeActor->GetXMinusFaceProperty()->SetInterpolationToFlat();
annotatedCubeActor->GetYPlusFaceProperty()->SetColor(0.0,1.0,0.0);
annotatedCubeActor->GetYPlusFaceProperty()->SetInterpolationToFlat();
annotatedCubeActor->GetYMinusFaceProperty()->SetColor(0.0,1.0,0.0);
annotatedCubeActor->GetYMinusFaceProperty()->SetInterpolationToFlat();
annotatedCubeActor->GetZPlusFaceProperty()->SetColor(1.0,0.0,0.0);
annotatedCubeActor->GetZPlusFaceProperty()->SetInterpolationToFlat();
annotatedCubeActor->GetZMinusFaceProperty()->SetColor(1.0,0.0,0.0);
annotatedCubeActor->GetZMinusFaceProperty()->SetInterpolationToFlat();
vtkSmartPointer<vtkOrientationMarkerWidget> widget =
vtkSmartPointer<vtkOrientationMarkerWidget>::New();
widget->SetOutlineColor( 0.9300, 0.5700, 0.1300 );
widget->SetOrientationMarker( annotatedCubeActor );
widget->SetInteractor( renWinInt );
widget->SetViewport( 0.0, 0.0, 0.2, 0.2 );
widget->SetEnabled( 1 );
widget->InteractiveOn();
//
--------------------------------------------------------------------------------------------
//
//
// Render
//
//
--------------------------------------------------------------------------------------------
//
renderer->ResetCamera();
renWin->Render();
renWinInt->Start();
return 0;
}
--
View this message in context: http://vtk.1045678.n5.nabble.com/vtkDICOMReader-and-vtkGDCMImageReader-orientation-tp5737481.html
Sent from the VTK - Users mailing list archive at Nabble.com.
More information about the vtkusers
mailing list