Ed Brignole asked about setting a volume contouring level so that the resulting surface encloses a desired volume (cubic Angstroms). Attached is a Python script that does this in two ways. The first way is fast but approximate finding the level such that the number of enclosed grid points times the voxel volume equals the requested volume. Estimating the enclosed volume of the surface this way can give a much different result than actually computing analytically the volume within the surface. The second approach computes surfaces at different contour levels and calculates analytically the volume. It uses bisection to find the correct contour level. This approach can be slow, about 30 seconds for a 256^3 volume (sfv.mrc virus map) with ~2.5e6 surface triangles. Most of that time (80%) is actually computing the enclosed volume. That should be very fast but the code is checking for holes in the surface in a rather slow way. Both methods use the step size (subsampling) set in the volume dialog. Tom # ----------------------------------------------------------------------------- # Find threshold level so that a contour surface encloses a specified volume. # # Method: Estimates how many grid points should lie within # surface then finds the level having that many grid points with higher values. # This method can produce a level where the enclosed volume is very # different from the requested volume. # def level_enclosing_grid_points(v, vol): (x0,y0,z0), (x1,y1,z1) = v.xyz_bounds() tvol = (x1-x0)*(y1-y0)*(z1-z0) r = 1.0 - vol/tvol s = v.matrix_value_statistics() level = s.rank_data_value(r) return level # ----------------------------------------------------------------------------- # Find threshold level so that a contour surface encloses a specified volume. # # Method: Computes surfaces for different levels and resulting volumes # using bisection to determine the correct level. # def level_by_bisection(v, vol, iterations = 10): m = v.matrix() l0, l1 = m.min(), m.max() ijk_step = v.region[2] from numpy import product voxel_volume = product(ijk_step) * product(v.data.step) mvol = vol / voxel_volume return find_level(m, l0, l1, mvol, iterations) # ----------------------------------------------------------------------------- # def find_level(m, l0, l1, mvol, iterations): l = 0.5*(l0+l1) if iterations <= 1: return l from _contour import surface varray, tarray = surface(m, l, cap_faces = True) from _surface import enclosed_volume vl, hole_count = enclosed_volume(varray, tarray) varray = tarray = None if vl == mvol: return l if vl > mvol: l0 = l else: l1 = l return find_level(m, l0, l1, mvol, iterations - 1) # ----------------------------------------------------------------------------- # vol = 3.227e6 # sfv.mrc, threshold level 1000. from chimera import openModels from VolumeViewer import volume_list for v in volume_list(): level1 = level_enclosing_grid_points(v, vol) level2 = level_by_bisection(v, vol) print v.name, level1, level2