MassProperties

# !/usr/bin/env python3

from pathlib import Path

# noinspection PyUnresolvedReferences
import vtkmodules.vtkInteractionStyle
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkFiltersCore import (
    vtkMassProperties,
    vtkPolyDataNormals,
    vtkTriangleFilter
)
from vtkmodules.vtkFiltersModeling import vtkFillHolesFilter
from vtkmodules.vtkFiltersSources import vtkSphereSource
from vtkmodules.vtkIOGeometry import (
    vtkBYUReader,
    vtkOBJReader,
    vtkSTLReader
)
from vtkmodules.vtkIOLegacy import vtkPolyDataReader
from vtkmodules.vtkIOPLY import vtkPLYReader
from vtkmodules.vtkIOXML import vtkXMLPolyDataReader


def get_program_parameters():
    import argparse
    description = 'Compute volume and surface area of a closed, triangulated mesh.'
    epilogue = '''
    '''
    parser = argparse.ArgumentParser(description=description, epilog=epilogue,
                                     formatter_class=argparse.RawTextHelpFormatter)
    parser.add_argument('file_name', nargs='?', default=None,
                        help='The polydata source file name,e.g. Torso.vtp.')
    args = parser.parse_args()

    return args.file_name


def main():
    file_name = get_program_parameters()
    poly_data = None
    if file_name:
        if Path(file_name).is_file():
            poly_data = read_poly_data(file_name)
        else:
            print(f'{file_name} not found.\nUsing a sphere instead.')
    if file_name is None or poly_data is None:
        source = vtkSphereSource(radius=0.5, phi_resolution=51, theta_resolution=51)
        poly_data = source.update().output

    fill_holes_filter = vtkFillHolesFilter(input_data=poly_data, hole_size=1000.0)
    triangle_filter = vtkTriangleFilter()

    # Make the triangle winding order consistent.
    normals = vtkPolyDataNormals(consistency=True, splitting=False)

    mass_properties = vtkMassProperties()
    fill_holes_filter >> triangle_filter >> normals >> mass_properties
    mass_properties.update()
    fmt = '8.6f'
    print(
        f'Volume: {mass_properties.volume:{fmt}}\n'
        f'    VolumeX: {mass_properties.volume_x:{fmt}}\n '
        f'   VolumeY: {mass_properties.volume_y:{fmt}}\n'
        f'    VolumeZ: {mass_properties.volume_z:{fmt}}\n'
        f'Area: {mass_properties.surface_area:{fmt}}\n'
        f'    MinCellArea: {mass_properties.min_cell_area:{fmt}}\n'
        f'    MaxCellArea: {mass_properties.max_cell_area:{fmt}}\n'
        f'NormalizedShapeIndex: {mass_properties.normalized_shape_index:{fmt}}')


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:
        return reader.update().output
    else:
        return None


if __name__ == '__main__':
    main()