fractal dimension computing implemented and tested, but is very slow;…

… fix branching angle bug to deal with special cases.

REQUIRE

@@ -2,3 +2,4 @@ julia 0.5
 LightGraphs
 BigArrays
 Libz
+LsqFit

src/BranchNets.jl

@@ -3,6 +3,7 @@ include("Branches.jl")
 using .Branches 
 using ..RealNeuralNetworks.NodeNets
 using ..RealNeuralNetworks.SWCs
+using LsqFit
 
 const ONE_UINT32 = UInt32(1)
 const EXPANSION = (ONE_UINT32, ONE_UINT32, ONE_UINT32)
@@ -467,6 +468,11 @@ function get_branching_angle( self::BranchNet, branchIndex::Integer; nodeDistanc
             break 
         end 
     end
+    if parentNode == branchingNode
+        warn("parent node is the same with branching node: $(branchingNode)")
+        return 0.0
+    end 
+
     childNode = branch[1]
     for index in length(branch):1
         node = branch[index]
@@ -475,14 +481,24 @@ function get_branching_angle( self::BranchNet, branchIndex::Integer; nodeDistanc
             break 
         end 
     end 
+    if childNode == branchingNode 
+        warn("child node is the same with branching node: $(branchingNode)")
+        return 0.0 
+    end 
+
     # compute the angle among the three nodes using the definition of dot product
     # get the two vectors
     v1 = map(-, parentNode[1:3], branchingNode[1:3] )
     v2 = map(-, branchingNode[1:3], childNode[1:3] )
     # normalize the vector
     nv1 = normalize([v1...]) 
     nv2 = normalize([v2...])
-    return acos( dot(nv1, nv2) )
+    #@show nv1, nv2
+    #@show childNode, branchingNode, parentNode
+    dotProduct = dot(nv1, nv2)
+    # tolerate some numerical varition. the dot product could go greater than 1.
+    @assert dotProduct < 1.001 "impossible dotProduct: $(dotProduct)"
+    return acos( min(1.0, dotProduct) )
 end 
 
 """
@@ -591,7 +607,7 @@ function get_fractal_dimension( self::BranchNet )
     # iterate all the nodes inside gyration radius as the center of scanning disks
     averageMassList = zeros( length(radiusList) )
     for (centerIndex, center) in enumerate(diskCenterList)
-        for (radiusIndex, radius) in enumerate()
+        for (radiusIndex, radius) in enumerate(radiusList)
             for node in nodeList 
                 if Branches.get_nodes_distance( center,  node) < radius
                     averageMassList[radiusIndex] += 1.0 / length(diskCenterList)
@@ -601,7 +617,11 @@ function get_fractal_dimension( self::BranchNet )
     end 
 
     # fit the curve and get slop as fractal dimension
-    error("imcomplete implementation...")
+    model(x,p) = p[1]*x + p[2]
+    p0 = [1.0, 0]
+    fit = curve_fit(model, log(radiusList), log(averageMassList), p0)
+    fractalDimension = fit.param[1]
+    return fractalDimension, radiusList, averageMassList 
 end 
 
 ############################### Base functions ###################################

src/Branches.jl

@@ -17,7 +17,7 @@ type Branch
     boundingBox ::BoundingBox
 end 
 
-function Branch(nodeList::Vector, class=CLASS)
+function Branch(nodeList::Vector; class=CLASS)
     Branch(nodeList, class, BoundingBox(nodeList))
 end 
 
@@ -26,7 +26,7 @@ end
     get_nodes_distance(self::Node, other::Node)
 compute the euclidean distance between two nodes 
 """
-function get_nodes_distance(self::Node, other::Node)
+function get_nodes_distance(self::Union{Vector,Tuple}, other::Union{Vector,Tuple})
     norm( [map((x,y)->x-y, self[1:3], other[1:3]) ...])
 end 
 function get_node_list(self::Branch) self.nodeList end 
@@ -175,7 +175,7 @@ function remove_node(self::Branch, removeNodeIndex::Integer)
             push!(newNodeList, node)
         end 
     end 
-    Branch(newNodeList, get_class(self))
+    Branch(newNodeList; class = get_class(self))
 end 
 
 function remove_redundent_nodes!(self::Branch)

test/BranchNets.jl

@@ -4,31 +4,9 @@ using RealNeuralNetworks.BranchNets
 using RealNeuralNetworks.NodeNets
 using RealNeuralNetworks.SWCs
 
-@testset "test BranchNet IO" begin 
-    branchNet = BranchNets.load_swc( joinpath(dirname(@__FILE__), "../assert/example.swc" ))
-    BranchNets.save(branchNet, "/tmp/branchNet.swc")
-    rm("/tmp/branchNet.swc")
-end 
-
 @testset "test BranchNets" begin
-    #println("create fake cylinder segmentation...")
-    #@time seg = FakeSegmentations.broken_cylinder()
-    #println("skeletonization to build a BranchNet ...")
-    #@time branchNet = BranchNet(seg)
-    println("create fake ring segmentation ...")
-    seg = FakeSegmentations.broken_ring()
-    branchNet = BranchNet(seg)
-    println("transform to SWC structure ...")
-    @time swc = SWCs.SWC( branchNet )
-    tempFile = tempname() * ".swc"
-    SWCs.save(swc, tempFile)
-    rm(tempFile)
-
     println("load swc of a real neuron...")
-    @time swc = SWCs.load_gzip_swc("../assert/76869.swc.gz")
-    println("save as gzip compressed file ...")
-    @time SWCs.save_gzip_swc(swc, "/tmp/76869.swc.gz")
-    rm("/tmp/76869.swc.gz")
+    @time swc = SWCs.load_swc_bin("../assert/76869.swc.bin")
 
     branchNet = BranchNet( swc )
     println("get node list ...")
@@ -38,8 +16,16 @@ end
     println("get branch order list...")
     @time branchOrderList = BranchNets.get_branch_order_list( branchNet )
 
-    println("remove subtree in soma...")
-    @time newBranchNet = BranchNets.remove_subtree_in_soma(branchNet)
+    println("clean up the neuron ...")
+    branchNet = BranchNets.remove_subtree_in_soma(branchNet)
+    branchNet = BranchNets.remove_hair(branchNet)
+    branchNet = BranchNets.remove_subtree_in_soma(branchNet)
+    branchNet = BranchNets.remove_terminal_blobs(branchNet)
+    branchNet = BranchNets.remove_redundent_nodes(branchNet)
+
+    println("get fractal dimension ...")
+    @time fractalDimension, _,_ = BranchNets.get_fractal_dimension( branchNet )
+    @show fractalDimension 
 
     println("get typical radius ...")
     @show BranchNets.get_typical_radius( branchNet )
@@ -78,3 +64,27 @@ end
     @test !isempty(tree1)
     @test !isempty(tree2)
 end 
+
+@testset "test BranchNet IO" begin 
+    branchNet = BranchNets.load_swc( joinpath(dirname(@__FILE__), "../assert/example.swc" ))
+    BranchNets.save(branchNet, "/tmp/branchNet.swc")
+    rm("/tmp/branchNet.swc")
+end 
+
+@testset "test fake segmentation skeletonization" begin 
+    println("create fake cylinder segmentation...")
+    @time seg = FakeSegmentations.broken_cylinder()
+    println("skeletonization to build a BranchNet ...")
+    @time branchNet = BranchNet(seg)
+
+    println("create fake ring segmentation ...")
+    seg = FakeSegmentations.broken_ring()
+    branchNet = BranchNet(seg)
+    println("transform to SWC structure ...")
+    @time swc = SWCs.SWC( branchNet )
+    tempFile = tempname() * ".swc"
+    SWCs.save(swc, tempFile)
+    rm(tempFile)
+end 
+
+