A toolbox providing methods for data-acquisition, visualisation and statistical methods related to Geometric Morphometrics and shape analysis

Description

A toolbox for Morphometric calculations. Including sliding operations for Semilandmarks, importing, exporting and manipulating of 3D-surface meshes and semi-automated placement of surface landmarks.

Details

Note

The pdf-version of Morpho-help can be obtained from CRAN on https://cran.r-project.org/package=Morpho

For more advanced operations on triangular surface meshes, check out my package Rvcg: https://cran.r-project.org/package=Rvcg or the code repository on github https://github.com/zarquon42b/Rvcg

Author(s)

Stefan Schlager zarquon42@gmail.com

Maintainer: Stefan Schlager zarquon42@gmail.com

References

Schlager S. 2013. Soft-tissue reconstruction of the human nose: population differences and sexual dimorphism. PhD thesis, Universitätsbibliothek Freiburg. URL: http://www.freidok.uni-freiburg.de/volltexte/9181/.

calculate common allometric component

Description

calculate common allometric component

Usage

CAC(x, size, groups = NULL, log = FALSE)

Arguments

Value

References

Mitteroecker P, Gunz P, Bernhard M, Schaefer K, Bookstein FL. 2004. Comparison of cranial ontogenetic trajectories among great apes and humans. Journal of Human Evolution 46(6):679-97.

Examples

data(boneData)
proc <- procSym(boneLM)
pop.sex <- name2factor(boneLM,which=3:4)
cac <- CAC(proc$rotated,proc$size,pop.sex)
plot(cac$CACscores,cac$size)#plot scores against Centroid size
cor.test(cac$CACscores,cac$size)#check for correlation
#visualize differences between large and small on the sample's consensus
## Not run: 
large <- restoreShapes(max(cac$CACscores),cac$CAC,proc$mshape)
small <- restoreShapes(min(cac$CACscores),cac$CAC,proc$mshape)
deformGrid3d(small,large,ngrid=0)

## End(Not run)

Canonical Variate Analysis

Description

performs a Canonical Variate Analysis.

Usage

CVA(
  dataarray,
  groups,
  weighting = TRUE,
  tolinv = 1e-10,
  plot = TRUE,
  rounds = 0,
  cv = FALSE,
  p.adjust.method = "none",
  robust = c("classical", "mve", "mcd"),
  prior = NULL,
  ...
)

Arguments

Value

Author(s)

Stefan Schlager

References

Cambell, N. A. & Atchley, W. R.. 1981 The Geometry of Canonical Variate Analysis: Syst. Zool., 30(3), 268-280.

Klingenberg, C. P. & Monteiro, L. R. 2005 Distances and directions in multidimensional shape spaces: implications for morphometric applications. Systematic Biology 54, 678-688.

See Also

groupPCA

Examples

## all examples are kindly provided by Marta Rufino

if (require(shapes)) {
# perform procrustes fit on raw data
alldat<-procSym(abind(gorf.dat,gorm.dat))
# create factors
groups<-as.factor(c(rep("female",30),rep("male",29)))
# perform CVA and test Mahalanobis distance
# between groups with permutation test by 100 rounds)            
cvall<-CVA(alldat$orpdata,groups,rounds=10000)     
## visualize a shape change from score -5 to 5:
cvvis5 <- 5*matrix(cvall$CVvis[,1],nrow(cvall$Grandm),ncol(cvall$Grandm))+cvall$Grandm
cvvisNeg5 <- -5*matrix(cvall$CVvis[,1],nrow(cvall$Grandm),ncol(cvall$Grandm))+cvall$Grandm
plot(cvvis5,asp=1)
points(cvvisNeg5,col=2)
for (i in 1:nrow(cvvisNeg5))
  lines(rbind(cvvis5[i,],cvvisNeg5[i,]))
}
### Morpho CVA
data(iris)
vari <- iris[,1:4]
facto <- iris[,5]

cva.1=CVA(vari, groups=facto)
## get the typicality probabilities and resulting classifications - tagging
## all specimens with a probability of < 0.01 as outliers (assigned to no class)
typprobs <- typprobClass(cva.1$CVscores,groups=facto)
print(typprobs)
## visualize the CV scores by their groups estimated from (cross-validated)
## typicality probabilities:
if (require(car)) {
scatterplot(cva.1$CVscores[,1],cva.1$CVscores[,2],groups=typprobs$groupaffinCV,
                  smooth=FALSE,reg.line=FALSE)
}
# plot the CVA
plot(cva.1$CVscores, col=facto, pch=as.numeric(facto), typ="n",asp=1,
   xlab=paste("1st canonical axis", paste(round(cva.1$Var[1,2],1),"%")),
   ylab=paste("2nd canonical axis", paste(round(cva.1$Var[2,2],1),"%")))
  
  text(cva.1$CVscores, as.character(facto), col=as.numeric(facto), cex=.7)

  # add chull (merge groups)
  for(jj in 1:length(levels(facto))){
        ii=levels(facto)[jj]
    kk=chull(cva.1$CVscores[facto==ii,1:2])
    lines(cva.1$CVscores[facto==ii,1][c(kk, kk[1])],
    cva.1$CVscores[facto==ii,2][c(kk, kk[1])], col=jj)
    }

  # add 80% ellipses
  if (require(car)) {
  for(ii in 1:length(levels(facto))){
    dataEllipse(cva.1$CVscores[facto==levels(facto)[ii],1],
    cva.1$CVscores[facto==levels(facto)[ii],2], 
                    add=TRUE,levels=.80, col=c(1:7)[ii])}
  }
  # histogram per group
  if (require(lattice)) {
  histogram(~cva.1$CVscores[,1]|facto,
  layout=c(1,length(levels(facto))),
          xlab=paste("1st canonical axis", paste(round(cva.1$Var[1,2],1),"%")))
  histogram(~cva.1$CVscores[,2]|facto, layout=c(1,length(levels(facto))),
          xlab=paste("2nd canonical axis", paste(round(cva.1$Var[2,2],1),"%")))
  } 
  # plot Mahalahobis
  dendroS=hclust(cva.1$Dist$GroupdistMaha)
  dendroS$labels=levels(facto)
  par(mar=c(4,4.5,1,1))
  dendroS=as.dendrogram(dendroS)
  plot(dendroS, main='',sub='', xlab="Geographic areas",
          ylab='Mahalahobis distance')

 
   # Variance explained by the canonical roots:
   cva.1$Var
   # or plot it:
   barplot(cva.1$Var[,2])

# another landmark based example in 3D: 
data(boneData)
groups <- name2factor(boneLM,which=3:4)
proc <- procSym(boneLM)
cvall<-CVA(proc$orpdata,groups)    
#' ## visualize a shape change from score -5 to 5:
cvvis5 <- 5*matrix(cvall$CVvis[,1],nrow(cvall$Grandm),ncol(cvall$Grandm))+cvall$Grandm
cvvisNeg5 <- -5*matrix(cvall$CVvis[,1],nrow(cvall$Grandm),ncol(cvall$Grandm))+cvall$Grandm
## Not run: 
#visualize it
deformGrid3d(cvvis5,cvvisNeg5,ngrid = 0)

## End(Not run)

#for using (e.g. the first 5) PCscores, one will do:
cvall <- CVA(proc$PCscores[,1:5],groups)    
#' ## visualize a shape change from score -5 to 5:
cvvis5 <- 5*cvall$CVvis[,1]+cvall$Grandm
cvvisNeg5 <- -5*cvall$CVvis[,1]+cvall$Grandm
cvvis5 <- restoreShapes(cvvis5,proc$PCs[,1:5],proc$mshape)
cvvisNeg5 <- restoreShapes(cvvisNeg5,proc$PCs[,1:5],proc$mshape)
## Not run: 
#visualize it
deformGrid3d(cvvis5,cvvisNeg5,ngrid = 0)

## End(Not run)

Create Matrices necessary for Thin-Plate Spline

Description

Create (Bending Engergy) Matrices necessary for Thin-Plate Spline, and sliding of Semilandmarks

Usage

CreateL(
  matrix,
  lambda = 1e-08,
  output = c("K", "L", "Linv", "Lsubk", "Lsubk3"),
  threads = 1
)

Arguments

Value

depending on the choices in output:

.

Note

This function is not intended to be called directly - except for playing around to grasp the mechansims of the Thin-Plate Spline.

References

Gunz, P., P. Mitteroecker, and F. L. Bookstein. 2005. Semilandmarks in Three Dimensions, in Modern Morphometrics in Physical Anthropology. Edited by D. E. Slice, pp. 73-98. New York: Kluwer Academic/Plenum Publishers.

Bookstein FL. 1989. Principal Warps: Thin-plate splines and the decomposition of deformations. IEEE Transactions on pattern analysis and machine intelligence 11(6).

See Also

tps3d

Examples

data(boneData)
L <- CreateL(boneLM[,,1])
## calculate Bending energy between first and second specimen:
be <- t(boneLM[,,2])%*%L$Lsubk%*%boneLM[,,2]
## calculate Frobenius norm 
sqrt(sum(be^2))
## the amount is dependant on on the squared scaling factor
# scale landmarks by factor 5 and compute bending energy matrix
be2 <- t(boneLM[,,2]*5)%*%L$Lsubk%*%(boneLM[,,2]*5)
sqrt(sum(be2^2)) # exactly 25 times the result from above
## also this value is not symmetric:
L2 <- CreateL(boneLM[,,2])
be3 <- t(boneLM[,,1])%*%L2$Lsubk%*%boneLM[,,1]
sqrt(sum(be3^2))

convert data from LPS to RAS space and back

Description

convert data from LPS to RAS space and back

Usage

LPS2RAS(x)

Arguments

Details

As e.g. the Slicer versions >= 4.11 are using LPS space, it might be needed to transform data like fiducials and surface models from and back to that space.

Value

returns the rotated data

Estimate the shape by averaging the shape of the nearest neighbours.

Description

Estimate the shape of one set of landmarks by averaging the shape of the nearest neighbours obtained by a second set of landmarks. Weights are calculated either form Mahalanobis or Procrustes distances. This can be useful for data with missing landmarks.

Usage

NNshapeReg(
  x,
  y = NULL,
  n = 3,
  mahalanobis = FALSE,
  mc.cores = parallel::detectCores()
)

Arguments

Details

This function calculates weights from one set of shape data and then estimates the shape of another (or same) set of landmarks. CAUTION: landmark data has to be registered beforehand.

Value

matrix or array of estimates.

See Also

proc.weight, fixLMtps

Examples

if (require(shapes)) {
proc <- procSym(gorf.dat)
#use the closest 3 specimen based on the first 4 landmarks
#to estimate the shape
estim <- NNshapeReg(proc$rotated[1:4,,],proc$rotated,n=3,mc.cores=1)
#compare estimation and true config
plot(proc$rotated[,,1],asp=1)
points(estim[,,1],col=2)
}

correlation between a reduced space and the original space

Description

Calculates the correlation between distances in a reduced space and the original space

Usage

PCdist(PCs, PCscores, x = 5, plot.type = "b")

Arguments

Value

a vector of R-squared values between subspace and fullspace distances and a barplot depicting the correlations belonging to the subspace.

Author(s)

Stefan Schlager

Examples

if (require(shapes)) {
a <- procSym(gorf.dat)
PCdist(a$PCs, a$PCscores, x = 2)
}

Workhorse function for procSym, responsible for Procrustes registration

Description

Workhorse function for procSym, responsible for Procrustes registration

Usage

ProcGPA(
  dat.array,
  tol = 1e-05,
  scale = TRUE,
  CSinit = FALSE,
  silent = TRUE,
  weights = NULL,
  centerweight = FALSE,
  reflection = TRUE,
  pcAlign = TRUE
)

Arguments

Value

returns a list with

Author(s)

Stefan Schlager

References

Goodall C. 1991. Procrustes methods in the statistical analysis of shape. Journal of the Royal Statistical Society. Series B. Statistical Methodology 53:285-239.

Dryden IL, Mardia KV. 1998. Statistical shape analysis. John Wiley and sons, Chichester.

See Also

procSym, rotonto

Examples

data(boneData)
proc <- ProcGPA(boneLM, CSinit=TRUE, silent=TRUE)
#now we landmarks 5 - 9 double the weight as  the others
weights <- c(rep(1,4),rep(2,5),1)
proc.wt <- ProcGPA(boneLM, CSinit=TRUE, weights=weights, silent=TRUE)

calulate regression scores for linear model

Description

calulate regression scores for linear model as specified in Drake & Klingenberg(2008)

Usage

RegScore(model, x = NULL)

Arguments

Details

the data are orthogonally projected onto the regression lines associated with each factor.

Value

returns a n x m matrix containing the regression scores for each specimen.

Warning

if model contains factors with more than 2 levels, R calculates one regression line per 2 factors. Check the colnames of the returned matrix to select the appropriate one. See examples for details.

References

Drake, AG. & Klingenberg, CP. The pace of morphological change: historical transformation of skull shape in St Bernard dogs. Proceedings of the Royal Society B: Biological Sciences, The Royal Society, 2008, 275, 71-76.

Examples

model <- lm(as.matrix(iris[,1:3]) ~ iris[,4])
rs <- RegScore(model)
plot(rs,iris[,4])

##now use a random subsample for model fitting
rand <- sample(nrow(iris))
x <- iris[rand[1:100],4]
newmod <- lm(as.matrix(iris[rand[1:100],1:3]) ~ x)
##predict the rest of data and get their regression scores
rsPred <- RegScore(newmod,as.matrix(iris[rand[101:150],1:3]))
plot(rsPred,iris[rand[101:150],4])
## Not run: 
data(boneData)
proc <- procSym(boneLM)
pop.sex <- name2factor(boneLM,which=3:4) # generate a factor with 4 levels
lm.ps.size <- lm(proc$PCscores ~ pop.sex+proc$size)
rs <- RegScore(lm.ps.size)
colnames(rs) # in this case, the last column contains the regression
# scores associated with proc$size
## validate by using a subsample for fitting
rand <- sample(dim(boneLM)[3])
lm.ps.size0 <- lm(proc$PCscores[rand[1:50],] ~ proc$size[rand[1:50]])
rs0 <- RegScore(lm.ps.size0,proc$PCscores[rand[-c(1:50)],] )
plot(rs0,proc$size[rand[-c(1:50)]])

## End(Not run)

align new data to an existing Procrustes registration

Description

align new data to an existing Procrustes registration

Usage

align2procSym(x, newdata, orp = TRUE)

Arguments

Value

an array with data aligned to the mean shape in x (and projected into tangent space)

Note

this will never yield the same result as a pooled Procrustes analysis because the sample mean is iteratively updated and new data would change the mean.

Examples

require(Morpho)
data(boneData)
# run procSym on entire data set
proc <- procSym(boneLM)
# this is the training data
array1 <- boneLM[,,1:60]
newdata <- boneLM[,,61:80]
proc1 <- procSym(array1)
newalign <- align2procSym(proc1,newdata)
## compare alignment for one specimen to Proc. registration using all data
## Not run: 
deformGrid3d(newalign[,,1],proc$orpdata[,,61])

## End(Not run)

calculate angle between two vectors

Description

calculates unsigned angle between two vectors

Usage

angle.calc(x, y)

Arguments

Value

angle between x and y in radians.

Examples

#calculate angle between two centered and
# superimposed landmark configuration
data(boneData)
opa <- rotonto(boneLM[,,1],boneLM[,,2])
angle.calc(opa$X, opa$Y)

