Generating and reading .vtk files from .dcm files, and modifying model

Question

I would like to have a better understanding about dicom volume rendering.

I have a set of dicom images, from which I've been able to extract axial, coronal and sagittal cuts, as follow :

I first wanted to generate a 3D model from scratch, but it seems to be way too hard.

So I heard about VTK/ITK, and I've been using this code to generate a .vtk file from my set of images :

http://www.itk.org/Doxygen46/html/IO_2DicomSeriesReadImageWrite2_8cxx-example.html

It works, but I need some explanations :

When I open this file with ParaView, I get the following result :

First, it might be a stupid question, but why is it blue?

Is there a way to cut and see the inside of the model?

My aim is not to use ParaView, and I'd like to make my own .vtk reader, I found this code I don't remember where, which I think is suppose to work, but all I get with it is the green background with nothing more :

#include <vtkPolyDataReader.h>
#include <vtkSmartPointer.h>
#include <vtkPolyDataMapper.h>
#include <vtkActor.h>
#include <vtkRenderWindow.h>
#include <vtkRenderer.h>
#include <vtkRenderWindowInteractor.h>

int main ( int argc, char *argv[] ) {

  // Parse command line arguments                                                                     
  if (argc != 2) {
    std::cerr << "Usage: " << argv[0] << " Filename(.vtk)" << std::endl;
    return EXIT_FAILURE;
  }

  std::string filename = argv[1];

  // Read all the data from the file                                                                  
  vtkSmartPointer<vtkPolyDataReader> reader = vtkSmartPointer<vtkPolyDataReader>::New();
  reader->SetFileName(filename.c_str());
  reader->Update();

  // Visualize                                                                                        
  vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
  mapper->SetInputConnection(reader->GetOutputPort());

  vtkSmartPointer<vtkActor> actor = vtkSmartPointer<vtkActor>::New();
  actor->SetMapper(mapper);

  vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer<vtkRenderer>::New();
  vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
  renderWindow->AddRenderer(renderer);

  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
    vtkSmartPointer<vtkRenderWindowInteractor>::New();
  renderWindowInteractor->SetRenderWindow(renderWindow);

  renderer->AddActor(actor);
  renderer->SetBackground(.3, .6, .3); // Background color green                                      

  renderWindow->Render();
  renderWindowInteractor->Start();

  return EXIT_SUCCESS;
}

any idea why? I saw in ParaView that I had to activate the "Volume" mode to be able to see my model, is there something similar to handle here?

Last thing, which is very important : Is it possible to modify the 3D volume inside a .vtk file? For instance, if I want to change the color of a specific part of the model, does VTK provides tools allowing that?

Answer 1

A lot of questions here ! Here are some answers.

• The rendering is blue because the lookup table is the default one (from blue to red) in paraview. You can edit it using the color map editor
• There are indeed some way to "cut" inside the data, please look at the example posted after.
• The example doesn't work because you are trying to load a vtkImageData using a polydata (mesh) reader
• That is exactly why you have to select "Volume" in Paraview, it's because your data is a volume (3D array of voxels). You can do the same using any of the volume mappers available in VTK.
• Yes you can edit the value of the volume, but that would lead us a little to far for now, we'll come to that later if needed ;)

Here is a full example that reads all the DICOM files under a directory, build a volume, render it using a volume render and enables a box widget to clip the volume interactively.

#include "vtkBoxRepresentation.h"
#include "vtkBoxWidget2.h"
#include "vtkCamera.h"
#include "vtkColorTransferFunction.h"
#include "vtkCommand.h"
#include "vtkDICOMImageReader.h"
#include "vtkGPUVolumeRaycastMapper.h"
#include "vtkImageData.h"
#include "vtkInteractorStyle.h"
#include "vtkInteractorStyleTrackballCamera.h"
#include "vtkMath.h"
#include "vtkPiecewiseFunction.h"
#include "vtkPlanes.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkRenderer.h"
#include "vtkVolume.h"
#include "vtkVolumeProperty.h"

// Box interaction callback
class vtkBoxCallback : public vtkCommand
{
public:
    static vtkBoxCallback *New(){ return new vtkBoxCallback; }
    vtkGPUVolumeRayCastMapper* m_mapper;
    vtkPlanes* m_planes;

    virtual void Execute( vtkObject* a_caller, unsigned long, void* ){
        vtkBoxWidget2* l_box_wdget = vtkBoxWidget2::SafeDownCast( a_caller );
        ( (vtkBoxRepresentation*)l_box_wdget->GetRepresentation() )->GetPlanes( m_planes );
        this->m_mapper->SetClippingPlanes( m_planes );
    }

    vtkBoxCallback(){}
};



