In [1]:
%pylab inline
import skimage.external.tifffile as tif
# import our ability to read and write MRC files
import Mrc
In [2]:
data_path = "../fixtures/Test SIM Data"
In [3]:
ls "$data_path"
In [4]:
import glob
In [175]:
data = np.array([tif.imread(path) for path in glob.glob(data_path + "/RAWYFP*1_ch*.tif")])
In [9]:
# orientation, z * phase, y, x
data.shape
Out[9]:
In [10]:
# split up axes to be o, z, p, y, x
o, zp, y, x = data.shape
data.shape = o, -1, 5, y, x
data.shape
Out[10]:
In [11]:
# rearrange so now its z, o, p, y, x
data = np.rollaxis(data, 1)
In [13]:
# collapse the z, o, p dimensions
data = data.reshape(-1, y, x)
In [15]:
tif.imsave(data_path + "/rearranged_data.tif", data)
Everything looks good from a data point of view
In [13]:
slice(None, None, -1)
Out[13]:
In [19]:
def rearrange_spim_data(data, dz=1):
"""Assumes data with 4 dimensions"""
# data is orienation, z and p, y, x
o, zp, y, x = data.shape
# reshape so it is o, z, p, y, x
data.shape = o, -1, 5, y, x
# rearrange so now its z, o, p, y, x
data = np.rollaxis(data, 1)
# flip z, collapse the z, o, p dimensions and return
if dz > 0:
# acquired data is in reverse of PSFs
s = slice(None, None, -1)
else:
s = slice(None)
return data[s].reshape(-1, y, x)
In [172]:
a = np.ones((2, 3, 4)) * np.arange(1,3)[:, None, None] * np.arange(1,4)[None, :, None] * np.arange(1,5)[None, None]
In [173]:
a.reshape(2, -1)
Out[173]:
In [6]:
import re
In [7]:
def get_center_wavlength(emission_filter):
razor_edge_wl_re = re.compile(r"(\d+)(?=(?:\s*nm)?(?:\s*)?Razor Edge)")
band_pass_re = re.compile(r"(\d+)\/(\d+)")
short_pass_re = re.compile(r"(\d+)(?:\s*nm\s*)")
razor_edge = razor_edge_wl_re.findall(emission_filter)
band_pass = band_pass_re.findall(emission_filter)
short_pass = short_pass_re.findall(emission_filter)
try:
short_side = int(razor_edge[0])
except IndexError:
short_side = np.nan
try:
bandcenter, bandwidth = np.array(band_pass[0], float)
except IndexError:
bandcenter, bandwidth = np.nan, np.nan
band_short = bandcenter - bandwidth / 2
band_long = bandcenter + bandwidth / 2
try:
long_side = int(short_pass[0])
except IndexError:
long_side = np.nan
short_side = np.nanmin((short_side, band_short))
long_side = np.nanmin((long_side, band_long))
return np.nanmean((short_side, long_side))
In [8]:
def parse_settings(path):
re_dz = re.compile(r"(?:Z Obj Offset\D*\d[^-\d]*)((?:\s*-?\d+(?:\.\d+)?)+)")
re_dz = re.compile(r"Z Obj Offset.*")
re_wl = re.compile(r"Excitation Filter.*")
with open(path) as f:
buffer = f.readlines()
dz_str = re_dz.findall("\n".join(buffer))
wl_str = re_wl.findall("\n".join(buffer))
_, z0, dz, nz = dz_str[0].split("\t")
_, emission_filter, laser, power, exposure = wl_str[0].split("\t")
emission_filter, laser, power, exposure
cwl = get_center_wavlength(emission_filter)
return dict(z0=float(z0), dz=float(dz), nz=int(nz), exposure=float(exposure),
power=float(power), laser=int(laser), center_wl=cwl)
In [9]:
def save_mrc(path, data, wl, dz, dr=0.13, **kwargs):
header = Mrc.makeHdrArray()
# initialize it
# set type and shape
Mrc.init_simple(header, 4, data.shape)
# set wavelength in nm
header.wave = wl
# set number of wavelengths
header.NumWaves = 1
# set dimensions
header.d = dr, dr, dz
Mrc.save(data, path, hdr=header, ifExists='overwrite', **kwargs)
In [139]:
testd = np.array(Mrc.Mrc(data_path + "/VSIM_488_test_data.mrc").data)
testd.shape
Out[139]:
In [148]:
save_mrc(data_path + "/VSIM_488_test_data2.mrc", testd, 519, 0.25)
In [20]:
data = np.array([tif.imread(path) for path in glob.glob(data_path + "/RAWYFP*1_ch*.tif")])
settings = parse_settings(data_path + "/YFP Cell 1_Settings.txt")
save_mrc(data_path + "/rearranged_data.mrc", rearrange_spim_data(data, -1), settings["center_wl"], settings["dz"])
In [5]:
import os
import simrecon_utils as su
In [6]:
for raw_sub in ("Combined/YFP 1.mrc", "Combined/YFP 3.mrc"):
base_kwargs = dict(
nphases=5,
ndirs=3,
angle0= -1.292,
negDangle=False,
na= 0.85,
nimm= 1.0,
zoomfact= 2.0,
background= 100.0,
wiener= 0.001,
fastSIM=True,
otfRA= True,
dampenOrder0=True,
k0searchall=True,
equalizez=True,
preciseapo=True,
gammaApo=0.1,
suppressR=1.5
)
base_kwargs.update(dict(gammaApo=0.3, suppressR=1, wiener=0.005))
raw = data_path + "/" + raw_sub
otf = data_path + "/488 nm Bead 21.25_20170410_202122.mrc"
sim_kwargs = dict(
input_file=raw,
otf_file=otf,
ls=(488/1000)/2/0.81
)
sim_kwargs.update(base_kwargs)
filename = os.path.split(otf)[1]
out_name = sim_kwargs["output_file"] = sim_kwargs["input_file"].replace(".mrc", "_proc.mrc")
%time sim_output = su.simrecon(**sim_kwargs)
with open(out_name.replace(".mrc", ".txt"), "w") as myfile:
myfile.write(str(sim_kwargs))
myfile.write("".join(sim_output))
In [22]:
ls
In [26]:
os.path.isdir(".ipynb_checkpoints")
Out[26]:
In [52]:
os.path.dirname(os.path.join(".ipynb_checkpoints/", ""))
Out[52]:
In [ ]: