#
#	weights.S
#
#	Utilities for computing quadrature weights
#
#	$Revision: 4.27 $	$Date: 2011/05/18 09:24:01 $
#
#
# Main functions:
		
#	gridweights()	    Divide the window frame into a regular nx * ny
#			    grid of rectangular tiles. Given an arbitrary
#			    pattern of (data + dummy) points derive the
#			    'counting weights'.
#
#	dirichlet.weights() Compute the areas of the tiles of the
#			    Dirichlet tessellation generated by the 
#			    given pattern of (data+dummy) points,
#			    restricted to the window.
#	
# Auxiliary functions:	
#			
#       countingweights()   compute the counting weights
#                           for a GENERIC tiling scheme and an arbitrary
#			    pattern of (data + dummy) points,
#			    given the tile areas and the information
#			    that point number k belongs to tile number id[k]. 
#
#
#	gridindex()	    Divide the window frame into a regular nx * ny
#			    grid of rectangular tiles. 
#			    Compute tile membership for arbitrary x,y.
#				    
#       grid1index()        1-dimensional analogue of gridindex()
#
#
#-------------------------------------------------------------------
	
countingweights <- function(id, areas, check=TRUE) {
	#
	# id:        cell indices of n points
	#                     (length n, values in 1:k)
	#
	# areas:     areas of k cells 
	#                     (length k)
	#
    id <- as.integer(id)
    fid <- factor(id, levels=seq_along(areas))
    counts <- table(fid)
    w <- areas[id] / counts[id]     # ensures denominator > 0
    w <- as.vector(w)
#	
# that's it; but check for funny business
#
    zerocount <- (counts == 0)
    zeroarea <- (areas == 0)
    if(any(!zeroarea & zerocount))
	warning("some tiles with positive area do not contain any points")
    if(any(!zerocount & zeroarea)) {
	warning("Some tiles with zero area contain points")
	warning("Some weights are zero")
	attr(w, "zeroes") <- zeroarea[id]
    }
#
    names(w) <- NULL
    w
}

gridindex <- function(x, y, xrange, yrange, nx, ny) {
	#
	# The box with dimensions xrange, yrange is divided
	# into nx * ny cells.
	#
	# For each point (x[i], y[i]) compute the index (ix, iy)
	# of the cell containing the point.
	# 
	ix <- grid1index(x, xrange, nx)
	iy <- grid1index(y, yrange, ny)
	#
	return(list(ix=ix, iy=iy, index=as.integer((iy-1) * nx + ix)))
}

grid1index <- function(x, xrange, nx) {
	i <- ceiling( nx * (x - xrange[1])/diff(xrange))
	i <- pmax(1, i)
	i <- pmin(i, nx)
	i
}

gridweights <- function(X, ntile=NULL, ..., window=NULL, verbose=FALSE,
                        npix=NULL, areas=NULL) {
	#
	# Compute counting weights based on a regular tessellation of the
	# window frame into ntile[1] * ntile[2] rectangular tiles.
	#
	# Arguments X and (optionally) 'window' are interpreted as a
	# point pattern.
	#
	# The window frame is divided into a regular ntile[1] * ntile[2] grid
	# of rectangular tiles. The counting weights based on this tessellation
	# are computed for the points (x, y) of the pattern.
	#
        # 'npix' determines the dimensions of the pixel raster used to
        # approximate tile areas.
	
	X <- as.ppp(X, window)
	x <- X$x
	y <- X$y
        win <- X$window

        # determine number of tiles
        if(is.null(ntile))
          ntile <- default.ntile(X)
        if(length(ntile) == 1)
          ntile <- rep(ntile, 2)
        nx <- ntile[1]
        ny <- ntile[2]

        if(verbose) 
          cat(paste("grid weights for a", nx, "x", ny, "grid of tiles\n"))
        
        if(is.null(areas)) {
  	  # compute tile areas
          if(win$type == "rectangle") {

            tilearea <- area.owin(win)/(nx * ny)
            areas <- rep(tilearea, nx * ny)

          } else {
            # convert window to mask

            win <- as.mask(win, dimyx=rev(npix))

            if(verbose) {
              if(!is.null(npix)) 
                np <- npix
              else {
                np <- rev(spatstat.options("npixel"))
                if(length(np) == 1) np <- rep(np, 2)
              }
              cat(paste("Approximating window by mask (",
                        np[1], " x ", np[2], " pixels)\n", sep=""))
            }
                
            # extract pixel coordinates inside window
            xx <- as.vector(raster.x(win)[win$m])
            yy <- as.vector(raster.y(win)[win$m])
                                
	    # classify all pixels into tiles
            pixelid <- gridindex(xx, yy, 
                                 win$xrange, win$yrange, nx, ny)$index
            pixelid <- factor(pixelid, levels=seq_len(nx * ny))
                                
	    # compute digital areas of tiles
            tilepixels <- as.vector(table(pixelid))
            pixelarea <- win$xstep * win$ystep
            areas <- tilepixels * pixelarea
          } 
        }
        
	# classify each point according	to its tile

        if(win$type == "mask") {
          # first move each data point to nearest pixel
          Xapprox <- nearest.raster.point(X$x, X$y, win, indices=FALSE)
          x <- Xapprox$x
          y <- Xapprox$y
        }
        id <- gridindex(x, y, win$xrange, win$yrange, nx, ny)$index
        
	# compute counting weights 
	w <- countingweights(id, areas)

        # attach information about weight construction parameters
        attr(w, "weight.parameters") <-
          list(method="grid", ntile=ntile, npix=npix, areas=areas)
        
	return(w)
}


dirichlet.weights <- function(X, window = NULL, exact=TRUE, ...) {
	#
	# Compute weights based on Dirichlet tessellation of the window 
	# induced by the point pattern X. 
	# The weights are just the tile areas.
	#
	# NOTE:	X should contain both data and dummy points,
	# if you need these weights for the B-T-B method.
	#
	# Arguments X and (optionally) 'window' are interpreted as a
	# point pattern.
	#
	# If the window is a rectangle, we invoke Rolf Turner's "deldir"
	# package to compute the areas of the tiles of the Dirichlet
	# tessellation of the window frame induced by the points.
	# [NOTE: the functionality of deldir to create dummy points
	# is NOT used. ]
	#	if exact=TRUE	compute the exact areas, using "deldir"
	#	if exact=FALSE      compute the digital areas using exactdt()
	# 
	# If the window is a mask, we compute the digital area of
	# each tile of the Dirichlet tessellation by counting pixels.
	#
	#
	# 
	#
	
	X <- as.ppp(X, window)
	x <- X$x
	y <- X$y
	win <- X$window

	if(exact && (win$type == "rectangle")) {
		rw <- c(win$xrange, win$yrange)
	        # invoke deldir() with NO DUMMY POINTS
		tessellation <- deldir(x, y, dpl=NULL, rw=rw)
	        # extract tile areas
	        w <- tessellation$summary[, 'dir.area']
	} else {
		# Compute digital areas of Dirichlet tiles.
                win <- as.mask(win)
                X$window <- win
		#
                # Nearest data point to each pixel:
                tileid <- exactdt(X)$i
                # 
		if(win$type == "mask") 
			# Restrict to window (result is a vector - OK)
			tileid <- tileid[win$m]
		# Count pixels in each tile
		id <- factor(tileid, levels=seq_len(X$n))
		counts <- table(id)
                # turn off the christmas lights
                class(counts) <- NULL
                names(counts) <- NULL
                dimnames(counts) <- NULL
		# Convert to digital area
		pixelarea <- win$xstep * win$ystep
		w <- pixelarea * counts
		# Check for zero pixel counts
		zeroes <- (counts == 0)
		if(any(zeroes)) {
			warning("some Dirichlet tiles have zero digital area")
			attr(w, "zeroes") <- zeroes
		}
	} 
        # attach information about weight construction parameters
        attr(w, "weight.parameters") <- list(method="dirichlet", exact=exact)

        return(w)
}

default.ntile <- function(X) { 
	# default number of tiles (n x n) for tile weights
        # when data and dummy points are X
  X <- as.ppp(X)
  guess.ngrid <- 10 * floor(sqrt(X$n)/10)
  max(5, guess.ngrid/2)
}