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.

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.

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!

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.

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).

Advertisement

Visualization of 3D tumor using isosurfaces and simple blur

Today, I’ll show how to simply visualize 3D tumor using MATLAB’s isosurface function combined with a simple blurring technique.

First, we will modify the code from the previous post in order to simulate tumor in 3D. There are only a few changes: redefinition of the domain and the neighborhood; generating permutations on the fly instead of storing them in Pms varaible; and the way in which the initial cell is placed.

N = 200; %cubic domain dimension
nSteps = 30*24; %number of simulation steps

pprol = 1/24; %probability of proliferating
pmig = 10/24; %probability of migrating
pdeath = 1/100; %probability of dying
pmax = 10; %proliferation capacity
ps = 3/10; %probability of symmetric division

L = false(N,N,N); %domain definition
L(1,:,:) = true; L(end,:,:) = true; %filling boundary
L(:,1,:) = true; L(:,end,:) = true; %filling boundary
L(:,:,1) = true; L(:,:,end) = true; %filling boundary

L(N*N*round(N/2)+N*round(N/2)+round(N/2)) = true;
cells = int32(N*N*round(N/2)+N*round(N/2)+round(N/2));
cellsIsStem = true;
cellsPmax = uint8(pmax);

aux = int32([[-N-1 -N -N+1 -1 1 N-1 N N+1] ...
             [-N-1 -N -N+1 0 -1 1 N-1 N N+1]-N*N ...
             [-N-1 -N -N+1 0 -1 1 N-1 N N+1]+N*N])'; %indices to heighborhood

for i = 1:nSteps
    sh = randperm(length(cells));
    cells = cells(sh);
    cellsIsStem = cellsIsStem(sh);
    cellsPmax = cellsPmax(sh);

    Pms = cell2mat(arrayfun(@(x)randperm(26),(1:length(cells))','UniformOutput',0))';
    S = bsxfun(@plus,cells,aux(Pms));
    S(L(S)) = 0; %setting occupied spots to 0
    indxF = find(any(S)); %selecting cells with at least one free spot
    nC = length(indxF); %number of cells with free spot
    
    P = rand(1,nC)<pprol; %proliferation
    Ps = P & rand(1,nC)<ps & cellsIsStem(indxF);%symmetric division
    De = P & (cellsPmax(indxF) == 0);%proliferation exhaution
    D = P & (rand(1,nC)<pdeath) & ~cellsIsStem(indxF); %death at proliferation attempt
    M = ~P & (rand(1,nC)<pmig); %go when no grow
    
    del = D | De; %cells to delete
    act = find((P | M) & ~del); %indices to the cells that will do something
    for ii = act %only for those that will do anything
        ngh = S(:,indxF(ii)); 
        ngh(ngh==0) = [];
        indO = find(~L(ngh),1,'first'); %selecting free spot
        if ~isempty(indO) %if there is still a free spot
            L(ngh(indO)) = true;
            if P(ii) %proliferation
                cells = [cells uint32(ngh(indO))];
                if Ps(ii) %symmetric division
                   cellsIsStem = [cellsIsStem true];
                   cellsPmax = [cellsPmax cellsPmax(indxF(ii))];  
                else
                   cellsIsStem = [cellsIsStem false];
                   cellsPmax = [cellsPmax cellsPmax(indxF(ii))-1];
                   if ~cellsIsStem(indxF(ii))
                    cellsPmax(indxF(ii)) = cellsPmax(indxF(ii))-1;
                   end
                end
            else %migration
                L(cells(indxF(ii))) = false;
                cells(indxF(ii)) = uint32(ngh(indO));
            end
        end
    end
    
    if ~isempty(del) %updating death
        L(cells(indxF(del))) = false;
        cells(indxF(del)) = [];
        cellsIsStem(indxF(del)) = [];
        cellsPmax(indxF(del)) = [];
    end 
end

Using the above code I simulated two tumors that have different growth characteristics, both with pmig = 5/24 and one with ps = 3/10 and the other one with ps = 5/100. Now we need to visualize them.

Having the whole tumor stored in the L variable we can use isosurface function straight away and plot the tumor as it is. The only thing that we need to do is to remove the cells from the L boundary.

function visualize( N, L)

%clearing cells from the boundary
L(1,:,:) = false; L(end,:,:) = false;
L(:,1,:) = false; L(:,end,:) = false;
L(:,:,1) = false; L(:,:,end) = false;

%calculating isosurfaces and plotting
p = patch(isosurface(1:N,1:N,1:N,L,0.25));
isonormals(1:N,1:N,1:N,L,p)
set(p,'FaceColor','red','EdgeColor','none');

xlim([1 N]);
ylim([1 N]);
zlim([1 N]);
view([90 0]);
camlight 
lighting gouraud    
end

The results of the above function are plotted below.

As you can see the visualization doesn’t really allow to see the 3D structure of the tumors, as individual separated cells are plotted and there is no clear 3D surface. We will fix that by using a simple blurring technique. The basic idea is that we will replace each site value by the average value from its neighborhood (several times in a loop). That will allow to delete separated cells and make 3D surface smoother. We add an input variable blLevels to the function in order to define how may iterations of that procedure should be performed.

function visualize( N, L, cells, blLevels )

%clearing cells from the boundary
L(1,:,:) = false; L(end,:,:) = false;
L(:,1,:) = false; L(:,end,:) = false;
L(:,:,1) = false; L(:,:,end) = false;

if blLevels %if perform basic blur

%auxilary variable with indices to the cell negiborhood
aux = int32([[-N-1 -N -N+1 -1 1 N-1 N N+1] ...
             [-N-1 -N -N+1 0 -1 1 N-1 N N+1]-N*N ...
             [-N-1 -N -N+1 0 -1 1 N-1 N N+1]+N*N])';

%creating indices to cells and their beignorhoods
S = [cells unique(reshape(bsxfun(@plus,cells,aux),1,[]))];
%creating indices to neigborhood of indices in S
S2 = bsxfun(@plus,S,aux);
%making sure that indices are still within the lattice
S2(S2<1) = []; S2(S2>N*N*N) = [];

    %changing lattice from logical variable to float
    L = single(L);
    for i = 1:blLevels %for number of blurs
        L(S) = mean(L(S2)); %taking the average of neighborhood
    end
end

%calculating isosurfaces and plotting
p = patch(isosurface(1:N,1:N,1:N,L,0.25));
isonormals(1:N,1:N,1:N,L,p)
set(p,'FaceColor','red','EdgeColor','none');

xlim([1 N]);
ylim([1 N]);
zlim([1 N]);
view([90 0]);
camlight 
lighting gouraud    
end

If we apply blLevels = 1 we get the following, nicer plot of the tumor.

For blLevels = 2 we get the following plots.

Finally, when applying three iterations of blur (blLevels = 3) we get 3D plots nicely showing the tumor shape.

What is also great it takes only about 3 seconds to plot a tumor with more than quarter of a million of cells (left one).