Merge pull request #58 from JuliaGeometry/sjk/clean3

Sjk/clean3

Author: sjkelly

Project.toml

@@ -12,6 +12,7 @@ julia = "0.7, 1"
 [extras]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [targets]
-test = ["LinearAlgebra", "Test"]
+test = ["LinearAlgebra", "StatsBase", "Test"]

src/Contour.jl

@@ -2,6 +2,8 @@ module Contour
 
 using StaticArrays
 
+include("interpolate.jl")
+
 export
     ContourLevel,
     Curve2,
@@ -144,42 +146,48 @@ const WN, WS, WE = W|N, W|S, W|E
 const NWSE = NW | 0x10 # special (ambiguous case)
 const NESW = NE | 0x10 # special (ambiguous case)
 
-const dirStr = ["N", "S", "NS", "E", "NE", "NS", "Invalid crossing",
-                "W", "NW", "SW", "Invalid crossing", "WE"]
+# Maps cell type to crossing types for non-ambiguous cells
+const edge_LUT = (SW, SE, EW, NE, 0x0, NS, NW, NW, NS, 0x0, NE, EW, SE, SW)
 
 # The way a contour crossing goes through a cell is labeled
 # by combining compass directions (e.g. a NW crossing connects
 # the N edge and W edges of the cell).  The Cell type records
 # the type of crossing that a cell contains.  While most
 # cells will have only one crossing, cell type 5 and 10 will
 # have two crossings.
-function get_next_edge!(cells::Dict, xi, yi, entry_edge::UInt8)
-    key = (xi,yi)
-    cell = cells[key]
+function get_next_edge!(cells::Dict, key, entry_edge::UInt8)
+    cell = pop!(cells, key)
     if cell == NWSE
-        if !iszero(NW & entry_edge)
+        if entry_edge == N || entry_edge == W
             cells[key] = SE
-            return NW ⊻ entry_edge
-        else #Nw
+            cell = NW
+        else #SE
             cells[key] = NW
-            return SE ⊻ entry_edge
+            cell = SE
         end
     elseif cell == NESW
-        if !iszero(NE & entry_edge)
+        if entry_edge == N || entry_edge == E
             cells[key] = SW
-            return NE ⊻ entry_edge
+            cell = NE
         else #SW
             cells[key] = NE
-            return SW ⊻ entry_edge
+            cell = SW
         end
-    else
-        next_edge = cell ⊻ entry_edge
-        delete!(cells, key)
-        return next_edge
     end
+    return cell ⊻ entry_edge
 end
 
-function get_first_crossing(cell)
+# N, S, E, W
+const next_map = ((0,1), (0,-1), (1,0), (-1,0))
+const next_edge = (S,N,W,E)
+
+@inline function advance_edge(ind, edge)
+    n = trailing_zeros(edge) + 1
+    nt = ind .+ next_map[n]
+    return nt, next_edge[n]
+end
+
+@inline function get_first_crossing(cell)
     if cell == NWSE
         return NW
     elseif cell == NESW
@@ -189,10 +197,7 @@ function get_first_crossing(cell)
     end
 end
 
-# Maps cell type to crossing types for non-ambiguous cells
-const edge_LUT = (SW, SE, EW, NE, 0x0, NS, NW, NW, NS, 0x0, NE, EW, SE, SW)
-
-function _get_case(z, h)
+@inline function _get_case(z, h)
     case = z[1] > h ? 0x01 : 0x00
     z[2] > h && (case |= 0x02)
     z[3] > h && (case |= 0x04)
@@ -229,67 +234,12 @@ function get_level_cells(z, h::Number)
     return cells
 end
 
-# Some constants used by trace_contour
-
-const fwd, rev = (false, true)
-
-function add_vertex!(curve::Curve2{T}, pos::Tuple{T,T}, dir) where {T}
-    if dir == fwd
-        push!(curve.vertices, SVector{2,T}(pos...))
-    else
-        pushfirst!(curve.vertices, SVector{2,T}(pos...))
-    end
-end
-
-# Given the row and column indices of the lower left
-# vertex, add the location where the contour level
-# crosses the specified edge.
-function interpolate(x, y, z::AbstractMatrix{T}, h::Number, xi::Int, yi::Int, edge::UInt8) where {T <: AbstractFloat}
-    @inbounds if edge == W
-        y_interp = y[yi] + (y[yi + 1] - y[yi]) * (h - z[xi, yi]) / (z[xi, yi + 1] - z[xi, yi])
-        x_interp = x[xi]
-    elseif edge == E
-        y_interp = y[yi] + (y[yi + 1] - y[yi]) * (h - z[xi + 1, yi]) / (z[xi + 1, yi + 1] - z[xi + 1, yi])
-        x_interp = x[xi + 1]
-    elseif edge == N
-        y_interp = y[yi + 1]
-        x_interp = x[xi] + (x[xi + 1] - x[xi]) * (h - z[xi, yi + 1]) / (z[xi + 1, yi + 1] - z[xi, yi + 1])
-    elseif edge == S
-        y_interp = y[yi]
-        x_interp = x[xi] + (x[xi + 1] - x[xi]) * (h - z[xi, yi]) / (z[xi + 1, yi] - z[xi, yi])
-    end
-
-    return x_interp, y_interp
-end
-
-function interpolate(x::AbstractMatrix{T}, y::AbstractMatrix{T}, z::AbstractMatrix{T}, h::Number, xi::Int, yi::Int, edge::UInt8) where {T <: AbstractFloat}
-    if edge == W
-        Δ = [y[xi,  yi+1] - y[xi,  yi  ], x[xi,  yi+1] - x[xi,  yi  ]].*(h - z[xi,  yi  ])/(z[xi,  yi+1] - z[xi,  yi  ])
-        y_interp = y[xi,yi] + Δ[1]
-        x_interp = x[xi,yi] + Δ[2]
-    elseif edge == E
-        Δ = [y[xi+1,yi+1] - y[xi+1,yi  ], x[xi+1,yi+1] - x[xi+1,yi  ]].*(h - z[xi+1,yi  ])/(z[xi+1,yi+1] - z[xi+1,yi  ])
-        y_interp = y[xi+1,yi] + Δ[1]
-        x_interp = x[xi+1,yi] + Δ[2]
-    elseif edge == N
-        Δ = [y[xi+1,yi+1] - y[xi,  yi+1], x[xi+1,yi+1] - x[xi,  yi+1]].*(h - z[xi,  yi+1])/(z[xi+1,yi+1] - z[xi,  yi+1])
-        y_interp = y[xi,yi+1] + Δ[1]
-        x_interp = x[xi,yi+1] + Δ[2]
-    elseif edge == S
-        Δ = [y[xi+1,yi  ] - y[xi,  yi  ], x[xi+1,yi  ] - x[xi,  yi  ]].*(h - z[xi,  yi  ])/(z[xi+1,yi  ] - z[xi,  yi  ])
-        y_interp = y[xi,yi] + Δ[1]
-        x_interp = x[xi,yi] + Δ[2]
-    end
-
-    return x_interp, y_interp
-end
-
 
 # Given a cell and a starting edge, we follow the contour line until we either
 # hit the boundary of the input data, or we form a closed contour.
-function chase!(cells, curve, x, y, z, h, xi_start, yi_start, entry_edge, xi_max, yi_max, dir)
+function chase!(cells, curve, x, y, z, h, start, entry_edge, xi_max, yi_max, ::Type{VT}) where VT
 
-    xi, yi = xi_start, yi_start
+    ind = start
 
     # When the contour loops back to the starting cell, it is possible
     # for it to not intersect with itself.  This happens if the starting
