[vtkusers] vtkDICOMReader and vtkGDCMImageReader orientation
David Gobbi
david.gobbi at gmail.com
Fri Apr 1 16:32:03 EDT 2016
Hi Richard,
I'm only going to touch on one part of your email now... hopefully I'll be
able to
address some of the other points later.
You cannot generate the correct 3D surface from the image data alone. The
DICOM images are oriented images, which means that, if you want to know
where things are, you need not only the image data, but also "orientation"
and "position" information for the images. And of course you are already
familiar with this... so let's jump ahead to what you actually need to _do_
in
order to generate the correct surface.
For the vtkDICOMReader, it should be easy:
1) read the image (obviously!)
2) call the GetPatientMatrix() method to get position/orientation info
3) run marching cubes on the image to get a surface
4) correct the surface by using vtkTransformPolyDataFilter to apply the
aforementioned PatientMatrix to the surface
These same basic steps can be used with any of the DICOM readers,
but of course each reader is slightly different and the devil is in the
details.
For example, here are some tricky details for vtkDICOMReader (and it is
the least "tricky" of the three readers, IMHO):
1) If you look at the orientation returned by
vtkDICOMReader::GetPatientMatrix()
and compare it to the DICOM file metadata, you will notice that it is
different!
This is because the vtkDICOMReader (like other VTK readers) rasterizes its
output from bottom to top, to historically match the glDrawPixels function.
Since native DICOM is rasterized top-to-bottom, this requires a vertical
flip of
the orientation matrix.
For this reason, I strongly recommend that anyone using vtkDICOMReader in
a "serious" application call the SetMemoryRowOrderToFileNative() method,
which will cause the reader to rasterize the data from top-to-bottom and
will
cause the matrix returned by GetPatientMatrix() to exactly match the
original
DICOM meta data.
2) The vtkDICOMReader always uses IPP sorting for a volumetric series,
unless you specifically call SortingOff(), in which case it uses the order
of
the files. I do not recommend calling SortingOff(). Relying on IPP is
best.
3) For isosurfaces, note that RescaleIntercept/RescaleSlope can vary from
slice-to-slice in CT and PET. For such images, I wrote a filter called
vtkDICOMApplyRescale that can be applied to the data immediately after
the reader (but please read the doxygen docs for this filter before using
it).
4) An even more pernicious problem for CTs is that, if the scan is done with
a tilted gantry, the resulting series does not form a rectangular "volume".
For such images, the vtkDICOMCTRectifier can resample the data so that
it can be properly used within VTK.
Hopefully this info will be of some help. Let me know if you want
clarification on any of the above points.
Cheers,
- David
On Fri, Apr 1, 2016 at 6:52 AM, mbcx9rb9 <richard.j.brown at live.co.uk> wrote:
> 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;
> }
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/vtkusers/attachments/20160401/52605864/attachment.html>
More information about the vtkusers
mailing list