MedicalDemo3

#!/usr/bin/env python

import sys
from dataclasses import dataclass
from pathlib import Path

# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.vtkCommonCore import (
    vtkLookupTable
)
from vtkmodules.vtkCommonDataModel import (
    vtkPlane
)
from vtkmodules.vtkFiltersCore import (
    vtkFlyingEdges3D,
    vtkMarchingCubes,
    vtkStripper
)
from vtkmodules.vtkFiltersModeling import vtkOutlineFilter
from vtkmodules.vtkIOImage import vtkMetaImageReader
from vtkmodules.vtkImagingCore import vtkImageMapToColors
from vtkmodules.vtkInteractionWidgets import vtkCameraOrientationWidget
from vtkmodules.vtkRenderingCore import (
    vtkActor,
    vtkCamera,
    vtkImageSlice,
    vtkPolyDataMapper,
    vtkProperty,
    vtkRenderWindow,
    vtkRenderWindowInteractor,
    vtkRenderer
)
from vtkmodules.vtkRenderingImage import vtkImageResliceMapper


def main():
    colors = vtkNamedColors()

    file_name, use_flying_edges = get_program_parameters()

    colors.SetColor('SkinColor', 240, 184, 160, 255)
    colors.SetColor('BkgColor', 51, 77, 102, 255)

    # Create the renderer, the render window, and the interactor. The
    # renderer draws into the render window, the interactor enables
    # mouse- and keyboard-based interaction with the data within the
    # render window.
    #
    ren = vtkRenderer()
    ren_win = vtkRenderWindow()
    ren_win.AddRenderer(ren)
    # Set a background color for the renderer and set the name and
    # size of the render window (expressed in pixels).
    ren = vtkRenderer(background=colors.GetColor3d('BkgColor'))
    ren_win = vtkRenderWindow(size=(640, 480),
                              window_name=f'{Path(sys.argv[0]).stem:s}')
    ren_win.AddRenderer(ren)

    iren = vtkRenderWindowInteractor()
    iren.SetRenderWindow(ren_win)
    # Since we import vtkmodules.vtkInteractionStyle we can do this
    # because vtkInteractorStyleSwitch is automatically imported:
    iren.GetInteractorStyle().SetCurrentStyleToTrackballCamera()

    # The following reader is used to read a series of 2D slices (images)
    # that compose the volume. The slice dimensions are set, and the
    # pixel spacing. The data Endianness must also be specified. The
    # reader uses the FilePrefix in combination with the slice number to
    # construct filenames using the format FilePrefix.%d. (In this case
    # the FilePrefix is the root name of the file: quarter.)
    reader = vtkMetaImageReader(file_name=file_name)

    # An isosurface, or contour value of 500 is known to correspond to
    # the skin of the patient.
    # The triangle stripper is used to create triangle
    # strips from the isosurface these render much faster on may
    # systems.
    if use_flying_edges:
        try:
            skin_extractor = vtkFlyingEdges3D()
        except AttributeError:
            skin_extractor = vtkMarchingCubes()
    else:
        skin_extractor = vtkMarchingCubes()
    skin_extractor.SetValue(0, 500)

    skin_stripper = vtkStripper()

    skin_mapper = vtkPolyDataMapper(scalar_visibility=False)

    reader >> skin_extractor >> skin_stripper >> skin_mapper

    bounds = skin_mapper.GetBounds()
    centroid = list(map(lambda x, y: x + (y - x) / 2, bounds[0::2], bounds[1::2]))

    # Set skin to semi-transparent.
    skin_prop = vtkProperty(diffuse_color=colors.GetColor3d('SkinColor'), specular=0.3, specular_power=20, opacity=0.5)
    skin = vtkActor(mapper=skin_mapper, property=skin_prop)

    # An isosurface, or contour value of 1150 is known to correspond to the
    # bone of the patient.
    # The triangle stripper is used to create triangle strips from the
    # isosurface these render much faster on may systems.
    if use_flying_edges:
        bone_extractor = vtkFlyingEdges3D()
    else:
        bone_extractor = vtkMarchingCubes()
    bone_extractor.SetValue(0, 1150)

    bone_stripper = vtkStripper()
    bone_mapper = vtkPolyDataMapper(scalar_visibility=False)
    reader >> bone_extractor >> bone_stripper >> bone_mapper

    bone_prop = vtkProperty(diffuse_color=colors.GetColor3d('Ivory'))
    # Turn off bone for this example.
    bone = vtkActor(mapper=bone_mapper, property=bone_prop, visibility=False)

    # An outline provides context around the data.
    #
    outline_data = vtkOutlineFilter()

    outline_mapper = vtkPolyDataMapper()
    reader >> outline_data >> outline_mapper

    outline_prop = vtkProperty(diffuse_color=colors.GetColor3d('Black'))
    outline = vtkActor(mapper=outline_mapper, property=outline_prop)

    # Now we are creating three orthogonal planes passing through the
    # volume. Each plane uses a different texture map and therefore has
    # different coloration.

    # Start by creating a black/white lookup table.
    bw_lut = vtkLookupTable(table_range=(0, 2000), hue_range=(0, 0), saturation_range=(0, 0), value_range=(0, 1))
    bw_lut.Build()

    # Now create a lookup table that consists of the full hue circle
    # (from HSV).
    hue_lut = vtkLookupTable(table_range=(0, 2000), hue_range=(0, 1), saturation_range=(1, 1), value_range=(1, 1))
    hue_lut.Build()

    # Finally, create a lookup table with a single hue but having a range
    # in the saturation of the hue.
    sat_lut = vtkLookupTable(table_range=(0, 2000), hue_range=(0.6, 0.6), saturation_range=(0, 1), value_range=(1, 1))
    sat_lut.Build()

    # Use vtkImageMapToColors to map the scalar components of an input image
    # through a lookup table to produce an RGBA or RGB output imag
    # Then create a slice planes through a 3D image volume using
    # vtkImageSlice, combining it with a vtkImageResliceMapper and a vtkPlane.

    # Create the first (sagittal) plane of the three planes.
    sagittal_colors = vtkImageMapToColors(lookup_table=bw_lut)
    sp = vtkPlane(origin=centroid, normal=(1, 0, 0))
    sagittal_mapper = vtkImageResliceMapper(slice_plane=sp)
    reader >> sagittal_colors >> sagittal_mapper
    sagittal_slice = vtkImageSlice(mapper=sagittal_mapper)

    # Create the second (axial) plane of the three planes. We use the
    # same approach as before except that the extent differs.
    axial_colors = vtkImageMapToColors(lookup_table=hue_lut)
    ap = vtkPlane(origin=centroid, normal=(0, 0, 1))
    axial_mapper = vtkImageResliceMapper(slice_plane=ap)
    reader >> axial_colors >> axial_mapper
    axial_slice = vtkImageSlice(mapper=axial_mapper)

    # Create the third (coronal) plane of the three planes. We use
    # the same approach as before except that the extent differs.
    coronal_colors = vtkImageMapToColors(lookup_table=sat_lut)
    cp = vtkPlane(origin=centroid, normal=(0, 1, 0))
    coronal_mapper = vtkImageResliceMapper(slice_plane=cp)
    reader >> coronal_colors >> coronal_mapper
    coronal_slice = vtkImageSlice(mapper=coronal_mapper)

    # It is convenient to create an initial view of the data. The
    # FocalPoint and Position form a vector direction. Later on
    # (ResetCamera() method) this vector is used to position the camera
    # to look at the data in this direction.
    camera = vtkCamera(view_up=(0, 0, 1), position=(0, -1, 0), focal_point=(0, 0, 0))
    camera.ComputeViewPlaneNormal()
    camera.Azimuth(30.0)
    camera.Elevation(30.0)

    # Actors are added to the renderer.
    ren.AddActor(outline)
    ren.AddActor(sagittal_slice)
    ren.AddActor(axial_slice)
    ren.AddActor(coronal_slice)
    ren.AddActor(skin)
    ren.AddActor(bone)

    # An initial camera view is created.  The Dolly() method moves
    # the camera towards the FocalPoint, thereby enlarging the image.
    ren.SetActiveCamera(camera)
    ren.ResetCamera()
    camera.Dolly(1.5)

    # Calling Render() directly on a vtkRenderer is strictly forbidden.
    # Only calling Render() on the vtkRenderWindow is a valid call.
    # ren_win.Render()

    # Note that when camera movement occurs (as it does in the Dolly()
    # method), the clipping planes often need adjusting. Clipping planes
    # consist of two planes: near and far along the view direction. The
    # near plane clips out objects in front of the plane; the far plane
    # clips out objects behind the plane. This way only what is drawn
    # between the planes is actually rendered.
    ren.ResetCameraClippingRange()

    cow = make_camera_orientation_widget(ren)
    cow.On()

    # The labels and colors for the various axes.
    category = 'lrpasi'
    cow1 = make_camera_orientation_widget(ren, category, 1)
    cow1.On()

    # Interact with the data.
    iren.Initialize()
    iren.Start()


