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