compute arbor density map

src/BoundingBoxes.jl

@@ -9,7 +9,14 @@ mutable struct BoundingBox
     maxCorner   ::NTuple{3,Float32}
 end 
 
-function BoundingBox(minCorner::Vector, maxCorner::Vector)
+function BoundingBox()
+    minCorner = (Inf32, Inf32, Inf32)
+    maxCorner = (ZERO_FLOAT32, ZERO_FLOAT32, ZERO_FLOAT32)
+    BoundingBox( minCorner, maxCorner )
+end 
+
+function BoundingBox(   minCorner::Vector = [Inf32, Inf32, Inf32],
+                        maxCorner::Vector = [ZERO_FLOAT32, ZERO_FLOAT32, ZERO_FLOAT32])
     @assert length(minCorner) == 3
     @assert length(maxCorner) == 3
     BoundingBox((minCorner...), (maxCorner...))
@@ -28,14 +35,21 @@ function BoundingBox(nodeList::Vector)
     BoundingBox(minCorner, maxCorner)
 end 
 
+function Base.size(self::BoundingBox)
+    map((x,y)->ceil(Int, x-y), self.maxCorner, self.minCorner)
+end 
+
 function Base.isequal(self::BoundingBox, other::BoundingBox)
     self.minCorner==other.minCorner && self.maxCorner==other.maxCorner 
 end
 function Base.:(==)(self::BoundingBox, other::BoundingBox)
     isequal(self, other) 
 end 
 
-
+"""
+    Base.union(self::BoundingBox, other::BoundingBox)
+get the bounding box containing these two bounding boxes 
+"""
 function Base.union(self::BoundingBox, other::BoundingBox)
     minCorner = map(min, self.minCorner, other.minCorner)
     maxCorner = map(max, self.maxCorner, other.maxCorner)
@@ -47,6 +61,16 @@ function isinside(self::BoundingBox, point::Union{Tuple, Vector})
     all( map((x,y)->x<y, self.minCorner, point[1:3]) )
 end 
 
+"""
+    to_voxel_space(self::BoundingBox, voxelSize::Union{Tuple, Vector})
+
+"""
+function to_voxel_space(self::BoundingBox, voxelSize::Union{Tuple, Vector})
+    minCorner = map(/, self.minCorner, voxelSize)
+    maxCorner = map(/, self.maxCorner, voxelSize)
+    BoundingBox(minCorner, maxCorner)
+end 
+
 """
     distance_from(self::BoundingBox, point::Union{Tuple, Vector})
 

src/BranchNets.jl

@@ -1,10 +1,12 @@
 module BranchNets
 include("Branches.jl")
-using .Branches 
+using .Branches
+using .Branches.BoundingBoxes 
 using ..RealNeuralNetworks.NodeNets
 using ..RealNeuralNetworks.SWCs
 using LsqFit
 #using JLD2
+using ImageFiltering
 
 const ONE_UINT32 = UInt32(1)
 const ZERO_FLOAT32 = Float32(0)
@@ -628,6 +630,47 @@ function get_fractal_dimension( self::BranchNet )
     return fractalDimension, radiusList, averageMassList 
 end 
 
+function Base.BitArray(self::BranchNet, voxelSize::Union{Tuple, Vector})
+    nodeList = get_node_list(self)
+    voxelList = Vector{NTuple{3, Float32}}()
+    for node in nodeList 
+        voxelCoordinate = map(/, node[1:3], voxelSize)
+        push!(voxelList, (voxelCoordinate...))
+    end 
+    boundingBox = BoundingBox( voxelList )
+    @show boundingBox 
+    sz = size(boundingBox)
+    # initialize the map 
+    bitMask = falses(map(x->x+1, sz))
+    @inbounds for voxel in voxelList
+        # add a small number to make sure that there is no 0
+        idx = [voxel...] .- boundingBox.minCorner + 0.0000001
+        idx = map(x->ceil(Int, x), idx)
+        bitMask[idx...] = true 
+    end 
+    bitMask
+end 
+
+"""
+    get_arbor_density_map(self::BranchNet, 
+                        voxelSize::Union{Tuple, Vector},
+                        gaussianFilterStd ::AbstractFloat)
+compute the arbor density map
+
+Return:
+    * densityMap::Array{Float64, 3} 
+"""
+function get_arbor_density_map(self::BranchNet, voxelSize::Union{Tuple, Vector},
+                                            gaussianFilterStd::AbstractFloat)
+    bitMask = BitArray(self, voxelSize)
+    densityMap = Array{Float64,3}(bitMask)
+    # gaussion convolution
+    kernel = fill(gaussianFilterStd, 3) |> Kernel.gaussian
+    @time densityMap = imfilter(densityMap, kernel)
+    densityMap
+end 
+
+
 ############################### Base functions ###################################
 function Base.getindex(self::BranchNet, index::Integer)
     get_branch_list(self)[index]
@@ -691,8 +734,8 @@ function Base.merge(self::BranchNet, other::BranchNet,
                                     1:size(self.connectivityMatrix,2)] = 
                                                         self.connectivityMatrix 
                 # do not include the connection of root in net2
-                mergedConnectivityMatrix[num_branches1+1+1 : end, 
-                                         num_branches1+1+1 : end] = 
+                mergedConnectivityMatrix[num_branches1+1 : end, 
+                                         num_branches1+1 : end] = 
                                                 other.connectivityMatrix[2:end, 2:end]
                 # reestablish the connection of root2
                 childrenBranchIndexList2 = get_children_branch_index_list(other, 1)
@@ -704,7 +747,7 @@ function Base.merge(self::BranchNet, other::BranchNet,
             end 
         end 
     else 
-        # need to break the nearest branch and rebuild connectivity matrix
+        #println("need to break the nearest branch and rebuild connectivity matrix")
         total_num_branches = num_branches1 + 1 + num_branches2 
         mergedBranchList = branchList1 
         mergedConnectivityMatrix = spzeros(Bool, total_num_branches, total_num_branches)

src/Branches.jl

@@ -1,7 +1,6 @@
 module Branches
 
 include("BoundingBoxes.jl")
-# include(joinpath(dirname(@__FILE__), "BoundingBoxes.jl"))
 using .BoundingBoxes 
 
 const Node = NTuple{4,Float32}

test/BranchNets.jl

@@ -10,9 +10,10 @@ using RealNeuralNetworks.SWCs
     branchNet = BranchNet( swc )
 
     BranchNets.save(branchNet, "/tmp/branchNet.swc")
-    branchNet2 = BranchNets.resample(branchNet, Float32(10))
+    branchNet2 = BranchNets.resample(branchNet, Float32(40))
     BranchNets.save_swc(branchNet2, "/tmp/branchNet2.swc")
-    #rm("/tmp/branchNet.swc")
+    rm("/tmp/branchNet.swc")
+    rm("/tmp/branchNet2.swc")
 end 
 
 @testset "test BranchNets" begin