R/Amberg2011.r

In zarquon42b/mesheR: Meshing Operations on Triangular Meshes

#  List all imports necessary in this file
#' @importFrom Matrix sparseMatrix Matrix
#' @importClassesFrom Matrix dgCMatrix sparseMatrix
#' @importMethodsFrom Matrix t
#' @importFrom spam as.spam.dgCMatrix rbind.spam solve.spam diag.spam diag.spam<-
#' @importClassesFrom spam spam spam.chol.NgPeyton
#' @importMethodsFrom spam t chol rbind diag 
NULL

#' subfunction of AmbergDeformSpam
#'
#' subroutine preparing vertex assortment in order to create S using an 'arcface' matrix
#' @param mesh triangular mesh
#' @export buildAffineAmberg

buildAffineAmberg <- function(mesh) 
    {
        verts <- vert2points(mesh)
        norms <- t(facenormals(mesh)$normals[1:3,])##calculate face normals
        allcoo <- rbind(verts,norms)##combine vertices and face normals
        rowind <- (1:ncol(mesh$it))
        rowind3 <- rowind*3
        rowind2 <- rowind3-1
        rowind1 <- rowind3-2
        whichsel <- cbind(rowind1,rowind2,rowind3,t(mesh$it),nrow(verts)+rowind)
        ## setup sparse selection matrix to build Affine per-face trafo-matrix 
        sel <- Matrix(0,3*ncol(mesh$it),nrow(allcoo))
        sel[whichsel[,c(1,4)]] <- -1
        sel[whichsel[,c(1,5)]] <- 1
        sel[whichsel[,c(2,4)]] <- -1
        sel[whichsel[,c(2,6)]] <- 1
        sel[whichsel[,c(3,7)]] <- 1
        return(list(sel=sel,allcoo=allcoo,verts=verts,normals=norms))
    }

#' subfunction of AmbergDeformSpam
#'
#' convert affine per-face trafo-matrix
#' @param mat quadratic matrix
#' @export mat2sparseBlock

mat2sparseBlock <- function (mat) 
{
    mat <- as.matrix(mat)
    nvb <- nrow(mat)
    data <- t(mat)
    data <- as.vector(data)
    j <- matrix(1:(nvb),3,nvb/3)
    j1 <- cbind(j,j,j)*0
    j1[,c(1:ncol(j))*3] <- j
    j1[,(c(1:ncol(j))*3)-1] <- j
    j1[,(c(1:ncol(j))*3)-2] <- j
    j1 <- as.vector(j1)
    i1 <- as.vector(rbind(1:(nvb),1:(nvb),1:(nvb))[,1:nvb])
    out <- sparseMatrix(j = j1, i = i1, x = data)
    invisible((out))
}

#' subfunction of AmbergDeformSpam
#'
#' create per face arcnode matrix to differentiate per-face affine trafo based on per edge approximation
#' @param mesh triangular mesh
#' @export createArcNode
#' 
createArcNode <- function(mesh)
    {
        edges <- vcgNonBorderEdge(mesh, silent=TRUE)##get all non-border edges
        n <- dim(edges)[1]
        dimit <- ncol(mesh$it)
        arcnode <- Matrix(0,n,dimit)
        edges$ind <- 1:n
        ## select adjacent faces later used for approx differentiation
        arcnode[cbind(edges$ind,edges$face1)] <- -1 
        arcnode[cbind(edges$ind,edges$face2)] <- 1
        out <- list()
        out$arcnode <- arcnode
        out$edges <- edges
        invisible(out)
    }

#' subfunction of AmbergDeformSpam
#'
#' calculate area between edge and barycenters of adjacent faces needed in AijkWeights
#' @param x n x 12  matrix containing cornerstones of area 
#' @export areafun
areafun <- function(x)
    {
        out <- as.vector(.Call("areafunCpp",x))
        return(out)
    }

