How to improve vtkDICOMImageReader output

Michele · March 7, 2019, 3:04pm

Hi everyone, I am VTK user and recently I needed to start visualizing some DICOM in the 3D space. The idea is to shows some stl from segmentation together with some relative DICOM.

Essentially, I follow the vtkDICOMImageReader example (https://vtk.org/Wiki/VTK/Examples/Cxx/IO/ReadDICOM) which works fine except for the quality of the images I get. Attached you may find an example: on the left the result I get with VTK, on the right the same DICOM read with MicroDicom.

Browsing the net I saw that the expected outcome should be much clearer, however I couldn’t figure it out how to adjust brightness and contrast.

Any hint?

Thank you in advance,

Michele

here is my code

    // Read the DICOM 
vtkSmartPointer<vtkDICOMImageReader> reader =    
  vtkSmartPointer<vtkDICOMImageReader>::New();

reader->SetFileName("myDICOM");
reader->Update();

// Mapper
vtkSmartPointer<vtkImageSliceMapper > imageMapper =
  vtkSmartPointer<vtkImageSliceMapper>::New();

imageMapper->SetInput(reader->GetOutput());


//VISUALIZATION
vtkImageActor *Actor	= vtkImageActor::New();

vtkRenderer			*render			= vtkRenderer::New();
vtkRenderWindow		*renWin			= vtkRenderWindow::New();
vtkRenderWindowInteractor *iren		= vtkRenderWindowInteractor::New();

Actor->SetMapper( imageMapper );
Actor->InterpolateOn();
render->AddActor( Actor );

render->SetBackground(255,255,255);
renWin->AddRenderer( render );
iren->SetRenderWindow(renWin);

iren->Initialize();
iren->Start();

banesullivan · March 8, 2019, 6:26am

Can you share that data file?

Does the sample dataset in the example you shared also look different?

Michele · March 8, 2019, 9:09am

Hi Bane,

thank you for your prompt reply.
You may find one of my DICOM at the following link (sorry, as new user I am not allowed to attach files).

I tried the dataset of the example, and it works fine. However these data are in .img format. Do you suggest I should convert my DICOM?

Michele

banesullivan · March 9, 2019, 6:16am

Here’s a way to render the image with a bit more clarity. Wrote the code in Python and used matplotlib to get a nice colormap but this could be translated to C++:

import vtk
from vtk.util import numpy_support
import numpy as np
from matplotlib.cm import get_cmap

reader = vtk.vtkDICOMImageReader()

reader.SetFileName('DICOM_test.dcm')
reader.Update()

mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(reader.GetOutput())
mapper.SetArrayName('DICOMImage')


actor = vtk.vtkActor()
actor.SetMapper( mapper );

mapper.SetScalarModeToUsePointData()
mapper.GetLookupTable().SetNumberOfTableValues(256)
mapper.InterpolateScalarsBeforeMappingOn()
mapper.SetScalarRange(0, 6930)


cmap = get_cmap('bone')
ctable = cmap(np.linspace(0, 1, 256))*255
ctable = ctable.astype(np.uint8)
mapper.GetLookupTable().SetTable(numpy_support.numpy_to_vtk(ctable))

actor.GetProperty().SetRepresentationToSurface()
actor.GetProperty().LightingOff()

render = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
iren = vtk.vtkRenderWindowInteractor()

render.AddActor( actor )

render.SetBackground(255,255,255)
renWin.AddRenderer( render )
iren.SetRenderWindow(renWin)
istyle = vtk.vtkInteractorStyleTrackballCamera()
iren.SetInteractorStyle(istyle)


iren.Initialize()
iren.Start()

Note that rendering images like this can be quite cumbersome and if you’re new to VTK, I’d recommend using Python for rapid prototyping. Also, check out vtki. I’ll add a DICOM reader to vtki soon but you could get a similar rendering with a bit more intuitve controls using vtki:

import vtk
import vtki

filename = 'DICOM_test.dcm'

# read using vtk
reader = vtk.vtkDICOMImageReader()
reader.SetFileName(filename)
reader.Update()
data = vtki.wrap(reader.GetOutput())
# OR when I add the DICOM reader to `vtki`:
# data = vtki.read(filename)

# and plot it with a bone colormap
data.plot(cmap='bone', cpos='xy')

banesullivan · March 9, 2019, 8:01pm

Also is it okay to distribute this data as an example for reading DICOM images?

Michele · March 10, 2019, 9:09am

Hi Bane,

thank you for your help. First of all, there is no problem to use this data for an example.

Coming to your suggestion, if I understood correctly you built a color map outside VTK and than pass it to the vtkDataSetMapper to forcing the right visualizaiton. I have not yet tested your solution since I am coding in c++ and at the moment I have not the time to switch to pyton. However I should be able to do the same using Matlab and C++. I will keep you updated on the results.

However, what it is not clear to me is why the DICOMReader does behave this way. I had the same results with all the DICOM I tested, coming from very different MR scanners. So it seems that the DICOMReader actually has problem in reading standard DICOM format. Do you know if DICOMReader needs a specified data format? Has said, data from the example are in .img.

Thank you once again,

Michele

banesullivan · March 10, 2019, 11:38pm

Yes, you should be able to copy most of that code and switch a few syntax things. You could also try simply number of colors and scalar range instead of making a colormap (I do this below).

I’m not sure there is anything wrong with the DICOM reader as the data values being read in seem correct, it’s just a bit tricky to get the rendering/color mapping to properly display those values (as is the case with any data). To check, when I read in that DICOM file, it sees a single vtkShortArray of data values (i.e. int16). These values range at (0, 6930) for the specific file you shared with me. Are those results not expected?

I think the reader is behaving correctly, you’re just experiencing trouble with your mapper. Perhaps this code works best for you to avoid building a colormap (this is pure VTK so it can be translated to C++ easily):

import vtk

reader = vtk.vtkDICOMImageReader()
reader.SetFileName('DICOM_test.dcm')
reader.Update()

dataset = reader.GetOutput()
array = dataset.GetPointData().GetScalars()
mini, maxi = array.GetFiniteRange()

mapper = vtk.vtkDataSetMapper()
mapper.SetInputData(reader.GetOutput())
table = mapper.GetLookupTable()
table.SetNumberOfTableValues(256)
mapper.SetScalarRange(mini, maxi)

actor = vtk.vtkActor()
actor.SetMapper(mapper);

actor.GetProperty().SetRepresentationToSurface()
actor.GetProperty().LightingOff()

render = vtk.vtkRenderer()
renWin = vtk.vtkRenderWindow()
iren = vtk.vtkRenderWindowInteractor()

render.AddActor(actor)

render.SetBackground(255, 255, 255)
renWin.AddRenderer( render )
iren.SetRenderWindow(renWin)
istyle = vtk.vtkInteractorStyleTrackballCamera()
iren.SetInteractorStyle(istyle)

iren.Initialize()
iren.Start()

Kenichiro-Yoshimi · March 11, 2019, 3:12am

Hi Michele,

Here is another way to improve the output of the first vtk example. This is a rough method to additionally set the color level and window. For now, I don’t know if it is the usual way.

#include <vtkSmartPointer.h>
#include <vtkDataArray.h>
#include <vtkImageData.h>
#include <vtkImageViewer2.h>
#include <vtkDICOMImageReader.h>
#include <vtkPointData.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderer.h>

int main(int argc, char* argv[])
{
  // Verify input arguments
  if ( argc != 2 )
  {
    std::cout << "Usage: " << argv[0]
              << " Filename(.img)" << std::endl;
    return EXIT_FAILURE;
  }

  std::string inputFilename = argv[1];

  // Read all the DICOM files in the specified directory.
  auto reader = vtkSmartPointer<vtkDICOMImageReader>::New();
  reader->SetFileName(inputFilename.c_str());
  reader->Update();

  auto image = reader->GetOutput();
  double *range = image->GetPointData()->GetScalars()->GetRange();

  // Visualize
  auto imageViewer = vtkSmartPointer<vtkImageViewer2>::New();
  imageViewer->SetInputConnection(reader->GetOutputPort());
  auto renderWindowInteractor =
    vtkSmartPointer<vtkRenderWindowInteractor>::New();
  imageViewer->SetColorLevel(0.5 * (range[0] + range[1]));
  imageViewer->SetColorWindow(range[1] - range[0]);
  imageViewer->SetupInteractor(renderWindowInteractor);
  imageViewer->GetRenderer()->ResetCamera();
  imageViewer->Render();

  renderWindowInter
  actor->Start();

  return EXIT_SUCCESS;
}

Michele · March 11, 2019, 9:34am

Thank you all for the help and support. Actually, the solution proposed by Kenichiro Yoshimi works fine with the minimum amount of additional code so I will go for that.

Best

Michele