Render 3D DICOM object with marching cube?

[truongleit]

Hi, i am trying to visualize a set of 2D DICOM slices into 3D object. But the important thing is I have to apply the marching cube algorithm. I have been searching for solutions on many sites but found nothing, also, I could see any clear tutorials about my problem on VTK.js website. Can you guys give me some suggestion? @thewtex

[jourdain]

Here is an example of doing marching cube of an ImageData. You just need to load your DICOM like ParaView Glance is doing using itk.js and you should be good to go.

[thewtex]

Adding the example that @jourdain points to which shows how to execute marching cubes, the other calls you need to make are readImageDICOMFileSeries and convertItkToVtkImage.

[truongleit]

Hi again, thanks for the clear explanation. The example is really clear and useful. So in general (considering the example attached by @jourdain ), I need to read the DICOM series, then convert them into a variable by using convertItkToVtkImage. Afterward, the file with VTI format in the example will be replace by assigned variable, then I will get what I need, is it right?

[thewtex]

@truongleit you got it :+1:

[truongleit]

Hi again, need some helps from you guys! I followed the code from itk-vtk-viewer to see how you read and convert the input files into the VTK image, then tried to do the apply it to the VolumeContour example. But in the final step, I receive an error which is:

"Uncaught (in promise) TypeError: a.getOrigin is not a function"

Here is my code:

function render(files) {
    let reader = itkreadImageDICOMFileSeries;
    let arg = files;
    reader(null, arg).then(({
        image: itkImage,
        webWorker
    }) => {
        webWorker.terminate()
        console.log(itkImage)

        const imageData = vtkITKHelper.convertItkToVtkImage(itkImage);
        console.log(imageData)

        const dataRange = imageData
            .getPointData()
            .getScalars()
            .getRange();
        const firstIsoValue = (dataRange[0] + dataRange[1]) / 3;

        const el = document.querySelector('.isoValue');
        el.setAttribute('min', dataRange[0]);
        el.setAttribute('max', dataRange[1]);
        el.setAttribute('value', firstIsoValue);
        el.addEventListener('input', updateIsoValue);

        marchingCube.setContourValue(firstIsoValue);
        renderer.addActor(actor);
        renderer.getActiveCamera().set({
            position: [1, 1, 0],
            viewUp: [0, 0, -1]
        });
        renderer.resetCamera();
        renderWindow.render();
    })
}

The imageData variable is the one that will be rendered. I also compared the imageData in my code with the imageData in itk-vtk-viewer. They are exactly similar when loading the same DICOM serie. Note that the VolumeContour example is running perfectly on my device. @thewtex @jourdain

[jourdain]

Where do you call marchingCube.setInputData(imageData)?

[truongleit]

HUGE thanks to you! It is working now 👍 @jourdain

[truongleit]

Hi again! I would like to render the 3D object into a specific container. In my case, the div class is 'threeD', how can I do that? I read the code and found out that it's now using vtkFullScreenRenderWindow to visualize object in fullscreen mode. But I don't know what I should do next, read some more examples on the site but most of them are using fullscreen mode as well. Hope to reiceive your help soon, thanks again. @jourdain

Here is my code:

const fullScreenRenderWindow = vtkFullScreenRenderWindow.newInstance({
    background: [0, 0, 0],
});
const renderWindow = fullScreenRenderWindow.getRenderWindow();
const renderer = fullScreenRenderWindow.getRenderer();

// fullScreenRenderWindow.addController(controlPanel);

const actor = vtkActor.newInstance();
const mapper = vtkMapper.newInstance();
const marchingCube = vtkImageMarchingCubes.newInstance({
    contourValue: 0.0,
    computeNormals: true,
    mergePoints: true,
});

actor.setMapper(mapper);
mapper.setInputConnection(marchingCube.getOutputPort());

function updateIsoValue(e) {
    const isoValue = Number(e.target.value);
    marchingCube.setContourValue(isoValue);
    renderWindow.render();
}

const reader = vtkHttpDataSetReader.newInstance({
    fetchGzip: true
});
marchingCube.setInputConnection(reader.getOutputPort());

function render(files) {
    let reader = itkreadImageDICOMFileSeries;
    let arg = files;
    reader(null, arg).then(({
        image: itkImage,
        webWorker
    }) => {
        webWorker.terminate()
        console.log(itkImage)

        const imageData = vtkITKHelper.convertItkToVtkImage(itkImage);
        console.log(imageData)

        const dataRange = imageData
            .getPointData()
            .getScalars()
            .getRange();
        const firstIsoValue = (dataRange[0] + dataRange[1]) / 3;

        marchingCube.setInputData(imageData);
        marchingCube.setContourValue(firstIsoValue);
        renderer.addActor(actor);
        renderer.getActiveCamera().set({
            position: [1, 1, 0],
            viewUp: [0, 0, -1]
        });
        renderer.resetCamera();
        renderWindow.render();
    })
}

var fileInput = document.getElementById('file-input');
fileInput.addEventListener('change', function(event) {
    var input = event.target;
    render(input.files);
});

global.fullScreen = fullScreenRenderWindow;
global.actor = actor;
global.mapper = mapper;
global.marchingCube = marchingCube;

[jourdain]

https://kitware.github.io/vtk-js/examples/SimpleCone.html

[truongleit]

Hi again, is there any configs except iso value for marching cube rendering? Is it possible to add the control UI from vtkVolumeController to the render so I can change the 3D object's colors? @jourdain

[jourdain]

You can do what you want, but the generated polydata won't have any field to color by since you contour (same value) it.

You can set the color to its actor using actor.getProperty().setColor(r, g, b) where r,g,b are float between [0,1].

[truongleit]

Thanks for suggestion! Also, I am building a color-changing function based on your instruction, so how can I make the object's color change immediately right after I choose new color? Currently I am testing it by typing directly the code in Chrome DevTools, so after I modify the color's value of actor, I need to do something (ex: rotate) with the object's canvas in order to see the change. @jourdain

[jourdain]

You need to call render() on the vtkRenderWindow instance.

[truongleit]

Hi again, I am thinking of creating 4 view modes for my application which contains both 3D and 2D (X, Y, Z) objects. Thanks to your clear explanation, the 3D mode is now performing pretty well. Now I am moving to visualize the input collection of DICOM slices into 3 axes (X, Y and Z). However, the only example I found on VTK.JS site isMultiSliceImageMapper. The problem of mine is I need to display them in 3 different DOM elements, not like the given example. Can you instruct me a little bit more? Thanks @jourdain

[jourdain]

You can look at ParaView Glance which already do all of that.

[truongleit]

Thanks a lot! My app also has other rendering algorithms which are implemented by using other JavaScript library. When switching from VTK's Marching Cube to another, how can I destroy the render for memory-saving purposes? @jourdain

[SaneSoup]

Hello there,

I know this issue was resolved 5 years ago, however I also have a problem regarding the vtk library.

My goal is it to have a python code that reads DICOM files and then turns them into a 3D .obj file, the problem is that I don't know how the values of vtkMarchingCubes work because I can't seem to get all the data but only parts of it. So to say, sometimes I only get the Skull as a model or just the skin, the important thing for me would be the brain tho.

I have tried to render everything in python before, where nothing went missing and at least this part worked well but because of other problems with vtk not liking tk inter I made up my mind to use Godot for the UI stuff. Now using Godot I need the before mentioned .obj file to kind of import it into the 3D environment as a node or something like that.

