liviaerxin/Dicom3dCube.py

## README.md

      
    Raw
  

              README.md
            
          
    DICOM Files Read, 3D Reconstruction and Visualization


vtk read --> vtk 3D Visualization
pydicom read --> numpy 3D Reconstruction --> vtk 3D Visualization

DICOM Processing and Segmentation in Python – Radiology Data Quest
Advanced DICOM-CT 3D visualizations with VTK | Kaggle
Pulmonary Dicom Preprocessing | Kaggle
How to actually to convert 3d numpy array of DICOM images into vtkImageData and vtkAlgorithmOutput? - #2 by dgobbi - Support - VTK

  
## Dicom3dCube.py
"""
DICOM files --> vtk array --> vtk 3D visualization

1. read a series of DICOM files by vtkDICOMImageReader(which will do HU and don't resample)
2. process vtk 3D array with marching cubes algorithm
3. render vtk 3D array
"""
import vtk
from utils import CreateTissue


def main():
    dicom_dir = get_program_parameters()

    reader = vtk.vtkDICOMImageReader()
    reader.SetDirectoryName(dicom_dir)
    reader.Update()

    # Setup render
    renderer = vtk.vtkRenderer()
    renderer.AddActor(CreateTissue(reader, -900, -400, "lung"))
    renderer.AddActor(CreateTissue(reader, 0, 120, "blood"))
    renderer.AddActor(CreateTissue(reader, 100, 1000, "skeleton"))
    renderer.SetBackground(1.0, 1.0, 1.0)
    renderer.ResetCamera()

    # Setup render window
    renWin = vtk.vtkRenderWindow()
    renWin.AddRenderer(renderer)
    renWin.SetSize(800, 800)
    renWin.Render()

    # Setup render window interactor
    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)

    # Setup coordinate axes helper
    axesActor = vtk.vtkAxesActor()
    widget = vtk.vtkOrientationMarkerWidget()
    rgba = [0] * 4
    colors = vtk.vtkNamedColors()
    colors.GetColor("Carrot", rgba)
    widget.SetOutlineColor(rgba[0], rgba[1], rgba[2])
    widget.SetOrientationMarker(axesActor)
    widget.SetInteractor(iren)
    widget.SetViewport(0.0, 0.0, 0.4, 0.4)
    widget.SetEnabled(1)
    widget.InteractiveOn()

    iren.Initialize()
    iren.Start()


def get_program_parameters():
    import argparse

    description = "DICOM Directory including `*.dcm` files"
    epilogue = """
    Derived from VTK/Examples/Cxx/Medical1.cxx
    This example reads a volume dataset, extracts an isosurface that
     represents the skin and displays it.
    """
    parser = argparse.ArgumentParser(
        description=description,
        epilog=epilogue,
        formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    parser.add_argument("dicom_dir", help="patients/series/")
    args = parser.parse_args()
    return args.dicom_dir


if __name__ == "__main__":
    main()

## Dicom3dCubeWithNumpy.py
"""
DICOM files --> 3D numpy array --> vtk array --> vtk 3D visualization

1. read a series of DICOM files by pydicom
2. create a 3D numpy array(NHW) by stacking a series of 2D DICOM image data(HW) and preprocessing(hu, resample,..etc) below
    1. HU
    2. resample
3. convert the 3D numpy to vtk 3D array(vtkImageData)
4. bridge vtk 3D array `vtkImageData` to `vtkAlgorigthmOutput` by `vtkTrivialProducer`, then process with marching cubes algorithm
5. render vtk 3D array
"""


import numpy as np

import vtk
from vtk.util import numpy_support
from utils import CreateTissue, get_pixels_hu, load_scan, resample


def main():
    dicom_dir = get_program_parameters()

    # Load the DICOM files
    slices = load_scan(dicom_dir)

    # pixel aspects, assuming all slices are the same
    ps = slices[0].PixelSpacing
    ss = slices[0].SliceThickness
    ax_aspect = ps[1] / ps[0]
    sag_aspect = ps[1] / ss
    cor_aspect = ss / ps[0]

    print(
        f"PixelSpacing[{ps}] SliceThickness[{ss}] AxialAspect[{ax_aspect}] SagittalAspect[{sag_aspect}] CoronalAspect[{cor_aspect}]"
    )

    # Convert raw data into Houndsfeld units
    img3d = get_pixels_hu(slices)

    print(f"img3d [{img3d.shape}] [{img3d.dtype}] [{img3d.nbytes}] [{np.mean(img3d)}] ")

    # Resample
    img3d_after_resamp, spacing = resample(img3d, slices)

    print(
        f"img3d_after_resamp [{img3d_after_resamp.shape}] [{img3d_after_resamp.dtype}] [{img3d_after_resamp.nbytes}] [{np.mean(img3d_after_resamp)}] "
    )

    # Prepare 3D visualization
    print(f"3D visualization...")

    # Convert numpy array to vtk array
    img3d_64f = img3d_after_resamp.astype(np.float64)
    num, h, w = img3d_64f.shape
    print(f"img3d_64f [{img3d_64f.shape}] [{np.max(img3d_64f[1])}]")

    imageData = vtk.vtkImageData()
    depthArray = numpy_support.numpy_to_vtk(
        num_array=img3d_64f.ravel(), deep=True, array_type=vtk.VTK_FLOAT
    )
    imageData.SetDimensions((w, h, num))  # x,y,z coordination system
    imageData.SetSpacing([1, 1, 1])
    # imageData.SetSpacing([ps[0], ps[1], ss])
    imageData.SetOrigin([0, 0, 0])
    imageData.GetPointData().SetScalars(depthArray)

    image_producer: vtk.vtkAlgorithm = vtk.vtkTrivialProducer()  # data provider
    image_producer.SetOutput(imageData)

    # Rendering

    # Setup render
    renderer = vtk.vtkRenderer()
    renderer.AddActor(CreateTissue(image_producer, -900, -400, "lung"))
    renderer.AddActor(CreateTissue(image_producer, 0, 120, "blood"))
    renderer.AddActor(CreateTissue(image_producer, 100, 2000, "skeleton"))
    renderer.SetBackground(1.0, 1.0, 1.0)
    renderer.ResetCamera()

    # Setup render window
    renWin = vtk.vtkRenderWindow()
    renWin.AddRenderer(renderer)
    renWin.SetSize(800, 800)
    renWin.Render()

    # Setup render window interaction
    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)

    # Setup coordinate axes helper
    axesActor = vtk.vtkAxesActor()
    widget = vtk.vtkOrientationMarkerWidget()
    rgba = [0] * 4
    colors = vtk.vtkNamedColors()
    colors.GetColor("Carrot", rgba)
    widget.SetOutlineColor(rgba[0], rgba[1], rgba[2])
    widget.SetOrientationMarker(axesActor)
    widget.SetInteractor(iren)
    widget.SetViewport(0.0, 0.0, 0.4, 0.4)
    widget.SetEnabled(1)
    widget.InteractiveOn()

    iren.Initialize()
    iren.Start()


