Guide for running Davinci commands from Python and passing image data to / from Davinci to Python

run_davinci.py

# Built-in dependencies
import os
import subprocess

# External dependencies (requires install, e.g. with conda)
import h5py
import numpy as np

# Uncomment and edit this if your Davinci interpreter is not at '/usr/local/bin/davinci'
# DAVINCI_PATH = "/path/to/bin/davinci"

def run_davinci(command_list, verbosity=False, logfile=None):
    """
    Run Davinci commands in a cludgey way. Writes a list of Davinci 
    commands to a tmp file and runs it as a .dav script using
    Python subprocess. After script terminates, print stdout line by line, 
    skipping Davinci setup logs.
    Note: subprocess gets no information until the process terminates,
    so you won't get any feedback about the state of your Davinci call
    until it completes all commands (or crashes). 
    
    Parameters
    ----------
    command_list: list of str 
        The Davinci commands to run. Each list element is one "command"
        or line of Davinci code to run
    verbosity: bool
        Print explicit steps and tmp script (useful for debugging).
    logfile: str
        Save the Davinci stdout logs to the file specified

    Return
    ------
    log: str
        The Davinci output
    """
    # DAVINCI_PATH is a global variable I set to specify the exact Davinci
    #  interpreter. On most machines it is: /usr/local/bin/davinci
    if 'DAVINCI_PATH' not in globals():
        DAVINCI_PATH = '/usr/local/bin/davinci'
    else:
        DAVINCI_PATH = globals()['DAVINCI_PATH']
    
    # Build the temp .dav script
    # See http://davinci.asu.edu/index.php?title=User_Guide#Writing_Scripts
    lines = ["#!" + DAVINCI_PATH + " -f"]  # Header for the .dav script
    lines.append("'===START OUTPUT==='")  # Lets us skip davinci setup logs

    lines.extend(command_list)  # Davinci commands, passed as list
    lines.append("\n")  # Always end davinci script with newline
    fstring = "\n".join(lines) # put each element on its own line

    # Write commands to temp davinci script
    fhash = hash(tuple(lines))  # Unique hash to avoid name conflicts
    fname = "tmp_{}.dav".format(fhash)
    with open(fname, "w+") as f:
        if verbosity:
            print("Writing following commands to " + fname + ":")
            print(fstring)
        f.write(fstring)

    # Run .dav script with Python subprocess
    if verbosity:
        print("Running " + fname + "...")
    result = subprocess.check_output(['./' + fname], shell=True)

    # Parse log output from result, skipping the davinci setup junk
    log = ""
    start_output = False
    for line in result.splitlines():
        if start_output:
            log += str(line).strip("bt'\\") + "\n"
        elif line == b'"===START OUTPUT==="':
            start_output = True

    # Write log output to file if one was specified
    if logfile:
        if verbosity:
            print("Writing logs to " + logfile)
        with open(logfile, "w+") as f:
            f.write(log)

    # Remove tmp file
    if verbosity:
        print("Removing file " + fname)
    os.remove("./"+fname)

    return log
  

def read_hdf_data(data, squeeze=False):
    """
    Return data if it is a dataset. If data is a group, recurse.
    """
    d = {}
    for key in list(data.keys()):
        if type(data[key]) == h5py._hl.group.Group:
            # Recurse if data[key] is still a group (nested group)
            d[key] = read_hdf_data(data[key])
        elif type(data[key][()]) == np.ndarray:
            # Swap xy because Davinci writes images as (rows,cols,z)
            # Squeeze to remove empty axes (MxNx1) -> (MxN)
            if squeeze:
                d[key] = np.squeeze(swap_xy(data[key][()]))
            else:
                d[key] = swap_xy(data[key][()])
        else:
            # Save data to dict
            d[key] = data[key][()]
    return d


def read_hdf(fname, squeeze=False):
    """
    Return data in hdf file fname as dict. Recursively fetches nested data.
    """
    with h5py.File(fname, 'r') as f:
        d = read_hdf_data(f, squeeze=squeeze)
    return d


def write_hdf_data(f, data):
    """
    Write data to f if it is a dataset. If data is a dict, recurse.
    """
    for key in list(data.keys()):
        if type(data[key]) == dict:
            ds = f.create_group(key)
            write_hdf_data(ds, data[key])
        elif type(data[key]) == np.ndarray:
            # reverse axes because Davinci reads (z,y,x) as (y,x,z)
            # for some dumb reason
            ds = f.create_dataset(key, data=reverse_axes(data[key]))
        elif type(data[key]) == str:
            dt = "S" + str(len(data[key]))
            ds = f.create_dataset(data[key], (100,), dtype=dt)
        else:
            ds = f.create_dataset(key, data=data[key])
    
    return ds


def write_hdf(fname, data):
    """
    Write dictionary data to hdf file fname as dict.
    """
    with h5py.File(fname, 'w') as f:
        d = write_hdf_data(f, data)
    return d

def swap_xy(arr):
    """
    Returns the source array reshaped from (x, y, z) => (y, x, z).
    Used to swap numpy axes order to Davinci and vice versa since
    numpy expects (x, y, bands) and Davinci expects (rows, cols, bands)
    
    Parameters
    ----------
    arr : numpy 3D array
        image to reshape
    """
    im = np.ma.transpose(arr, [1, 0, 2])
    return im

def reverse_axes(arr):
    """Reverse axis order of nd array arr"""
    axes = np.arange(len(arr.shape))
    return np.ma.transpose(arr, axes[::-1])

run_davinci_tutorial.ipynb