Test whether the direction of two vectors is similar

Description

Test whether the direction of two vectors is similar

Usage

angleTest(x, y)

Arguments

Details

Under the assumption of all (normalized) n-vectors being represented by an n-dimensional hypersphere, the probability of the angle between two vectors is <= the measured values can be estimated as the area of a cap defined by that angle and divided by the hypersphere's complete surface area.

Value

a list with

References

S. Li , 2011. Concise Formulas for the Area and Volume of a Hyperspherical Cap. Asian Journal of Mathematics & Statistics, 4: 66-70.

Examples

x <- c(1,0); y <- c(1,1) # for a circle this should give us p = 0.25 as the angle between vectors
##  is pi/4 and for any vector the segment +-pi/4 covers a quarter of the circle 
angleTest(x,y)

Replace ID-strings of data and associated files.

Description

Replace ID-strings with for digits - e.g. for blind observer error testing.

Usage

anonymize(
  data,
  remove,
  path = NULL,
  dest.path = NULL,
  ext = ".ply",
  split = "_",
  levels = TRUE,
  prefix = NULL,
  suffix = NULL,
  sample = TRUE
)

Arguments

Value

Examples

anonymize(iris,remove=1)

apply affine transformation to data

Description

apply affine transformation to data

Usage

applyTransform(x, trafo, ...)

## S3 method for class 'matrix'
applyTransform(x, trafo, inverse = FALSE, threads = 1, ...)

## S3 method for class 'mesh3d'
applyTransform(x, trafo, inverse = FALSE, threads = 1, ...)

## Default S3 method:
applyTransform(x, trafo, inverse = FALSE, threads = 1, ...)

Arguments

Value

the transformed object

See Also

rotonto, link{rotmesh.onto}, tps3d, computeTransform

Examples

data(boneData)
rot <- rotonto(boneLM[,,1],boneLM[,,2])
trafo <- getTrafo4x4(rot)
boneLM2trafo <- applyTransform(boneLM[,,2],trafo)

compute the area of an n-dimensional hypersphere

Description

compute the area of an n-dimensional hypersphere

Usage

areaSphere(n, r = 1)

Arguments

Value

returns the area

Examples

areaSphere(2) #gives us the circumference of a circle of radius 1

compute the area of an n-dimensional hypersphere cap

Description

compute the area of an n-dimensional hypersphere cap

Usage

areaSpherePart(n, phi, r = 1)

Arguments

Value

returns the area of the hypersphere cap

Examples

areaSpherePart(2,pi/2) # covers half the area of a circle

calculate Pseudo-inverse of a Matrix using RcppArmadillo

Description

a simple wrapper to call Armadillo's pinv function

Usage

armaGinv(x, tol = NULL)

Arguments

Value

Pseudo-inverse

Examples

mat <- matrix(rnorm(12),3,4)
pinvmat <- armaGinv(mat)

calculate mean of an array

Description

calculate mean of a 3D-array (e.g. containing landmarks) (fast) using the Armadillo C++ Backend

Usage

arrMean3(arr)

Arguments

Value

matrix of dimensions k x m.

Note

this is the same as apply(arr, 1:2, mean), only faster for large configurations.

Examples

data(boneData)
proc <- ProcGPA(boneLM, silent = TRUE)
mshape <- arrMean3(proc$rotated)

reverts list2array, converting an array to a list of matrices

Description

reverts list2array, converting an array to a list of matrices

Usage

array2list(x)

Arguments

Value

returns a list containing the matrices

Assess differences in amount and direction of asymmetric variation (only object symmetry)

Description

Assess differences in amount and direction of asymmetric variation (only object symmetry)

Usage

asymPermute(x, groups, rounds = 1000, which = NULL)

Arguments

Value

Note