def get_program_parameters():
    import argparse

    description = "DICOM Directory including `*.dcm` files"
    epilogue = """
    Derived from VTK/Examples/Cxx/Medical1.cxx
    This example reads a volume dataset, extracts an isosurface that
     represents the skin and displays it.
    """
    parser = argparse.ArgumentParser(
        description=description,
        epilog=epilogue,
        formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    parser.add_argument("dicom_dir", help="patients/series/")
    args = parser.parse_args()
    return args.dicom_dir


if __name__ == "__main__":
    main()

## Dicom3dVolume.py
"""
DICOM files --> vtk array --> vtk 3D visualization

1. read a series of DICOM files by vtkDICOMImageReader(which will do HU)
2. process vtk 3D array with volume
3. render vtk 3D array
"""

import vtkmodules.all as vtk


def main():
    dicom_dir = get_program_parameters()

    reader = vtk.vtkDICOMImageReader()
    reader.SetDirectoryName(dicom_dir)
    reader.Update()

    # The volume will be displayed by ray-cast alpha compositing.
    # A ray-cast mapper is needed to do the ray-casting.
    volume_mapper = vtk.vtkFixedPointVolumeRayCastMapper()
    volume_mapper.SetInputConnection(reader.GetOutputPort())

    # The color transfer function maps voxel intensities to colors.
    # It is modality-specific, and often anatomy-specific as well.
    # The goal is to one color for flesh (between 500 and 1000)
    # and another color for bone (1150 and over).
    volume_color = vtk.vtkColorTransferFunction()
    volume_color.AddRGBPoint(0, 0.0, 0.0, 0.0)
    volume_color.AddRGBPoint(500, 240.0 / 255.0, 184.0 / 255.0, 160.0 / 255.0)
    volume_color.AddRGBPoint(1000, 240.0 / 255.0, 184.0 / 255.0, 160.0 / 255.0)
    volume_color.AddRGBPoint(1150, 1.0, 1.0, 240.0 / 255.0)  # Ivory

    # The opacity transfer function is used to control the opacity
    # of different tissue types.
    volume_scalar_opacity = vtk.vtkPiecewiseFunction()
    volume_scalar_opacity.AddPoint(0, 0.00)
    volume_scalar_opacity.AddPoint(500, 0.15)
    volume_scalar_opacity.AddPoint(1000, 0.15)
    volume_scalar_opacity.AddPoint(1150, 0.85)

    # The gradient opacity function is used to decrease the opacity
    # in the 'flat' regions of the volume while maintaining the opacity
    # at the boundaries between tissue types.  The gradient is measured
    # as the amount by which the intensity changes over unit distance.
    # For most medical data, the unit distance is 1mm.
    volume_gradient_opacity = vtk.vtkPiecewiseFunction()
    volume_gradient_opacity.AddPoint(0, 0.0)
    volume_gradient_opacity.AddPoint(90, 0.5)
    volume_gradient_opacity.AddPoint(100, 1.0)

    # The VolumeProperty attaches the color and opacity functions to the
    # volume, and sets other volume properties.  The interpolation should
    # be set to linear to do a high-quality rendering.  The ShadeOn option
    # turns on directional lighting, which will usually enhance the
    # appearance of the volume and make it look more '3D'.  However,
    # the quality of the shading depends on how accurately the gradient
    # of the volume can be calculated, and for noisy data the gradient
    # estimation will be very poor.  The impact of the shading can be
    # decreased by increasing the Ambient coefficient while decreasing
    # the Diffuse and Specular coefficient.  To increase the impact
    # of shading, decrease the Ambient and increase the Diffuse and Specular.
    volume_property = vtk.vtkVolumeProperty()
    volume_property.SetColor(volume_color)
    volume_property.SetScalarOpacity(volume_scalar_opacity)
    volume_property.SetGradientOpacity(volume_gradient_opacity)
    volume_property.SetInterpolationTypeToLinear()
    volume_property.ShadeOn()
    volume_property.SetAmbient(0.4)
    volume_property.SetDiffuse(0.6)
    volume_property.SetSpecular(0.2)

    # The vtkVolume is a vtkProp3D (like a vtkActor) and controls the position
    # and orientation of the volume in world coordinates.
    volume = vtk.vtkVolume()
    volume.SetMapper(volume_mapper)
    volume.SetProperty(volume_property)

    # Setup render
    renderer = vtk.vtkRenderer()
    renderer.AddViewProp(volume)
    renderer.SetBackground(1.0, 1.0, 1.0)
    renderer.ResetCamera()

    # Set up an initial view of the volume.  The focal point will be the
    # center of the volume, and the camera position will be 400mm to the
    # patient's left (which is our right).
    camera = renderer.GetActiveCamera()
    c = volume.GetCenter()
    # camera.SetViewUp(0, 0, -1)
    # camera.SetPosition(c[0], c[1] - 400, c[2])
    # camera.SetFocalPoint(c[0], c[1], c[2])
    camera.Azimuth(0.0)
    camera.Elevation(80.0)

    # Setup render window
    renWin = vtk.vtkRenderWindow()
    renWin.AddRenderer(renderer)
    renWin.SetSize(640, 480)

    renWin.Render()

    # Setup render window interactor
    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)

    iren.Initialize()
    iren.Start()


