[vtkusers] Correct anatomical orientation of volumes from Niftii and DICOM

[ochampao]

Hi vtkUsers,

I am developing an application for displaying medical images (four pane
viewer with Sagittal, Coronal and Axial views). The volumes are loaded from
different files that can be DICOM, Niftii or HDF5. Multiple volumes can be
loaded and displayed simultaneously (by overlaying 2D slices). 

Until recently I have been using vtkDICOMImageReader for DICOM (see
loadDicom()) and vtkNIFTIImageReader for Niftii (see loadNiftii()). Using
these two functions I was able to display the 3 views with the correct
anatomical orientations (i.e. Sagittal, Coronal and Axial). After reading
the very nice guide by David Gobbi about Medical Image Orientation (see
http://calgaryimageanalysis.ca/wiki/images/5/52/Image-orientation.pdf )
I followed the recommendation therein and decided to switch to
vtkGDCMImageReader (see loadDicomGDCM() ). Despite using FileLowerLeftOn()
and SetOutputOrigin() (as recommended in the guide) after switching to the
vtkGDCMImageReader I am unable to display datasets coming from both DICOM
and Niftii with the correct anatomical orientation. 

Using the pipeline as configured in the code segments below, I am able to
correctly display volumes loaded from Niftii files (using loadNiftii()) and
incorrectly display volumes loaded from DICOM (using loadDicomGDCM()). The
2D slices loaded using loadDicomGDCM() appear flipped both in the X and Y
axes. The X axis flip can be avoided by skipping the call to
FileLowerLeftOn().

So, my question is: how can I setup/correct my pipeline below to correctly
display volumes loaded from both Niftii and DICOM (using
vtkGDCMImageReader)? 

Below you can see relevant segments of my code.

Thanks for your help
Panos.

================================================================

// setups empty view with a vtkImageStack for overalying slices from
multiple datasets
void QuadView::initStaticPipelineForView(const IQuadViewUI::VTKDockWidget
view)
{
    // create an object that will hold the pointers to the static objects of
the pipeline
    PipelineObjects2D pipelineObject2D;

    // setup image stack for overalying the multiple slices
    pipelineObject2D.imageStack = vtkSmartPointer<vtkImageStack>::New();
    
    // Setup renderers
    pipelineObject2D.renderer = vtkSmartPointer<vtkRenderer>::New();
    pipelineObject2D.renderer->AddViewProp(pipelineObject2D.imageStack);
    pipelineObject2D.renderer->GetActiveCamera()->ParallelProjectionOn();
    pipelineObject2D.renderer->ResetCameraClippingRange();
    pipelineObject2D.renderer->ResetCamera();

    // Setup camera for each view
    if (view == IQuadViewUI::VTKDockWidget::TopLeft) // Axial
    {
        // do nothing
    }
    
    if (view == IQuadViewUI::VTKDockWidget::TopRight) // Sagittal
    {
        pipelineObject2D.renderer->GetActiveCamera()->Azimuth(90); 
        pipelineObject2D.renderer->GetActiveCamera()->Roll(90); 
        pipelineObject2D.renderer->GetActiveCamera()->OrthogonalizeViewUp();
    }
    
    if (view == IQuadViewUI::VTKDockWidget::BottomLeft) // Coronal
    {
        pipelineObject2D.renderer->GetActiveCamera()->Elevation(90);
        pipelineObject2D.renderer->GetActiveCamera()->OrthogonalizeViewUp();
    }

    // Setup render window
    pipelineObject2D.renderWindow =
vtkSmartPointer<vtkGenericOpenGLRenderWindow>::New();
    pipelineObject2D.renderWindow->AddRenderer(pipelineObject2D.renderer);

    // Attach render window to QVTK widget
    QVTKOpenGLWidget& qVtkWidget = mQuadViewUI.getQVTKOpenGLWidget(view); 
    qVtkWidget.SetRenderWindow(pipelineObject2D.renderWindow);

    //setup interaction style
    pipelineObject2D.interactorStyle =
vtkSmartPointer<vtkInteractorStyleImage>::New();
    pipelineObject2D.interactorStyle->SetInteractionModeToImageSlicing();

    // Setup render window interactor
    pipelineObject2D.renderWindowInteractor =
vtkSmartPointer<vtkRenderWindowInteractor>::New();
   
pipelineObject2D.renderWindowInteractor->SetInteractorStyle(pipelineObject2D.interactorStyle);
   
pipelineObject2D.renderWindowInteractor->SetRenderWindow(pipelineObject2D.renderWindow);
    pipelineObject2D.renderWindowInteractor->Initialize();

    mPipelineObject2DContainer.insert(std::pair<IQuadViewUI::VTKDockWidget,
PipelineObjects2D>(view, pipelineObject2D));

    // render
    pipelineObject2D.renderWindow->Render();
}

// load a new volume
unsigned long QuadView::addImage(const char* filename)
{
    ImageDataStruct imageDataStruct;

    // *** SELECT ONE OF READERS
    imageDataStruct.imageData = loadDicomGDCM(filename);
    //imageDataStruct.imageData = loadDicomGDCM(filename);
    //imageDataStruct.imageData = loadNiftii(filename);
    //imageDataStruct.imageData = loadH5(filename);

    // setup lookup for the image data. 
    imageDataStruct.lookupTable = setupLookupTableInternal();

    // iterate through the active 2D views
    for (const auto &view : mQuadViewUI.getQVTKOpenGLWidget2DKeyList())
    {
        // get a reference to the 
        PipelineObjects2D &pipelineObjects2D =
mPipelineObject2DContainer.at(view);

        // setup slice mapper for the current view
        imageDataStruct.resliceMapperContainer[view] =
vtkSmartPointer<vtkImageResliceMapper>::New();
       
imageDataStruct.resliceMapperContainer[view]->SetInputData(imageDataStruct.imageData);
        imageDataStruct.resliceMapperContainer[view]->SliceFacesCameraOn();
        imageDataStruct.resliceMapperContainer[view]->SliceAtFocalPointOn();

        // setup prop holding the slice
        imageDataStruct.imageSliceContainer[view] =
vtkSmartPointer<vtkImageSlice>::New();
       
imageDataStruct.imageSliceContainer[view]->SetMapper(imageDataStruct.resliceMapperContainer[view]);
       
imageDataStruct.imageSliceContainer[view]->GetProperty()->SetLookupTable(imageDataStruct.lookupTable);
       
imageDataStruct.imageSliceContainer[view]->GetProperty()->UseLookupTableScalarRangeOn(); 

        // add to image stack
       
pipelineObjects2D.imageStack->AddImage(imageDataStruct.imageSliceContainer[view]);
    }

    // reset camera and clipping range
    resetCameras2D();

    // render views
    renderViews2D();
    
    // save object holding the pointers to the dynamic part of pipeline
associated with the current image data
    mImageDataContainer.insert(
        std::pair<unsigned long, ImageDataStruct>(imageDataStruct.getUID(),
imageDataStruct));

    return imageDataStruct.getUID();
}

// load volume from dicom using vtkDICOMImageReader
vtkSmartPointer<vtkImageData> QuadView::loadDicom(const char* filename)
const
{
    vtkNew<vtkDICOMImageReader> dicomReader;
    dicomReader->SetDirectoryName(filename);
    dicomReader->Update();

    return dicomReader->GetOutput();
}

// load volume from dicom using vtkDICOMImageReader
vtkSmartPointer<vtkGDCMImageReader> QuadView::loadDicomGDCM(const char*
finename) const
{
    gdcm::Directory gdcmDir;
    gdcmDir.Load(finename, false);

    gdcm::IPPSorter ippSorter;
    ippSorter.SetComputeZSpacing(true);
    ippSorter.SetZSpacingTolerance(1e-3);
    ippSorter.Sort(gdcmDir.GetFilenames());

    vtkNew<vtkStringArray> sortedFiles;
    for (const auto file : ippSorter.GetFilenames())
        sortedFiles->InsertNextValue(file);

    vtkNew<vtkGDCMImageReader> reader;
    reader->SetFileNames(sortedFiles);
    reader->FileLowerLeftOn();
    reader->Update();

    vtkNew<vtkImageChangeInformation> imageInfo;
    imageInfo->SetOutputOrigin(0.0, 0.0, 0.0);
    imageInfo->SetOutputSpacing(
        reader->GetOutput()->GetSpacing()[0], 
        reader->GetOutput()->GetSpacing()[1], 
        ippSorter.GetZSpacing());
    imageInfo->SetInputConnection(reader->GetOutputPort());
    imageInfo->Update();

    return imageInfo->GetOutput();
}

// load volume from Niftii
vtkSmartPointer<vtkImageData> QuadView::loadNiftii(const char* filename)
const
{
    vtkNew<vtkNIFTIImageReader> niftiiReader;
    niftiiReader->SetFileName(filename);
    niftiiReader->Update();

    return niftiiReader->GetOutput();
}

void QuadView::renderViews2D() const
{
    // Render changes in 2D views
    for (const auto view : mQuadViewUI.getQVTKOpenGLWidget2DKeyList())
        mPipelineObject2DContainer.at(view).renderWindow->Render();
}

void QuadView::resetCameras2D() const
{
    // iterate through the active 2D views.
    for (const auto &view : mQuadViewUI.getQVTKOpenGLWidget2DKeyList())
    {
        // Reset camera & clipping range
        mPipelineObject2DContainer.at(view).renderer->ResetCamera();

        switch (view)
        {
            case IQuadViewUI::VTKDockWidget::TopLeft: //Axial
               
mPipelineObject2DContainer.at(view).renderer->GetActiveCamera()->SetViewUp(0,
1, 0);
               
mPipelineObject2DContainer.at(view).renderer->GetActiveCamera()->OrthogonalizeViewUp();
                break;
            case IQuadViewUI::VTKDockWidget::TopRight: //Sagittal
               
mPipelineObject2DContainer.at(view).renderer->GetActiveCamera()->SetViewUp(0,
0, -1);
               
mPipelineObject2DContainer.at(view).renderer->GetActiveCamera()->OrthogonalizeViewUp();
                break;
            case IQuadViewUI::VTKDockWidget::BottomLeft: //Coronal
               
mPipelineObject2DContainer.at(view).renderer->GetActiveCamera()->SetViewUp(0,
0, -1);
               
mPipelineObject2DContainer.at(view).renderer->GetActiveCamera()->OrthogonalizeViewUp();
                break;
            default:
                break;
        }
    }
}
================================================================




--
Sent from: http://vtk.1045678.n5.nabble.com/VTK-Users-f1224199.html