Some tips&tricks for 3D tumor visualization

There are several reasons why we rarely see 3D agent based models of tumor growth, with I think the two most important ones being: 1) 3D simulations are way more computationally expensive than 2D simulations; 2) it’s harder to visualize results of 3D simulations to present all necessary information. The other frequently used argument is that there is no need for 3D, because the crucial phenomena can be observed as well in 2D. However, when one considers e.g. stochastic models there might be a huge difference in the system behavior between 2D and 3D case (see very nice post Why is that a 2D random walk is recurrent, while a 3D random walk is transient?).

Possible improvements in the computational speed are highly dependent on the given model setting and typically require sophisticated techniques. Problem with nice visualizations, however, just requires proper tools. Today, I will show my attempts to make a good-looking 3D tumor visualizations, which I have taken over the last several years.

Several years ago when working with @heikman on the project about the impact of cellular senescence  I was using The Persistence of Vision Raytracer (POV-Ray) in order to obtain high quality 3D tumor visualization. After simulating the tumor up to a given cell number I was parsing the information about cells locations and their basic properties into a POV-Ray specific format. In my approach, each cell was separate actor in the scene with possible different surface/texture/color properties. As it can be expected renderings with raytracing software are very good (see Fig.1 below), but the whole procedure to prepare input files to POV-Ray, to set-up properly camera and light was really time-consuming. In addition, rendering takes a long time, so there is no chance for “live” rendering during the simulation, which would allow to observe simulation behavior.

Figure 1. Exemplary 3D tumor renderings obtained using POV-Ray software. Bottom-left picture shows tumor consisted of 9.4 million of cells.

I really wanted to have a tool that would allow me to seamlessly simulate and visualize 3D tumor at the same time. I quickly realized that rendering each cell as a separate object is not a good approach, when you simulate more than 1 million cells. Hence, I started looking around for an algorithm that would allow to extract tumor surface from the cloud of cells. I quickly found well known Marching Cubes algorithm  that does the trick very well.  I’ve implemented the algorithm and visualized resulting surface using Visualization Library in C++. I was quite happy with the results (see Fig. 2 below), however, at some point code became so complicated and so model specific that I assumed that it doesn’t have any future.

Untitled-2

Figure 2. “Live” 3D tumor visualization during the simulation. Code developed in C++ using Visualization Library.

Quite recently I started using VTK library from within C++ (you can also use it form within Python) and I quickly realized that it might be the ultimate tool for 3D visualization. Now I will show how easy it is to produce great visualizations that can be used in “live” renderings during 3D simulations using VTK library.

In the following example I will use glyph mapper that basically takes my cloud of points (cells), creates cubes around them, and extracts the surface. First, let us include all necessary headers and define necessary variables.

#include "vtkCubeSource.h"
#include "vtkPolyDataMapper.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkInteractorStyleTrackballCamera.h"
#include "vtkPolyData.h"
#include "vtkGlyph3DMapper.h"

vtkSmartPointer<vtkRenderer> renderer;
vtkSmartPointer<vtkRenderWindow> renderWindow;
vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor;
vtkSmartPointer<vtkPoints> points;

Now we can write the function that will initialize the renderer and rendering window together with its interactor object.

void initializeRendering() {
    renderer = vtkSmartPointer<vtkRenderer>::New(); //creating new renderer
    renderer->SetBackground(1., 1., 1.); //setting background color to white

    renderWindow = vtkSmartPointer<vtkRenderWindow>::New(); //creating renderer window
    renderWindow->SetSize(1000, 1000); //setting size of the window
    renderWindow->AddRenderer(renderer); //adding previously created renderer

    renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New(); //creating render
    renderWindowInteractor->SetRenderWindow(renderWindow); //adding previously created renderer window

    points = vtkSmartPointer<vtkPoints>::New();         //vector with points
}

I assume that we have our cell objects are stored in STL vector structure and that each cell contains information about its (x,y,z) position in 3D space. Having that we can write a really short function that adds our cells to the scene.

void addCells(vector<Cell>::iterator begin, vector<Cell>::iterator end) {
    vector<Cell>::iterator cell = begin;
    for (cell; cell<end; cell++)
points->InsertNextPoint(cell->x,cell->y,cell->z);
}

Finally, we can write the rendering function that is invoked in the main part of the code.

void visualizeLesion(vector<Cell>::iterator begin, vector<Cell>::iterator end) {

    initializeRendering();

    addCells(begin, end);

    // Combine into a polydata
    vtkSmartPointer<vtkPolyData> polydata = vtkSmartPointer<vtkPolyData>::New();
    polydata->SetPoints(points);

    vtkSmartPointer<vtkCubeSource> cubeSource = vtkSmartPointer<vtkCubeSource>::New();//you can use sphere if you want

    vtkSmartPointer<vtkGlyph3D> glyph3D = vtkSmartPointer<vtkGlyph3D>::New();
    glyph3D->SetSourceConnection(cubeSource->GetOutputPort());
    glyph3D->ScalingOff();
    glyph3D->SetInputData(polydata);
    glyph3D->Update();

    // Create a mapper and actor
    vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    mapper->SetInputConnection(glyph3D->GetOutputPort());

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

    // Add the actor to the scene
    renderer->AddActor(actor);

    // Render
    renderWindow->Render();

    vtkSmartPointer<vtkInteractorStyleTrackballCamera> style = vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New(); //like paraview

    renderWindowInteractor->SetInteractorStyle( style );

    renderWindowInteractor->Start();
}