def get_program_parameters():
    import argparse

    description = "DICOM Directory including `*.dcm` files"
    epilogue = """
    This example reads a series of DICOM files, and render 3d volume.
    """
    parser = argparse.ArgumentParser(
        description=description,
        epilog=epilogue,
        formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    parser.add_argument(
        "dicom_dir", help="a DICOM files directory, like: patients/series/"
    )
    args = parser.parse_args()
    return args.dicom_dir


if __name__ == "__main__":
    main()


## Dicom3dVolumeWithNumpy.py
"""
DICOM files --> vtk array --> vtk 3D visualization

1. read a series of DICOM files by vtkDICOMImageReader(which will do HU)
2. process vtk 3D array with volume
3. render vtk 3D array
"""

import vtkmodules.all as vtk


def main():
    dicom_dir = get_program_parameters()

    reader = vtk.vtkDICOMImageReader()
    reader.SetDirectoryName(dicom_dir)
    reader.Update()

    # The volume will be displayed by ray-cast alpha compositing.
    # A ray-cast mapper is needed to do the ray-casting.
    volume_mapper = vtk.vtkFixedPointVolumeRayCastMapper()
    volume_mapper.SetInputConnection(reader.GetOutputPort())

    # The color transfer function maps voxel intensities to colors.
    # It is modality-specific, and often anatomy-specific as well.
    # The goal is to one color for flesh (between 500 and 1000)
    # and another color for bone (1150 and over).
    volume_color = vtk.vtkColorTransferFunction()
    volume_color.AddRGBPoint(0, 0.0, 0.0, 0.0)
    volume_color.AddRGBPoint(500, 240.0 / 255.0, 184.0 / 255.0, 160.0 / 255.0)
    volume_color.AddRGBPoint(1000, 240.0 / 255.0, 184.0 / 255.0, 160.0 / 255.0)
    volume_color.AddRGBPoint(1150, 1.0, 1.0, 240.0 / 255.0)  # Ivory

    # The opacity transfer function is used to control the opacity
    # of different tissue types.
    volume_scalar_opacity = vtk.vtkPiecewiseFunction()
    volume_scalar_opacity.AddPoint(0, 0.00)
    volume_scalar_opacity.AddPoint(500, 0.15)
    volume_scalar_opacity.AddPoint(1000, 0.15)
    volume_scalar_opacity.AddPoint(1150, 0.85)

    # The gradient opacity function is used to decrease the opacity
    # in the 'flat' regions of the volume while maintaining the opacity
    # at the boundaries between tissue types.  The gradient is measured
    # as the amount by which the intensity changes over unit distance.
    # For most medical data, the unit distance is 1mm.
    volume_gradient_opacity = vtk.vtkPiecewiseFunction()
    volume_gradient_opacity.AddPoint(0, 0.0)
    volume_gradient_opacity.AddPoint(90, 0.5)
    volume_gradient_opacity.AddPoint(100, 1.0)

    # The VolumeProperty attaches the color and opacity functions to the
    # volume, and sets other volume properties.  The interpolation should
    # be set to linear to do a high-quality rendering.  The ShadeOn option
    # turns on directional lighting, which will usually enhance the
    # appearance of the volume and make it look more '3D'.  However,
    # the quality of the shading depends on how accurately the gradient
    # of the volume can be calculated, and for noisy data the gradient
    # estimation will be very poor.  The impact of the shading can be
    # decreased by increasing the Ambient coefficient while decreasing
    # the Diffuse and Specular coefficient.  To increase the impact
    # of shading, decrease the Ambient and increase the Diffuse and Specular.
    volume_property = vtk.vtkVolumeProperty()
    volume_property.SetColor(volume_color)
    volume_property.SetScalarOpacity(volume_scalar_opacity)
    volume_property.SetGradientOpacity(volume_gradient_opacity)
    volume_property.SetInterpolationTypeToLinear()
    volume_property.ShadeOn()
    volume_property.SetAmbient(0.4)
    volume_property.SetDiffuse(0.6)
    volume_property.SetSpecular(0.2)

    # The vtkVolume is a vtkProp3D (like a vtkActor) and controls the position
    # and orientation of the volume in world coordinates.
    volume = vtk.vtkVolume()
    volume.SetMapper(volume_mapper)
    volume.SetProperty(volume_property)

    # Setup render
    renderer = vtk.vtkRenderer()
    # renderer.AddActor(axesActor)
    renderer.AddViewProp(volume)
    renderer.SetBackground(1.0, 1.0, 1.0)
    renderer.ResetCamera()

    # Setup render window
    renWin = vtk.vtkRenderWindow()
    renWin.AddRenderer(renderer)
    renWin.SetSize(800, 800)

    renWin.Render()

    # Setup render window interactor
    iren = vtk.vtkRenderWindowInteractor()
    iren.SetRenderWindow(renWin)

    iren.Initialize()
    iren.Start()


def get_program_parameters():
    import argparse

    description = "DICOM Directory including `*.dcm` files"
    epilogue = """
    Derived from VTK/Examples/Cxx/Medical1.cxx
    This example reads a volume dataset, extracts an isosurface that
     represents the skin and displays it.
    """
    parser = argparse.ArgumentParser(
        description=description,
        epilog=epilogue,
        formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    parser.add_argument("dicom_dir", help="patients/series/")
    args = parser.parse_args()
    return args.dicom_dir


if __name__ == "__main__":
    main()

## DicomReadPydicom.py
"""
1. read DICOM files via pydicom
2. convert vtk 3D array to numpy 3D array
"""

import numpy as np
import vtk
from utils import get_pixels_hu, load_scan


def main():
    dicom_dir = get_program_parameters()

    # Load the DICOM files
    slices = load_scan(dicom_dir)

    # pixel aspects, assuming all slices are the same
    ps = slices[0].PixelSpacing
    ss = slices[9].SliceThickness
    ax_aspect = ps[1] / ps[0]
    sag_aspect = ps[1] / ss
    cor_aspect = ss / ps[0]

    # NOTE: `SliceThickness` from DICOM tag may be different with `slice_thickness` here
    try:
        slice_thickness = np.abs(
            slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2]
        )
    except:
        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)

    print(f"PixelSpacing[{ps}] SliceThickness[{ss}] [{slice_thickness}]")

    # Convert raw data into Houndsfeld units
    img3d = get_pixels_hu(slices)

    print(f"img3d [{img3d.shape}] [{img3d.dtype}] [{img3d.nbytes}] [{np.mean(img3d)}] ")


