https://github.com/pierre-guillou/pdiags_bench
Tip revision: 8a2de479abf7bcbf2002b4f8f620c635217cc088 authored by Julien J Tierny on 30 January 2023, 16:30:54 UTC
added reference to the arxiv repo
added reference to the arxiv repo
Tip revision: 8a2de47
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)