CellsInsideObject

#!/usr/bin/env python3

from pathlib import Path

# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.vtkCommonDataModel import vtkDataObject
from vtkmodules.vtkCommonTransforms import vtkTransform
from vtkmodules.vtkFiltersGeneral import (
    vtkMultiThreshold,
    vtkTransformPolyDataFilter
)
from vtkmodules.vtkFiltersModeling import vtkSelectEnclosedPoints
from vtkmodules.vtkIOGeometry import (
    vtkBYUReader,
    vtkOBJReader,
    vtkSTLReader
)
from vtkmodules.vtkIOLegacy import vtkPolyDataReader
from vtkmodules.vtkIOPLY import vtkPLYReader
from vtkmodules.vtkIOXML import vtkXMLPolyDataReader
from vtkmodules.vtkRenderingCore import (
    vtkActor,
    vtkDataSetMapper,
    vtkRenderWindow,
    vtkRenderWindowInteractor,
    vtkRenderer
)


def get_program_parameters():
    import argparse
    description = 'Read a polydata file of a surface and determine if it is a closed surface.'
    epilogue = '''
This example illustrates how to extract the cells that exist inside a closed surface.
It uses vtkSelectEnclosedPoints to mark points that are inside and outside the surface.
vtkMultiThreshold is used to extract the three meshes into three vtkMultiBlockDataSet's.
The cells completely outside are shown in crimson, completely inside are yellow and
 border cells are green.
A translucent copy of the closed surface helps illustrate the selection process.

If two polydata datasets are provided, the example uses the second as the closed surface.
If only one dataset is provided, the closed surface is generated by rotating the
 first dataset by 90 degrees around its Y axis.

   '''
    parser = argparse.ArgumentParser(description=description, epilog=epilogue,
                                     formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument('filename1', help='Enter a polydata file e.g cow.g.')
    parser.add_argument('filename2', default=None, nargs='?', help='Enter another polydata file e.g cow.g.')
    args = parser.parse_args()
    return args.filename1, args.filename2


def main():
    fn1, fn2 = get_program_parameters()
    poly_data1 = read_poly_data(fn1)
    if fn2:
        poly_data2 = read_poly_data(fn2)
    else:
        # If only one polydata is present, generate a second polydata by
        # rotating the original about its center.
        print('Generating modified poly_data1')
        center = poly_data1.GetCenter()
        transform = vtkTransform()
        transform.Translate(center[0], center[1], center[2])
        transform.RotateY(90.0)
        transform.Translate(-center[0], -center[1], -center[2])
        transform_pd = vtkTransformPolyDataFilter(input_data=poly_data1, transform=transform)
        poly_data2 = transform_pd.update().output

    # Mark points inside with 1 and outside with a 0.
    select = vtkSelectEnclosedPoints(input_data=poly_data1, surface_data=poly_data2)

    # Extract three meshes, one completely inside, one completely
    # outside and on the border between the inside and outside.

    threshold = vtkMultiThreshold()
    # Outside points have a 0 value in ALL points of a cell
    outside_id = threshold.AddBandpassIntervalSet(
        0, 0,
        vtkDataObject.FIELD_ASSOCIATION_POINTS, 'SelectedPoints',
        0, 1)
    # Inside points have a 1 value in ALL points of a cell.
    inside_id = threshold.AddBandpassIntervalSet(
        1, 1,
        vtkDataObject.FIELD_ASSOCIATION_POINTS, 'SelectedPoints',
        0, 1)
    # Border points have a 0 or a 1 in at least one point of a cell.
    border_id = threshold.AddIntervalSet(
        0, 1,
        vtkMultiThreshold.OPEN, vtkMultiThreshold.OPEN,
        vtkDataObject.FIELD_ASSOCIATION_POINTS, 'SelectedPoints',
        0, 0)

    # Select the intervals to be output
    threshold.OutputSet(outside_id)
    threshold.OutputSet(inside_id)
    threshold.OutputSet(border_id)
    (select >> threshold).update()

    # Visualize
    colors = vtkNamedColors()
    outside_color = colors.GetColor3d('Crimson')
    inside_color = colors.GetColor3d('Banana')
    border_color = colors.GetColor3d('Mint')
    surface_color = colors.GetColor3d('Peacock')
    background_color = colors.GetColor3d('Silver')

    # Outside
    outside_mapper = vtkDataSetMapper(input_data=threshold.output.GetBlock(outside_id).GetBlock(0),
                                      scalar_visibility=False)

    outside_actor = vtkActor(mapper=outside_mapper)
    outside_actor.property.diffuse_color = outside_color
    outside_actor.property.specular = 0.6
    outside_actor.property.specular_power = 30

    # Inside
    inside_mapper = vtkDataSetMapper(input_data=threshold.output.GetBlock(inside_id).GetBlock(0),
                                     scalar_visibility=False)

    inside_actor = vtkActor(mapper=inside_mapper)
    inside_actor.property.diffuse_color = inside_color
    inside_actor.property.specular = 0.6
    inside_actor.property.specular_power = 30
    inside_actor.property.edge_visibility = True

    # Border
    border_mapper = vtkDataSetMapper(input_data=threshold.output.GetBlock(border_id).GetBlock(0),
                                     scalar_visibility=False)

    border_actor = vtkActor(mapper=border_mapper)
    border_actor.property.diffuse_color = border_color
    border_actor.property.specular = 0.6
    border_actor.property.specular_power = 30
    border_actor.property.edge_visibility = True

    surface_mapper = vtkDataSetMapper(input_data=poly_data2, scalar_visibility=False)

    # Surface of the object containing the cells.
    surface_actor = vtkActor(mapper=surface_mapper)
    surface_actor.property.diffuse_color = surface_color
    surface_actor.property.opacity = 0.1

    ren = vtkRenderer(background=background_color, use_hidden_line_removal=True)
    ren_win = vtkRenderWindow(size=(640, 480), window_name='CellsInsideObject')
    ren_win.AddRenderer(ren)

    iren = vtkRenderWindowInteractor()
    iren.render_window = ren_win

    ren.AddActor(surface_actor)
    ren.AddActor(outside_actor)
    ren.AddActor(inside_actor)
    ren.AddActor(border_actor)

    ren_win.Render()
    ren.active_camera.Azimuth(30)
    ren.active_camera.Elevation(30)
    ren.active_camera.Dolly(1.25)
    ren_win.Render()

    iren.Start()


def read_poly_data(file_name):
    if not file_name:
        print(f'No file name.')
        return None

    valid_suffixes = ['.g', '.obj', '.stl', '.ply', '.vtk', '.vtp']
    path = Path(file_name)
    ext = None
    if path.suffix:
        ext = path.suffix.lower()
    if path.suffix not in valid_suffixes:
        print(f'No reader for this file suffix: {ext}')
        return None

    reader = None
    if ext == '.ply':
        reader = vtkPLYReader(file_name=file_name)
    elif ext == '.vtp':
        reader = vtkXMLPolyDataReader(file_name=file_name)
    elif ext == '.obj':
        reader = vtkOBJReader(file_name=file_name)
    elif ext == '.stl':
        reader = vtkSTLReader(file_name=file_name)
    elif ext == '.vtk':
        reader = vtkPolyDataReader(file_name=file_name)
    elif ext == '.g':
        reader = vtkBYUReader(file_name=file_name)

    if reader:
        reader.update()
        poly_data = reader.output
        return poly_data
    else:
        return None


if __name__ == '__main__':
    main()