-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathslide_semiLMs.Rmd
86 lines (58 loc) · 3 KB
/
slide_semiLMs.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
---
title: "Sliding semiLM via GeoMorph"
author: "SlicerMorph"
date: "4/02/2021"
output:
html_document: default
pdf_document: default
---
This script is based on geomorph 3.3.2.
In this example we are using the data from the *Gorilla patch semiLMs* dataset distributed by SlicerMorph (available via `File->Download Sample Data`). This particular dataset contains four Gorilla specimens, each of which has a fixed and semiLM file pair. First use the `MergeMarkups` module in batch mode to combine fixed and semiLMs pairs into a single file, and then run the SlicerMorph GPA module on them. More information about `MergeMarkups` module can be found [here](https://github.com/SlicerMorph/SlicerMorph/tree/master/Docs/MergeMarkups).
After this the first task is to load the geomorph library and the parser convenience function to pull all the files and analytical settings.
These functions eventually will be turned into an R package.
```{r setup}
library(geomorph)
source("https://raw.githubusercontent.com/muratmaga/SlicerMorph_Rexamples/main/log_parser.R")
```
Next we need to point to the location of the analysis.log file that was saved by SlicerMorph's GPA module either via coding the path OR interactively.
If you are knitting the document, make sure to hard code the path.
```{r load log file}
# interactively choose log file
# SM.log.file <- file.choose()
# coding the path to log file example
SM.log.file <- "c:/temp/RemoteIO/gorilla_patches/merged/2021-04-05_20_51_17/analysis.log"
SM.log <- parser(SM.log.file)
```
Check if there are any LMs that are tagged as semiLMs
```{r check semiLM indices}
if (SM.log$semi) print (SM.log$semiLMs)
```
Here, as an example we fit three different GPA superimposition: Two with sliding (via bending energy and procrustes distance) and one without sliding.
Note that depending on the sample size and LM numbers, bending energy option may take a few minutes to complete.
```{r GPA with and without sliding}
gpa.bending = gpagen(A=SM.log$LM,
surfaces = as.numeric(SM.log$semiLMs),
ProcD = FALSE,
print.progress = FALSE)
# Sliding semiLMs via optimizing procrustes distance
gpa.procD = gpagen(A=SM.log$LM,
surfaces = as.numeric(SM.log$semiLMs),
ProcD = TRUE,
print.progress = FALSE)
# no sliding
gpa = gpagen(A=SM.log$LM,
print.progress = FALSE)
```
Lets plot the mean shapes of each of these to visualize the differences
```{r plot consensus shapes}
par(mfrow=c(2,2))
plot(gpa$consensus[,c(1,2)], pch=20, col='black')
points(gpa.procD$consensus[,c(1,2)], pch=20, col='red')
points(gpa.bending$consensus[,c(1,2)], pch=20, col='green')
plot(gpa$consensus[,c(1,3)], pch=20, col='black')
points(gpa.procD$consensus[,c(1,3)], pch=20, col='red')
points(gpa.bending$consensus[,c(1,3)], pch=20, col='green')
plot(gpa$consensus[,c(2,3)], pch=20, col='black')
points(gpa.procD$consensus[,c(2,3)], pch=20, col='red')
points(gpa.bending$consensus[,c(2,3)], pch=20, col='green')
```