## Copyright 2017 Irene Epifanio

## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.


#removing titles

plotPDMr=function (mean, pc, sd, xl, xu, yl, yu, lineorder) 
{
    for (i in c(-3, 0, 3)) {
        fig <- mean + i * pc * sd
        k <- length(mean)/2
        figx <- fig[1:k]
        figy <- fig[(k + 1):(2 * k)]
      # plot(figx, figy, axes = FALSE, xlab = " ", ylab = " ",   ylim = c(yl, yu), xlim = c(xl, xu))
        plot(figx, figy, axes = TRUE, xlab = " ", ylab = " ",   ylim = c(yl, yu), xlim = c(xl, xu),cex=1.5,cex.axis=1.2)
        lines(figx[lineorder], figy[lineorder])
        if (i == -3) 
        #    title(sub = "mean - c sd")
        if (i == 0) 
        #    title(sub = "mean")
        if (i == 3) 
         #   title(sub = "mean + c sd")
        par(lty = 1)
    }
}

 plotpcar=
function (proc, pcno, type, mag, xl, yl, width, joinline = c(1, 
    1), project = c(1, 2)) 
{
    k <- proc$k
    zero <- matrix(0, k - 1, k)
    h <- defh(k - 1)
    H <- cbind(h, zero)
    H1 <- cbind(zero, h)
    H <- rbind(H, H1)
    if (project[1] == 1) {
        select1 <- 1:k
    }
    if (project[1] == 2) {
        select1 <- (k + 1):(2 * k)
    }
    if (project[1] == 3) {
        select1 <- (2 * k + 1):(3 * k)
    }
    if (project[2] == 1) {
        select2 <- 1:k
    }
    if (project[2] == 2) {
        select2 <- (k + 1):(2 * k)
    }
    if (project[2] == 3) {
        select2 <- (2 * k + 1):(3 * k)
    }
    select <- c(select1, select2)
    meanxy <- c(proc$mshape[, project[1]], proc$mshape[, project[2]])
    if (dim(proc$pcar)[1] == (2 * (k - 1))) {
        pcarot <- (t(H) %*% proc$pcar)[select, ]
    }
    if (dim(proc$pcar)[1] != (2 * (k - 1))) {
        pcarot <- proc$pcar[select, ]
    }
    par(pty = "s")
    par(lty = 1)
    np <- length(pcno)
    nr <- trunc((length(pcno) + 1)/2)
    if (type == "g") {
        par(mfrow = c(nr, 2))
        if (np == 1) {
            par(mfrow = c(1, 1))
        }
        for (i in 1:np) {
            j <- pcno[i]
            fig <- meanxy + pcarot[, j] * 3 * mag * proc$pcasd[j]
            figx <- fig[1:k]
            figy <- fig[(k + 1):(2 * k)]
            YY <- cbind(figx, figy)
            tpsgrid(cbind(proc$mshape[, project[1]], proc$mshape[, 
                project[2]]), YY, xl, yl, width, 1, 0.1, 22)
        }
    }
    else {
        if (type == "r") {
            par(mfrow = c(np, 3))
            op <- par(mfrow = c(np,3),oma = c(5,4,0,0) + 0.1,
          mar = c(1,1,1,1) + 0.1)
            for (i in 1:np) {
                j <- pcno[i]
                plotPDMr(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                  xl, xl + width, yl, yl + width, joinline)
               # title(as.character(paste("PC ", as.character(pcno[i]), ": ", as.character(round(proc$percent[i], 1)), "%")))
            }
            par(op)
        }
        else {
            if (type == "v") {
                par(mfrow = c(nr, 2))
                if (np == 1) {
                  par(mfrow = c(1, 1))
                }
                for (i in 1:np) {
                  j <- pcno[i]
                  plotPDM3(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                    xl, xl + width, yl, yl + width, joinline)
                #  title(as.character(paste("PC ", as.character(pcno[i]), ": ", as.character(round(proc$percent[i],   1)), "%")))
                }
            }
            else {
                if (type == "b") {
                  par(mfrow = c(nr, 2))
                  if (np == 1) {
                    par(mfrow = c(1, 1))
                  }
                  for (i in 1:np) {
                    j <- pcno[i]
                    plotPDMbook(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                      -0.6, 0.6, -0.6, 0.6, joinline)
                    title(as.character(paste("PC ", as.character(pcno[i]), 
                      ": ", as.character(round(proc$percent[i], 
                        1)), "%")))
                  }
                }
                else {
                  if (type == "s") {
                    par(mfrow = c(nr, 2))
                    if (np == 1) {
                      par(mfrow = c(1, 1))
                    }
                    for (i in 1:np) {
                      j <- pcno[i]
                      plotPDM2(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                        xl, xl + width, yl, yl + width, joinline)
                      title(as.character(paste("PC ", as.character(pcno[i]), 
                        ": ", as.character(round(proc$percent[i], 
                          1)), "%")))
                    }
                  }
                  else {
                    if (type == "m") {
                      par(mfrow = c(1, 1))
                      for (i in 1:np) {
                        j <- pcno[i]
                        cat(paste("PC ", pcno[i], " \n"))
                        movie(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                          xl, xl + width, yl, yl + width, joinline, 
                          20)
                      }
                    }
                  }
                }
            }
        }
    }
    par(mfrow = c(1, 1))
}

 shapepcar=
