VTKを用いて3次元Numpy配列をボクセルボリュームとして可視化

Numpyを用いて生成した3次元行列をそれぞれボクセルとみなして,ボリュームの可視化をしたかったため実装した.

vtkHexahedronでの実装やmatplotで実装しようとしたところ,激重だったため,一度vtkArrayにコンバートしている.いろいろなサイトを見ながらコードを書いたけど,どこを参考にしたかは忘れてしまいました.すんません.

   

以下,必須インポート

import vtk
import numpy as np

以下,使用例

hoge = vtkNumpyVolumeRender(ndarrayVolume)
hoge.visualize()

で使えるはず.HSVαのレンジを指定することでスカラ値からどのように着色するかを示すカラーテーブルを変更できる.コンストラクタで,SpacingとOriginも変更可能.IS_PHOTO_MODE=Trueにすると,スクショとかを取りやすいビューになる.

以下,クラス.

class vtkNumpyVolumeRender():

    def numpyArray2vtkArray(self, numpyArray, maxValue):
        # np.uint8での可視化は255が上限
        if maxValue > 255:
            raise ValueError('Only values between 0 and 255 allowed.')
        # nanを1に置換
        ndarray[np.isnan(numpyArray)] = 0
        # [0-1]の値を[0-maxValue]に拡張
        ndarray = ndarray.astype(np.float) / np.max(ndarray) * maxValue
        ndarray = ndarray.astype(np.uint8)
        # needs to be done due to different array ordering in numpy and vtk
        ndarray = ndarray.T
        return ndarray

    def initDataImporter(self, vtkArray, origin, spacing):
        dataImporter = vtk.vtkImageImport()
        data_string = vtkArray.tostring()
        dataImporter.CopyImportVoidPointer(data_string, len(data_string))
        dataImporter.SetDataScalarTypeToUnsignedChar()
        dataImporter.SetNumberOfScalarComponents(1)
        w, d, h = vtkArray.shape
        dataImporter.SetDataExtent(0, h - 1, 0, d - 1, 0, w - 1)
        dataImporter.SetWholeExtent(0, h - 1, 0, d - 1, 0, w - 1)
        dataImporter.SetDataOrigin(origin[0], origin[1], origin[2])
        dataImporter.SetDataSpacing(spacing[0], spacing[1], spacing[2])
        return dataImporter

    def initColorTransferFunction(self, componentValueRange, hueRange, saturationRange, valueRange, steps=100):
        volumeColor = vtk.vtkColorTransferFunction()
        if self.IS_PHOTO_MODE:
            valueRange = (0,0)
        for i in range(steps + 1):
            componentValue = i / steps * (componentValueRange[1]-componentValueRange[0]) + componentValueRange[0]
            h = i / steps * (hueRange[1]-hueRange[0]) + hueRange[0]
            s = i / steps * (saturationRange[1]-saturationRange[0]) + saturationRange[0]
            v = i / steps * (valueRange[1]-valueRange[0]) + valueRange[0]
            volumeColor.AddHSVPoint(componentValue, h, s, v)
        volumeColor.SetNanColor(0, 0, 0)
        return volumeColor

    def initScalarOpacity(self, componentValueRange, opacityRange):
        volumeScalarOpacity = vtk.vtkPiecewiseFunction()
        steps = 100
        for i in range(steps + 1):
            componentValue = i / steps * (componentValueRange[1] - componentValueRange[0]) + componentValueRange[0]
            alpha = i / steps * (opacityRange[1]-opacityRange[0]) + opacityRange[0]
            volumeScalarOpacity.AddPoint(componentValue, alpha)
        return volumeScalarOpacity

    def initGradientOpacity(self, componentValueRange):
        # 輝度値によるグラデーション具合を操作する
        volumeGradientOpacity = vtk.vtkPiecewiseFunction()
        steps = 100
        for i in range(steps + 1):
            componentValue = i / steps * (componentValueRange[1] - componentValueRange[0]) + componentValueRange[0]
            gradiendt = i / steps
            volumeGradientOpacity.AddPoint(componentValue, gradiendt)
        return volumeGradientOpacity

    def initVolumeProperty(self, colorTransferFunction, scalarOpacity):
        # ボリュームレンダリングのパラメータ
        volumeProperty = vtk.vtkVolumeProperty()
        volumeProperty.SetInterpolationTypeToLinear()  # 欠損部分の線形補完
        volumeProperty.SetColor(colorTransferFunction)
        volumeProperty.SetScalarOpacity(scalarOpacity)
        # volumeProperty.SetGradientOpacity(volumeGradientOpacity)
        # volumeProperty.ShadeOn() # 指向性照明と影の描画
        # volumeProperty.SetAmbient(0.3) # sharderのパラメータ
        # volumeProperty.SetDiffuse(0.2) # sharderのパラメータ
        # volumeProperty.SetSpecular(0.2) # sharderのパラメータ
        return volumeProperty

    def initVolumeMapper(self, inputData):
        volumeMapper = vtk.vtkSmartVolumeMapper()
        volumeMapper.SetInputConnection(inputData.GetOutputPort())
        return volumeMapper

    def initVolume(self, volumeMapper, volumeProperty):
        volume = vtk.vtkVolume()
        volume.SetMapper(volumeMapper)
        volume.SetProperty(volumeProperty)
        return volume

    def initRenderer(self, volume, SET_AXIS=True):
        ren = vtk.vtkRenderer()
        ren.AddVolume(volume)
        ren.SetBackground(1, 1, 1)
        if not self.IS_PHOTO_MODE:
            ren.GradientBackgroundOn()
            ren.SetBackground2(0.4, 0.55, .75)
            if SET_AXIS:
                axes = vtk.vtkAxesActor()
                ren.AddActor(axes)
        return ren

    def initCamera(self, renderer, volume, zoom_rate = 20):
        srideLen = self.slice_num * zoom_rate
        camera = renderer.GetActiveCamera()
        volumeCenter = volume.GetCenter()
        camera.SetPosition(volumeCenter[0]+self.slice_num/2+srideLen, volumeCenter[1], volumeCenter[2])
        camera.SetFocalPoint(volumeCenter[0], volumeCenter[1], volumeCenter[2])

    def initRenderingWindow(self, renderer, windowSize = (960, 720)):
        ren_win = vtk.vtkRenderWindow()
        ren_win.AddRenderer(renderer)
        ren_win.SetSize(windowSize[0], windowSize[1])
        return ren_win

    def initInteracter(self, ren_win):
        inter = vtk.vtkRenderWindowInteractor()
        inter.SetRenderWindow(ren_win)
        inter.SetInteractorStyle(MouseInteractorStyle())
        return inter

    def initScalarBar(self, volumeColor):
        scalar_bar = vtk.vtkScalarBarActor()
        scalar_bar.SetLookupTable(volumeColor)
        scalar_bar.SetOrientationToHorizontal()
        return scalar_bar

    def initScalarBarWidget(self, interacter, scalar_bar):
        scalar_bar_widget = vtk.vtkScalarBarWidget()
        scalar_bar_widget.SetInteractor(interacter)
        scalar_bar_widget.SetScalarBarActor(scalar_bar)
        scalar_bar_widget.On()
        return scalar_bar_widget

    def __init__(self, numpy3dArray, spacing=(1,1,1), origin=(0,0,0), maxValue = 255, IS_PHOTO_MODE = False):
        self.IS_PHOTO_MODE = IS_PHOTO_MODE
        self.slice_num, self.height_num, self.width_num = numpy3dArray.shape[:3]
        vtkArray = self.numpyArray2vtkArray(numpy3dArray, maxValue)
        minComponentValue = np.min(vtkArray)
        maxComponentValue = np.max(vtkArray)

        dataImporter = self.initDataImporter(vtkArray, origin, spacing)


        hueRange = (0, 0.55) # range[0, 1]
        saturationRange = (0.8, 0.8) # range[0, 1]
        valueRange = (1, 1) # range[0, 1]
        componentValueRange = (minComponentValue, maxComponentValue)
        volumeColor = self.initColorTransferFunction(componentValueRange, hueRange, saturationRange, valueRange)

        alphaRange = (1, 0)
        componentValueRange = (0, maxComponentValue)
        volumeScalarOpacity = self.initScalarOpacity(componentValueRange, alphaRange)

        volumeGradientOpacity = self.initGradientOpacity(componentValueRange)

        volumeProperty = self.initVolumeProperty(volumeColor, volumeScalarOpacity)

        volumeMapper = self.initVolumeMapper(dataImporter)

        self.volume = self.initVolume(volumeMapper, volumeProperty)

        # Renderer
        self.ren = self.initRenderer(self.volume, SET_AXIS=True)

        # camera
        self.initCamera(self.ren, self.volume)

        # Rendering Window
        self.ren_win = self.initRenderingWindow(self.ren)

        # Interactor
        self.inter = self.initInteracter(self.ren_win)

        # Scalar bar
        if not self.IS_PHOTO_MODE:
            scalarBar = self.initScalarBar(volumeColor)
            self.scalarBarWidget = self.initScalarBarWidget(self.inter, scalarBar)

        self.ren_win.Render()

    def setFocalPoint(self, divCenter=(0,0,0)):
        volumeCenter = self.volume.GetCenter()
        camera = self.ren.GetActiveCamera()
        x = volumeCenter[0] + divCenter[0]
        y = volumeCenter[1] + divCenter[1]
        z = volumeCenter[2] + divCenter[2]
        camera.SetFocalPoint(x, y, z)

    def visualize(self):
        self.inter.Initialize()
        self.inter.Start()
   
カテゴリー: プログラミング タグ: , , パーマリンク