Localize small particles and determine how the presence of cells affects their position.

analyze_stack.py

#!/usr/bin/env python

import argparse as ap

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

VOXEL_SIZE = (0.2059, 0.2059, 0.5669)

def process_single_stack(stack):
    fluor = stack[:,0,:,:]
    bf = stack[:,1,:,:]

    # Get cell locations in XY
    from skimage.filters import threshold_otsu
    fluor_mp = fluor.max(axis=0)
    cell_map = fluor_mp > threshold_otsu(fluor_mp)

    # Get particle locations in XY
    min_proj = bf.mean(axis=0)

    from skimage.morphology import black_tophat, disk
    bth = black_tophat(min_proj, selem=disk(6))

    from skimage.feature import blob_log
    blobs = blob_log(-min_proj, min_sigma=2, max_sigma=5, threshold=10, overlap=0.2)

    # Use Z profile to localize particles in Z
    zpros = []

    for blob in blobs:
        y, x, r = blob.astype('int')
        ys, xs = np.where(disk(r))
        #mask = np.zeros_like(bth)
        #mask[(y-r):(y+r+1),(x-r):(x+r+1)] = disk(r)

        ys, xs = zip(*[(a, b) for (a, b) in zip(ys+y, xs+x) if (-1 < min(a, b)) and (max(a, b) < 1024)])
        zpro = bf[:,ys,xs].mean(axis=-1)
        zpros.append(zpro)

    zpros = np.vstack(zpros)

    from scipy.ndimage.filters import gaussian_filter1d
    smoothpros = gaussian_filter1d(zpros, sigma=1, axis=-1)
    depth = smoothpros.argmin(axis=-1)

    on_cell = []
    off_cell = []

    for blob, d in zip(blobs, depth):
        y, x, r = blob.astype(int)
        if cell_map[y,x]:
            on_cell.append(d * VOXEL_SIZE[2])
        else:
            off_cell.append(d * VOXEL_SIZE[2])

    # Setup plots
    sns.set_color_codes()
    sns.set_style('white')
    sns.set_context('poster')

    import matplotlib.gridspec as gridspec
    fig = plt.figure()
    gs = gridspec.GridSpec(4, 4)

    # Plot: BF + particle detection
    ax_bf = plt.subplot(gs[0:2,0:2])
    ax_bf.imshow(min_proj, cmap='gray')

    for blob in blobs:
        y, x, r = blob
        c = plt.Circle((x, y), r, color='r', lw=1, fill=False, alpha=0.7)
        ax_bf.add_patch(c)

    ax_bf.set_axis_off()

    # Plot: Fluor + depth color map
    ax_fl = plt.subplot(gs[0:2, 2:4])
    cm = plt.get_cmap('inferno')
    H = zpros.shape[-1]

    def map_depth(depth):
        return cm(int((H - 1 - depth) / (H - 1) * 256))

    ax_fl.imshow(fluor_mp, cmap='gray')

    for blob, d in zip(blobs, depth):
        y, x, r = blob
        c = plt.Circle((x, y), r, color=map_depth(d), lw=2, fill=True, alpha=1)
        ax_fl.add_patch(c)

    ax_fl.set_axis_off()

    bins = np.linspace(0, 40 * VOXEL_SIZE[2], 10)
    sm = plt.cm.ScalarMappable(cmap=cm, norm=plt.Normalize(vmin=bins[0], vmax=bins[-1]))
    sm._A = []
    plt.colorbar(sm, ax=ax_fl, label=r"Height from lowest plane ($\mu m$)")

    # Plot: Normalized Histogram
    from matplotlib.ticker import FuncFormatter
    ax_hist = plt.subplot(gs[2:4, 0:3])
    ax_hist.hist(off_cell, bins=bins, normed=True, color='r', alpha=0.7, label="Away from cells")
    ax_hist.hist(on_cell, bins=bins, normed=True, color='b', alpha=0.7, label="On cells")
    ax_hist.yaxis.set_major_formatter(FuncFormatter(lambda y, pos: r'${:0.0f}\%$'.format(y*400)))
    ax_hist.set_xlabel(r"Depth in sample ($\mu m$)")
    ax_hist.set_xlim(bins[[0, -1]])

    # Plot: Relative count
    ax_comp = plt.subplot(gs[2:4, 3])
    ax_comp.bar([0, 1], [len(on_cell), len(off_cell)], color=['b', 'r'], alpha=0.7)
    ax_comp.set_xticks([0, 1])
    ax_comp.set_xticklabels(['On cells', 'Off cells'])
    ax_comp.set_ylabel("Particle count")

    # Subplot lettering
    import string
    for idx, ax in enumerate([ax_bf, ax_fl, ax_hist, ax_comp]):
        ax.text(-0.1, 1.05, string.ascii_uppercase[idx], transform=ax.transAxes,
                size=20, weight='bold')

    plt.tight_layout()

    return fig


if __name__ == "__main__":
     parser = ap.ArgumentParser()
     parser.add_argument('file_paths', type=str, nargs='+',
                         help="List of files to process")
     args = parser.parse_args()

     from os.path import splitext
     from skimage.io import imread
     from tqdm import tqdm

     for path in tqdm(args.file_paths):
         stack = imread(path)
         fig = process_single_stack(stack)
         base, ext = splitext(path)
         figpath = base + '.png'
         fig.savefig(figpath, dpi=300)