per partes node aggregation function

lisa/skeleton_analyser.py

@@ -327,14 +327,67 @@ def __skeleton_nodes_aggregation(self, data3d_skel):
 
         aggregate near nodes
         """
+
+        method = 'auto'
         if self.aggregate_near_nodes_distance >= 0:
-#          TODO doplnit zzávislost na voxelsize 
+            print 'generate structure'
             structure = generate_binary_elipsoid(
                     self.aggregate_near_nodes_distance / np.asarray(self.voxelsize_mm))
-            nd_dil = scipy.ndimage.binary_dilation(data3d_skel==2, structure)
+            print 'perform dilation ', data3d_skel.shape
+            # import ipdb; ipdb.set_trace() #  noqa BREAKPOINT
+
+            # TODO select best method
+            # old simple method
+            # nd_dil = scipy.ndimage.binary_dilation(data3d_skel==2, structure)
+
+            # per partes method even slower
+            nd_dil = self.__skeleton_nodes_aggregation_per_each_node(data3d_skel==2, structure)
+            print "dilatation completed"
             data3d_skel[nd_dil] = 2
+            print 'set to 2'
         return data3d_skel
 
+    def __skeleton_nodes_aggregation_per_each_node(self, data3d_skel2, structure):
+        node_list = np.nonzero(data3d_skel2)
+        nlz = zip(node_list[0], node_list[1], node_list[2])
+
+        print 'zip finished'
+
+        for node_xyz in nlz:
+            self.__node_dilatation(data3d_skel2, node_xyz, structure)
+
+        return data3d_skel2
+
+
+    def __node_dilatation(self, data3d_skel2, node_xyz, structure):
+        """
+        this function is called for each node
+        """
+        border = structure.shape
+
+        print 'border ', border
+        print structure.shape
+        xlim = [max(0, node_xyz[0] - border[0]), min(data3d_skel2.shape[0], node_xyz[0] + border[0])]
+        ylim = [max(0, node_xyz[1] - border[1]), min(data3d_skel2.shape[1], node_xyz[1] + border[1])]
+        zlim = [max(0, node_xyz[2] - border[2]), min(data3d_skel2.shape[2], node_xyz[2] + border[2])]
+
+
+        # dilation on small box
+        nd_dil = scipy.ndimage.binary_dilation(
+            data3d_skel2[xlim[0]:xlim[1],
+                         ylim[0]:ylim[1],
+                         zlim[0]:zlim[1]]==2, structure)
+
+        # nd_dil = nd_dil * 2
+
+        data3d_skel2[xlim[0]:xlim[1],
+                        ylim[0]:ylim[1],
+                        zlim[0]:zlim[1]]=nd_dil
+
+
+
+
+
 
     def __label_edge_by_its_terminal(self, labeled_terminals):
         import functools

tests/skeleton_test.py

@@ -19,20 +19,49 @@
 class TemplateTest(unittest.TestCase):
 
     @attr('actual')
+    @attr('slow')
+    def test_nodes_aggregation_big_data(self):
+
+        data = np.zeros([1000, 1000, 100], dtype=np.int8)
+        voxelsize_mm = [14, 10, 6]
+
+        # snake
+        # data[15:17, 13, 13] = 1
+        data[18, 3:17, 12] = 1
+        data[18, 17, 13:17] = 1
+        data[18, 9, 4:12] = 1
+        data[18, 11, 12:19] = 1
+        # data[18, 18, 15:17] = 1
+
+        # T-junction on the left
+        data[18, 4:16, 3] = 1
+        import sed3
+        # ed = sed3.sed3(data)
+        # ed.show()
+
+        skel = data
+
+        skan = sk.SkeletonAnalyser(
+                copy.copy(skel), 
+                volume_data=data,
+                voxelsize_mm=voxelsize_mm, 
+                cut_wrong_skeleton=False, 
+                aggregate_near_nodes_distance=20)
+        vessel_tree = skan.skeleton_analysis()
+
+        # ed = sed3.sed3(skan.sklabel, contour=data)
+        # ed.show()
+        # import ipdb; ipdb.set_trace() #  noqa BREAKPOINT
+
+        # number of terminals + inner nodes
+        self.assertEqual(np.min(skan.sklabel), -7)
+        # number of edges
+        self.assertEqual(np.max(skan.sklabel), 6)
+
     def test_nodes_aggregation(self):
 
         data = np.zeros([20, 20, 20], dtype=np.int8)
-
-        # data[2:13, 4, 4] = 1
-        # data[6, 2:13, 6] = 1
-        # data[8, 8, 2:13] = 1
-        #
-        # data[6, 8, 6:17] = 1
-        #
-        # # diagonala
-        # data[11, 11, 11] = 1
-        # data[12, 12, 12] = 1
-        # data[13, 13, 13] = 1
+        voxelsize_mm = [14, 10, 6]
 
         # snake
         # data[15:17, 13, 13] = 1
@@ -50,8 +79,12 @@ def test_nodes_aggregation(self):
 
         skel = data
 
-        skan = sk.SkeletonAnalyser(copy.copy(skel), volume_data=data,
-                                   voxelsize_mm=[14, 10, 6], cut_wrong_skeleton=False, aggregate_near_nodes_distance=20)
+        skan = sk.SkeletonAnalyser(
+                copy.copy(skel), 
+                volume_data=data,
+                voxelsize_mm=voxelsize_mm, 
+                cut_wrong_skeleton=False, 
+                aggregate_near_nodes_distance=20)
         vessel_tree = skan.skeleton_analysis()
 
         # ed = sed3.sed3(skan.sklabel, contour=data)