(I'm a Godot noob so I also have no idea on how to implement the .obj but first I gotta find out how to solve my python problem)

My first full python try looked like this (the # are in german):

import tkinter as tk
from tkinter import filedialog
import vtk

class DICOMViewer:
    def __init__(self, root):
        self.root = root
        self.root.title("DICOM Viewer")

        # Initialisieren Sie die VTK-Variable
        self.vtk_window = None

        # Menü erstellen
        menu_bar = tk.Menu(root)
        root.config(menu=menu_bar)

        file_menu = tk.Menu(menu_bar, tearoff=0)
        menu_bar.add_cascade(label="File", menu=file_menu)
        file_menu.add_command(label="Open DICOM", command=self.open_dicom_folder)

    def open_dicom_folder(self):
        folder_path = filedialog.askdirectory(title="Open DICOM Folder")

        if folder_path:
            # Verzögere das Laden der VTK-Operationen im Hauptthread
            self.root.after(0, self.load_dicom_and_render, folder_path)

    def load_dicom_and_render(self, folder_path):
        # Laden aller DICOM-Dateien im ausgewählten Ordner
        dicom_reader = vtk.vtkDICOMImageReader()
        dicom_reader.SetDirectoryName(folder_path)
        dicom_reader.Update()
        dicom_volume = dicom_reader.GetOutput()

        # 3D-Rekonstruktion
        volume_mapper = vtk.vtkGPUVolumeRayCastMapper()
        volume_mapper.SetInputData(dicom_volume)

        volume_property = vtk.vtkVolumeProperty()
        volume_property.ShadeOn()
        volume_property.SetInterpolationTypeToLinear()

        volume_color = vtk.vtkColorTransferFunction()
        volume_color.AddRGBPoint(-1024, 0.0, 0.0, 0.0)
        volume_color.AddRGBPoint(0, 1.0, 0.5, 0.5)
        volume_color.AddRGBPoint(3071, 1.0, 1.0, 1.0)

        volume_scalar_opacity = vtk.vtkPiecewiseFunction()
        volume_scalar_opacity.AddPoint(-1024, 0.0)
        volume_scalar_opacity.AddPoint(0, 0.2)
        volume_scalar_opacity.AddPoint(3071, 0.2)

        volume_property.SetColor(volume_color)
        volume_property.SetScalarOpacity(volume_scalar_opacity)

        volume_actor = vtk.vtkVolume()
        volume_actor.SetMapper(volume_mapper)
        volume_actor.SetProperty(volume_property)

        # Rendering
        renderer = vtk.vtkRenderer()
        renderer.SetBackground(0.0, 0.0, 0.0)

        render_window = vtk.vtkRenderWindow()
        render_window.SetWindowName("DICOM 3D Viewer")
        render_window.SetSize(self.root.winfo_screenwidth() // 2, self.root.winfo_screenheight())
        render_window.AddRenderer(renderer)

        interactor = vtk.vtkRenderWindowInteractor()
        interactor.SetRenderWindow(render_window)

        renderer.AddVolume(volume_actor)
        renderer.ResetCamera()

        # Interaktionsstil erstellen
        interactor_style = vtk.vtkInteractorStyleTrackballCamera()
        interactor.SetInteractorStyle(interactor_style)

        # Ändere die Position des VTK-Fensters
        self.root.update_idletasks()
        render_window.SetPosition(self.root.winfo_width(), 0)

        # Speichern Sie die Referenz auf das VTK-Fenster
        self.vtk_window = render_window

        # Starte die Interaktion
        interactor.Start()

if __name__ == "__main__":
    root = tk.Tk()
    app = DICOMViewer(root)

    # Setze die Tkinter-GUI auf die gesamte linke Seite
    root.geometry("{}x{}+0+0".format(root.winfo_screenwidth() // 2, root.winfo_screenheight()))

    root.mainloop()

my current attempt looks more like this:

import vtk

def load_dicom_series(a):
    reader = vtk.vtkDICOMImageReader()
    reader.SetDirectoryName(a)
    reader.Update()
    return reader

def create_3d_model(b):
    # Apply Marching Cubes to extract the surface
    marching_cubes = vtk.vtkMarchingCubes()
    marching_cubes.SetInputConnection(b.GetOutputPort())
    marching_cubes.SetValue(10, 1500)  # Adjust this value to include more/less of the data
    marching_cubes.Update()
    return marching_cubes

def save_obj(c, d):
    # Convert to PolyData
    poly_data = d.GetOutput()

    # Write to OBJ file
    obj_writer = vtk.vtkOBJWriter()
    obj_writer.SetFileName(c)
    obj_writer.SetInputData(poly_data)
    obj_writer.Write()

# Pfad zur DICOM-Serie
dicom_path = 'C:\\Users\\User\\Desktop\\python\\privat\\Germanische\\Scans'
output_path = 'C:\\Users\\User\\Desktop\\model.obj'

# Laden der DICOM-Daten
dicom_reader = load_dicom_series(dicom_path)

# Erzeugung des 3D-Modells
e = create_3d_model(dicom_reader)

# Speichern des 3D-Modells in einer .obj-Datei
save_obj(output_path, e)

[finetjul]

It's not clear what your question/problem is.

Is your issue to pick the marching cube "value" ? if so, it is expected. DICOMs typically have a different voxel values for "brain", you can't hardcode it. You would need to let the user "find it" interactively.

(you might want to use trame instead of godot)

[SaneSoup]

Hey @finetjul thanks for your response, my issue is indeed picking the marching cube value. It doesn't really matter that I can't hardcode it for just the brains but there surely has to be a way for everything to be in my 3D model at once (meaning that all of the skin, bones, and organs show up and not just parts of them) And while we are at it, how do these voxel values even work?

For example the following .obj files were done with the values (0, 100), (0,500), (0, 60), (0,10): (this is made with the second code from before + I viewed the files with an online viewer)

It should look more like this with a filled and not empty head: (this is from my pure python code, which was the first one from before)

[finetjul]

vtkMarchingCube creates a "surface" mesh, that's why it will ALWAYS be empty. Your python results show "volume rendering". It's not clear why you want to extract a mesh (triangles) from an existing 3D image (i.e. volume).

[SaneSoup]

Thats not good :( is there any way to get the whole thing and not just the surface?

And the "volume rendering" stuff was just a different approach and a test to see if the whole program would be possible in just python, which turned out to not work, however it looks good. So it has nothing to do with the MarchingCubes approach, it's a different code.

[jadh4v]

Looks like you are trying to re-create volume rendering style visualization using polygonal data. This is not a trivial problem to solve. You could get a rough visual approximate of the volume rendering through multiple translucent surfaces. I'd start by extracting multiple iso-contours. You can set the number of contours first and then set multiple values through calls to marching_cubes.SetValue as below:

marching_cubes.SetNumberOfContours(4)
marching_cubes.SetValue(0, 10)
marching_cubes.SetValue(1, 60)
marching_cubes.SetValue(2, 100)
marching_cubes.SetValue(3, 500)

Rather than doing this programmatically, a better approach could be to do this through Paraview's user interface. This way, you can tune your values and also tweak rendering properties. You might also be able to export the mesh directly to obj from Paraview.

[SaneSoup]

@jadh4v I wasn't trying to do that, I just used different values before to explain my standpoint + I didn't even know that it'll just be a surface mesh. I actually just want a method to get a full 3D Model of combined DICOM files. It's not obligatory to be a .obj file afterwards it just has to be compatible with Godot.

[finetjul]

You might want to 3D Slicer to extract the brain from your image. You can look at those extensions:

• https://github.com/lassoan/SlicerSwissSkullStripper
• https://extensions.slicer.org/view/HDBrainExtraction/32919/win

[SaneSoup]

Thank you @finetjul the source code of these extensions might be usefull in a bit.