analyse muvicyte datasets

segment/muvicyte.py

+from glob import glob
+from droplet_growth import multiwell
+import matplotlib.pyplot as plt
+from multiprocessing import Pool, cpu_count
+import numpy as np
+from nd2tif.transform import Well
+import os
+import pandas as pd
+from PIL import Image
+from scipy.ndimage import gaussian_filter1d
+from segment import seg
+from skimage.color import label2rgb
+from skimage.transform import rotate as _rotate
+
+
+def process_dataset(path, px=.453, interval_h=3, labels=['SOX1', 'BRA']):
+    ch = get_channels(path)
+    print(ch)
+    grouped_paths = get_grouped_paths(path, ch)
+    print(f'{len(grouped_paths)} timepoints')
+    try:
+        p = Pool(cpu_count()-1)
+        rotated_orgs = p.map(rotate_organoid, grouped_paths)
+        w = Well(np.array(rotated_orgs, dtype='float32'), order='tcyx', calibration_um=px)
+        w.save_tif(movie := os.path.join(path, 'aligned.tif'))
+        print(f'saved movie to {movie}')
+        peaks = np.array([get_maxima(r) for r in rotated_orgs[:]]) * px
+        df = pd.DataFrame(data=peaks, columns=labels)
+        df.loc[:, 'time, h'] = (time := np.arange(len(peaks)) * interval_h)
+        df.to_csv(csv := os.path.join(path, 'peaks.csv'))
+        print(f'saved peak coordiantes to {csv}')
+
+        fig, ax = plt.subplots(dpi=150)
+        plt.plot(time, peaks[:,0], label=labels[0])
+        plt.plot(time, peaks[:,1], label=labels[1])
+        plt.legend()
+        plt.ylabel('peak distance from center, μm')
+        plt.xlabel('time, h')
+        plt.savefig(plot := os.path.join(path, 'plot.png'))
+
+        print(f'saved plot to {plot}')
+    except Exception as e:
+        print('ERROR', e.args)
+    finally:
+        p.close()
+    
+    return 
+
+def get_channels(path):
+    return sorted([os.path.split(os.path.split(p)[0])[-1] for p in glob(os.path.join(path, '*/'))])
+
+
+def get_grouped_paths(path, channels):
+    f_list = {ch: sorted(glob(os.path.join(path, ch, '*.JPG'))) for ch in channels}
+    grouped_list = [{ch: f_list[ch][i] for ch in channels} for i in range(len(f_list[channels[0]]))]
+    return grouped_list
+
+def get_maxima(rotated_stack, channels=[1,2]):
+    mask = rotated_stack[3]
+    ch1,ch2 = [rotated_stack[i] * mask for i in channels]
+    cx = ch1.shape[1] // 2
+    max_proj_x = [ch.max(axis=0) for ch in (ch1, ch2)]
+    peaks_x = [_get_peak_position(mp, min_val=2, default_position=cx, smooth=10) - cx for mp in max_proj_x]
+    return peaks_x
+
+def _get_peak_position(proj, min_val=2, default_position=0, smooth=10):
+    if proj.max() < min_val:
+        return default_position
+    return np.argmax(gaussian_filter1d(proj, smooth))
+
+def segment_bf(well, thr=0.5, smooth=10, erode=10, fill=True, plot=False):
+    '''
+    Serments input 2d array using thresholded gradient with filling
+    Returns SegmentedImage object
+    '''
+    grad = multiwell.get_2d_gradient(well)
+    sm = multiwell.gaussian_filter(grad, smooth)
+#     sm = multiwell.gaussian_filter(well, smooth)
+    
+    regions = sm > thr
+    
+    if fill:
+        regions = multiwell.binary_fill_holes(regions)
+    
+    if erode:
+        regions = multiwell.binary_erosion(regions, iterations=erode)
+    
+    labels, n_labels = multiwell.label(regions)
+    # print(f'{n_labels} regions')
+    if plot:
+        fig, ax = multiwell.plt.subplots(1,2)
+        ax[0].imshow(sm, cmap='gray')
+        ax[1].imshow(labels)
+        plt.show()
+        
+    return labels
+
+def get_biggest_region(labels):
+    props = seg.regionprops(labels)
+    return labels == props[np.argmax([p.area for p in props])].label
+
+
+def get_contour(mask, size=1, overlay_image=None):
+    contour = np.logical_xor(mask, seg.binary_erosion(mask, np.ones((size, size))))
+    if overlay_image is not None:
+        plt.imshow(label2rgb(contour, overlay_image, bg_label=0))
+        plt.show()
+    return contour
+
+def segment_jpg(path, thr=.5, smooth=10, plot=True, contour=20):
+    img = np.array(Image.open(path))
+    if smooth > 0:
+        sm_img = seg.gaussian_filter(img, smooth)
+    else:
+        sm_img = img.copy()
+    mask = segment_bf(sm_img, thr=thr, plot=plot)
+    big = get_biggest_region(mask)
+    if plot:
+        plt.imshow(big)
+        plt.show()
+        contour = get_contour(big, size=contour, overlay_image=img)
+    return img, big
+
+def sub_bg(img:np.ndarray, smooth=50):
+    img = img.astype('float32')
+    if smooth:
+        return img - seg.gaussian_filter(img, smooth)
+    return img - img.min()
+
+
+def get_max_pos2D(a:np.ndarray):
+    return np.unravel_index(np.argmax(a, axis=None), a.shape)
+
+
+def rotate_organoid(grouped_path, seg_ch='BRIGHT', orient='GFP', plot=False):
+    bf, mask = segment_jpg(grouped_path[seg_ch], plot=plot)
+    props = seg.regionprops(seg.label(mask))
+    rotated = [rotate(sub_bg(np.array(Image.open(path)), 150), props[0]) for path in grouped_path.values()] + \
+                [rotate(mask, props[0]).astype('uint8')]
+    
+    # flip 180 deg if needed to get `orient` channel to the left
+    index = list(grouped_path.keys()).index(orient)
+    max_pos = get_max_pos2D(rotated[index])
+    if max_pos[1] > rotated[0].shape[1] // 2:
+        rotated1 = [_rotate(r, 180, preserve_range=True) for r in rotated]
+        rotated = rotated1
+    
+    if plot:
+        fig, ax = plt.subplots(ncols=len(rotated), figsize=(10,3), dpi=150)
+        for a, img in zip(ax, rotated):
+            im1 = a.imshow(img)
+            fig.colorbar(im1, ax=a)
+    return np.array(rotated)
+
+
+def move_to_center(img, centroid):
+    shape = img.shape
+    cyx = np.array(shape, dtype='int') // 2
+    dy, dx = np.array(centroid, dtype='int') - cyx
+    padded = np.pad(img, ((abs(dy), abs(dy)), (abs(dx), abs(dx))), mode='edge')
+    return padded[abs(dy) + dy : -1 + dy - abs(dy), abs(dx) + dx : -1 + dx - abs(dx)]
+
+
+def rotate(img, prop, add=0):
+    centered = move_to_center(img, prop.centroid)
+    angle = -prop.orientation / np.pi * 180 - 90 + add
+    return _rotate(centered, angle, preserve_range=True, resize=False, center=None)