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