Larson Hogstrom
Jan 2015
SPIM samples are typically embedded with micrometer-sized fluorescent beads in a 1 % Low Melting Point (LMP) Agarose in deionised water. The fluorescent beads serve as fuducial markers during registration of multiview data. Fluorescent beads are purchased commercially at a concentration of around 5*10^13 particles/ml and are diluted at a ratio of 1:10^6 in the agrose. This will yield approximately 400 beads per sample.
In this simulation, the microscope's light sheet is represented by a plane through the sample space. A series of planes are generated at multiple angles to reflect SPIM data acquisition. Tiff images will be generated along each of these planes. This will serve as simulated data to be processed in by existing bead-based registration systems such as those available for Fiji.
In [212]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pylab import *
import Image
#from tifffile import imsave
%matplotlib inline
In [4]:
def addBead(r,xf,yf,zf):
'''
r - sphere size
xf - x offset
yf - y offset
zf - z offset
'''
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
#construct bead
x = r * np.outer(np.cos(u), np.sin(v)) + xf
y = r * np.outer(np.sin(u), np.sin(v)) + yf
z = r * np.outer(np.ones(np.size(u)), np.cos(v)) +zf
return x, y, z
In [41]:
def make_plane(xoff,yoff,zoff,normal,grid_range,ref_plane='z'):
'''
xoff - x offset
yoff - y offset
zoff - z offset
normal - define plane, eg: np.array([0, 0, 1])
'''
point = np.array([0, 0, 0])
# a plane is a*x+b*y+c*z+d=0
# [a,b,c] is the normal. Thus, we have to calculate
# d and we're set
d = -point.dot(normal)
### x_ref
if ref_plane =='x':
# xy are the base grid
z2, yy = np.meshgrid(grid_range, grid_range)
xx = (-normal[2]*z2 - normal[1]*yy - d) * 1. /normal[0] # convert to real
### y_ref
if ref_plane =='y':
# xz are the base grid
xx, z2 = np.meshgrid(grid_range, grid_range)
yy = (-normal[0]*xx - normal[2]*z2 - d) * 1. /normal[1] # convert to real
### z_ref
if ref_plane =='z':
# xy are the base grid
xx, yy = np.meshgrid(grid_range, grid_range)
z2 = (-normal[0]*xx - normal[1]*yy - d) * 1. /normal[2] # convert to real
#apply offset
xx = xx + xoff
yy = yy + yoff
zz = z2 + zoff
return xx, yy, zz
In [213]:
beadSize = .1
nBeads = 50
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# add multiple beads
for i in range(nBeads):
irand = np.random.random_integers(-3,3,size=3)
x,y,z = addBead(beadSize, irand[0], irand[1], irand[2])
ax.plot_surface(x, y, z, rstride=4, cstride=4, color='b')
grid_range = np.arange(-5,5)
for j in [-3,-1.5,0,1.5,3]:
#for j in [0]:
### horizontal
#normal = np.array([0, 0, 1])
#xx, yy, z2 = make_plane(0,0,j,normal,grid_range,ref_plane='z')
### orthogonal to x
#normal = np.array([1, 0, 0])
#xx, yy, z2 = make_plane(j,0,0,normal,grid_range,ref_plane='x')
### orthogonal to y
normal = np.array([0, 1, 0])
xx, yy, z2 = make_plane(0,j,0,normal,grid_range,ref_plane='y')
### orthogonal to (1,1,0)
#normal = np.array([1, 1, 0])
#xx, yy, z2 = make_plane(0,j,0,normal,grid_range,ref_plane='y')
ax.plot_surface(xx, yy, z2, color='r',alpha=0.5)
show()
In [57]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
grid_range = np.arange(-5,5)
for j in [-3,-1.5,0,1.5,3]:
#for j in [0]:
### orthogonal to x
normal = np.array([1, 0, 0])
xx, yy, z2 = make_plane(j,0,0,normal,grid_range,ref_plane='x')
ax.plot_surface(xx, yy, z2, color='r',alpha=0.5)
show()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for j in [-3,-1.5,0,1.5,3]:
#for j in [0]:
### orthogonal to (2,1,0)
normal = np.array([2, 1, 0])
xx, yy, z2 = make_plane(2*j,j,0,normal,grid_range,ref_plane='y')
ax.plot_surface(xx, yy, z2, color='r',alpha=0.5)
show()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for j in [-3,-1.5,0,1.5,3]:
#for j in [0]:
### orthogonal to (1,1,0)
normal = np.array([1, 1, 0])
xx, yy, z2 = make_plane(j,j,0,normal,grid_range,ref_plane='y')
ax.plot_surface(xx, yy, z2, color='r',alpha=0.5)
show()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for j in [-3,-1.5,0,1.5,3]:
#for j in [0]:
### orthogonal to (1,2,0)
normal = np.array([1, 2, 0])
xx, yy, z2 = make_plane(j,2*j,0,normal,grid_range,ref_plane='y')
ax.plot_surface(xx, yy, z2, color='r',alpha=0.5)
show()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
for j in [-3,-1.5,0,1.5,3]:
#for j in [0]:
### orthogonal to y
normal = np.array([0, 1, 0])
xx, yy, z2 = make_plane(0,j,0,normal,grid_range,ref_plane='y')
ax.plot_surface(xx, yy, z2, color='r',alpha=0.5)
show()
In [150]:
fig = plt.figure(facecolor="k")
ax = fig.add_subplot(1,1,1, axisbg="k")
ax.plot(x[:,99],y[:,99],'.',c='w') #, backgroundcolor='b'
ax.set_axis_bgcolor('r')
plt.axis('off')
fig.savefig(fname,facecolor='k',dpi=200)
#show()
In [184]:
wkdir = '/Users/hogstrom/Dropbox (MIT)/Neuron_data/SPIM_simulation/tiff'
beadSize = .1
nBeads = 10
bcMtrx = np.random.random_integers(-3,3,size=(nBeads,3)) #bead coordinate matrix
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# add multiple beads
x = np.zeros((0,100))
y = np.zeros((0,100))
z = np.zeros((0,100))
for i in range(nBeads):
x1,y1,z1 = addBead(beadSize, bcMtrx[i,0], bcMtrx[i,1], bcMtrx[i,2])
x = np.concatenate((x,x1))
y = np.concatenate((y,y1))
z = np.concatenate((z,z1))
ax.plot_surface(x1, y1, z1, rstride=4, cstride=4, color='b')
show()
In [183]:
#pick z slice: 0-.1
slice_thickness = .05
for j in [-3,-1.5,0,1.5,3]:
igr = z > j
ilt = z < j + slice_thickness
irange = igr * ilt
#z[irange]
fig = plt.figure(facecolor="k")
ax = fig.add_subplot(1,1,1, axisbg="k")
ax.plot(x[irange],y[irange],'.',c='w') #, backgroundcolor='b'
ax.set_axis_bgcolor('r')
plt.axis('off')
show()
#fig.savefig(fname,facecolor='k',dpi=200)
#fig = plt.figure()
#ax = fig.add_subplot(111, projection='3d')
#ax.plot_surface(x, y, z, color='r',alpha=0.5)
#show()
In [182]:
j=3
igr = z > j
ilt = z < j + slice_thickness
irange = igr * ilt
#z[irange]
fig = plt.figure(facecolor="k")
ax = fig.add_subplot(1,1,1, axisbg="k")
ax.plot(x[irange],y[irange],'.',c='w') #, backgroundcolor='b'
ax.set_axis_bgcolor('r')
plt.axis('off')
show()
close()
len(irange)
Out[182]:
In [ ]:
mtrx = np.zeros((3,3))
mtrx[1,2] = 1
igr = mtrx > .5
ilt = mtrx < 1.5
In [210]:
wkdir = '/Users/hogstrom/Dropbox (MIT)/Neuron_data/SPIM_simulation/tiff'
beadSize = .1
nBeads = 10
bcMtrx = np.random.random_integers(-3,3,size=(nBeads,3)) #bead coordinate matrix
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# add multiple beads
x = np.zeros((0,100))
y = np.zeros((0,100))
z = np.zeros((0,100))
for i in range(nBeads):
x1,y1,z1 = addBead(beadSize, bcMtrx[i,0], bcMtrx[i,1], bcMtrx[i,2])
x = np.concatenate((x,x1))
y = np.concatenate((y,y1))
z = np.concatenate((z,z1))
ax.plot_surface(x1, y1, z1, rstride=4, cstride=4, color='b')
show()
print bcMtrx
slice_thickness = .05
for j in [-3,-2,-1,0,1,2,3]:
igr = bcMtrx[:,2] >= j
ilt = bcMtrx[:,2] <= j + slice_thickness
irange = igr * ilt
#z[irange]
fig = plt.figure(facecolor="k")
ax = fig.add_subplot(1,1,1, axisbg="k")
ax.plot(bcMtrx[irange,0],bcMtrx[irange,1],'o',c='w') #, backgroundcolor='b'
ax.axis([-4, 4, -4, 4])
ax.set_axis_bgcolor('r')
plt.axis('off')
fname = wkdir + '/test_' + str(j) + '.tif'
fig.savefig(fname,facecolor='k',dpi=200)
show()
plt.close()
In [ ]:
In [ ]: