.packageName <- "geometry"
"Unique" <-
function (X, rows.are.sets = FALSE) 
{
    if (rows.are.sets) 
        X = matsort(X)
    X = X[matorder(X), ]
    dX = apply(X, 2, diff)
    uniq = c(TRUE, ((dX^2) %*% rep(1, ncol(dX))) > 0)
    X = X[uniq, ]
    return(X)
}
"convhulln" <-
function (p, options = " ") 
.Call("convhulln", as.matrix(p), as.character(options))
"delaunayn" <-
function (p, options = "QJ") 
.Call("delaunayn", p, options)
"distmesh2d" <-
function(fd, fh, h0, bbox, p=NULL, pfix=array(0,dim=c(0,2)),
  ..., dptol=.001, ttol=.1, Fscale=1.2, deltat=.2,
  geps=.001*h0, deps=sqrt(.Machine$double.eps)*h0,
  maxiter=1000){

  rownorm2 = function(x) drop(sqrt((x^2)%*%c(1,1)))

  if(any(apply(bbox,2,diff)==0))
     stop("Supplied bounding box has zero area.")

  if(is.null(p)){
    #%1 generate initial grid
    y = seq(bbox[1,2],bbox[2,2], by=h0*sqrt(3)/2)
    x = seq(bbox[1,1], bbox[2,1], by=h0)
    x = matrix(x,length(y),length(x),byrow=TRUE)
    x[seq(2,length(y),by=2),] = x[seq(2,length(y),by=2),] + h0/2
    p = cbind(c(x),y)

    #%2 remove nodes outside boundary specified by fd (points evaluated negative are
    # considered to lie inside the boundary)
    p = p[fd(p,...)<geps,]
  }

  r0 = 1 / fh(p,...)^2                           # acceptance probability
  p = rbind(pfix, p[runif(nrow(p))<r0/max(r0),]) # rejection sampling
  N = nrow(p)
  if(N<=3)
    stop("Not enough starting points inside boundary (is h0 too large?).")

  on.exit(return(invisible(p)));                 # in case we need to stop earlier
  cat("Press esc if the mesh seems fine but the algorithm hasn't converged.\n")
  flush.console();

  #%3 main loop: iterative improvement of points
  pold = 1.0/.Machine$double.eps;
  iter = 0
  while(TRUE){
    if( max( rownorm2(p-pold)/h0 )>ttol ){
        pold = p;

        T = delaunayn(p)                           # generate a Delaunay triangulation

        pmid = (p[T[,1],] + p[T[,2],] + p[T[,3],])/3          # calculate average of node locations as centers
        T  = T[fd(pmid, ...) < (-geps),1:3];                  # remove triangles with center outside region

        #%4 describe edges by uniqe pairs of nodes
        #bars = unique(rbind(T[,-1],T[,-2],T[,-3]),MARGIN=1);
        #bars = bars[order(bars[,1],bars[,2]),];

        bars = rbind(T[,-3],T[,-2],T[,-1])                    # select unique edges
        #too slow: bars = unique(matsort(bars), MARGIN=1)
        bars = Unique(matsort(bars))                          # order the edges according to the node indices

        #%5 Graphical display
        trimesh(T,p)         # a la Matlab
    }

    #%6 compute force F on the basis of edge lenghts
    barvec = p[bars[,1],] - p[bars[,2],]                      # bar vectors
    L = rownorm2(barvec)                                      # their lengths

    # calculate desired lengths L0 by use of fh
    hbars = fh((p[bars[,1],]+p[bars[,2],])/2, ...)
    L0 = hbars * Fscale * sqrt(sum(L^2)/sum(hbars^2));
    F = drop(L0-L); F[F<0] = 0;                               # the forces on the edges = max(L0-L,0)


    Fvec = barvec * (F/L)
    Ftot = matrix(0,N,2);
    ii = bars[,c(1,1,2,2)]
    jj = rep(1,length(F)) %o% c(1,2,1,2)
    s  = c(cbind(Fvec,-Fvec))
#    for(k in 1:length(s))                                     # sum all forces on each node
#        Ftot[ii[k],jj[k]] = Ftot[ii[k],jj[k]] + s[k];
    ns = length(s)
    Ftot[1:(2*N)] = rowsum(s,ii[1:ns]+ns*(jj[1:ns]-1))  # sum all forces on each node
    if(nrow(pfix)>0) Ftot[1:nrow(pfix),] = 0;                 # Force = 0 at fixed points
    p = p + deltat*Ftot;

    #%7 excercise normal force at boundary: move overshoot points to nearest boundary point
    d = fd(p);
    ix= d>0;                                                  # find points outside
    dgradx= (fd(cbind(p[ix,1]+deps, p[ix,2]),...) - d[ix])/deps;  # Numerical
    dgrady= (fd(cbind(p[ix,1], p[ix,2]+deps),...) - d[ix])/deps;  #    gradient
    p[ix,] = p[ix,] - cbind(d[ix]*dgradx, d[ix]*dgrady);      # Project back to boundary

    #%8 test for convergence
    if(max(rownorm2(deltat*Ftot[d < (-geps),])/h0) < dptol | iter>=maxiter) break;
    iter = iter + 1
  }
  if(iter>=maxiter)
    warning(" Maximum iterations reached. Relaxation process not \n completed")
  return(p);
}
"distmeshnd"  <-
function (fdist, fh, h, box, pfix = array(dim = c(0, ncol(box))),
    ..., ptol = 0.001, ttol = 0.1, deltat = 0.1, geps = 0.1 *
        h, deps = sqrt(.Machine$double.eps) * h)
{
# %DISTMESHND N-D Mesh Generator using Distance Functions.
    dim = ncol(as.matrix(box))
    L0mult = 1 + 0.4/2^(dim - 1)
    rownorm2 = function(x) drop(sqrt((x^2) %*% rep(1, ncol(x))))

    # %1. Create initial distribution in bounding box
    if (dim == 1) {
        p = seq(box[1], box[2], by = h)
    }
    else {
        cbox = lapply(1:dim, function(ii) seq(box[1, ii], box[2,
            ii], by = h))
        p = do.call("expand.grid", cbox)
        p = as.matrix(p)
    }

    # %2. Remove points outside the region, apply the rejection method
    p = p[fdist(p, ...) < geps, ]
    r0 = fh(p, ...)
    p = rbind(pfix, p[runif(nrow(p)) < min(r0)^dim/r0^dim, ])
    N = nrow(p)
    if (N <= dim + 1)
        stop("Not enough starting points inside boundary (is h0 too large?).")
    on.exit(return(invisible(p)))

    cat("Press esc if the mesh seems fine but the algorithm hasn't converged.\n")
    flush.console()
    count = 0

    p0 = 1/.Machine$double.eps

    # mimick Matlab call ``localpairs=nchoosek(1:dim+1,2)'':
    localpairs = as.matrix(expand.grid(1:(dim + 1), 1:(dim + 1)))
    localpairs = localpairs[lower.tri(matrix(TRUE, dim + 1, dim + 1)), 2:1]

    while (TRUE) {
        if (max(rownorm2(p - p0)) > ttol * h) {
            # %3. Retriangulation by Delaunay:

            p0 = p
            t = delaunayn(p)
            pmid = matrix(0, nrow(t), dim)
            for (ii in 1:(dim + 1)) pmid = pmid + p[t[, ii],
                ]/(dim + 1)
            t = t[fdist(pmid, ...) < (-geps), ]
            pair = array(dim = c(0, 2))
            for (ii in 1:nrow(localpairs)) {
                pair = rbind(pair, t[, localpairs[ii, ]])
            }

            # %4. Describe each edge by a unique pair of nodes
            pair = Unique(pair, TRUE); # base-function `unique' is way too slow
            if (dim == 2) {
                trimesh(t, p[, 1:2])
            }
            else if (dim == 3) {
                if (count%%5 == 0) {
                  tetramesh(t, p)
                }
            }
            else {
                cat("Retriangulation #", 15, "\n")
                flush.console()
            }
            count = count + 1
        }
        bars = p[pair[, 1], ] - p[pair[, 2], ]
        L = rownorm2(bars)
        L0 = fh((p[pair[, 1], ] + p[pair[, 2], ])/2, ...)
        L0 = L0 * L0mult * (sum(L^dim)/sum(L0^dim))^(1/dim)
        F = L0 - L
        F[F < 0] = 0
        Fbar = cbind(bars, -bars) * matrix(F/L, nrow = nrow(bars),
            ncol = 2 * dim)
        ii = pair[, t(matrix(1:2, 2, dim))]
        jj = rep(1, nrow(pair)) %o% c(1:dim, 1:dim)
        s = c(Fbar)
        ns = length(s)
        dp = matrix(0, N, dim)
        dp[1:(dim * N)] = rowsum(s, ii[1:ns] + ns * (jj[1:ns] -
            1))
        if (nrow(pfix) > 0)
            dp[1:nrow(pfix), ] = 0
        p = p + deltat * dp
        d = fdist(p, ...)
        ix = d > 0
        gradd = matrix(0, sum(ix), dim)
        for (ii in 1:dim) {
            a = rep(0, dim)
            a[ii] = deps
            d1x = fdist(p[ix, ] + rep(1, sum(ix)) %o% a, ...)
            gradd[, ii] = (d1x - d[ix])/deps
        }
        p[ix, ] = p[ix, ] - (d[ix] %o% rep(1, dim)) * gradd
        maxdp = max(deltat * rownorm2(dp[d < (-geps), ]))
        if (maxdp < ptol * h)
            break
    }
}
"entry.value" <-
function (a, idx) 
{
    if (!is.array(a)) 
        stop(paste("First argument `", deparse(substitute(a)), 
            "' should be an array.", sep = ""))
    if (!is.matrix(idx)) 
        stop(paste("Second argument `", substitute(idx), "' should be a matrix.", 
            sep = ""))
    n <- length(dim(a))
    if (n != ncol(idx)) 
        stop(paste("Number of columns in", deparse(substitute(idx)), 
            "is incompatible is dimension of", deparse(substitute(a))))
    a[(idx - 1) %*% c(1, cumprod(dim(a))[-n]) + 1]
}
"entry.value<-" <-
function (a, idx, value) 
{
    if (!is.array(a)) 
        stop(paste("First argument `", deparse(substitute(a)), 
            "' should be an array.", sep = ""))
    if (!is.matrix(idx)) 
        stop(paste("Second argument `", substitute(idx), "' should be a matrix.", 
            sep = ""))
    n <- length(dim(a))
    if (n != ncol(idx)) 
        stop(paste("Number of columns in", deparse(substitute(idx)), 
            "is incompatible is dimension of", deparse(substitute(a))))
    a[(idx - 1) %*% c(1, cumprod(dim(a))[-n]) + 1] <- value
    return(a)
}