def get_program_parameters():
    import argparse
    description = 'A Composite image of three planes and translucent skin extracted from a CT dataset of the head.'
    epilogue = '''
    Derived from VTK/Examples/Cxx/Medical3.cxx
    This example reads a volume dataset, extracts two isosurfaces that
     represent the skin and bone, creates three orthogonal planes
     (sagittal, axial, coronal), and displays them.
    '''
    parser = argparse.ArgumentParser(description=description, epilog=epilogue,
                                     formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument('filename', help='FullHead.mhd.')
    parser.add_argument('-m', '--marching_cubes', action='store_false',
                        help='Use Marching Cubes instead of Flying Edges.')
    args = parser.parse_args()
    return args.filename, args.marching_cubes


def set_axes_labels():
    """
    Define the axes labels.

    :return: The axes labels.
    """
    return {
        # Labels are: Anterior, Posterior, Dorsal, Ventral, Left, Right
        'apdvlr': {'+X': 'A', '-X': 'P', '+Y': 'D', '-Y': 'V', '+Z': 'L', '-Z': 'R'},
        'apdvrl': {'+X': 'A', '-X': 'P', '+Y': 'D', '-Y': 'V', '+Z': 'R', '-Z': 'L'},
        'padvlr': {'+X': 'P', '-X': 'A', '+Y': 'D', '-Y': 'V', '+Z': 'L', '-Z': 'R'},
        # Labels are: Left, Right, Superior, Inferior, Anterior, Posterior
        'lrsiap': {'+X': 'L', '-X': 'R', '+Y': 'S', '-Y': 'I', '+Z': 'A', '-Z': 'P'},
        'lrpasi': {'+X': 'L', '-X': 'R', '+Y': 'P', '-Y': 'A', '+Z': 'S', '-Z': 'I'},
        'rlpais': {'+X': 'R', '-X': 'L', '+Y': 'P', '-Y': 'A', '+Z': 'I', '-Z': 'S'},
        # Default labels
        'xyz': None
    }


def set_axes_colors():
    """
    Define the axes colors.

    :return: The axes colors.
    """
    color1 = {'+X': 'IndianRed', '-X': 'FireBrick',
              '+Y': 'LimeGreen', '-Y': 'DarkGreen',
              '+Z': 'Blue', '-Z': 'SteelBlue'}
    return {
        'apdvlr': color1,
        'apdvrl': color1,
        'padvlr': color1,
        'lrsiap': color1,
        'lrpasi': None,
        'rlpais': None,
        # Default colors.
        'xyz': None
    }


def get_axes_params():
    """
    Gather the defined axes labels and colors into a dictionary.

    :return: The dictionary of axes labels and colors.
    """

    # The keys must be the same.
    axes_labels = set_axes_labels()
    axes_colors = set_axes_colors()
    label_keys = set(axes_labels.keys())
    color_keys = set(axes_colors.keys())
    common_keys = label_keys.intersection(color_keys)
    alc = dict()
    for k in common_keys:
        alc[k] = (axes_labels[k], axes_colors[k])
    return alc


def make_camera_orientation_widget(ren, alc_key='xyz', position=3):  # , reposition):
    """
    Make a camera orientation widget for a given renderer.

    position has these values 0: LowerLeft, 1: UpperLeft, 2: LowerRight, 3: UpperRight

    :param ren: The renderer.
    :param alc_key: The key specifying the desired labels and colors for the axes.
    :param position: Position the camera orientation widget.
    :return: The camera orientation widget.
    """

    cow = vtkCameraOrientationWidget(parent_renderer=ren, enabled=True)

    rep = cow.representation

    axes_parameters = get_axes_params()

    if alc_key not in axes_parameters:
        print(
            f'Invalid key for axes labels and colors.'
            f'\nValid keys are: {sorted(axes_parameters.keys())}'
            f'\nUsing the key: xyz.')
        alc_key = 'xyz'

    alc = axes_parameters[alc_key]

    match position:
        case 0:
            rep.AnchorToLowerLeft()
        case 1:
            rep.AnchorToUpperLeft()
        case 2:
            rep.AnchorToLowerRight()
        case _:
            rep.AnchorToUpperRight()

    if not alc[0] is None:
        rep.x_plus_label_text = alc[0]['+X']
        rep.x_minus_label_text = alc[0]['-X']
        rep.y_plus_label_text = alc[0]['+Y']
        rep.y_minus_label_text = alc[0]['-Y']
        rep.z_plus_label_text = alc[0]['+Z']
        rep.z_minus_label_text = alc[0]['-Z']

    if not alc[1] is None:
        colors = vtkNamedColors()
        rep.x_axis_color = colors.GetColor3d(alc[1]['+X'])
        rep.y_axis_color = colors.GetColor3d(alc[1]['+Y'])
        rep.z_axis_color = colors.GetColor3d(alc[1]['+Z'])

    cow.SetRepresentation(rep)

    cow.Off()

    return cow


if __name__ == '__main__':
    main()