Snapshot of resulting rendering window is presented in Fig. 3 – and yes, I agree, it’s not very impressive. However, we can easily tweak it and make it significantly better.

plain3D

Figure 3. Exemplary 3D tumor visualized using 3D glyph mapper in VTK.

First of all to have better a better perspective we can slice the tumor in half and take resulting parts slightly apart. This can be achieved by modification of the addCells function (note that in my setting center of the tumor is at (x,y,z) = (250,250,250)).

void addCells(vector<Cell>::iterator begin, vector<Cell>::iterator end) {
    vector<Cell>::iterator cell = begin;
    for (cell; cell<end; cell++) {         if (cell->z < 250) {             points->InsertNextPoint(cell->x,cell->y,cell->z);
        } else {
            points->InsertNextPoint(cell->x,cell->y,cell->z+50);
        }
    }
}

Fig. 4 shows the snapshot of the rendering window after slicing the tumor – much better, but there is still place for easy improvement. Let’s add some color!

sliced
Figure 4. Exemplary 3D tumor visualized using 3D glyph mapper in VTK.

I assume that each cell has also information about the number of free spots in its direct neighborhood. We will color the tumor based on that information.

First, we need to add additional headers and variables.

#include "vtkCubeSource.h"
#include "vtkPolyDataMapper.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"
#include "vtkInteractorStyleTrackballCamera.h"
#include "vtkPolyData.h"
#include "vtkGlyph3DMapper.h"
#include "vtkPointData.h"
#include "vtkFloatArray.h"
#include "vtkLookupTable.h"

vtkSmartPointer<vtkRenderer> renderer;
vtkSmartPointer<vtkRenderWindow> renderWindow;
vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor;

vtkSmartPointer<vtkPoints> points;

vtkSmartPointer<vtkFloatArray> scalar;

Then we need to modify the function that initializes rendering

void initializeRendering() {
    renderer = vtkSmartPointer<vtkRenderer>::New(); //creating new renderer
    renderer->SetBackground(1., 1., 1.); //setting background color to white

    renderWindow = vtkSmartPointer<vtkRenderWindow>::New(); //creating renderer window
    renderWindow->SetSize(1000, 1000); //setting size of the window
    renderWindow->AddRenderer(renderer); //adding previously created renderer

    renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New(); //creating render
    renderWindowInteractor->SetRenderWindow(renderWindow); //adding previously created renderer window

    points = vtkSmartPointer<vtkPoints>::New();         //vector with points
    scalar = vtkSmartPointer<vtkFloatArray>::New(); //vector with scalar values of the points
    scalar->SetNumberOfComponents(1); //one value per point
}

and addCells function

void addCells(vector<Cell>::iterator begin, vector<Cell>::iterator end) {
    vector<Cell>::iterator cell = begin;
    for (cell; cell<end; cell++) {         if (cell->z < 250) {             points->InsertNextPoint(cell->x,cell->y,cell->z);
        } else {
            points->InsertNextPoint(cell->x,cell->y,cell->z+50);
        }

        scalar->InsertNextValue((double)cell->nSpots/26.);
    }
}

Finally, we can modify the visualization function.

void visualizeLesion(vector<Cell>::iterator begin, vector<Cell>::iterator end) {

    initializeRendering();

    addCells(begin, end);

    //creating new color table
    vtkSmartPointer<vtkLookupTable> colorLookupTable = vtkSmartPointer<vtkLookupTable>::New();
    colorLookupTable->SetNumberOfColors(256);
    colorLookupTable->SetHueRange( 0.667, 0.0);
    colorLookupTable->Build();

    // Combine into a polydata
    vtkSmartPointer<vtkPolyData> polydata = vtkSmartPointer<vtkPolyData>::New();
    polydata->SetPoints(points);
    polydata->GetPointData()->SetScalars(scalar);

    vtkSmartPointer<vtkCubeSource> cubeSource = vtkSmartPointer<vtkCubeSource>::New();//you can use sphere if you want

    vtkSmartPointer<vtkGlyph3D> glyph3D = vtkSmartPointer<vtkGlyph3D>::New();
    glyph3D->SetColorModeToColorByScalar();
    glyph3D->SetSourceConnection(cubeSource->GetOutputPort());
    glyph3D->ScalingOff();
    glyph3D->SetInputData(polydata);
    glyph3D->Update();

    // Create a mapper and actor
    vtkSmartPointer<vtkPolyDataMapper> mapper = vtkSmartPointer<vtkPolyDataMapper>::New();
    mapper->SetInputConnection(glyph3D->GetOutputPort());
    mapper->SetLookupTable(colorLookupTable);

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

    // Add the actor to the scene
    renderer->AddActor(actor);

    // Render
    renderWindow->Render();

    vtkSmartPointer<vtkInteractorStyleTrackballCamera> style = vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New(); //like paraview

    renderWindowInteractor->SetInteractorStyle( style );

    renderWindowInteractor->Start();
}

Snapshot of sliced and colored tumor is presented in Fig. 5 – fast, good-looking and easy.
slicedColor
Figure 5. Exemplary 3D tumor visualized using 3D glyph mapper in VTK. Color represents the number of free spots in the direct neighborhood of each cell (jet colormap).

Advertisements

One thought on “Some tips&tricks for 3D tumor visualization

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s