function (proc, pcno = c(1, 2, 3), type = "r", mag = 1, joinline = c(1, 
    1), project = c(1, 2), scores3d = FALSE, color = 2, axes3 = FALSE, 
    rglopen = TRUE, zslice = 0) 
{
    if (scores3d == TRUE) {
        axes3 <- TRUE
        sz <- max(proc$rawscores[, max(pcno)]) - min(proc$rawscores[, 
            max(pcno)])
        spheres3d(proc$rawscores[, pcno], radius = sz/30, col = color)
        if (axes3) {
            axes3d()
        }
    }
    m <- dim(proc$mshape)[2]
    k <- dim(proc$mshape)[1]
    if (scores3d == FALSE) {
        if ((m == 2)) {
            out <- defplotsize2(proc$rotated, project = project)
            xl <- out$xl
            yl <- out$yl
            width <- out$width
            plotpcar(proc, pcno, type, mag, xl, yl, width, joinline, 
                project)
        }
        if ((m == 3) && (type == "m")) {
            for (i in 1:length(pcno)) {
                cat("PC ", pcno[i], " \n")
                plotpca3d(proc, pcno[i])
            }
        }
        h <- defh(k - 1)
        zero <- matrix(0, k - 1, k)
        H <- cbind(h, zero, zero)
        H1 <- cbind(zero, h, zero)
        H2 <- cbind(zero, zero, h)
        H <- rbind(H, H1, H2)
        if (dim(proc$pcar)[1] == (3 * (k - 1))) {
            pcarot <- (t(H) %*% proc$pcar)
            proc$pcar <- pcarot
        }
        if (((m == 3) && (type != "m")) && (type != "g")) {
            if (rglopen) {
                rgl.open()
                rgl.bg(color = "white")
            }
            sz <- centroid.size(proc$mshape)/sqrt(k)/30
            spheres3d(proc$mshape, radius = sz, col = color)
            if (axes3) {
                axes3d()
            }
            for (i in pcno) {
                pc <- proc$mshape + 3 * mag * proc$pcasd[i] * 
                  cbind(proc$pcar[1:k, i], proc$pcar[(k + 1):(2 * 
                    k), i], proc$pcar[(2 * k + 1):(3 * k), i])
                for (j in 1:k) {
                  lines3d(rbind(proc$mshape[j, ], pc[j, ]), col = i)
                }
            }
        }
    }
    if ((m == 3) && (type == "g")) {
        if (rglopen) {
            rgl.open()
            rgl.bg(color = "white")
        }
        for (i in pcno) {
            pc <- proc$mshape + 3 * mag * proc$pcasd[i] * cbind(proc$pcar[1:k, 
                i], proc$pcar[(k + 1):(2 * k), i], proc$pcar[(2 * 
                k + 1):(3 * k), i])
            tpsgrid(proc$mshape, pc, zslice = zslice)
        }
    }
}




 plotpcaro=
