atlas_report.py

#!/usr/bin/env python3
"""
Create a report of intra and inter-observer atlas label statistics
- requires that atlas.py has been run previously on the labels directory
- generates HTML report pages in subdirectory of atlas directory

Usage
----
atlas_report.py -a <atlas directory created by atlas.py>
atlas_report.py -h

Authors
----
Mike Tyszka, Caltech Brain Imaging Center

Dates
----
2017-02-21 JMT Split from atlas.py

License
----
This file is part of atlaskit.

    atlaskit is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    atlaskit is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with atlaskit.  If not, see <http://www.gnu.org/licenses/>.

Copyright
----
2017 California Institute of Technology.
"""

import os
import sys
import argparse
import jinja2
import colorsys
import numpy as np
import nibabel as nib
import matplotlib.pyplot as plt
from datetime import datetime
from skimage.util.montage import montage2d
from skimage import color
from atlas import get_label_name
__version__ = '1.1'


def main():

    # Parse command line arguments
    parser = argparse.ArgumentParser(description='Create labeling report for a probabilistic atlas')
    parser.add_argument('-a', '--atlasdir', required=True, help='Directory containing probabilistic atlas')
    parser.add_argument('--strip', dest='strip', action='store_true', help='Strep prefixes from label names')
    
    # Parse command line arguments
    args = parser.parse_args()
    atlas_dir = args.atlasdir
    strip_prefix = args.strip
    
    print('')
    print('-----------------------------')
    print('Atlas label similarity report')
    print('-----------------------------')

    # Check for atlas directory existence
    if not os.path.isdir(atlas_dir):
        print('Atlas directory does not exist (%s) - exiting' % atlas_dir)
        sys.exit(1)

    # Create report directory within atlas directory
    report_dir = os.path.join(atlas_dir, 'report')
    if not os.path.isdir(report_dir):
        os.mkdir(report_dir)

    print('Atlas directory  : %s' % atlas_dir)
    print('Report directory : %s' % report_dir)
    print('')

    print('Loading similarity metrics')
    intra_stats, inter_stats = load_metrics(atlas_dir)

    # Intra-observer reports (one per observer)
    print('')
    print('Generating intra-observer reports')
    obs_reports = intra_observer_reports(atlas_dir, report_dir, intra_stats, strip_prefix)

    # Inter-observer report
    print('')
    print('Generating inter-observer report')
    inter_observer_report(report_dir, inter_stats, strip_prefix)

    # Summary report page
    print('')
    print('Writing report summary page')
    summary_report(atlas_dir, report_dir, intra_stats, inter_stats, obs_reports, strip_prefix)

    # Clean exit
    sys.exit(0)


def summary_report(atlas_dir, report_dir, intra_metrics, inter_metrics, obs_reports, strip_prefix):
    """
    Summary report for the entire atlas
    - maximum probability projections for all labels

    Parameters
    ----------
    atlas_dir: atlas directory path
    report_dir: report directory path
    intra_metrics: intra-observer metrics tuple
    inter_metrics: inter-observer metrics tuple
    obs_reports: list of intra-observer report tuples (obs, fname)

    Returns
    -------

    """

    # Setup Jinja2 template
    html_loader = jinja2.FileSystemLoader(searchpath=sys.path[0])
    html_env = jinja2.Environment(loader=html_loader)
    html_fname = "atlas_summary.jinja"
    html = html_env.get_template(html_fname)

    # Parse metrics tuples
    label_names, label_nos, observers, templates, intra_dice, intra_haus = intra_metrics
    _, _, _, _, inter_dice, inter_haus = inter_metrics

    # Create grand prob label overlays on bg image
    print('  Generating probability montages')
    montage_fname = overlay_montage(atlas_dir, report_dir, 'prob_atlas.nii.gz')
    colorkey_fname = create_colorkey(atlas_dir, report_dir, 'prob_atlas.nii.gz', strip_prefix)

    # Template variables
    template_vars = {
        "obs_reports": obs_reports,
        "montage_fname": montage_fname,
        "colorkey_fname": colorkey_fname,
        "report_time": datetime.now().strftime('%Y-%m-%d %H:%M')}

    # Finally, process the template to produce our final text.
    output_text = html.render(template_vars)

    # Write page to report directory
    with open(os.path.join(os.path.join(atlas_dir), 'report', 'index.html'), "w") as f:
        f.write(output_text)


