ms_rectify_maps with type = "lines"

[mpadge]

This is much tricker than polygons or points, but likely to also be useful, so should be done ...

[mpadge]

@mdsumner I've just wasted couple of hours trying to crack this to no avail. Any insights on how to either:

1. Convert a binary raster into some kind of vector "midline"? or ...
2. Get the midline of an elongated polygon? (I can easily trace polygons around the perimeter of the binary raster).

The only approaches I've found use voronoi tesselations, remove all edge triangles, and then somehow "trace" along a central ridge - but I can't figure out how to get something like that to work. I also played with decido a bit here, but couldn't get anywhere with that either ... Any help greatly appreciated!

Sample traced polygon here:

[mdsumner]

Post the data vectors!

This is a common ask that no one ever solves, it's called "skeletonization" - but I've had some wins doing it before, the problem is having some kind of idea of the start and end

[mdsumner]

If you try a rasterized version then isoband::isolines might do it

[mdsumner]

Nope.

[mpadge]

Skeletonization as I understand is this:

I've coded up this algorithm to do that job, but those are still pixels that have no idea of order ...

[mpadge]

I've got an idea for how to identify start/end points. As you state, once those have been identified, the rest is relatively easy.

[mdsumner]

This is the result from sfdct::ct_triangulate on my own snaky polygon:

Those internal edges don't occur in the input edges, and their mid-points might be enough?

Still, order is always tough!

[mdsumner]

Oh, it was something about being able to order the internal triangles, so when they are dense they can be classified by distance to the edge:

lighter colour is further from the original boundary, then maybe traverse the triangles along that peak, or even rasterize it

[mdsumner]

Getting pretty heavy on triangles!

I dunno ...

[mpadge]

That looks good! Here's some code to extract one line in both raster and polygon form:

devtools::load_all (".", export_all = TRUE)
# code to extract the drawn-over stuff:
original <- system.file ("extdata", "omaha.png", package = "mapscanner")
f_mod <- system.file ("extdata", "omaha-lines.png", package = "mapscanner")
map <- png::readPNG (original)
map_scanned <- png::readPNG (f_mod)
res <- m_niftyreg (map_scanned, map)
channel <- extract_channel (res)
image (channel)

# code to separate the objects:
img <- matrix (as.logical (channel), nrow = nrow (channel))
cmat <- rcpp_components (img)
comps <- seq (1:max (cmat))

imgs <- list ()
for (ci in comps)
{
    cmat_i <- cmat
    cmat_i [cmat_i != ci] <- 0
    i1 <- get_shape_index (cmat_i, x = TRUE)
    i2 <- get_shape_index (cmat_i, x = FALSE)

    imgs [[length (imgs) + 1]] <- cmat_i [i1, i2]
}

# then code to analyse one of those:
i <- 2
im <- array (FALSE, dim = dim (imgs [[i]]))
im [which (imgs [[i]] > 0)] <- TRUE
image (im)
b <- rcpp_boundary (im)
#b <- rcpp_edge_thin (im) # get the pixel skeleton

That b object should give the above. Any code you could jig together to get a centre line would be awesome!!!!

[mdsumner]

Ok, this is horrific but looks promising.

#a <- do.call(cbind, locator())
library(sf)
library(dplyr)
a <- structure(c(116.931396841307, 152.008842201758, 173.445058810923, 
            174.614306989605, 152.008842201758, 115.762148662625, 67.4332239437809, 
            117.710895627095, 42.879012191465, -25.3271315649682, -1.55241859844007, 
            56.9099903356455, 54.1817445853882, 52.6227470138126, 67.8229733366748, 
            1.6816106511235, 30.5230657252724, 74.5647471222835, 125.232168198491, 
            227.346509136694, 241.377487280874, 237.479993351935, 170.443097774184, 
            172.391844738653, 175.899589274699, 124.062920019809, 97.5599613030239, 
            27.4050705821212, -2.99538206360334, 0.902111865335695), .Dim = c(15L, 
                                                                              2L), .Dimnames = list(NULL, c("x", "y")))

## we segmentize to get a dense Voronoi
sf <- sf::st_sfc(sf::st_segmentize(sf::st_buffer(sf::st_linestring(a), 3), 1)); plot(sf)
x <- silicate::SC(sf); max_area <- NULL

## converte SC edges to RTriangle (Delaunay + Voronoi)
edge_RTriangle <- anglr:::edge_RTriangle
edge_per_object_lt <- x$object["object_"] %>% dplyr::inner_join(x$object_link_edge,
"object_") %>% dplyr::inner_join(x$object_link_edge,
"object_") %>% dplyr::group_by(.data$object_) %>% dplyr::tally(dplyr::n()) %>%
dplyr::filter(n < 3)
if (nrow(edge_per_object_lt) > 0) {
 message("dropping untriangulatable objects")
 x <- dplyr::filter(x, !.data$object_ %in% edge_per_object_lt$object_)
 drop <- x$vertex$vertex_ %in% c(x$edge$.vx0, x$edge$.vx1)
 if (any(drop)) x$vertex <- x$vertex[drop, ]
}
objs <- vector("list", nrow(x$object))
for (i in seq_along(objs)) {
  x1 <- x
  x1$object <- x1$object[i, ]
  x1$edge <- x1$object %>% dplyr::inner_join(x1$object_link_edge,
  "object_") %>% dplyr::inner_join(x1$edge, "edge_")
  ordered_verts <- t(apply(as.matrix(x1$edge[c(".vx0",
  ".vx1")]), 1, sort))
  x1$vertex <- x$vertex[x$vertex$vertex_ %in% c(ordered_verts),
  ]
  dots <- list()
  dots[["a"]] <- max_area
  dots[["x"]] <- x1
  RTri <- do.call(edge_RTriangle, dots)
  objs[[i]] <- RTri
}
plot(RTri)

plot(sf, reset = F)
drop_inf <- function(x) {
  ## infinite rays in Voronoi
  x[x[,1] > 0 & x[,2] > 0, ]
}
VP <- RTri$VP
VE <- drop_inf(RTri$VE)
segments_mesh <- function(P, I, ..., add = TRUE) {
  if (!add) plot(P)
  segments(P[I[,1], 1], P[I[,1], 2], 
           P[I[,2], 1], P[I[,2], 2], ...)
}

segments_mesh(VP, VE, add = TRUE, col = "firebrick")

I don't think I have the sf->RTriangle pathway separated out anywhere, this is is just ripped out of anglr, and assuming a single sf polygon df.

In the end it's just the VE Voronoi Edges without the infinite rays, built from a "dense-enough" segmentized polygon.

Just want to record this for now and come back later ...

[mpadge]

Wow, yeah, a royal nightmare, but looking strong! I'm going to leave this issue and #1 for sometime in the hopefully not-too-distant futures. Functionality is plenty enough to plough on ahead regardless, and these are just desirable embellishments. I've got a reasonably clear C++ idea for end-point identification - let me know if you need that, and I'll code it up, otherwise I'll also put that aside for the moment. The one issue I'd like done before submission is #3 - I'll ping you there. Thanks loads!