#!/usr/bin/python3

"""
Display volumetric data from the Stanford Volume Data Archive. The
data (layers) are expected to be concatenated together.
"""

from mayavi import mlab
import os
import sys

def usage():
    binary = sys.argv[0]
    print(f"Usage: {binary} <data file>")

if len(sys.argv) < 2:
    usage()
    raise SystemExit


def cube_root(n):
    """
    For integers, this probably even works.
    """
    for x in range(0,n+1):
        if x**3 == n:
            return x

    return None


# Read the data in a numpy 3D array
import numpy as np
data = np.fromfile(sys.argv[1], dtype='>u2')

original_data_length = 109*256*256
multiplier = cube_root(len(data) // original_data_length)
data.shape = (109*multiplier, 256*multiplier, 256*multiplier)
data = data.T

# Display the data
mlab.figure( bgcolor=(0,0,0), size=(1000,1000) )

# Our data is scalar (grayscale), but not equally spaced in all
# directions.
src = mlab.pipeline.scalar_field(data)
src.spacing = [1, 1, 1.5]
src.update_image_data = True

mlab.pipeline.volume(src, vmin=1400, vmax=2500)
mlab.show()