text: import numpy as np
  import tensorflow as tf
  import tensorflow.keras.backend as K
  import itertools
  
  
  def face_normals(vertices, faces):
      edges1 = tf.gather(vertices, faces[:, 1]) - tf.gather(vertices, faces[:, 0])
      edges2 = tf.gather(vertices, faces[:, 2]) - tf.gather(vertices, faces[:, 0])
      return tf.linalg.l2_normalize(tf.cross(edges1, edges2), axis=1)
  
  
  def vertex_normals(vertices, faces, vfaces):
      fnorm = face_normals(vertices, faces)
  
      max_faces = max([len(x) for x in vfaces])
      nvertices = len(vfaces)
  
      inds = np.zeros((nvertices, max_faces), dtype='int32')
      mask = np.zeros((nvertices, max_faces, 3), dtype='int16')
      for v, f in enumerate(vfaces):
          inds[v, :len(f)] = f
          mask[v, :len(f), :] = 1
  
      vfnorm = tf.gather(fnorm, inds) * mask
      avg = tf.reduce_sum(vfnorm, axis=1) / tf.cast(tf.reduce_sum(mask, axis=1), dtype=tf.float32)
      return tf.linalg.l2_normalize(avg, axis=1)
  
  
  def prod_n(lst):
      '''From Adrian...'''
      prod = lst[0]
      for p in lst[1:]:
          prod *= p
      return prod
  
  
  def sub2ind(siz, subs, **kwargs):
      '''From Adrian...'''
      # subs is a list
      assert len(siz) == len(subs), 'found inconsistent siz and subs: %d %d' % (len(siz), len(subs))
  
      k = np.cumprod(siz[::-1])
  
      ndx = subs[-1]
      for i, v in enumerate(subs[:-1][::-1]):
          ndx = ndx + v * k[i]
  
      return ndx
  
  
  def interp(vol, loc, method='linear'):
      '''From Adrian...'''
  
      if isinstance(loc, (list, tuple)):
          loc = tf.stack(loc, -1)
      nb_dims = loc.shape[-1]
  
      if not nb_dims.value:
          raise Exception("Loc dimension is None")
  
      if len(vol.shape) not in [nb_dims, nb_dims+1]:
          raise Exception("Number of loc Tensors %d does not match volume dimension %d" % (int(nb_dims), len(vol.shape[:-1])))
  
      if nb_dims > len(vol.shape):
          raise Exception("Loc dimension %d does not match volume dimension %d" % (int(nb_dims), len(vol.shape)))
  
      if vol.shape.ndims == nb_dims:
          vol = K.expand_dims(vol, -1)
  
      # flatten and float location Tensors
      loc = tf.cast(loc, 'float32')
      
      if isinstance(vol.shape, (tf.Dimension, tf.TensorShape)):
          volshape = vol.shape.as_list()
      else:
          volshape = vol.shape
  
      # interpolate
      if method == 'linear':
          loc0 = tf.floor(loc)
  
          # clip values
          max_loc = [d - 1 for d in vol.get_shape().as_list()]
          clipped_loc = [tf.clip_by_value(loc[...,d], 0, max_loc[d]) for d in range(nb_dims)]
          loc0lst = [tf.clip_by_value(loc0[...,d], 0, max_loc[d]) for d in range(nb_dims)]
  
          # get other end of point cube
          loc1 = [tf.clip_by_value(loc0lst[d] + 1, 0, max_loc[d]) for d in range(nb_dims)]
          locs = [[tf.cast(f, 'int32') for f in loc0lst], [tf.cast(f, 'int32') for f in loc1]]
  
          # compute the difference between the upper value and the original value
          # differences are basically 1 - (pt - floor(pt))
          #   because: floor(pt) + 1 - pt = 1 + (floor(pt) - pt) = 1 - (pt - floor(pt))
          diff_loc1 = [loc1[d] - clipped_loc[d] for d in range(nb_dims)]
          diff_loc0 = [1 - d for d in diff_loc1]
          weights_loc = [diff_loc1, diff_loc0] # note reverse ordering since weights are inverse of diff.
  
          # go through all the cube corners, indexed by a ND binary vector 
          # e.g. [0, 0] means this "first" corner in a 2-D "cube"
          cube_pts = list(itertools.product([0, 1], repeat=nb_dims))
          interp_vol = 0
          
          for c in cube_pts:
              
              # get nd values
              # note re: indices above volumes via https://github.com/tensorflow/tensorflow/issues/15091
              #   It works on GPU because we do not perform index validation checking on GPU -- it's too
              #   expensive. Instead we fill the output with zero for the corresponding value. The CPU
              #   version caught the bad index and returned the appropriate error.
              subs = [locs[c[d]][d] for d in range(nb_dims)]
  
              # tf stacking is slow for large volumes, so we will use sub2ind and use single indexing.
              # indices = tf.stack(subs, axis=-1)
              # vol_val = tf.gather_nd(vol, indices)
              # faster way to gather than gather_nd, because the latter needs tf.stack which is slow :(
              idx = sub2ind(vol.shape[:-1], subs)
              vol_val = tf.gather(tf.reshape(vol, [-1, volshape[-1]]), idx)
  
              # get the weight of this cube_pt based on the distance
              # if c[d] is 0 --> want weight = 1 - (pt - floor[pt]) = diff_loc1
              # if c[d] is 1 --> want weight = pt - floor[pt] = diff_loc0
              wts_lst = [weights_loc[c[d]][d] for d in range(nb_dims)]
              # tf stacking is slow, we we will use prod_n()
              # wlm = tf.stack(wts_lst, axis=0)
              # wt = tf.reduce_prod(wlm, axis=0)
              wt = prod_n(wts_lst)
              wt = K.expand_dims(wt, -1)
              
              # compute final weighted value for each cube corner
              interp_vol += wt * vol_val
          
      else:
          assert method == 'nearest'
          roundloc = tf.cast(tf.round(loc), 'int32')
  
          # clip values
          max_loc = [tf.cast(d - 1, 'int32') for d in vol.shape]
          roundloc = [tf.clip_by_value(roundloc[...,d], 0, max_loc[d]) for d in range(nb_dims)]
  
          # get values
          # tf stacking is slow. replace with gather
          # roundloc = tf.stack(roundloc, axis=-1)
          # interp_vol = tf.gather_nd(vol, roundloc)
          idx = sub2ind(vol.shape[:-1], roundloc)
          interp_vol = tf.gather(tf.reshape(vol, [-1, vol.shape[-1]]), idx) 
  
      return interp_vol
  