library(geomorph) strike_table= read.table("~/strike info.txt", header=T) # info associated with 3260 motions shapes (326 feeding srikes*10 frames per strike) ######################################################################## ### MOTION SHAPES: scaling and alignment of shape data ######################################################################## coord_data= readland.tps(file="~/cichlid feeding LM.TPS") # landmark coordinate data semiland= as.matrix(read.table("~/sliding_semi.txt", header=T)) # file indicating semilandmarks along ventral profile of head unaligned_2d= two.d.array(coord_data) # creates a 2D table of landmarks n_specimen=326 # number of feeding strikes/motions scaled= unaligned_2d/as.numeric(strike_table[,6]) # scales coordinate data p=18 # number of landmarks k=2 # dimensionality of shape data scaled_3d= arrayspecs(scaled, p=p, k=k) # 3D array of shape data, needed for alignment with gpagen function, below y= gpagen(scaled_3d, curves=semiland, ProcD=T) # generalized Procrustes analysis (i.e. "shape alignment") y2= two.d.array(y$coords) # 2D table of scaled and aligned landmarks ############################################################################## ## TRAJECTORY LENGTHS & KINESIS CALCULATION ############################################################################## kinetic_procrustes= as.data.frame(matrix(ncol=n_specimen, nrow=11)) # table of Procrustes distances b/t successive motion points (columns represent individual motions) for(f in 1:n_specimen){ spec= y2[10*(f-1) +c(1:10),] # selecting the 10 shapes for motion "f" of 326 motions for(i in 1:9){ proc= sqrt(sum((spec[i,]-spec[i+1,])^2)) # Procrustes distance between shape "i" along motion trajectory and next shape (i+1) kinetic_procrustes[i,f]= proc # successive Procrustes distances b/t motion shapes kinetic_procrustes[10,f]= sum(kinetic_procrustes[1:9,f]) # total path length (i.e., "kinesis") kinetic_procrustes[11,f]= sqrt(sum((spec[1,]-spec[10,])^2)) # Procrustes distance between start and end shapes (i.e., length of linear trajectory) }} row.names(kinetic_procrustes)= c(1:9, "total", "linear") ############################################################################## ## ESTIMATING LINEAR TRAJECTORY SHAPES ############################################################################## kinetic_proportional= as.data.frame(matrix(ncol=n_specimen, nrow=9)) # Table of successive Procrustes distances, expressed as proportions of total path length (i.e., "spacing regime" from materials and methods) for(t in 1:n_specimen){ kinetic_proportional[,t]= kinetic_procrustes[1:9,t]/kinetic_procrustes[10,t] } linear_traj= as.data.frame(matrix(ncol=36, nrow=3260)) # new shape data estimated along a linear trajectory between start and end motion shapes for(f in 1:n_specimen){ spec= y2[10*(f-1) + c(1:10),] # selecting the 10 shapes for motion "f" out of 326 motions start_motion=spec[1,] # first motion shape end_motion= spec[10,] # last motion shape diff_start_end= end_motion-start_motion for(i in 1:8){ linear_traj[rownames(spec)[1],]= start_motion linear_traj[rownames(spec)[i+1],]= linear_traj[rownames(spec)[i],] + ((diff_start_end)*kinetic_proportional[i,f]) # calculates linear shapes between start & end shapes linear_traj[rownames(spec)[10],]= end_motion }} ############################################################################## ## TRAJECTORY NONLINEARITY & KINEMATIC EFFICIENCY ############################################################################## linear_v_actual= as.data.frame(matrix(ncol=5, nrow=n_specimen*10)) # table of procrustes distances b/t matching shapes on linear vs. actual motion trajectories (see ms. fig 1c) for(g in 1:3260){ p_dist= sqrt(sum((y2[g,]-linear_traj[g,])^2)) linear_v_actual[g,5]= p_dist # distances from linear at each motion stage linear_v_actual[,1]= strike_table[1] # strike number linear_v_actual[,2]= strike_table[3] # species linear_v_actual[,3]= as.factor(rep(1:10,n_specimen)) # sequence (1:10) linear_v_actual[,4]= strike_table[7] # individual specimen } linear_distances= as.data.frame(t(kinetic_procrustes[11,])) # linear path length (b/t start & end shapes) for each motion rownames(linear_distances)= c(1:326) linear_distances_all= as.data.frame(matrix(ncol=1, nrow=3260)) # expands linear path length data to match table of linear deviations (i.e., "linear_v_actual") for(f in 1:n_specimen){ list_linear= rep(linear_distances[f,], 10) linear_distances_all[10*(f-1) +c(1:10),1]= list_linear } noncorrected_deviation= linear_v_actual[,5] deviations_table_raw= as.data.frame(matrix(nrow=n_specimen, ncol=10)) # generates a table of raw deviations from linear (NOT CORRECTED FOR LINEAR DISTANCE) for(f in 1:n_specimen){ devs_b= noncorrected_deviation[10*(f-1) + c(1:10)] deviations_table_raw[f,]= devs_b names(deviations_table_raw)= c("dev_start","dev_2","dev_3","dev_4","dev_5","dev_6","dev_7","dev_8","dev_9","dev_end") } corrected_deviation= linear_v_actual[,5]/linear_distances_all deviations_table= as.data.frame(matrix(nrow=n_specimen, ncol=10)) # generates table of deviations, adjusted by length of linear trajectory (USED AS "KINEMATIC EFFICIENCY") for(f in 1:n_specimen){ devs_b= corrected_deviation[10*(f-1) +c(1:10), 1] deviations_table[f,]= devs_b names(deviations_table)= c("dev_start","dev_2","dev_3","dev_4","dev_5","dev_6","dev_7","dev_8","dev_9","dev_end") }