In [1]:
%matplotlib inline
In [2]:
from numpy import array
import matplotlib.pyplot as plt
import seaborn as sns
In [3]:
import thunder
from showit import image, tile
import matplotlib.animation as animation
In [4]:
from os.path import join, exists
from os import mkdir, makedirs
In [5]:
from numpy import save
In [6]:
from skimage.io import imsave, imread
In [7]:
from regional import many
from numpy import random
In [8]:
def normalize(oim):
# normalizes 3D image across first axis
assert oim.ndim == 3
means = oim.mean(axis=(1, 2), dtype='float32')
return array([oim[i]/means[i]/4 for i in range(oim.shape[0])]).clip(0, 1)
In [9]:
directory = '/tier2/freeman/Nick/lfov.calibration'
In [10]:
key = '2016-04-16-L23'
name = 'anm-0330549'
In [11]:
path = join(directory, 'reprocessed', name, key)
print exists(path)
In [12]:
savepath = join(path, 'sources')
if not exists(savepath):
makedirs(savepath)
In [15]:
data = thunder.images.frombinary(join(path, 'patch'), engine=sc, npartitions = 1000)
In [21]:
data = data[:,:200,:200]
In [23]:
data.compute()
In [16]:
from numpy import arange, polyfit, polyval
In [24]:
def detrend(y, order=5):
"""
Detrend series data with linear or nonlinear detrending.
Preserve intercept so that subsequent operations can adjust the baseline.
Parameters
----------
method : str, optional, default = 'linear'
Detrending method
order : int, optional, default = 5
Order of polynomial, for non-linear detrending only
"""
x = arange(len(y))
p = polyfit(x, y, order)
p[-1] = 0
yy = polyval(p, x)
return y - yy
In [25]:
detrended = data.map_as_series(lambda x: detrend(x))
In [26]:
detrended.cache()
detrended.compute()
In [29]:
#detrended.tobinary(join(path, 'patch'))
In [30]:
detrended
Out[30]:
In [32]:
mean = detrended.mean().toarray()
In [34]:
image(mean, clim=(0, 3.5*mean.mean()), size=12)
Out[34]:
In [57]:
meanG = detrended.uniform_filter(2).mean().toarray()
In [58]:
image(meanG, clim=(0, 3.5*meanG.mean()), size=12)
Out[58]:
In [16]:
localcorr2 = data.localcorr(2).astype('float32')
In [18]:
image(localcorr2, clim=(0, 2*localcorr2.mean()), size=12)
Out[18]:
In [ ]:
In [35]:
localcorr4 = data.localcorr(4).astype('float32')
In [36]:
image(localcorr4, clim=(0, 3.5*localcorr4.mean()), size=12)
Out[36]:
In [38]:
localcorr8 = detrended.localcorr(8).astype('float32')
In [39]:
image(localcorr8, clim=(0, 3.5*localcorr8.mean()), size=12)
Out[39]:
In [59]:
localcorr1G8 = detrended.uniform_filter(2).localcorr(8).astype('float32')
In [60]:
image(localcorr1G8, clim=(0, 3.5*localcorr1G8.mean()), size=12)
Out[60]:
In [52]:
localcorr16 = detrended.localcorr(16).astype('float32')
In [53]:
image(localcorr16, clim=(0, 3.5*localcorr16.mean()), size=12)
Out[53]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [17]:
data = imread(join(path, 'summary', 'localcorr.tif'), plugin='tifffile')
#data = data.astype('float32')/255
In [16]:
#norm = normalize(data)
In [32]:
img = data[:200,:200]
In [ ]:
In [33]:
image(img, size=12)
Out[33]:
In [19]:
from skimage.filters import gaussian_filter
from skimage.feature import peak_local_max
In [20]:
from scipy.optimize import curve_fit
from scipy.stats import multivariate_normal
In [21]:
def inside(im, center, window):
inside = [(center[i]-window[i] >= 0) & (center[i]+window[i] <= im.shape[i]) for i in range(0,2)]
inside = np.all(inside)
return inside
In [22]:
def volume(im, center, window):
if np.all(inside(im, center, window)):
volume = im[(center[0]-window[0]):(center[0]+window[0]), (center[1]-window[1]):(center[1]+window[1])]
volume = volume.astype('float64')
baseline = volume[[0,-1],[0,-1]].mean()
volume = volume - baseline
volume = volume/volume.max()
return volume
In [24]:
imgT = thunder.images.fromarray(localcorr2)
In [25]:
smoothed = imgT.uniform_filter(4)
In [26]:
image(smoothed, size=12)
Out[26]:
In [27]:
import numpy as np
In [28]:
threshold_abs = np.percentile(smoothed.toarray().flatten(),70)
print threshold_abs
threshold_rel = np.percentile(smoothed.toarray().flatten(),80) - threshold_abs
print threshold_rel
In [29]:
from skimage.feature import blob_log, blob_doh, blob_dog
In [31]:
centers = blob_log(localcorr2, min_sigma=.2, max_sigma=4, num_sigma=20, threshold=.05)
len(centers)
Out[31]:
In [215]:
centers = blob_dog(localcorr16, min_sigma=.3, max_sigma=3, sigma_ratio=1.05, threshold=.01)
len(centers)
Out[215]:
In [140]:
def findBeads(img):
centers = peak_local_max(img, min_distance=1, threshold_rel=threshold_rel, threshold_abs=threshold_abs, exclude_border=False)
return centers
In [141]:
centers = findBeads(smoothed.toarray())
len(centers)
Out[141]:
In [34]:
image(localcorr2, size=12)
plt.plot(centers[:,1], centers[:,0], 'g.', ms=20);
plt.xlim([0, localcorr2.shape[1]])
plt.ylim([localcorr2.shape[0], 0])
Out[34]:
In [219]:
[centers[20,1], centers[20,0]]
Out[219]:
In [239]:
blurred = detrended.uniform_filter(16)
In [242]:
trace = blurred[:,50,33].toarray()
In [243]:
plt.plot(trace)
Out[243]:
In [232]:
records = detrended.toseries()
In [245]:
corIndvB = records.correlate(trace).toarray()
In [244]:
corIndv = X
In [240]:
image(corIndv, size=12)
Out[240]:
In [246]:
image(corIndvB, size=12)
Out[246]:
In [ ]:
In [ ]:
In [ ]:
In [48]:
def getCenters(im, options):
window = [options['windowUm'][0]*options['pxPerUmAx'], options['windowUm'][1]*options['pxPerUmLat']]
window = [round(x) for x in window]
centers, smoothed = findBeads(im, window)
print len(centers)
beads = [volume(im, x, window) for x in centers]
maxima = [im[x[0], x[1]] for x in centers]
return beads, maxima, centers, smoothed
In [49]:
FOVumLat = 590.0
FOVpxLat = 512.0 # 512
pxPerUmLat = FOVpxLat/FOVumLat
pxPerUmAx = .05 # 2.0
wavelength = 970.0
NA = 0.6
windowUm = [2, 2]
options = {'FOVumLat':FOVumLat, 'FOVpxLat':FOVpxLat, 'pxPerUmLat':FOVpxLat/FOVumLat, 'pxPerUmAx':pxPerUmAx, 'wavelength':970.0, 'NA':0.6, 'windowUm':windowUm}
In [50]:
import numpy as np
In [51]:
beads, maxima, centers, smoothed = getCenters(patch, options)
In [ ]:
In [ ]:
In [18]:
img = data
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [56]:
plane = 0
img = norm[plane]
In [205]:
with open(join(path, 'sources', 'sources-%04d.tif' % plane)) as fid:
sources = many([x['coordinates'] for x in json.load(fid)])
In [19]:
viz = lgn.imagepoly(img)
viz
In [ ]:
X = viz.points(z=plane)
In [ ]:
X
In [ ]:
In [58]:
sources = many(viz.points(z=plane))
In [65]:
from numpy import maximum, tile
In [77]:
base.shape
Out[77]:
In [74]:
sources[0]
Out[74]:
In [78]:
base = tile(norm,(3,1,1,1)).transpose(1,2,3,0).transpose(1,2,0,3)
masks = sources.mask((512,512,4), background='black', fill='blue', stroke='orange')
blend = maximum(base, masks)
In [81]:
fig = plt.figure(figsize=[12,12])
ax = plt.axes()
im = image(blend[:,:,0], ax=ax)
plt.xlim([0, blend.shape[1]]);
plt.ylim([blend.shape[0], 0]);
#for s in range(sources.count):
# plt.annotate(s=str(s), xy=[sources.center[s][1],sources.center[s][0]], color='w', size = 20);
In [82]:
imsave(join(path, 'sources', 'sources.tif'), (255*blend).astype('uint8'), plugin='tifffile', photometric='rgb')
#imsave(join(path, 'sources', 'sources.tif'), (255*blend).astype('uint8'), plugin='tifffile', photometric='rgb')
In [203]:
imsave(join(path, 'sources', 'sources-%04d.tif' % plane), (255*blend).astype('uint8'), plugin='tifffile', photometric='rgb')
In [207]:
###multi dim nature of this ...
In [ ]:
foo = sources.masks((512,512,4), color='blue', base=X)
image(foo[:,:,1,:])
In [227]:
from numpy import concatenate
In [231]:
x = many([[concatenate((x, [plane])) for x in roi.coordinates] for roi in sources])
In [235]:
if norm.ndim == 2:
base = tile(norm,(3,1,1)).transpose(1,2,0)
masks = sources.mask(norm.shape, background='black', fill='blue', stroke='orange')
blend = maximum(base, masks)
else:
base = tile(norm,(3,1,1,1)).transpose(1,2,3,0)
masks = [sources.mask(norm.shape, background='black', fill='blue', stroke='orange')
In [ ]:
def mask(sources, shape, **kwargs):
if len(shape) == 3:
else:
return sources.mask(shape, **kwargs)
In [1]:
import regional
In [2]:
regional.__version__
Out[2]:
In [19]:
a = [1, 18, 35]
b = [2, 20, 40]
In [5]:
import numpy as np
In [20]:
x = np.stack((a, b), axis = 1)
In [21]:
x
Out[21]:
In [23]:
x[:,0]
Out[23]:
In [ ]:
In [21]:
from sima.misc.imagej import read_imagej_roi_zip
In [25]:
X = read_imagej_roi_zip('/Users/sofroniewn/github/regional/RoiSet.zip')
In [33]:
X[0]
Out[33]:
In [38]:
Y[0]
Out[38]:
In [35]:
Y = read_roi_zip('/Users/sofroniewn/github/regional/RoiSet.zip')
In [34]:
# Copyright: Luis Pedro Coelho <luis@luispedro.org>, 2012
# License: MIT
import numpy as np
def read_roi(fileobj):
'''
points = read_roi(fileobj)
Read ImageJ's ROI format
'''
# This is based on:
# http://rsbweb.nih.gov/ij/developer/source/ij/io/RoiDecoder.java.html
# http://rsbweb.nih.gov/ij/developer/source/ij/io/RoiEncoder.java.html
SPLINE_FIT = 1
DOUBLE_HEADED = 2
OUTLINE = 4
OVERLAY_LABELS = 8
OVERLAY_NAMES = 16
OVERLAY_BACKGROUNDS = 32
OVERLAY_BOLD = 64
SUB_PIXEL_RESOLUTION = 128
DRAW_OFFSET = 256
pos = [4]
def get8():
pos[0] += 1
s = fileobj.read(1)
if not s:
raise IOError('readroi: Unexpected EOF')
return ord(s)
def get16():
b0 = get8()
b1 = get8()
return (b0 << 8) | b1
def get32():
s0 = get16()
s1 = get16()
return (s0 << 16) | s1
def getfloat():
v = np.int32(get32())
return v.view(np.float32)
magic = fileobj.read(4)
if magic != 'Iout':
raise IOError('Magic number not found')
version = get16()
# It seems that the roi type field occupies 2 Bytes, but only one is used
roi_type = get8()
# Discard second Byte:
get8()
if not (0 <= roi_type < 11):
raise ValueError('roireader: ROI type %s not supported' % roi_type)
if roi_type != 7:
raise ValueError('roireader: ROI type %s not supported (!= 7)' % roi_type)
top = get16()
left = get16()
bottom = get16()
right = get16()
n_coordinates = get16()
x1 = getfloat()
y1 = getfloat()
x2 = getfloat()
y2 = getfloat()
stroke_width = get16()
shape_roi_size = get32()
stroke_color = get32()
fill_color = get32()
subtype = get16()
if subtype != 0:
raise ValueError('roireader: ROI subtype %s not supported (!= 0)' % subtype)
options = get16()
arrow_style = get8()
arrow_head_size = get8()
rect_arc_size = get16()
position = get32()
header2offset = get32()
if options & SUB_PIXEL_RESOLUTION:
getc = getfloat
points = np.empty((n_coordinates, 2), dtype=np.float32)
else:
getc = get16
points = np.empty((n_coordinates, 2), dtype=np.int16)
points[:,1] = [getc() for i in xrange(n_coordinates)]
points[:,0] = [getc() for i in xrange(n_coordinates)]
points[:,1] += left
points[:,0] += top
points -= 1
return points
def read_roi_zip(fname):
import zipfile
with zipfile.ZipFile(fname) as zf:
return [read_roi(zf.open(n))
for n in zf.namelist()]
In [ ]: