Example
    -------
    >>> import ants
    >>> mni = ants.image_read(ants.get_data('mni'))
    >>> ch2 = ants.image_read(ants.get_data('ch2'))
    >>> ants.plot_ortho_stack([mni,mni,mni])

def plot_ortho_stack(images, overlays=None, reorient=True,
    # xyz arguments
    xyz=None, xyz_lines=False, xyz_color='red', xyz_alpha=0.6, xyz_linewidth=2, xyz_pad=5,
    # base image arguments
    cmap='Greys_r', alpha=1,
    # overlay arguments
    overlay_cmap='jet', overlay_alpha=0.9, 
    # background arguments
    black_bg=True, bg_thresh_quant=0.01, bg_val_quant=0.99, 
    # scale/crop/domain arguments
    crop=False, scale=False, domain_image_map=None, 
    # title arguments
    title=None, titlefontsize=24, title_dx=0, title_dy=0, 
    # 4th panel text arguemnts
    text=None, textfontsize=24, textfontcolor='white', text_dx=0, text_dy=0,
    # save & size arguments
    filename=None, dpi=500, figsize=1., colpad=0, rowpad=0,
    transpose=False, transparent=True):
    """
    Example
    -------
    >>> import ants
    >>> mni = ants.image_read(ants.get_data('mni'))
    >>> ch2 = ants.image_read(ants.get_data('ch2'))
    >>> ants.plot_ortho_stack([mni,mni,mni])
    """
    def mirror_matrix(x):
        return x[::-1,:]
    def rotate270_matrix(x):
        return mirror_matrix(x.T)
    def reorient_slice(x, axis):
        return rotate270_matrix(x)

    # need this hack because of a weird NaN warning from matplotlib with overlays
    warnings.simplefilter('ignore')

    n_images = len(images)

    # handle `image` argument
    for i in range(n_images):
        if isinstance(images[i], str):
            images[i] = iio2.image_read(images[i])
        if not isinstance(images[i], iio.ANTsImage):
            raise ValueError('image argument must be an ANTsImage')
        if images[i].dimension != 3:
            raise ValueError('Input image must have 3 dimensions!')

    if overlays is None:
        overlays = [None]*n_images
    # handle `overlay` argument
    for i in range(n_images):
        if overlays[i] is not None:
            if isinstance(overlays[i], str):
                overlays[i] = iio2.image_read(overlays[i])
            if not isinstance(overlays[i], iio.ANTsImage):
                raise ValueError('overlay argument must be an ANTsImage')
            if overlays[i].dimension != 3:
                raise ValueError('Overlay image must have 3 dimensions!')

            if not iio.image_physical_space_consistency(images[i], overlays[i]):
                overlays[i] = reg.resample_image_to_target(overlays[i], images[i], interp_type='linear')

    for i in range(1,n_images):
        if not iio.image_physical_space_consistency(images[0], images[i]):
            images[i] = reg.resample_image_to_target(images[0], images[i], interp_type='linear')

    # reorient images
    if reorient != False:
        if reorient == True:
            reorient = 'RPI'
        
        for i in range(n_images):
            images[i] = images[i].reorient_image2(reorient)
        
        if overlays[i] is not None:
            overlays[i] = overlays[i].reorient_image2(reorient)

    # handle `slices` argument
    if xyz is None:
        xyz = [int(s/2) for s in images[0].shape]
    for i in range(3):
        if xyz[i] is None:
            xyz[i] = int(images[0].shape[i]/2)

    # resample image if spacing is very unbalanced
    spacing = [s for i,s in enumerate(images[0].spacing)]
    if (max(spacing) / min(spacing)) > 3.:
        new_spacing = (1,1,1)
        for i in range(n_images):
            images[i] = images[i].resample_image(tuple(new_spacing))
        if overlays[i] is not None:
            overlays[i] = overlays[i].resample_image(tuple(new_spacing))
        xyz = [int(sl*(sold/snew)) for sl,sold,snew in zip(xyz,spacing,new_spacing)]

    # potentially crop image
    if crop:
        for i in range(n_images):
            plotmask = images[i].get_mask(cleanup=0)
            if plotmask.max() == 0:
                plotmask += 1
            images[i] = images[i].crop_image(plotmask)
            if overlays[i] is not None:
                overlays[i] = overlays[i].crop_image(plotmask)
                
    # pad images
    for i in range(n_images):
        if i == 0:
            images[i], lowpad, uppad = images[i].pad_image(return_padvals=True)
        else:
            images[i] = images[i].pad_image()
        if overlays[i] is not None:
            overlays[i] = overlays[i].pad_image()
    xyz = [v+l for v,l in zip(xyz,lowpad)]

    # handle `domain_image_map` argument
    if domain_image_map is not None:
        if isinstance(domain_image_map, iio.ANTsImage):
            tx = tio2.new_ants_transform(precision='float', transform_type='AffineTransform',
                                         dimension=3)
            for i in range(n_images):
                images[i] = tio.apply_ants_transform_to_image(tx, images[i], domain_image_map)
            
                if overlays[i] is not None:
                    overlays[i] = tio.apply_ants_transform_to_image(tx, overlays[i], 
                                                            domain_image_map, 
                                                            interpolation='linear')       
        elif isinstance(domain_image_map, (list, tuple)):
            # expect an image and transformation
            if len(domain_image_map) != 2:
                raise ValueError('domain_image_map list or tuple must have length == 2')
            
            dimg = domain_image_map[0]
            if not isinstance(dimg, iio.ANTsImage):
                raise ValueError('domain_image_map first entry should be ANTsImage')

            tx = domain_image_map[1]
            for i in range(n_images):
                images[i] = reg.apply_transforms(dimg, images[i], transform_list=tx)
                if overlays[i] is not None:
                    overlays[i] = reg.apply_transforms(dimg, overlays[i], transform_list=tx,
                                                   interpolator='linear')

    # potentially find dynamic range
    if scale == True:
        vmins = []
        vmaxs = []
        for i in range(n_images):
            vmin, vmax = images[i].quantile((0.05,0.95))
            vmins.append(vmin)
            vmaxs.append(vmax)
    elif isinstance(scale, (list,tuple)):
        if len(scale) != 2:
            raise ValueError('scale argument must be boolean or list/tuple with two values')
        vmins = []
        vmaxs = []
        for i in range(n_images):
            vmin, vmax = images[i].quantile(scale)
            vmins.append(vmin)
            vmaxs.append(vmax)
    else:
        vmin = None
        vmax = None

    if not transpose:
        nrow = n_images
        ncol = 3
    else:
        nrow = 3
        ncol = n_images

    fig = plt.figure(figsize=((ncol+1)*2.5*figsize, (nrow+1)*2.5*figsize))
    if title is not None:
        basey = 0.93
        basex = 0.5
        fig.suptitle(title, fontsize=titlefontsize, x=basex+title_dx, y=basey+title_dy)

    if (colpad > 0) and (rowpad > 0):
        bothgridpad = max(colpad, rowpad)
        colpad = 0
        rowpad = 0
    else:
        bothgridpad = 0.0

    gs = gridspec.GridSpec(nrow, ncol, wspace=bothgridpad, hspace=0.0,
                 top=1.-0.5/(nrow+1), bottom=0.5/(nrow+1) + colpad, 
                 left=0.5/(ncol+1) + rowpad, right=1-0.5/(ncol+1))

    # pad image to have isotropic array dimensions
    for i in range(n_images):
        images[i] = images[i].numpy()
        if overlays[i] is not None:
            overlays[i] = overlays[i].numpy()
            overlays[i][np.abs(overlays[i]) == 0] = np.nan

    ####################
    ####################
    for i in range(n_images):
        yz_slice = reorient_slice(images[i][xyz[0],:,:],0)
        if not transpose:
            ax = plt.subplot(gs[i,0])
        else:
            ax = plt.subplot(gs[0,i])
        ax.imshow(yz_slice, cmap=cmap, vmin=vmin, vmax=vmax)
        if overlays[i] is not None:
            yz_overlay = reorient_slice(overlays[i][xyz[0],:,:],0)
            ax.imshow(yz_overlay, alpha=overlay_alpha, cmap=overlay_cmap)
        if xyz_lines:
            # add lines
            l = mlines.Line2D([yz_slice.shape[0]-xyz[1],yz_slice.shape[0]-xyz[1]], 
                              [xyz_pad,yz_slice.shape[0]-xyz_pad],
                              color=xyz_color, alpha=xyz_alpha, linewidth=xyz_linewidth)
            ax.add_line(l)
            l = mlines.Line2D([xyz_pad,yz_slice.shape[1]-xyz_pad], 
                              [yz_slice.shape[1]-xyz[2],yz_slice.shape[1]-xyz[2]],
                color=xyz_color, alpha=xyz_alpha, linewidth=xyz_linewidth)
            ax.add_line(l)
        ax.axis('off')
        ####################
        ####################

        xz_slice = reorient_slice(images[i][:,xyz[1],:],1)
        if not transpose:
            ax = plt.subplot(gs[i,1])
        else:
            ax = plt.subplot(gs[1,i])
        ax.imshow(xz_slice, cmap=cmap, vmin=vmin, vmax=vmax)
        if overlays[i] is not None:
            xz_overlay = reorient_slice(overlays[i][:,xyz[1],:],1)
            ax.imshow(xz_overlay, alpha=overlay_alpha, cmap=overlay_cmap)
        if xyz_lines:
            # add lines
            l = mlines.Line2D([xz_slice.shape[0]-xyz[0],xz_slice.shape[0]-xyz[0]], 
                              [xyz_pad,xz_slice.shape[0]-xyz_pad],
                               color=xyz_color, alpha=xyz_alpha, linewidth=xyz_linewidth)
            ax.add_line(l)
            l = mlines.Line2D([xyz_pad,xz_slice.shape[1]-xyz_pad], 
                              [xz_slice.shape[1]-xyz[2],xz_slice.shape[1]-xyz[2]],
                              color=xyz_color, alpha=xyz_alpha, linewidth=xyz_linewidth)
            ax.add_line(l)
        ax.axis('off')

        ####################
        ####################
        xy_slice = reorient_slice(images[i][:,:,xyz[2]],2)
        if not transpose:
            ax = plt.subplot(gs[i,2])
        else:
            ax = plt.subplot(gs[2,i])
        ax.imshow(xy_slice, cmap=cmap, vmin=vmin, vmax=vmax)
        if overlays[i] is not None:
            xy_overlay = reorient_slice(overlays[i][:,:,xyz[2]],2)
            ax.imshow(xy_overlay, alpha=overlay_alpha, cmap=overlay_cmap)
        if xyz_lines:
            # add lines
            l = mlines.Line2D([xy_slice.shape[0]-xyz[0],xy_slice.shape[0]-xyz[0]], 
                              [xyz_pad,xy_slice.shape[0]-xyz_pad],
                              color=xyz_color, alpha=xyz_alpha, linewidth=xyz_linewidth)
            ax.add_line(l)
            l = mlines.Line2D([xyz_pad,xy_slice.shape[1]-xyz_pad], 
                              [xy_slice.shape[1]-xyz[1],xy_slice.shape[1]-xyz[1]],
                              color=xyz_color, alpha=xyz_alpha, linewidth=xyz_linewidth)
            ax.add_line(l)
        ax.axis('off')

    ####################
    ####################

    if filename is not None:
        plt.savefig(filename, dpi=dpi, transparent=transparent)
        plt.close(fig)
    else:
        plt.show()

    # turn warnings back to default
    warnings.simplefilter('default')