Use arrays rather than Dicts for Contour Chasing

src/Contour.jl

@@ -62,9 +62,17 @@ over the result.
 """
 function contour(x, y, z, level::Number)
     # Todo: size checking on x,y,z
-    trace_contour(x, y, z, level, get_level_cells(z, level))
+    cells, cell_pop = get_level_cells(z, level)
+    trace_contour(x, y, z, level, cells, cell_pop)
 end
 
+function contour(cells, x, y, z, level::Number)
+    # Todo: size checking on x,y,z
+    cell_pop = get_level_cells!(cells, z, level)
+    trace_contour(x, y, z, level, cells, cell_pop)
+end
+
+
 """
 `contours` returns a set of isolines.
 
@@ -76,7 +84,10 @@ contours(::Any...) = error("This method exists only for documentation purposes")
 `contours(x,y,z,levels)` Trace the contour levels indicated by the `levels`
 argument.
 """
-contours(x, y, z, levels) = ContourCollection([contour(x, y, z, l) for l in levels])
+function contours(x, y, z, levels)
+    cells = cell_matrix(z)
+    ContourCollection([contour(cells, x, y, z, l) for l in levels])
+end
 
 """
 `contours(x,y,z,Nlevels::Integer)` Trace `Nlevels` contour levels at heights
@@ -141,6 +152,8 @@ const NS, NE, NW = N|S, N|E, N|W
 const SN, SE, SW = S|N, S|E, S|W
 const EN, ES, EW = E|N, E|S, E|W
 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"]
@@ -151,24 +164,40 @@ const dirStr = ["N", "S", "NS", "E", "NE", "NS", "Invalid crossing",
 # 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.
-struct Cell
-    crossings::Vector{UInt8}
-end
-
-function get_next_edge!(cell::Cell, entry_edge::UInt8)
-    for (i, edge) in enumerate(cell.crossings)
-        if !iszero(edge & entry_edge)
-            next_edge = edge ⊻ entry_edge
-            deleteat!(cell.crossings, i)
-
-            return next_edge
+function get_next_edge!(cells::Array, xi, yi, entry_edge::UInt8, cell_pop)
+    cell = cells[xi,yi]
+    if cell != NWSE && cell != NESW
+        next_edge = cell ⊻ entry_edge
+        cells[xi,yi] = 0x0
+        return next_edge, cell_pop - 1
+    else  # ambiguous case flag
+        if cell == NWSE
+            if !iszero(NW & entry_edge)
+                cells[xi,yi] = SE
+                return NW ⊻ entry_edge, cell_pop
+            elseif !iszero(SE & entry_edge)
+                cells[xi,yi] = NW
+                return SE ⊻ entry_edge, cell_pop
+            end
+        elseif cell == NESW
+            if !iszero(NE & entry_edge)
+                cells[xi,yi] = SW
+                return NE ⊻ entry_edge, cell_pop
+            elseif !iszero(SW & entry_edge)
+                cells[xi,yi] = NE
+                return SW ⊻ entry_edge, cell_pop
+            end
         end
     end
-    error("There is no edge containing ", entry_edge)
+end
+
+function get_first_crossing(cell)
+    cell & 0x0f
 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)
+const start_edge_LUT = (0x00, 0x00, 0x01, 0x00, 0x01, 0x02, 0x00, 0x00, 0x01, 0x02, 0x00, 0x04, 0x00, 0x00)
 
 function _get_case(z, h)
     case = z[1] > h ? 0x01 : 0x00
@@ -178,12 +207,24 @@ function _get_case(z, h)
     case
 end
 
-function get_level_cells(z, h::Number)
-    cells = Dict{(Tuple{Int,Int}),Cell}()
+function cell_matrix(z)
+    xi_max, yi_max = size(z)
+    cells = zeros(UInt8, (xi_max-1,yi_max-1))
+end
+
+function get_level_cells(z,h::Number)
+    cells = cell_matrix(z)
+    cell_pop = get_level_cells!(cells, z, h)
+    cells, cell_pop
+end
+
+function get_level_cells!(cells, z, h::Number)
     xi_max, yi_max = size(z)
 
-    @inbounds for xi in 1:xi_max - 1
-        for yi in 1:yi_max - 1
+    cell_pop = 0
+
+    @inbounds for yi in 1:yi_max - 1
+        for xi in 1:xi_max - 1
             elts = (z[xi, yi], z[xi + 1, yi], z[xi + 1, yi + 1], z[xi, yi + 1])
             case = _get_case(elts, h)
 
@@ -192,34 +233,44 @@ function get_level_cells(z, h::Number)
                 continue
             end
 
+            cell_pop += 1 # number of cells in the array populated
+
             # Process ambiguous cells (case 5 and 10) using
             # a bilinear interpolation of the cell-center value.
             if case == 0x05
                 if 0.25*sum(elts) >= h
