R/selfIntersections.R

#' Find if two line segments intersect
#' @title Check if line segments intersect
#' @param a first coordinate of first line segment
#' @param b second coordinate of first line segment
#' @param c first coordinate of second line segment
#' @param d second coordinate of second line segment
#' @param seq indicator for whether the two line segments are intersecting
#' @return boolean indicating if there is an intersection
#' @export
#' @examples
#'  check_intersect(c(0, 0), c(1, 1), c(2, 2), c(3, 3))
#'  check_intersect(c(0, 0), c(1, 1), c(0.5, 0.5), c(2, 2))
check_intersect <- function(a, b, c, d, seq = FALSE) {
  ##typical lines refers to lines that are not horizontal or vertical

  ##inequality functions that correctly account for floating point issues.
  ##(Care with equality of floating points is needed througout. This is why all.equal is used)
  less_eq <- function(a, b) {
    ((a < b) || isTRUE(all.equal(as.numeric(a), as.numeric(b))))
  }
  greater_eq <- function(a, b) {
    ((a > b) || isTRUE(all.equal(as.numeric(a), as.numeric(b))))
  }

  ##return true for identical points
  if (isTRUE(all.equal(a, c)) & isTRUE(all.equal(b, d))) {
    return(TRUE)

    ##two vertical lines
  } else if (isTRUE(all.equal(a[1], b[1])) & isTRUE(all.equal(c[1], d[1]))) {
    if (!isTRUE(all.equal(a[1], c[1]))) { #check if the two lines are not on the same x-coord
      return(FALSE)
      #if on same x-coord; y-coords need to overlap
    } else if ((max(a[2], b[2]) < min(c[2], d[2])) ||
               (max(c[2], d[2]) < min(a[2], b[2]))) {
      return(FALSE) #if sequential equality means not an overlap
    } else if ((less_eq(max(a[2], b[2]), min(c[2], d[2])) ||
                less_eq(max(c[2], d[2]), min(a[2], b[2]))) & seq == TRUE) {
      return(FALSE)
    } else {
      return(TRUE)
    }

    ##two horizontal lines
  } else if (isTRUE(all.equal(a[2], b[2])) & isTRUE(all.equal(c[2], d[2]))) {
    if (!(isTRUE(all.equal(a[2], c[2])))) { #check if the two lines are on the same y-coord
      return(FALSE)
      #if on same y-coord; x-coords need to overlap
    } else if ((max(a[1], b[1]) < min(c[1], d[1])) ||
               (max(c[1], d[1]) < min(a[1], b[1]))) {
      return(FALSE) #if sequential equality means not an overlap
    } else if ((less_eq(max(a[1], b[1]), min(c[1], d[1])) ||
                less_eq(max(c[1], d[1]), min(a[1], b[1]))) & seq == TRUE) {
      return(FALSE)
    } else {
      return(TRUE)
    }

    ##first line horizontal, second line vertical
  } else if (isTRUE(all.equal(a[2], b[2])) & isTRUE(all.equal(c[1], d[1]))) {
    if (greater_eq(c[1], min(a[1], b[1])) & less_eq(c[1], max(a[1], b[1])) &
        greater_eq(a[2], min(c[2], d[2])) & less_eq(a[2], max(c[2], d[2]))) {
      inter <- c(c[1], a[2])
      if (seq & isTRUE(all.equal(inter, b))) { #okay if line intersects only at end points for sequential coords
        return(FALSE)
      } else {
        return(TRUE)
      }
    } else {
      return(FALSE)
    }

    ##first line vertical, second line horizontal
  } else if (isTRUE(all.equal(c[2], d[2])) & isTRUE(all.equal(a[1], b[1]))) {
    if (greater_eq(a[1], min(c[1], d[1])) & less_eq(a[1], max(c[1], d[1])) &
        greater_eq(c[2], min(a[2], b[2])) & less_eq(c[2], max(a[2], b[2]))) {
      inter <- c(a[1], c[2])
      if (seq & isTRUE(all.equal(inter, d))) { #okay if line intersects only at end points for sequential coords
        return(FALSE)
      } else {
        return(TRUE)
      }
    } else {
      return(FALSE)
    }

    ##first line horizontal, second typical
  } else if (isTRUE(all.equal(a[2], b[2]))) {
    #find typical line
    m2 <- (d[2] - c[2])/(d[1] - c[1])
    b2 <- c[2] - m2*c[1]
    #find possible intersection (ya = m2*xi + b2)
    xi <- (a[2] - b2)/m2
    yi <- m2*xi + b2
    #check if intersect point is in region; if so, return intersect point
    if (!(less_eq(xi, max(a[1], b[1])) & greater_eq(xi, min(a[1], b[1])) &
          less_eq(yi, max(a[2], b[2])) & greater_eq(yi, min(a[2], b[2])) &
          less_eq(xi, max(c[1], d[1])) & greater_eq(xi, min(c[1], d[1])) &
          less_eq(yi, max(c[2], d[2])) & greater_eq(yi, min(c[2], d[2])))) {
      return(FALSE)
    } else {
      inter <- c(xi, yi)
      if (seq & isTRUE(all.equal(inter, b))) { #okay if line intersects only at end points for sequential coords
        return(FALSE)
      } else {
        return(TRUE)
      }
    }

    ##first line typical, second horizontal
  } else if (isTRUE(all.equal(c[2], d[2]))) {
    #find typical line
    m2 <- (b[2] - a[2])/(b[1] - a[1])
    b2 <- a[2] - m2*a[1]
    #find possible intersection (ya = m2*xi + b2)
    xi <- (c[2] - b2)/m2
    yi <- m2*xi + b2
    #check if intersect point is in region; if so, return intersect point
    if (!(less_eq(xi, max(c[1], d[1])) & greater_eq(xi, min(c[1], d[1])) &
          less_eq(yi, max(c[2], d[2])) & greater_eq(yi, min(c[2], d[2])) &
          less_eq(xi, max(a[1], b[1])) & greater_eq(xi, min(a[1], b[1])) &
          less_eq(yi, max(a[2], b[2])) & greater_eq(yi, min(a[2], b[2])))) {
      return(FALSE)
    } else {
      inter <- c(xi, yi)
      if (seq & isTRUE(all.equal(inter, c))) { #okay if line intersects only at end points for sequential coords
        return(FALSE)
      } else {
        return(TRUE)
      }
    }

    ##first line typical, second vertical
  } else if (isTRUE(all.equal(c[1], d[1]))) {
    #find other line
    m1 <- (b[2] - a[2])/(b[1] - a[1])
    b1 <- a[2] - m1*a[1]
    #find possible intersection (yi = m1*xc + b1)
    yi <- m1*c[1] + b1
    #check if intersect point is in region; if so, return intersect point
    if (!(less_eq(yi, max(a[2], b[2])) & greater_eq(yi, min(a[2], b[2])) &
          less_eq(yi, max(c[2], d[2])) & greater_eq(yi, min(c[2], d[2])))) {
      return(FALSE)
    } else {
      inter <- c(c[1], yi)
      if (seq & isTRUE(all.equal(inter, b))) { #okay if line intersects only at end points for sequential coords
        return(FALSE)
      } else {
        return(TRUE)
      }
    }

    ##first line vertical, second line typical
  } else if (isTRUE(all.equal(a[1], b[1]))) {
    #find other line
    m2 <- (d[2] - c[2])/(d[1] - c[1])
    b2 <- c[2] - m2*c[1]
    #find possible intersection (yi = m2*xa + b2)
    yi <- m2*a[1] + b2
    #check if intersect point is in region; if so, return intersect point
    if (!(less_eq(yi, max(a[2], b[2])) & greater_eq(yi, min(a[2], b[2])) &
          less_eq(yi, max(c[2], d[2])) & greater_eq(yi, min(c[2], d[2])))) {
      return(FALSE)
    } else {
      inter <- c(a[1], yi)
      if (seq & isTRUE(all.equal(inter, b))) { #okay if line intersects only at end points for sequential coords
        return(FALSE)
      } else {
        return(TRUE)
      }
    }

    ##two typical lines
  } else {
    #find eq's of both lines
    m1 <- (b[2] - a[2])/(b[1] - a[1])
    m2 <- (d[2] - c[2])/(d[1] - c[1])
    b1 <- a[2] - m1*a[1]
    b2 <- c[2] - m2*c[1]

    #points are on the same line
    if (isTRUE(all.equal(m1, m2)) & isTRUE(all.equal(b1, b2))) {
      #points overlap when min(a, b) > max(c, d) or vice versa
      #It's sufficient to test x or y coords only, since coords are only the same line
      if (!((greater_eq(min(a[2], b[2]), max(c[2], d[2])) || greater_eq(min(c[2], d[2]), max(a[2], b[2]))))) {
        return(TRUE)
      } else if (isTRUE(all.equal(b, c)) & seq == FALSE) {
        return(TRUE)
      } else {
        return(FALSE)
      }

    } else {
      #find possible point of intersection (Yi = m1*xi + b1 & Yi = m2*xi + b2, solve for xi)
      xi <- (b2 - b1)/(m1 - m2)
      #check if possible intersect point is in region; if so, return intersect point
      if (!((less_eq(xi, max(a[1], b[1])) & greater_eq(xi, min(a[1], b[1]))) &
            (less_eq(xi, max(c[1], d[1])) & greater_eq(xi, min(c[1], d[1]))))) {
        return(FALSE)
      } else {
        inter <- c(xi, m1*xi + b1)
        if (seq & isTRUE(all.equal(inter, b))) { #okay if line intersects only at end points for sequential coords
          return(FALSE)
        } else {
          return(TRUE)
        }
      }
    }
  }
}

#' Determine if there are any intersecting line segments in a matrix of coordinates representing a line
#' @title Check if a line has intersecting segments
#' @param line matrix of coordinates corresponding to the line of interest
#' @return list where \code{list$any} is a boolean indicating if there are any intersections
#' and \code{list$val} is an index corresponding to the first intersection found
#' @importFrom utils combn
#' @examples
#' check_results <- any_intersect(currSecEx)
#' check_results$any #true/false
#' check_results$val #indices of first intersection found
any_intersect <- function(line) {
  ##line segments are referred to by the index of their first point
  n_points <- nrow(line)

  ##only two points, no way to have an intersection
  if(n_points == 2) {
    return(list("any" = FALSE, "val" = NA))
  }

  ##check all pairs of line segments for intersections
  pairs <- combn(1:(n_points - 1), 2) #indices of all the combinations of ways pairs could be lined up
  for (i in 1:ncol(pairs)) { #loop through line seqment pairs, stop if an intersection is found
    curr_pair <- pairs[,i]
    seq <- ((curr_pair[2] - curr_pair[1]) == 1) #boolean to indicate if points are sequential (meaning they should have an intersection)
    if (check_intersect(a = line[curr_pair[1], ], b = line[curr_pair[1] + 1,],
                       c = line[curr_pair[2],], d = line[curr_pair[2] + 1,], seq)) {
      #return indices of which line segments are the source of the intersection
      return(list("any" = TRUE, "val" = pairs[, i]))
    }
  }

  #return FALSE and NA if no intersections are found
  return(list("any" = FALSE, "val" = NA))
}


