<jupyter_start><jupyter_text>GoalRotate raw images of embryos so that orientation is consistent before vector flow calculation. Process collection from 11-30-2018<jupyter_code>import glob
import tqdm

import tifffile

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage

import bokeh.io
bokeh.io.output_notebook()
import bebi103

from imp import reload
import gbeflow<jupyter_output><empty_output><jupyter_text>Import raw data, first timepoint only<jupyter_code>files = glob.glob('../data/original/20181130-gbe_mutants_brt/*.tif')
files
raw = {}
for f in tqdm.tqdm(files):
    raw[f] = tifffile.imread(f,key=(0,1))<jupyter_output><empty_output><jupyter_text>Load previously calculated alignment parametersSee 20181203-gbe_mutant_brt.ipynb for original calculations<jupyter_code>line = pd.read_csv('20181204-alignment.csv')<jupyter_output><empty_output><jupyter_text>Calculate rotationThe angle of rotation is calculated as follows$$ \theta = \arctan\bigg(\frac{y_2-y_1}{x_2-x_1}\bigg)$$This calculation can be coded using the `np.arctan2`, which has two arguments that correspond to $\Delta y$ and $\Delta x$.<jupyter_code>line = gbeflow.calc_embryo_theta(line)
# line.head()<jupyter_output><empty_output><jupyter_text>Extract the filenumber from the file path so that we can directly compare to the shorter paths present in `files`<jupyter_code>line['fnum'] = line['f'].str.slice(-6,-4).astype('int')
line['fnum']
fnums = [int(f.split('_')[-2]) for f in files]
fnums<jupyter_output><empty_output><jupyter_text>Apply rotation based on $\theta$With $\theta$ calculated, we are now ready to rotate each sample accordingly. Since we cannot load all of the data into memory at the same time, we will currently only rotate the first timepoint to check that it worked. After we have determined all necessary manipulations for each embryo, we will run the actual rotation.<jupyter_code># Dataframe to save first timepoint from each rotate embryo
rot = {}

# List to save bokeh plots
Lp = []

for f,n in tqdm.tqdm(zip(files,fnums)):
    # Extract the theta value for this sample
    theta = line[line['fnum']==n]['theta'].values[0]
    
    # Rotate single image
    rimg = scipy.ndimage.rotate(raw[f][1],theta)
    
    # Save and plot first timepoint
    rot[f] = rimg
    p = bebi103.viz.imshow(rimg,title=f)
    Lp.append(p)
    