int main( int argc, char *argv[] ){

    // Read volume
    vtkDICOMImageReader* l_reader = vtkDICOMImageReader::New();
    l_reader->SetDirectoryName( "C:/PathToDicomFiles/" );
    l_reader->Update();

    // Setup rendering stuff
    vtkRenderer* l_renderer = vtkRenderer::New();
    l_renderer->SetBackground( 0.3, 0.3, 0.3 );

    vtkRenderWindow* l_render_windows = vtkRenderWindow::New();
    l_render_windows->AddRenderer( l_renderer );
    l_render_windows->SetSize( 900, 900 );

    vtkInteractorStyleTrackballCamera* l_trackball = vtkInteractorStyleTrackballCamera::New();

    vtkRenderWindowInteractor* l_iren = vtkRenderWindowInteractor::New();
    l_iren->SetInteractorStyle( l_trackball );
    l_iren->SetRenderWindow( l_render_windows );
    l_iren->GetInteractorStyle()->SetDefaultRenderer( l_renderer );
    l_iren->SetDesiredUpdateRate( 15 );

    // Make sure we have an opengl context
    l_render_windows->Render();


    // Setup GPU volume raycast mapper
    vtkGPUVolumeRayCastMapper* l_gpu_mapper = vtkGPUVolumeRayCastMapper::New();
    l_gpu_mapper->SetInputConnection( l_reader->GetOutputPort() );

    // Setup Volume property
    // Window/Level
    double wl = 260;
    double ww = 270;

    // Color function
    vtkColorTransferFunction* l_color = vtkColorTransferFunction::New();
    l_color->SetColorSpaceToRGB();
    l_color->AddRGBPoint( wl - ww / 2, 0, 0, 0 );
    l_color->AddRGBPoint( wl - ww / 2 + 94 * ( ww / 255.0 ), 1., 21. / 255.0, 27. / 255.0 );
    l_color->AddRGBPoint( wl - ww / 2 + 147 * ( ww / 255.0 ), 1., 176. / 255.0, 9. / 255.0 );
    l_color->AddRGBPoint( wl - ww / 2 + 201 * ( ww / 255.0 ), 1., 241. / 255.0, 39. / 255.0 );
    l_color->AddRGBPoint( wl - ww / 2 + 255 * ( ww / 255.0 ), 1, 1, 1. );
    l_color->Build();

    // Opacity function
    vtkPiecewiseFunction* l_opacity = vtkPiecewiseFunction::New();
    l_opacity->AddPoint( wl - ww / 2, 0 );
    l_opacity->AddPoint( wl + ww / 2, 1 );

    // Volume property, light, shading
    vtkVolumeProperty* l_volume_property = vtkVolumeProperty::New();
    l_volume_property->SetColor( l_color );
    l_volume_property->SetScalarOpacity( l_opacity );
    l_volume_property->SetInterpolationTypeToLinear();
    l_volume_property->ShadeOn();
    l_volume_property->SetAmbient( 0.15 );
    l_volume_property->SetDiffuse( 0.8 );
    l_volume_property->SetSpecular( 0.25 );
    l_volume_property->SetSpecularPower( 40 );

    // Put everything together
    vtkVolume* l_volume = vtkVolume::New();
    l_volume->SetProperty( l_volume_property );
    l_volume->SetMapper( l_gpu_mapper );
    l_renderer->AddVolume( l_volume );
    l_renderer->ResetCamera();

    // setup Box interactive widget
    vtkBoxRepresentation* l_box_rep = vtkBoxRepresentation::New();
    l_box_rep->SetInsideOut( true );

    vtkBoxWidget2* l_voi_widget = vtkBoxWidget2::New();
    l_voi_widget->SetRepresentation( l_box_rep );
    l_voi_widget->SetInteractor( l_iren );
    l_voi_widget->GetRepresentation()->SetPlaceFactor( 1. );
    l_voi_widget->GetRepresentation()->PlaceWidget( l_reader->GetOutput()->GetBounds() );
    l_voi_widget->SetEnabled( true );

    vtkPlanes* l_planes = vtkPlanes::New();

    vtkBoxCallback* l_callback = vtkBoxCallback::New();
    l_callback->m_mapper = l_gpu_mapper;
    l_callback->m_planes = l_planes;
    l_voi_widget->AddObserver( vtkCommand::InteractionEvent, l_callback );

    // Go rendering !
    l_iren->Start();

    // Memory cleanup
    l_reader->Delete();
    l_renderer->Delete();
    l_render_windows->Delete();
    l_trackball->Delete();
    l_iren->Delete();
    l_gpu_mapper->Delete();
    l_color->Delete();
    l_opacity->Delete();
    l_volume_property->Delete();
    l_volume->Delete();
    l_voi_widget->Delete();
    l_planes->Delete();
    l_callback->Delete();
}

As a general advice, I suggest you read VTK examples that should helps you to understand all VTK capabilities.

Hope that helps :)