#' Function to remove all self-intersections from a contour.
#' @title Remove self-intersections
#' @param my_poly \code{SpatialPolygons} object from which self-intersections need to be removed
#' @param plotting boolean indicating if results should be plotted
#' @param poly_name name for \code{SpatialPolygons} object to return (defaults to "unspecified")
#' @param min_area minimum area for any individual polygon
#' @return \code{SpatialPolygons} object with self-intersections removed
#' @importFrom rgeos gIsValid
#' @importFrom sp Polygon
#' @export
#' @examples
#' \dontrun{
#' par(mfrow = c(1, 2))
#' plot(interEx, main = "Original Contour")
#' noInter <- untwist(interEx, poly_name = "interEx")
#' plot(noInter, main = "Final Contour")
#' }
untwist <- function(my_poly, plotting = FALSE, poly_name = "unspecified",
                    min_area = 12.5) {
  if (is.null(my_poly)) {
    return(NULL)
  }

  if (get_area(my_poly) < min_area) {
    return(NULL)
  }

  while (suppressWarnings(!gIsValid(my_poly))) {

    ##pull out coordinates
    coords <- my_poly@polygons[[1]]@Polygons[[1]]@coords
    coords <- matrix(as.numeric(coords), ncol = 2)

    ##return NULL if you start with just a horizontal or vertinal line
    if ((length(unique(coords[, 1])) == 1) || (length(unique(coords[, 2])) == 1)) {
      return (NULL)
    }

    ##find point at or near self intersection point
    prob <- suppressWarnings(gIsValid(my_poly, reason = T))
    temp <- strsplit(prob, " ")
    if (length(temp[[1]]) == 3) { #special cases with more text (e.g. "ring)
      prob <- temp[[1]][2:3]
    }
    prob <- as.numeric(unlist(strsplit(gsub("Self-intersection\\[|\\]", "", prob), " ")))
    prob <- prob[!is.na(prob)]

    ##find minimal region around intersection
    does_intersect <- FALSE
    n_neigh <- 5
    curr_sec <- c()
    while (!does_intersect) {
      index <- order(apply(coords, 1, function(x){  sqrt((x[1] - prob[1])^2 + (x[2] - prob[2])^2)}),
                     decreasing = FALSE)[1:n_neigh]
      index <- index[!is.na(index)] #remove NA's, which are put in place if more neightbors requested than points
      if (((nrow(coords) - max(index)) + min(index) <= max(index) - min(index)) & #Check if looping around end (BUT not needing all points)
          ((max(index) - min(index) + 1) !=  nrow(coords))) {
        cat <-  ((nrow(coords) - index) <= index) #categorize each point as closer to the end of the contour or the beginning
        curr <- c(min(index[cat == TRUE]):nrow(coords), 1:max(index[cat == FALSE]))
      } else { #standard case (no looping)
        curr <- min(index):max(index)
      }
      curr_sec <- coords[curr,]
      temp <- any_intersect(curr_sec)
      if (temp$any) {

        ##if we've found the intersection, keep the minimal number of line segments
        ep <- rep(NA, 4)
        ep[1] <- curr[temp$val[1]] #first coord of intersecting line segment 1
        ep[3] <- curr[temp$val[2]] #first coord of intersecting line segment 2
        if (ep[1] != nrow(coords)) {#use max + 1, since line segments are indexed by their first point (checking for looping over)
          ep[2] <- ep[1] + 1
        } else {
          ep[2] <- 1
        }
        if (ep[3] != nrow(coords)) {
          ep[4] <- ep[3] + 1
        } else {
          ep[4] <- 1
        }
        if ((nrow(coords) - max(ep)) + min(ep) < max(ep) - min(ep) & (nrow(coords) > length(ep) + 1)) {
          #Check if indices loop back to beginning. The second section condition ensures that there is actually
          #a line segment to add back not just a point, in which case all points should be replaced
          cat <-  ((nrow(coords) - ep) <= ep) #categorize each point as closer to the end of the contour or the beginning
          curr <- c(min(ep[cat]):nrow(coords), 1:max(ep[!cat]))
        } else { #standard case (no looping)
          curr <- min(ep):max(ep)
        }
        curr_sec <- coords[curr, ]
        stopifnot(any_intersect(curr_sec)$any) #you should always have an intersection at this point
        does_intersect <- TRUE
      } else { #otherwise expand the region
        n_neigh <- n_neigh + 5
        if (n_neigh > nrow(coords)) {
          n_neigh <- nrow(coords)
        }
      }
    }
    new <- untwist_sec(curr_sec)
    if (curr[1] > curr[length(curr)]) {
      #remove section being adjusted and put new points on the end
      cat <-  ((nrow(coords) - curr) <= curr) #categorize each point as closer to the end of the contour or the beginning
      coords <- rbind(coords[(max(curr[cat == FALSE]) + 1):(min((curr[cat == TRUE])) - 1),],  new)
    } else if (curr[length(curr)] == nrow(coords)) { #last point
      coords <- rbind(coords[1:(min(curr) - 1),], new)
    } else if (curr[1] == 1) { #check if indices start at one
      coords <- rbind(new, coords[(max(curr) + 1):nrow(coords),])
    } else { #no looping
      coords <- rbind(coords[1:(min(curr) - 1),], new, coords[(max(curr) + 1):nrow(coords), ])
    }

    ##if all the algorithm leads to is a single vertical line or horiztonal line,  assume no polygon
    if (length(unique(coords[,1])) == 1 || length(unique(coords[,2])) == 1) {
      return(NULL)
    }

    #If all the algorithm, leads to is a single set of two connected points, assume no polygon
    if (nrow(unique(round(coords, 2))) <= 3) {
      return(NULL)
    }

    ##make new myPoly and start over
    my_poly <- SpatialPolygons(list(
      Polygons(list(Polygon(coords)), ID = poly_name))
    )
  }

  if (!is.null(my_poly)) {
    if (get_area(my_poly) < min_area) {
      return(NULL)
    }
  }

  return(my_poly)
}

#' Function to correct self-intersections in a section of a line.
#' @title Remove self-intersections from one section of a contour
#' @param line  N x 2 matrix of coordinates
#' @param tol how much of a difference between the original line and the simplified line is allowed
#' @param eps how much to increase \code{tol} by on each iteration
#' @return n x 2 matrix of the new coordinates with self-intersections removed
#' @importFrom sp Line Lines SpatialLines
#' @importFrom rgeos gSimplify
#' @export
#' @examples
#' par(mfrow = c(1, 2))
#' plot(currSecEx, type = "l", main = "Original Line Section", xlab = "", ylab = "")
#' new_sec <- untwist_sec(currSecEx)
#' plot(new_sec, type = "l", main = "New Line Section", xlab = "", ylab =  "")
untwist_sec <- function(line, tol = 0, eps = .25) {

  l <- SpatialLines(list(Lines(Line(line), ID = "temp")))
  while (any_intersect(line)$any) {
    tol <- tol + eps
    l <- gSimplify(l, tol)
    line <- l@lines[[1]]@Lines[[1]]@coords
  }

  return(line)
}

Try the IceCast package in your browser

Any scripts or data that you put into this service are public.

IceCast documentation built on June 24, 2019, 9:03 a.m.