Update 3D voxelization. (#892)

* Update 3D voxelization.
Adds;
    ics_distance_threshold option to avoid spurious voxel being included due to floating point precision when calculating a voxels position in the grid.
    ics_partial_volume_resolution option for how much to sub-sample surface voxels when calculating the volume.
Updated rxd pytests, so the tests are not performed when saving a new copy of the test data

Co-authored-by: Pramod Kumbhar <pramod.kumbhar@epfl.ch>

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

@@ -7,6 +7,7 @@
 from .GeneralizedVoxelization import voxelize
 from .simplevolume_helper import simplevolume
 from .surface_a import surface_area
+from ..options import ics_distance_threshold
 import warnings
 from neuron import h
 
@@ -35,7 +36,7 @@ def find_parent_seg(join, sdict, objects):
 
 def all_in(dist):
     for i in dist:
-        if i > 0:
+        if i > ics_distance_threshold:
             return False 
     return True
 

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

@@ -1,4 +1,5 @@
 from . import graphicsPrimitives as graphics
+from .. import options
 
 def find_voxel(x,y,z,g):
     """returns (i,j,k) of voxel containing point x,y,z"""
@@ -47,7 +48,7 @@ def verts_in(f,voxel,surf,g):
     for v in verts:
         dist = f.distance(v[0],v[1],v[2])
         distlist.append(dist)
-        if dist <= 0:
+        if dist <= options.ics_distance_threshold:
             ins+=1
     if 1 <= ins <= 7:
         surf[voxel] = distlist
@@ -184,7 +185,7 @@ def voxelize(grid, Object, corners=None, include_ga=False):
         (i0,j0,k0) = find_voxel(x0,y0,z0,grid)
         # find the contained endpoints and start the set with initial row and initial endpoints
         s = set()
-        ends = find_endpoints(Object,surface,include_ga,(j0,k0),(i0,i0),grid)
+        ends = find_endpoints(Object,surface,include_ga,(j0,k0),(i0-1,i0+1),grid)
 
     # the given starting voxel is not actually found
     possibly_missed = False

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

@@ -1,4 +1,5 @@
 from . import graphicsPrimitives as graphics
+from .. import options
 import math
 import random
 
@@ -15,7 +16,8 @@ def get_verts(voxel,g):
                 (v1_0+dx,v1_1,v1_2+dz),
                 (v1_0+dx,v1_1+dy,v1_2+dz),
                 (v1_0,v1_1+dy,v1_2+dz)]
-    return vertices 
+    return vertices
+
 
 def get_subverts(voxel,g,step):
     """return list (len=8) of point coordinates (x,y,z) that are vertices of the voxel (i,j,k)"""
@@ -33,64 +35,60 @@ def get_subverts(voxel,g,step):
                 (v1_0,v1_1+DY,v1_2+DZ)]
     return vertices 
 
-def add_res(flist, voxel, verts_in, res, g):        
-    (i, j, k) = voxel
-    dx,dy,dz = g['dx'],g['dy'],g['dz']
-    bit = (dx * dy * dz)/(res**3)
+def add_res(flist, voxel, verts_in, res, g):
+    dx, dy, dz = g['dx'],g['dy'],g['dz']
+    Sx, Sy, Sz = dx/res, dy/res, dz/res      
+    bit = Sx * Sy * Sz
 
     step = [dx/2, dy/2, dz/2]
     subverts = get_subverts(voxel, g, step)               # the 'voxel positions' for the new subvoxels
 
-    #verts_in = [0,1,2,3,4,5,6,7]
     count = 0
     #select only the subvoxels of the vertices that are in
     for i in verts_in:
         v0 = subverts[i]