function (proc, pcno, type, mag, xl, yl, width, joinline = c(1, 
    1), project = c(1, 2)) 
{
    k <- proc$k
    zero <- matrix(0, k - 1, k)
    h <- defh(k - 1)
    H <- cbind(h, zero)
    H1 <- cbind(zero, h)
    H <- rbind(H, H1)
    if (project[1] == 1) {
        select1 <- 1:k
    }
    if (project[1] == 2) {
        select1 <- (k + 1):(2 * k)
    }
    if (project[1] == 3) {
        select1 <- (2 * k + 1):(3 * k)
    }
    if (project[2] == 1) {
        select2 <- 1:k
    }
    if (project[2] == 2) {
        select2 <- (k + 1):(2 * k)
    }
    if (project[2] == 3) {
        select2 <- (2 * k + 1):(3 * k)
    }
    select <- c(select1, select2)
    meanxy <- c(proc$mshape[, project[1]], proc$mshape[, project[2]])
    if (dim(proc$pcar)[1] == (2 * (k - 1))) {
        pcarot <- (t(H) %*% proc$pcar)[select, ]
    }
    if (dim(proc$pcar)[1] != (2 * (k - 1))) {
        pcarot <- proc$pcar[select, ]
    }
    par(pty = "s")
    par(lty = 1)
    np <- length(pcno)
    nr <- trunc((length(pcno) + 1)/2)
    if (type == "g") {
        par(mfrow = c(nr, 2))
        if (np == 1) {
            par(mfrow = c(1, 1))
        }
        for (i in 1:np) {
            j <- pcno[i]
            fig <- meanxy + pcarot[, j] * 3 * mag * proc$pcasd[j]
            figx <- fig[1:k]
            figy <- fig[(k + 1):(2 * k)]
            YY <- cbind(figx, figy)
            tpsgrid(cbind(proc$mshape[, project[1]], proc$mshape[, 
                project[2]]), YY, xl, yl, width, 1, 0.1, 22)
        }
    }
    else {
        if (type == "r") {
            
            for (i in 1:np) {
                j <- pcno[i]
                plotPDMro(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                  xl, xl + width, yl, yl + width, joinline)
               # title(as.character(paste("PC ", as.character(pcno[i]), ": ", as.character(round(proc$percent[i], 1)), "%")))
            }
            par(op)
        }
        else {
            if (type == "v") {
                par(mfrow = c(nr, 2))
                if (np == 1) {
                  par(mfrow = c(1, 1))
                }
                for (i in 1:np) {
                  j <- pcno[i]
                  plotPDM3(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                    xl, xl + width, yl, yl + width, joinline)
                #  title(as.character(paste("PC ", as.character(pcno[i]), ": ", as.character(round(proc$percent[i],   1)), "%")))
                }
            }
            else {
                if (type == "b") {
                  par(mfrow = c(nr, 2))
                  if (np == 1) {
                    par(mfrow = c(1, 1))
                  }
                  for (i in 1:np) {
                    j <- pcno[i]
                    plotPDMbook(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                      -0.6, 0.6, -0.6, 0.6, joinline)
                    title(as.character(paste("PC ", as.character(pcno[i]), 
                      ": ", as.character(round(proc$percent[i], 
                        1)), "%")))
                  }
                }
                else {
                  if (type == "s") {
                    par(mfrow = c(nr, 2))
                    if (np == 1) {
                      par(mfrow = c(1, 1))
                    }
                    for (i in 1:np) {
                      j <- pcno[i]
                      plotPDM2(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                        xl, xl + width, yl, yl + width, joinline)
                      title(as.character(paste("PC ", as.character(pcno[i]), 
                        ": ", as.character(round(proc$percent[i], 
                          1)), "%")))
                    }
                  }
                  else {
                    if (type == "m") {
                      par(mfrow = c(1, 1))
                      for (i in 1:np) {
                        j <- pcno[i]
                        cat(paste("PC ", pcno[i], " \n"))
                        movie(meanxy, pcarot[, j], mag * proc$pcasd[j], 
                          xl, xl + width, yl, yl + width, joinline, 
                          20)
                      }
                    }
                  }
                }
            }
        }
    }
    par(mfrow = c(1, 1))
}

#overlaid
 shapepcaro=
