RFC: Interpolate Highways

src/crop.jl

@@ -69,9 +69,7 @@ end
 function crop!(nodes::Dict, bounds::Bounds, highways::Dict{Int,Highway})
     missing_nodes = Int[]
 
-    for key in keys(highways)
-        highway = highways[key]
-
+    for (key, highway) in highways
         valid = falses(length(highway.nodes))
         #println(highway.nodes)
         #for n = 1:length(highway.nodes)
@@ -92,6 +90,7 @@ function crop!(nodes::Dict, bounds::Bounds, highways::Dict{Int,Highway})
         end
 
         nodes_in_bounds = sum(valid)
+
         if nodes_in_bounds == 0
             delete!(highways,key)   # Remove highway from list
         elseif nodes_in_bounds < length(valid)
@@ -144,26 +143,32 @@ function inBounds(loc::LLA, bounds::Bounds)
     lat = loc.lat
     lon = loc.lon
 
-    if lat < bounds.min_lat || lat > bounds.max_lat
-        return false
-    elseif lon < bounds.min_lon || lon > bounds.max_lon
-        return false
-    end
-
-    return true
+    bounds.min_lat <= lat <= bounds.max_lat &&
+    bounds.min_lon <= lon <= bounds.max_lon
 end
 
 function inBounds(loc::ENU, bounds::Bounds)
     north = loc.north
     east = loc.east
 
-    if north < bounds.min_lat || north > bounds.max_lat
-        return false
-    elseif east < bounds.min_lon || east > bounds.max_lon
-        return false
-    end
+    bounds.min_lat <= north <= bounds.max_lat &&
+    bounds.min_lon <= east <= bounds.max_lon
+end
+
+function onBounds(loc::LLA, bounds::Bounds)
+    lat = loc.lat
+    lon = loc.lon
 
-    return true
+    lat == bounds.min_lat || lat == bounds.max_lat ||
+    lon == bounds.min_lon || lon == bounds.max_lon
+end
+
+function onBounds(loc::ENU, bounds::Bounds)
+    north = loc.north
+    east = loc.east
+
+    north == bounds.min_lat || north == bounds.max_lat ||
+    east == bounds.min_lon || east == bounds.max_lon
 end
 
 ### Remove specified items from an array ###
@@ -179,138 +184,62 @@ function cropList!(list::Array, crop_list::BitArray{1})
     return nothing
 end
 
-### Crop highway to fit within bounds, interpolating to place ###
-### new nodes on the boundary as necessary.                   ###
-function cropHighway!(nodes::Dict, bounds::Bounds, highway::Highway, valid::BitArray{1})
-    inside = find(valid)
-    first_inside = inside[1]
-    last_inside = inside[end]
-
-    # Remove bad nodes at end of highway node list
-    if last_inside+1 < length(highway.nodes)
-        for k = (last_inside+2):length(highway.nodes)
-            pop!(highway.nodes)
-            pop!(valid)
-        end
-    end
+function boundaryPoint{T}(p1::T, p2::T, bounds::Bounds)
+    x1, y1 = getX(p1), getY(p1)
+    x2, y2 = getX(p2), getY(p2)
 
-    # Remove bad nodes at start of highway node list
-    if first_inside > 2
-        ind = first_inside - 2
-        for k = 1:(first_inside-2)
-            splice!(highway.nodes,ind)
-            splice!(valid,ind)
-            ind -= 1
-        end
-    end
+    x, y = x1, y1
 
-    interpolate_start = !valid[1]
-    interpolate_end = !valid[end]
-
-    if interpolate_end
-        last_inside = find(valid)[end]
-        const node0 = highway.nodes[last_inside]
-        const x0 = getX( nodes[node0] )
-        const y0 = getY( nodes[node0] )
-        node1 = highway.nodes[last_inside+1]
-        x1 = getX( nodes[node1] )
-        y1 = getY( nodes[node1] )
-
-        if x1 < bounds.min_lon || x1 > bounds.max_lon
-            if x1 < bounds.min_lon
-                x = bounds.min_lon
-            else
-                x = bounds.max_lon
-            end
-            y = y0 + (y1 - y0) * (x - x0) / (x1 - x0)
+    # checks assume inBounds(p1) != inBounds(p2)
+    if x1 < bounds.min_lon < x2 || x1 > bounds.min_lon > x2
+        x = bounds.min_lon
+        y = y1 + (y2 - y1) * (bounds.min_lon - x1) / (x2 - x1)
+    elseif x1 < bounds.max_lon < x2 || x1 > bounds.max_lon > x2
+        x = bounds.max_lon
+        y = y1 + (y2 - y1) * (bounds.max_lon - x1) / (x2 - x1)
+    end
 
-            # Add a new node to nodes list
-            if typeof(nodes[node0]) == LLA
-                new_id = addNewNode(nodes,LLA(y,x))
-            else # ENU
-                new_id = addNewNode(nodes,ENU(x,y))
-            end
-            highway.nodes[last_inside+1] = new_id
-            valid[last_inside+1] = inBounds(nodes[new_id],bounds)
-        end
+    p3 = T == LLA ? T(y, x) : T(x, y)
+    inBounds(p3, bounds) && return p3
 