"extprod3d" <-
function (x, y) 
{
    x = matrix(x, ncol = 3)
    y = matrix(y, ncol = 3)
    drop(cbind(x[, 2] * y[, 3] - x[, 3] * y[, 2], x[, 3] * y[, 
        1] - x[, 1] * y[, 3], x[, 1] * y[, 2] - x[, 2] * y[, 
        1]))
}
"matmax" <-
function (...)
{
    x = cbind(...)
    if(!is.numeric(x))
        stop("Input should by numeric.")
    if (!is.matrix(drop(x)))
        x = t(x)
    x[1:nrow(x) + nrow(x) * (max.col(x) - 1)]
}
"matmin" <-
function (...)
{
    x = cbind(...)
    if(!is.numeric(x))
        stop("Input should by numeric.")
    if (!is.matrix(drop(x)))
        x = t(x)
    x[1:nrow(x) + nrow(x) * (max.col(-x) - 1)]
}
"matorder" <-
function (...)
{
    x = cbind(...)
    if(!is.numeric(x))
        stop("Input should by numeric.")
    do.call("order", lapply(1:ncol(x), function(i) x[, i]))
}
"matsort" <-
function (...) {
    x = cbind(...)
    if(!is.numeric(x))
        stop("Input should by numeric.")
    res = array(dim = c(nrow(x), 0))
    if (!is.matrix(drop(x)))
        return(x)
    else if (ncol(x) > 30)
        return(t(apply(x, 1, sort)))
    else while (is.matrix(drop(x))) {
        imc = max.col(x)
        x = t(x)
        imx = nrow(x) * (1:ncol(x) - 1) + imc
        xmax = x[imx]
        x = t(matrix(x[-imx], ncol = ncol(x)))
        res = cbind(res, xmax)
    }
    res = cbind(res, x)
    colnames(res) = NULL
    rownames(res) = NULL
    return(res)
}
"mesh.dcircle" <-
function (p, radius = 1, ...)
{
    if (!is.matrix(p))
        p = t(as.matrix(p))
    sqrt((p^2) %*% c(1, 1))-radius
}
"mesh.diff" <-
function (p, regionA, regionB, ...) 
matmax(regionA(p, ...), -regionB(p, ...))
"mesh.drectangle" <-
function (p, x1 = -1/2, y1 = -1/2, x2 = 1/2, y2 = 1/2, ...) 
{
    if (!is.matrix(p)) 
        p = t(as.matrix(p))
    d1 = y1 - p[, 2]
    d2 = -y2 + p[, 2]
    d3 = x1 - p[, 1]
    d4 = -x2 + p[, 1]
    d5 = sqrt(d1^2 + d3^2)
    d6 = sqrt(d1^2 + d4^2)
    d7 = sqrt(d2^2 + d3^2)
    d8 = sqrt(d2^2 + d4^2)
    matmin = function(...) apply(cbind(...), 1, min)
    d = -matmin(matmin(matmin(-d1, -d2), -d3), -d4)
    ix = d1 > 0 & d3 > 0
    d[ix] = d5[ix]
    ix = d1 > 0 & d4 > 0
    d[ix] = d6[ix]
    ix = d2 > 0 & d3 > 0
    d[ix] = d7[ix]
    ix = d2 > 0 & d4 > 0
    d[ix] = d8[ix]
    d
}
"mesh.dsphere" <-
function (p, radius = 1, ...) 
{
    if (!is.matrix(p)) 
        p = t(as.matrix(p))
    sqrt((p^2) %*% rep(1, ncol(p))) - radius
}
"mesh.hunif" <-
function (p, ...) 
{
    if (!is.matrix(p)) 
        stop("Input `p' should be matrix.")
    rep(1, nrow(p))
}
"mesh.intersect" <-
function (p, regionA, regionB, ...) 
matmax(regionA(p, ...), regionB(p, ...))
"mesh.union" <-
function (p, regionA, regionB, ...) 
matmin(regionA(p, ...), regionB(p, ...))
"paint.fill" <-
function (m, startloc = NULL, plot = TRUE) 
{
    k = list()
    if (plot) {
        k = lapply(1:length(dim(a)), function(i) if (i == 1 || 
            i == 2) 
            1:dim(a)[i]
        else round(dim(a)[i]/2))
        image(1:length(k[[1]]), 1:length(k[[2]]), do.call("[", 
            append(list(m), k)))
    }
    if (is.null(startloc) && plot) {
        startloc = lapply(locator(), round)
        k[[1]] = startloc[[1]]
        k[[2]] = startloc[[2]]
        startloc = k
    }
    if (is.null(startloc)) 
        stop("No starting element provided.")
    if (length(dim(m)) != length(startloc)) 
        stop(paste("Dimension of first argument `", deparse(substitute(m)), 
            "' incompatible with length of \nsecond argument '", 
            deparse(substitute(startloc)), "'", sep = ""))
    cl = do.call("cbind", startloc)
    testvalue = entry.value(m, cl)[1]
    neighb = as.matrix(do.call("expand.grid", lapply(dim(m), 
        function(b) -1:1)))
    checked = array(FALSE, dim = dim(m))
    boundary = array(dim = c(0, length(dim(m))))
    on.exit(return(invisible(list(boundary = boundary, painted = checked))))
    while (nrow(cl) > 0) {
        if (plot) 
            points(cl, pch = 20)

        # generate new elements
        nl = rep(1, nrow(neighb)) %x% cl + neighb %x% rep(1, 
            nrow(cl))
        nl = Unique(nl)
        outside = ((nl < 1) | t(t(nl) > dim(m))) %*% rep(1, ncol(nl))
        nl = nl[!outside, , drop = FALSE]                   # remove invalid elements
        nl = nl[!entry.value(checked, nl), , drop = FALSE]  # remove already checked elements
        if (nrow(nl) <= 0) 
            break

        # check image for boundary at new elements
        is.bound.elem = entry.value(m, nl) != testvalue
        if (any(is.bound.elem)) 
            boundary = rbind(boundary, nl[is.bound.elem, ])

        # make new elements the current elements
        entry.value(checked, nl) = TRUE                     # these elements have been checked
        cl = nl[!is.bound.elem, , drop = FALSE]
    }
}