function (proc, pcno = c(1, 2, 3), type = "r", mag = 1, joinline = c(1, 
    1), project = c(1, 2), scores3d = FALSE, color = 2, axes3 = FALSE, 
    rglopen = TRUE, zslice = 0) 
{
    if (scores3d == TRUE) {
        axes3 <- TRUE
        sz <- max(proc$rawscores[, max(pcno)]) - min(proc$rawscores[, 
            max(pcno)])
        spheres3d(proc$rawscores[, pcno], radius = sz/30, col = color)
        if (axes3) {
            axes3d()
        }
    }
    m <- dim(proc$mshape)[2]
    k <- dim(proc$mshape)[1]
    if (scores3d == FALSE) {
        if ((m == 2)) {
            out <- defplotsize2(proc$rotated, project = project)
            xl <- out$xl
            yl <- out$yl
            width <- out$width
            plotpcaro(proc, pcno, type, mag, xl, yl, width, joinline, 
                project)
        }
        if ((m == 3) && (type == "m")) {
            for (i in 1:length(pcno)) {
                cat("PC ", pcno[i], " \n")
                plotpca3d(proc, pcno[i])
            }
        }
        h <- defh(k - 1)
        zero <- matrix(0, k - 1, k)
        H <- cbind(h, zero, zero)
        H1 <- cbind(zero, h, zero)
        H2 <- cbind(zero, zero, h)
        H <- rbind(H, H1, H2)
        if (dim(proc$pcar)[1] == (3 * (k - 1))) {
            pcarot <- (t(H) %*% proc$pcar)
            proc$pcar <- pcarot
        }
        if (((m == 3) && (type != "m")) && (type != "g")) {
            if (rglopen) {
                rgl.open()
                rgl.bg(color = "white")
            }
            sz <- centroid.size(proc$mshape)/sqrt(k)/30
            spheres3d(proc$mshape, radius = sz, col = color)
            if (axes3) {
                axes3d()
            }
            for (i in pcno) {
                pc <- proc$mshape + 3 * mag * proc$pcasd[i] * 
                  cbind(proc$pcar[1:k, i], proc$pcar[(k + 1):(2 * 
                    k), i], proc$pcar[(2 * k + 1):(3 * k), i])
                for (j in 1:k) {
                  lines3d(rbind(proc$mshape[j, ], pc[j, ]), col = i)
                }
            }
        }
    }
    if ((m == 3) && (type == "g")) {
        if (rglopen) {
            rgl.open()
            rgl.bg(color = "white")
        }
        for (i in pcno) {
            pc <- proc$mshape + 3 * mag * proc$pcasd[i] * cbind(proc$pcar[1:k, 
                i], proc$pcar[(k + 1):(2 * k), i], proc$pcar[(2 * 
                k + 1):(3 * k), i])
            tpsgrid(proc$mshape, pc, zslice = zslice)
        }
    }
}



plotPDMro=function (mean, pc, sd, xl, xu, yl, yu, lineorder) 
{
par(new=FALSE)
    for (i in c(-3, 0, 3)) {
        fig <- mean + i * pc * sd
        k <- length(mean)/2
        figx <- fig[1:k]
        figy <- fig[(k + 1):(2 * k)]
        par(new=TRUE)
       plot(figx, figy, axes = FALSE, xlab = " ", ylab = " ",   ylim = c(yl, yu), xlim = c(xl, xu))
       par(new=TRUE)
       # plot(figx, figy, axes = TRUE, xlab = " ", ylab = " ",   ylim = c(yl, yu), xlim = c(xl, xu))
        lines(figx[lineorder], figy[lineorder])
        if (i == -3) 
        #    title(sub = "mean - c sd")
        if (i == 0) 
        #    title(sub = "mean")
        if (i == 3) 
         #   title(sub = "mean + c sd")
        par(lty = 1)
    }
    par(new=FALSE)
}