#' subfunction of AmbergDeformSpam
#'
#' calculate per face-deform weights for approx. derivation
#' @param mesh triangular mesh
#' @param arcnode output from createArcNode(mesh)
#' @export AijklWeights
AijklWeights <- function(mesh,arcnode)
    {
        verts <- vert2points(mesh)
        ## get facenormals and barycenters
        bary <-t(facenormals(mesh)$vb[1:3,])
        xijk <- bary[arcnode$edges$face1,]-bary[arcnode$edges$face2,]
        xijk <- apply(xijk,1,function(x){x <- sqrt(sum(x^2))})
        nedge <- length(xijk)
        arraw <- cbind(verts[arcnode$edges$vert1,],verts[arcnode$edges$vert2,],bary[arcnode$edges$face1,],bary[arcnode$edges$face2,])

        aijkl <- sqrt(areafun(arraw))
        diagmat <- aijkl/xijk
        diatmp <- sparseMatrix(i=1:nedge,j=1:nedge,x=diagmat)
        dia3 <- sparseMatrix(i=1:3,j=1:3,x=rep(1,3))
        diatmp <- kronecker(diatmp,dia3)
        return(diatmp)
    }


### 
#' subfunction of AmbergDeformSpam
#'
#' create S - upper block of Jacobian matrix
#' @param mesh triangular mesh
#' @export createS
createS <- function(mesh)
    {
        out <- list()
        sel <- buildAffineAmberg(mesh)
        transel <- sel$sel %*% sel$allcoo
        invtransel <- multisolve3(transel,trans=TRUE)
        arcnodes <- createArcNode(mesh) 
        weight <- AijklWeights(mesh,arcnodes)
        out$sparseAijk <- t(mat2sparseBlock(invtransel))
        out$sel <- sel
        dia3 <- sparseMatrix(i=1:3,j=1:3,x=rep(1,3))
        arcnodes3 <- kronecker(arcnodes$arcnode,dia3)
        out$process <- weight%*%arcnodes3%*%out$sparseAijk%*%sel$sel
        
        return(out)
    }

#' subfunction of AmbergDeformSpam
#'
#' create C - middle block of Jacobian matrix
#' @param lm1 landmarks
#' @param mesh triangular mesh
#' @export createC
createC <- function(lm1,mesh)
    {
        proj <- vcgClost(lm1,mesh,barycentric=TRUE)
        C <- Matrix(0,nrow(lm1),ncol(mesh$vb))
        vertptr <- t(mesh$it[,proj$faceptr])
        faceptr <- cbind(1:nrow(lm1),vertptr)
        for(i in 2:4)
            C[faceptr[,c(1,i)]] <- proj$barycoords[i-1,]
        return(C)
    }

#' subfunction of AmbergDeformSpam
#'
#' function to setup Jc Jacobian submatrix
#' 
#' @param lm1 landmarks
#' @param ncol column numbers of complete Jacbian
#' @param mesh triangular mesh
#' @export createJc

createJc <- function(lm1,ncol,mesh)
    {
        C <- createC(lm1,mesh)
        Jc <- Matrix(0,nrow(lm1),ncol-ncol(C))
        Jc <- cbind(C,Jc)
        Jc <- as.spam.dgCMatrix(Jc)
        return(Jc)
    }

### deformation of a mesh based on correspondences
### converts all sparse matrices to class spam and let spam handle it (much faster)