-                    cells[(xi, yi)] = Cell([NW, SE])
+                    cells[xi, yi] = NWSE
                 else
-                    cells[(xi, yi)] = Cell([NE, SW])
+                    cells[xi, yi] = NESW
                 end
             elseif case == 0x0a
                 if 0.25*sum(elts) >= h
-                    cells[(xi, yi)] = Cell([NE, SW])
+                    cells[xi, yi] = NESW
                 else
-                    cells[(xi, yi)] = Cell([NW, SE])
+                    cells[xi, yi] = NWSE
                 end
             else
-                cells[(xi, yi)] = Cell([edge_LUT[case]])
+                cells[xi, yi] = edge_LUT[case]
             end
         end
     end
 
-    return cells
+    return cell_pop
+end
+
+function findfirst_cell(m, from_x, from_y)
+    s = size(m)
+    @inbounds for xi = from_x:s[1], yi = from_y:s[2]
+        !iszero(m[xi,yi]) && return xi,yi
+    end
+    return 0,0
 end
 
 # Some constants used by trace_contour
 
 const fwd, rev = (UInt8(0)), (UInt8(1))
 
-@inline function add_vertex!(curve::Curve2{T}, pos::Tuple{T,T}, dir::UInt8) where {T}
+function add_vertex!(curve::Curve2{T}, pos::Tuple{T,T}, dir::UInt8) where {T}
     if dir == fwd
         push!(curve.vertices, SVector{2,T}(pos...))
     else
@@ -273,7 +324,7 @@ 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, xi_start, yi_start, entry_edge, xi_max, yi_max, cell_pop, dir)
 
     xi, yi = xi_start, yi_start
 
@@ -284,11 +335,7 @@ function chase!(cells, curve, x, y, z, h, xi_start, yi_start, entry_edge, xi_max
     loopback_edge = entry_edge
 
     while true
-        cell = cells[(xi, yi)]
-        exit_edge = get_next_edge!(cell, entry_edge)
-        if length(cell.crossings) == 0
-            delete!(cells, (xi, yi))
-        end
+        exit_edge, cell_pop = get_next_edge!(cells, xi, yi, entry_edge, cell_pop)
 
         add_vertex!(curve, interpolate(x, y, z, h, xi, yi, exit_edge), dir)
 
@@ -309,51 +356,43 @@ function chase!(cells, curve, x, y, z, h, xi_start, yi_start, entry_edge, xi_max
            0 < yi < yi_max && 0 < xi < xi_max) && break
     end
 
-    return xi, yi
+    return xi, yi, cell_pop
 end
 
 
-function trace_contour(x, y, z, h::Number, cells::Dict{(Tuple{Int,Int}),Cell})
+function trace_contour(x, y, z, h::Number, cells::Array, cell_pop)
 
     contours = ContourLevel(h)
 
-    local yi::Int
-    local xi::Int
-    local xi_0::Int
-    local yi_0::Int
-
-    local xi_max::Int
-    local yi_max::Int
-
     (xi_max, yi_max) = size(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.
 
-    while length(cells) > 0
+    nonempty_cells = cell_pop
+    (xi_0, yi_0) = (1,1)
+
+    @inbounds while nonempty_cells > 0
         contour = Curve2(promote_type(map(eltype, (x, y, z))...))
 
         # Pick initial box
-        (xi_0, yi_0), cell = first(cells)
+        (xi_0, yi_0) = findfirst_cell(cells, xi_0, yi_0)
+        iszero(xi_0) && iszero(yi_0) && break
         (xi, yi) = (xi_0, yi_0)
+        cell = cells[xi,yi]
 
         # Pick a starting edge
-        crossing = first(cell.crossings)
-        starting_edge = UInt8(0)
-        for edge in (N, S, E, W)
-            if !iszero(edge & crossing)
-                starting_edge = edge
-                break
-            end
-        end
+        crossing = get_first_crossing(cell)
+        starting_edge = start_edge_LUT[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)
 
         # 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)
+        xi_end, yi_end, cp = chase!(cells, contour, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max, nonempty_cells, fwd)
+        nonempty_cells = cp
 
         if (xi_end, yi_end) == (xi_0, yi_0)
             push!(contours.lines, contour)
@@ -380,7 +419,8 @@ function trace_contour(x, y, z, h::Number, cells::Dict{(Tuple{Int,Int}),Cell})
         end
 
         # Start trace in reverse direction
-        (xi, yi) = chase!(cells, contour, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max, rev)
+        xi, yi, cp = chase!(cells, contour, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max, nonempty_cells, rev)
+        nonempty_cells = cp
         push!(contours.lines, contour)
     end