-
         startpt = (v0[0]+dx/(2*res), v0[1]+dy/(2*res), v0[2]+dx/(2*res))
         for x in range(res//2):
             for y in range(res//2):
                 for z in range(res//2):
                     v = (startpt[0] + x*(dx/res), startpt[1] + y*(dy/res), startpt[2] + z*(dz/res))
-                    if  min([f.distance(v[0],v[1],v[2]) for f in flist]) <= 0:
+                    if  min([f.distance(v[0],v[1],v[2]) for f in flist]) <= options.ics_distance_threshold:
                         count += 1
 
-    vol = count*bit 
-    # add in partials if still 0
-    if vol == 0:
-        vol = bit/2
-
-    return vol
+    if count > 0:
+        return count * bit 
+    return bit / 8
+
 
 def Put(flist, voxel, v0, verts_in, res, g):
     """ add voxel key with partial volume value to dict of surface voxels"""
     # v0 is the start coordinates of voxel (verts[0] feed in)
     # res is resolution of sampling points (only works for even values of res!)
-    dx,dy,dz = g['dx'],g['dy'],g['dz']
-    bit = (dx * dy * dz)/(res**3)
+    dx, dy, dz = g['dx'], g['dy'], g['dz']
+    Sx, Sy, Sz = dx/res, dy/res, dz/res
 
     count = 0
-    startpt = (v0[0]+dx/(2*res), v0[1]+dy/(2*res), v0[2]+dx/(2*res))
-    for x in range(res):
-        for y in range(res):
-            for z in range(res):
-                v = (startpt[0] + x*(dx/res), startpt[1] + y*(dy/res), startpt[2] + z*(dz/res))
-                if min([f.distance(v[0],v[1],v[2]) for f in flist]) <= 0:
+    startpt = v0[0] + Sx/2.0, v0[1] + Sy/2.0, v0[2] + Sz/2.0
+    for i in range(res):
+        for j in range(res):
+            for k in range(res):
+                v = startpt[0] + i*Sx, startpt[1] + j*Sy, startpt[2] + k*Sz
+                if min([f.distance(v[0],v[1],v[2]) for f in flist]) <= options.ics_distance_threshold:
                     count += 1
-    bitvol = count*bit             
-    if bitvol == 0:
-        bitvol = add_res(flist, voxel, verts_in, res*2, g)
-
-    return bitvol
+    if count > 0:
+        return Sx * Sy * Sz * count
+    return add_res(flist, voxel, verts_in, 2 * res, g)
 
 def simplevolume(flist,distances,voxel,g):
     """return the number of vertices of this voxel that are contained within the surface"""
-    verts = get_verts(voxel,g)
+
+    res = options.ics_partial_volume_resolution
+    verts = (voxel,g)
     verts_in = []
     for i in range(8):
-        if distances[i] <= 0:
+        if distances[i] <= options.ics_distance_threshold:
             verts_in.append(i)
-    Vol = Put(flist, voxel, verts[0], verts_in, 2, g)
+    Vol = Put(flist, voxel, verts[0], verts_in, res, g)
     return Vol
 
 

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

@@ -7,6 +7,14 @@
 
 concentration_nodes_3d = "surface"
 
+# how far inside must a voxel be to be consider inside
+# the value is necessary to account for floating point precision
+ics_distance_threshold = -1e-12
+
+# resolution (relative to dx) of sampling points use to determine surface
+# volume fractions.
+ics_partial_volume_resolution = 2
+
 # 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

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

@@ -26,9 +26,9 @@ def _sort_secs(secs):
     for root in root_secs:
         all_sorted.wholetree(sec=root)
     secs_names = dict([(sec.hoc_internal_name(),sec) for sec in secs])
-    for sec in secs:
-        if sec.orientation():
-            raise RxDException('still need to deal with backwards sections')
+    #for sec in secs:
+    #    if sec.orientation():
+    #        raise RxDException('still need to deal with backwards sections')
     return [secs_names[sec.hoc_internal_name()] for sec in all_sorted if sec.hoc_internal_name() in secs_names]
 
 
@@ -414,6 +414,7 @@ def _indices_from_sec_x(self, sec, position):
         assert(position in (0, 1))
         # NOTE: some care is necessary in constructing normal vector... must be
         #       based on end frusta, not on vector between end points
+        dx = self.dx
         if position == 0:
             x = sec.x3d(0)
             y = sec.y3d(0)
@@ -426,10 +427,13 @@ def _indices_from_sec_x(self, sec, position):
             x = sec.x3d(n - 1)
             y = sec.y3d(n - 1)
             z = sec.z3d(n - 1)
-            # NOTE: sign of the normal is irrelevant
             nx = x - sec.x3d(n - 2)
             ny = y - sec.y3d(n - 2)
             nz = z - sec.z3d(n - 2)
+            x -= dx * nx / (nx**2 + ny**2 + nz**2)**0.5
+            y -= dx * ny / (nx**2 + ny**2 + nz**2)**0.5
+            z -= dx * nz / (nx**2 + ny**2 + nz**2)**0.5
+
         else:
             raise RxDException('should never get here')
         #dn = (nx**2 + ny**2 + nz**2)**0.5
@@ -514,7 +518,6 @@ def __init__(self, secs=None, nrn_region=None, geometry=None, dimension=None, dx
                 dx = -1
             if dx <= 0:
                 raise RxDException("dx must be a positive real number or None")
-
         self.dx = dx
         _all_regions.append(weakref.ref(self))
 

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

@@ -680,30 +680,20 @@ def _setup_matrices():
                                     hybrid_index1d_grid_ids[index1d] = grid_id
                                     index1d_sec1d[index1d] = parent_sec
                                     hybrid_vols1d[index1d] = vols1d
-                                else:
-                                    for sec1d in r._secs1d:
-                                        parent_1d_seg = sec1d.trueparentseg()
-                                        parent_1d = None if not parent_1d_seg else parent_1d_seg.sec 
-                                        if parent_1d == sec:
-                                            # it is the parent of a 1d section
-                                            index1d, indices3d, vols1d, vols3d = _get_node_indices(s, r, sec, parent_1d_seg.x , sec1d, sec1d.orientation())
-                                            hybrid_neighbors[index1d] += indices3d
-                                            hybrid_vols[index1d] += vols3d
-                                            hybrid_diams[index1d] = sec1d(h.section_orientation(sec=sec1d)).diam
-                                            hybrid_index1d_grid_ids[index1d] = grid_id
-                                            index1d_sec1d[index1d] = sec1d
-                                            hybrid_vols1d[index1d] = vols1d
-
-
-                                        elif parent_1d == parent_sec and parent_1d is not None:
-                                            # it connects to the parent of a 1d section
-                                            index1d, indices3d, vols1d, vols3d = _get_node_indices(s, r, sec, sec.orientation(), sec1d, sec1d.orientation())
-                                            hybrid_neighbors[index1d] += indices3d
-                                            hybrid_vols[index1d] += vols3d
-                                            hybrid_diams[index1d] = sec1d(h.section_orientation(sec=sec1d)).diam
-                                            hybrid_index1d_grid_ids[index1d] = grid_id
-                                            index1d_sec1d[index1d] = sec1d
-                                            hybrid_vols1d[index1d] = vols1d
+
+                                for sec1d in r._secs1d:
+                                    parent_1d_seg = sec1d.trueparentseg()
+                                    parent_1d = None if not parent_1d_seg else parent_1d_seg.sec 
+                                    if parent_1d == sec:
+                                        # it is the parent of a 1d section
+                                        index1d, indices3d, vols1d, vols3d = _get_node_indices(s, r, sec, parent_1d_seg.x , sec1d, sec1d.orientation())
+                                        hybrid_neighbors[index1d] += indices3d
+                                        hybrid_vols[index1d] += vols3d
+                                        hybrid_diams[index1d] = sec1d(h.section_orientation(sec=sec1d)).diam
+                                        hybrid_index1d_grid_ids[index1d] = grid_id
+                                        index1d_sec1d[index1d] = sec1d
+                                        hybrid_vols1d[index1d] = vols1d
+
 
 
         if len(dxs) > 1:

test/rxd/3d/test_ics_currents.py

@@ -11,7 +11,7 @@ def ics_example(neuron_instance):
        currents.
     """
 
-    h, rxd, data = neuron_instance
+    h, rxd, data, save_path = neuron_instance
     rxd.set_solve_type(dimension=3)
     # create cell1 where `x` will be created and leak out
     cell1 = h.Section(name='cell1')
@@ -44,23 +44,25 @@ def ics_example(neuron_instance):
 def test_ics_currents(ics_example):
     """Test ics_example with fixed step methods"""
 
-    (h, rxd, data), model = ics_example
+    (h, rxd, data, save_path), model = ics_example
     h.finitialize(-65)
     h.continuerun(100)
 
-    max_err = compare_data(data)
-    assert max_err < tol
+    if not save_path: 
+        max_err = compare_data(data)
+        assert max_err < tol
 
 
 def test_ics_currents_cvode(ics_example):
     """Test ics_example with variable step methods"""
 
-    (h, rxd, data), model = ics_example
+    (h, rxd, data, save_path), model = ics_example
 
     h.CVode().active(True)
 
     h.finitialize(-65)
     h.continuerun(100)
 
-    max_err = compare_data(data)
-    assert max_err < tol
+    if not save_path: 
+        max_err = compare_data(data)
+        assert max_err < tol

test/rxd/3d/test_include_flux3d.py

@@ -11,7 +11,7 @@ def ics_include_flux(neuron_instance):
     pointer. All three are tested here for 3d intracellular rxd.
     """
 
-    h, rxd, data = neuron_instance
+    h, rxd, data, save_path = neuron_instance
     sec = h.Section(name="sec")
     sec.L = 1 
     sec.nseg = 11 
@@ -41,21 +41,23 @@ def test_include_flux3d(ics_include_flux):
     """Test ics_include_flux with fixed step methods"""
 
     neuron_instance, model = ics_include_flux
-    h, rxd, data = neuron_instance
+    h, rxd, data, save_path = neuron_instance
     h.finitialize(-70)
     h.continuerun(10)
-    max_err = compare_data(data)
-    assert max_err < tol
+    if not save_path:
+        max_err = compare_data(data)
+        assert max_err < tol
 
 
 def test_include_flux3d_cvode(ics_include_flux):
     """Test ics_include_flux with variable step methods"""
 
     neuron_instance, model = ics_include_flux
-    h, rxd, data = neuron_instance
+    h, rxd, data, save_path = neuron_instance
     h.CVode().active(True)
     h.finitialize(-70)
     h.continuerun(10)
 
-    max_err = compare_data(data)
-    assert max_err < tol
+    if not save_path:
+        max_err = compare_data(data)
+        assert max_err < tol