def get_program_parameters():
    import argparse

    description = "DICOM Directory including `*.dcm` files"
    epilogue = """
    Derived from VTK/Examples/Cxx/Medical1.cxx
    This example reads a volume dataset, extracts an isosurface that
     represents the skin and displays it.
    """
    parser = argparse.ArgumentParser(
        description=description,
        epilog=epilogue,
        formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    parser.add_argument("dicom_dir", help="patients/series/")
    args = parser.parse_args()
    return args.dicom_dir


if __name__ == "__main__":
    main()

## DicomReadVtk.py
"""
1. read DICOM files via VTK `vtkDICOMImageReader` class
2. convert vtk 3D array to numpy 3D array
"""
from vtk.util import numpy_support
import numpy as np
import vtk


def main():
    dicom_dir = get_program_parameters()

    reader = vtk.vtkDICOMImageReader()
    reader.SetDirectoryName(dicom_dir)
    reader.Update()

    # Load dimensions using `GetDataExtent`
    _extent = reader.GetDataExtent()
    ConstPixelDims = [
        _extent[1] - _extent[0] + 1,
        _extent[3] - _extent[2] + 1,
        _extent[5] - _extent[4] + 1,
    ]  # x,y,z
    w = ConstPixelDims[0]  #
    h = ConstPixelDims[1]  #
    num = ConstPixelDims[2]  #

    # Load spacing values
    ConstPixelSpacing = reader.GetPixelSpacing()  # x, y, z

    print(
        f"PixelSpacing[{ConstPixelSpacing[0]}, {ConstPixelSpacing[1]}] SliceThickness[{ConstPixelSpacing[2]}]"
    )

    # Get the 'vtkImageData' object from the reader
    imageData = reader.GetOutput()
    # Get the 'vtkPointData' object from the 'vtkImageData' object
    pointData = imageData.GetPointData()
    # Ensure that only one array exists within the 'vtkPointData' object
    assert pointData.GetNumberOfArrays() == 1
    # Get the `vtkArray` (or whatever derived type) which is needed for the `numpy_support.vtk_to_numpy` function
    arrayData = pointData.GetArray(0)

    # Convert the `vtkArray` to a NumPy array
    ArrayDicom = numpy_support.vtk_to_numpy(arrayData)
    # Reshape the NumPy array to 3D using 'ConstPixelDims' as a 'shape'
    ArrayDicom = ArrayDicom.reshape((num, h, w), order="F")

    print(
        f"ArrayDicom [{ArrayDicom.shape}] [{ArrayDicom.dtype}] [{ArrayDicom.nbytes}] [{np.mean(ArrayDicom)}] "
    )


def get_program_parameters():
    import argparse

    description = "DICOM Directory including `*.dcm` files"
    epilogue = """
    Derived from VTK/Examples/Cxx/Medical1.cxx
    This example reads a volume dataset, extracts an isosurface that
     represents the skin and displays it.
    """
    parser = argparse.ArgumentParser(
        description=description,
        epilog=epilogue,
        formatter_class=argparse.RawDescriptionHelpFormatter,
    )
    parser.add_argument("dicom_dir", help="patients/series/")
    args = parser.parse_args()
    return args.dicom_dir


if __name__ == "__main__":
    main()

## utils.py
import vtk
import numpy as np


def CreateLut():
    colors = vtk.vtkNamedColors()

    colorLut = vtk.vtkLookupTable()
    colorLut.SetNumberOfColors(17)
    colorLut.SetTableRange(0, 16)
    colorLut.Build()

    colorLut.SetTableValue(0, 0, 0, 0, 0)
    colorLut.SetTableValue(1, colors.GetColor4d("salmon"))  # blood
    colorLut.SetTableValue(2, colors.GetColor4d("beige"))  # brain
    colorLut.SetTableValue(3, colors.GetColor4d("orange"))  # duodenum
    colorLut.SetTableValue(4, colors.GetColor4d("misty_rose"))  # eye_retina
    colorLut.SetTableValue(5, colors.GetColor4d("white"))  # eye_white
    colorLut.SetTableValue(6, colors.GetColor4d("tomato"))  # heart
    colorLut.SetTableValue(7, colors.GetColor4d("raspberry"))  # ileum
    colorLut.SetTableValue(8, colors.GetColor4d("banana"))  # kidney
    colorLut.SetTableValue(9, colors.GetColor4d("peru"))  # l_intestine
    colorLut.SetTableValue(10, colors.GetColor4d("pink"))  # liver
    colorLut.SetTableValue(11, colors.GetColor4d("powder_blue"))  # lung
    colorLut.SetTableValue(12, colors.GetColor4d("carrot"))  # nerve
    colorLut.SetTableValue(13, colors.GetColor4d("wheat"))  # skeleton
    colorLut.SetTableValue(14, colors.GetColor4d("violet"))  # spleen
    colorLut.SetTableValue(15, colors.GetColor4d("plum"))  # stomach

    return colorLut


def CreateTissueMap():
    tissueMap = dict()
    tissueMap["blood"] = 1
    tissueMap["brain"] = 2
    tissueMap["duodenum"] = 3
    tissueMap["eyeRetina"] = 4
    tissueMap["eyeWhite"] = 5
    tissueMap["heart"] = 6
    tissueMap["ileum"] = 7
    tissueMap["kidney"] = 8
    tissueMap["intestine"] = 9
    tissueMap["liver"] = 10
    tissueMap["lung"] = 11
    tissueMap["nerve"] = 12
    tissueMap["skeleton"] = 13
    tissueMap["spleen"] = 14
    tissueMap["stomach"] = 15

    return tissueMap


tissueMap = CreateTissueMap()
colorLut = CreateLut()


# Generate a 3D mesh with marching cubes algorithm
def CreateTissue(
    image_producer: vtk.vtkAlgorithm, ThrIn, ThrOut, color="skeleton", isoValue=127.5
):
    selectTissue = vtk.vtkImageThreshold()
    selectTissue.ThresholdBetween(ThrIn, ThrOut)
    selectTissue.ReplaceInOn()
    selectTissue.SetInValue(255)
    selectTissue.ReplaceOutOn()
    selectTissue.SetOutValue(0)
    # selectTissue.SetInputData(imageData)
    selectTissue.SetInputConnection(image_producer.GetOutputPort())
    # selectTissue.Update()

    gaussianRadius = 5
    gaussianStandardDeviation = 2.0
    gaussian = vtk.vtkImageGaussianSmooth()
    gaussian.SetStandardDeviations(
        gaussianStandardDeviation, gaussianStandardDeviation, gaussianStandardDeviation
    )
    gaussian.SetRadiusFactors(gaussianRadius, gaussianRadius, gaussianRadius)
    gaussian.SetInputConnection(selectTissue.GetOutputPort())

    # isoValue = 127.5
    mcubes = vtk.vtkMarchingCubes()
    mcubes.SetInputConnection(gaussian.GetOutputPort())
    mcubes.ComputeScalarsOff()
    mcubes.ComputeGradientsOff()
    mcubes.ComputeNormalsOff()
    mcubes.SetValue(0, isoValue)

    smoothingIterations = 5
    passBand = 0.001
    featureAngle = 60.0
    smoother = vtk.vtkWindowedSincPolyDataFilter()
    smoother.SetInputConnection(mcubes.GetOutputPort())
    smoother.SetNumberOfIterations(smoothingIterations)
    smoother.BoundarySmoothingOff()
    smoother.FeatureEdgeSmoothingOff()
    smoother.SetFeatureAngle(featureAngle)
    smoother.SetPassBand(passBand)
    smoother.NonManifoldSmoothingOn()
    smoother.NormalizeCoordinatesOn()
    smoother.Update()

    normals = vtk.vtkPolyDataNormals()
    normals.SetInputConnection(smoother.GetOutputPort())
    normals.SetFeatureAngle(featureAngle)

    stripper = vtk.vtkStripper()
    stripper.SetInputConnection(normals.GetOutputPort())

    mapper = vtk.vtkPolyDataMapper()
    mapper.SetInputConnection(stripper.GetOutputPort())

    actor = vtk.vtkActor()
    actor.SetMapper(mapper)
    actor.GetProperty().SetColor(colorLut.GetTableValue(tissueMap[color])[:3])
    actor.GetProperty().SetSpecular(0.5)
    actor.GetProperty().SetSpecularPower(10)

    return actor


import pathlib
import pydicom
from typing import List


def load_scan(dicom_dir: str) -> List[pydicom.FileDataset]:
    slices = [pydicom.dcmread(f) for f in pathlib.Path(dicom_dir).glob("*.dcm")]
    print(f"file count: {len(slices)}")

    # skip files with no SliceLocation or InstanceNumber(eg scout views)
    slices = list(filter(lambda s: hasattr(s, "SliceLocation"), slices))

    print(f"file count with SliceLocation: {len(slices)}")
    slices = sorted(slices, key=lambda s: s.SliceLocation)

    return slices


def get_pixels_hu(slices: List[pydicom.FileDataset]) -> np.ndarray:
    # s.pixel_array: 2D array with HW format
    # img3d: 3D array with NHW format
    img3d = np.stack([s.pixel_array for s in slices])

    # Convert to int16 (from sometimes int16),
    # should be possible as values should always be low enough (<32k)
    img3d = img3d.astype(np.int16)

    # Set outside-of-scan pixels to 1
    # The intercept is usually -1024, so air is approximately 0
    # img3d[img3d == -1000] = 0
    img3d[img3d == -2000] = 0

    # Convert to Hounsfield units (HU)

    # Operate on whole array
    intercept = slices[0].RescaleIntercept
    slope = slices[0].RescaleSlope

    if slope != 1:
        img3d = slope * img3d.astype(np.float64)
        img3d = img3d.astype(np.int16)

    img3d += np.int16(intercept)

    # # Operate on individual array
    # for n in range(len(slices)):
    #     intercept = slices[n].RescaleIntercept
    #     slope = slices[n].RescaleSlope

    #     if slope != 1:
    #         img3d[n] = slope * img3d[n].astype(np.float64)
    #         img3d[n] = img3d[n].astype(np.int16)

    #     img3d[n] += np.int16(intercept)

    return np.array(img3d, dtype=np.int16)


def get_pixels(scans: List[pydicom.FileDataset]) -> np.ndarray:
    # s.pixel_array: 2D array with HW format
    # img3d: 3D array with NHW format
    img3d = np.stack([s.pixel_array for s in scans])

    return img3d


def plot_stack_sample(stack, rows=6, cols=6, start_with=10, show_every=3):
    import matplotlib.pyplot as plt

    fig, ax = plt.subplots(rows, cols, figsize=[12, 12])
    for i in range(rows * cols):
        ind = start_with + i * show_every
        ax[int(i / rows), int(i % rows)].set_title("slice %d" % ind)
        ax[int(i / rows), int(i % rows)].imshow(stack[ind], cmap="gray")
        ax[int(i / rows), int(i % rows)].axis("off")
    plt.show()


def resample(image, scan, new_spacing=[1, 1, 1]):
    import scipy.ndimage

    # Determine current pixel spacing
    spacing = np.array(
        [
            float(scan[0].SliceThickness),
            scan[0].PixelSpacing[0],
            scan[0].PixelSpacing[1],
        ]
    )

    resize_factor = spacing / new_spacing
    new_real_shape = image.shape * resize_factor
    new_shape = np.round(new_real_shape)
    real_resize_factor = new_shape / image.shape
    new_spacing = spacing / real_resize_factor

    image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)

    return image, new_spacing
	"""
	DICOM files --> vtk array --> vtk 3D visualization

	1. read a series of DICOM files by vtkDICOMImageReader(which will do HU and don't resample)
	2. process vtk 3D array with marching cubes algorithm
	3. render vtk 3D array
	"""
	import vtk
	from utils import CreateTissue


	def main():
	dicom_dir = get_program_parameters()

	reader = vtk.vtkDICOMImageReader()
	reader.SetDirectoryName(dicom_dir)
	reader.Update()

	# Setup render
	renderer = vtk.vtkRenderer()
	renderer.AddActor(CreateTissue(reader, -900, -400, "lung"))
	renderer.AddActor(CreateTissue(reader, 0, 120, "blood"))
	renderer.AddActor(CreateTissue(reader, 100, 1000, "skeleton"))
	renderer.SetBackground(1.0, 1.0, 1.0)
	renderer.ResetCamera()

	# Setup render window
	renWin = vtk.vtkRenderWindow()
	renWin.AddRenderer(renderer)
	renWin.SetSize(800, 800)
	renWin.Render()

	# Setup render window interactor
	iren = vtk.vtkRenderWindowInteractor()
	iren.SetRenderWindow(renWin)

	# Setup coordinate axes helper
	axesActor = vtk.vtkAxesActor()
	widget = vtk.vtkOrientationMarkerWidget()
	rgba = [0] * 4
	colors = vtk.vtkNamedColors()
	colors.GetColor("Carrot", rgba)
	widget.SetOutlineColor(rgba[0], rgba[1], rgba[2])
	widget.SetOrientationMarker(axesActor)
	widget.SetInteractor(iren)
	widget.SetViewport(0.0, 0.0, 0.4, 0.4)
	widget.SetEnabled(1)
	widget.InteractiveOn()

	iren.Initialize()
	iren.Start()


	def get_program_parameters():
	import argparse

	description = "DICOM Directory including `*.dcm` files"
	epilogue = """
	Derived from VTK/Examples/Cxx/Medical1.cxx
	This example reads a volume dataset, extracts an isosurface that
	represents the skin and displays it.
	"""
	parser = argparse.ArgumentParser(
	description=description,
	epilog=epilogue,
	formatter_class=argparse.RawDescriptionHelpFormatter,
	)
	parser.add_argument("dicom_dir", help="patients/series/")
	args = parser.parse_args()
	return args.dicom_dir


	if __name__ == "__main__":
	main()
	"""
	DICOM files --> 3D numpy array --> vtk array --> vtk 3D visualization

	1. read a series of DICOM files by pydicom
	2. create a 3D numpy array(NHW) by stacking a series of 2D DICOM image data(HW) and preprocessing(hu, resample,..etc) below
	1. HU
	2. resample
	3. convert the 3D numpy to vtk 3D array(vtkImageData)
	4. bridge vtk 3D array `vtkImageData` to `vtkAlgorigthmOutput` by `vtkTrivialProducer`, then process with marching cubes algorithm
	5. render vtk 3D array
	"""


	import numpy as np

	import vtk
	from vtk.util import numpy_support
	from utils import CreateTissue, get_pixels_hu, load_scan, resample


	def main():
	dicom_dir = get_program_parameters()

	# Load the DICOM files
	slices = load_scan(dicom_dir)

	# pixel aspects, assuming all slices are the same
	ps = slices[0].PixelSpacing
	ss = slices[0].SliceThickness
	ax_aspect = ps[1] / ps[0]
	sag_aspect = ps[1] / ss
	cor_aspect = ss / ps[0]

	print(
	f"PixelSpacing[{ps}] SliceThickness[{ss}] AxialAspect[{ax_aspect}] SagittalAspect[{sag_aspect}] CoronalAspect[{cor_aspect}]"
	)

	# Convert raw data into Houndsfeld units
	img3d = get_pixels_hu(slices)

	print(f"img3d [{img3d.shape}] [{img3d.dtype}] [{img3d.nbytes}] [{np.mean(img3d)}] ")

	# Resample
	img3d_after_resamp, spacing = resample(img3d, slices)

	print(
	f"img3d_after_resamp [{img3d_after_resamp.shape}] [{img3d_after_resamp.dtype}] [{img3d_after_resamp.nbytes}] [{np.mean(img3d_after_resamp)}] "
	)

	# Prepare 3D visualization
	print(f"3D visualization...")

	# Convert numpy array to vtk array
	img3d_64f = img3d_after_resamp.astype(np.float64)
	num, h, w = img3d_64f.shape
	print(f"img3d_64f [{img3d_64f.shape}] [{np.max(img3d_64f[1])}]")

	imageData = vtk.vtkImageData()
	depthArray = numpy_support.numpy_to_vtk(
	num_array=img3d_64f.ravel(), deep=True, array_type=vtk.VTK_FLOAT
	)
	imageData.SetDimensions((w, h, num)) # x,y,z coordination system
	imageData.SetSpacing([1, 1, 1])
	# imageData.SetSpacing([ps[0], ps[1], ss])
	imageData.SetOrigin([0, 0, 0])
	imageData.GetPointData().SetScalars(depthArray)

	image_producer: vtk.vtkAlgorithm = vtk.vtkTrivialProducer() # data provider
	image_producer.SetOutput(imageData)

	# Rendering

	# Setup render
	renderer = vtk.vtkRenderer()
	renderer.AddActor(CreateTissue(image_producer, -900, -400, "lung"))
	renderer.AddActor(CreateTissue(image_producer, 0, 120, "blood"))
	renderer.AddActor(CreateTissue(image_producer, 100, 2000, "skeleton"))
	renderer.SetBackground(1.0, 1.0, 1.0)
	renderer.ResetCamera()

	# Setup render window
	renWin = vtk.vtkRenderWindow()
	renWin.AddRenderer(renderer)
	renWin.SetSize(800, 800)
	renWin.Render()

	# Setup render window interaction
	iren = vtk.vtkRenderWindowInteractor()
	iren.SetRenderWindow(renWin)

	# Setup coordinate axes helper
	axesActor = vtk.vtkAxesActor()
	widget = vtk.vtkOrientationMarkerWidget()
	rgba = [0] * 4
	colors = vtk.vtkNamedColors()
	colors.GetColor("Carrot", rgba)
	widget.SetOutlineColor(rgba[0], rgba[1], rgba[2])
	widget.SetOrientationMarker(axesActor)
	widget.SetInteractor(iren)
	widget.SetViewport(0.0, 0.0, 0.4, 0.4)
	widget.SetEnabled(1)
	widget.InteractiveOn()

	iren.Initialize()
	iren.Start()


	def get_program_parameters():
	import argparse

	description = "DICOM Directory including `*.dcm` files"
	epilogue = """
	Derived from VTK/Examples/Cxx/Medical1.cxx
	This example reads a volume dataset, extracts an isosurface that
	represents the skin and displays it.
	"""
	parser = argparse.ArgumentParser(
	description=description,
	epilog=epilogue,
	formatter_class=argparse.RawDescriptionHelpFormatter,
	)
	parser.add_argument("dicom_dir", help="patients/series/")
	args = parser.parse_args()
	return args.dicom_dir


	if __name__ == "__main__":
	main()
	"""
	DICOM files --> vtk array --> vtk 3D visualization

	1. read a series of DICOM files by vtkDICOMImageReader(which will do HU)
	2. process vtk 3D array with volume
	3. render vtk 3D array
	"""

	import vtkmodules.all as vtk


	def main():
	dicom_dir = get_program_parameters()

	reader = vtk.vtkDICOMImageReader()
	reader.SetDirectoryName(dicom_dir)
	reader.Update()

	# The volume will be displayed by ray-cast alpha compositing.
	# A ray-cast mapper is needed to do the ray-casting.
	volume_mapper = vtk.vtkFixedPointVolumeRayCastMapper()
	volume_mapper.SetInputConnection(reader.GetOutputPort())

	# The color transfer function maps voxel intensities to colors.
	# It is modality-specific, and often anatomy-specific as well.
	# The goal is to one color for flesh (between 500 and 1000)
	# and another color for bone (1150 and over).
	volume_color = vtk.vtkColorTransferFunction()
	volume_color.AddRGBPoint(0, 0.0, 0.0, 0.0)
	volume_color.AddRGBPoint(500, 240.0 / 255.0, 184.0 / 255.0, 160.0 / 255.0)
	volume_color.AddRGBPoint(1000, 240.0 / 255.0, 184.0 / 255.0, 160.0 / 255.0)
	volume_color.AddRGBPoint(1150, 1.0, 1.0, 240.0 / 255.0) # Ivory

	# The opacity transfer function is used to control the opacity
	# of different tissue types.
	volume_scalar_opacity = vtk.vtkPiecewiseFunction()
	volume_scalar_opacity.AddPoint(0, 0.00)
	volume_scalar_opacity.AddPoint(500, 0.15)
	volume_scalar_opacity.AddPoint(1000, 0.15)
	volume_scalar_opacity.AddPoint(1150, 0.85)

	# The gradient opacity function is used to decrease the opacity
	# in the 'flat' regions of the volume while maintaining the opacity
	# at the boundaries between tissue types. The gradient is measured
	# as the amount by which the intensity changes over unit distance.
	# For most medical data, the unit distance is 1mm.
	volume_gradient_opacity = vtk.vtkPiecewiseFunction()
	volume_gradient_opacity.AddPoint(0, 0.0)
	volume_gradient_opacity.AddPoint(90, 0.5)
	volume_gradient_opacity.AddPoint(100, 1.0)

	# The VolumeProperty attaches the color and opacity functions to the
	# volume, and sets other volume properties. The interpolation should
	# be set to linear to do a high-quality rendering. The ShadeOn option
	# turns on directional lighting, which will usually enhance the
	# appearance of the volume and make it look more '3D'. However,
	# the quality of the shading depends on how accurately the gradient
	# of the volume can be calculated, and for noisy data the gradient
	# estimation will be very poor. The impact of the shading can be
	# decreased by increasing the Ambient coefficient while decreasing
	# the Diffuse and Specular coefficient. To increase the impact
	# of shading, decrease the Ambient and increase the Diffuse and Specular.
	volume_property = vtk.vtkVolumeProperty()
	volume_property.SetColor(volume_color)
	volume_property.SetScalarOpacity(volume_scalar_opacity)
	volume_property.SetGradientOpacity(volume_gradient_opacity)
	volume_property.SetInterpolationTypeToLinear()
	volume_property.ShadeOn()
	volume_property.SetAmbient(0.4)
	volume_property.SetDiffuse(0.6)
	volume_property.SetSpecular(0.2)

	# The vtkVolume is a vtkProp3D (like a vtkActor) and controls the position
	# and orientation of the volume in world coordinates.
	volume = vtk.vtkVolume()
	volume.SetMapper(volume_mapper)
	volume.SetProperty(volume_property)

	# Setup render
	renderer = vtk.vtkRenderer()
	renderer.AddViewProp(volume)
	renderer.SetBackground(1.0, 1.0, 1.0)
	renderer.ResetCamera()

	# Set up an initial view of the volume. The focal point will be the
	# center of the volume, and the camera position will be 400mm to the
	# patient's left (which is our right).
	camera = renderer.GetActiveCamera()
	c = volume.GetCenter()
	# camera.SetViewUp(0, 0, -1)
	# camera.SetPosition(c[0], c[1] - 400, c[2])
	# camera.SetFocalPoint(c[0], c[1], c[2])
	camera.Azimuth(0.0)
	camera.Elevation(80.0)

	# Setup render window
	renWin = vtk.vtkRenderWindow()
	renWin.AddRenderer(renderer)
	renWin.SetSize(640, 480)

	renWin.Render()

	# Setup render window interactor
	iren = vtk.vtkRenderWindowInteractor()
	iren.SetRenderWindow(renWin)

	iren.Initialize()
	iren.Start()


	def get_program_parameters():
	import argparse

	description = "DICOM Directory including `*.dcm` files"
	epilogue = """
	This example reads a series of DICOM files, and render 3d volume.
	"""
	parser = argparse.ArgumentParser(
	description=description,
	epilog=epilogue,
	formatter_class=argparse.RawDescriptionHelpFormatter,
	)
	parser.add_argument(
	"dicom_dir", help="a DICOM files directory, like: patients/series/"
	)
	args = parser.parse_args()
	return args.dicom_dir


	if __name__ == "__main__":
	main()
	"""
	1. read DICOM files via pydicom
	2. convert vtk 3D array to numpy 3D array
	"""

	import numpy as np
	import vtk
	from utils import get_pixels_hu, load_scan


	def main():
	dicom_dir = get_program_parameters()

	# Load the DICOM files
	slices = load_scan(dicom_dir)

	# pixel aspects, assuming all slices are the same
	ps = slices[0].PixelSpacing
	ss = slices[9].SliceThickness
	ax_aspect = ps[1] / ps[0]
	sag_aspect = ps[1] / ss
	cor_aspect = ss / ps[0]

	# NOTE: `SliceThickness` from DICOM tag may be different with `slice_thickness` here
	try:
	slice_thickness = np.abs(
	slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2]
	)
	except:
	slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)

	print(f"PixelSpacing[{ps}] SliceThickness[{ss}] [{slice_thickness}]")

	# Convert raw data into Houndsfeld units
	img3d = get_pixels_hu(slices)

	print(f"img3d [{img3d.shape}] [{img3d.dtype}] [{img3d.nbytes}] [{np.mean(img3d)}] ")


	def get_program_parameters():
	import argparse

	description = "DICOM Directory including `*.dcm` files"
	epilogue = """
	Derived from VTK/Examples/Cxx/Medical1.cxx
	This example reads a volume dataset, extracts an isosurface that
	represents the skin and displays it.
	"""
	parser = argparse.ArgumentParser(
	description=description,
	epilog=epilogue,
	formatter_class=argparse.RawDescriptionHelpFormatter,
	)
	parser.add_argument("dicom_dir", help="patients/series/")
	args = parser.parse_args()
	return args.dicom_dir


	if __name__ == "__main__":
	main()
	"""
	1. read DICOM files via VTK `vtkDICOMImageReader` class
	2. convert vtk 3D array to numpy 3D array
	"""
	from vtk.util import numpy_support
	import numpy as np
	import vtk


	def main():
	dicom_dir = get_program_parameters()

	reader = vtk.vtkDICOMImageReader()
	reader.SetDirectoryName(dicom_dir)
	reader.Update()

	# Load dimensions using `GetDataExtent`
	_extent = reader.GetDataExtent()
	ConstPixelDims = [
	_extent[1] - _extent[0] + 1,
	_extent[3] - _extent[2] + 1,
	_extent[5] - _extent[4] + 1,
	] # x,y,z
	w = ConstPixelDims[0] #
	h = ConstPixelDims[1] #
	num = ConstPixelDims[2] #

	# Load spacing values
	ConstPixelSpacing = reader.GetPixelSpacing() # x, y, z

	print(
	f"PixelSpacing[{ConstPixelSpacing[0]}, {ConstPixelSpacing[1]}] SliceThickness[{ConstPixelSpacing[2]}]"
	)

	# Get the 'vtkImageData' object from the reader
	imageData = reader.GetOutput()
	# Get the 'vtkPointData' object from the 'vtkImageData' object
	pointData = imageData.GetPointData()
	# Ensure that only one array exists within the 'vtkPointData' object
	assert pointData.GetNumberOfArrays() == 1
	# Get the `vtkArray` (or whatever derived type) which is needed for the `numpy_support.vtk_to_numpy` function
	arrayData = pointData.GetArray(0)

	# Convert the `vtkArray` to a NumPy array
	ArrayDicom = numpy_support.vtk_to_numpy(arrayData)
	# Reshape the NumPy array to 3D using 'ConstPixelDims' as a 'shape'
	ArrayDicom = ArrayDicom.reshape((num, h, w), order="F")

	print(
	f"ArrayDicom [{ArrayDicom.shape}] [{ArrayDicom.dtype}] [{ArrayDicom.nbytes}] [{np.mean(ArrayDicom)}] "
	)


	def get_program_parameters():
	import argparse

	description = "DICOM Directory including `*.dcm` files"
	epilogue = """
	Derived from VTK/Examples/Cxx/Medical1.cxx
	This example reads a volume dataset, extracts an isosurface that
	represents the skin and displays it.
	"""
	parser = argparse.ArgumentParser(
	description=description,
	epilog=epilogue,
	formatter_class=argparse.RawDescriptionHelpFormatter,
	)
	parser.add_argument("dicom_dir", help="patients/series/")
	args = parser.parse_args()
	return args.dicom_dir


	if __name__ == "__main__":
	main()
	import vtk
	import numpy as np


	def CreateLut():
	colors = vtk.vtkNamedColors()

	colorLut = vtk.vtkLookupTable()
	colorLut.SetNumberOfColors(17)
	colorLut.SetTableRange(0, 16)
	colorLut.Build()

	colorLut.SetTableValue(0, 0, 0, 0, 0)
	colorLut.SetTableValue(1, colors.GetColor4d("salmon")) # blood
	colorLut.SetTableValue(2, colors.GetColor4d("beige")) # brain
	colorLut.SetTableValue(3, colors.GetColor4d("orange")) # duodenum
	colorLut.SetTableValue(4, colors.GetColor4d("misty_rose")) # eye_retina
	colorLut.SetTableValue(5, colors.GetColor4d("white")) # eye_white
	colorLut.SetTableValue(6, colors.GetColor4d("tomato")) # heart
	colorLut.SetTableValue(7, colors.GetColor4d("raspberry")) # ileum
	colorLut.SetTableValue(8, colors.GetColor4d("banana")) # kidney
	colorLut.SetTableValue(9, colors.GetColor4d("peru")) # l_intestine
	colorLut.SetTableValue(10, colors.GetColor4d("pink")) # liver
	colorLut.SetTableValue(11, colors.GetColor4d("powder_blue")) # lung
	colorLut.SetTableValue(12, colors.GetColor4d("carrot")) # nerve
	colorLut.SetTableValue(13, colors.GetColor4d("wheat")) # skeleton
	colorLut.SetTableValue(14, colors.GetColor4d("violet")) # spleen
	colorLut.SetTableValue(15, colors.GetColor4d("plum")) # stomach

	return colorLut


	def CreateTissueMap():
	tissueMap = dict()
	tissueMap["blood"] = 1
	tissueMap["brain"] = 2
	tissueMap["duodenum"] = 3
	tissueMap["eyeRetina"] = 4
	tissueMap["eyeWhite"] = 5
	tissueMap["heart"] = 6
	tissueMap["ileum"] = 7
	tissueMap["kidney"] = 8
	tissueMap["intestine"] = 9
	tissueMap["liver"] = 10
	tissueMap["lung"] = 11
	tissueMap["nerve"] = 12
	tissueMap["skeleton"] = 13
	tissueMap["spleen"] = 14
	tissueMap["stomach"] = 15

	return tissueMap


	tissueMap = CreateTissueMap()
	colorLut = CreateLut()


	# Generate a 3D mesh with marching cubes algorithm
	def CreateTissue(
	image_producer: vtk.vtkAlgorithm, ThrIn, ThrOut, color="skeleton", isoValue=127.5
	):
	selectTissue = vtk.vtkImageThreshold()
	selectTissue.ThresholdBetween(ThrIn, ThrOut)
	selectTissue.ReplaceInOn()
	selectTissue.SetInValue(255)
	selectTissue.ReplaceOutOn()
	selectTissue.SetOutValue(0)
	# selectTissue.SetInputData(imageData)
	selectTissue.SetInputConnection(image_producer.GetOutputPort())
	# selectTissue.Update()

	gaussianRadius = 5
	gaussianStandardDeviation = 2.0
	gaussian = vtk.vtkImageGaussianSmooth()
	gaussian.SetStandardDeviations(
	gaussianStandardDeviation, gaussianStandardDeviation, gaussianStandardDeviation
	)
	gaussian.SetRadiusFactors(gaussianRadius, gaussianRadius, gaussianRadius)
	gaussian.SetInputConnection(selectTissue.GetOutputPort())

	# isoValue = 127.5
	mcubes = vtk.vtkMarchingCubes()
	mcubes.SetInputConnection(gaussian.GetOutputPort())
	mcubes.ComputeScalarsOff()
	mcubes.ComputeGradientsOff()
	mcubes.ComputeNormalsOff()
	mcubes.SetValue(0, isoValue)

	smoothingIterations = 5
	passBand = 0.001
	featureAngle = 60.0
	smoother = vtk.vtkWindowedSincPolyDataFilter()
	smoother.SetInputConnection(mcubes.GetOutputPort())
	smoother.SetNumberOfIterations(smoothingIterations)
	smoother.BoundarySmoothingOff()
	smoother.FeatureEdgeSmoothingOff()
	smoother.SetFeatureAngle(featureAngle)
	smoother.SetPassBand(passBand)
	smoother.NonManifoldSmoothingOn()
	smoother.NormalizeCoordinatesOn()
	smoother.Update()

	normals = vtk.vtkPolyDataNormals()
	normals.SetInputConnection(smoother.GetOutputPort())
	normals.SetFeatureAngle(featureAngle)

	stripper = vtk.vtkStripper()
	stripper.SetInputConnection(normals.GetOutputPort())

	mapper = vtk.vtkPolyDataMapper()
	mapper.SetInputConnection(stripper.GetOutputPort())

	actor = vtk.vtkActor()
	actor.SetMapper(mapper)
	actor.GetProperty().SetColor(colorLut.GetTableValue(tissueMap[color])[:3])
	actor.GetProperty().SetSpecular(0.5)
	actor.GetProperty().SetSpecularPower(10)

	return actor


	import pathlib
	import pydicom
	from typing import List


	def load_scan(dicom_dir: str) -> List[pydicom.FileDataset]:
	slices = [pydicom.dcmread(f) for f in pathlib.Path(dicom_dir).glob("*.dcm")]
	print(f"file count: {len(slices)}")

	# skip files with no SliceLocation or InstanceNumber(eg scout views)
	slices = list(filter(lambda s: hasattr(s, "SliceLocation"), slices))

	print(f"file count with SliceLocation: {len(slices)}")
	slices = sorted(slices, key=lambda s: s.SliceLocation)

	return slices


	def get_pixels_hu(slices: List[pydicom.FileDataset]) -> np.ndarray:
	# s.pixel_array: 2D array with HW format
	# img3d: 3D array with NHW format
	img3d = np.stack([s.pixel_array for s in slices])

	# Convert to int16 (from sometimes int16),
	# should be possible as values should always be low enough (<32k)
	img3d = img3d.astype(np.int16)

	# Set outside-of-scan pixels to 1
	# The intercept is usually -1024, so air is approximately 0
	# img3d[img3d == -1000] = 0
	img3d[img3d == -2000] = 0

	# Convert to Hounsfield units (HU)

	# Operate on whole array
	intercept = slices[0].RescaleIntercept
	slope = slices[0].RescaleSlope

	if slope != 1:
	img3d = slope * img3d.astype(np.float64)
	img3d = img3d.astype(np.int16)

	img3d += np.int16(intercept)

	# # Operate on individual array
	# for n in range(len(slices)):
	# intercept = slices[n].RescaleIntercept
	# slope = slices[n].RescaleSlope

	# if slope != 1:
	# img3d[n] = slope * img3d[n].astype(np.float64)
	# img3d[n] = img3d[n].astype(np.int16)

	# img3d[n] += np.int16(intercept)

	return np.array(img3d, dtype=np.int16)


	def get_pixels(scans: List[pydicom.FileDataset]) -> np.ndarray:
	# s.pixel_array: 2D array with HW format
	# img3d: 3D array with NHW format
	img3d = np.stack([s.pixel_array for s in scans])

	return img3d


	def plot_stack_sample(stack, rows=6, cols=6, start_with=10, show_every=3):
	import matplotlib.pyplot as plt

	fig, ax = plt.subplots(rows, cols, figsize=[12, 12])
	for i in range(rows * cols):
	ind = start_with + i * show_every
	ax[int(i / rows), int(i % rows)].set_title("slice %d" % ind)
	ax[int(i / rows), int(i % rows)].imshow(stack[ind], cmap="gray")
	ax[int(i / rows), int(i % rows)].axis("off")
	plt.show()


	def resample(image, scan, new_spacing=[1, 1, 1]):
	import scipy.ndimage

	# Determine current pixel spacing
	spacing = np.array(
	[
	float(scan[0].SliceThickness),
	scan[0].PixelSpacing[0],
	scan[0].PixelSpacing[1],
	]
	)

	resize_factor = spacing / new_spacing
	new_real_shape = image.shape * resize_factor
	new_shape = np.round(new_real_shape)
	real_resize_factor = new_shape / image.shape
	new_spacing = spacing / real_resize_factor

	image = scipy.ndimage.interpolation.zoom(image, real_resize_factor)

	return image, new_spacing