[vtkusers] vtkDICOMReader and vtkGDCMImageReader orientation
Richard Brown
richard.j.brown at live.co.uk
Mon Apr 4 06:09:55 EDT 2016
Dear David,
Thanks for the help. I applied the suggested changes to vtkDICOMReader and I seem to be able to consistently produce 3D surfaces that agree with Slicer, which is always reassuring.
> 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.
I read that scaling, reflection/rotation, shear and translation are included in the 4x4 matrix. Given that shear is included in this matrix, shouldn’t gantry tilt be accounted for when I use GetPatientMatrix()?
And then moving on from 3D surfaces…
If I didn’t want to produce 3D surfaces; if I wanted to manipulate vtkImageData, is there any precautions I should take?
For example, if I want to superpose DICOM RT-Structures or other surfaces (notably .STL in my case) onto my vtkImageData, do I need to apply any sort of matrix to either the DICOM vtkImageData or the superposed structures to make sure that they superpose in the correct position with the correct orientation?
And lastly, to get the world coordinates (x,y,z) from my vtkImageData coordinates (i,j,k), do I simply apply the GetPatientMatrix()? Is there a method for this? And the same question applies for world coordinate —> index coordinate conversions - I assume I calculate the inverse matrix and apply that in a similar fashion?
Many thanks,
Richard
> On 01 Apr 2016, at 22:32, David Gobbi <david.gobbi at gmail.com> wrote:
>
> 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;
>> }
>>
More information about the vtkusers
mailing list