Add an option to allow higher resolution calculation of surface areas.

share/lib/python/neuron/rxd/geometry3d/FullJoinMorph.py

@@ -178,7 +178,7 @@ def fullmorph(source, dx, soma_step=100, mesh_grid=None, relevant_pts=None):
             # no need to delete from surface voxels, it is never added to final surface voxels
         else:
             V = simplevolume(itemlist, distances, vox, grid)
-            A = surface_area(distances[0], distances[1], distances[2], distances[3], distances[4], distances[5], distances[6], distances[7], 0,dx,0,dy,0,dz)
+            A = surface_area(itemlist, vox, grid)
             final_surface_voxels[vox] = [V, A, seg]
 
     def has_vox(*vox):

share/lib/python/neuron/rxd/geometry3d/surface_a.py

@@ -1,13 +1,49 @@
 import ctypes
-from ctypes import c_double, POINTER
 import numpy
+from .. import options
 from neuron import nrn_dll_sym
+from ctypes import c_double, POINTER
+
 
 """Function to take any surface vox with information on distance of each vertex from surface 
 and output approx. surface area contained in that vox"""
 
+def surface_area(itemlist, vox, grid):
+    res = options.ics_partial_surface_resolution
+    i0, j0, k0 = vox
+    rng = range(res+1)
+    rng0 = range(res)
+    dx, dy, dz = grid['dx'], grid['dy'], grid['dz']
+    Sx, Sy, Sz = dx/res, dy/res, dz/res
+    x0 = grid['xlo'] + i0 * dx
+    y0 = grid['ylo'] + j0 * dy
+    z0 = grid['zlo'] + k0 * dz
+    dist = {}
+    for i in rng:
+        x = x0 + i * Sx
+        for j in rng:
+            y = y0 + j * Sy
+            for k in rng:
+                z =  z0 + k * Sz
+                dist[(i,j,k)] = min([f.distance(x,y,z) for f in itemlist])
+    area = 0
+    for i in rng0:
+        for j in rng0:
+            for k in rng0:
+                area += sub_surface_area(dist[i,j,k], dist[i+1,j,k],
+                                         dist[i+1,j+1,k], dist[i,j+1,k],
+                                         dist[i,j,k+1], dist[i+1,j,k+1],
+                                         dist[i+1,j+1,k+1], dist[i,j+1,k+1],
+                                         0, Sx, 0, Sy, 0, Sz)
+    return area
+
+def sub_surface_area(v0, v1, v2, v3, v4, v5, v6, v7, x0, x1, y0, y1, z0, z1):
+
+    def allornone(lst): return all(lst) or not any(lst)
+    # if all corners are inside -- return 0
+    if allornone([x <= options.ics_distance_threshold for x in [v0, v1, v2, v3, v4, v5, v6, v7]]):
+        return 0
 
-def surface_area(v0, v1, v2, v3, v4, v5, v6, v7, x0, x1, y0, y1, z0, z1):
     results = numpy.array([0.] * 48)
 
     find_triangles = nrn_dll_sym('find_triangles')

share/lib/python/neuron/rxd/options.py

@@ -15,6 +15,10 @@
 # volume fractions.
 ics_partial_volume_resolution = 2
 
+# resolution (relative to dx) of sampling points use to determine surface
+# areas of surface voxels. 
+ics_partial_surface_resolution = 1
+
 # the number of electrophysiology fixed steps per rxd step
 # WARNING: setting this to anything other than 1 is probably a very bad
 #          idea, numerically speaking, at least for now