surface_spheres_3D

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 [110]:

    
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

#from itertools import product, combinations

from pylab import *
%matplotlib inline

create a single sphere to represent a fluorescent bead



In [3]:

    
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

# bead offset 
xf = 10
yf = 10
zf = 10

#construct bead
x = 3 * np.outer(np.cos(u), np.sin(v)) + xf
y = 3 * np.outer(np.sin(u), np.sin(v)) + yf
z = 3 * np.outer(np.ones(np.size(u)), np.cos(v)) +zf
# instantiate plot object
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z,  rstride=4, cstride=4, color='b')
# ax.plot_surface(x, y, z,  rstride=1, cstride=1, color='b')
show()
#plt.show()

plot crosssection of two spheres



In [17]:

    
# bead offset 
xf1 = 0
yf1 = 0
#construct bead
x1 = 3 * np.outer(np.cos(u), np.sin(v)) + xf1
y1 = 3 * np.outer(np.sin(u), np.sin(v)) + yf1

xf2 = 3
yf2 = 3
#construct bead
x2 = 3 * np.outer(np.cos(u), np.sin(v)) + xf2
y2 = 3 * np.outer(np.sin(u), np.sin(v)) + yf2

x = np.concatenate((x1[:,1],x2[:,1]))
y = np.concatenate((y1[:,1],y2[:,1]))

plt.plot(x,y,'.')
#plt.plot(x2[:,1],y2[:,1])
show()

plot two spheres



In [70]:

    
#construct bead
x1 = 3 * np.outer(np.cos(u), np.sin(v)) + 0
y1 = 3 * np.outer(np.sin(u), np.sin(v)) + 0
z1 = 3 * np.outer(np.ones(np.size(u)), np.cos(v)) + 0

x2 = 3 * np.outer(np.cos(u), np.sin(v)) + 10
y2 = 3 * np.outer(np.sin(u), np.sin(v)) + 10
z2 = 3 * np.outer(np.ones(np.size(u)), np.cos(v)) + 10

xc = np.concatenate((x1.T,x2.T))
yc = np.concatenate((y1.T,y2.T))
zc = np.concatenate((z1.T,z2.T))


# instantiate plot object
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(xc, yc, zc,  rstride=4, cstride=4, color='b')
show()

plot to discontinuous spheres



In [71]:

    
# instantiate plot object
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x1, y1, z1,  rstride=4, cstride=4, color='b')
ax.plot_surface(x2, y2, z2,  rstride=4, cstride=4, color='b')
show()

define function for creating many spheres



In [76]:

    
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 [165]:

    
beadSize = .1
nBeads = 20
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')
show()

show plane of image

The plane here will represent the light sheet of the SPIM microscope



In [51]:

    
point  = np.array([1, 2, 3])
#normal = np.array([1, 1, 2])
normal = np.array([1, 0, 1])

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

# create x,y
xx, yy = np.meshgrid(range(10), range(10))
#xx, yy = np.meshgrid(range(100), range(100))

# calculate corresponding z
#z = (-normal[0] * xx - normal[1] * yy - d) * 1. /normal[2]
z2 = (-normal[0]*xx - normal[1]*yy - d) * 1. /normal[2] # convert to real

# plot the surface
plt3d = plt.figure().gca(projection='3d')
plt3d.plot_surface(xx, yy, z2)
show()

show sphere + plane of image



In [87]:

    
### make plane
point  = np.array([1, 2, 3])
#normal = np.array([1, 1, 2])
normal = np.array([0, 1, 1])

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

# create x,y
xx, yy = np.meshgrid(range(10), range(10))
#xx, yy = np.meshgrid(range(100), range(100))

# calculate corresponding z
#z = (-normal[0] * xx - normal[1] * yy - d) * 1. /normal[2]
z2 = (-normal[0]*xx - normal[1]*yy - d) * 1. /normal[2] # convert to real
z2 = z2 + 5

### sphere
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)
# bead offset 
xf = 5
yf = 5
zf = 5

#construct bead
x = 3 * np.outer(np.cos(u), np.sin(v)) + xf
y = 3 * np.outer(np.sin(u), np.sin(v)) + yf
z = 3 * np.outer(np.ones(np.size(u)), np.cos(v)) +zf

# instantiate plot object
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z,  rstride=4, cstride=4, color='b')
ax.plot_surface(xx, yy, z2,color='r',alpha=0.5)
show()
#plt.show()

Plot multiple beads with a plane



In [119]:

    
def make_plane(xoff,yoff,zoff,normal):
    '''
    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)
    # create x,y
    xx, yy = np.meshgrid(range(8), range(8))
    # calculate corresponding z
    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 [150]:

    
beadSize = .1
nBeads = 20
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')

for j in [-3,-1.5,0,1.5,3]:
    normal = np.array([0, 0, 1])
    xx, yy, z2 = make_plane(-4,-4,j,normal)    
    ax.plot_surface(xx, yy, z2,color='r',alpha=0.5)

show()

make verticle planes traversing the beads



In [154]:

    
beadSize = .1
nBeads = 20
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')

for j in [-3,-1.5,0,1.5,3]:
    normal = np.array([0, 0, 1])
    xx, yy, z2 = make_plane(-4,-4,j,normal)    
    ax.plot_surface(xx, z2, yy,color='r',alpha=0.5)

show()



In [160]:

    
### 
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')

for j in [-3,-1.5,0,1.5,3]:
    normal = np.array([1, 0, 1])
    xx, yy, z2 = make_plane(-4,-4,j,normal)    
    ax.plot_surface(xx, z2, yy,color='r',alpha=0.5)
show()

### 
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')

for j in [-3,-1.5,0,1.5,3]:
    normal = np.array([1, 0, 10])
    xx, yy, z2 = make_plane(-4,-4,j,normal)    
    ax.plot_surface(xx, z2, yy,color='r',alpha=0.5)
show()

plot hundreds of beads



In [163]:

    
beadSize = .1
nBeads = 200
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
# add multiple beads
for i in range(nBeads):
    irand = np.random.random_integers(-30,30,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')
show()

sample cube



In [117]:

    
#draw cube
#fig = plt.figure()
#ax = fig.gca(projection='3d')
#r = [-1, 4]
#for s, e in combinations(np.array(list(product(r,r,r))), 2):
#    if np.sum(np.abs(s-e)) == r[1]-r[0]:
#        ax.plot3D(*zip(s,e), color="b")  
#show()
vs = reshape(mgrid[-1:2:2,-1:2:2,-1:2:2].T, (8,3))
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(vs[:,0], vs[:,1], vs[:,2],color='r',alpha=0.5)
show()