def intra_observer_reports(atlas_dir, report_dir, intra_metrics, strip_prefix):
    """
    Generate intra-observer report for each observer

    Parameters
    ----------
    atlas_dir: string
        atlas directory path
    report_dir: string
        report directory path
    intra_metrics: tuple
        containing labelNames, labelNos, observers, templates, dice and haussdorff metrics

    Returns
    -------
    obs_reports: list of tuples (obs, fname)
        List of observer numbers and report filenames
    """

    # Setup Jinja2 template
    html_loader = jinja2.FileSystemLoader(searchpath=sys.path[0])
    html_env = jinja2.Environment(loader=html_loader)
    html_fname = "atlas_intra_observer.jinja"
    html = html_env.get_template(html_fname)

    # Parse metrics tuple
    label_names, label_nos, observers, templates, dice, haus = intra_metrics

    # Determine montage size from number of labels
    ncols = 4
    nrows = np.ceil(len(label_names)/ncols).astype(int)

    # Metric limits
    dlims = 0.0, 1.0
    hlims = 0.0, 10.0

    # Init image filename and stats lists
    intra_dice_imgs = []
    intra_haus_imgs = []
    obs_reports = []

    for obs in observers:

        print('')
        print('Observer %02d' % obs)

        # Generate Dice and Hausdorf similarity matrix figures

        dice_fname = "intra_obs_%02d_dice.png" % obs
        similarity_figure(dice[:,obs,:,:],
                          "Observer %02d Dice Coefficient" % obs,
                          dice_fname,
                          report_dir, label_names, dlims, nrows, ncols, 0.0, 12, strip_prefix)
        intra_dice_imgs.append(dice_fname)

        haus_fname = "intra_obs_%02d_haus.png" % obs
        similarity_figure(haus[:,obs,:,:],
                          "Observer %02d Hausdorff Distance (mm)" % obs,
                          haus_fname,
                          report_dir, label_names, hlims, nrows, ncols, 1e6, 12, strip_prefix)
        intra_haus_imgs.append(haus_fname)

        # Compile stats results for each label for this observer

        obs_stats = []

        for ll, label_name in enumerate(label_names):

            this_intra_dice = dice[ll, obs, :, :]
            this_intra_haus = haus[ll, obs, :, :]

            # Similarity matrices are upper triangle symmetric
            # so calculate upper triangle mean, excluding leading diagonal
            # Returns a string (to allow for '-')
            intra_dice_mean = mean_triu_str(this_intra_dice)
            intra_haus_mean = mean_triu_str(this_intra_haus)

            # Find unfinished template labels
            # Search for NaNs on leading diagonals in intra dice data
            unfinished = str(np.where(np.isnan(np.diagonal(this_intra_dice)))[0])

            label_dict = dict([("label_name", label_name),
                               ("label_no", label_nos[ll]),
                               ("intra_dice_mean", intra_dice_mean),
                               ("intra_haus_mean", intra_haus_mean),
                               ("unfinished", unfinished)])

            obs_stats.append(label_dict)

        # Mean label overlay montage
        print('  Generating mean label montage')
        montage_fname = overlay_montage(atlas_dir, report_dir, 'obs-{0:02d}_label_mean.nii.gz'.format(obs))

        # Template variables
        html_vars = {
            "obs": "{0:02d}".format(obs),
            "montage_fname": montage_fname,
            "dice_fname": dice_fname,
            "haus_fname": haus_fname,
            "obs_stats": obs_stats,
            "report_time": datetime.now().strftime('%Y-%m-%d %H:%M')
        }

        # Render page
        html_text = html.render(html_vars)

        # Write report
        obs_html = "observer_%02d_report.html" % obs
        with open(os.path.join(report_dir, obs_html), "w") as f:
            f.write(html_text)
        obs_reports.append(dict(fname=obs_html, obs="{0:02d}".format(obs)))

    return obs_reports


