SPIM_bead_simulation

Larson Hogstrom

Jan 2015

The goal of this code is to construct a system to simulate beads in a 3D SPIM sample and simulate multiple angles of image acquisition

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

define function for creating many spheres and imaging planes



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

Plot multiple beads with horizontal plane



In [48]:

    
beadSize = .1
nBeads = 10
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()

plot planes at multiple angles



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()

Need to match the following physical parameters

Sample dimensions
bead size
image resolution (1344 x 1024)
light-sheet width
spacing of imaging planes - e.g., 51 slices along the sample spaced 6 microns apart
multiple angles of acquisition
scanning to 2/3 of the sample

save tif image of a single plane



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()

save tif image of a single plane



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)



In [ ]:

    
mtrx = np.zeros((3,3))
mtrx[1,2] = 1

igr = mtrx > .5
ilt = mtrx < 1.5

plot beads based on coordinates

assume that bead is only detected in plane if the bead center is in the plane



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()









    












    



[[ 2  1  0]
 [-2  1  0]
 [-2  0 -2]
 [ 1  1 -2]
 [-1 -1 -3]
 [ 1 -3  1]
 [-3  2 -3]
 [ 2  2  2]
 [ 2 -2 -2]
 [ 0  3  2]]



In [ ]:



In [ ]: