From 5714816534c26a4bea61c8217e81cd8d895bf470 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sun, 12 Apr 2020 19:33:38 -0400 Subject: [PATCH 01/10] cleanup forward/reverse logic --- src/Contour.jl | 49 +++++++++++++++++-------------------------------- 1 file changed, 17 insertions(+), 32 deletions(-) diff --git a/src/Contour.jl b/src/Contour.jl index 7cdfd72..5414b64 100644 --- a/src/Contour.jl +++ b/src/Contour.jl @@ -144,9 +144,6 @@ 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"] - # 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 @@ -192,7 +189,7 @@ 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,18 +226,6 @@ 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. @@ -259,7 +244,7 @@ function interpolate(x, y, z::AbstractMatrix{T}, h::Number, xi::Int, yi::Int, ed 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 + return SVector{2,T}(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} @@ -281,13 +266,13 @@ function interpolate(x::AbstractMatrix{T}, y::AbstractMatrix{T}, z::AbstractMatr x_interp = x[xi,yi] + Δ[2] end - return x_interp, y_interp + return SVector{2,T}(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, xi_start, yi_start, entry_edge, xi_max, yi_max) xi, yi = xi_start, yi_start @@ -300,7 +285,7 @@ function chase!(cells, curve, x, y, z, h, xi_start, yi_start, entry_edge, xi_max @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) + push!(curve, interpolate(x, y, z, h, xi, yi, exit_edge)) if exit_edge == N yi += 1 @@ -327,7 +312,9 @@ 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 @@ -335,7 +322,7 @@ function trace_contour(x, y, z, h::Number, cells::Dict) # 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) @@ -352,13 +339,13 @@ function trace_contour(x, y, z, h::Number, cells::Dict) end # 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, xi_0, yi_0, starting_edge)) # 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) = chase!(cells, contour_arr, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max) - if (xi_end, yi_end) == (xi_0, yi_0) - push!(contours.lines, contour) + if xi_end == xi_0 && yi_end == yi_0 + push!(contours.lines, Curve2(contour_arr)) continue end @@ -376,14 +363,12 @@ function trace_contour(x, y, z, h::Number, cells::Dict) starting_edge = E end - if !(0 < yi < yi_max && 0 < xi < xi_max) - push!(contours.lines, contour) - continue + if 0 < yi < yi_max && 0 < xi < xi_max + # Start trace in reverse direction + (xi, yi) = chase!(cells, reverse!(contour_arr), x, y, z, h, xi, yi, starting_edge, xi_max, yi_max) 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 From fcf3e0e4204fa36a03f6b639e5c8225ed81b3c59 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sun, 12 Apr 2020 21:18:49 -0400 Subject: [PATCH 02/10] refactor edge advancing and more cleanup --- src/Contour.jl | 60 ++++++++++++++++++++------------------------------ 1 file changed, 24 insertions(+), 36 deletions(-) diff --git a/src/Contour.jl b/src/Contour.jl index 5414b64..33515df 100644 --- a/src/Contour.jl +++ b/src/Contour.jl @@ -144,6 +144,9 @@ const WN, WS, WE = W|N, W|S, W|E const NWSE = NW | 0x10 # special (ambiguous case) const NESW = NE | 0x10 # special (ambiguous case) +# 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 @@ -176,7 +179,19 @@ function get_next_edge!(cells::Dict, xi, yi, entry_edge::UInt8) end end -function get_first_crossing(cell) +@inline function advance_edge(xi::T,yi::T,edge) where T + if edge == N + return xi, yi+one(T), S + elseif edge == S + return xi, yi-one(T), N + elseif edge == E + return xi+one(T), yi, W + elseif edge == W + return xi-one(T), yi, E + end +end + +@inline function get_first_crossing(cell) if cell == NWSE return NW elseif cell == NESW @@ -186,9 +201,6 @@ 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) - @inline function _get_case(z, h) case = z[1] > h ? 0x01 : 0x00 z[2] > h && (case |= 0x02) @@ -287,19 +299,8 @@ function chase!(cells, curve, x, y, z, h, xi_start, yi_start, entry_edge, xi_max push!(curve, interpolate(x, y, z, h, xi, yi, exit_edge)) - 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 = advance_edge(xi, yi, exit_edge) + !((xi, yi, entry_edge) != (xi_start, yi_start, loopback_edge) && 0 < yi < yi_max && 0 < xi < xi_max) && break end @@ -325,8 +326,7 @@ function trace_contour(x, y, z, h::Number, cells::Dict) contour_arr = VT[] # Pick initial box - (xi_0, yi_0), cell = first(cells) - (xi, yi) = (xi_0, yi_0) + (xi, yi), cell = first(cells) # Pick a starting edge crossing = get_first_crossing(cell) @@ -339,33 +339,21 @@ function trace_contour(x, y, z, h::Number, cells::Dict) end # Add the contour entry location for cell (xi_0,yi_0) - push!(contour_arr, interpolate(x, y, z, h, xi_0, yi_0, starting_edge)) + push!(contour_arr, interpolate(x, y, z, h, xi, yi, starting_edge)) # Start trace in forward direction - (xi_end, yi_end) = chase!(cells, contour_arr, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max) + xi_end, yi_end = chase!(cells, contour_arr, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max) - if xi_end == xi_0 && yi_end == yi_0 + if xi_end == xi && yi_end == yi 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 + xi, yi, starting_edge = advance_edge(xi, yi, starting_edge) if 0 < yi < yi_max && 0 < xi < xi_max # Start trace in reverse direction - (xi, yi) = chase!(cells, reverse!(contour_arr), x, y, z, h, xi, yi, starting_edge, xi_max, yi_max) + chase!(cells, reverse!(contour_arr), x, y, z, h, xi, yi, starting_edge, xi_max, yi_max) end push!(contours.lines, Curve2(contour_arr)) From bd373350de34829da51ab9739f0097db7f246221 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sun, 12 Apr 2020 23:49:34 -0400 Subject: [PATCH 03/10] reuse cell dictionary to avoid excessive GC. Reduces memory roughly 25% --- src/Contour.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/Contour.jl b/src/Contour.jl index 33515df..10ecf65 100644 --- a/src/Contour.jl +++ b/src/Contour.jl @@ -60,9 +60,9 @@ argument `level`. You'll usually call [`lines`](@ref) on the output of `contour`, and then iterate over the result. """ -function contour(x, y, z, level::Number) +function contour(x, y, z, level::Number, cells = Dict{Tuple{Int,Int},UInt8}()) # Todo: size checking on x,y,z - trace_contour(x, y, z, level, get_level_cells(z, level)) + trace_contour(x, y, z, level, get_level_cells(z, level, cells)) end """ @@ -76,7 +76,11 @@ 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) + # reuse lookup dictionary as it is run multiple times and empty after each run + cells = Dict{Tuple{Int,Int},UInt8}() + ContourCollection([contour(x, y, z, l, cells) for l in levels]) +end """ `contours(x,y,z,Nlevels::Integer)` Trace `Nlevels` contour levels at heights @@ -209,8 +213,8 @@ end case end -function get_level_cells(z, h::Number) - cells = Dict{Tuple{Int,Int},UInt8}() +function get_level_cells(z, h::Number, cells = Dict{Tuple{Int,Int},UInt8}()) + xi_max, yi_max = size(z) @inbounds for xi in 1:xi_max - 1 From 4915fced329847faa994e61be71d88d1cd13a286 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Mon, 13 Apr 2020 15:25:30 -0400 Subject: [PATCH 04/10] add length checks for test data --- test/verify_vertices.jl | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/test/verify_vertices.jl b/test/verify_vertices.jl index ae4649a..e61c3f5 100644 --- a/test/verify_vertices.jl +++ b/test/verify_vertices.jl @@ -198,6 +198,21 @@ contourlevels = Contour.contour(X, Y, Z, h) # Test contour location on a realistic dataset include("testdata.jl") -Contour.contours(x, y, z) - +cts = Contour.contours(x, y, z) +@test length(cts.contours) == 10 +cts_ct = (8, 8, 8, 8, 126, 7, 5, 5, 5, 4) +lines_ct = ([138, 220, 469, 138, 469, 208, 143, 143], + [220, 475, 140, 210, 146, 475, 140, 146], + [222, 481, 140, 214, 481, 145, 140, 145], + [228, 485, 214, 142, 485, 142, 146, 146], + [9, 7, 515, 501, 9, 9, 9, 7, 7, 9, 7, 15, 7, 18, 39, 9, 7, 7, 9, 9, 7, 7, 7, 7, 9, 7, 18, 9, 7, 7, 9, 9, 35, 9, 9, 9, 7, 7, 9, 35, 9, 7, 9, 7, 9, 9, 7, 7, 7, 9, 7, 7, 7, 7, 9, 7, 9, 7, 9, 7, 7, 9, 9, 9, 7, 9, 7, 7, 9, 7, 7, 9, 7, 9, 7, 7, 7, 7, 7, 9, 7, 7, 7, 9, 7, 7, 7, 7, 9, 7, 7, 7, 7, 7, 7, 7, 9, 7, 9, 9, 7, 7, 7, 7, 7, 7, 9, 7, 7, 7, 7, 7, 7, 7, 9, 7, 9, 7, 7, 7, 7, 9, 7, 7, 7, 7], + [17, 34, 14, 34, 14, 5, 5], + [15, 29, 29, 12, 12], + [13, 26, 10, 10, 26], + [11, 23, 21, 11, 11], + [18, 7, 19, 7]) +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 end From e7df228cd7817a87ae4b0a1f18e0aa8ebff99307 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Mon, 13 Apr 2020 15:49:31 -0400 Subject: [PATCH 05/10] reorg interpolation and propogate vector type, closes #41 --- Project.toml | 3 ++- src/Contour.jl | 55 +++++-------------------------------- src/interpolate.jl | 60 +++++++++++++++++++++++++++++++++++++++++ test/verify_vertices.jl | 11 ++++++++ 4 files changed, 80 insertions(+), 49 deletions(-) create mode 100644 src/interpolate.jl diff --git a/Project.toml b/Project.toml index 5aa343c..5b8ccee 100644 --- a/Project.toml +++ b/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"] diff --git a/src/Contour.jl b/src/Contour.jl index 10ecf65..2aa07d1 100644 --- a/src/Contour.jl +++ b/src/Contour.jl @@ -2,6 +2,8 @@ module Contour using StaticArrays +include("interpolate.jl") + export ContourLevel, Curve2, @@ -242,53 +244,10 @@ function get_level_cells(z, h::Number, cells = Dict{Tuple{Int,Int},UInt8}()) return cells 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 SVector{2,T}(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 SVector{2,T}(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) +function chase!(cells, curve, x, y, z, h, xi_start, yi_start, entry_edge, xi_max, yi_max, ::Type{VT}) where VT xi, yi = xi_start, yi_start @@ -301,7 +260,7 @@ function chase!(cells, curve, x, y, z, h, xi_start, yi_start, entry_edge, xi_max @inbounds while true exit_edge = get_next_edge!(cells, xi, yi, entry_edge) - push!(curve, interpolate(x, y, z, h, xi, yi, exit_edge)) + push!(curve, interpolate(x, y, z, h, xi, yi, exit_edge, VT)) xi, yi, entry_edge = advance_edge(xi, yi, exit_edge) @@ -343,10 +302,10 @@ function trace_contour(x, y, z, h::Number, cells::Dict) end # Add the contour entry location for cell (xi_0,yi_0) - push!(contour_arr, interpolate(x, y, z, h, xi, yi, starting_edge)) + push!(contour_arr, interpolate(x, y, z, h, xi, yi, starting_edge, VT)) # Start trace in forward direction - xi_end, yi_end = chase!(cells, contour_arr, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max) + xi_end, yi_end = chase!(cells, contour_arr, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max, VT) if xi_end == xi && yi_end == yi push!(contours.lines, Curve2(contour_arr)) @@ -357,7 +316,7 @@ function trace_contour(x, y, z, h::Number, cells::Dict) if 0 < yi < yi_max && 0 < xi < xi_max # Start trace in reverse direction - chase!(cells, reverse!(contour_arr), x, y, z, h, xi, yi, starting_edge, xi_max, yi_max) + chase!(cells, reverse!(contour_arr), x, y, z, h, xi, yi, starting_edge, xi_max, yi_max, VT) end push!(contours.lines, Curve2(contour_arr)) diff --git a/src/interpolate.jl b/src/interpolate.jl new file mode 100644 index 0000000..091f786 --- /dev/null +++ b/src/interpolate.jl @@ -0,0 +1,60 @@ +# 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, xi, yi, edge::UInt8, ::Type{VT}) where {VT} + @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, xi, yi, edge::UInt8, ::Type{VT}) where {VT} + @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, xi, yi, edge::UInt8, ::Type{VT}) where {VT} + 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 diff --git a/test/verify_vertices.jl b/test/verify_vertices.jl index e61c3f5..ba5f91e 100644 --- a/test/verify_vertices.jl +++ b/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 From a12a018706c2ea42663f97b4d6aa7aa8e7d93065 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Mon, 13 Apr 2020 21:03:00 -0400 Subject: [PATCH 06/10] cleanup get next edge to use pop! --- src/Contour.jl | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/src/Contour.jl b/src/Contour.jl index 2aa07d1..5b0b821 100644 --- a/src/Contour.jl +++ b/src/Contour.jl @@ -161,28 +161,25 @@ const edge_LUT = (SW, SE, EW, NE, 0x0, NS, NW, NW, NS, 0x0, NE, EW, SE, SW) # have two crossings. function get_next_edge!(cells::Dict, xi, yi, entry_edge::UInt8) key = (xi,yi) - cell = cells[key] + 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 @inline function advance_edge(xi::T,yi::T,edge) where T From 4d067a8902e2cbbf9b1abc62b5b673c74ebff83f Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Mon, 13 Apr 2020 02:27:10 -0400 Subject: [PATCH 07/10] wip --- src/Contour.jl | 61 +++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 60 insertions(+), 1 deletion(-) diff --git a/src/Contour.jl b/src/Contour.jl index 5b0b821..2bc3402 100644 --- a/src/Contour.jl +++ b/src/Contour.jl @@ -152,6 +152,8 @@ const NESW = NE | 0x10 # special (ambiguous case) # 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 edge_pairs = ((S,W),(S,E),(E,W),(N,E),(0x0,0x0),(N,S),(N,W), + (W,N),(S,N),(0x0,0x0),(E,N),(W,E),(E,S),(W,S)) # The way a contour crossing goes through a cell is labeled # by combining compass directions (e.g. a NW crossing connects @@ -189,7 +191,7 @@ end return xi, yi-one(T), N elseif edge == E return xi+one(T), yi, W - elseif edge == W + else # W return xi-one(T), yi, E end end @@ -241,6 +243,63 @@ function get_level_cells(z, h::Number, cells = Dict{Tuple{Int,Int},UInt8}()) return cells end +struct MSEdges + +end + +function contours(x,y,z, levels::Integer, T::MSEdges) + contours(x,y,z,contourlevels(z, levels), T) +end + +function contours(x, y, z, levels, ::MSEdges) + + xi_max, yi_max = size(z) + + VT = SVector{2,eltype(z)} + + edge_list = [Vector{Pair{VT,VT}}() for _ in levels] + + @inbounds for xi in 1:xi_max - 1 + for yi in 1:yi_max - 1 + elts = (z[xi, yi], z[xi + 1, yi], z[xi + 1, yi + 1], z[xi, yi + 1]) + for i in eachindex(levels) + h = levels[i] + case = _get_case(elts, h) + + # Contour does not go through these cells + if iszero(case) || case == 0x0f + continue + end + + # 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 #? NWSE : NESW + push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N) => interpolate(x,y,z,h,xi,yi,W), + interpolate(x,y,z,h,xi,yi,S) => interpolate(x,y,z,h,xi,yi,E)) + else + push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N) => interpolate(x,y,z,h,xi,yi,E), + interpolate(x,y,z,h,xi,yi,S) => interpolate(x,y,z,h,xi,yi,W)) + end + elseif case == 0x0a + if 0.25*sum(elts) >= h #? NESW : NWSE + push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N) => interpolate(x,y,z,h,xi,yi,E), + interpolate(x,y,z,h,xi,yi,S) => interpolate(x,y,z,h,xi,yi,W)) + else + push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N) => interpolate(x,y,z,h,xi,yi,W), + interpolate(x,y,z,h,xi,yi,S) => interpolate(x,y,z,h,xi,yi,E)) + end + else + ep = edge_pairs[case] + push!(edge_list[i], interpolate(x,y,z,h,xi,yi,ep[1]) => interpolate(x,y,z,h,xi,yi,ep[2])) + end + end + end + end + + return edge_list +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. From cf8740c54e1a48560f5f28dd1b79940a3b16c497 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Mon, 13 Apr 2020 23:26:53 -0400 Subject: [PATCH 08/10] wip2 --- src/Contour.jl | 30 +++++++++++++++++++++--------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/src/Contour.jl b/src/Contour.jl index 2bc3402..b11bc70 100644 --- a/src/Contour.jl +++ b/src/Contour.jl @@ -259,9 +259,21 @@ function contours(x, y, z, levels, ::MSEdges) edge_list = [Vector{Pair{VT,VT}}() for _ in levels] + for e in edge_list + sizehint!(e, (2*xi_max+2*yi_max)) + end + @inbounds for xi in 1:xi_max - 1 for yi in 1:yi_max - 1 + elts = (z[xi, yi], z[xi + 1, yi], z[xi + 1, yi + 1], z[xi, yi + 1]) + + # check upper and lower bounds first to see if we can skip + lb_h = first(levels); lb = _get_case(elts, lb_h) + ub_h = last(levels); ub = _get_case(elts, ub_h) + iszero(lb) && iszero(ub) && continue + lb == 0x0f && ub == 0x0f && continue + for i in eachindex(levels) h = levels[i] case = _get_case(elts, h) @@ -275,23 +287,23 @@ function contours(x, y, z, levels, ::MSEdges) # a bilinear interpolation of the cell-center value. if case == 0x05 if 0.25*sum(elts) >= h #? NWSE : NESW - push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N) => interpolate(x,y,z,h,xi,yi,W), - interpolate(x,y,z,h,xi,yi,S) => interpolate(x,y,z,h,xi,yi,E)) + push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,W,VT), + interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,E,VT)) else - push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N) => interpolate(x,y,z,h,xi,yi,E), - interpolate(x,y,z,h,xi,yi,S) => interpolate(x,y,z,h,xi,yi,W)) + push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,E,VT), + interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,W,VT)) end elseif case == 0x0a if 0.25*sum(elts) >= h #? NESW : NWSE - push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N) => interpolate(x,y,z,h,xi,yi,E), - interpolate(x,y,z,h,xi,yi,S) => interpolate(x,y,z,h,xi,yi,W)) + push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,E,VT), + interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,W,VT)) else - push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N) => interpolate(x,y,z,h,xi,yi,W), - interpolate(x,y,z,h,xi,yi,S) => interpolate(x,y,z,h,xi,yi,E)) + push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,W,VT), + interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,E,VT)) end else ep = edge_pairs[case] - push!(edge_list[i], interpolate(x,y,z,h,xi,yi,ep[1]) => interpolate(x,y,z,h,xi,yi,ep[2])) + push!(edge_list[i], Pair{VT,VT}(interpolate(x,y,z,h,xi,yi,ep[1],VT),interpolate(x,y,z,h,xi,yi,ep[2],VT)) ) end end end From ff75d4af7f95a2b29913967cac6d48a2a09a4f13 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Tue, 14 Apr 2020 00:48:43 -0400 Subject: [PATCH 09/10] better mem --- bench/perf.jl | 3 ++- src/Contour.jl | 68 +++++++++++++++++++++----------------------------- 2 files changed, 31 insertions(+), 40 deletions(-) diff --git a/bench/perf.jl b/bench/perf.jl index 39ca5f9..7d187aa 100644 --- a/bench/perf.jl +++ b/bench/perf.jl @@ -13,7 +13,8 @@ include(joinpath(@__DIR__,"../test/testdata.jl")) v = Contour.contours(x, y, z) @show typeof(v) -suite["contour"]["testdata"] = @benchmarkable Contour.contours($x,$y,$z) +suite["contour"]["chase"] = @benchmarkable Contour.contours($x,$y,$z, 10) +suite["contour"]["edge_list"] = @benchmarkable Contour.contours($x,$y,$z, 10, Contour.MSEdges()) results = run(suite) diff --git a/src/Contour.jl b/src/Contour.jl index b11bc70..266d187 100644 --- a/src/Contour.jl +++ b/src/Contour.jl @@ -251,60 +251,50 @@ function contours(x,y,z, levels::Integer, T::MSEdges) contours(x,y,z,contourlevels(z, levels), T) end -function contours(x, y, z, levels, ::MSEdges) +function contours(x,y,z, levels::Union{AbstractArray,AbstractRange}, T::MSEdges) + [contour(x,y,z,h,T) for h in levels] +end + +function contour(x, y, z, h::Number, ::MSEdges) xi_max, yi_max = size(z) VT = SVector{2,eltype(z)} - edge_list = [Vector{Pair{VT,VT}}() for _ in levels] - - for e in edge_list - sizehint!(e, (2*xi_max+2*yi_max)) - end + edge_list = Pair{VT,VT}[] @inbounds for xi in 1:xi_max - 1 for yi in 1:yi_max - 1 elts = (z[xi, yi], z[xi + 1, yi], z[xi + 1, yi + 1], z[xi, yi + 1]) + case = _get_case(elts, h) - # check upper and lower bounds first to see if we can skip - lb_h = first(levels); lb = _get_case(elts, lb_h) - ub_h = last(levels); ub = _get_case(elts, ub_h) - iszero(lb) && iszero(ub) && continue - lb == 0x0f && ub == 0x0f && continue - - for i in eachindex(levels) - h = levels[i] - case = _get_case(elts, h) + # Contour does not go through these cells + if iszero(case) || case == 0x0f + continue + end - # Contour does not go through these cells - if iszero(case) || case == 0x0f - continue + # 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 #? NWSE : NESW + push!(edge_list, interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,W,VT), + interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,E,VT)) + else + push!(edge_list, interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,E,VT), + interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,W,VT)) end - - # 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 #? NWSE : NESW - push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,W,VT), - interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,E,VT)) - else - push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,E,VT), - interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,W,VT)) - end - elseif case == 0x0a - if 0.25*sum(elts) >= h #? NESW : NWSE - push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,E,VT), - interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,W,VT)) - else - push!(edge_list[i], interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,W,VT), - interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,E,VT)) - end + elseif case == 0x0a + if 0.25*sum(elts) >= h #? NESW : NWSE + push!(edge_list, interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,E,VT), + interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,W,VT)) else - ep = edge_pairs[case] - push!(edge_list[i], Pair{VT,VT}(interpolate(x,y,z,h,xi,yi,ep[1],VT),interpolate(x,y,z,h,xi,yi,ep[2],VT)) ) + push!(edge_list, interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,W,VT), + interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,E,VT)) end + else + ep = edge_pairs[case] + push!(edge_list, Pair{VT,VT}(interpolate(x,y,z,h,xi,yi,ep[1],VT),interpolate(x,y,z,h,xi,yi,ep[2],VT)) ) end end end From 211adb77c8b0e23ba40e2c2ecf665d35dfd5e7cc Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Tue, 14 Apr 2020 14:47:15 -0400 Subject: [PATCH 10/10] non exhaustive interpolation cache --- src/Contour.jl | 32 ++++++++++++++++++++++---------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/src/Contour.jl b/src/Contour.jl index 266d187..73b6342 100644 --- a/src/Contour.jl +++ b/src/Contour.jl @@ -152,7 +152,7 @@ const NESW = NE | 0x10 # special (ambiguous case) # 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 edge_pairs = ((S,W),(S,E),(E,W),(N,E),(0x0,0x0),(N,S),(N,W), +const edge_pairs = ((W,S),(E,S),(E,W),(E,N),(0x0,0x0),(N,S),(W,N), (W,N),(S,N),(0x0,0x0),(E,N),(W,E),(E,S),(W,S)) # The way a contour crossing goes through a cell is labeled @@ -263,6 +263,9 @@ function contour(x, y, z, h::Number, ::MSEdges) edge_list = Pair{VT,VT}[] + # save vertices + prior_vts = Vector{VT}(undef, yi_max-1) + @inbounds for xi in 1:xi_max - 1 for yi in 1:yi_max - 1 @@ -278,23 +281,23 @@ function contour(x, y, z, h::Number, ::MSEdges) # a bilinear interpolation of the cell-center value. if case == 0x05 if 0.25*sum(elts) >= h #? NWSE : NESW - push!(edge_list, interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,W,VT), - interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,E,VT)) + add_edge!(x,y,z,h,xi,yi,(N,W),edge_list,prior_vts,VT) + add_edge!(x,y,z,h,xi,yi,(S,E),edge_list,prior_vts,VT) else - push!(edge_list, interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,E,VT), - interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,W,VT)) + add_edge!(x,y,z,h,xi,yi,(N,E),edge_list,prior_vts,VT) + add_edge!(x,y,z,h,xi,yi,(S,W),edge_list,prior_vts,VT) end elseif case == 0x0a if 0.25*sum(elts) >= h #? NESW : NWSE - push!(edge_list, interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,E,VT), - interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,W,VT)) + add_edge!(x,y,z,h,xi,yi,(N,E),edge_list,prior_vts,VT) + add_edge!(x,y,z,h,xi,yi,(S,W),edge_list,prior_vts,VT) else - push!(edge_list, interpolate(x,y,z,h,xi,yi,N,VT) => interpolate(x,y,z,h,xi,yi,W,VT), - interpolate(x,y,z,h,xi,yi,S,VT) => interpolate(x,y,z,h,xi,yi,E,VT)) + add_edge!(x,y,z,h,xi,yi,(N,W),edge_list,prior_vts,VT) + add_edge!(x,y,z,h,xi,yi,(S,E),edge_list,prior_vts,VT) end else ep = edge_pairs[case] - push!(edge_list, Pair{VT,VT}(interpolate(x,y,z,h,xi,yi,ep[1],VT),interpolate(x,y,z,h,xi,yi,ep[2],VT)) ) + add_edge!(x,y,z,h,xi,yi,ep,edge_list,prior_vts,VT) end end end @@ -302,6 +305,15 @@ function contour(x, y, z, h::Number, ::MSEdges) return edge_list end +function add_edge!(x,y,z,h,xi,yi,ep,edge_list,prior_vts,::Type{VT}) where VT + if ep[1] == W + fv = interpolate(x,y,z,h,xi,yi,ep[1],VT) + prior_vts[yi] = fv + push!(edge_list, fv => interpolate(x,y,z,h,xi,yi,ep[2],VT)) + elseif ep[1] == E && xi != 1 + push!(edge_list, prior_vts[yi] => interpolate(x,y,z,h,xi,yi,ep[2],VT)) + end +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.