This test is only sensible if between-group differences concerning directional asymmetry have been established (e.g. by applying a MANOVA on the "asymmetric" PCscores (see also procSym) and one wants to test whether these can be attributed to differences in amount and/or direction of asymmetric displacement. Careful interpretation for very small amounts of directional asymmetry is advised. The Null-Hypothesis is that we have the same directional asymmetry in both groups. If you want to test whether the angle between groups is similar, please use angleTest.

See Also

procSym

calculates the barycenters for all faces of a triangular mesh

Description

calculates the barycenters for all faces of a triangular mesh

Usage

barycenter(mesh)

Arguments

Value

k x 3 matrix of barycenters for all k faces of input mesh.

See Also

closemeshKD

Examples

data(nose)
bary <- barycenter(shortnose.mesh)
## Not run: 
require(rgl)
##visualize mesh
wire3d(shortnose.mesh)
# visualize barycenters
points3d(bary, col=2)
## now each triangle is equipped with a point in its barycenter

## End(Not run)

concatenate multiple arrays/matrices

Description

concatenate multiple 3-dimensional arrays and/or 2-dimensional matrices to one big array

Usage

bindArr(..., along = 1, collapse = FALSE)

Arguments

Details

dimnames, if present and if differing between entries, will be concatenated, separated by a "_".

Value

returns array of combined matrices/arrays

See Also

cbind, rbind, array

Examples

A <- matrix(rnorm(18),6,3)
B <- matrix(rnorm(18),6,3)
C <- matrix(rnorm(18),6,3)

#combine to 3D-array
newArr <- bindArr(A,B,C,along=3)
#combine along first dimension
newArr2 <- bindArr(newArr,newArr,along=1)

Landmarks and a triangular mesh

Description

Landmarks on the osseous human nose and a triangular mesh representing this structure.

Format

boneLM: A 10x3x80 array containing 80 sets of 3D-landmarks placed on the human osseous nose.

skull_0144_ch_fe.mesh: The mesh representing the area of the first individual of boneLM

extract information about fixed landmarks, curves and patches from and atlas generated by "landmark"

Description

After exporting the pts file of the atlas from "landmark" and importing it into R via "read.pts" cExtract gets information which rows of the landmark datasets belong to curves or patches.

Usage

cExtract(pts.file)

Arguments

Value

returns a list containing the vectors with the indices of matrix rows belonging to the in "landmark" defined curves, patches and fix landmarks and a matrix containing landmark coordinates.

Author(s)

Stefan Schlager

See Also

read.lmdta ,read.pts

calculate Centroid Size for a landmark configuration

Description

calculate Centroid Size for a landmark configuration

Usage

cSize(x)

Arguments

Value

returns Centroid size

Examples

data(boneData)
cSize(boneLM[,,1])

Visually browse through a sample rendering its landmarks and corresponding surfaces.

Description

Browse through a sample rendering its landmarks and corresponding surfaces. This is handy e.g. to check if the landmark projection using placePatch was successful, and to mark specific specimen.

Usage

checkLM(
  dat.array,
  path = NULL,
  prefix = "",
  suffix = ".ply",
  col = "white",
  pt.size = NULL,
  alpha = 1,
  begin = 1,
  render = c("w", "s"),
  point = c("s", "p"),
  add = FALSE,
  meshlist = NULL,
  Rdata = FALSE,
  atlas = NULL,
  text.lm = FALSE
)

Arguments

Value

returns an invisible vector of indices of marked specimen.

See Also

placePatch, createAtlas, plotAtlas, file2mesh

Examples

data(nose)
###create mesh for longnose
longnose.mesh <- tps3d(shortnose.mesh,shortnose.lm,longnose.lm,threads=1)
### write meshes to disk
save(shortnose.mesh, file="shortnose")
save(longnose.mesh, file="longnose")

## create landmark array
data <- bindArr(shortnose.lm, longnose.lm, along=3)
dimnames(data)[[3]] <- c("shortnose", "longnose")
## Not run: 
checkLM(data, path="./",Rdata=TRUE, suffix="")

## End(Not run)

## now visualize by using an atlas:
atlas <- createAtlas(shortnose.mesh, landmarks =
           shortnose.lm[c(1:5,20:21),],
patch=shortnose.lm[-c(1:5,20:21),])
if (interactive()){
checkLM(data, path="./",Rdata=TRUE, suffix="", atlas=atlas)
}
## remove data from disk
unlink("shortnose")
unlink("longnose")

check for NA values in a matrix (of landmarks)

Description

check for NA values in a matrix (of landmarks)

Usage

checkNA(x)

Arguments

Value

returns a vector with missin landmarks and a vector of length=0 if none are missing

classify specimen based on between-group PCA or CVA or typprobClass

Description

classify specimen based on between-group PCA, CVA or typprobClass

Usage

classify(x, cv = TRUE, ...)

## S3 method for class 'bgPCA'
classify(x, cv = TRUE, newdata = NULL, ...)

## S3 method for class 'CVA'
classify(x, cv = T, newdata = NULL, prior = NULL, ...)

## S3 method for class 'typprob'
classify(x, cv = TRUE, ...)

Arguments

Value

See Also

CVA,groupPCA, typprobClass

Project coordinates onto a target triangular surface mesh.

Description

For a set of 3D-coordinates the closest matches on a target surface are determined and normals at as well as distances to that point are calculated.

Usage

closemeshKD(
  x,
  mesh,
  k = 50,
  sign = FALSE,
  barycoords = FALSE,
  cores = 1,
  method = 0,
  ...
)

Arguments

Details

The search for the clostest point is designed as follows: Calculate the barycenter of each target face. For each coordinate of x, determine the k closest barycenters and calculate the distances to the closest point on these faces.

Value

returns an object of class mesh3d. with:

Author(s)

Stefan Schlager

References

Baerentzen, Jakob Andreas. & Aanaes, H., 2002. Generating Signed Distance Fields From Triangle Meshes. Informatics and Mathematical Modelling.

Moshfeghi M, Ranganath S, Nawyn K. 1994. Three-dimensional elastic matching of volumes IEEE Transactions on Image Processing: A Publication of the IEEE Signal Processing Society 3:128-138.

See Also

ply2mesh

Examples

data(nose)
out <- closemeshKD(longnose.lm,shortnose.mesh,sign=TRUE)
### show distances - they are very small because
###longnose.lm is scaled to unit centroid size.
hist(out$quality)

predefined colors for bone and skin

Description

predefined colors for bone and skin

Details

available colors are:

bone1

bone2

bone3

skin1

skin2

skin3

skin4

Compute area enclosed within an irregular polygon

Description

Compute area enclosed within an irregular polygon - i.e. defined by curves

Usage

computeArea(x)

Arguments

Details

For 3D coordinates, a PCA is computed and only the first two PCs are used to compute the area. This is a projection of the coordinates onto a 2D plane spanned by those PCs.

Value

returns a list containing

Note

in case custom planes are preferred, the data can first be projected onto such a custom defined plane via points2plane first.

Examples

require(shapes)
require(sf)
myarea <- computeArea(gorf.dat[c(1,6:8,2:5),,1])
myarea$area
plot(myarea$poly)


## 3D example
data(boneData)
myarea3D <- computeArea(boneLM[c(4,2,3,7,5,6,8),,1])
plot(myarea3D$poly)
cent <- colMeans(myarea3D$xpro2D)
text(cent[1],cent[2],labels=paste0("Area=",round(myarea3D$area,digits=2)))

calculate an affine transformation matrix

Description

calculate an affine transformation matrix

Usage

computeTransform(
  x,
  y,
  type = c("rigid", "similarity", "affine", "tps"),
  reflection = FALSE,
  lambda = 1e-08,
  weights = NULL,
  centerweight = FALSE,
  threads = 1
)

Arguments

Details

x and y can also be a pair of meshes with corresponding vertices.

Value

returns a 4x4 (3x3 in 2D case) transformation matrix or an object of class "tpsCoeff" in case of type="tps".

Note

all lines containing NA, or NaN are ignored in computing the transformation.

See Also

rotonto, link{rotmesh.onto}, tps3d

Examples

data(boneData)
trafo <- computeTransform(boneLM[,,1],boneLM[,,2])
transLM <- applyTransform(boneLM[,,2],trafo)

calculates distances and PC-coordinates of covariance matrices

Description

calculates PC-coordinates of covariance matrices by using the Riemannian metric in their respective space.

Usage

covDist(s1, s2)

covPCA(
  data,
  groups,
  rounds = 1000,
  bootrounds = 0,
  lower.bound = 0.05,
  upper.bound = 0.95
)

Arguments

Details

covDist calculates the Distance between covariance matrices while covPCA uses a MDS (multidimensional scaling) approach to obtain PC-coordinates from a distance matrix derived from multiple groups. P-values for pairwise distances can be computed by permuting group membership and comparing actual distances to those obtained from random resampling. To calculate confidence intervals for PC-scores, within-group bootstrapping can be performed.

Value

covDist returns the distance between s1 and s2

covPCA returns a list containing:

if scores = TRUE

if rounds > 0

if bootrounds > 0

Author(s)

Stefan Schlager

References

Mitteroecker P, Bookstein F. 2009. The ontogenetic trajectory of the phenotypic covariance matrix, with examples from craniofacial shape in rats and humans. Evolution 63:727-737.

Hastie T, Tibshirani R, Friedman JJH. 2013. The elements of statistical learning. Springer New York.

See Also

prcomp

Examples

cpca <- covPCA(iris[,1:4],iris[,5])
cpca$p.matrix #show pairwise p-values for equal covariance matrices
## Not run: 
require(car)
sp(cpca$PCscores[,1],cpca$PCscores[,2],groups=levels(iris[,5]),
   smooth=FALSE,xlim=range(cpca$PCscores),ylim=range(cpca$PCscores))

data(boneData)
proc <- procSym(boneLM)
pop <- name2factor(boneLM, which=3)
## compare covariance matrices for PCscores of Procrustes fitted data
cpca1 <- covPCA(proc$PCscores, groups=pop, rounds = 1000)
## view p-values:
cpca1$p.matrix # differences between covariance matrices
# are significant
## visualize covariance ellipses of first 5 PCs of shape
spm(proc$PCscores[,1:5], groups=pop, smooth=FALSE,ellipse=TRUE, by.groups=TRUE)
## covariance seems to differ between 1st and 5th PC
## for demonstration purposes, try only first 4 PCs
cpca2 <- covPCA(proc$PCscores[,1:4], groups=pop, rounds = 1000)
## view p-values:
cpca2$p.matrix # significance is gone

## End(Not run)

#do some bootstrapping 1000 rounds
cpca <- covPCA(iris[,1:4],iris[,5],rounds=0, bootrounds=1000)
#plot bootstrapped data of PC1 and PC2 for first group
plot(t(cpca$boot.data[1,1:2,]),xlim=range(cpca$boot.data[,1,]),
                               ylim=range(cpca$boot.data[,2,]))
points(t(cpca$PCscores[1,]),col="white",pch=8,cex=1.5)##plot actual values
                      
for (i in 2:3) {
  points(t(cpca$boot.data[i,1:2,]),col=i)##plot other groups
  points(t(cpca$PCscores[i,]),col=1,pch=8,cex=1.5)##plot actual values
}

calculate the pooled within groups covariance matrix

Description

calculate the pooled within groups covariance matrix

Usage

covW(data, groups, robust = c("classical", "mve", "mcd"), ...)

Arguments

Value

Returns the pooled within group covariance matrix. The attributes contain the entry means, containing the respective group means.

Author(s)

Stefan Schlager

See Also

cov, typprobClass

Examples

data(iris)
poolCov <- covW(iris[,1:4],iris[,5])

Create an atlas needed in placePatch

Description

Create an atlas needed in placePatch

Usage

createAtlas(
  mesh,
  landmarks,
  patch,
  corrCurves = NULL,
  patchCurves = NULL,
  keep.fix = NULL
)

Arguments

Value

Returns a list of class "atlas". Its content is corresponding to argument names.

Note

This is a helper function of placePatch.

See Also

placePatch, plotAtlas

Examples

data(nose)
atlas <- createAtlas(shortnose.mesh, landmarks =
            shortnose.lm[c(1:5,20:21),], patch=shortnose.lm[-c(1:5,20:21),])

create a list with empty entries to be used as missingList in slider3d

Description

create a list with empty entries to be used as missingList in slider3d

Usage

createMissingList(x)

Arguments

Value

returns a list of length x filled with numerics of length zero.

See Also

fixLMtps,fixLMmirror, slider3d

Examples

## Assume in a sample of 10, the 9th individual has (semi-)landmarks 10:50
#   hanging in thin air (e.g. estimated using fixLMtps)
#   while the others are complete.
## create empty list
missingList <- createMissingList(10)
missingList[[9]] <- 10:50

calculate the orthogonal complement of a 3D-vector

Description

calculate the orthogonal complement of a 3D-vector

Usage

crossProduct(x, y, normalize = TRUE)

tangentPlane(x)

Arguments

Details

calculate the orthogonal complement of a 3D-vector or the 3D-crossproduct, finding an orthogonal vector to a plane in 3D.

Value

tangentPlane:

crossProduct: returns a vector of length 3.

Author(s)

Stefan Schlager

Examples

require(rgl)

x <- c(1,0,0)
y <- c(0,1,0)

#example tangentPlane
z <- tangentPlane(x)
#visualize result
## Not run: 
lines3d(rbind(0, x), col=2, lwd=2)
## show complement
lines3d(rbind(z$y, 0, z$z), col=3, lwd=2)

## End(Not run)
# example crossProduct
z <- crossProduct(x, y)
# show x and y
## Not run: 
lines3d(rbind(x, 0, y), col=2, lwd=2)
# show z
lines3d(rbind(0, z), col=3, lwd=2)

## End(Not run)

cut a mesh by a hyperplane and remove parts above/below that plane

Description

cut a mesh by a hyperplane and remove parts above/below that plane

Usage

cutMeshPlane(mesh, v1, v2 = NULL, v3 = NULL, normal = NULL, keep.upper = TRUE)

Arguments

Details

see cutSpace for more details.

Value

mesh with part above/below hyperplane removed

separate a 3D-pointcloud by a hyperplane

Description

separate a 3D-pointcloud by a hyperplane

Usage

cutSpace(pointcloud, v1, v2 = NULL, v3 = NULL, normal = NULL, upper = TRUE)

Arguments

Details

As above and below are specified by the normal calculated from (v2-v1) \times (v3-v1), where \times denotes the vector crossproduct. This means the normal points "upward" when viewed from the positon where v1, v2 and v3 are arranged counter-clockwise. Thus, which side is "up" depends on the ordering of v1, v2 and v3.

Value

logical vector of length n. Reporting for each point if it is above or below the hyperplane

Examples

data(nose)
v1 <- shortnose.lm[1,]
v2 <- shortnose.lm[2,]
v3 <- shortnose.lm[3,]
pointcloud <- vert2points(shortnose.mesh)
upper <- cutSpace(pointcloud, v1, v2, v3)
## Not run: 
require(rgl)
normal <- crossProduct(v2-v1,v3-v1)
zeroPro <- points2plane(rep(0,3),v1,normal)
## get sign of normal displacement from zero
sig <- sign(crossprod(-zeroPro,normal))
d <- sig*norm(zeroPro,"2")
planes3d(normal[1],normal[2],normal[3],d=d)
points3d(pointcloud[upper,])

## End(Not run)

creates 3D shapes from data to be saved as triangular meshes

Description

creates 3D shapes from 3-dimensional data that can be saved as triangular meshes

Usage

data2platonic(
  datamatrix,
  shape = Rvcg::vcgSphere(),
  col = "red",
  scale = FALSE,
  scalefactor = 1
)

Arguments

Value

returns all shapes merged into a single mesh

Examples

mymesh <- data2platonic(iris[iris$Species=="setosa",1:3],scalefactor=0.1)
mymesh <- mergeMeshes(mymesh,data2platonic(iris[iris$Species=="versicolor",1:3],
                      shape=Rvcg::vcgIcosahedron(),scalefactor=0.1,col="green"))
mymesh <- mergeMeshes(mymesh,data2platonic(iris[iris$Species=="virginica",1:3],
                      shape=Rvcg::vcgTetrahedron(),scalefactor=0.1,col="blue"))
## Not run: 
rgl::shade3d(mymesh)
## save to disk
Rvcg::vcgPlyWrite(mymesh,filename="3D_Data.ply")

## End(Not run)

visualise differences between two superimposed sets of 2D landmarks

Description

visualise differences between two superimposed sets of 2D landmarks by deforming a square grid based on a thin-plate spline interpolation

Usage

deformGrid2d(
  matrix,
  tarmatrix,
  ngrid = 0,
  lwd = 1,
  show = c(1:2),
  lines = TRUE,
  lcol = 1,
  lty = 2,
  col1 = 2,
  col2 = 3,
  pcaxis = FALSE,
  add = FALSE,
  wireframe = NULL,
  margin = 0.2,
  gridcol = "grey",
  gridlty = 1,
  cex1 = 1,
  cex2 = 1,
  ...
)

Arguments

Value

if ngrid > 1 the coordinates of the displaced grid knots are returned.

Author(s)

Stefan Schlager

See Also

tps3d

Examples

if (require(shapes)) {
proc <- procSym(gorf.dat)
deformGrid2d(proc$mshape,proc$rotated[,,1],ngrid=5,pch=19)
}

visualise differences between two superimposed sets of 3D landmarks

Description

visualise differences between two superimposed sets of 3D landmarks by deforming a cubic grid based on a thin-plate spline interpolation

Usage

deformGrid3d(
  matrix,
  tarmatrix,
  ngrid = 0,
  align = FALSE,
  lwd = 1,
  showaxis = c(1, 2),
  show = c(1, 2),
  lines = TRUE,
  lcol = 1,
  add = FALSE,
  col1 = 2,
  col2 = 3,
  type = c("s", "p"),
  size = NULL,
  pcaxis = FALSE,
  ask = TRUE,
  margin = 0.2,
  createMesh = FALSE,
  slice1 = NULL,
  slice2 = NULL,
  slice3 = NULL,
  gridcol = 1,
  gridwidth = 1,
  ...
)

Arguments

Value

if createMesh=TRUE, a mesh containing spheres of reference and target as well as the displacement vectors is returned. Otherwise the knots of the displaced grid is returned.

Author(s)

Stefan Schlager

See Also

tps3d

Examples

if (interactive()){
data(nose)
deformGrid3d(shortnose.lm,longnose.lm,ngrid=10)

## select some slices
deformGrid3d(shortnose.lm,longnose.lm,showaxis=1:3,ngrid=10,slice1=2,slice2=5,slice3=7)
}

deprecated functions of Morpho

Description

document deprecated functions

Usage

deform.grid(...)

showPC(scores, PC, mshape)

adnormals(...)

regdist.raw(
  dataarray,
  plot = TRUE,
  main = "",
  rho = "angle",
  dist.mat.out = FALSE,
  ...
)

crossp(...)

tanplan(...)

conv2backf(...)

warp.mesh(
  mesh,
  matr,
  matt,
  lambda = 1e-08,
  updateNormals = TRUE,
  silent = FALSE
)

showPC(scores, PC, mshape)

make a curve equidistant (optionally up/downsampling)

Description

make a curve equidistant (optionally up/downsampling)

Usage

equidistantCurve(
  x,
  n = NULL,
  open = TRUE,
  subsample = 0,
  increment = 2,
  smoothit = 0,
  mesh = NULL,
  iterations = 1
)

Arguments

Details

Equidistancy is reached by iteratively deforming (using TPS) a straight line with equidistantly placed points to the target using control points with the same spacing as the actual curve. To avoid singularity, the straight line containes a small amount of noise, which can (optionally) be accounted for by smoothing the line by its neighbours.

Value

matrix containing equidistantly placed points

Note

if n >> number of original points, the resulting curves can show unwanted distortions.

Examples

## Not run: 
data(nose)
x <- shortnose.lm[c(304:323),]
xsample <- equidistantCurve(x,n=50,iterations=10,increment=2)

require(rgl)
points3d(xsample,size=5)
spheres3d(x,col=2,radius=0.3,alpha=0.5)

## End(Not run)

calculate variance of a distribution stemming from prediction models

Description

calculates a quotient of the overall varriance within a predicted distribution to that from the original one. This function calculates a naive extension of the univariate R^2-value by dividing the variance in the predicted dat by the variance of the original data. No additional adjustments are made!!

Usage

exVar(model, ...)

## S3 method for class 'lm'
exVar(model, ...)

## S3 method for class 'mvr'
exVar(model, ncomp, val = FALSE, ...)

Arguments

Value

returns the quotient.

Note

The result is only!! a rough estimate of the variance explained by a multivariate model. And the result can be misleading - especially when there are many predictor variables involved. If one is interested in the value each factor/covariate explains, we recommend a 50-50 MANOVA perfomed by the R-package "ffmanova", which reports this value factor-wise.

Author(s)

Stefan Schlager

References

Langsrud O, Juergensen K, Ofstad R, Naes T. 2007. Analyzing Designed Experiments with Multiple Responses Journal of Applied Statistics 34:1275-1296.

Examples

lm1 <- lm(as.matrix(iris[,1:4]) ~ iris[,5])
exVar(lm1)

fast kmeans clustering for 2D or 3D point clouds

Description

fast kmeans clustering for 2D or 3D point clouds - with the primary purpose to get a spatially equally distributed samples

Usage

fastKmeans(x, k, iter.max = 10, project = TRUE, threads = 0)

Arguments

Value

returns a list containing

Examples

require(Rvcg)
data(humface)
set.seed(42)
clust <- fastKmeans(humface,k=1000,threads=1)
## Not run: 
require(rgl)

## plot the cluster centers
spheres3d(clust$centers)

## now look at the vertices closest to the centers
wire3d(humface)
spheres3d(vert2points(humface)[clust$selected,],col=2)

## End(Not run)

Import 3D surface mesh files

Description

Import 3D surface mesh files

Usage

file2mesh(filename, clean = TRUE, readcol = FALSE)

obj2mesh(filename, adnormals = TRUE)

ply2mesh(
  filename,
  adnormals = TRUE,
  readnormals = FALSE,
  readcol = FALSE,
  silent = FALSE
)

Arguments

Details

imports 3D mesh files and store them as an R .object of class mesh3d

Value

Examples

data(nose)
mesh2ply(shortnose.mesh)
mesh <- ply2mesh("shortnose.mesh.ply")

mesh2obj(shortnose.mesh)
mesh2 <- obj2mesh("shortnose.mesh.obj")
## cleanup
unlink(c("shortnose.mesh.obj","shortnose.mesh.ply"))

Graphical interface to find outliers and/or to switch mislabeld landmarks

Description

Graphical interface to find outliers and/or to switch mislabeld landmarks

Usage

find.outliers(
  A,
  color = 4,
  lwd = 1,
  lcol = 2,
  mahalanobis = FALSE,
  PCuse = NULL,
  text = TRUE,
  reflection = FALSE
)

Arguments

Details

This function performs a procrustes fit and sorts all specimen according to their distances (either Procrustes or Mahalanobis-distance) to the sample's consensus. It provides visual help for rearranging landmarks and/or excluding outliers.

Value

Author(s)

Stefan Schlager

See Also

typprob,typprobClass

Examples

data(boneData)
## look for outliers using the mahalanobis distance based on the first
# 10 PCscores
# to perform the example below, you need,of course, uncomment the answers
if (interactive()){
outliers <- find.outliers(boneLM, mahalanobis= TRUE, PCuse=10)
# n # everything is fine 
# n # proceed to next
# s # let's switch some landmarks (3 and 4)
# 3
# 4
# n # we are done
# y # yes, because now it is an outlier
# s #enough for now
}

estimate missing landmarks from their bilateral counterparts

Description

estimate missing landmarks from their bilateral counterparts

Usage

fixLMmirror(x, pairedLM, ...)

## S3 method for class 'array'
fixLMmirror(x, pairedLM, ...)

## S3 method for class 'matrix'
fixLMmirror(x, pairedLM, ...)

Arguments

Details

the configurations are mirrored and the relabled version is matched onto the original using a thin-plate spline deformation. The missing landmark is now estimated using its bilateral counterpart. If one side is completely missing, the landmarks will be mirrored and aligned by the unilateral landmarks.

Value

a matrix or array with fixed missing bilateral landmarks.

Note

in case both landmarks of a bilateral pair are missing a message will be issued. As well if there are missing landmarks on the midsaggital plane are detected.

Examples

data(boneData)
left <- c(4,6,8)
## determine corresponding Landmarks on the right side:
# important: keep same order
right <- c(3,5,7)
pairedLM <- cbind(left, right)
exampmat <- boneLM[,,1]
exampmat[4,] <- NA #set 4th landmark to be NA
fixed <- fixLMmirror(exampmat, pairedLM=pairedLM)
## Not run: 
deformGrid3d(fixed, boneLM[,,1],ngrid=0)
## result is a bit off due to actual asymmetry

## End(Not run)
## example with one side completely missing
oneside <- boneLM[,,1]
oneside[pairedLM[,1],] <- NA
onesidefixed <- fixLMmirror(oneside,pairedLM)
## Not run: 
deformGrid3d(onesidefixed, boneLM[,,1],ngrid=0)
## result is a bit off due to actual asymmetry

## End(Not run)

estimate missing landmarks

Description

Missing landmarks are estimated by deforming a sample average or a weighted estimate of the configurations most similar onto the deficient configuration. The deformation is performed by a Thin-plate-spline interpolation calculated by the available landmarks.

Usage

fixLMtps(data, comp = 3, weight = TRUE, weightfun = NULL)

Arguments

Details

This function tries to estimate missing landmark data by mapping weighted averages from complete datasets onto the missing specimen. The weights are the inverted Procrustes (see proc.weight) distances between the 'comp' closest specimen (using the available landmark configuration).

Value

Note

Be aware that these estimates might be grossly wrong when the missing landmark is quite far off the rest of the landmarks (due to the radial basis function used in the Thin-plate spline interpolation.

Author(s)

Stefan Schlager

References

Bookstein FL. 1989. Principal Warps: Thin-plate splines and the decomposition of deformations IEEE Transactions on pattern analysis and machine intelligence 11.

See Also

proc.weight, tps3d

Examples

if (require(shapes)) {
data <- gorf.dat
### set first landmark of first specimen to NA
data[1,,1] <- NA
repair <- fixLMtps(data,comp=5)
### view difference between estimated and actual landmark
plot(repair$out[,,1],asp=1,pch=21,cex=0.7,col=2)#estimated landmark
points(gorf.dat[,,1],col=3,pch=20)#actual landmark
}
## 3D-example:
data(boneData)
data <- boneLM
### set first and 5th landmark of first specimen to NA
data[c(1,5),,1] <- NA
repair <- fixLMtps(data,comp=10)
##  view difference between estimated and actual landmark
## Not run: 
deformGrid3d(repair$out[,,1], boneLM[,,1],ngrid=0)

## End(Not run)

## Now use a gaussian kernel to compute the weights and use all other configs
gaussWeight <- function(r,sigma=0.05) {
   sigma <- 2*sigma^2
   return(exp(-r^2/ sigma))
}
repair <- fixLMtps(data,comp=79,weightfun=gaussWeight)

find indices of faces that contain specified vertices

Description

find indices of faces that contain specified vertices

Usage

getFaces(mesh, index)

Arguments

Value

vector of face indices

get number of meaningful Principal components

Description

get number of meaningful Principal components

Usage

getMeaningfulPCs(values, n, expect = 2, sdev = FALSE)

Arguments

Details

This implements the method suggested by Bookstein (2014, pp. 324), to determine whether a PC is entitled to interpretation. I.e. a PC is regarded meaningful (its direction) if the ratio of this PC and its successor is above a threshold based on a log-likelihood ratio (and dependend on sample size).

Value

References

Bookstein, F. L. Measuring and reasoning: numerical inference in the sciences. Cambridge University Press, 2014

See Also

getPCtol

Examples

data(boneData)
proc <- procSym(boneLM)
getMeaningfulPCs(proc$eigenvalues,n=nrow(proc$PCscores))
## the first 3 PCs are reported as meaningful
## show barplot that seem to fit the bill
barplot(proc$eigenvalues)

Get viewpoints on a sphere around a 3D mesh

Description

Get viewpoints on a sphere around a 3D mesh to be used with virtualMeshScan

Usage

getOuterViewpoints(
  x,
  n,
  inflate = 1.5,
  radius = NULL,
  subdivision = 3,
  PCA = FALSE
)

Arguments

Value

a list containing

Examples

## Not run: 
data(boneData)
vp <- getOuterViewpoints(skull_0144_ch_fe.mesh,n=100)

require(rgl)
shade3d(skull_0144_ch_fe.mesh,col="white")
spheres3d(vp$viewpoints)
wire3d(vp$sphere)

### Fit to principal axes
vppca <- getOuterViewpoints(skull_0144_ch_fe.mesh,n=100,PCA=TRUE,inflate=1.5)

require(rgl)
shade3d(skull_0144_ch_fe.mesh,col="white")
spheres3d(vppca$viewpoints)
wire3d(vppca$sphere)

## End(Not run)

Obtain PC-scores for new landmark data

Description

Obtain PC-scores for new landmark data

Usage

getPCscores(x, PC, mshape)

Arguments

Value

returns a matrix containing the PC scores

See Also

restoreShapes

Examples

## Not run: 
data(boneData)
proc <- procSym(boneLM[,,-c(1:2)])
newdata <- boneLM[,,c(1:2)]
newdataAlign <- align2procSym(proc,newdata)
scores <- getPCscores(newdataAlign,proc$PCs,proc$mshape)

## End(Not run)

determine the minimum ratio for two subsequent eigenvalues to be considered different

Description

determine the minimum ratio for two subsequent eigenvalues to be considered different

Usage

getPCtol(n, expect = 2)

Arguments

Value

returns the minimum ratio between two subsequent subsequent eigenvalues to be considered different.

References

Bookstein, F. L. Measuring and reasoning: numerical inference in the sciences. Cambridge University Press, 2014

See Also

getMeaningfulPCs

Examples

## reproduce the graph from Bookstein (2014, p. 324)
## and then compare it to ratios for values to be considered
## statistically significant
myseq <- seq(from=10,to = 50, by = 2)
myseq <- c(myseq,seq(from=50,to=1000, by =20))
ratios <- getPCtol(myseq)
plot(log(myseq),ratios,cex=0,xaxt = "n",ylim=c(1,5.2))
ticks <- c(10,20,50,100,200,300,400,500,600,700,800,900,1000)
axis(1,at=log(ticks),labels=ticks)
lines(log(myseq),ratios)
abline(v=log(ticks), col="lightgray", lty="dotted")
abline(h=seq(from=1.2,to=5, by = 0.2), col="lightgray", lty="dotted")

## now we raise the bar and compute the ratios for values
## to be beyond the 95th percentile of
## the corresponding chi-square distribution:
ratiosSig <- getPCtol(myseq,expect=qchisq(0.95,df=2))
lines(log(myseq),ratiosSig,col=2)

Get the linear combinations associated with the common shape change in each latent dimension of a pls2B

Description

Get the linear combinations associated with the common shape change in each latent dimension of a pls2B

Usage

getPLSCommonShape(pls)

Arguments

Value

returns a list containing

References

Mitteroecker P, Bookstein F. 2007. The conceptual and statistical relationship between modularity and morphological integration. Systematic Biology 56(5):818-836.

See Also

plsCoVarCommonShape

Examples

data(boneData)
proc <- procSym(boneLM)
pls <- pls2B(proc$orpdata[1:4,,],proc$orpdata[5:10,,])
commShape <- getPLSCommonShape(pls)
## get common shape for first latent dimension at +-2 sd of the scores
## (you can do this much more convenient using \code{\link{plsCoVarCommonShape}}
scores <- c(-2,2) * sd(c(commShape$XscoresScaled[,1],commShape$YscoresScaled[,1]))
pred <- restoreShapes(scores,commShape$shapevectors[,1],matrix(commShape$commoncenter,10,3))
## Not run: 
deformGrid3d(pred[,,1],pred[,,2])

## End(Not run)

compute changes associated with 2-Block PLS-scores

Description

compute changes associated with 2-Block PLS-scores

Usage

getPLSfromScores(pls, x, y)

Arguments

Details

other than predictPLSfromScores, providing Xscores will not compute predictions of y, but the changes in the original data x that is associated with the specific scores

Value

returns data in the original space associated with the specified values.

compute 2-Block PLS scores for new data

Description

compute 2-Block PLS scores for new data from an existing pls2B

Usage

getPLSscores(pls, x, y)

Arguments

Value

returns a vector of pls-scores

Note

either x or y must be missing

See Also

pls2B, predictPLSfromScores,predictPLSfromData

Get a point along a line with a given distance from the start of the line

Description

Get a point along a line with a given distance from the start of the line

Usage

getPointAlongOutline(mat, dist = 15, reverse = FALSE)

Arguments

Value

returns a vector containing the resulting coordinate

try to identify bilateral landmarks and sort them by side

Description

try to identify bilateral landmarks and sort them by side

Usage

getSides(x, tol = 3, pcAlign = TRUE, icpiter = 100, ...)

Arguments

Details

This function mirrors the landmark set and aligns it to the original. Then it tries to find pairs. If you have a sample, run a Procrustes registration first (without scaling to unit centroid size, or you later have to adapt tol - see examples) and then use the mean as it is usually more symmetrical.

Value

returns a list containing

Examples

data(boneData)
proc <- procSym(boneLM,CSinit=FALSE)
mysides <- getSides(proc$mshape)
if (interactive()){
#visualize bilateral landmarks
deformGrid3d(boneLM[mysides$side1,,1],boneLM[mysides$side2,,1])
## visualize unilateral landmarks
rgl::spheres3d(boneLM[mysides$unilat,,1],radius=0.5)
}

get 4x4 Transformation matrix

Description

get 4x4 Transformation matrix

Usage

getTrafo4x4(x)

## S3 method for class 'rotonto'
getTrafo4x4(x)

Arguments

Value

returns a 4x4 transformation matrix

Examples

data(boneData)
rot <- rotonto(boneLM[,,1],boneLM[,,2])
trafo <- getTrafo4x4(rot)

compute a 4x4 Transformation matrix for rotation around an arbitrary axis

Description

compute a 4x4 Transformation matrix for rotation around an arbitrary axis

Usage

getTrafoRotaxis(pt1, pt2, theta)

Arguments

Note

the resulting matrix can be used in applyTransform

find vertices visible from a given viewpoints

Description

find vertices visible from a given viewpoints

Usage

getVisibleVertices(mesh, viewpoints, offset = 0.001, cores = 1)

Arguments

Value

a vector with (1-based) indices of points visible from at least one of the viewpoints

Note

The function tries to filter out all vertices where the line connecting each vertex with the viewpoints intersects with the mesh itself. As, technically speaking this always occurs at a distance of value=0, a mesh with a tiny offset is generated to avoid these false hits.

Examples

SCP1 <- file2mesh(system.file("extdata","SCP1.ply",package="Morpho"))
viewpoints <- read.fcsv(system.file("extdata","SCP1_Endo.fcsv",package="Morpho"))
visivert <- getVisibleVertices(SCP1,viewpoints)

Perform PCA based of the group means' covariance matrix

Description

Calculate covariance matrix of the groupmeans and project all observations into the eigenspace of this covariance matrix. This displays a low dimensional between group structure of a high dimensional problem.

Usage

groupPCA(
  dataarray,
  groups,
  rounds = 10000,
  tol = 1e-10,
  cv = TRUE,
  mc.cores = parallel::detectCores(),
  weighting = TRUE
)

Arguments

Value

Author(s)

Stefan Schlager

References

Mitteroecker P, Bookstein F 2011. Linear Discrimination, Ordination, and the Visualization of Selection Gradients in Modern Morphometrics. Evolutionary Biology 38:100-114.

Boulesteix, A. L. 2005: A note on between-group PCA, International Journal of Pure and Applied Mathematics 19, 359-366.

See Also

CVA

Examples

data(iris)
vari <- iris[,1:4]
facto <- iris[,5]
pca.1 <-groupPCA(vari,groups=facto,rounds=100,mc.cores=1)

### plot scores
if (require(car)) {
scatterplotMatrix(pca.1$Scores,groups=facto, ellipse=TRUE,
        by.groups=TRUE,var.labels=c("PC1","PC2","PC3"))
}
## example with shape data
data(boneData)
proc <- procSym(boneLM)
pop_sex <- name2factor(boneLM, which=3:4)
gpca <- groupPCA(proc$orpdata, groups=pop_sex, rounds=0, mc.cores=2)
## Not run: 
## visualize shape associated with first between group PC
dims <- dim(proc$mshape)
## calculate matrix containing landmarks of grandmean
grandmean <-gpca$Grandmean
## calculate landmarks from first between-group PC
#                   (+2 and -2 standard deviations)
gpcavis2sd<- restoreShapes(c(-2,2)*sd(gpca$Scores[,1]), gpca$groupPCs[,1], grandmean)
deformGrid3d(gpcavis2sd[,,1], gpcavis2sd[,,2], ngrid = 0,size=0.01)
require(rgl)
## visualize grandmean mesh

grandm.mesh <- tps3d(skull_0144_ch_fe.mesh, boneLM[,,1],grandmean,threads=1)
wire3d(grandm.mesh, col="white")
spheres3d(grandmean, radius=0.01)

## End(Not run)

plot histogram for multiple groups.

Description

plot a histogram for multiple groups, each group colored individually

Usage

histGroup(
  data,
  groups,
  main = paste("Histogram of", dataname),
  xlab = dataname,
  ylab,
  col = NULL,
  alpha = 0.5,
  breaks = "Sturges",
  legend = TRUE,
  legend.x = 80,
  legend.y = 80,
  legend.pch = 15,
  freq = TRUE
)

Arguments

Details

Just a wrapper for the function hist from the "graphics" package

Author(s)

Stefan Schlager

See Also

hist

Examples

data(iris)
histGroup(iris$Petal.Length,iris$Species)

match two landmark configurations using iteratively closest point search

Description

match two landmark configurations using iteratively closest point search

Usage

icpmat(
  x,
  y,
  iterations,
  mindist = 1e+15,
  subsample = NULL,
  type = c("rigid", "similarity", "affine"),
  weights = NULL,
  threads = 1,
  centerweight = FALSE
)

Arguments

Value

returns the rotated landmarks

Examples

data(nose)
icp <- icpmat(shortnose.lm,longnose.lm,iterations=10)

## example using weights
## we want to assign high weights to the first three cordinates
icpw <- icpmat(shortnose.lm,longnose.lm,iterations=10,
               weights=c(rep(100,3),rep(1,620)),centerweight = TRUE)
## the RMSE between those four points and the target is now smaller:
require(Rvcg)
RMSE <- sqrt(sum(vcgKDtree(longnose.lm,icp[1:3,],k=1)$distance^2))
RMSEW<- sqrt(sum(vcgKDtree(longnose.lm,icpw[1:3,],k=1)$distance^2))
barplot(c(RMSE,RMSEW),names.arg=c("RMSE weighted","RMSE unweighted"))
## Not run: 
## plot the differences between unweighted and weighted icp
deformGrid3d(icp,icpw)
## plot the first four coordinates from the icps:
spheres3d(icp[1:3,],col="red",radius = 0.5)
spheres3d(icpw[1:3,],col="green",radius = 0.5)
## plot the target
spheres3d(longnose.lm,col="yellow",radius = 0.2)

## End(Not run)
##2D example  using icpmat to determine point correspondences
if (require(shapes)) {
## we scramble rows to show that this is independent of point order
moving <- gorf.dat[sample(1:8),,1]
plot(moving,asp=1) ## starting config
icpgorf <- icpmat(moving,gorf.dat[,,2],iterations = 20)
points(icpgorf,asp = 1,col=2)
points(gorf.dat[,,2],col=3)## target

## get correspondences using nearest neighbour search
index <- mcNNindex(icpgorf,gorf.dat[,,2],k=1,cores=1)
icpsort <- icpgorf[index,]
for (i in 1:8)
lines(rbind(icpsort[i,],gorf.dat[i,,2]))
}

invert faces' orientation of triangular mesh

Description

inverts faces' orientation of triangular mesh and recomputes vertex normals

Usage

invertFaces(mesh)

Arguments

Value

returns resulting mesh

Author(s)

Stefan Schlager

See Also

updateNormals

Examples

data(nose)
## Not run: 
rgl::shade3d(shortnose.mesh,col=3)

## End(Not run)
noseinvert <- invertFaces(shortnose.mesh)
## show normals
## Not run: 
plotNormals(noseinvert,long=0.01)

## End(Not run)

Calculates the Riemannian distance between two superimposed landmark configs.

Description

Calculates the Riemannian distance between two superimposed landmark configs.

Usage

kendalldist(x, y)

Arguments

Value

returns Riemannian distance

Examples

if(require(shapes)) {
OPA <- rotonto(gorf.dat[,,1],gorf.dat[,,2])
kendalldist(OPA$X,OPA$Y)
}

get intersection between a line and a plane

Description

get intersection between a line and a plane

Usage

line2plane(ptLine, ptDir, planePt, planeNorm)

Arguments

Value

hit point

Note

in case you only have three points on a plane (named pt1, pt2, pt3 you can get the plane's normal by calling crossProduct(pt1-pt2,pt1-pt3).

plot lines between landmarks

Description

add lines connecting landmarks to visualise a sort of wireframe

Usage

lineplot(
  x,
  point,
  col = 1,
  lwd = 1,
  line_antialias = FALSE,
  lty = 1,
  add = TRUE
)

Arguments

Note

works with 2D and 3D configurations

Author(s)

Stefan Schlager

See Also

pcaplot3d

Examples

if (require(shapes)) {
##2D example
plot(gorf.dat[,,1],asp=1)
lineplot(gorf.dat[,,1],point=c(1,5:2,8:6,1),col=2)
}
##3D example
## Not run: 
require(rgl)
data(nose)
points3d(shortnose.lm[1:9,])
lineplot(shortnose.lm[1:9,],point=list(c(1,3,2),c(3,4,5),c(8,6,5,7,9)),col=2)

## End(Not run)

converts a list of matrices to an array

Description

converts a list of matrices to an array

Usage

list2array(x)

Arguments

Value

returns an array concatenating all matrices

find nearest neighbours for 2D and 3D point clouds

Description

find nearest neighbours for point clouds using a kd-tree search. This is just a wrapper of the function vcgKDtree from package Rvcg. Wwraps the function vcgKDtree from package 'Rvcg' (for backward compatibility )

Usage

mcNNindex(target, query, cores = parallel::detectCores(), k = k, ...)

Arguments

Value

l x k matrix containing indices of closest points.

See Also

closemeshKD

Examples

require(rgl)
data(nose)
# find closest vertex on surface for each landmark
clost <- mcNNindex(vert2points(shortnose.mesh),shortnose.lm, k=1,
mc.cores=1)
## Not run: 
spheres3d(vert2points(shortnose.mesh)[clost,],col=2,radius=0.3)
spheres3d(shortnose.lm,radius=0.3)
wire3d(shortnose.mesh)

## End(Not run)

merge multiple triangular meshes into a single one

Description

merge multiple triangular meshes into a single one, preserving color and vertex normals.

Usage

mergeMeshes(...)

Arguments

Value

returns the meshes merged into a single one.

See Also

mesh2ply, file2mesh, ply2mesh

Examples

## Not run: 
require(rgl)
data(boneData)
data(nose)
mergedMesh <- mergeMeshes(shortnose.mesh, skull_0144_ch_fe.mesh)
require(rgl)
shade3d(mergedMesh, col=3)

## End(Not run)

convert a colored mesh to greyscale.

Description

convert the colors of a colored mesh to greyscale values

Usage

mesh2grey(mesh)

Arguments

Value

returns a mesh with material$color replaced by greyscale rgb values.

Author(s)

Stefan Schlager

See Also

ply2mesh,file2mesh

export mesh objects to disk

Description

export mesh objects to disk.

Usage

mesh2obj(x, filename = dataname, writeNormals = TRUE)

mesh2ply(x, filename = dataname, col = NULL, writeNormals = FALSE)

Arguments

Details

export an object of class mesh3d or a set of coordinates to a common mesh file.

Note

meshes containing quadrangular faces will be converted to triangular meshes by splitting the faces. Additionally, mesh2obj is now simply a wrapper of Rvcg::vcgObjWrite.

Author(s)

Stefan Schlager

See Also

ply2mesh, quad2trimesh

Examples

require(rgl)
vb <- c(-1.8,-1.8,-1.8,1.0,1.8,-1.8,-1.8,1.0,-1.8,1.8,-1.8,1.0,1.8,
1.8,-1.8,1.0,-1.8,-1.8,1.8,1.0,1.8,
-1.8,1.8,1.0,-1.8,1.8,1.8,1.0,1.8,1.8,1.8,1.0)
it <- c(2,1,3,3,4,2,3,1,5,5,7,3,5,1,2,2,6,5,6,8,7,7,5,6,7,8,4,4,3,7,4,8,6,6,2,4)
vb <- matrix(vb,4,8) ##create vertex matrix
it <- matrix(it,3,12) ## create face matrix
cube<-list(vb=vb,it=it)
class(cube) <- "mesh3d"
## Not run: 
shade3d(cube,col=3) ## view the green cube

## End(Not run)
mesh2ply(cube,filename="cube") # write cube to a file called cube.ply
unlink("cube.ply")

calculates and visualises distances between surface meshes or 3D coordinates and a surface mesh.

Description

calculates and visualises distances between surface meshes or 3D coordinates and a surface mesh.

Usage

meshDist(x, ...)

## S3 method for class 'mesh3d'
meshDist(
  x,
  mesh2 = NULL,
  distvec = NULL,
  from = NULL,
  to = NULL,
  steps = 20,
  ceiling = FALSE,
  rampcolors = colorRamps::blue2green2red(steps - 1),
  NAcol = "white",
  file = "default",
  imagedim = "100x800",
  uprange = 1,
  ray = FALSE,
  raytol = 50,
  raystrict = FALSE,
  save = FALSE,
  plot = TRUE,
  sign = TRUE,
  tol = NULL,
  tolcol = "green",
  displace = FALSE,
  shade = TRUE,
  method = c("vcglib", "morpho"),
  add = FALSE,
  scaleramp = TRUE,
  threads = 1,
  titleplot = "Distance in mm",
  ...
)

## S3 method for class 'matrix'
meshDist(
  x,
  mesh2 = NULL,
  distvec = NULL,
  from = NULL,
  to = NULL,
  steps = 20,
  ceiling = FALSE,
  rampcolors = colorRamps::blue2green2red(steps - 1),
  NAcol = "white",
  uprange = 1,
  plot = TRUE,
  sign = TRUE,
  tol = NULL,
  tolcol = "green",
  type = c("s", "p"),
  radius = NULL,
  displace = FALSE,
  add = FALSE,
  scaleramp = FALSE,
  titleplot = "Distance in mm",
  ...
)

Arguments

Details

calculates the distances from a mesh or a set of 3D coordinates to another at each vertex; either closest point or along the normals

Value

Returns an object of class "meshDist" if the input is a surface mesh and one of class "matrixDist" if input is a matrix containing 3D coordinates.

Author(s)

Stefan Schlager

References

Detection of inside/outside uses the algorithm proposed in:

Baerentzen, Jakob Andreas. & Aanaes, H., 2002. Generating Signed Distance Fields From Triangle Meshes. Informatics and Mathematical Modelling, .

See Also

render.meshDist, , export.meshDist, shade3d

Examples

data(nose)##load data
##warp a mesh onto another landmark configuration:
longnose.mesh <- tps3d(shortnose.mesh, shortnose.lm, longnose.lm,threads=1)
## Not run: 
mD <- meshDist(longnose.mesh, shortnose.mesh)
##now change the color ramp
render(mD,rampcolors = c("white","red"))

## End(Not run)
#use unsigned distances and a ramp from blue to red
#color distances < 0.01 green:
## Not run: 
meshDist(longnose.mesh, shortnose.mesh, rampcolors = c("blue", "red"),sign=FALSE, tol=0.5)

## End(Not run)

get intersections between mesh and a plane

Description

get intersections between mesh and a plane

Usage

meshPlaneIntersect(mesh, v1, v2 = NULL, v3 = NULL, normal = NULL)

Arguments

Value

returns the intersections of edges and the plane

Examples

data(nose)
v1 <- shortnose.lm[1,]
v2 <- shortnose.lm[2,]
v3 <- shortnose.lm[3,]
intersect <- meshPlaneIntersect(shortnose.mesh,v1,v2,v3)
## Not run: 
require(rgl)
wire3d(shortnose.mesh)
spheres3d(shortnose.lm[1:3,],col=2)#the plane
spheres3d(intersect,col=3,radius = 0.2)#intersections

## End(Not run)

calculate the corners of a mesh's bouning box

Description

calculate the corners of a mesh's bouning box

Usage

meshcube(x)

Arguments

Value

returns a 8 x 3 matrix with the coordinates of the corners of the bounding box.

Examples

require(rgl)
data(boneData)
mc <- meshcube(skull_0144_ch_fe.mesh)
## Not run: 
spheres3d(mc)
wire3d(skull_0144_ch_fe.mesh)

## End(Not run)

calculate average edge length of a triangular mesh

Description

calculate average edge length of a triangular mesh, by iterating over all faces.

Usage

meshres(mesh)

Arguments

Value

returns average edge length (a.k.a. mesh resolution)

Author(s)

Stefan Schlager

Examples

data(boneData)
mres <- meshres(skull_0144_ch_fe.mesh)

mirror landmarks or triangular mesh in place

Description

mirror landmarks or triangular mesh in place

Usage

mirror(
  x,
  icpiter = 50,
  subsample = NULL,
  pcAlign = FALSE,
  mirroraxis = 1,
  initPC = TRUE,
  initCenter = TRUE,
  v1 = NULL,
  v2 = NULL,
  v3 = NULL,
  normal = NULL,
  mc.cores = 2
)

## S3 method for class 'matrix'
mirror(
  x,
  icpiter = 50,
  subsample = NULL,
  pcAlign = FALSE,
  mirroraxis = 1,
  initPC = TRUE,
  initCenter = TRUE,
  v1 = NULL,
  v2 = NULL,
  v3 = NULL,
  normal = NULL,
  mc.cores = 2
)

## S3 method for class 'mesh3d'
mirror(
  x,
  icpiter = 50,
  subsample = NULL,
  pcAlign = FALSE,
  mirroraxis = 1,
  initPC = TRUE,
  initCenter = TRUE,
  v1 = NULL,
  v2 = NULL,
  v3 = NULL,
  normal = NULL,
  mc.cores = 2
)

Arguments

Details

reflect a mesh configuration at the plane spanned by its first 2 principal axis, then try to rigidily register the reflected configuration onto the original one using iterative closest point search to establish correspondences. Also, if a reflection plane is defined, pcAlign, initPC, initCenter and mirroraxis will be ignored and the object will be mirrored on the defined plane (and optionally aligned using an ICP approach).

Value

returns the reflected object

Examples

data(boneData)
boneMir <- mirror(boneLM[,,1],icpiter=50,mc.cores=2,mirroraxis=3)

### mirror on 3 midsaggital landmarks and then optimize it with an ICP
boneMirPlane <- mirror(boneLM[,,1],v1=boneLM[1,,1],v2=boneLM[2,,1],v3=boneLM[9,,1])

## 2D Example:
if (require(shapes)) {
gorfMir <- mirror(gorf.dat[,,1],mirroraxis=2,pcAlign=TRUE,icpiter = 0)
plot(gorfMir,asp = 1)
points(gorf.dat[,,1],col=3)
}
## Not run: 
## now mirror a complete mesh
require(rgl)
skullMir <- mirror(skull_0144_ch_fe.mesh,icpiter=10,subsample = 30,
                   mc.cores=2,mirroraxis=3,pcAlign=TRUE)
###compare result to original
wire3d(skull_0144_ch_fe.mesh,col=3)
wire3d(skullMir,col=2)

## End(Not run)

mirror points or mesh on an arbitrary plane

Description

mirror points or mesh on an arbitrary plane

Usage

mirror2plane(x, v1, normal = NULL, v2 = NULL, v3 = NULL)

## S3 method for class 'matrix'
mirror2plane(x, v1, normal = NULL, v2 = NULL, v3 = NULL)

## S3 method for class 'mesh3d'
mirror2plane(x, v1, normal = NULL, v2 = NULL, v3 = NULL)

Arguments

Value

mirrored coordinates mesh

Examples

# mirror mesh on plane spanned by 3 midsagital landmarks
data(boneData)
mirrmesh <- mirror2plane(skull_0144_ch_fe.mesh,v1=boneLM[1,,1],v2=boneLM[9,,1],v3=boneLM[10,,1])

extract data from array names

Description

extract data from array names

Usage

name2factor(x, sep = "_", which, collapse = sep, as.factor = TRUE)

name2num(x, sep = "_", which, collapse = sep, dif = TRUE)

Arguments

Details

extract data from array names and convert to factors or numbers

If an array is used as input, the data info is expected to be in the 3rd dimension, for a matrix, rownames are used.

Value

returns a vector containing factors or numbers

Author(s)

Stefan Schlager

Examples

data <- matrix(rnorm(200),100,2)
id <- paste("id",1:100,sep="")
pop <- c(rep("pop1",50),rep("pop2",50))
sex <- c(rep("male",50),rep("female",50))
age <- floor(rnorm(100,mean=50,sd=10))
rownames(data) <- paste(id,pop,sex,age,sep="_")
infos <- data.frame(pop=name2factor(data,which=2))
infos$age <- name2num(data,which=4)
infos$pop.sex <- name2factor(data,which=2:3)

landmarks and a triangular mesh representing a human nose

Description

triangular mesh representing a human nose and two matrices containing landmark data

Format

shortnose.mesh: A triangular mesh of class 'mesh3d'.

shortnose.lm: matrix containing example landmark data placed on shortnose.mesh.

longnose.lm: matrix containing example landmark data representing a caricaturesquely deformed human nose.

align two 3D-pointclouds/meshes by their principal axes

Description

align two 3D-pointclouds/meshes by their principal axes

Usage

pcAlign(x, y, optim = TRUE, subsample = NULL, iterations = 10, mc.cores = 2)

## S3 method for class 'matrix'
pcAlign(x, y, optim = TRUE, subsample = NULL, iterations = 10, mc.cores = 2)

## S3 method for class 'mesh3d'
pcAlign(x, y, optim = TRUE, subsample = NULL, iterations = 10, mc.cores = 2)

Arguments

Details

x and y will first be centered and aligned by their PC-axes. If optim=TRUE,all possible 8 ordinations of PC-axes will be tested and the one with the smallest RMSE between the transformed version of x and the closest points on y will be used. Then the rotated version of x is translated to the original center of mass of y.

Value

rotated and translated version of x to the center and principal axes of y.

Examples

data(boneData)
blm1 <- pcAlign(boneLM[,,1],boneLM[,,2])
## Not run: 
require(rgl)
spheres3d(boneLM[,,1])#original position
spheres3d(blm1,col=2)#aligned configuration
spheres3d(boneLM[,,2],col=3)#target

## End(Not run)

visualization of shape variation

Description

visualization of shape change

Usage

pcaplot3d(x, ...)

## S3 method for class 'symproc'
pcaplot3d(
  x,
  pcshow = c(1, 2, 3),
  mag = 3,
  color = 4,
  lwd = 1,
  sym = TRUE,
  legend = TRUE,
  type = c("spheres", "points"),
  ...
)

## S3 method for class 'nosymproc'
pcaplot3d(
  x,
  pcshow = c(1, 2, 3),
  mag = 3,
  color = 4,
  lwd = 1,
  legend = TRUE,
  type = c("spheres", "points"),
  ...
)

Arguments

Details

visualization of the shape changes explained by Principal components

Value

returns an invisible array containing the shapes associated with the Principal components selected.

See Also

procSym

Examples

## Not run: 
data(boneData)
proc <- procSym(boneLM)
pcaplot3d(proc,pcshow=1:3,mag=-3)#only one PC available

## End(Not run)

performs permutation testing for group differences.

Description

This function compares the distance between two groupmeans to the distances obtained by random assignment of observations to this groups.

Usage

permudist(
  data,
  groups,
  rounds = 1000,
  which = NULL,
  p.adjust.method = "none",
  median = FALSE
)

Arguments

Value

Examples

data(boneData)
proc <- procSym(boneLM)
groups <- name2factor(boneLM,which=3)
perm <- permudist(proc$PCscores[,1:10], groups=groups, rounds=10000)

## now we concentrate only on sex dimorphism between Europeans
groups <- name2factor(boneLM,which=3:4)
levels(groups)
perm1 <- permudist(proc$PCscores, groups=groups,which=3:4, rounds=10000)

perfom permutation testing on angles and distances between subgroups of two major groups.

Description

perform permutation test on length and angle of the vectors connecting the subgroup means of two groups: e.g. compare if length and angle between sex related differences in two populations differ significantly.

Usage

permuvec(
  data,
  groups,
  subgroups = NULL,
  rounds = 9999,
  scale = TRUE,
  tol = 1e-10,
  mc.cores = parallel::detectCores()
)

Arguments

Details

This function calculates means of all four subgroups and compares the residual vectors of the major grouping variables by angle and distance.

Value

means.

Examples

data(boneData)
proc <- procSym(boneLM)
pop <- name2factor(boneLM,which=3)
sex <- name2factor(boneLM,which=4)
## use non scaled distances by setting \code{scale = FALSE}
## and only use first 10 PCs
perm <- permuvec(proc$PCscores[,1:10], groups=pop, subgroups=sex,
                 scale=FALSE, rounds=100, mc.cores=2)


## visualize if the amount of sexual dimorphism differs between
# (lenghts of vectors connecting population specific sex's averages)
# differs between European and Chines
hist(perm$permudist, xlim=c(0,0.1),main="measured vs. random distances",
     xlab="distances")
points(perm$dist,10,col=2,pch=19)#actual distance
text(perm$dist,15,label=paste("actual distance\n
     (p=",perm$p.dist,")"))
## not significant!!

## visualize if the direction of sexual dimorphism
# (angle between vectors connecting population specific sex's averages)
# differs between European and Chines
hist(perm$permutangles, main="measured vs. random angles",
     xlab="angles")
points(perm$angle,10,col=2,pch=19)#actual distance
text(perm$angle,15,label=paste("actual distance\n
    (p=",perm$p.angle,")"))
## also non-significant

Project semi-landmarks from a predefined atlas onto all specimen in a sample

Description

Project semi-landmarks from a predefined atlas onto all specimen in a sample. Various mechanisms are implemented to avoid errorneous placement on the wrong surface layer (e.g. inside the bone).

Usage

placePatch(
  atlas,
  dat.array,
  path,
  prefix = NULL,
  fileext = ".ply",
  ray = TRUE,
  inflate = NULL,
  tol = inflate,
  relax.patch = TRUE,
  keep.fix = NULL,
  rhotol = NULL,
  silent = FALSE,
  mc.cores = 1
)

Arguments

Details

This function allows the (relatively) easy projection of surface points defined on an atlas onto all surface of a given sample by Thin-Plate Spline deformation and additional mechanisms to avoid distortions. The algorithm can be outlined as followed.

1. relax curves (if specified) against atlas.

2. deform atlas onto targets by TPS based on predefined landmarks (and curves).

3. project coordinates on deformed atlas onto target mesh

4. 'inflate' or 'deflate' configuration along their normals to make sure all coordinates are on the outside/inside

5. Project inflated points back onto surface along these normals.

6. Check if normals are roughly pointing into the same direction as those on the (deformed) atlas.

7. Relax all points against atlas.

8. the predefined coordinates will note change afterwards!

Value

array containing the projected coordinates appended to the data.array specified in the input. In case dat.array is a matrix only a matrix is returned.

Author(s)

Stefan Schlager

References

Schlager S. 2013. Soft-tissue reconstruction of the human nose: population differences and sexual dimorphism. PhD thesis, Universitätsbibliothek Freiburg. URL: http://www.freidok.uni-freiburg.de/volltexte/9181/.

See Also

createAtlas, relaxLM, checkLM, slider3d, tps3d

Examples

## Not run: 
data(nose)
require(rgl)
###create mesh for longnose
longnose.mesh <- tps3d(shortnose.mesh,shortnose.lm,longnose.lm,threads=1)
## create atlas
fix <- c(1:5,20:21)
atlas <- createAtlas(shortnose.mesh, landmarks =
           shortnose.lm[fix,], patch=shortnose.lm[-c(1:5,20:21),])
## view atlas

plotAtlas(atlas)

## create landmark array with only fix landmarks
data <- bindArr(shortnose.lm[fix,], longnose.lm[fix,], along=3)
dimnames(data)[[3]] <- c("shortnose", "longnose")

### write meshes to disk
mesh2ply(shortnose.mesh, filename="shortnose")
mesh2ply(longnose.mesh, filename="longnose")

patched <- placePatch(atlas, data, path="./", inflate=5)
## now browse through placed patches
checkLM(patched, path="./", atlas=atlas)

## same example with only one target specimen
data <- longnose.lm[fix, ]

patched <- placePatch(atlas, data, prefix="longnose", path="./", inflate=5)
wire3d(longnose.mesh,col=3)
spheres3d(patched)

## End(Not run)

plot the result of slider3d

Description

plot the result of slider3d

Usage

## S3 method for class 'slider3d'
plot(
  x,
  cols = 2:4,
  pt.size = NULL,
  point = c("sphere", "point"),
  specimen = 1,
  add = TRUE,
  ...
)

Arguments

visualize an atlas defined by createAtlas

Description

visualize an atlas defined by createAtlas

Usage

plotAtlas(
  atlas,
  pt.size = NULL,
  alpha = 1,
  render = c("w", "s"),
  point = c("s", "p"),
  meshcol = "white",
  add = TRUE,
  legend = TRUE,
  cols = 2:5
)

Arguments

Details

If legend=TRUE, a plot with a legend will open where coloring of the 3D-spheres is specified.

Value

returns invisible vector containing rgl.id of rendered objects.

See Also

placePatch, createAtlas

Examples

data(nose)
atlas <- createAtlas(shortnose.mesh, landmarks =
           shortnose.lm[c(1:5,20:21),], patch=shortnose.lm[-c(1:5,20:21),])
## Not run: 
plotAtlas(atlas)

## End(Not run)

plots the normals of a triangular surface mesh.

Description

visualises the vertex normals of a triangular surface mesh of class mesh3d. If no normals are contained, they are computed.

Usage

plotNormals(x, length = 1, lwd = 1, col = 1, ...)

Arguments

Author(s)

Stefan Schlager

Examples

## Not run: 
require(rgl)
data(nose)
plotNormals(shortnose.mesh,col=4,length=0.01)
shade3d(shortnose.mesh,col=3)

## End(Not run)

Two-Block partial least square regression.

Description

Performs a Two-Block PLS on two sets of data and assesses the significance of each score by permutation testing

Usage

pls2B(
  x,
  y,
  tol = 1e-12,
  same.config = FALSE,
  rounds = 0,
  useCor = FALSE,
  cv = FALSE,
  cvlv = NULL,
  mc.cores = parallel::detectCores(),
  ...
)

Arguments

Details

The Two-Block PLS tries to find those linear combinations in each block maximising the covariance between blocks. The significance of each linear combination is assessed by comparing the singular value to those obtained from permuted blocks. If both blocks contain landmarks superimposed TOGETHER, the option same.config=TRUE requests superimposition of the permuted configurations (i.e. where the the landmarks of block x are replaced by corresponding landmarks of other specimen.

Value

Author(s)

Stefan Schlager

References

Rohlf FJ, Corti M. 2000. Use of two-block partial least-squares to study covariation in shape. Systematic Biology 49:740-753.

See Also

plsCoVar, getPLSfromScores, predictPLSfromScores, getPLSscores, predictPLSfromData,svd , plsCoVarCommonShape, getPLSCommonShape

Examples

if (require(shapes)) {
### very arbitrary test:
### check if first 4 landmarks covaries with the second 4
proc <- procSym(gorf.dat)
## we do only 50 rounds to minimize computation time
## Not run: #same.config takes too long for CRAN check
pls1 <- pls2B(proc$rotated[1:4,,],proc$rotated[5:8,,],
              same.config=TRUE,rounds=50,mc.cores=2)

## End(Not run)
pls1 <- pls2B(proc$rotated[1:4,,],proc$rotated[5:8,,],
              same.config=FALSE,rounds=50,mc.cores=1)
pls1
layout(matrix(1:4,2,2,byrow=TRUE))
for(i in 1:4)
 plot(pls1$Xscores[,i]~pls1$Yscores[,i])


## predict first 4 landmarks from second 4 for first config
layout(1)
predPLS <- predictPLSfromData(pls1,y=proc$rotated[5:8,,1])
## show differences between prediction and original
deformGrid2d(predPLS,proc$rotated[1:4,,1],pch=19)
##plot the complete first config
points(proc$rotated[,,1])

##show effects of first latent variable
plsEffects <- plsCoVar(pls1,i=1)
deformGrid2d(plsEffects$x[,,1],plsEffects$x[,,2])##show on x
deformGrid2d(plsEffects$y[,,1],plsEffects$y[,,2],add=TRUE,pch=19)##show on y

##show effects of 2nd latent variable
plsEffects2 <- plsCoVar(pls1,i=2)
deformGrid2d(plsEffects2$x[,,1],plsEffects2$x[,,2])##show on x
deformGrid2d(plsEffects2$y[,,1],plsEffects2$y[,,2],add=TRUE,pch=19)##show on y
}

Get the shape changes from pls2B associated with each latent variable

Description

Get the shape changes from pls2B associated with each latent variable

Usage

plsCoVar(pls, i, sdx = 3, sdy = 3)

Arguments

Value

.

See Also

pls2B, getPLSfromScores, predictPLSfromScores, getPLSscores, predictPLSfromData,svd, plsCoVarCommonShape

Compute the shape changes along the common axis of deformations

Description

Compute the shape changes between two blocks of 2D or 3D shape coordiantes along the common axis of deformations defined by each dimension of the latent space

Usage

plsCoVarCommonShape(pls, i, sdcommon = 1)

Arguments

Value

returns an k x m x 2 array with the common shape changes associated with +-sdcommon SD of the i-th latent dimension

Note

this give the same results as plsCoVar, however, using common shape vectors as suggested by Mitteroecker and Bookstein (2007)

References

Mitteroecker P, Bookstein F. 2007. The conceptual and statistical relationship between modularity and morphological integration. Systematic Biology 56(5):818-836.

See Also

pls2B, getPLSfromScores, predictPLSfromScores, getPLSscores, predictPLSfromData,svd, plsCoVar, getPLSCommonShape

Examples

data(boneData)
proc <- procSym(boneLM)
pls <- pls2B(proc$orpdata[1:4,,],proc$orpdata[5:10,,])
commShape <- getPLSCommonShape(pls)
## get common shape for first latent dimension at +-2 sd of the scores
pred <- plsCoVarCommonShape(pls,1,2)
## Not run: 
deformGrid3d(pred[,,1],pred[,,2])

## End(Not run)

projects a 3D coordinate orthogonally onto a plane

Description

projects a 3D coordinate orthogonally onto a plane

Usage

points2plane(x, v1, normal = NULL, v2 = NULL, v3 = NULL)

Arguments

Value

projected point

Examples

data(boneData)
##project rhinion onto plane spanned by Nasion and both Nariales
rpro <- points2plane(boneLM[10,,1],v1=boneLM[9,,1],v2=boneLM[3,,1],v3=boneLM[4,,1])

## Not run: 
require(rgl)
#visualize
wire3d(skull_0144_ch_fe.mesh,col="white")
##get plane normal
normal <- crossProduct(boneLM[3,,1]-boneLM[9,,1],boneLM[4,,1]-boneLM[9,,1])
#' ## get plane offset
d <- norm(points2plane(c(0,0,0),v1=boneLM[9,,1],normal=normal),"2")
spheres3d(boneLM[,,1],radius=0.5)
spheres3d(boneLM[c(3,4,9),,1],radius=0.6,col=3)
##original position of Rhinion
spheres3d(boneLM[10,,1],radius=0.6,col=2)
##projected onto plane
spheres3d(rpro,radius=0.9,col=6)
lines3d(rbind(rpro,boneLM[10,,1]),lwd=3)
##plot plane
planes3d(normal[1],normal[2],normal[3],d=d,col=2,alpha=0.5)

##now we project all points onto that plane:
spheres3d(points2plane(boneLM[,,1],v1=boneLM[9,,1],v2=boneLM[3,,1],v3=boneLM[4,,1]),col=3)

## and finally project the vertices of the mesh onto the plane
meshpro <- points2plane(vert2points(skull_0144_ch_fe.mesh),v1=boneLM[9,,1],normal=normal)
points3d(meshpro,col=2)

## End(Not run)

fast Principal Component Analysis (PCA)

Description

fast Principal Component Analysis (PCA)

Usage

prcompfast(x, retx = TRUE, center = TRUE, scale. = FALSE, tol = NULL, ...)

Arguments

Value

prcomp returns a list with class prcomp containing the followin components:

. pcafast <- prcompfast(iris[,1:4]) pcadefault <- prcompfast(iris[,1:4]) ## check if both results are idential (ignoring the sign) all.equal(lapply(pcafast,abs),lapply(pcadefault,abs))

Note

this function returns the same results as prcomp (apart from sign differences) but uses smarter matrix decompositions making it faster for nrow(x) >> ncol(x) and nrow(x) << ncol(x).

Compute CV-scores from new data

Description

Compute CV-scores from new data

Usage

## S3 method for class 'CVA'
predict(object, newdata, ...)

Arguments

Value

returns the CV-scores for new data

Compute between-group-PC scores from new data

Description

Compute between-group-PC scores from new data

Usage

## S3 method for class 'bgPCA'
predict(object, newdata, ...)

Arguments

Value

returns the between-group-PC scores for new data

Examples

data(boneData)

boneLMPart <- boneLM[,,-(1:2)]
procPart <- procSym(boneLMPart)
pop_sex <- name2factor(boneLMPart, which=3:4)
## compute group PCA without first 2 specimens
gpcaPart <- groupPCA(procPart$orpdata, groups=pop_sex, rounds=0, mc.cores=2,cv=FALSE)
## align new data to Procrustes analysis
newdata <- align2procSym(procPart,boneLM[,,1:2])
## get scores for new data
newscores <- predict(gpcaPart,newdata)

predict 2 Block-PLS from new data

Description

predict 2 Block-PLS from new data

Usage

predictPLSfromData(pls, x, y, ncomp = NULL)

Arguments

Value

returns an array/matrix/vector of predictions - depending on input for computing pls

Note

either x or y must be missing

See Also

pls2B, getPLSscores,predictPLSfromScores

Examples

##see examples in pls2B

predict data from 2-Block PLS-scores

Description

predict data from 2-Block PLS-scores

Usage

predictPLSfromScores(pls, x, y)

Arguments

Value

returns an array/matrix of landmarks or original values, depending on input for computing pls

Note

either x or y must be missing. If x-scores are provided, the yscores will be estimated and the predictions calculated.

See Also

pls2B, getPLSscores,predictPLSfromData

predict relative warps for data not included in the training data set

Description

predict relative warps for data not included in the training data set

Usage

predictRelWarps(x, newdata, noalign = FALSE)

Arguments

Details

This function aligns the new data to the mean from x and transforms it into the relative warp space computed from the training data.

Value

returns a list containing

Examples

data(boneData)
set.seed(42)
training <- sample(1:80,size=60)
rW1 <- relWarps(boneLM[,,training], alpha = -1)
## predict scores for the entire sample
predAll <- predictRelWarps(rW1,boneLM)

## now compare the scores predicted scores to the original ones
layout(matrix(1:4,2,2))
for (i in 1:2) {
  plot(rW1$bescores[,i],predAll$bescores[training,i],main=paste("RW",i))
  plot(rW1$uniscores[,i],predAll$uniscores[training,i],main=paste("UC",i))
}

Predict shapes based on linear models calculated from PCscores

Description

Predict shapes based on linear models calculated from PCscores.

Usage

predictShape.lm(fit, datamod, PC, mshape)

Arguments

Details

This function predicts the landmarks based on models calculated from PCscores.

Value

Warning

Make sure that the levels of the variables used in datamod correspond exactly to those used in fit. Otherwise model matrix will be calculated erroneous.

See Also

model.matrix, lm, formula

Examples

data(boneData)
proc <- procSym(boneLM)
pop <- name2factor(boneLM,which=3)
##easy model with only one factor based on the first four PCs
fit <- lm(proc$PCscores[,1:4] ~ pop)
## get shape for Europeans only
datamod <- ~as.factor(levels(pop))[2]
Eu <- predictShape.lm(fit,datamod, proc$PCs[,1:4],proc$mshape)

## get shape for Europeans and Chinese
datamod <- ~as.factor(levels(pop))
pred <- predictShape.lm(fit,datamod, proc$PCs[,1:4],proc$mshape)
## Not run: 
deformGrid3d(pred$predicted[,,1], pred$predicted[,,2], ngrid = 0)

## End(Not run)

## more complicated model

sex <- name2factor(boneLM,which=4)
fit <- lm(proc$PCscores[,1:4] ~ pop*sex)
## predict female for chinese and European
datamod <- ~(as.factor(levels(pop))*rep(as.factor(levels(sex))[1],2))
pred <- predictShape.lm(fit,datamod, proc$PCs[,1:4],proc$mshape)

## predict female and malefor chinese and European
popmod <- factor(c(rep("eu",2),rep("ch",2)))
sexmod <- rep(as.factor(levels(sex)),2)
datamod <- ~(popmod*sexmod)
pred <- predictShape.lm(fit,datamod, proc$PCs[,1:4],proc$mshape)

## add some (randomly generated) numeric covariate
somevalue <- rnorm(80,sd=10)
fit <- lm(proc$PCscores[,1:4] ~ pop+somevalue)
probs <- quantile(somevalue, probs=c(0.05, 0.95))
## make model for European at 5% and 95% quantile
popmod <- rep(factor(levels(pop))[2],2)
datamod <- ~(popmod+probs)
pred <- predictShape.lm(fit,datamod, proc$PCs[,1:4],proc$mshape)

calculate weights inverse to the distances from the specified observation.

Description

for calculation of a shape model by averaging the observations neighbouring the configuration in question, it is necessary to calculate weights by similarity.

Usage

proc.weight(
  data,
  number,
  ref,
  report = TRUE,
  reg = 0,
  log = FALSE,
  mahalanobis = FALSE,
  weightfun = NULL
)

Arguments

Details

distances of zero will get a weight of 1e12 (this is scaled to all weights summing to one), thus weights for observations further away are converging to zero.

Value

Examples

if (require(shapes)) {
proc <- procSym(gorf.dat)
##get weights for the the four specimen closest to the first observation.
weights <- proc.weight(proc$rotated,4,1)

##estimate the first specimen by weighted neighbour shapes.
estim <- proc$mshape*0;
for (i in 1:4)
{estim <-estim+proc$rotated[,,weights$data$nr[i]]*weights$data$weight[i]}

### visualise
plot(estim,asp=1)## show estimation
points(proc$rotated[,,1],col=3)##show original

## use a gaussian smoother to compute weights using a bandwidth of 0.05
gaussWeight <- function(r,sigma=0.05) {
   sigma <- 2*sigma^2
   return(exp(-r^2/ sigma))
}
weights <- proc.weight(proc$rotated,4,1,weightfun=gaussWeight)
}

Procrustes ANOVA for structures with object symmetry

Description

Procrustes ANOVA for structures with object symmetry, currently only supporting the factors 'specimen', 'side' and the interaction term.

Usage

procAOVsym(symproc, indnames = NULL)

Arguments

Details

performs a Procrustes ANOVA for configurations with object symmetry (as described in Klingenberg et al. 2002).

Value

returns a dataframe containing Sums of Squares for each factor.

Note

In future releases the implementation of support for bilateral symmetry and more factors is intended.

Author(s)

Stefan Schlager

References

Klingenberg CP, Barluenga M, Meyer A. 2002. Shape analysis of symmetric structures: quantifying variation among individuals and asymmetry. Evolution 56:1909-20.

See Also

procSym

Examples

data(boneData)
left <- c(4,6,8)
## determine corresponding Landmarks on the right side:
# important: keep same order
right <- c(3,5,7)
pairedLM <- cbind(left,right)
symproc <- procSym(boneLM, pairedLM=pairedLM)
procAOVsym(symproc)

Procrustes registration

Description

procSym performs Procrustes superimposition including sliding of semi-landmarks on curves/outlines in 2D and 3D.

Usage

procSym(
  dataarray,
  scale = TRUE,
  reflect = TRUE,
  CSinit = TRUE,
  orp = TRUE,
  proctol = 1e-05,
  tol = 1e-05,
  pairedLM = NULL,
  sizeshape = FALSE,
  use.lm = NULL,
  center.part = FALSE,
  weights = NULL,
  centerweight = FALSE,
  pcAlign = TRUE,
  distfun = c("angle", "riemann"),
  SMvector = NULL,
  outlines = NULL,
  deselect = FALSE,
  recursive = TRUE,
  iterations = 0,
  initproc = FALSE,
  bending = TRUE,
  stepsize = 1
)

Arguments

Details

This function performs Procrustes registration, allowing a variety of options, including scaling, orthogonal projection into tangentspace and relaxation of semi-landmarks on curves (without reprojection onto the surface/actual outline). It also allows the superimpositioning to be performed using only a subset of the available landmark. For taking into account object symmetry, pairedLM needs to be set. This generates an object of class "symproc". Otherwise an object of class "nosymproc".

Value

Note

For processing of surface landmarks or including the reprojection of slid landmarks back onto 3D-surface representations, the usage of slider3d is recommended.

Author(s)

Stefan Schlager

References

Dryden IL, and Mardia KV. 1998. Statistical shape analysis. Chichester.

Klingenberg CP, Barluenga M, and Meyer A. 2002. Shape analysis of symmetric structures: quantifying variation among individuals and asymmetry. Evolution 56(10):1909-1920.

Gunz, P., P. Mitteroecker, and F. L. Bookstein. 2005. Semilandmarks in Three Dimensions, in Modern Morphometrics in Physical Anthropology. Edited by D. E. Slice, pp. 73-98. New York: Kluwer Academic/Plenum Publishers.

See Also

slider3d

Examples

require(rgl)
data(boneData)

### do an analysis of symmetric landmarks
## visualize landmarks on surface
## Not run: 
 spheres3d(boneLM[,,1])
wire3d(skull_0144_ch_fe.mesh,col=3)
## get landmark numbers
text3d(boneLM[,,1],text=paste(1:10),adj = 1, cex=3)

## End(Not run)
## determine paired Landmarks left side:
left <- c(4,6,8)
## determine corresponding Landmarks on the right side:
# important: keep same order
right <- c(3,5,7)
pairedLM <- cbind(left,right)
symproc <- procSym(boneLM, pairedLM=pairedLM)
## Not run: 
## visualize first 3 PCs of symmetric shape
pcaplot3d(symproc, sym=TRUE)
## visualize first 3 PCs of asymmetric shape
pcaplot3d(symproc, sym=FALSE)

## visualze distribution of symmetric PCscores population
pop <- name2factor(boneLM, which=3)
if (require(car)) {
spm(~symproc$PCscore_sym[,1:5], groups=pop)
## visualze distribution of asymmetric PCscores population
spm(~symproc$PCscore_asym[,1:5], groups=pop)
}

## End(Not run)

Project points onto the closest point on a mesh

Description

project points onto a given surface and return projected points and normals.

Usage

projRead(lm, mesh, readnormals = TRUE, smooth = FALSE, sign = TRUE, ...)

Arguments

Value

if readnormals = FALSE, a m x 3 matrix containing projected points is returned, otherwise a list, where

Author(s)

Stefan Schlager

References

Detection of inside/outside uses the algorithm proposed in:

Baerentzen, Jakob Andreas. & Aanaes, H., 2002. Generating Signed Distance Fields From Triangle Meshes. Informatics and Mathematical Modelling.

See Also

closemeshKD

Examples

data(nose)
## Not run: 
repro <- projRead(shortnose.lm,shortnose.mesh)

## End(Not run)

Q-Q plot to assess normality of data

Description

qqmat plots Mahalanobisdistances of a given sample against those expected from a Gaussian distribution

Usage

qqmat(x, output = FALSE, square = FALSE)

Arguments

Value

if output=TRUE, the following values are returned

Author(s)

Stefan Schlager

See Also

qqplot

Examples

require(MASS)
### create normally distributed data
data <- mvrnorm(100,mu=rep(0,5),Sigma = diag(5:1))
qqmat(data)

###create non normally distributed data
data1 <- rchisq(100,df=3)
qqmat(data1,square=FALSE)

converts a mesh containing quadrangular faces into one only consisting of triangles

Description

converts a mesh containing quadrangular faces into one only consisting of triangles

Usage

quad2trimesh(mesh, updateNormals = TRUE)

Arguments

Value

triangular mesh with updated normals

Examples

Sigma <- diag(3:1) #create a 3D-covariance matrix
require(rgl)
quadmesh <- ellipse3d(Sigma)##create quadmesh
trimesh <- quad2trimesh(quadmesh)# convert to trimesh

Export data to MorphoJ and Morphologika

Description

Export data to MorphoJ and Morphologika

Usage

r2morphoj(x, file, id.string = NULL)

r2morphologika(x, file = file, labels = NULL, labelname = NULL, ...)

Arguments

Details

Export data to MorphoJ and Morphologika

Examples

if (require(shapes)) {
r2morphoj(gorf.dat,file="gorf.dat")

data <- bindArr(gorf.dat, gorm.dat, along=3)
datalabels <- c(rep("female",dim(gorf.dat)[3]),
rep("male",dim(gorm.dat)[3]))
labelname <- "sex"
r2morphologika(data, labels=datalabels, labelname= labelname, file="data.dat")
## cleanup
unlink(c("gorf.dat","data.dat"))
}

projects the vertices of a mesh along its normals onto the surface of another one.

Description

projects the vertices of a mesh onto the surface of another one by searching for the closest point along vertex normals on the target by for each vertex.

Usage

ray2mesh(mesh1, tarmesh, tol = 1e+12, inbound = FALSE, mindist = FALSE, ...)

Arguments

Value

returns projected mesh with additional list entries:

Author(s)

Stefan Schlager

See Also

ply2mesh, closemeshKD

batch import data from files

Description

imports all data files contained in a specified folder.

Usage

read.csv.folder(
  folder,
  x,
  y = 2:4,
  rownames = NULL,
  header = TRUE,
  dec = ".",
  sep = ";",
  pattern = "csv",
  addSpec = NULL,
  back = TRUE
)

Arguments

Value

Author(s)

Stefan Schlager

See Also

read.table

read fiducials from slicer4

Description

read fiducials from slicer4

Usage

read.fcsv(x, na = NULL, lps2ras = FALSE)

Arguments

Value

a k x 3 matrix with landmarks

read dta files

Description

reads .dta files created by the software Landmark http://graphics.idav.ucdavis.edu/research/EvoMorph

Usage

read.lmdta(file = "x", na = 9999)

Arguments

Value

Read saved pick-points from meshlab

Description

imports pick points selected with meshlab

Usage

read.mpp(file, info = FALSE)

Arguments

Value

if info=FALSE:

a matrix containing picked-points coordinates (only those tagged as active).

if info=TRUE: a list containing

Author(s)

Stefan Schlager

See Also

read.pts

reads pts files

Description

reads Landmark data exported from the software Landmark from http://graphics.idav.ucdavis.edu/research/EvoMorph

Usage

read.pts(file = "x", na = 9999)

Arguments

Value

See Also

read.pts

Examples

data(nose)
write.pts(shortnose.lm, filename="shortnose")
data <- read.pts("shortnose.pts")

read Landmarks from Slicer in Json format

Description

read Landmarks from Slicer in Json format

Usage

read.slicerjson(x, lps2ras = FALSE, na = NULL)

Arguments

Value

returns matrix or list of matrices with imported landmark coordinates

import landmark data from csv files

Description

import landmark data from csv files

Usage

readLandmarks.csv(
  file,
  x,
  y = 2:4,
  rownames = NULL,
  header = TRUE,
  dec = ".",
  sep = ";"
)

Arguments

Value

Author(s)

Stefan Schlager

See Also

read.table

Import landmarks and outlines from TPS files

Description

Imports outlines and landmarks from files generated by tpsdig2

Usage

readallTPS(file, scale = TRUE)

Arguments

Value

Note

currently only landmarks, ID and outlines are read from the TPS-file

Author(s)

Stefan Schlager

References

http://life.bio.sunysb.edu/ee/rohlf/software.html

See Also

read.lmdta, read.pts

correlation between shape space and tangent space

Description

performs a partial Procrustes superimposition of landmark data and calculates the correlation between tangent and shape space.

Usage

regdist(
  dataarray,
  plot = TRUE,
  main = "",
  rho = "angle",
  dist.mat.out = FALSE,
  ...
)

Arguments

Value

Author(s)

Stefan Schlager

See Also

regdist

Examples

if (require(shapes)) {
regdist(gorf.dat)
}

calculate relative Warp analysis

Description

After Procrustes registration the data is scaled by the bending energy or its inverse to emphasize global/local differences when exploring a sample's shape.

Usage

relWarps(
  data,
  scale = TRUE,
  CSinit = TRUE,
  alpha = 1,
  tol = 1e-10,
  orp = TRUE,
  pcAlign = TRUE,
  computeBasis = TRUE,
  noalign = FALSE
)

Arguments

Value

Author(s)

Stefan Schlager

References

Bookstein FL 1989. Principal Warps: Thin-plate splines and the decomposition of deformations. IEEE Transactions on pattern analysis and machine intelligence 11.

Bookstein FL, 1991. Morphometric tools for landmark data. Geometry and biology. Cambridge Univ. Press, Cambridge.

Rohlf FJ, Bookstein FL 2003. Computing the Uniform Component of Shape Variation. Systematic Biology 52:66-69.

Examples

data(boneData)
pop <- name2factor(boneLM,which=3)
rW <- relWarps(boneLM, alpha = -1)
## Not run: 
if (require(car)) {
# plot first 5 relative warps scores grouped by population
spm(rW$bescores[,1:5],group=pop)
# plot uniform component scores grouped by population
spm(rW$uniscores[,1:5],group=pop)
}
##plot non-affine variance associated with each relative warp
barplot(rW$Var[,2], xlab="relative Warps")
## visualize first relative warp +-3 sd of the scores
rw1 <- restoreShapes(as.matrix(c(-3,3)*sd(rW$bescores[,1])),rW$bePCs[,1,drop=FALSE],rW$mshape)
deformGrid3d(rw1[,,1],rw1[,,2],ngrid=5)

## 2D example:
if (require(shapes)) {
data <- bindArr(gorf.dat, gorm.dat, along=3)
sex <- factor(c(rep("fem", dim(gorf.dat)[3]), rep("male",dim(gorm.dat)[3])))
rW <- relWarps(data, alpha = -1)
if (require(car)) {
# plot first 3 relative warps scores grouped by population
spm(rW$bescores[,1:3],group=sex)
# plot uniform component scores grouped by population
spm(rW$uniscores[,1:2],group=sex)
}
##plot non-affine variance associated with each relative warp
barplot(rW$Var[,2], xlab="relative Warps")
## visualize first relative warp +-3 sd of the scores
rw1 <- restoreShapes(as.matrix(c(-3,3)*sd(rW$bescores[,1])),rW$bePCs[,1,drop=FALSE],rW$mshape)
deformGrid2d(rw1[,,1],rw1[,,2],ngrid=10)
}
## End(Not run)

relax one specific landmark configuration against a reference

Description

relax one specific landmark configuration against a reference (e.g. a sample mean)

Usage

relaxLM(lm, ...)

## S3 method for class 'matrix'
relaxLM(
  lm,
  reference,
  SMvector,
  outlines = NULL,
  surp = NULL,
  sur.name = NULL,
  mesh = NULL,
  tol = 1e-05,
  deselect = FALSE,
  inc.check = TRUE,
  iterations = 0,
  fixRepro = TRUE,
  missing = NULL,
  bending = TRUE,
  stepsize = ifelse(bending, 1, 0.5),
  use.lm = NULL,
  silent = FALSE,
  ...
)

## S3 method for class 'mesh3d'
relaxLM(
  lm,
  reference,
  tol = 1e-05,
  deselect = FALSE,
  inc.check = TRUE,
  iterations = 0,
  fixRepro = TRUE,
  missing = NULL,
  bending = FALSE,
  stepsize = ifelse(bending, 1, 0.5),
  use.lm = NULL,
  silent = FALSE,
  ...
)

Arguments

Details

if lm is a surface mesh, all vertices will be treated as semilandmarks and a allowed to freely slide along the surface.

Value

returns kx3 matrix of slidden landmarks

Author(s)

Stefan Schlager

References

Gunz, P., P. Mitteroecker, and F. L. Bookstein. 2005. Semilandmarks in Three Dimensions, in Modern Morphometrics in Physical Anthropology. Edited by D. E. Slice, pp. 73-98. New York: Kluwer Academic/Plenum Publishers.

See Also

slider3d

Examples

require(rgl)
data(nose)
### relax shornose against longnose

# define fix landmarks
fix <- c(1:5,20:21)
# define surface patch by specifying row indices of matrices
# all except those defined as fix
surp <- c(1:dim(shortnose.lm)[1])[-fix]
 
relax <- relaxLM(shortnose.lm,
         longnose.lm, mesh=shortnose.mesh, iterations=1,
         SMvector=fix, deselect=TRUE, surp=surp)

## example minimizing Procrustes distance when displacement is not
## dampened by stepsize
relaxProcD <- relaxLM(shortnose.lm,
         longnose.lm, mesh=shortnose.mesh, iterations=1,
         SMvector=fix, deselect=TRUE, surp=c(1:623)[-fix],bending=FALSE,stepsize=1)

## Not run: 
# visualize differences red=before and green=after sliding
deformGrid3d(shortnose.lm, relax, ngrid=0)

 
# visualize differences minimizing Procrusted distances red=before and green=after sliding

deformGrid3d(shortnose.lm, relaxProcD, ngrid=0)
## no smooth displacement, now let's check the distances:
rot2ref <- rotonto(relaxProcD,longnose.lm)
angle.calc(rot2ref$X,rot2ref$Y)
# 0.2492027 Procrustes distance between reference and slided shape
# (minimizing Procrustes distance)
rot2refBend <- rotonto(relax,longnose.lm)
angle.calc(rot2refBend$X,rot2refBend$Y)
# 0.2861322 Procrustes distance between reference and slided shape
# (minimizing bending energy)

rot2refOrig <- rotonto(shortnose.lm,longnose.lm)
angle.calc(rot2refOrig$X,rot2refOrig$Y)
# 0.3014957 Procrustes distance between reference and original shape
##result: while minimizing Procrustes distance, displacement is not
##guaranteed to be smooth

# add surface
wire3d(shortnose.mesh, col="white")


## finally relax two meshes with corresponding vertices:

mediumnose.mesh <- tps3d(shortnose.mesh,shortnose.lm, (shortnose.lm+longnose.lm)/2,threads=1)
## we use Procrustes distance as criterion as bending energy is pretty slow because
## of too many coordinates (more than 3000 is very unreasonable).
relaxMesh <- relaxLM(shortnose.mesh,mediumnose.mesh,iterations=2,bending=FALSE,stepsize=0.05)

## End(Not run)

plot or save the results of meshDist

Description

plot or save the results of meshDist

Usage

render(x, ...)

## S3 method for class 'meshDist'
render(
  x,
  from = NULL,
  to = NULL,
  steps = NULL,
  ceiling = NULL,
  uprange = NULL,
  tol = NULL,
  tolcol = NULL,
  rampcolors = NULL,
  NAcol = NULL,
  displace = FALSE,
  shade = TRUE,
  sign = NULL,
  add = FALSE,
  scaleramp = NULL,
  titleplot = "Distance in mm",
  ...
)

## S3 method for class 'matrixDist'
render(
  x,
  from = NULL,
  to = NULL,
  steps = NULL,
  ceiling = NULL,
  uprange = NULL,
  tol = NULL,
  tolcol = NULL,
  type = c("s", "p"),
  radius = NULL,
  rampcolors = NULL,
  NAcol = NULL,
  displace = FALSE,
  sign = NULL,
  add = FALSE,
  scaleramp = FALSE,
  titleplot = "Distance in mm",
  ...
)

export(x, ...)

## S3 method for class 'meshDist'
export(
  x,
  file = "default",
  imagedim = "100x800",
  titleplot = "Distance in mm",
  ...
)

Arguments

Details

Visualise or save the results of meshDist to disk.

render.meshDist renders the colored mesh and displays the color ramp and returns an object of class "meshDist". export.meshDist exports the colored mesh as ply file and the color chart as png file.

Author(s)

Stefan Schlager

See Also

meshDist, shade3d

Resample a curve equidistantly

Description

Resample a curve equidistantly (optionally with smoothing)

Usage

resampleCurve(x, n, smooth = FALSE, smoothn = n, open = TRUE)

Arguments

Value

returns a matrix containing the resampled curve

Examples

data(nose)
x <- shortnose.lm[c(304:323),]
xsample <- resampleCurve(x,n=50)

restore original data from PCA

Description

restore original data from PCA by reverting rotation and centering

Usage

restoreFromPCA(scores, rotation, center)

Arguments

Examples

myirispca <- prcomp(iris[,1:4])
myirisRecovered <- restoreFromPCA(myirispca$x,myirispca$rotation,myirispca$center)
all.equal(myirisRecovered,as.matrix(iris[,1:4]))

restore shapes from PC-Scores or similar projections

Description

restore shapes from PC-Scores or similar projections

Usage

restoreShapes(
  scores,
  PC,
  mshape,
  sizeshape = FALSE,
  origsize = FALSE,
  meanlogCS
)

Arguments

Details

Rotates and translates PC-scores (or similar) derived from shape data back into configuration space.

Value

returns matrix or array containing landmarks

Author(s)

Stefan Schlager

See Also

prcomp, procSym

getPCscores

Examples

if (require(shapes)) {
## generate landmarks using
##the first PC-score of the first specimen

proc <- procSym(gorf.dat)
lm <- restoreShapes(proc$PCscores[1,1],proc$PCs[,1],proc$mshape)
plot(lm,asp=1)

##now the first 3 scores
lm2 <- restoreShapes(proc$PCscores[1,1:3],proc$PCs[,1:3],proc$mshape)
points(lm2,col=2)

## Now restore some sizeshape data
procSize <- procSym(gorf.dat,sizeshape=TRUE)
est1 <- restoreShapes(range(procSize$PCscores[,1]),procSize$PCs[,1],procSize$mshape,
                      sizeshape=TRUE,origsize=TRUE,meanlogCS=procSize$meanlogCS)
}

symmetrize a bilateral landmark configuration

Description

symmetrize a bilateral landmark configuration by removing bending and stretching

Usage

retroDeform3d(mat, pairedLM, hmult = 5, alpha = 0.01)

Arguments

Value

References

Ghosh, D.; Amenta, N. & Kazhdan, M. Closed-form Blending of Local Symmetries. Computer Graphics Forum, Wiley-Blackwell, 2010, 29, 1681-1688

symmetrize a triangular mesh

Description

symmetrize a triangular mesh

Usage

retroDeformMesh(
  mesh,
  mat,
  pairedLM,
  hmult = 5,
  alpha = 0.01,
  rot = TRUE,
  lambda = 1e-08,
  threads = 0
)

Arguments

Details

this function performs retroDeform3d and deforms the mesh accordingly using the function tps3d.

Value

Rotate an object (matrix or mesh) around an arbitrary axis in 3D

Description

Rotate an object around an arbitrary axis in 3D

Usage

rotaxis3d(x, pt1, pt2 = c(0, 0, 0), theta)

## S3 method for class 'matrix'
rotaxis3d(x, pt1, pt2 = c(0, 0, 0), theta)

## S3 method for class 'mesh3d'
rotaxis3d(x, pt1, pt2 = c(0, 0, 0), theta)

Arguments

Details

Rotate an object (matrix or triangular mesh) around an 3D-axis defined by two points.

Value

returns rotated object (including updated normals for mesh3d objects)

Author(s)

Stefan Schlager

References

http://en.wikipedia.org/wiki/Rotation_matrix

See Also

rotonto, rotmesh.onto

Examples

require(rgl)
data(nose)
shrot.rot <- rotaxis3d(shortnose.mesh,pt1=c(1,1,1),theta=pi)
## Not run: 
shade3d(shortnose.mesh,col=3,specular=1)
shade3d(shrot.rot,col=2)

###print rotation axis
#' lines3d(rbind(rep(-0.1,3),rep(0.1,3)))

## End(Not run)

calculate a rotation matrix around an arbitrary axis through the origin in 3D

Description

calculate a rotation matrix around an arbitrary axis in 3D

Usage

rotaxisMat(u, theta, homogeneous = FALSE)

Arguments

Value

returns 3x3 rotation matrix

References

http://en.wikipedia.org/wiki/Rotation_matrix

See Also

rotaxis3d

rotate ,scale and translate a mesh based on landmark information.

Description

rotates and reflects a mesh onto by calculating the transformation from two sets of referenced landmarks.

Usage

rotmesh.onto(
  mesh,
  refmat,
  tarmat,
  adnormals = FALSE,
  scale = FALSE,
  reflection = FALSE,
  ...
)

Arguments

Value

Author(s)

Stefan Schlager

See Also

file2mesh,tps3d ,rotonto,mesh2ply

Examples

require(rgl)
data(boneData)
## rotate, translate and scale the mesh belonging to the first specimen
## onto the landmark configuration of the 10th specimen
rotmesh <- rotmesh.onto(skull_0144_ch_fe.mesh,boneLM[,,1],
                        boneLM[,,10], scale=TRUE)
## Not run: 
## render rotated mesh and landmarks
shade3d(rotmesh$mesh, col=2, specular=1)
spheres3d(boneLM[,,1])
## render original mesh
shade3d(skull_0144_ch_fe.mesh, col=3, specular=1)
spheres3d(boneLM[,,10])

## End(Not run)

rotate matrix of landmarks

Description

rotate matrix of landmarks by using a rotation determined by two matrices.

Usage

rotonmat(
  X,
  refmat,
  tarmat,
  scale = TRUE,
  reflection = FALSE,
  weights = NULL,
  centerweight = FALSE,
  getTrafo = FALSE
)

Arguments

Details

A matrix is rotated by rotation parameters determined by two different matrices. This is usefull, if the rotation parameters are to be estimated by a subset of landmark coordinates.

Value

if getTrafo=FALSE the transformed X will be returned, else alist containing:

Author(s)

Stefan Schlager

See Also

rotonto,rotmesh.onto

Examples

data(nose)
shortnose.rot <-
rotonmat(shortnose.lm,shortnose.lm[1:9,],longnose.lm[1:9,])

##view result
## Not run: 
deformGrid3d(shortnose.rot,shortnose.lm,ngrid=0)

## End(Not run)

rotates, translates and scales one matrix onto an other using Procrustes fitting

Description

rotates, translates and scales one matrix onto an other using Procrustes fitting

Usage

rotonto(
  x,
  y,
  scale = FALSE,
  signref = TRUE,
  reflection = TRUE,
  weights = NULL,
  centerweight = FALSE,
  ...
)

rotreverse(mat, rot)

## S3 method for class 'matrix'
rotreverse(mat, rot)

## S3 method for class 'mesh3d'
rotreverse(mat, rot)

Arguments

Details

rotate a matrix onto an other without loosing information about the location of the targetmatrix and reverse this transformations using rotreverse

Value

Note

all lines containing NA, or NaN are ignored in computing the transformation.

Author(s)

Stefan Schlager

References

Lissitz, R. W., Schoenemann, P. H., & Lingoes, J. C. (1976). A solution to the weighted Procrustes problem in which the transformation is in agreement with the loss function. Psychometrika, 41,547-550.

See Also

rotmesh.onto

Examples

if (require(shapes)) {
lims <- c(min(gorf.dat[,,1:2]),max(gorf.dat[,,1:2]))
rot <- rotonto(gorf.dat[,,1],gorf.dat[,,2]) ### rotate the second onto the first config
plot(rot$yrot,pch=19,xlim=lims,ylim=lims) ## view result
points(gorf.dat [,,2],pch=19,col=2) ## view original config
rev1 <- rotreverse(rot$yrot,rot)
points(rev1,cex=2) ### show inversion by larger circles around original configuration
}

scale a mesh of class "mesh3d"

Description

scales (the vertices of a mesh by a scalar

Usage

scalemesh(mesh, size, center = c("bbox", "mean", "none"))

Arguments

Details

The mesh's center is determined either as mean of the bounding box (center="bbox") or mean of vertex coordinates (center="mean") and then scaled according to the scaling factor. If center="none", vertex coordinates will simply be multiplied by "size".

Value

returns a scaled mesh

Author(s)

Stefan Schlager

See Also

rotmesh.onto

Examples

data(nose)
#inflate mesh by factor 4
largenose <- scalemesh(shortnose.mesh,4)

slides Semilandmarks along curves 2D by minimising bending energy of a thin-plate spline deformation.

Description

slides Semilandmarks along curves 2D. The positions are sought by minimising bending energy (of a thin-plate spline deformation) or Procrustes distance

Usage

slider2d(
  dataframe,
  SMvector,
  outlines,
  tol = 1e-05,
  deselect = FALSE,
  recursive = TRUE,
  iterations = 0,
  initproc = FALSE,
  pairedLM = NULL,
  bending = TRUE,
  stepsize = 1,
  silent = FALSE
)

Arguments

Value

returns an array containing slided coorndinates in the original space - not yet processed by a Procrustes analysis.

Warning

Depending on the amount of landmarks this can use an extensive amount of your PC's resources, especially when running in parallel. As the computation time and RAM usage of matrix algebra involved is quadratic to the amount of landmarks used, doubling the amount of semi-landmarks will quadruple computation time and system resource usage. You can easily stall you computer with this function with inappropriate data.

Author(s)

Stefan Schlager

See Also

relaxLM, slider3d

slides Semilandmarks along curves and surfaces in 3D by minimising bending energy of a thin-plate spline deformation.

Description

slides Semilandmarks along curves and surfaces in 3D. The positions on the surface are sought which minimise bending energy (of a thin-plate spline deformation)

Usage

slider3d(
  dat.array,
  SMvector,
  outlines = NULL,
  surp = NULL,
  sur.path = NULL,
  sur.name = NULL,
  meshlist = NULL,
  ignore = NULL,
  sur.type = "ply",
  tol = 1e-05,
  deselect = FALSE,
  inc.check = TRUE,
  recursive = TRUE,
  iterations = 0,
  initproc = TRUE,
  fullGPA = FALSE,
  pairedLM = 0,
  bending = TRUE,
  stepsize = ifelse(bending, 1, 0.5),
  mc.cores = parallel::detectCores(),
  fixRepro = TRUE,
  missingList = NULL,
  use.lm = NULL,
  smoothnormals = FALSE,
  silent = FALSE
)

Arguments