convert_datasets.py

import argparse
import enum
import logging
import pathlib
import time

from paraview import simple

import vti2nc3

RESAMPL_3D = 192
RESAMPL_2D = 4096
RESAMPL_1D = 1024 ** 2
logging.basicConfig(format="%(asctime)s %(levelname)s %(message)s", level=logging.INFO)


def write_output(outp, fname, out_dir, explicit):
    if out_dir:
        fname = out_dir + "/" + fname

    partial = pathlib.Path(".not_all_apps").exists()
    if partial and not explicit:
        return

    # Dipha Explicit Complex (Dipha) or Image Data (Dipha, CubicalRipser)
    simple.SaveData(fname + ".dipha", proxy=outp)

    if explicit:
        # vtkUnstructuredGrid (TTK)
        simple.SaveData(fname + ".vtu", proxy=outp)
        # TTK Simplicial Complex (Gudhi, Dionysus, Ripser)
        simple.SaveData(fname + ".tsc", proxy=outp)
        # PHAT ASCII boundary_matrix file format
        simple.SaveData(fname + ".phat", proxy=outp)
        if not partial:
            # Perseus Uniform Triangulation (Perseus)
            simple.SaveData(fname + ".pers", proxy=outp)
            # Eirene.jl Sparse Column Format CSV
            simple.SaveData(fname + ".eirene", proxy=outp)
            # Oineus Custom Simplicial Complex Format
            simple.SaveData(fname + ".oin", proxy=outp)
    else:
        # vtkImageData (TTK)
        simple.SaveData(fname + ".vti", proxy=outp)
        # Perseus Cubical Grid (Perseus, Gudhi)
        simple.SaveData(fname + ".pers", proxy=outp)
        # NetCDF3 (Diamorse)
        vti2nc3.main(fname + ".vti")


def read_file(input_file):
    extension = input_file.split(".")[-1]
    if extension == "vti":
        return simple.XMLImageDataReader(FileName=input_file)

    if extension == "raw":
        extent, dtype = input_file.split(".")[0].split("_")[-2:]
        extent = [int(dim) for dim in extent.split("x")]

        dtype_pv = {
            "uint8": "unsigned char",
            "int16": "short",
            "uint16": "unsigned short",
            "float32": "float",
            "float64": "double",
        }

        raw = simple.ImageReader(FileNames=[input_file])
        raw.DataScalarType = dtype_pv[dtype]
        raw.DataByteOrder = "LittleEndian"
        raw.DataExtent = [0, extent[0] - 1, 0, extent[1] - 1, 0, extent[2] - 1]
        return raw

    return None


class SliceType(enum.Enum):
    VOL = 0
    SURF = 1
    LINE = 2

    @classmethod
    def from_filename(cls, fname):
        if "x1x1_" in fname:
            return cls.LINE
        if "x1_" in fname:
            return cls.SURF
        return cls.VOL


def slice_data(input_dataset, slice_type, dims):
    if slice_type in (SliceType.SURF, SliceType.LINE):
        # force generation of input_vti (otherwise, slice is empty)
        simple.Show(input_dataset)
        # slice along depth/z axis
        sl0 = simple.Slice(Input=input_dataset)
        sl0.SliceType.Normal = [0.0, 0.0, 1.0]

        if slice_type == SliceType.LINE:
            # resample to a 2D strip before slicing again
            rsi = simple.ResampleToImage(Input=sl0)
            rsi.SamplingDimensions = [dims[0], 3, 1]
            simple.Show(rsi)  # same here...

            # slice along vertical/y axis
            sl1 = simple.Slice(Input=rsi)
            sl1.SliceType.Normal = [0.0, 1.0, 0.0]
            return sl1

        # resample to something like 4096x4096x1
        rsi = simple.ResampleToImage(Input=sl0)
        rsi.SamplingDimensions = dims
        return rsi

    # resample to something like 192^3
    rsi = simple.ResampleToImage(Input=input_dataset)
    rsi.SamplingDimensions = dims
    return rsi


def pipeline(raw_file, raw_stem, dims, slice_type, out_dir):
    reader = read_file(raw_file)
    # convert input scalar field to float
    calc = simple.Calculator(Input=reader)
    calc.Function = "ImageFile"
    calc.ResultArrayType = "Float"
    calc.ResultArrayName = "ImageFile"

    # get a slice
    cut = slice_data(calc, slice_type, dims)

    # compute order field
    arrprec = simple.TTKArrayPreconditioning(Input=cut)
    arrprec.PointDataArrays = ["ImageFile"]

    # trash input scalar field, save order field
    pa = simple.PassArrays(Input=arrprec)
    pa.PointDataArrays = ["ImageFile_Order"]

    # save implicit mesh
    if slice_type != SliceType.LINE:
        write_output(pa, raw_stem + "_order_impl", out_dir, False)

    # tetrahedralize grid
    tetrah = simple.Tetrahedralize(Input=pa)
    # remove vtkGhostType arrays (only applies on vtu & vtp)
    rgi = simple.RemoveGhostInformation(Input=tetrah)
    # save explicit mesh
    write_output(rgi, raw_stem + "_order_expl", out_dir, True)


def main(raw_file, out_dir="", resampl_size=RESAMPL_3D, slice_type=SliceType.VOL):
    if raw_file == "":
        return

    if slice_type == SliceType.VOL:
        dims = [resampl_size] * 3
    elif slice_type == SliceType.SURF:
        dims = [resampl_size] * 2 + [1]
    elif slice_type == SliceType.LINE:
        dims = [resampl_size] + [1] * 2
    extent_s = "x".join([str(d) for d in dims])

    raw_stem = raw_file.split(".")[0].split("/")[-1]
    try:
        raw_stem_parts = raw_stem.split("_")
        # update extent
        raw_stem_parts[-2] = extent_s
        # remove data type in file name
        raw_stem_parts.pop()
        raw_stem = "_".join(raw_stem_parts)
    except IndexError:
        # not an Open-Scivis-Datasets raw file (elevation or random)
        raw_stem = f"{raw_stem}_{extent_s}"

    logging.info("Converting %s to input formats (resampled to %s)", raw_file, extent_s)
    beg = time.time()

    pipeline(raw_file, raw_stem, dims, slice_type, out_dir)

    end = time.time()
    logging.info("Converted %s (took %ss)", raw_file, round(end - beg, 3))


if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        description="Generate input files from Open-Scivis Raw format"
    )

    parser.add_argument("raw_file", type=str, help="Path to input Raw file")
    parser.add_argument(
        "-d", "--dest_dir", type=str, help="Destination directory", default="datasets"
    )
    parser.add_argument(
        "-s",
        "--resampling_size",
        type=int,
        help="Resampling to a cube of given vertices edge",
    )
    parser.add_argument(
        "-2",
        "--slice",
        action="store_true",
        help="Generate a 2D slice",
    )
    parser.add_argument(
        "-1",
        "--line",
        action="store_true",
        help="Generate a 1D line",
    )
    args = parser.parse_args()

    if args.line and args.slice:
        raise argparse.ArgumentError

    if args.slice:
        stype = SliceType.SURF
        if args.resampling_size is None:
            args.resampling_size = RESAMPL_2D
    elif args.line:
        stype = SliceType.LINE
        if args.resampling_size is None:
            args.resampling_size = RESAMPL_1D
    else:
        stype = SliceType.VOL
        if args.resampling_size is None:
            args.resampling_size = RESAMPL_3D

    main(args.raw_file, args.dest_dir, args.resampling_size, stype)