plotshapesr=
function (A, B = 0, joinline = c(1, 1), orthproj = c(1, 2), color = 1, 
    symbol = 1,mainm=c()) 
{
    CHECKOK <- TRUE
    if (is.array(A) == FALSE) {
        if (is.matrix(A) == FALSE) {
            cat("Error !! argument should be an array or matrix \n")
            CHECKOK <- FALSE
        }
    }
    if (CHECKOK) {
        k <- dim(A)[1]
        m <- dim(A)[2]
        kk <- k
        if (k >= 15) {
            kk <- 1
        }
        par(pty = "s")
        if (length(c(B)) != 1) {
            par(mfrow = c(1, 2))
        }
        if (length(dim(A)) == 3) {
            A <- A[, orthproj, ]
        }
        if (is.matrix(A) == TRUE) {
            a <- array(0, c(k, 2, 1))
            a[, , 1] <- A[, orthproj]
            A <- a
        }
        out <- defplotsize2(A)
        width <- out$width
        if (length(c(B)) != 1) {
            if (length(dim(B)) == 3) {
                B <- B[, orthproj, ]
            }
            if (is.matrix(B) == TRUE) {
                a <- array(0, c(k, 2, 1))
                a[, , 1] <- B[, orthproj]
                B <- a
            }
            ans <- defplotsize2(B)
            width <- max(out$width, ans$width)
        }
        n <- dim(A)[3]
        lc <- length(color)
        lt <- k * m * n/lc
        color <- rep(color, times = lt)
        lc <- length(symbol)
        lt <- k * m * n/lc
        symbol <- rep(symbol, times = lt)
        plot(A[, , 1], xlim = c(out$xl, out$xl + width), ylim = c(out$yl, 
            out$yl + width), type = "n", xlab = " ", ylab = " ",axes=F,main=mainm)
        for (i in 1:n) {
            select <- ((i - 1) * k * m + 1):(i * k * m)
            points(A[, , i], pch = symbol[select], col = color[select])
            lines(A[joinline, , i])
        }
        if (length(c(B)) != 1) {
            A <- B
            if (is.matrix(A) == TRUE) {
                a <- array(0, c(k, 2, 1))
                a[, , 1] <- A
                A <- a
            }
            out <- defplotsize2(A)
            n <- dim(A)[3]
            plot(A[, , 1], xlim = c(ans$xl, ans$xl + width), 
                ylim = c(ans$yl, ans$yl + width), type = "n", 
                xlab = " ", ylab = " ",axes=F,main=mainm)
            for (i in 1:n) {
                points(A[, , i], pch = symbol[select], col = color[select])
                lines(A[joinline, , i])
            }
        }
    }
}




 plotshapesa3=
function (A, B = 0, joinline = c(1, 1), orthproj = c(1, 2), color = 1, 
    symbol = 1) 
{
    CHECKOK <- TRUE
    if (is.array(A) == FALSE) {
        if (is.matrix(A) == FALSE) {
            cat("Error !! argument should be an array or matrix \n")
            CHECKOK <- FALSE
        }
    }
    if (CHECKOK) {
        k <- dim(A)[1]
        m <- dim(A)[2]
        kk <- k
        if (k >= 15) {
            kk <- 1
        }
        par(pty = "s")
        if (length(c(B)) != 1) {
            par(mfrow = c(1, 2))
        }
        if (length(dim(A)) == 3) {
            A <- A[, orthproj, ]
        }
        if (is.matrix(A) == TRUE) {
            a <- array(0, c(k, 2, 1))
            a[, , 1] <- A[, orthproj]
            A <- a
        }
        out <- defplotsize2(A)
        width <- out$width
        if (length(c(B)) != 1) {
            if (length(dim(B)) == 3) {
                B <- B[, orthproj, ]
            }
            if (is.matrix(B) == TRUE) {
                a <- array(0, c(k, 2, 1))
                a[, , 1] <- B[, orthproj]
                B <- a
            }
            ans <- defplotsize2(B)
            width <- max(out$width, ans$width)
        }
        n <- dim(A)[3]
        lc <- length(color)
        lt <- k * m * n/lc
        color <- rep(color, times = lt)
        lc <- length(symbol)
        lt <- k * m * n/lc
        symbol <- rep(symbol, times = lt)
        plot(A[, , 1], xlim = c(-.6, .4), ylim = c(-.6, 
            .4), type = "n", xlab = " ", ylab = " ")
        for (i in 1:n) {
            select <- ((i - 1) * k * m + 1):(i * k * m)
            points(A[, , i], pch = symbol[select], col = color[select])
            lines(A[joinline, , i])
        }
        if (length(c(B)) != 1) {
            A <- B
            if (is.matrix(A) == TRUE) {
                a <- array(0, c(k, 2, 1))
                a[, , 1] <- A
                A <- a
            }
            out <- defplotsize2(A)
            n <- dim(A)[3]
            plot(A[, , 1], xlim = c(ans$xl, ans$xl + width), 
                ylim = c(ans$yl, ans$yl + width), type = "n", 
                xlab = " ", ylab = " ",cex.axis=1.2)
            for (i in 1:n) {
                points(A[, , i], pch = symbol[select], col = color[select])
                lines(A[joinline, , i])
            }
        }
    }
}