head_stretching_recycle.R

tfile <- sprintf("%s.rda", tempfile())
download.file("https://github.com/mdsumner/scsf/raw/master/data/minimal_mesh.rda", tfile, mode = "wb")
load(tfile)
library(sf)

## Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.3

inner_ring_touching <- st_sf(a = 1, geometry = st_sfc(st_multipolygon(unlist(lapply(st_geometry(minimal_mesh), unclass), recursive = FALSE))))

#plot(inner_ring_touching)
st_is_valid(inner_ring_touching)

## Warning in evalq((function (..., call. = TRUE, immediate. = FALSE,
## noBreaks. = FALSE, : Self-intersection at or near point 0.68999999999999995
## 0

## [1] FALSE

# devtools::install_github("mdsumner/scsf")
library(scsf)

## Loading required package: sc

library(sc)
prim <- PRIMITIVE(inner_ring_touching)

## which segments have to go? find where interchanging vertex pairs is not unique
prim$segment$useg <-  apply(cbind(prim$segment$.vertex0, prim$segment$.vertex1), 1, function(x) paste(sort(x), collapse = "-"))
prim$segment$bad <- duplicated(prim$segment$useg)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

### aa is a matrix of vertex id
## this is slow and a bit bogus, note that there's a problem with 
## multiple matches - but they should be covered by filter out those
## segments in the first place as per "useg" anti join
cycles <- function(aa) {
  ii <- 1
  set0 <- ii
  visited <- logical(nrow(aa))
  while(!all(visited)) {
    i0 <- ii
    repeat {
      ii <- which(aa[,1] == aa[ii, 2])
      if (length(ii) < 1 | ii[1] == i0) {
        set0 <- c(set0, NA_integer_)
        break; 
      }
      set0 <- c(set0, ii)
    }
    visited <- seq(nrow(aa)) %in% na.omit(set0)
    ii <- which(!visited)[1L]
    if (!is.na(ii)) set0 <- c(set0, ii)
  }
  l <- split(set0, c(0, cumsum(abs(diff(is.na(set0))))))
  bind_rows(lapply(l[!unlist(lapply(l, function(x) all(is.na(x))))], function(x) tibble(row = x)), .id = "cycle")
}

## these are the ids of the vertex pairs to be kept, one row per segment
#bb <- prim$segment %>% dp %>% dplyr::select(.vertex0, .vertex1) %>% as.matrix()
## I'm failing to get anti_join to work today for some reason
anti <- prim$segment %>% filter(bad)  %>% inner_join(prim$segment %>% select(useg))

## Joining, by = "useg"

segs <- prim$segment %>% filter( !prim$segment$useg %in% anti$useg) 
bb <- prim$segment %>% filter( !prim$segment$useg %in% anti$useg) %>% dplyr::select(.vertex0, .vertex1) %>% as.matrix() 
## these are the indices from bb of each valid
ind <- cycles(bb)

#segs <- prim$segment %>% dplyr::filter(!bad) 

new_paths_list <- lapply(split(ind$row, ind$cycle), function(index) tibble(vertex_path = segs[index, c(".vertex0", ".vertex1")] %>% as.matrix() %>% t() %>% as.vector()))
plot(inner_ring_touching)

new_paths <- bind_rows(lapply(new_paths_list, 
                              function(vertex_tab) prim$vertex[match(vertex_tab$vertex_path[!duplicated(vertex_tab$vertex_path)], prim$vertex$vertex_), c("x_", "y_")]), .id = "newpath")
library(ggplot2)
#ggplot(new_paths, aes(x = x_, y = y_, group = newpath)) + geom_polygon()

library(ggpolypath)
ggplot(new_paths, aes(x = x_, y = y_, group = newpath, col = newpath)) + ggpolypath::geom_polypath()

#polygon(prim$vertex[match(t1$vertex_path, prim$vertex$vertex_), c("x_", "y_")], lwd = 4)
#polygon(prim$vertex[match(t2$vertex_path, prim$vertex$vertex_), c("x_", "y_")], lwd = 4)
#polygon(prim$vertex[match(t4$vertex_path, prim$vertex$vertex_), c("x_", "y_")], lwd = 4)