def inter_observer_report(report_dir, inter_metrics, strip_prefix):
    """
    Generate inter-observer report for each template

    Parameters
    ----------
    report_dir: report directory path
    inter_metrics: tuple containing labelNames, labelNos, observers, templates, dice and haussdorff metrics

    Returns
    -------

    """

    # Setup Jinja2 template
    html_loader = jinja2.FileSystemLoader(searchpath=sys.path[0])
    html_env = jinja2.Environment(loader=html_loader)
    html_fname = "atlas_inter_observer.jinja"
    html = html_env.get_template(html_fname)

    # Parse metrics tuple
    label_names, label_nos, observers, templates, dice, haus = inter_metrics

    # Determine subplot matrix dimensions from number of labels
    # nrows x ncols where nrows = ceil(n_labels/ncols)
    ncols = 4
    nrows = np.ceil(len(label_names)/ncols).astype(int)

    # Metric limits
    dlims = 0.0, 1.0
    hlims = 0.0, 10.0

    # Init image filename lists for HTML template
    inter_dice_imgs = []
    inter_haus_imgs = []

    # Loop over all templates, constructing dice and haus matrix images
    for tt in templates:

        # Create similarity figures over all labels and observers
        dice_fname = "inter_tmp_%02d_dice.png" % tt
        similarity_figure(dice[:,tt,:,:],
                          "Template %02d : Dice Coefficient" % tt,
                          dice_fname,
                          report_dir, label_names, dlims, nrows, ncols, 0.0, 12, strip_prefix)
        inter_dice_imgs.append(dice_fname)

        haus_fname = "inter_tmp_%02d_haus.png" % tt
        similarity_figure(haus[:,tt,:,:],
                          "Template %02d Hausdorff Distance (mm)" % tt,
                          haus_fname,
                          report_dir, label_names, hlims, nrows, ncols, 1e6, 12, strip_prefix)
        inter_haus_imgs.append(haus_fname)

    # Composite all images into a single dictionary list
    inter_imgs = []
    for i, dimg in enumerate(inter_dice_imgs):
        himg = inter_haus_imgs[i]
        inter_imgs.append(dict(dimg=dimg, himg=himg))

    # Template variables
    html_vars = {"inter_imgs": inter_imgs,
                 "report_time": datetime.now().strftime('%Y-%m-%d %H:%M')}

    # Render page
    html_text = html.render(html_vars)

    # Write report
    obs_html = os.path.join(report_dir, "inter_report.html")
    with open(obs_html, "w") as f:
        f.write(html_text)

        
def do_strip_prefixes(label_names):
    # if desired, remove prefixes from label names 
    stripped_label_names = []
    for l in label_names:
        stripped_label_names.append(do_strip_prefix(l))
    return stripped_label_names


def do_strip_prefix(label_name):
    # if desired, remove prefixes from label names 
    idx = label_name.rfind('_') + 1
    stripped_label_name = label_name[idx:]
    # print("Stripping atlas label name, from %s to %s" % (label_name, stripped_label_name))
    return stripped_label_name


