#
#        edgeRipley.R
#
#    $Revision: 1.5 $    $Date: 2008/04/02 13:42:05 $
#
#    Ripley isotropic edge correction weights
#
#  edge.Ripley(X, r, W)      compute isotropic correction weights
#                            for centres X[i], radii r[i,j], window W
#
#  To estimate the K-function see the idiom in "Kest.S"
#
#######################################################################

edge.Ripley <- function(X, r, W=X$window, method="C") {
  # X is a point pattern, or equivalent
  X <- as.ppp(X, W)
  W <- X$window

  switch(W$type,
         rectangle={},
         polygonal={
           if(method != "C")
             stop(paste("Ripley isotropic correction for polygonal windows",
                        "requires method = ", dQuote("C")))
         },
         mask={
           stop(paste("sorry, Ripley isotropic correction",
                      "is not implemented for binary masks"))
         }
         )

  n <- X$n

  if(is.matrix(r) && nrow(r) != n)
    stop("the number of rows of r should match the number of points in X")
  if(!is.matrix(r)) {
    if(length(r) != n)
      stop("length of r is incompatible with the number of points in X")
    r <- matrix(r, nrow=n)
  }

  stopifnot(method %in% c("interpreted", "C"))

  ##########
  
  x <- X$x
  y <- X$y
  
  ############### interpreted R code for rectangular case ##############
  
  if(method == "interpreted") {

    # perpendicular distance from point to each edge of rectangle
    # L = left, R = right, D = down, U = up
    dL  <- x - W$xrange[1]
    dR  <- W$xrange[2] - x
    dD  <- y - W$yrange[1]
    dU  <- W$yrange[2] - y

    # detect whether any points are corners of the rectangle
    small <- function(x) { abs(x) < .Machine$double.eps }
    corner <- (small(dL) + small(dR) + small(dD) + small(dU) >= 2)
  
    # angle between (a) perpendicular to edge of rectangle
    # and (b) line from point to corner of rectangle
    bLU <- atan2(dU, dL)
    bLD <- atan2(dD, dL)
    bRU <- atan2(dU, dR)
    bRD <- atan2(dD, dR)
    bUL <- atan2(dL, dU)
    bUR <- atan2(dR, dU)
    bDL <- atan2(dL, dD)
    bDR <- atan2(dR, dD)

    # The above are all vectors [i]
    # Now we compute matrices [i,j]

    # half the angle subtended by the intersection between
    # the circle of radius r[i,j] centred on point i
    # and each edge of the rectangle (prolonged to an infinite line)

    hang <- function(d, r) {
      answer <- matrix(0, nrow=nrow(r), ncol=ncol(r))
      # replicate d[i] over j index
      d <- matrix(d, nrow=nrow(r), ncol=ncol(r))
      hit <- (d < r)
      answer[hit] <- acos(d[hit]/r[hit])
      answer
    }

    aL <- hang(dL, r)
    aR <- hang(dR, r)
    aD <- hang(dD, r) 
    aU <- hang(dU, r)

    # apply maxima
    # note: a* are matrices; b** are vectors;
    # b** are implicitly replicated over j index
    cL <- pmin(aL, bLU) + pmin(aL, bLD)
    cR <- pmin(aR, bRU) + pmin(aR, bRD)
    cU <- pmin(aU, bUL) + pmin(aU, bUR)
    cD <- pmin(aD, bDL) + pmin(aD, bDR)

    # total exterior angle
    ext <- cL + cR + cU + cD

    # add pi/2 for corners 
    if(any(corner))
    ext[corner,] <- ext[corner,] + pi/2

    # OK, now compute weight
    weight <- 1 / (1 - ext/(2 * pi))
    return(weight)
  }
  ############ C code #############################

  DUP <- spatstat.options("dupC")
  
  switch(W$type,
         rectangle={
           z <- .C("ripleybox",
                   nx=as.integer(n),
                   x=as.double(x),
                   y=as.double(y),
                   rmat=as.double(r),
                   nr=as.integer(ncol(r)),
                   xmin=as.double(W$xrange[1]),
                   ymin=as.double(W$yrange[1]),
                   xmax=as.double(W$xrange[2]),
                   ymax=as.double(W$yrange[2]),
                   epsilon=as.double(.Machine$double.eps),
                   out=as.double(numeric(length(r))),
                   DUP=DUP,
                   PACKAGE="spatstat")
           weight <- matrix(z$out, nrow(r), ncol(r))
         },
         polygonal={
           Y <- as.psp(W)
           z <- .C("ripleypoly",
                   nc=as.integer(n),
                   xc=as.double(x),
                   yc=as.double(y),
                   nr=as.integer(ncol(r)),
                   rmat=as.double(r),
                   nseg=as.integer(Y$n),
                   x0=as.double(Y$ends$x0),
                   y0=as.double(Y$ends$y0),
                   x1=as.double(Y$ends$x1),
                   y1=as.double(Y$ends$y1),
                   out=as.double(numeric(length(r))),
                   DUP=DUP,
                   PACKAGE="spatstat")
           angles <- matrix(z$out, nrow(r), ncol(r))
           weight <- 2 * pi/angles
         }
         )
  return(weight)
}