-        if !valid[last_inside+1]
-            node1 = highway.nodes[last_inside+1]
-            x1 = getX( nodes[node1] )
-            y1 = getY( nodes[node1] )
-
-            if y1 < bounds.min_lat || y1 > bounds.max_lat
-                if y1 < bounds.min_lat
-                    y = bounds.min_lat
-                else
-                    y = bounds.max_lat
-                end
-                x = x0 + (x1-x0) * (y - y0) / (y1 - y0)
-
-                # Add a new node to nodes list
-                if typeof(nodes[node0]) == LLA
-                    new_id = addNewNode(nodes,LLA(y,x))
-                else # ENU
-                    new_id = addNewNode(nodes,ENU(x,y))
-                end
-                highway.nodes[last_inside+1] = new_id
-                valid[last_inside+1] = inBounds(nodes[new_id],bounds)
-            end
-        end
+    if y1 < bounds.min_lat < y2 || y1 > bounds.min_lat > y2
+        x = x1 + (x2 - x1) * (bounds.min_lat - y1) / (y2 - y1)
+        y = bounds.min_lat
+    elseif y1 < bounds.max_lat < y2 || y1 > bounds.max_lat > y2
+        x = x1 + (x2 - x1) * (bounds.max_lat - y1) / (y2 - y1)
+        y = bounds.max_lat
     end
 
-    if interpolate_start
-        first_inside = find(valid)[1]
-        const node0 = highway.nodes[first_inside]
-        const x0 = getX( nodes[node0] )
-        const y0 = getY( nodes[node0] )
-        node1 = highway.nodes[first_inside-1]
-        x1 = getX( nodes[node1] )
-        y1 = getY( nodes[node1] )
-
-        if x1 < bounds.min_lon || x1 > bounds.max_lon
-            if x1 < bounds.min_lon
-                x = bounds.min_lon
-            else
-                x = bounds.max_lon
-            end
-            y = y0 + (y1 - y0) * (x - x0) / (x1 - x0);
+    p3 = T == LLA ? T(y, x) : T(x, y)
+    inBounds(p3, bounds) && return p3
 
-            # Add a new node to nodes list
-            if typeof(nodes[node0]) == LLA
-                new_id = addNewNode(nodes,LLA(y,x))
-            else # ENU
-                new_id = addNewNode(nodes,ENU(x,y))
-            end
-            highway.nodes[first_inside-1] = new_id
-            valid[first_inside-1] = inBounds(nodes[new_id],bounds)
-        end
+    error("Failed to find boundary point.")
+end
 
-        if !valid[first_inside-1]
-            node1 = highway.nodes[first_inside-1]
-            x1 = getX( nodes[node1] )
-            y1 = getY( nodes[node1] )
-
-            if y1 < bounds.min_lat || y1 > bounds.max_lat
-                if y1 < bounds.min_lat
-                    y = bounds.min_lat
-                else
-                    y = bounds.max_lat
-                end
-                x = x0 + (x1-x0) * (y - y0) / (y1 - y0)
-
-                # Add a new node to nodes list
-                if typeof(nodes[node0]) == LLA
-                    new_id = addNewNode(nodes,LLA(y,x))
-                else # ENU
-                    new_id = addNewNode(nodes,ENU(x,y))
-                end
-                highway.nodes[first_inside-1] = new_id
-                valid[first_inside-1] = inBounds(nodes[new_id],bounds)
+function cropHighway!(nodes::Dict, bounds::Bounds, highway::Highway, valids::BitArray{1})
+    prev_id, prev_valid = highway.nodes[1], valids[1]
+    ni = 1
+    for i in 1:length(valids)
+        id, valid = highway.nodes[ni], valids[i]
+
+        if !valid
+            deleteat!(highway.nodes, ni)
+            ni -= 1
+        end
+        if valid != prev_valid
+            prev_node, node = nodes[prev_id], nodes[id]
+            if !(onBounds(prev_node, bounds) || onBounds(node, bounds))
+                new_node = boundaryPoint(prev_node, node, bounds)
+                new_id = addNewNode(nodes, new_node)
+                insert!(highway.nodes, ni + !valid, new_id)
+                ni += 1
             end
         end
+
+        ni += 1
+
+        prev_id, prev_valid = id, valid
     end
 
     return nothing
 end
-

test/classes.jl

@@ -1,4 +1,5 @@
 # Test extracting way, building, and feature classes from OSM data
+module TestClasses
 
 using OpenStreetMap
 using Base.Test
@@ -25,3 +26,5 @@ feat_classes = classify( feats )
 for key in keys(feat_classes)
     @test feat_classes[key] == 1
 end
+
+end # module TestClasses

test/coordinates.jl

@@ -1,4 +1,5 @@
 # Test Coordinate Transforms
+module TestCoordinates
 
 using OpenStreetMap
 using Base.Test