def overlay_montage(atlas_dir, report_dir, overlay_fname):
    """
    Construct a montage of colored label overlays on a T1w background
    - Each label is colored according to the ITK-SNAP label key
    - Calculate coronal slice skip from minimum BB for 4 x 4 montage (16 slices)

    Parameters
    ----------
    atlas_dir: string
        atlas directory path
    report_dir: string
        report directory path
    overlay_fname: string
        4D overlay image filename (within atlas_dir)

    Returns
    -------
    montage_png: prob label montage
    """

    # Use ITK-SNAP label key colors
    atlas_color = False

    # CIT atlas directory from shell environment
    cit_dir = os.environ['CIT168_DIR']

    # Load label key from atlas directory
    label_key = load_key(os.path.join(atlas_dir, 'labels.txt'))

    # Extract HSV label colors (n_labels x 3 array)
    hsv = label_rgb2hsv(label_key)

    # Probability threshold for minimum BB
    p_thresh = 0.25

    # Size of coronal section montage
    n_rows, n_cols = 6, 6

    # Load background image
    print('  Loading background image')
    bg_fname = os.path.join(cit_dir, 'CIT168_700um', 'CIT168_T1w_700um.nii.gz')
    bg_nii = nib.load(bg_fname)
    bg_img = bg_nii.get_data()

    # Normalize background intensity range to [0,1]
    bg_img = bg_img / np.max(bg_img)

    # Load the 4D probabilistic atlas
    print('  Loading probabilistic image')
    p_nii = nib.load(os.path.join(atlas_dir, overlay_fname))
    p_atlas = p_nii.get_data()

    # Count prob labels
    n_labels = p_atlas.shape[3]

    # Find minimum bounding box for all prob labels > 0.25
    # x0, y0, z0 : minimum corner of BB (closest to origin)
    print('  Determining minimum isotropic bounding box')
    p_all = np.sum(p_atlas, axis=3)
    x0, x1, y0, y1, z0, z1 = bb(p_all > p_thresh, padding=4)

    # Crop bg image and prob atlas
    bg_crop = bg_img[x0:x1, y0:y1, z0:z1]
    p_crop = p_atlas[x0:x1, y0:y1, z0:z1, :]

    # Create montage of coronal sections through cropped bg image
    bg_mont = coronal_montage(bg_crop, n_rows, n_cols)

    bg_mont_rgb = tint(bg_mont, hue=0.0, saturation=0.0)

    # Initialize the all-label overlay
    overlay_mont_rgb = np.zeros_like(bg_mont_rgb)

    # Create equivalent montage for all prob labels with varying hues
    for lc in range(0, n_labels):

        # Construct prob label montage
        p_mont = coronal_montage(p_crop[:,:,:,lc], n_rows, n_cols)

        # Hue and saturation for label overlay
        if atlas_color:
            # Pull HSV from ITK-SNAP label key
            hue, sat, val = hsv[lc, 0], hsv[lc, 1], hsv[lc,2]
        else:
            # Calculate rotating hue
            hue = float(np.mod(lc * 3, n_labels)) / n_labels
            sat, val = 1.0, 1.0

        # Tint the montage
        p_mont_rgb = tint(p_mont, hue=hue, saturation=sat, value=val)

        # Add tinted overlay to running total
        overlay_mont_rgb += p_mont_rgb

    # Composite prob atlas overlay on bg image
    mont_rgb = composite(overlay_mont_rgb, bg_mont_rgb)

    # Create figure and render montage
    fig = plt.figure(figsize=(15,10), dpi=100)
    plt.imshow(mont_rgb, interpolation='none')
    plt.axis('off')
    plt.legend()
    
    # Save figure to PNG
    montage_fname = overlay_fname.replace('.nii.gz', '_montage.png')
    print('  Saving image to %s' % montage_fname)
    plt.savefig(os.path.join(report_dir, montage_fname), bbox_inches='tight')

    return montage_fname


