From 4274a149c572d269786113aea40997d08794b7f0 Mon Sep 17 00:00:00 2001 From: jingpengw Date: Tue, 30 Apr 2019 10:42:08 -0400 Subject: [PATCH] change to GeometryTypes for better meaning --- REQUIRE | 1 + src/Neurons.jl | 112 ++++++++++++------------------------ src/NodeNets.jl | 43 +++++++------- src/PointArrays.jl | 8 ++- src/Segments.jl | 48 +++++++--------- src/Synapses.jl | 11 ++-- src/Utils/BoundingBoxes.jl | 38 ++++++------ test/Neurons.jl | 5 +- test/NodeNets.jl | 8 ++- test/PointArrays.jl | 4 +- test/Utils/BoundingBoxes.jl | 5 +- 11 files changed, 121 insertions(+), 162 deletions(-) diff --git a/REQUIRE b/REQUIRE index ab06eda3..ff7d8b8e 100644 --- a/REQUIRE +++ b/REQUIRE @@ -12,3 +12,4 @@ Query JSON NearestNeighbors ProgressMeter +GeometryTypes diff --git a/src/Neurons.jl b/src/Neurons.jl index 837407ac..b54be981 100644 --- a/src/Neurons.jl +++ b/src/Neurons.jl @@ -16,9 +16,10 @@ using LinearAlgebra using Serialization using Statistics using NearestNeighbors +import GeometryTypes: Vec, Vec4, Vec4f0, Vec3, Vec3f0 -const EXPANSION = (one(UInt32), one(UInt32), one(UInt32)) -const VOXEL_SIZE = (500, 500, 500) +const EXPANSION = Vec3(one(UInt32)) +const VOXEL_SIZE = Vec3(500) export Neuron @@ -110,7 +111,7 @@ function Neuron!(seedNodeId::Integer, nodeNet::NodeNet{T}, collectedFlagVec::Bit while !isempty(segmentSeedList) seedNodeId, segmentParentId = pop!(segmentSeedList) # grow a segment from seed - nodeListInSegment = Vector{NTuple{4, T}}() + nodeListInSegment = Vec4f0[] seedNodeIdList = [seedNodeId] while true # construct this segment @@ -169,8 +170,8 @@ function Neuron!(seedNodeId::Integer, nodeNet::NodeNet{T}, collectedFlagVec::Bit Neuron{T}(segmentList, connectivityMatrix) end -function Neuron( seg::Array{ST,3}; obj_id::ST = convert(ST,1), - expansion::NTuple{3,UInt32}=EXPANSION ) where ST +function Neuron(seg::Array{ST,3}; obj_id::ST = convert(ST,1), + expansion::Vec3{TE}=EXPANSION ) where {ST, TE} nodeNet = NodeNet( seg; obj_id = obj_id, expansion = expansion ) Neuron( nodeNet ) end @@ -360,7 +361,7 @@ end get the node list. the first one is root node. """ function get_node_list( self::Neuron{T} ) where T - nodeList = Vector{NTuple{4, T}}(undef, get_node_num(self)) + nodeList = Vector{Vec4f0}(undef, get_node_num(self)) i = 0 for segment in get_segment_list(self) @@ -583,7 +584,7 @@ Return: get_furthest_terminal_node_pair_direction(terminalNodeList) end -function get_furthest_terminal_node_pair_direction(terminalNodeList::Vector{NTuple{4,T}}) where T +function get_furthest_terminal_node_pair_direction(terminalNodeList::Vector{Vec4{T}}) where T # find the farthest node pair vec = zeros(T, 2) maxNodeDistance = zero(T) @@ -612,7 +613,7 @@ function get_mass_center( self::Neuron ) x = sum(map(n->n[1], nodeList)) / length(nodeList) y = sum(map(n->n[2], nodeList)) / length(nodeList) z = sum(map(n->n[3], nodeList)) / length(nodeList) - (x,y,z) + Vec3f0(x,y,z) end """ @@ -624,7 +625,7 @@ function get_typical_radius( self::Neuron ) nodeList = get_node_list( self ) typicalRadius = 0.0 for node in nodeList - typicalRadius += norm([massCenter...] .- [node[1:3]...]) / length(nodeList) + typicalRadius += norm(massCenter .- node[1:3]) / length(nodeList) end typicalRadius end @@ -636,7 +637,7 @@ asymmetry was measured by the euclidean distance between root node and mass cent function get_asymmetry( self::Neuron ) massCenter = get_mass_center( self ) root = get_root_node( self ) - norm( [massCenter...] .- [root[1:3]...] ) + norm( massCenter .- root[1:3] ) end """ @@ -770,7 +771,7 @@ function get_sholl_number(self::Neuron, shollRadius::AbstractFloat) get_sholl_number( root, nodeList, edgeList, shollRadius ) end -function get_sholl_number(root::NTuple{4, Float32}, nodeList::Vector{NTuple{4,Float32}}, +function get_sholl_number(root::Vec4f0, nodeList::Vector{Vec4f0}, edgeList::Vector{NTuple{2,Int}}, shollRadius::AbstractFloat) shollNum = 0 for edge in edgeList @@ -842,7 +843,7 @@ function get_surface_area(self::Neuron) parentSegmentId = get_parent_segment_id(self, segmentId) if parentSegmentId > 0 parentNode = segmentList[parentSegmentId][end] - h = Segments.euclidean_distance(parentNode[1:3], segment[1][1:3]) + h = norm(parentNode[1:3] .- segment[1][1:3]) # average diameter r1 = parentNode[4] r2 = segment[1][4] @@ -869,7 +870,7 @@ function get_volume(self::Neuron) node = segment[1] r1 = node[4] r2 = parentNode[4] - h = Segments.euclidean_distance(node[1:3], parentNode[1:3]) + h = norm(node[1:3] .- parentNode[1:3]) ret += pi * h * (r1*r1 + r1*r2 + r2*r2) / Float32(3) end end @@ -892,7 +893,7 @@ function get_fractal_dimension( self::Neuron ) gyrationRadius = get_gyration_radius( self; nodeList=nodeList, massCenter=massCenter ) # find the nodes inside the gyration radius - diskCenterList = Vector{NTuple{4,Float32}}() + diskCenterList = Vector{Vec4f0}() for node in nodeList if Segments.get_nodes_distance( massCenter, node ) < gyrationRadius push!(diskCenterList, node) @@ -931,15 +932,15 @@ end """ - get_mask(self::Neuron, voxelSize::Union{Tuple, Vector}) + get_mask(self::Neuron, voxelSize::Vec3f0) compute binary representation of the neuron skeleton """ -function get_mask(self::Neuron, voxelSize::Union{Tuple, Vector}) +function get_mask(self::Neuron{T}, voxelSize::Vec3{T}) where T nodeList = get_node_list(self) - voxelSet = Set{NTuple{3, Int}}() + voxelSet = Set{Vec3{Int}}() for node in nodeList voxelCoordinate = map((x,y)->round(Int, x/y), node[1:3], voxelSize) - push!(voxelSet, (voxelCoordinate...,)) + push!(voxelSet, voxelCoordinate) end boundingBox = Segments.BoundingBox( voxelSet ) range = BoundingBoxes.get_unit_range(boundingBox) @@ -987,15 +988,15 @@ end """ get_arbor_density_map(self::Neuron, - voxelSize::Union{Tuple, Vector}, + voxelSize::Vec3f0, gaussianFilterStd ::AbstractFloat) compute the arbor density map Return: * densityMap::OffsetArray """ -function get_arbor_density_map(self::Neuron, voxelSize::Union{Tuple, Vector}, - gaussianFilterStd::AbstractFloat) +function get_arbor_density_map(self::Neuron, voxelSize::Vec3f0, + gaussianFilterStd::AbstractFloat) densityMap = get_mask(self, voxelSize) # gaussion convolution kernel = Kernel.gaussian((3,3,3)) @@ -1014,7 +1015,7 @@ end translate the soma to coordinate (0,0,0) """ function translate_soma_to_coordinate_origin(self::Neuron, densityMap::OffsetArray, - voxelSize::Tuple) + voxelSize::Vec3) # assume that the first node is the soma node with largest radius somaCorrdinate = get_root_node(self)[1:3] newRange = map((x,y,s)->(x .- round(Int,y/s)), @@ -1043,7 +1044,7 @@ function get_arbor_density_map_distance(self::OffsetArray{T,N,Array{T,N}}, end function _find_correlation_peak_along_axis(self::Array{T,3}, other::Array{T,3}, - axis::Int) where {T} + axis::Int) where {T} crossCorrelation = _cross_correlation_along_axis(self, other, axis) maximum(crossCorrelation) end @@ -1823,8 +1824,6 @@ end remove_axons(self::Neuron) = remove_segments_with_class(self, Segments.AXON_CLASS) remove_dendrites(self::Neuron) = remove_segments_with_class(self, Segments.DENDRITE_CLASS) - - function downsample_nodes(self::Neuron{T}; nodeNumStep::Int=24) where T @assert nodeNumStep > 1 segmentList = get_segment_list(self) @@ -1832,10 +1831,10 @@ function downsample_nodes(self::Neuron{T}; nodeNumStep::Int=24) where T sizehint!(newSegmentList, length(segmentList) ) for segment in segmentList nodeList = Segments.get_node_list(segment) - newNodeList = Vector{Segments.Node}() + newNodeList = Vector{Vec4{T}}() for i in 1:nodeNumStep:length(nodeList) center = Segments.get_center(segment, - i: min(i+nodeNumStep-1,length(nodeList))) + i: min(i+nodeNumStep-1,length(nodeList))) push!(newNodeList, center) end newSegment = Segment(newNodeList, class=Segments.get_class(segment)) @@ -1846,13 +1845,13 @@ end """ - smooth(nodeList::Vector{NTuple{4,T}}, k::Integer) where T + smooth(nodeList::Vector{Vec4{T}, k::Integer) where T """ -function smooth(nodeList::Vector{NTuple{4,T}}) where T +function smooth(nodeList::Vector{Vec4{T}}) where T if length(nodeList) < 3 return nodeList end - newNodeList = Vector{NTuple{4,T}}() + newNodeList = Vec4f0[] push!(newNodeList, nodeList[1]) for i in 2:length(nodeList)-1 newNode = map((x,y,z)-> (x+y+z) / T(3), nodeList[i-1], nodeList[i], nodeList[i+1]) @@ -1938,7 +1937,7 @@ function resample(self::Neuron{T}, resampleDistance::T) where T neuron end -function resample(nodeList::Vector{NTuple{4,T}}, resampleDistance::T) where T +function resample(nodeList::Vector{Vec4{T}}, resampleDistance::T) where T @assert !isempty(nodeList) # always keep the first node ret = [nodeList[1]] @@ -2001,66 +2000,27 @@ function build_tree_with_position(self::Neuron{T}; leafSize::Integer=1) where T return kdTree, nodePosition end -function find_closest_node(self::Neuron{T}, seedNode::NTuple{3,T}) where T +function find_closest_node(self::Neuron{T}, seedNode::Union{Vec{N,T}, Vector{T}, NTuple{N,T}}) where {N,T} kdTree, nodePosition = build_tree_with_position(self) find_closest_node(kdTree, nodePosition, seedNode) end function find_closest_node(kdTree::KDTree, nodePosition::Matrix{Int}, - seedNode::Union{Tuple, Vector}) - idxs, dists = knn(kdTree, [seedNode[1:3]...], 1, false) + seedNode::Union{Vec{N,T}, Vector{T}, NTuple{N,T}}) where {N,T} + idxs, dists = knn(kdTree, seedNode[1:3], 1, false) segmentId, nodeIdInSegment = nodePosition[:, idxs[1]] return segmentId, nodeIdInSegment, dists[1] end """ - find_closest_node_V1(self::Neuron{T}, seedNode::NTuple{3,T}) where T - -Return: - closestSegmentId::Int, the closest segment Id - closestNodeIdInSegment::Int, the closest node id in the segment - distance::T, the closest distance -""" -function find_closest_node_V1(self::Neuron{T}, seedNode::NTuple{3,T}) where T - segmentList = get_segment_list(self) - - # initialization - closestSegmentId = 0 - closestNodeIdInSegment = 0 - distance = typemax(T) - - @assert !isempty(segmentList) - - for (segmentId, segment) in enumerate(segmentList) - # distance from bounding box give the maximum bound of closest distance - # this was used to filter out far away segments quickly - # no need to compute all the distances between nodes - # this is the lower bound of the distance - bbox_distance = Segments.get_bounding_box_distance(segment, seedNode) - if bbox_distance < distance - d, nodeIdInSegment = Segments.distance_from(segment, seedNode ) - if d < distance - distance = d - closestSegmentId = segmentId - closestNodeIdInSegment = nodeIdInSegment - end - end - end - @assert closestNodeIdInSegment > 0 - @assert closestNodeIdInSegment <= length(segmentList[closestSegmentId]) - closestSegmentId, closestNodeIdInSegment, distance -end - - -""" - find_merging_terminal_node_id(self::Neuron, seedNode2::NTuple{4, Float32}) + find_merging_terminal_node_id(self::Neuron, seedNode2::Vec4{Float32}) # Note that the weighted version was commented out since it seems not working well #the distance was weighted according to angle θ: d' = d * tan(θ) #this function has a nice property that shrinking the distance within 45 degrees #and enlarge it with angle greater than 45 degrees """ -function find_merging_terminal_node_id(self::Neuron{T}, seedNode2::NTuple{4, T}) where T +function find_merging_terminal_node_id(self::Neuron{T}, seedNode2::Vec4{T}) where T # initialization closestTerminalSegmentId1 = 0 terminalNodeIdInSegment1 = 0 @@ -2202,7 +2162,7 @@ function find_closest_node(neuron::Neuron, nodeNet::NodeNet, collectedFlagVec::B return closestNodeId, mergingSegmentId, mergingNodeIdInSegment end -function add_offset(self::Neuron{T}, offset::Union{Vector, Tuple}) where T +function add_offset(self::Neuron{T}, offset::Vec3f0) where T segmentList = Vector{Segment{T}}() for segment in self.segmentList newSegment = Segments.add_offset(segment, offset) diff --git a/src/NodeNets.jl b/src/NodeNets.jl index 76cc71a5..9b282c45 100644 --- a/src/NodeNets.jl +++ b/src/NodeNets.jl @@ -3,6 +3,7 @@ include("DBFs.jl"); using .DBFs; include("PointArrays.jl"); using .PointArrays; using ..RealNeuralNetworks.SWCs import LinearAlgebra: norm +import GeometryTypes: Vec2, Vec4, Vec4f0, Vec3, Vec3f0 export NodeNet @@ -11,10 +12,10 @@ import BigArrays using SparseArrays import DataStructures: IntSet -const OFFSET = (zero(UInt32), zero(UInt32), zero(UInt32)) +const OFFSET = Vec3(zero(UInt32)) # rescale the skeleton -const EXPANSION = (one(UInt32), one(UInt32), one(UInt32)) +const EXPANSION = Vec3(one(UInt32)) # control the removing points around path based on DBF const REMOVE_PATH_SCALE = 3 @@ -23,7 +24,7 @@ const REMOVE_PATH_CONST = 4 mutable struct NodeNet{T} # x,y,z,r - nodeList :: Vector{NTuple{4,T}} + nodeList :: Vector{Vec4{T}} nodeClassList :: Vector{UInt8} # connectivity matrix to represent edges # conn[2,3]=true means node 2 and 3 connected with each other @@ -31,7 +32,7 @@ mutable struct NodeNet{T} connectivityMatrix :: SparseMatrixCSC{Bool,UInt32} end -@inline function NodeNet(nodeList::Vector{NTuple{4,T}}, +@inline function NodeNet(nodeList::Vector{Vec4{T}}, connectivityMatrix::SparseMatrixCSC{Bool, UInt32}) where {T} nodeClassList = zeros(UInt8, length(nodeList)) NodeNet{T}(nodeList, nodeClassList, connectivityMatrix) @@ -43,7 +44,7 @@ Perform the teasar algorithm on the passed binary array. """ function NodeNet( seg::Array{T,3}; obj_id::T = convert(T,1), - expansion::NTuple{3, UInt32} = EXPANSION, + expansion::Vec3{UInt32} = EXPANSION, penalty_fn::Function = alexs_penalty ) where T # note that the object voxels are false and non-object voxels are true! # bin_im = DBFs.create_binary_image( seg, obj_id ) @@ -59,8 +60,8 @@ Return: nodeNet object """ function NodeNet(bin_im::Union{BitArray, Array{Bool,3}}; - offset::NTuple{3, UInt32} = OFFSET, - expansion::NTuple{3, UInt32} = EXPANSION, + offset::Vec3{UInt32} = OFFSET, + expansion::Vec3{UInt32} = EXPANSION, penalty_fn::Function = alexs_penalty) # transform segmentation to points points = PointArrays.from_binary_image(bin_im) @@ -82,7 +83,7 @@ end """ function NodeNet( points::Matrix{T}; dbf::DBF=DBFs.compute_DBF(points), penalty_fn::Function = alexs_penalty, - expansion::NTuple{3, UInt32} = EXPANSION) where T + expansion::Vec3{UInt32} = EXPANSION) where T @assert length(dbf) == size(points, 1) println("total number of points: $(size(points,1))") points, bbox_offset = translate_to_origin!( points ); @@ -159,22 +160,20 @@ function NodeNet( points::Matrix{T}; dbf::DBF=DBFs.compute_DBF(points), nodes, edges = distill!(points, path_nodes, path_edges) conn = get_connectivity_matrix(edges) - nodeList = Vector{NTuple{4,Float32}}() + nodeList = Vector{Vec4f0}() sizehint!(nodeList, length(node_radii)) for i in 1:length(node_radii) - push!(nodeList, (map(Float32,nodes[i,:])..., node_radii[i])) + push!(nodeList, Vec4f0(nodes[i,:]..., node_radii[i])) end nodeNet = NodeNet(nodeList, conn) # add the offset from shift bounding box function - bbox_offset = map(Float32, bbox_offset) - @show bbox_offset add_offset!(nodeNet, bbox_offset) return nodeNet end function NodeNet(swc::SWC) - nodeList = Vector{NTuple{4, Float32}}() + nodeList = Vector{Vec4f0}() nodeClassList = Vector{UInt8}() I = Vector{UInt32}() J = Vector{UInt32}() @@ -231,7 +230,7 @@ the connectivity matrix is symmetric, so the connection is undirected @inline function get_edge_num(self::NodeNet) div(nnz(self.connectivityMatrix), 2) end function get_edges(self::NodeNet) - edges = Vector{Tuple{UInt32,UInt32}}() + edges = Vector{NTuple{2,UInt32}}() conn = get_connectivity_matrix(self) I,J,V = findnz(conn) # only record the triangular part of the connectivity matrix @@ -314,8 +313,8 @@ end end function Base.UnitRange(self::NodeNet) - minCoordinates = [typemax(UInt32), typemax(UInt32), typemax(UInt32)] - maxCoordinates = [zero(UInt32), zero(UInt32), zero(UInt32)] + minCoordinates = Vec3{UInt32}[typemax(UInt32)] + maxCoordinates = Vec3{UInt32}[zero(UInt32)] for node in get_node_list(self) minCoordinates = map(min, minCoordinates, node[1:3]) maxCoordinates = map(max, maxCoordinates, node[1:3]) @@ -412,11 +411,11 @@ function save_edges(self::NodeNet, fileName::String) end ##################### manipulate ############################ -@inline function add_offset!(self::NodeNet{T}, offset::Union{Vector{T},NTuple{3,T}} ) where T +@inline function add_offset!(self::NodeNet{T}, offset::Vec3{T}) where T @assert length(offset) == 3 for i in 1:get_node_num(self) - xyz = map((x,y)->x+y, self.nodeList[i][1:3],offset) - self.nodeList[i] = (xyz..., self.nodeList[i][4]) + xyz = self.nodeList[i][1:3] + offset + self.nodeList[i] = Vec4f0(xyz..., self.nodeList[i][4]) end end @@ -424,7 +423,7 @@ end expansion = (2^(mip-1), 2^(mip-1), 1) stretch_coordinates!(self, expansion) end -@inline function stretch_coordinates!(self::NodeNet, expansion::Union{Vector, Tuple}) +@inline function stretch_coordinates!(self::NodeNet, expansion::Vec3{T}) where T @assert length(expansion) == 3 nodeList = get_node_list(self) for (nodeIndex, node) in enumerate(nodeList) @@ -455,7 +454,7 @@ function translate_to_origin!( points::Matrix{T} ) where T points[:,2] .-= offset[2] points[:,3] .-= offset[3] - points, offset + points, Vec3f0(offset) end """ @@ -495,7 +494,7 @@ end be weighted or modified easily upon return. """ function make_neighbor_graph( points::Array{T,2}, volumeIndex2NodeIndex=nothing, max_dims=nothing; - expansion::NTuple{3, UInt32} = EXPANSION ) where T + expansion::Vec3{UInt32} = EXPANSION ) where T if volumeIndex2NodeIndex == nothing volumeIndex2NodeIndex, max_dims = create_node_lookup(points) end diff --git a/src/PointArrays.jl b/src/PointArrays.jl index 27d412a8..fbadf78b 100644 --- a/src/PointArrays.jl +++ b/src/PointArrays.jl @@ -1,12 +1,14 @@ module PointArrays +import GeometryTypes: Vec3f0, Vec3 + export PointArray const PointArray = Array{UInt32, 2} const ZERO_UINT32 = zero(UInt32) const ONE_UINT32 = one(UInt32) -const DEFAULT_OFFSET = (ZERO_UINT32, ZERO_UINT32, ZERO_UINT32) +const DEFAULT_OFFSET = Vec3(ZERO_UINT32) const MAX_BOUNDARY_DISTANCE = 100000 @@ -35,7 +37,7 @@ end find points inside an object from a segmentation array. """ function from_seg(seg::Array{T,3}; obj_id::T=convert(T,1), - offset::NTuple{3,UInt32}=DEFAULT_OFFSET) where T + offset::Vec3{UInt32}=DEFAULT_OFFSET) where T # compute the number of object voxels nov = 0 for i in eachindex(seg) @@ -93,7 +95,7 @@ end """ add offset to points """ -function add_offset!(self::Array{T,2}, offset::NTuple{3,T}) where T +function add_offset!(self::Array{T,2}, offset::Vec3{T}) where T self[:, 1] .+= offset[1] self[:, 2] .+= offset[2] self[:, 3] .+= offset[3] diff --git a/src/Segments.jl b/src/Segments.jl index 10615757..feaa6dca 100644 --- a/src/Segments.jl +++ b/src/Segments.jl @@ -2,12 +2,11 @@ module Segments import LinearAlgebra: norm, dot import Statistics: mean, std +import GeometryTypes: Vec, Vec4, Vec4f0, Vec3f0 using RealNeuralNetworks.Utils.BoundingBoxes include("Synapses.jl"); using .Synapses -const Node = NTuple{4,Float32} - # TO-DO: make the data type more general const SynapseList = Vector{Union{Missing, Vector{Synapse{Float32}}}} @@ -20,7 +19,7 @@ const UNDEFINED_CLASS = zero(UInt8) export Segment mutable struct Segment{T} # list of tuple (x,y,z,r) - nodeList ::Vector{NTuple{4,T}} + nodeList ::Vector{Vec4{T}} class ::UInt8 #boundingBox ::BoundingBox preSynapseList ::SynapseList @@ -28,11 +27,11 @@ mutable struct Segment{T} end function Segment() - nodeList = Vector{Node}() + nodeList = Vector{Vec4f0}() Segment(nodeList) end -function Segment(nodeList::Vector{Node}; +function Segment(nodeList::Vector{Vec4f0}; class::UInt8=UNDEFINED_CLASS, preSynapseList::Vector{X} = SynapseList(missing, length(nodeList)), postSynapseList::Vector{X} = SynapseList(missing, length(nodeList))) where {X<:Union{Missing, Vector{Synapse{Float32}}}} @@ -63,8 +62,8 @@ end get_nodes_distance(self::Node, other::Node) compute the euclidean distance between two nodes """ -@inline function get_nodes_distance(self::Union{Vector,Tuple}, other::Union{Vector,Tuple}) - norm( [map((x,y)->x-y, self[1:3], other[1:3]) ...]) +@inline function get_nodes_distance(self::Vec{N,T}, other::Vec{N,T}) where {N,T} + norm( self[1:3] .- other[1:3]) end @inline function get_node_list(self::Segment) self.nodeList end @@ -76,25 +75,21 @@ end @inline function get_pre_synapse( self::Segment, index::Int ) self.preSynapseList[index] end @inline function get_post_synapse( self::Segment, index::Int ) self.postSynapseList[index] end -@inline function get_bounding_box_distance(self::Segment, point::Union{Tuple, Vector}) - @assert length(point) >= 3 +@inline function get_bounding_box_distance(self::Segment{T}, point::Vec{N,T}) where {N,T} + @assert N >= 3 boundingBox = get_bounding_box(self) BoundingBoxes.distance_from(boundingBox, point) end -@inline function euclidean_distance( n1::NTuple{3,T}, n2::NTuple{3,T}) where T - norm([map((x,y)->x-y, n1, n2)...]) -end - """ get_path_length(self::Segment; nodeId::Int=length(self)) accumulate the euclidean distance between neighboring nodes """ -@inline function get_path_length(self::Segment; nodeId::Int=length(self)) +@inline function get_path_length(self::Segment; nodeId::Integer=length(self)) ret = 0.0 for i in 2:nodeId - ret += euclidean_distance(self[i][1:3], self[i-1][1:3]) + ret += norm(self[i][1:3] .- self[i-1][1:3]) end ret end @@ -157,7 +152,7 @@ function get_surface_area(self::Segment{T}) where T # average diameter r1 = self[i][4] r2 = self[i-1][4] - h = euclidean_distance(self[i][1:3], self[i-1][1:3]) + h = norm(self[i][1:3] .- self[i-1][1:3]) ret += pi* (r1+r2) * sqrt(h*h + (r1-r2)*(r1-r2)) end ret::T @@ -173,7 +168,7 @@ function get_volume(self::Segment{T}) where T for i in 2:length(self) r1 = self[i-1][4] r2 = self[i][4] - h = euclidean_distance(self[i-1][1:3], self[i][1:3]) + h = norm(self[i-1][1:3] .- self[i][1:3]) ret += pi * h * (r1*r1 + r1*r2 + r2*r2) / T(3) end ret @@ -208,7 +203,7 @@ the ratio of the actual path length to the euclidean distance between head and t pathLength / euclideanLength end -@inline function get_center(nodeList::Vector{Node}) +@inline function get_center(nodeList::Vector{Vec4f0}) (map(i->mean(y->y[i], nodeList), 1:4)...,) end @@ -321,18 +316,15 @@ end """ distance from a point """ -function distance_from(self::Segment, point::Tuple) - distance_from(self, [point[1:3]...]) -end -function distance_from(self::Segment{T}, point::Vector) where T +function distance_from(self::Segment{T}, point::Vec{N,T}) where {N,T} @assert !isempty(self) ret = (zero(T), zero(Int)) nodeList = get_node_list(self) @assert !isempty(nodeList) - @assert length(point) == 3 || length(point) == 4 + @assert N == 3 || N == 4 distance = typemax(T) for (index, node) in enumerate(nodeList) - d = norm( [node[1:3]...] .- [point[1:3]...] ) + d = norm( node[1:3] .- point[1:3] ) if d < distance distance = d ret = (d, index) @@ -389,9 +381,9 @@ function adjust_class!(self::Segment) end end -function add_offset(self::Segment{T}, offset::Union{Tuple, Vector}) where T +function add_offset(self::Segment{T}, offset::Vec3f0) where T @assert length(offset) == 3 - nodeList = Vector{NTuple{4,T}}() + nodeList = Vec4f0[] for node in self.nodeList newNode = map(+, node, [offset..., zero(T)]) push!(nodeList, newNode) @@ -413,7 +405,7 @@ function remove_nodes(self::Segment{T}, removeIdRange::UnitRange{Int}) where T if newLength == 0 return Segment() end - newNodeList = Vector{NTuple{4, T}}() + newNodeList = Vector{Vec4f0}() sizehint!(newNodeList, newLength) for (index,node) in enumerate(get_node_list(self)) @@ -443,7 +435,7 @@ remove neighboring nodes that is the same. """ function remove_redundent_nodes!(self::Segment{T}) where T nodeList = get_node_list(self) - newNodeList = Vector{NTuple{4, T}}() + newNodeList = Vector{Vec4f0}() newPreSynapseList = SynapseList() newPostSynapseList = SynapseList() for index in 1:length(nodeList)-1 diff --git a/src/Synapses.jl b/src/Synapses.jl index 82a94302..e4fdd8d5 100644 --- a/src/Synapses.jl +++ b/src/Synapses.jl @@ -2,20 +2,21 @@ module Synapses using RealNeuralNetworks.Utils.SynapseTables using RealNeuralNetworks.Utils.BoundingBoxes using DataFrames +import GeometryTypes: Vec3, Vec3f0 export Synapse mutable struct Synapse{T} psdSegmentationId ::Int - psdCoordinate ::NTuple{3,T} # nm + psdCoordinate ::Vec3{T} # nm psdBoundingBox ::BoundingBox psdSize ::Int preSynapticSegmentationId ::Int - preSynapticCoordinate ::NTuple{3,T} # nm - preSynapticWeight ::Float32 + preSynapticCoordinate ::Vec3{T} # nm + preSynapticWeight ::T postSynapticSegmentationId ::Int - postSynapticCoordinate ::NTuple{3,T} # nm - postSynapticWeight ::Float32 + postSynapticCoordinate ::Vec3{T} # nm + postSynapticWeight ::T end function Synapse( row::DataFrameRow ) diff --git a/src/Utils/BoundingBoxes.jl b/src/Utils/BoundingBoxes.jl index 7165bf78..6f177caf 100644 --- a/src/Utils/BoundingBoxes.jl +++ b/src/Utils/BoundingBoxes.jl @@ -1,37 +1,35 @@ module BoundingBoxes import LinearAlgebra: norm +import GeometryTypes: Vec, Vec3, Vec3f0, Vec4 const ZERO_FLOAT32 = zero(Float32) export BoundingBox -mutable struct BoundingBox - minCorner ::NTuple{3,Float32} - maxCorner ::NTuple{3,Float32} +mutable struct BoundingBox{T} + minCorner ::Vec3{T} + maxCorner ::Vec3{T} end function BoundingBox() - minCorner = (Inf32, Inf32, Inf32) - maxCorner = (ZERO_FLOAT32, ZERO_FLOAT32, ZERO_FLOAT32) + minCorner = Vec3f0(Inf32) + maxCorner = Vec3f0(ZERO_FLOAT32) BoundingBox( minCorner, maxCorner ) end -function BoundingBox( minCorner::NTuple{3,Int}, maxCorner::NTuple{3,Int} ) - BoundingBox( map(x->convert(Float32,x), minCorner), - map(x->convert(Float32,x), maxCorner) ) -end - -function BoundingBox( minCorner::Vector{Float32}, maxCorner::Vector{Float32} ) - BoundingBox((minCorner...,), (maxCorner...,)) +function BoundingBox(minCorner::Union{Vector{T}, NTuple{3,T}}, + maxCorner::Union{Vector{T}, NTuple{3,T}} ) where T + BoundingBox(Vec3{T}(minCorner), + Vec3{T}(maxCorner) ) end """ get bounding box from a node list """ function BoundingBox(nodeList::Union{Vector, Set}) - minCorner = fill(Inf32,3) - maxCorner = fill(ZERO_FLOAT32, 3) + minCorner = Vec3f0(Inf32) + maxCorner = Vec3f0(ZERO_FLOAT32) for node in nodeList minCorner = map(min, minCorner, node[1:3]) maxCorner = map(max, maxCorner, node[1:3]) @@ -44,7 +42,8 @@ function get_unit_range( self::BoundingBox ) end function Base.size(self::BoundingBox) - map(length, get_unit_range(self)) + sz = map(length, get_unit_range(self)) + return (sz...,) end function Base.isequal(self::BoundingBox, other::BoundingBox) @@ -64,14 +63,13 @@ function Base.union(self::BoundingBox, other::BoundingBox) BoundingBox(minCorner, maxCorner) end -function isinside(self::BoundingBox, point::Union{Tuple, Vector}) +function isinside(self::BoundingBox{T}, point::Vec{N,T}) where {N,T} all( map((x,y)->x>y, self.maxCorner, point[1:3]) ) && all( map((x,y)->x max(cmin-p, Float32(0), p-cmax), self.minCorner, point[1:3], self.maxCorner ) norm([d...]) diff --git a/test/Neurons.jl b/test/Neurons.jl index 14338bf9..45898287 100644 --- a/test/Neurons.jl +++ b/test/Neurons.jl @@ -8,12 +8,13 @@ using RealNeuralNetworks.Utils.SynapseTables using CSV using SparseArrays using LinearAlgebra -using DataFrames +using DataFrames +import GeometryTypes: Vec3f0 const NEURON_ID = 77625 const ASSERT_DIR = joinpath(@__DIR__, "../asset") const SWC_BIN_PATH = joinpath(ASSERT_DIR, "$(NEURON_ID).swc.bin") -const ARBOR_DENSITY_MAP_VOXEL_SIZE = (2000,2000,2000) +const ARBOR_DENSITY_MAP_VOXEL_SIZE = Vec3f0(2000) @testset "test basic functions of Neurons..." begin swc = SWCs.load(joinpath(ASSERT_DIR, "example.swc")) diff --git a/test/NodeNets.jl b/test/NodeNets.jl index 82e5c9c0..894d8daa 100644 --- a/test/NodeNets.jl +++ b/test/NodeNets.jl @@ -7,14 +7,16 @@ using RealNeuralNetworks.NodeNets using RealNeuralNetworks.SWCs using LinearAlgebra +import GeometryTypes: Vec3f0, Vec3 + const CELL_ID = UInt32(76880) -const EXPANSION= (UInt32(80), UInt32(80), UInt32(40)) +const EXPANSION= Vec3(UInt32(80), UInt32(80), UInt32(40)) const MIP = 4 function get_seg_from_h5() # read seg data - f = h5open("/usr/people/jingpeng/seungmount/research/Ashwin/Scripts/NG_scripts/77605.h5") + f = h5open("../asset/77605.h5") seg = f["main"][1,:,:,101:500] close(f) seg = reshape(seg, size(seg)[2:4]) @@ -49,7 +51,7 @@ end println("building nodeNet ...") #@time nodeNet = NodeNet( seg; obj_id = CELL_ID ) @time nodeNet = NodeNet( seg; obj_id = one(UInt32) ) - NodeNets.add_offset!(nodeNet, (-one(Float32),-one(Float32),-one(Float32))) + NodeNets.add_offset!(nodeNet, Vec3f0(-one(Float32))) bin = NodeNets.get_neuroglancer_precomputed(nodeNet) # open("/tmp/fake.bin", "w") do f write(f, bin) end @time swc = SWCs.SWC( nodeNet ) diff --git a/test/PointArrays.jl b/test/PointArrays.jl index 4b4ee3a7..ece975e2 100644 --- a/test/PointArrays.jl +++ b/test/PointArrays.jl @@ -1,6 +1,8 @@ using RealNeuralNetworks using RealNeuralNetworks.NodeNets.PointArrays +import GeometryTypes: Vec3 + using Test @testset "test PointArray add_offset! operation" begin @@ -9,7 +11,7 @@ using Test 0xa8e1555e 0x7fbeb45b 0xc786dc67 0x940649cd 0x20b148bc 0xc15152cb 0xc91b9db1 0x39fee982 0x7019e841] - offset = (UInt32(1), UInt32(2), UInt32(3)) + offset = Vec3(UInt32(1), UInt32(2), UInt32(3)) PointArrays.add_offset!(a, offset) b = [0xc95f39a8 0x21a2a09e 0xc9e447a5 0xa8e1555f 0x7fbeb45d 0xc786dc6a diff --git a/test/Utils/BoundingBoxes.jl b/test/Utils/BoundingBoxes.jl index 948b25ff..386da1de 100644 --- a/test/Utils/BoundingBoxes.jl +++ b/test/Utils/BoundingBoxes.jl @@ -1,5 +1,6 @@ using Test using RealNeuralNetworks.Utils.BoundingBoxes +import GeometryTypes: Vec3 function create_fake_bounding_box() bbox = BoundingBox((1,1,1), (204,204,300)) @@ -7,10 +8,10 @@ end @testset "test bounding box" begin bbox = BoundingBox((1,2,3), (2,3,4)) - d = BoundingBoxes.distance_from(bbox, (3,4,5)) + d = BoundingBoxes.distance_from(bbox, Vec3(3,4,5)) @test d ≈ sqrt(3) range = BoundingBoxes.get_unit_range(bbox) - @test range == (1:2, 2:3, 3:4) + @test range == UnitRange[1:2, 2:3, 3:4] @test size(bbox) == (2,2,2)