#' Deform triangular mesh based on correspondences
#' 
#' Perform smooth deformation of a triangular mesh, minimizing per-face
#' distortions.
#' 
#' Perform smooth deformation of a triangular mesh, minimizing per-face
#' distortions.No loose vertices, edges and degenerated faces are allowed, as
#' they lead to singular equation system.
#' 
#' @param mesh triangular mesh of class "mesh3d". No loose vertices, edges and
#' degenerated faces are allowed.
#' @param lm1 m x 3 matrix containing correspondences on "mesh"
#' @param lm2 m x 3 matrix containing target correspondences
#' @param k0 integer: parameter regularizing face normal distortion.
#' @param lambda numeric: parameter regularizing faces's distortion.  
#' @param S optional: object from function createS from previous calculation.
#' @param Hchol Cholesky decomposition of Hessian (obtained by a previous run of AmbergDeformSpam), speeds up things significantly.
#' @param clean logical: if TRUE, \code{vcgClean} from package Rvcg is run to remove duplicated and unreferenced vertices from the mesh and preventing segfaults.
#' @param ridgetol in case of a singular matrix, this value will be added to the diagonal (similar to a ridge regression) to avoid singularity.
#' @return
#' \item{mesh}{deformed mesh}
#' \item{Jn}{Jacobi submatrix Jn}
#' \item{Jc}{Jacobi submatrix Jc}
#' \item{J }{Jacobian matrix}
#' \item{H }{Hessian of J, class "spam"}
#' \item{Hchol}{Cholesky decomposition of H; class"spam"}
#' @author Stefan Schlager
#' @seealso \code{\link{gaussMatch}}
#' @references Amberg, B. 2011. Editing faces in videos, University of Basel.
#' 
#' @importFrom Rvcg vcgClean
#' @importFrom spam crossprod.spam
#' @export AmbergDeformSpam
AmbergDeformSpam <- function(mesh,lm1,lm2,k0=1,lambda=1,S=NULL,Hchol=NULL,clean=FALSE,ridgetol=1e-08)
    {
        if (clean) {
            mesh <- vcgClean(mesh,sel=c(0:1),silent=T)
        }
        options(spam.cholsymmetrycheck=FALSE,spam.eps=.Machine$double.eps)
        
        spam::powerboost("on")
        ## assign(".Spam",spamnosym,envir=asNamespace("spam")) #disable symmetry check
        t0 <- Sys.time()
        out <- list()
        if (is.null(S))
            S <- createS(mesh)
        ##convert S to class spam
        spamS <- as.spam.dgCMatrix(S$process)
        fn <- S$sel$normals
        ## setup Jc Jacobian submatrix
        Jc <- createJc(lm1,ncol(spamS),mesh)
        Jc <- lambda*Jc
        ## setup Jacobian J
        J <- rbind.spam(spamS,Jc)

        ## setup Jn Jacobian submatrix
        Jn0 <- sparseMatrix(i=1:ncol(mesh$it),j=1:ncol(mesh$it),x=rep(1,ncol(mesh$it)))
        Jn <- Matrix(0,ncol(mesh$it),ncol(J)-ncol(Jn0))
        Jn <- cbind(Jn,Jn0)
        Jn <- k0*Jn
        Jn <- as.spam.dgCMatrix(Jn)
        
        ## append Jn to Jacobian J
        J <- rbind.spam(J,Jn)
        Jtc <- t(Jc)%*%lm2
        ## calculate Hessian H
        if (is.null(Hchol)) {
            options(spam.cholsymmetrycheck=FALSE,spam.eps=.Machine$double.eps)
            H <- crossprod.spam(J)
            Hchk <- as.logical(length(which(!diff(H@rowpointers) > 0)))
            if (Hchk)
                diag(H) <- diag(H)
            
            ## Cholesky decomposition of Hessian H
            chk <- try(suppressWarnings(Hchol <- spam::chol(H)),silent = T)
            if (inherits(chk,"try-error")) {
                diag(H) <- diag(H)+ridgetol
                Hchol <- suppressWarnings(spam::chol(H))
                warning(paste("singularity avoided by adding",ridgetol," to diagonal, please check result"))
               
            }
        } else {
            H <- NULL
        }
        k <- solve.spam(Hchol,lambda*Jtc)
        v <- S$sel$allcoo
        v[-c(1:dim(vert2points(mesh))[1]),] <- 0
        v <- Matrix(v)
        k_v <- k-v
        Jtnv <- t(Jn)%*%(fn)
        deltav <- solve.spam(Hchol,k0*Jtnv)+k_v
        vert_new <- as.matrix(t((v+deltav)[1:ncol(mesh$vb),]))
        mesh_new <- mesh
        mesh_new$vb[1:3,] <- (vert_new)
        meshout <- vcgUpdateNormals(mesh_new)
        return(list(mesh=meshout,Jn=Jn,Jc=Jc,J=J,H=H,Hchol=Hchol,S=S))
    }