@@ -42,3 +43,5 @@ bounds_enu = lla2enu(bounds)
 @test_approx_eq bounds_enu.max_lat 249.9353534127322
 @test_approx_eq bounds_enu.min_lon -247.1091823451915
 @test_approx_eq bounds_enu.max_lon 247.1268196141778
+
+end # module TestCoordinates

test/crop_map.jl

@@ -1,16 +1,59 @@
 # Test cropping OSM data
+module TestCropMap
 
 using OpenStreetMap
 using Base.Test
 
+import OpenStreetMap: Bounds, inBounds, LLA, getX, getY
+
 MAP_FILENAME = "tech_square.osm"
 
-nodes, hwys, builds, feats = getOSMData( MAP_FILENAME, nodes=true, highways=true, buildings=true, features=true)
+nodes, hwys, builds, feats = getOSMData(MAP_FILENAME, nodes=true, highways=true, buildings=true, features=true)
 
-bounds = OpenStreetMap.Bounds(42.3637,42.3655,-71.0919,-71.0893)
+bounds = Bounds(42.3637, 42.3655, -71.0919, -71.0893)
 cropMap!(nodes, bounds, highways=hwys, buildings=builds, features=feats, delete_nodes=true)
 
 @test length(nodes) == 198
 @test length(hwys) == 16
 @test length(builds) == 1
 @test length(feats) == 0
+
+# cropHighways removes out-of-bounds points and interpolates points on boundary
+let bounds = Bounds(0, 1, 0, 1)
+    coords = [
+        LLA(-0.1, -0.02),
+        LLA(0.1, 0.1),
+        LLA(0.2, 0.0),
+        LLA(0.3, -0.1),
+        LLA(0.4, 0.1),
+        LLA(1.0, 0.3),
+        LLA(2.0, 0.4)
+    ]
+
+    hwy = first(values(hwys))
+    hwy.nodes = [1:length(coords)]
+    highways = [1 => hwy]
+
+    nodes = Dict(hwy.nodes, coords)
+    cropMap!(nodes, bounds, highways=highways)
+
+    hwy_nodes = highways[1].nodes
+
+    # same changes made to nodes and highways
+    @test isempty(symdiff(keys(nodes), hwy_nodes))
+
+    # everything ends up in bounds
+    @test all(map(x -> inBounds(x, bounds), values(nodes)))
+
+    # proper interpolation
+    @test getY(nodes[hwy_nodes[1]]) == 0.0
+    @test_approx_eq getX(nodes[hwy_nodes[1]]) 0.04
+
+    @test_approx_eq getY(nodes[hwy_nodes[4]]) 0.35
+    @test getX(nodes[hwy_nodes[4]]) == 0.0
+
+    @test getY(nodes[hwy_nodes[end]]) == 1.0
+    @test_approx_eq getX(nodes[hwy_nodes[end]]) 0.3
+end
+
+end # module TestCropMap

test/plots.jl

@@ -1,4 +1,5 @@
 # Test plotting
+module TestPlots
 
 using OpenStreetMap
 using Base.Test
@@ -19,3 +20,5 @@ feat_classes = classify( feats )
 fignum = plotMap(nodes, highways=hwys, buildings=builds, features=feats, bounds=bounds, width=500, feature_classes=feat_classes, building_classes=bldg_classes, cycleways=cycles, walkways=peds, roadways=roads)
 
 @test fignum == 1
+
+end # module TestPlots

test/read_data.jl

@@ -1,4 +1,5 @@
 # Test reading in OSM data
+module TestReadData
 
 using OpenStreetMap
 using Base.Test
@@ -24,3 +25,5 @@ bounds = getBounds(map)
 @test_approx_eq bounds.max_lat 42.3659000
 @test_approx_eq bounds.min_lon -71.0939000
 @test_approx_eq bounds.max_lon -71.0891000
+
+end # module TestReadData

test/routes.jl

@@ -1,7 +1,9 @@
 # Test route planning
+module TestRoutes
 
 using OpenStreetMap
 using Base.Test
+using Graphs
 
 MAP_FILENAME = "tech_square.osm"
 
@@ -14,6 +16,7 @@ cropMap!(nodes, bounds, highways=hwys, buildings=builds, features=feats, delete_
 nodesENU = lla2enu( nodes, OpenStreetMap.centerBounds(bounds) )
 
 # Form transportation network
+roads = roadways( hwys )
 network = createGraph( nodesENU, hwys, roads, Set(1:8...) )
 
 @test Graphs.num_vertices(network.g) == 155
@@ -58,3 +61,5 @@ shortest_edges = routeEdges( network, shortest_route )
 
 fastest_edges = routeEdges( network, fastest_route )
 @test length(fastest_edges) == 21
+
+end # module TestRoutes

test/simcity.jl

@@ -1,4 +1,5 @@
 # Test city simulation
+module TestSimCity
 
 using OpenStreetMap
 using Base.Test
@@ -31,3 +32,5 @@ nodes, highways, roads = simCityGrid(roads_north,roads_east)
 @test roads[3] == roads_north[3]
 @test roads[5] == roads_east[2]
 @test roads[7] == roads_east[4]
+
+end # module TestSimCity