bokeh.io.show(bokeh.layouts.gridplot(Lp,ncols=2))<jupyter_output><empty_output><jupyter_text>Correct orientation to be conventional AP DV Create dictionary to save rotation correction functions.<jupyter_code>flips = {}<jupyter_output><empty_output><jupyter_text>Create boolean variable to toggle plotting<jupyter_code>plot=True<jupyter_output><empty_output><jupyter_text>Anterior to the left, dorsal up In 3D, axis 0 flips in time, axis 1 flips vertically, and axis 2 flips horizontally? So a general solution should ge axis=-1 for horizontal and axis=-2 for vertical.<jupyter_code>if plot: bokeh.io.show(Lp[0])
flips[fnums[0]] = lambda img: scipy.ndimage.rotate(img,180,axes=(2,1))
if plot: bokeh.io.show(Lp[1])
flips[fnums[1]] = lambda img: np.flip(img,axis=-2)
if plot: bokeh.io.show(Lp[2])
flips[fnums[2]] = lambda img: np.flip(img,axis=-1)
if plot: bokeh.io.show(Lp[3])
flips[fnums[3]] = lambda img: np.flip(img,axis=-1)
if plot: bokeh.io.show(Lp[4])
flips[fnums[4]] = lambda img: np.flip(img,axis=-1)
if plot: bokeh.io.show(Lp[5])
flips[fnums[5]] = lambda img: np.flip(img,axis=-1)
if plot: bokeh.io.show(Lp[6])
flips[fnums[6]] = lambda img: np.flip(img,axis=-1)
if plot: bokeh.io.show(Lp[7])<jupyter_output><empty_output><jupyter_text>This embryo is oriented correctly, but we'll include a function that doesn't make any changes to avoid possible errors<jupyter_code>flips[fnums[7]] = lambda img: img
if plot: bokeh.io.show(Lp[8])
flips[fnums[8]] = lambda img: img
if plot: bokeh.io.show(Lp[9])
flips[fnums[9]] = lambda img: np.flip(img,axis=-1)
if plot: bokeh.io.show(Lp[10])
flips[fnums[10]] = lambda img: img
if plot: bokeh.io.show(Lp[11])
flips[fnums[11]] = lambda img: img
if plot: bokeh.io.show(Lp[12])
flips[fnums[12]] = lambda img: np.flip(img,axis=-2)
if plot: bokeh.io.show(Lp[13])
flips[fnums[13]] = lambda img: np.flip(img,axis=-2)
if plot: bokeh.io.show(Lp[14])
flips[fnums[14]] = lambda img: np.flip(img,axis=-1)
if plot: bokeh.io.show(Lp[15])
flips[fnums[15]] = lambda img: scipy.ndimage.rotate(img,180,axes=(2,1))
if plot: bokeh.io.show(Lp[16])
flips[fnums[16]] = lambda img: np.flip(img,axis=-2)
if plot: bokeh.io.show(Lp[17])
flips[fnums[17]] = lambda img: img
if plot: bokeh.io.show(Lp[18])
flips[fnums[18]] = lambda img: img
if plot: bokeh.io.show(Lp[19])
flips[fnums[19]] = lambda img: img
if plot: bokeh.io.show(Lp[20])
flips[fnums[20]] = lambda img: img
if plot: bokeh.io.show(Lp[21])
flips[fnums[21]] = lambda img: np.flip(img,axis=-1)
if plot: bokeh.io.show(Lp[22])
flips[fnums[22]] = lambda img: scipy.ndimage.rotate(img,180,axes=(2,1))
if plot: bokeh.io.show(Lp[23])<jupyter_output><empty_output><jupyter_text>Apply theta rotation and flips to each full 3D stackImport full stacks one at a time to avoid memory overload<jupyter_code>line.columns
f[:-4]+'_rot.tif'<jupyter_output><empty_output><jupyter_text>Apply rotation and flips to each full stack, one at a time<jupyter_code>for f,n in tqdm.tqdm(zip(files[19:],fnums[19:])):
    # Select just the brightfield channel, but all timepoints
    raw = tifffile.imread(f)[:,1]
    
    # Extract theta value
    theta = line[line['fnum']==n]['theta']
    
    # Apply theta rotation
    rot = scipy.ndimage.rotate(raw,theta,axes=(2,1))
    
    # Flip according to previously defined functions in flips
    img = flips[n](rot)
    
    tifffile.imsave(f[:-4]+'_rot.tif',data=img)<jupyter_output><empty_output><jupyter_text>Run optical flow function through matlab<jupyter_code>import matlab.engine
eng = matlab.engine.start_matlab()
eng.addpath(r'../matlab',nargout=0)<jupyter_output><empty_output><jupyter_text>Define parameters **Note**: Matlab only accepts complete paths<jupyter_code>movs = glob.glob('/Users/morganschwartz/Code/germband-extension/data/original/20181130-gbe_mutants_brt/*rot.tif')
names = [m[:-4] for m in movs]
BinaryMask = matlab.single([])
scale = 0.5
dt = 1.0
BoxSize = 30
BlurSTD = 1.0 # Must be float
ArrowSize= 5
for f,n in tqdm.tqdm(zip(movs[1:],names[1:])):
    try:
        eng.OpticalFlowOutput(n,f,BinaryMask,
                          BoxSize,BlurSTD,ArrowSize,
                          scale,dt,nargout=0)
    except:
        print('failed',f)<jupyter_output>0it [00:00, ?it/s][A[A[A[A



1it [03:33, 213.62s/it][A[A[A[A
{'fraction_non_alphanumeric': 0.14367816091954022, 'fraction_numerical': 0.031004234724742893, 'mean_word_length': 5.483333333333333, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '207608', 'license': 'MIT', 'n_tokens_mistral': 2720, 'n_tokens_neox': 2534, 'n_words': 585, 'path': 'notebooks/20181206-rotate_20181130.ipynb', 'repo_name': 'msschwartz21/germband-extension', 'size': '6612'}