"surf.tri" <-
function(p,t){
    # original by Per-Olof Persson (c) 2005 for MATLAB
    # ported to R and modified for efficiency by Raoul Grasman (c) 2005

    # construct all faces
    faces = rbind(t[,-4], t[,-3], t[,-2], t[,-1]);
    node4 = rbind(t[, 4], t[, 3], t[, 2], t[, 1]);

#    #original translated from MATLAB:
#    # select the faces that occur only once --> these are the surface boundary faces
#    faces = t(apply(faces,1,sort));                                         # sort each row
#    foo   = apply(faces,1,function(x) do.call("paste",as.list(x,sep=" "))); # makes a string from each row
#    vec   = table(foo);                                                     # tabulates the number of occurences of each string
#    ix    = sapply(names(vec[vec==1]),function(b) which(b==foo))            # obtain indices of faces with single occurence
#    tri   = faces[ix,];
#    node4 = node4[ix];


    # we wish to achieve
    #   > faces = t(apply(faces,1,sort));
    # but this is much too slow, we therefore use max.col and the fact
    # that there are only 3 columns in faces
    i.max = 3*(1:nrow(faces)-1) + max.col(faces)
    i.min = 3*(1:nrow(faces)-1) + max.col(-faces)
    faces = t(faces)
    faces = cbind(faces[i.min], faces[-c(i.max,i.min)], faces[i.max])
    ix = order(faces[,1], faces[,2], faces[,3])

    # Next, we wish to detect duplicated rows in faces, that is,
    #   > qx = duplicated(faces[ix,],MARGIN=1)              # logical indicating duplicates
    # but this is also much to slow, we therefore use the fact that
    # faces[ix,] has the duplicate rows ordered beneath each other
    # and the fact that each row occurs exactly once or twice
    fo = apply(faces[ix,],2,diff)
    dup = (abs(fo) %*% rep(1,3)) == 0        # a row of only zeros indicates duplicate
    dup = c(FALSE,dup)                       # first is never a duplicate
    qx = diff(dup)==0                        # only zero if two consecutive elems are not duplicates
    qx = c(qx, !dup[length(dup)])            # last row is either non-duplicate or should not be selected
    tri = faces[ix[qx],]                     # ix[qx] are indices of singly occuring faces
    node4 = node4[ix[qx]]

    # compute face orientations
    v1 = p[tri[,2],] - p[tri[,1],]; # edge vectors
    v2 = p[tri[,3],] - p[tri[,1],];
    v3 = p[node4,]   - p[tri[,1],];
    ix = which( apply(extprod3d(v1,v2) * v3, 1, sum) > 0 )
    tri[ix,c(2,3)] = tri[ix,c(3,2)]
    rownames(tri) = NULL
    tri
}
"tetramesh" <-
function (T, X, col = heat.colors(nrow(T)), clear = TRUE, ...)
{
    require(rgl)
    if (!is.numeric(T) | !is.numeric(T))
        stop("`T' and `X' should both be numeric.")
    if (ncol(T) != 4)
        stop("Expect first arg `T' to have 4 columns.")
    if (ncol(X) != 3)
        stop("Expect second arg `X' to have 3 columns.")
    t = t(rbind(T[, -1], T[, -2], T[, -3], T[, -4]))
    if (clear)
        rgl.clear()
    rgl.triangles(X[t, 1], X[t, 2], X[t, 3], col = col, ...)
}
trimesh <- function(T, p, p2, add=FALSE, axis=FALSE, boxed=FALSE, ...){
  if(!is.matrix(p)){
     p = cbind(p,p2) # automatically generates error if p2 not present
  }
  xlim = range(p[,1])
  ylim = range(p[,2])
  if(!add){
    plot.new()
    plot.window(xlim, ylim, ...)
  }
  if(boxed){
    box()
  }
  if(axis) {
    axis(1)
    axis(2)
  }
  m = rbind(T[,-1], T[, -2], T[, -3])
  segments(p[m[,1],1],p[m[,1],2],p[m[,2],1],p[m[,2],2], ...)
  return(invisible(list(T = T, p = p)))
}

"trisplinter" <-
function(T,p,threshold=sqrt(.Machine$double.eps)){
  rownorm2 = function(x) drop(sqrt((x^2)%*%c(1,1)))
  d1 = p[T[,1],] - p[T[,2],]
  d2 = p[T[,2],] - p[T[,3],]
  d1 = d1 / rownorm2(d1)
  d2 = d2 / rownorm2(d2)
  ar = d1[,1]*d2[,2] - d1[,2]*d2[,1]
  return(abs(ar) < threshold)
}
