diff --git a/btmorph2/structs/neuronMorph.py b/btmorph2/structs/neuronMorph.py
index 20626ff..2c18db5 100644
--- a/btmorph2/structs/neuronMorph.py
+++ b/btmorph2/structs/neuronMorph.py
@@ -1348,6 +1348,10 @@ def nodeXYZ(n):
                 nearerDist = nodeDist
                 nearerXYZ = nodeXYZ(node)
 
+            # Finding intersects with spheres can have slight error. To avoid inconsistencies, round distances
+            nearerDist = np.round(nearerDist, 3)
+            fartherDist = np.round(fartherDist, 3)
+
             if fartherDist > radiiSorted[0]:
                 radiiCrossedMask = np.logical_and(nearerDist < radiiSorted, radiiSorted <= fartherDist)
                 radiiCrossed = radiiSorted[radiiCrossedMask]
@@ -1424,7 +1428,11 @@ def nodeXYZ(n):
                     fartherPoint = nodeXYZ(node)
                     fartherDist = nodeDist
 
-                radiiCrossedMask = np.logical_and(nearerDist < radiiSorted, radiiSorted <= fartherDist)
+                # # Finding intersects with spheres can have slight error. To avoid inconsistencies, round distances
+                # nearerDist = np.round(nearerDist, 3)
+                # fartherDist = np.round(fartherDist, 3)
+
+                radiiCrossedMask = np.logical_and(nearerDist <= radiiSorted, radiiSorted <= fartherDist)
                 radiiCrossed = radiiSorted[radiiCrossedMask]
 
                 # line connecting the points are within one shell
@@ -1434,7 +1442,8 @@ def nodeXYZ(n):
                     intersects = getIntersectionXYZs(nearerPoint, fartherPoint, centeredAt, largestRLessNearerPoint)
 
                     # line joining the points does not intersect any sphere
-                    if len(intersects) == 0:
+                    # (len(intersects) == 1 can occur due to rounding error in finding intersects)
+                    if len(intersects) < 2:
                         shellIndex = radii.index(largestRLessNearerPoint) + 1
                         lengths[shellIndex] += np.linalg.norm(fartherPoint - nearerPoint)
                     # line joining the points intersects a sphere are two distinct points
@@ -1446,8 +1455,6 @@ def nodeXYZ(n):
                         lengths[outerShellIndex] += np.linalg.norm(fartherPoint - intersects[1])
                         lengths[innerShellIndex] += np.linalg.norm(intersects[1] - intersects[0])
 
-                    else:
-                        raise(ValueError("Impossible case! There has been a wrong assumption."))
 
                 # line connecting the points is contained in at least two shells
                 else:
@@ -1455,8 +1462,18 @@ def nodeXYZ(n):
                     for rad in radiiCrossed:
                         shellIndex = radii.index(rad)
                         intersects = getIntersectionXYZs(nearerPoint, fartherPoint, centeredAt, rad)
-                        assert len(intersects) == 1, "Impossible case! There has been a wrong assumption."
-                        intersect = np.array(intersects[0])
+                        if len(intersects) == 1:
+                            intersect = np.array(intersects[0])
+                        # len(intersects) == 2 can happen when the line joining the points crosses one sphere
+                        # two times and nearerPoint is on the sphere.
+                        elif len(intersects) == 2:
+                            intersect = np.array(intersects[1])
+                        else:
+                            # This case must theoretically not happen, but happens due to rounding errors
+                            pass
+                            # raise(ValueError("Impossible case! There has been a wrong assumption."))
+
+
                         lengths[shellIndex] += np.linalg.norm(intersect - tempNearestPoint)
                         tempNearestPoint = intersect
                     if shellIndex + 1 < len(radii):
diff --git a/tests/auxFuncs_test.py b/tests/auxFuncs_test.py
index e6cb2df..be56086 100644
--- a/tests/auxFuncs_test.py
+++ b/tests/auxFuncs_test.py
@@ -58,4 +58,5 @@ def getPointAtDistance_test():
                                                                                 atol=1e-2)
 
 if __name__ == "__main__":
-    intersects = getIntersectionXYZs(p1=[20, 30, 0], p2=[30, 40, 0], radius=50, centeredAt=[0, 0, 0])
\ No newline at end of file
+    intersects = getIntersectionXYZs(p1=[ 100.455,  159.113,  155.615], p2=[ 100.288,  158.798,  155.015],
+                                     radius=160, centeredAt=[ 254.79213333,  119.54293333,  168.7052])
\ No newline at end of file