@@ -298,92 +248,62 @@ function chase!(cells, curve, x, y, z, h, xi_start, yi_start, entry_edge, xi_max
     loopback_edge = entry_edge
 
     @inbounds while true
-        exit_edge = get_next_edge!(cells, xi, yi, entry_edge)
-
-        add_vertex!(curve, interpolate(x, y, z, h, xi, yi, exit_edge), dir)
-
-        if exit_edge == N
-            yi += 1
-            entry_edge = S
-        elseif exit_edge == S
-            yi -= 1
-            entry_edge = N
-        elseif exit_edge == E
-            xi += 1
-            entry_edge = W
-        elseif exit_edge == W
-            xi -= 1
-            entry_edge = E
-        end
-        !((xi, yi, entry_edge) != (xi_start, yi_start, loopback_edge) &&
-           0 < yi < yi_max && 0 < xi < xi_max) && break
+        exit_edge = get_next_edge!(cells, ind, entry_edge)
+
+        push!(curve, interpolate(x, y, z, h, ind, exit_edge, VT))
+
+        ind, entry_edge = advance_edge(ind, exit_edge)
+
+        !((ind[1], ind[2], entry_edge) != (start[1], start[2], loopback_edge) &&
+           0 < ind[2] < yi_max && 0 < ind[1] < xi_max) && break
     end
 
-    return xi, yi
+    return ind
 end
 
 
 function trace_contour(x, y, z, h::Number, cells::Dict)
 
     contours = ContourLevel(h)
 
-    (xi_max, yi_max) = size(z)
+    (xi_max, yi_max) = size(z)::Tuple{Int,Int}
+
+    VT = SVector{2,promote_type(map(eltype, (x, y, z))...)}
 
     # When tracing out contours, this algorithm picks an arbitrary
     # starting cell, then first follows the contour in one direction
     # until it either ends up where it started # or at one of the boundaries.
     # It then tries to trace the contour in the opposite direction.
 
     @inbounds while length(cells) > 0
-        contour = Curve2(promote_type(map(eltype, (x, y, z))...))
+        contour_arr = VT[]
 
         # Pick initial box
-        (xi_0, yi_0), cell = first(cells)
-        (xi, yi) = (xi_0, yi_0)
+        ind, cell = first(cells)
 
         # Pick a starting edge
         crossing = get_first_crossing(cell)
-        starting_edge = UInt8(0)
-        for edge in (N, S, E, W)
-            if !iszero(edge & crossing)
-                starting_edge = edge
-                break
-            end
-        end
+        starting_edge = 0x01 << trailing_zeros(crossing)
 
         # Add the contour entry location for cell (xi_0,yi_0)
-        add_vertex!(contour, interpolate(x, y, z, h, xi_0, yi_0, starting_edge), fwd)
+        push!(contour_arr, interpolate(x, y, z, h, ind, starting_edge, VT))
 
         # Start trace in forward direction
-        (xi_end, yi_end) = chase!(cells, contour, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max, fwd)
+        ind_end = chase!(cells, contour_arr, x, y, z, h, ind, starting_edge, xi_max, yi_max, VT)
 
-        if (xi_end, yi_end) == (xi_0, yi_0)
-            push!(contours.lines, contour)
+        if ind == ind_end
+            push!(contours.lines, Curve2(contour_arr))
             continue
         end
 
-        if starting_edge == N
-            yi = yi_0 + 1
-            starting_edge = S
-        elseif starting_edge == S
-            yi = yi_0 - 1
-            starting_edge = N
-        elseif starting_edge == E
-            xi = xi_0 + 1
-            starting_edge = W
-        elseif starting_edge == W
-            xi = xi_0 - 1
-            starting_edge = E
-        end
+        ind, starting_edge = advance_edge(ind, starting_edge)
 
-        if !(0 < yi < yi_max && 0 < xi < xi_max)
-            push!(contours.lines, contour)
-            continue
+        if 0 < ind[2] < yi_max && 0 < ind[1] < xi_max
+            # Start trace in reverse direction
+            chase!(cells, reverse!(contour_arr), x, y, z, h, ind, starting_edge, xi_max, yi_max, VT)
         end
 
-        # Start trace in reverse direction
-        (xi, yi) = chase!(cells, contour, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max, rev)
-        push!(contours.lines, contour)
+        push!(contours.lines, Curve2(contour_arr))
     end
 
     return contours

src/interpolate.jl

@@ -0,0 +1,63 @@
+# Given the row and column indices of the lower left
+# vertex, add the location where the contour level
+# crosses the specified edge.
+function interpolate(x, y, z::AbstractMatrix, h::Number, ind, edge::UInt8, ::Type{VT}) where {VT}
+    xi, yi = ind
+    @inbounds if edge == W
+        y_interp = y[yi] + (y[yi + 1] - y[yi]) * (h - z[xi, yi]) / (z[xi, yi + 1] - z[xi, yi])
+        x_interp = x[xi]
+    elseif edge == E
+        y_interp = y[yi] + (y[yi + 1] - y[yi]) * (h - z[xi + 1, yi]) / (z[xi + 1, yi + 1] - z[xi + 1, yi])
+        x_interp = x[xi + 1]
+    elseif edge == N
+        y_interp = y[yi + 1]
+        x_interp = x[xi] + (x[xi + 1] - x[xi]) * (h - z[xi, yi + 1]) / (z[xi + 1, yi + 1] - z[xi, yi + 1])
+    elseif edge == S
+        y_interp = y[yi]
+        x_interp = x[xi] + (x[xi + 1] - x[xi]) * (h - z[xi, yi]) / (z[xi + 1, yi] - z[xi, yi])
+    end
+
+    return VT(x_interp, y_interp)
+end
+
+function interpolate(x::AbstractRange, y::AbstractRange, z::AbstractMatrix, h::Number, ind, edge::UInt8, ::Type{VT}) where {VT}
+    xi, yi = ind
+    @inbounds if edge == W
+        y_interp = y[yi] + step(y) * (h - z[xi, yi]) / (z[xi, yi + 1] - z[xi, yi])
+        x_interp = x[xi]
+    elseif edge == E
+        y_interp = y[yi] + step(y) * (h - z[xi + 1, yi]) / (z[xi + 1, yi + 1] - z[xi + 1, yi])
+        x_interp = x[xi + 1]
+    elseif edge == N
+        y_interp = y[yi + 1]
+        x_interp = x[xi] + step(x) * (h - z[xi, yi + 1]) / (z[xi + 1, yi + 1] - z[xi, yi + 1])
+    elseif edge == S
+        y_interp = y[yi]
+        x_interp = x[xi] + step(x) * (h - z[xi, yi]) / (z[xi + 1, yi] - z[xi, yi])
+    end
+
+    return VT(x_interp, y_interp)
+end
+
+function interpolate(x::AbstractMatrix, y::AbstractMatrix, z::AbstractMatrix, h::Number, ind, edge::UInt8, ::Type{VT}) where {VT}
+    xi, yi = ind
+    @inbounds if edge == W
+        Δ = [y[xi,  yi+1] - y[xi,  yi  ], x[xi,  yi+1] - x[xi,  yi  ]].*(h - z[xi,  yi  ])/(z[xi,  yi+1] - z[xi,  yi  ])
+        y_interp = y[xi,yi] + Δ[1]
+        x_interp = x[xi,yi] + Δ[2]
+    elseif edge == E
+        Δ = [y[xi+1,yi+1] - y[xi+1,yi  ], x[xi+1,yi+1] - x[xi+1,yi  ]].*(h - z[xi+1,yi  ])/(z[xi+1,yi+1] - z[xi+1,yi  ])
+        y_interp = y[xi+1,yi] + Δ[1]
+        x_interp = x[xi+1,yi] + Δ[2]
+    elseif edge == N
+        Δ = [y[xi+1,yi+1] - y[xi,  yi+1], x[xi+1,yi+1] - x[xi,  yi+1]].*(h - z[xi,  yi+1])/(z[xi+1,yi+1] - z[xi,  yi+1])
+        y_interp = y[xi,yi+1] + Δ[1]
+        x_interp = x[xi,yi+1] + Δ[2]
+    elseif edge == S
+        Δ = [y[xi+1,yi  ] - y[xi,  yi  ], x[xi+1,yi  ] - x[xi,  yi  ]].*(h - z[xi,  yi  ])/(z[xi+1,yi  ] - z[xi,  yi  ])
+        y_interp = y[xi,yi] + Δ[1]
+        x_interp = x[xi,yi] + Δ[2]
+    end
+
+    return VT(x_interp, y_interp)
+end

test/verify_vertices.jl

@@ -215,4 +215,15 @@ for i in eachindex(cts_ct)
     @test length(cts.contours[i].lines) == cts_ct[i]
     @test all(lines_ct[i] .== [length(c.vertices) for c in cts.contours[i].lines])
 end
+
+
+# support non-float z
+using StatsBase
+
+N = 10000
+x = randn(N)
+y = randn(N)
+H = fit(Histogram, (x, y), closed = :left)
+contours(midpoints.(H.edges)..., H.weights)
+
 end