New fixed.angle() Function

Erik Otárola-Castillo

9 years ago

[This article was first published on geomorph, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

< !DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">< http-equiv="Content-Type" content="text/html; charset=utf-8" />

Hello morphometricians,

Below you can find a new fixed angle function addressing the problem discovered by Fabio Machado in the morphmet mail archive. We will include this function in our next schedule update to geomorph.

Cheers,

Erik

CODE:

fixed.angle<-function(A,art.pt=NULL,angle.pts=NULL,rot.pts=NULL,angle=0){
if (length(dim(A))!=3){
stop("Data matrix 1 not a 3D array (see ‘arrayspecs’).") }
if(length(grep("-999",A))!=0){
stop("Data matrix 1 contains missing values. Estimate these first(see ‘estimate.missing’).") }
n<-dim(A)[3]; k<-dim(A)[2]; p<-dim(A)[1]
if (k!=2){
stop("Method presently implemented for two-dimensional data only.")}
if (angle>pi*2){
stop("Additional angle must be specified in radians.")}
if (angle< -pi*2){
stop("Additional angle must be specified in radians.")}
angl<-array(NA,dim=n)
for (i in 1:n){
A[,,i]<-t(t(A[,,i])-A[art.pt,,i])
angl[i]<- acos((A[angle.pts[1],,i]/sqrt(sum(A[angle.pts[1],,i]^2)))%*%(A[angle.pts[2],,i]/sqrt(sum(A[angle.pts[2],,i]^2))))
}
dev.angle<- (angl-mean(angl))+angle
if(A[angle.pts[1],1,1]<0){dev.angle<- -1*dev.angle}
for (i in 1:n){
r = matrix(c(cos(dev.angle[i]),-sin(dev.angle[i]),sin(dev.angle[i]),cos(dev.angle[i])),2)
A[rot.pts,,i] = A[rot.pts,,i]%*%r
}
return(A)
}

To leave a comment for the author, please follow the link and comment on their blog: geomorph.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.