Procrustes

Library

library(shapes)

## 
## Attaching package: 'shapes'

## The following object is masked from 'package:stats':
## 
##     sigma

library(ggvis)

Set resporitory and read data

setwd("C:/Users/Admin/Documents/Intern/practice")
data <- read.csv("wings.csv")
data1 <- data.matrix(data[,c("ID","X","Y","isBoundary","isLandmark")])

Procrustes analysis

Load shapes

data_split <- split(data, data$ID)
wing1 <-  data.frame(data_split[[1]])
landmark1 <- split(wing1, wing1$isBoundary)
landmark1 <- data.matrix(landmark1[[2]])
rownames(landmark1)<- NULL

wing2 <-  data.frame(data_split[[2]])
landmark2 <- split(wing2, wing2$isBoundary)
landmark2 <- data.matrix(landmark2[[2]])
rownames(landmark2)<- NULL

standard = landmark1[,1:2]
data = landmark2[,1:2]

Procrustes

out <- procOPA(standard,data)

Plot to test

• Standard and Ahat

plot(standard,xlim=c(-1500,3000),ylim=c(-1500,3000),xlab=" ",ylab=" ")
lines(standard[c(1:19,1),])
points(out$Ahat)
lines(out$Ahat[c(1:19,1),],lty=2)

-Bhat and data

plot(data,xlim=c(-1500,3000),ylim=c(-1500,3000),xlab=" ",ylab=" ")
lines(data[c(1:19,1),])
points(out$Bhat,pch =2,col="red")
lines(out$Bhat[c(1:19,1),],lty=2, col="red")

• Standard and Bhat

plot(standard,xlim=c(-1500,3000),ylim=c(-1500,3000),xlab=" ",ylab=" ")
lines(standard[c(1:19,1),])
points(out$Bhat,pch =2, col = "red")
lines(out$Bhat[c(1:19,1),],lty=2,col="red")

Processing point Test for 1 point

temp2 <- split(wing2, wing2$isBoundary)

temp1 <- data.matrix(temp2[[1]])
rownames(temp1)<- NULL

boundary<- data.matrix(temp2[[2]])
boundary<- boundary[,1:2]
rownames(boundary)<- NULL

data1 = temp1[,1:2]
temp_point = matrix(data1[15,], nrow =1, ncol = 2)

1.Scaling

1.1 Scaling point first

point_scaling = out$s * temp_point
print(point_scaling)

##          [,1]     [,2]
## [1,] 1506.882 1032.281

1.2 Scaling boundary

boundary_scaling = out$s * boundary 
print(boundary_scaling)

##               X         Y
##  [1,]  237.8193 1151.7103
##  [2,]  393.5962 1203.6359
##  [3,]  733.1898 1263.8696
##  [4,] 1055.1286 1310.6027
##  [5,] 1523.4977 1414.4539
##  [6,] 1970.0581 1449.7634
##  [7,] 2337.6915 1395.7607
##  [8,] 2497.6224 1282.5629
##  [9,] 2572.3953 1198.4433
## [10,] 2617.0513 1111.2083
## [11,] 2516.3156 1023.9733
## [12,] 2363.6543  930.5071
## [13,] 2142.4511  837.0410
## [14,] 1925.4020  766.4222
## [15,] 1719.7766  735.2668
## [16,] 1340.7195  795.5005
## [17,]  794.4620  985.5483
## [18,]  482.9083 1070.7063
## [19,]  343.7476 1089.3995

2. Rotating 2.1 Rotating point

temp_rotate_x = point_scaling[1]* out$R[1,1] - point_scaling[2]* out$R[1,2]
temp_rotate_y = point_scaling[1]* out$R[1,2] + point_scaling[2]* out$R[1,1]
point_rotating = cbind(temp_rotate_x, temp_rotate_y)
#Points(point_rotating)

2.2 Rotating boundary

output_rotating = matrix(ncol=2)
output_rotating = output_rotating[-1,]
for (i in 1:nrow(boundary_scaling)){
    temp_x = boundary_scaling[i,1]* out$R[1,1] - boundary_scaling[i,2]* out$R[1,2]
    temp_y = boundary_scaling[i,1]* out$R[1,2] + boundary_scaling[i,2]* out$R[1,1]
    point_rotate_data = cbind(temp_x, temp_y)
    output_rotating <- rbind(output_rotating, point_rotate_data)    
}
rownames(output_rotating)<- NULL
#points(output_rotating)

3. Translating point

Value of vecctor of translation

vector = colMeans(out$Bhat - output_rotating )
print(vector)

##    temp_x    temp_y 
## -1566.599 -1091.324

Translation of point

point_translating = point_rotating + vector 
print(point_translating)

##      temp_rotate_x temp_rotate_y
## [1,]     -49.99977     -73.37304

Plot test

plot(data,xlim=c(-1500,3000),ylim=c(-1500,3000),xlab=" ",ylab=" ")
lines(data[c(1:19,1),])
points(out$Bhat,pch =2, col="red")
lines(out$Bhat[c(1:19,1),],lty=2, col="red")
points(temp_point,pch =2)
points(point_translating,pch =2,col="red")