Hi Jeff,
The Chimera mask command code can be adapted to calculate the average
thickness of a slab defined by two surfaces. I've attached the adapted
Python code. You select the surfaces then open the script and it
reports the mean thickness, standard deviation of thickness, min/max
thickness and number of thickness samples used. The script contains a
projection axis which should be perpendicular to the slab and a grid
step which defines the grid spacing. The thickness is measured on a
grid perpendicular to the axis and only grid points where the projection
axis intersects the surfaces at least twice are used. You have to edit
the script (at the bottom) to set the projection axis and grid step to
appropriate values. Those values could be automatically set using the
Python code of the inertia ellipsoid script I gave you earlier. Tested
in Chimera version 1.2498 on a spherical surface.
Tom
Jeff wrote:
> hi Tom,
>
> I'd like to compute the 'average distance' between two surfaces to
> calculate the width of a cellular compartment (have attached an image,
> the two surfaces delineating the compartment are colored blue).
>
> one way that has been suggested to me is to calculate the volume
> between the two surfaces, and then just divide by a mean path length
> along one of the surfaces. I know that chimera has a volume
> calculation tool for an isosurface, and was wondering if I could use
> the surfaces you see to do a volume calculation (and, if so, how I
> would go about capping those surfaces to create a closed surface). If
> able to calculate the volume that way, I could just use the path
> length of a contour (which I have) that generated one of the two
> surfaces to derive the "average width".
>
# -----------------------------------------------------------------------------
# Script to calculate average thickness of a slab defined by surfaces.
# Adapted from mask command code for Jeff Triffo, March 18, 2008.
#
# Adapted version of masked_volume() from depthmask.py.
#
# Only uses first two surface layers along projection axis.
#
def average_thickness(projection_axis = (0,0,1),
grid_step = 1.0,
surfaces = None):
# Use selected surfaces if none are passed in.
if surfaces is None:
from Surface import selected_surface_pieces
plist = selected_surface_pieces()
from Mask.depthmask import surface_geometry
surfaces = surface_geometry(plist, tf = None, pad = None)
# Determine transform from vertex coordinates to depth array indices
from Mask.depthmask import orthonormal_frame
fx,fy,fz = orthonormal_frame(projection_axis)
from numpy import array, float32, intc, zeros
tf = array(((fx[0], fx[1], fx[2], 0),
(fy[0], fy[1], fy[2], 0),
(fz[0], fz[1], fz[2], 0)), float32) / grid_step
# Transform vertices to depth array coordinates.
zsurf = []
tcount = 0
for vertices, triangles in surfaces:
varray = vertices.copy()
from Mask.depthmask import apply_transform
apply_transform(tf, varray)
zsurf.append((varray, triangles))
tcount += len(triangles)
if tcount == 0:
return 0,0,0,0,0
from Mask.depthmask import bounding_box
vmin, vmax = bounding_box(zsurf)
voffset = -(vmin - 0.5)
tf[:,3] += voffset
from _contour import shift_vertices
for varray, triangles in zsurf:
shift_vertices(varray, voffset)
from math import ceil, floor
dxsize = int(ceil(vmax[0] - vmin[0] + 1))
dysize = int(ceil(vmax[1] - vmin[1] + 1))
# Create depth arrays
depth = zeros((dysize,dxsize), float32)
tnum = zeros((dysize,dxsize), intc)
depth2 = zeros((dysize,dxsize), float32)
tnum2 = zeros((dysize,dxsize), intc)
# Copy volume to masked volume at masked depth intervals.
max_depth = 1e37
depth.fill(max_depth)
tnum.fill(-1)
from Mask.depthmask import surfaces_z_depth
any = surfaces_z_depth(zsurf, depth, tnum, None, None)
if not any:
return 0,0,0,0,0
depth2.fill(max_depth)
tnum2.fill(-1)
surfaces_z_depth(zsurf, depth2, tnum2, depth, tnum)
# Find thicknesses at grid points that intercept surfaces twice.
dlist = []
for i in range(dysize):
for j in range(dxsize):
if tnum[i,j] != -1 and tnum2[i,j] != -1:
d = (depth2[i,j] - depth[i,j]) * grid_step
dlist.append(d)
if len(dlist) == 0:
return 0,0,0,0,0
from numpy import mean, std
return mean(dlist), std(dlist), min(dlist), max(dlist), len(dlist)
# -----------------------------------------------------------------------------
#
projection_axis = (0,0,1)
grid_step = 1.0
m, sd, mn, mx, n = average_thickness(projection_axis, grid_step)
print ('Slab thickness along axis %s sampled at grid spacing %.3g'
% (str(projection_axis), grid_step))
print 'Mean %.5g, sdev %.5g, min %.5g, max %.5g, samples %d' % (m,sd,mn,mx,n)
from Accelerators.standard_accelerators import show_reply_log
show_reply_log()