def create_colorkey(atlas_dir, report_dir, overlay_fname, strip_prefix):
    """
    Construct an montage of colored label overlays on a T1w background
    - Each label is colored according to the ITK-SNAP label key
    - Calculate coronal slice skip from minimum BB for 4 x 4 montage (16 slices)

    Parameters
    ----------
    atlas_dir: string
        atlas directory path
    report_dir: string
        report directory path
    overlay_fname: string
        4D overlay image filename (within atlas_dir)

    Returns
    -------
    montage_png: prob label montage
    """

    # Use ITK-SNAP label key colors
    atlas_color = False

    # CIT atlas directory from shell environment
    cit_dir = os.environ['CIT168_DIR']

    # Load label key from atlas directory
    label_key = load_key(os.path.join(atlas_dir, 'labels.txt'))
        
    # Extract HSV label colors (n_labels x 3 array)
    hsv = label_rgb2hsv(label_key)

    # Load the 4D probabilistic atlas
    print('  Loading probabilistic image')
    p_nii = nib.load(os.path.join(atlas_dir, overlay_fname))
    p_atlas = p_nii.get_data()

    # Count prob labels
    n_labels = p_atlas.shape[3]

    rgb_colors = []
    labels = []
    # Create equivalent montage for all prob labels with varying hues
    for lc in range(0, n_labels):

        # Hue and saturation for label overlay
        if atlas_color:
            # Pull HSV from ITK-SNAP label key
            hue, sat, val = hsv[lc, 0], hsv[lc, 1], hsv[lc,2]
        else:
            # Calculate rotating hue
            hue = float(np.mod(lc * 3, n_labels)) / n_labels
            sat, val = 1.0, 1.0
        rgb = colorsys.hsv_to_rgb(hue, sat, val)
        rgb_colors.append(rgb)

    rgb_color_array = np.array(rgb_colors)
    
    x = np.array([0] * n_labels)
    y = np.linspace(1, n_labels, n_labels)

    fig, ax = plt.subplots()
    fig_size = fig.get_size_inches()
    fig.set_size_inches([1.0, 4.5]) # float(fig_size[0])/6, fig_size[1]])
    ax.scatter(x, y, c = rgb_color_array, edgecolors='none', s=25)
    plt.axis('off')

    for i in range(0, n_labels):
        if strip_prefix:
            label_name = do_strip_prefix(label_key['Name'][i])
        ax.annotate(label_name, (x[i],y[i]), xytext=(5,0), textcoords='offset points')
            
    
    # Save figure to PNG
    colorkey_fname = overlay_fname.replace('.nii.gz', '_colorkey.png')
    print('  Saving image to %s' % colorkey_fname)
    plt.savefig(os.path.join(report_dir, colorkey_fname), bbox_inches='tight', transparent = True, pad_inches=0)

    return colorkey_fname


def label_rgb2hsv(label_key):
    """
    Extract label RGB colors and convert to HSV

    Parameters
    ----------
    label_key: data frame

    Returns
    -------
    hsv: numpy array

    """

    rgb = np.array(label_key[['R','G','B']]) / 255.0
    rgb = rgb.reshape([rgb.shape[0], 1, 3])
    hsv = color.rgb2hsv(rgb)
    hsv = hsv.reshape([-1,3])

    return hsv


def coronal_montage(img, n_rows=4, n_cols=4, flip_x=False, flip_y=True, flip_z=True):
    """
    Create a montage of all coronal (XZ) slices from a 3D image

    Parameters
    ----------
    img: 3D image to montage
    n_rows: number of montage rows
    n_cols: number of montage columns
    rot: CCW 90deg rotations to apply to each section

    Returns
    -------

    cor_mont: coronal slice montage of img
    """

    # Total number of sections to extract
    n = n_rows * n_cols

    # Source image dimensions
    nx, ny, nz = img.shape

    # Coronal (XZ) sections
    yy = np.linspace(0, ny-1, n).astype(int)
    cors = img[:,yy,:]

    if flip_x:
        cors = np.flip(cors, axis=0)
    if flip_y:
        cors = np.flip(cors, axis=1)
    if flip_z:
        cors = np.flip(cors, axis=2)

    # Permute image axes for montage2d: original y becomes new x
    img = np.transpose(cors, (1,2,0))

    # Construct montage of coronal sections
    cor_mont = montage2d(img, fill=0, grid_shape=(n_rows, n_cols))

    return cor_mont


def tint(image, hue=0.0, saturation=1.0, value=1.0):
    """
    Add color of the given hue to an RGB image

    Parameters
    ----------
    image
    hue
    saturation

    Returns
    -------
    """

    hsv = np.zeros([image.shape[0], image.shape[1], 3])
    hsv[:, :, 0] = hue
    hsv[:, :, 1] = saturation
    hsv[:, :, 2] = image * value

    return color.hsv2rgb(hsv)


def composite(overlay_rgb, background_rgb):
    """
    Alpha composite RGB overlay on RGB background
    - derive alpha from HSV value of overlay

    Parameters
    ----------
    overlay_rgb:
    background_rgb:

    Returns
    -------

    """

    overlay_hsv = color.rgb2hsv(overlay_rgb)
    value = overlay_hsv[:,:,2]
    alpha_rgb = np.dstack((value, value, value))

    composite_rgb = overlay_rgb * alpha_rgb + background_rgb * (1.0 - alpha_rgb)

    return composite_rgb


def similarity_figure(metric, img_title, img_fname, report_dir, label_names, mlims, nrows, ncols, nansub=0.0, fontsize=8, strip_prefix=False):
    """
    Plot an array of similarity matrix figures for a given observer or template

    Parameters
    ----------
    metric: 3D numpy float array
        similarity metric array to plot
    img_title: string
        image title
    img_fname: string
        output image filename
    report_dir: string
        report directory
    label_names: string list
        list of label names
    mlims: float tuple
        scale limits for metric
    nrows: int
        plot grid rows
    ncols: int
        plot grid columns
    nansub: float
        value to replace NaNs in data

    Returns
    -------
    """

    # Create figure with subplot array
    fig, axs = plt.subplots(nrows, ncols)
    axs = np.array(axs).reshape(-1)

    im = []

    for aa, ax in enumerate(axs):

        if aa < len(label_names):
            mmaa = np.flipud(metric[aa, :, :]).copy()
            mmaa[np.isnan(mmaa)] = nansub
            im = ax.pcolor(mmaa, vmin=mlims[0], vmax=mlims[1], cmap='Spectral')
            if strip_prefix:
                label_name = do_strip_prefix(label_names[aa])
            else:
                label_name = label_names[aa]
            ax.set_title(label_name, fontsize=fontsize)
        else:
            ax.axis('off')

        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)
        ax.set(adjustable='box-forced', aspect='equal')

    # Tidy up spacing
    plt.tight_layout()

    # Make space for title and colorbar
    fig.subplots_adjust(bottom=0.1, top=0.9, left=0.1, right=0.8)
    plt.suptitle(img_title, x=0.5, y=0.99)
    cax = fig.add_axes([0.85, 0.1, 0.05, 0.8])  # [x0, y0, w, h]
    fig.colorbar(im, cax=cax)

    # Save figure to PNG
    print('  Saving image to %s' % img_fname)
    plt.savefig(os.path.join(report_dir, img_fname), bbox_inches='tight')

    # Clean up
    plt.close(fig)


def bb(mask, padding=8):
    """
    Determine minimum bounding box containing all non-zero voxels in mask

    Parameters
    ----------
    mask: 3D boolean array
        binary mask containing all regions
    padding: integer
        voxel padding around minimum BB

    Returns
    -------
    x0, x1, y0, y1, z0, z1: bounding box limits
    """

    # Mask dimensions
    nx, ny, nz = mask.shape

    # MIP in x, y, z
    xproj = np.max(np.max(mask, axis=2), axis=1)
    yproj = np.max(np.max(mask, axis=2), axis=0)
    zproj = np.max(np.max(mask, axis=1), axis=0)

    # Non-zero indices in each projection
    xnz = np.nonzero(xproj)[0]
    ynz = np.nonzero(yproj)[0]
    znz = np.nonzero(zproj)[0]

    # Min and max limits of non-zero projection
    x0, x1 = np.min(xnz), np.max(xnz)
    y0, y1 = np.min(ynz), np.max(ynz)
    z0, z1 = np.min(znz), np.max(znz)

    # Add padding then clip to image bounds
    x0 = np.clip(x0 - padding, 0, nx-1)
    x1 = np.clip(x1 + padding, 0, nx-1)
    y0 = np.clip(y0 - padding, 0, ny-1)
    y1 = np.clip(y1 + padding, 0, ny-1)
    z0 = np.clip(z0 - padding, 0, nz-1)
    z1 = np.clip(z1 + padding, 0, nz-1)

    return x0, x1, y0, y1, z0, z1


def load_metrics(atlas_dir):
    """
    Parse similarity metrics from CSV file

    Parameters
    ----------
    atlas_dir: atlas directory

    Returns
    -------
    m : numpy array containing label, observer and template indices and metrics
    """

    #
    # Load intra-observer metrics
    # Ignore number of voxels in each label (nA, nB) for now
    #

    intra_csv = os.path.join(atlas_dir, 'intra_observer_metrics.csv')
    m = np.genfromtxt(intra_csv,
                      dtype=None,
                      names=['labelName', 'labelNo', 'observer', 'tmpA', 'tmpB', 'dice', 'haus', 'nA', 'nB'],
                      delimiter=',', skip_header=1)

    # Find unique label numbers with initial row indices for each
    label_nos, idx = np.unique(m['labelNo'], return_index=True)

    # Extract corresponding label names to unique label numbers
    label_names = m['labelName'][idx].astype(str)
            
    # Unique template and observer lists can be sorted as usual
    observers = np.unique(m['observer'])
    templates = np.unique(m['tmpA'])

    # Count labels, templates and observers
    n_labels, n_templates, n_observers = len(label_names), len(templates), len(observers)

    # Cast to float and reshape metrics
    dice = m['dice'].reshape(n_labels, n_observers, n_templates, n_templates)
    haus = m['haus'].reshape(n_labels, n_observers, n_templates, n_templates)

    # Composite into intra_metrics tuple
    intra_metrics = label_names, label_nos, observers, templates, dice, haus

    #
    # Load inter-observer metrics
    # Ignore number of voxels in each label (nA, nB) for now
    #

    inter_csv = os.path.join(atlas_dir, 'inter_observer_metrics.csv')
    m = np.genfromtxt(inter_csv,
                      dtype=[('labelName', 'a32'), ('labelNo', 'u8'),
                             ('template', 'u8'), ('obsA', 'u8'), ('obsB', 'u8'),
                             ('dice', 'f8'), ('haus', 'f8'),
                             ('nA', 'u8'), ('nB', 'u8')],
                      delimiter=',', skip_header=1)

    # Find unique label numbers with initial row indices for each
    label_nos, idx = np.unique(m['labelNo'], return_index=True)

    # Extract corresponding label names to unique label numbers
    label_names = m['labelName'][idx].astype(str)

    # Unique template and observer lists can be sorted as usual
    templates = np.unique(m['template'])
    observers = np.unique(m['obsA'])

    # Count labels, templates and observers
    n_labels, n_templates, n_observers = len(label_names), len(templates), len(observers)

    # Cast to float and reshape metrics
    dice = m['dice'].reshape(n_labels, n_templates, n_observers, n_observers)
    haus = m['haus'].reshape(n_labels, n_templates, n_observers, n_observers)

    # Composite into inter_metrics tuple
    inter_metrics = label_names, label_nos, observers, templates, dice, haus

    return intra_metrics, inter_metrics


def mean_triu_str(x):
    """
    Calculate mean of upper triangle (excluding NaNs)
    Returns '-' when upper triangle is entirely NaNs

    Parameters
    ----------
    x: numpy float array

    Returns
    -------
    xms: formatted mean string

    """

    # Size of square matrix, x
    n = x.shape[0]

    # Upper triangle values, excluding leading diagonal
    xut = x[np.triu_indices(n, 1)]

    # Number of NaNs in upper triangle
    n_nans = np.sum(np.isnan(xut))

    if n_nans == xut.size:
        xms = "-"
    else:
        xms = "%0.3f" % np.nanmean(xut)

    return xms


def load_key(key_fname):
    """
    Parse an ITK-SNAP label key file

    Parameters
    ----------
    key_fname: ITK-SNAP label key filename

    Returns
    -------
    key: Data table containing ITK-SNAP style label key
    """

    import pandas as pd

    # Import key as a data table
    # Note the partially undocumented delim_whitespace flag
    key = pd.read_table(key_fname,
                        comment='#',
                        header=None,
                        names=['Index', 'R', 'G', 'B', 'A', 'Vis', 'Mesh', 'Name'],
                        delim_whitespace=True)

    return key


# def maxprob_projections(atlas_dir, report_dir, label_names, nrows, ncols):
#     """
#     *** CURRENTLY UNUSED ***
#
#     Construct an array of maximum probablity projections through each label
#     over all observers and templates
#
#     Parameters
#     ----------
#     atlas_dir: atlas directory path
#     report_dir: report directory path
#     label_names: list of unique label names (in label number order)
#     nrows, ncols: figure matrix size
#
#     Returns
#     -------
#     mpp_png: maxprob projection PNG filename
#     """
#
#     # Probability threshold for minimum BB
#     p_thresh = 0.25
#
#     # Load prob atlas
#     prob_nii = nib.load(os.path.join(atlas_dir, 'prob_atlas.nii.gz'))
#     prob_atlas = prob_nii.get_data()
#
#     # Create figure with subplot array
#     fig, axs = plt.subplots(nrows, ncols, figsize=(8,4))
#     axs = np.array(axs).reshape(-1)
#
#     # Loop over axes
#     for aa, ax in enumerate(axs):
#
#         if aa < len(label_names):
#
#             print('    %s' % label_names[aa])
#
#             # Current prob label
#             p = prob_atlas[:, :, :, aa]
#
#             # Create tryptic of central slices through ROI defined by p > p_thresh
#             tryptic = central_slices(p, isobb(p > p_thresh))
#
#             ax.pcolor(tryptic)
#             ax.set_title(label_names[aa], fontsize=8)
#
#         else:
#             ax.axis('off')
#
#         ax.get_xaxis().set_visible(False)
#         ax.get_yaxis().set_visible(False)
#         ax.set(adjustable='box-forced', aspect='equal')
#         ax.set(aspect='equal')
#
#     # Tidy up spacing
#     plt.tight_layout()
#
#     # Save figure to PNG
#     mpp_fname = 'mpp.png'
#     plt.savefig(os.path.join(report_dir, mpp_fname))
#
#     # Clean up
#     plt.close(fig)
#
#     return mpp_fname


# def central_slices(p, roi):
#     """
#     *** RETIRED ***
#
#     Create tryptic of central slices from ROI
#
#     Parameters
#     ----------
#     img: 3D numpy array
#     roi: tuple containing (x0, y0, z0, w) for ROI
#
#     Returns
#     -------
#     pp: numpy tryptic of central slices
#     """
#
#     # Base output image dimensions
#     base_dims = 64, 3 * 64
#
#     # Unpack ROI
#     x0, y0, z0, w = roi
#
#     # Half width of ROI
#     hw = np.ceil(w / 2.0).astype(int) + 1
#
#     if w > 0:
#
#         # Define central slices
#         xx = slice(x0 - hw, x0 + hw, 1)
#         yy = slice(y0 - hw, y0 + hw, 1)
#         zz = slice(z0 - hw, z0 + hw, 1)
#
#         # Extract slices
#         p_xy = p[xx, yy, z0]
#         p_xz = p[xx, y0, zz]
#         p_yz = p[x0, yy, zz]
#
#         # Create horizontal tryptic
#         pp = np.hstack([p_xy, p_xz, p_yz])
#
#         # Resize to base size
#         pp = imresize(pp, base_dims, interp='bicubic')
#
#     else:
#
#         # Blank tryptic
#         pp = np.zeros(base_dims)
#
#     return pp

# This is the standard boilerplate that calls the main() function.
if __name__ == '__main__':
    try:
        cit_dir = os.environ['CIT168_DIR']
    except KeyError:
        print('* Environmental variable CIT168_DIR not set - exiting')
        sys.exit(1)
        
    main()