#
# infline.R
#
# Infinite lines
#
# $Revision: 1.15 $ $Date: 2011/05/18 07:57:31 $
#

infline <- function(a=NULL, b=NULL, h=NULL, v=NULL, p=NULL, theta=NULL) {
  if(is.null(a) != is.null(b))
    stop("invalid specification of a,b")
  if(is.null(p) != is.null(theta))
    stop("invalid specification of p,theta")
  if(!is.null(h)) 
    out <- data.frame(a=h, b=0, h=h, v=NA, p=h, theta=pi/2)
  else if(!is.null(v)) 
    out <- data.frame(a=NA,b=NA,h=NA,v=v,p=v,theta=ifelse(v < 0, pi, 0))
  else if(!is.null(a)) {
    # a, b specified
    z <- data.frame(a=a,b=b)
    a <- z$a
    b <- z$b
    theta <- ifelse(b == 0, pi/2, atan(-1/b))
    theta <- theta %% (2*pi)
    theta <- ifelse(theta > pi, theta - pi, theta)
    p <- a * sin(theta)
    out <- data.frame(a=a, b=b,
                      h=ifelse(b==0, a, NA),
                      v=NA, p=p, theta=theta)
  } else if(!is.null(p)) {
    # p, theta specified
    z <- data.frame(p=p,theta=theta)
    p <- z$p
    theta <- z$theta
    theta <- theta %% (2*pi)
    if(any(reverse <- (theta >= pi))) {
      theta[reverse] <- theta[reverse] - pi
      p[reverse]     <- -p[reverse]
    }
    vert <- (theta == 0)
    horz <- (cos(theta) == 0)
    gene <- !(vert | horz)
    v <- ifelse(vert, p, NA)
    h <- ifelse(horz, p, NA)
    a <- ifelse(gene, p/sin(theta), NA)
    b <- ifelse(gene, -cos(theta)/sin(theta), NA)
    out <- data.frame(a=a,b=b,h=h,v=v,p=p,theta=theta)
  } else stop("No data given!")
  class(out) <- c("infline", class(out))
  return(out)
}

is.infline <- function(x) { inherits(x, "infline") }

plot.infline <- function(x, ...) {
  for(i in seq_len(nrow(x))) {
    xi <- x[i, 1:4]
    xi <- lapply(as.list(xi), function(z){if(is.na(z)) NULL else z})
    do.call("abline", append(xi, list(...)))
  }
  return(invisible(NULL))
}

print.infline <- function(x, ...) {
  n <- nrow(x)
  cat(paste(if(n > 1) n else NULL, "infinite ",
            ngettext(n, "line", "lines"), "\n"))
  print(as.data.frame(x), ...)
  return(invisible(NULL))
}

clip.infline <- function(L, win) {
  # clip a set of infinite straight lines to a window
  win <- as.owin(win)
  stopifnot(inherits(L, "infline"))
  # determine circumcircle of win
  xr <- win$xrange
  yr <- win$yrange
  xmid <- mean(xr)
  ymid <- mean(yr)
  width <- diff(xr)
  height <- diff(yr)
  rmax <- sqrt(width^2 + height^2)/2
  boundbox <- owin(xmid + c(-1,1) * rmax, ymid + c(-1,1) * rmax)
  # compute intersection points with circumcircle 
  p <- L$p
  theta <- L$theta
  hit <- (abs(p) < rmax)
  if(!any(hit)) 
    return(psp(numeric(0),numeric(0),numeric(0),numeric(0), window=win))
  p <- p[hit]
  theta <- theta[hit]
  q <- sqrt(rmax^2 - p^2)
  co <- cos(theta)
  si <- sin(theta)
  X <- psp(x0= xmid + p * co + q * si,
           y0= ymid + p * si - q * co,
           x1= xmid + p * co - q * si,
           y1= ymid + p * si + q * co,
           window=boundbox, check=FALSE)
  # clip to window
  X <- X[win]
  return(X)
}
  
chop.tess <- function(X, L) {
  stopifnot(is.infline(L))
  stopifnot(is.tess(X)||is.owin(X))
  X <- as.tess(X)

  if(!spatstat.options("gpclib")) {
    # polygonal computations unavailable
    # convert to image tessellation
    X <- tess(image=as.im(X))
  }
  
  if(X$type == "image") {
    Xim <- X$image
    xr <- Xim$xrange
    yr <- Xim$yrange
    # extract matrices of pixel values and x, y coordinates
    Zmat <- as.integer(as.matrix(Xim))
    xmat <- rasterx.im(Xim)
    ymat <- rastery.im(Xim)
    # process lines
    for(i in seq_len(nrow(L))) {
      # line i chops window into two pieces
      if(!is.na(h <- L[i, "h"])) {
        # horizontal line
        if(h > yr[1] && h < yr[2]) 
          Zmat <- 2 * Zmat + (ymat > h)
      } else if(!is.na(v <- L[i, "v"])) {
        # vertical line
        if(v > xr[1] && v < xr[2])
          Zmat <- 2 * Zmat + (xmat < h)
      } else {
        # generic line y = a + bx
        a <- L[i, "a"]
        b <- L[i, "b"]
        Zmat <- 2 * Zmat + (ymat > a + b * xmat)
      }
    }
    # Now just put back as factor image
    Zim <- im(Zmat, xcol=Xim$xcol, yrow=Xim$yrow, unitname=unitname(Xim))
    Z <- tess(image=Zim)
    return(Z)
  }

  #---- polygonal computation --------
  # get bounding box
  B <- as.rectangle(as.owin(X))
  xr <- B$xrange
  yr <- B$yrange

  # get coordinates
  for(i in seq_len(nrow(L))) {
    # line i chops box B into two pieces
    if(!is.na(h <- L[i, "h"])) {
      # horizontal line
      if(h < yr[1] || h > yr[2])
        Z <- NULL
      else {
        lower <- owin(xr, c(yr[1], h))
        upper <- owin(xr, c(h, yr[2]))
        Z <- tess(tiles=list(lower,upper), window=B)
      }
    } else if(!is.na(v <- L[i, "v"])) {
      # vertical line
      if(v < xr[1] || v > xr[2])
        Z <- NULL
      else {
        left <- owin(c(xr[1], v), yr)
        right <- owin(c(v, xr[2]), yr)
        Z <- tess(tiles=list(left,right), window=B)
      }
    } else {
      # generic line
      a <- L[i, "a"]
      b <- L[i, "b"]
      # Intersect with extended left and right sides of B
      yleft <- a + b * xr[1]
      yright <- a + b * xr[2]
      ylo <- min(yleft, yright, yr[1]) - 1
      yhi <- max(yleft, yright, yr[2]) + 1
      lower <- owin(poly=list(x=xr[c(1,1,2,2)],
                              y=c(yleft,ylo,ylo,yright)))
      upper <- owin(poly=list(x=xr[c(1,2,2,1)],
                              y=c(yleft,yright,yhi,yhi)))
      Bplus <- owin(xr, c(ylo, yhi))
      Z <- tess(tiles=list(lower,upper), window=Bplus)
    }
    # intersect this simple tessellation with X
    if(!is.null(Z))
      X <- intersect.tess(X, Z)
  }
  return(X)
}