Triangulation for Neuroimage Splines
En esencia trabajo con un paquete - Triangulation, de Funstats - que saca, a partir de una figura 2D, su Triangulacion de Delaunay, esto es, una red de triángulos conexa y convexa que cumple la condición de Delaunay. Esta condición dice que la circunferencia circunscrita de cada triángulo de la red no debe contener ningún vértice de otro triángulo. Las triangulaciones de Delaunay tienen importante relevancia en el campo de la geometría computacional).
De hecho ya llevo trabajando un tiempo con este paquete solo que de forma pasiva. En general, para obtener los Statistical Confidence Corridors (SCC’s) hace falta, antes de nada, dos objetos matriciales: (1) una matriz nx2 con las coordenadas de los vértices, y (2) una matriz nx3 donde, una vez definidos los vértices, se indica cuales van a conformar cada triángulo.
Lo que veo es que las funciones de este paquete inicial - que por cierto, sus outputs tienen un valor artístico innegable a mi gusto - son la base sobre la que realizar los splines bivariados penalizados. Voy a poner a continuación un par de ejemplos a ver qué tal se entiende:
Preambulo:
# Install packages from GitHub:
install.packages("devtools");library(devtools)
install.packages("remotes");library(remotes)
install.packages("readr")
install.packages("imager")
install.packages("itsadug")
Sys.setenv("R_REMOTES_NO_ERRORS_FROM_WARNINGS"=TRUE) # Forces installation to finish even with warning messages (EXTREMELY NECESSARY)
remotes::install_github("funstatpackages/BPST")
remotes::install_github("funstatpackages/Triangulation")
remotes::install_github("funstatpackages/ImageSCC")
library(BPST);library(Triangulation);library(ImageSCC);library(readr);library(imager);library(itsadug);library(fields)
Vale, con esto deberíamos tener todos los paquetes instalados directamente del GitHub de este grupo de chin@s de Ames (Iowa, que no Ames en Galicia). Se hace complicado descargar estos paquetes sin esa línea de código de Sys.setenv de modo que no se debe quitar.
Ejemplos de Triangulaciones:
Pues empezamos así. Tenemos una serie de datos llamados BMP que constan de varios objetos en la lista: los bound que no dejan de ser las fronteras o límites del polígono, los holes (1 y 2) que son fronteras internas que delimitan ‘agujeros’ dentro de la figura.
Tanto BMP$bound como BMP$H1 y BMP$H2 son data.frames con dos columnas (dos coordenadas) que ni siquiera están acotadas ni nada. Muy servicial el tema vaya.
library(Triangulation)
data("BMP") # Propios del paquete
head(BMP$bound,10) # Tenemos 3 datasets en BMP: los límites exteriores (bound) y dos huecos interiores (H1 y H2)
## V1 V2
## 1 -3.618052 14.35317
## 2 -3.617937 14.35378
## 3 -3.617673 14.35444
## 4 -3.617440 14.35449
## 5 -3.616956 14.35430
## 6 -3.615316 14.35396
## 7 -3.615193 14.35374
## 8 -3.614754 14.35452
## 9 -3.614028 14.35576
## 10 -3.613829 14.35578
head(BMP$H1,5)
## V1 V2
## 1 -3.615044 14.34641
## 2 -3.615346 14.34726
## 3 -3.614881 14.34736
## 4 -3.614328 14.34711
## 5 -3.613457 14.34677
head(BMP$H2,5)
## V1 V2
## 1 -3.605706 14.34656
## 2 -3.605935 14.34605
## 3 -3.606552 14.34550
## 4 -3.607993 14.34582
## 5 -3.608298 14.34564
str(BMP)
## List of 3
## $ bound:'data.frame': 95 obs. of 2 variables:
## ..$ V1: num [1:95] -3.62 -3.62 -3.62 -3.62 -3.62 ...
## ..$ V2: num [1:95] 14.4 14.4 14.4 14.4 14.4 ...
## $ H1 :'data.frame': 21 obs. of 2 variables:
## ..$ V1: num [1:21] -3.62 -3.62 -3.61 -3.61 -3.61 ...
## ..$ V2: num [1:21] 14.3 14.3 14.3 14.3 14.3 ...
## $ H2 :'data.frame': 11 obs. of 2 variables:
## ..$ V1: num [1:11] -3.61 -3.61 -3.61 -3.61 -3.61 ...
## ..$ V2: num [1:11] 14.3 14.3 14.3 14.3 14.3 ...
summary(BMP)
## Length Class Mode
## bound 2 data.frame list
## H1 2 data.frame list
## H2 2 data.frame list
plot(BMP$bound)
En lo único que parece que hay que tener un poco de ATENCIÓN es que los agujeros dentro del polígono tienen que ser previamente definidos como lista en caso de que sean múltiples, como es este caso.
A parte de eso, utilizando la función TriMesh se puede sacar la triangulación del polígono indicado (sería BMP$bound) con una n determinada. n es un argumento para la función que nos va a indicar la “finura de la triangulación”, literalmente esa sería la traducción. Se va a entender fácilmente de manera gráfica.
Nota: TriMesh muestra ya el resultado de la triangulación, pero con TriPlot se pueden hacer ciertas modificaciones como el color y el grosor de la línea.
Triangulación a n=8 con TriMesh y TriPlot:
VT8=TriMesh(BMP$bound,8,list(as.matrix(BMP$H1),as.matrix(BMP$H2)))
TriPlot(VT8$V, VT8$Tr, col = 3, lwd = 4)
Triangulación a n=25:
VT25=TriMesh(BMP$bound,25,list(as.matrix(BMP$H1),as.matrix(BMP$H2)))
TriPlot(VT25$V, VT25$Tr, col = 2, lwd = 3)
Lo mismo con otras formas del estilo:
data("shape")
plot(shape)
shape<-TriMesh(shape, n=8, H = NULL)
TriPlot(shape$V, shape$Tr, col = 5, lwd = 6)
La cuestión entonces, para poder tomar el código de los SCC’s y aplicarlo a datos reales de la base de datos ADNI, es obtener una lista de data.frames tal que estas con las que trabaja el paquete, pero para mi imagen PET. En un principio supongo que lo más práctico va a ser hacerlo para el corte en el que estoy trabajando (Data$z==30), aunque lo propio sería programar una función en loop que me sacase la triangulación para todos los niveles. De este modo tendría un data.frame con Vertices y otro con Triangulation vertices. Same as:
head(VT8$V) # Coordenadas de cada vertice
## X Y
## [1,] -3.619323 14.34922
## [2,] -3.619323 14.35120
## [3,] -3.619323 14.35318
## [4,] -3.617343 14.34922
## [5,] -3.617343 14.35120
## [6,] -3.617343 14.35318
head(VT8$Tr) # Combinacion de vertices que forma la triangulacion
## node1 node2 node3
## [1,] 1 125 124
## [2,] 1 124 122
## [3,] 1 122 121
## [4,] 1 121 4
## [5,] 1 4 2
## [6,] 1 2 125
Pero, para eso, antes necesito un data.frame de coordenadas tal que el de:
head(BMP$bound)
## V1 V2
## 1 -3.618052 14.35317
## 2 -3.617937 14.35378
## 3 -3.617673 14.35444
## 4 -3.617440 14.35449
## 5 -3.616956 14.35430
## 6 -3.615316 14.35396
head(BMP$H1)
## V1 V2
## 1 -3.615044 14.34641
## 2 -3.615346 14.34726
## 3 -3.614881 14.34736
## 4 -3.614328 14.34711
## 5 -3.613457 14.34677
## 6 -3.611361 14.34572
head(BMP$H2)
## V1 V2
## 1 -3.605706 14.34656
## 2 -3.605935 14.34605
## 3 -3.606552 14.34550
## 4 -3.607993 14.34582
## 5 -3.608298 14.34564
## 6 -3.608526 14.34607
¿Cómo puedo obtener yo esto a partir de una imagen PET? Pues eso es lo que habrá que descubrir. Con las funciones del paquete Triangulation puedo hacer la triangulación de Delaunay pero…NECESITO SABER CON LA MÁXIMA EXACTITUD CUALES SON LAS COORDENADAS DE MIS DATOS FRONTERA.
Obtener coordenadas para PETriangulation:
Antes de nada vamos a cargar una imagen que sirva de ejemplo. Cargo una que ya he utilizado antes que es, sin mas, una PET image con la máscara ya puesta. Por alguna razón no me deja cargar la máscara en si, aunque también debería ya que es un NIFTI. Sea como sea, la cuestión es que cualquiera de las imágenes vale ya que todas tienen la misma máscara (tridimensional). Con tal de sacarlo una vez, podría crearse una función y loopearla para sacarla de los niveles que se quiera. El corte z=30 se ve que no es el ideal ya que quisiera ver cómo funciona el tema para los huecos interiores. Este corte tiene 2 huecos interiores muy pequeños que deberían bastar, pero la verdad es que también habrá que comprobar cómo funciona en otras z’s con ‘holes’ más consistentes.
Nota de Juan del futuro: En realidad este corte es bastante bueno porque tiene dos agujeros interiores minúsculos, lo que pone a prueba el funcionamiento del código para ver si en efecto detecta o no detecta esos ‘holes’. A priori parece que sí que funciona bien, por lo que debería también funcionar correctamente para agujeros más grandes.
Como era de esperar, se ve en el histograma de los datos asociados a la imagen que una gran cantidad de celdas consisten en un valor de 0, mientras que el resto tienen otros valores acotados que, para el caso que nos atañe, serán 1’s para así sacar una imagen completamente binaria de la que obtener la triangulación. Esto sería lo mismo que trabajar con la máscara original.
Nota de Juan del futuro: Finalmente no merece la pena convertir la imagen a 0’s y 1’s ya que hay una gran diferencia entre los valores nulos (ceros) y los no nulos. La idea de convertir la imagen en binaria era para facilitar el proceso y que no hubiera píxeles “intermedios” que produjesen artefactos.
library(mgcv);library(gamair);library(oro.nifti);library(memisc)
#Lectura de la imagen NIFTI y muestra de datos generales de la Neuro-Imagen que vamos a utilizar
setwd("C:/Hugo/Sites/messy-dataset/content/post/triangulation_for_neuroimage")
sample <- readNIfTI("maskedPETdataEXAMPLE", verbose = FALSE, warn = -1, reorient = TRUE,
call = NULL, read_data = TRUE)
#Visualizacion de la imagen completa
image(sample,plane="axial",plot.type = "single",z=30,main="Neuro-imagen en plano axial Z=30, se aprecian los agujeros en los datos que han de ser captados por la triangulación") # single axial en Z=30
#Datos asociados a la imagen
img_data=img_data(sample) # Nos quedamos con los datos crudos como objeto R
hist(img_data) # Histograma mostrando que mayormente hay valores de 0 y luego valores de PET !=0
No obstante, aunque sabemos que hay celdas con un valor nulo y otras que no, no nos valdría de nada a estas alturas del proceso convertirlas ya a una imagen binaria, ya que no sabríamos las coordenadas que corresponden a dicho valor de 0 o de 1. Quizás haya una forma más rápida para acceder desde aquí, pero teniendo en cuenta la complejidad del data.frame a estas alturas, lo mejor será tirar de fórmulas comprobadas como f.clean, una función ya creada previamente por mí en 2019 que debería ahorrar mucho timepo. Voy a comprobar que es aplicable al 100% a los datos actuales pero debería serlo.
- FUNCIÓN PROPIA: f.clean
##########################################################
### f.clean (PPT,z,x,y,pet) ###
##########################################################
## Set as working directory -> directory with .hdr/.img NIfTI files and Demographics
f.clean <- function(name) { #### f.clean is meant for CLEANING ONE SINGLE PPT DATA
# Read the NIFTI image, transform it to dataframe, preserve slice Z and organize the table
## Load Data
# setwd("D:/Usuario/Desktop/T-tests")
file <- readNIfTI(fname = name, verbose = FALSE, warn = -1, reorient = TRUE, call = NULL, read_data = TRUE)
namex <- as.character(name)
n = img_data(file)
n = to.data.frame(n)
## Prepare data.frame base where surther data from the loop will be integrated
dataframe <- data.frame(z=integer(),x=integer(),y=integer(),pet=integer())
# Loop for 79 slices of Z in the NiFtI image -> move to dataframe
for (i in seq(1:79)) {
n_lim = n[n$Var2==i,] # Select just one Z slice
n_lim$Var1=NULL
n_lim$Var2=NULL
z <-rep(i,length.out=7505)
x <-rep(1:79, each=95, length.out = 7505)
y <-rep(1:95,length.out = 7505)
attach(n_lim)
pet<-c(`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`)
detach(n_lim)
temp0 = data.frame(z,x,y,pet) # temporal dataframe
temp1 <- print(temp0) # is this necessary?
dataframe <- rbind(dataframe,temp1) # sum new data with previous data
}
print(dataframe) # Necessary for asigning an object name
}
Ahora con f.clean puedo convertir la imagen PET en un data.frame del que exportar los ‘boundaries’ o fronteras.
#Example(s):
mask_data = f.clean("maskedPETdataEXAMPLE")
sample_30 <- mask_data[mask_data$z==30,-1]
head(sample_30)
De modo que ahora mismo tengo los valores ordenados en el z=30 de mi sample. Ahora a convertirlo a un dataset binario y sabremos simplemente dónde sí y dónde no hay datos en el z=30, a partir de ese dataset habrá que sacar las coordenadas de los valores frontera.
IMPORTANTE: Error cometido aquí. Hay que tener cuidado de, al hacer una substitución condicional en el data.frame, solo convierta los datos del PET que no sean iguales a 0 en un 1, no toda la fila incluyendo las coordenadas. Para eso parece que la función within es mejor que un condicional o un lapply.
Esto tarda la leche de tiempo (f.clean tiene que organizar medio millón de datos por participante) así que para resumir el Markdown voy a cargar directamente los datos de los que estamos hablando y que ya he exportado previamente.
load("C:/Hugo/Sites/messy-dataset/content/post/triangulation_for_neuroimage/sample_30.RData") # Lo dicho, por ahorrar tiempo
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.6.3
##
## Attaching package: 'ggplot2'
## The following object is masked from 'package:memisc':
##
## syms
## The following object is masked from 'package:plotfunctions':
##
## alpha
qplot(sample_30$pet,
geom="histogram",
main = "Histogram for PET activity in Sample_30",
xlab = "PET levels",
fill=I("blue"),
col=I("red"),
alpha=I(.2))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Bien, como se ve en el histograma, ahora mismo solamente hay valores de 0 o de no 0 pero separados en dos bloques bien distintos. A los efectos tampoco importa demasiado pero está bien ver que no hay datos intermedios que puedan llevar a confusión, o sea, artefactos.
El siguiente paso es más complicado, necesito las coordenadas de los valores frontera, o sea, aquellos valores que son los últimos ceros? o los últimos unos? He pensado en hacer una programación con for’s y if’s en la cual cada pixel “mira” encima, debajo, izquierda, y derecha y si está solamente rodeado de otros como él, no es pixel frontera; si está rodeado en parte por iguales y en parte por distintos sería frontera (?). La verdad es que esta clase de funciones son muy liosas así que voy a tratar de beber de funciones pre-programadas siempre que me sea posible:
En efecto, hacer este proceso mediante una serie de IF’s y de FOR’s es bastante problemático y más bien me da más preguntas que respuestas. ¿Qué hacer entonces? Por lo visto hay un paquete llamado contoureR que debería funcionar para mi caso aunque esa no sea su función habitual. En principio lo único que necesitaría sería un data.frame reducido en el que solo tengo tres columnas: una con las coordenadas de X, otra con las coordenadas de Y, y por último la de los valores de actividad cerebral. Sin embargo me ha dado muchísimos problemas inesperados que dejaré aquí registrados (los que recuerde) por si se repiten en un futuro no perder tanto tiempo en trouble-shooting.
- En primer lugar cualquier tipo de problemas relacionados con el rownames, por ejemplo, lleva a errores que a menudo hacen que Rstudio crashee. Por lo tanto lo primero que he hecho ha sido exportar de nuevo mis datos de prueba a un CSV e importarlos de vuelta. Mucho cuidado a la hora de hacer el import, ver anotaciones:
dat<-read.csv2("dat",
header=T, # 'Dat' tiene encabezados
stringsAsFactors=FALSE, # 'Dat' está formado por datos numéricos,
# si no se deja esto marcado lleva a muchos problemas con los decimales
sep=";", # El separador es el punto y coma
dec = "," # Los decimales se marcan con comas, no con puntos, lo cual me ha llevado a bastantes quebraderos
)
head(dat,10)
## X x y pet
## 1 1 1 1 0
## 2 2 1 2 0
## 3 3 1 3 0
## 4 4 1 4 0
## 5 5 1 5 0
## 6 6 1 6 0
## 7 7 1 7 0
## 8 8 1 8 0
## 9 9 1 9 0
## 10 10 1 10 0
- Se ve que hay que hacer algunos pequeños cambios extra sobre la estructura del data.frame, algunos ya se han obviado porque se hicieron antes de importar pero quedan aquí en el script para el futuro:
dat<-dat[,-1] # la primera columna no es útil
# names=c("x","y","pet")
# colnames(dat)=names
str(dat)
## 'data.frame': 7505 obs. of 3 variables:
## $ x : int 1 1 1 1 1 1 1 1 1 1 ...
## $ y : int 1 2 3 4 5 6 7 8 9 10 ...
## $ pet: int 0 0 0 0 0 0 0 0 0 0 ...
- La estructura todavía no es la adecuada para trabajar con el paquete ContoureR:
dat$x <- as.numeric(dat$x)
dat$y <- as.numeric(dat$y)
# dat$pet <- gsub(",",".",dat$pet) #Esto ya no hace falta tras el import detallado
dat$pet <- as.numeric(dat$pet)
str(dat)
## 'data.frame': 7505 obs. of 3 variables:
## $ x : num 1 1 1 1 1 1 1 1 1 1 ...
## $ y : num 1 2 3 4 5 6 7 8 9 10 ...
## $ pet: num 0 0 0 0 0 0 0 0 0 0 ...
- Otros detalles que hay que pulir normalmente:
# dat[is.na(dat)] <- 0
sum(is.na(dat$pet))
## [1] 0
# dat$pet<-round(dat$pet,0) # Para reducir carga en la memoria RAM
library(ggplot2)
library(cowplot)
## Warning: package 'cowplot' was built under R version 3.6.3
##
## ********************************************************
## Note: As of version 1.0.0, cowplot does not change the
## default ggplot2 theme anymore. To recover the previous
## behavior, execute:
## theme_set(theme_cowplot())
## ********************************************************
##
## Attaching package: 'cowplot'
## The following object is masked from 'package:imager':
##
## draw_text
p1 <- ggdraw() + draw_image("abort.png", scale = 1.5) + draw_label("(!!) Mensaje que he visto medio millón de veces para sacar la triangulación",x=0.5,y=0.95)
plot_grid(p1)
- Bueno, este error es muy recurrente y me ha tocado las narices durante días. Parece ser que, por alguna razón, este paquete tiene funciones que demandan mucha memoria RAM del ordenador y, cuando RStudio detecta que no tiene suficiente RAM simplemente aborta sesión y se cierra. Para solucionar eso, o al menos para garantizar que estamos usando el máximo de memoria RAM, hay que hacer:
memory.size(max = TRUE)
## [1] 327.94
- FINALMENTE: Tras una serie de pasos y transformaciones que parecen fáciles ahora pero que me han ido trayendo por el camino de la amargura a mi y a la comunidad de StackOverflow, esto debería funcionar:
library(ggplot2) # Para la visualizacion
library(contoureR)
rownames(dat) <- NULL # Quitar la indexacion a las filas, por alguna razón esto también me ha dado problemas
df = getContourLines(dat[1:7504,], # por último, seleccionar todas las filas menos una (!!!)
levels=c(0)) # Podría especificar otros saltos pero con c(0) se busca el salto entre 0 y el resto.
ggplot(df,aes(x,y,colour=z)) + geom_path()
#install.packages(c("cowplot","magick"))
library(magick)
library(cowplot)
library(ggplot2)
p1 <- ggdraw() + draw_image("brain_image.jpg")
p2 <- ggdraw() + draw_image("contour.png")
plot_grid(
p1, p2,
ncol = 2)
Ahora la cuestión es trabajar con las coordenadas que hemos obtenido a partir de la función getContourLines. Primero vamos a analizar el output obtenido, cambiarle el nombre a algo más intuitivo y entender el resultado que tenemos entre manos:
contour30=df
head(contour30)
## LID GID PID Group x y z
## 1 0 0 0 0-0 4 39 0
## 2 0 0 1 0-0 5 38 0
## 3 0 0 2 0-0 5 38 0
## 4 0 0 3 0-0 5 37 0
## 5 0 0 4 0-0 5 37 0
## 6 0 0 5 0-0 5 37 0
str(contour30)
## 'data.frame': 593 obs. of 7 variables:
## $ LID : num 0 0 0 0 0 0 0 0 0 0 ...
## $ GID : num 0 0 0 0 0 0 0 0 0 0 ...
## $ PID : num 0 1 2 3 4 5 6 7 8 9 ...
## $ Group: Factor w/ 3 levels "0-0","0-1","0-2": 1 1 1 1 1 1 1 1 1 1 ...
## $ x : num 4 5 5 5 5 5 5 5 5 6 ...
## $ y : num 39 38 38 37 37 37 36 36 35 34 ...
## $ z : num 0 0 0 0 0 0 0 0 0 0 ...
summary(contour30)
## LID GID PID Group x
## Min. :0 Min. :0.0000 Min. : 0.0 0-0:556 Min. : 4.00
## 1st Qu.:0 1st Qu.:0.0000 1st Qu.:111.0 0-1: 14 1st Qu.:15.00
## Median :0 Median :0.0000 Median :259.0 0-2: 23 Median :40.00
## Mean :0 Mean :0.1012 Mean :260.8 Mean :40.53
## 3rd Qu.:0 3rd Qu.:0.0000 3rd Qu.:407.0 3rd Qu.:65.00
## Max. :0 Max. :2.0000 Max. :555.0 Max. :76.00
## y z
## Min. : 5.00 Min. :0
## 1st Qu.:21.00 1st Qu.:0
## Median :52.00 Median :0
## Mean :48.78 Mean :0
## 3rd Qu.:75.00 3rd Qu.:0
## Max. :91.00 Max. :0
Finalmente tenemos el objeto data.frame contour30 que va a incluir los siguientes datos:
-
LID: Un número que representa el nivel. Entiendo que el valor en el nivel ya que para todo tenemos 0 así que básicamente nos indica que las coordenadas en esa fila corresponden a los valores de 0.
-
GID: Dentro de cada nivel, o sea, dentro del contour para 0, los diferentes grupos de contours. Es de esperar por el análisis de nuestra imagen que tengamos 3 grupos ya que hay un contour base y dos agujeros adicionales. Se comprueba que en efecto Contour30$GID presenta tres valores: 1, 2, y 3. Estos serán los valores a través de los que tenemos que filtrar coordenadas de nuestros contours.
-
PID: Segmento dentro del nivel GID. Esto en principio no no es útil…en principio.
-
x,y: Son las coordenadas del boundary en cuestión.
-
z: El valor en esa coordenada (0 en nuestro caso).
-
Group: Un identificador del salto que marca ese boundary.
Bueno, el caso es que necesito sacar las coordenadas de cada grupo o GID, hacer un data.frame, y luego una lista de data.frames. Para que sea como:
str(BMP)
## List of 3
## $ bound:'data.frame': 95 obs. of 2 variables:
## ..$ V1: num [1:95] -3.62 -3.62 -3.62 -3.62 -3.62 ...
## ..$ V2: num [1:95] 14.4 14.4 14.4 14.4 14.4 ...
## $ H1 :'data.frame': 21 obs. of 2 variables:
## ..$ V1: num [1:21] -3.62 -3.62 -3.61 -3.61 -3.61 ...
## ..$ V2: num [1:21] 14.3 14.3 14.3 14.3 14.3 ...
## $ H2 :'data.frame': 11 obs. of 2 variables:
## ..$ V1: num [1:11] -3.61 -3.61 -3.61 -3.61 -3.61 ...
## ..$ V2: num [1:11] 14.3 14.3 14.3 14.3 14.3 ...
Vamos entonces primero a hacer la función:
contour<-function(x){
aa<-contour30[contour30$GID==x,] # Nos quedamos con el GID==x (1,2,3)
a<-aa[,5:6] # De esas, nos quedamos solamente con las columnas de las coordenadas
print(a) # E imprimimos el resultado para luego hacer una lista
}
Y una ve definida la función hacemos un FOR LOOP:
head(coord[[1]],10)
## x y
## 1 4 39
## 2 5 38
## 3 5 38
## 4 5 37
## 5 5 37
## 6 5 37
## 7 5 36
## 8 5 36
## 9 5 35
## 10 6 34
head(coord[[2]])
## x y
## 1 40 56
## 2 40 56
## 3 40 56
## 4 40 57
## 5 40 57
## 6 40 57
head(coord[[3]])
## x y
## 1 62 59
## 2 62 59
## 3 62 59
## 4 62 60
## 5 62 60
## 6 63 61
Brain Triangulation:
Comprobamos que efectivamente los boundaries corresponden con lo esperable de un corte de un cerebro en axial.
plot(coord[[1]])
A continuación las triangulaciones propuestas. La n=8 propuesta por el paquete parece bastante deficiente, aunque habrá que evaluar luego sobre la marcha un equilibrio entre calidad y tiempo de procesamiento. Un ejemplo máximo sería el n=100 que sería probablemente excesivo.
IMPORTANTE: Se ve que los agujeros están siendo detectados.
# An integer parameter controlling the fineness of the triangulation and subsequent triangulation. As n increases the fineness increases. Usually, n = 8 seems to be a good choice.
VT8=TriMesh(coord[[1]],8,list(as.matrix(coord[[2]]),as.matrix(coord[[3]])))
head(VT8$V,10);head(VT8$Tr,10)
## X Y
## [1,] 13 32
## [2,] 13 41
## [3,] 13 50
## [4,] 13 59
## [5,] 22 14
## [6,] 22 23
## [7,] 22 32
## [8,] 22 41
## [9,] 22 50
## [10,] 22 59
## node1 node2 node3
## [1,] 1 54 57
## [2,] 1 57 60
## [3,] 1 60 62
## [4,] 1 62 63
## [5,] 1 63 64
## [6,] 1 64 66
## [7,] 1 66 6
## [8,] 1 6 7
## [9,] 1 7 2
## [10,] 1 2 54
VT15=TriMesh(coord[[1]],15,list(as.matrix(coord[[2]]),as.matrix(coord[[3]])))
head(VT15$V,10);head(VT15$Tr,10)
## X Y
## [1,] 8.8 33.8
## [2,] 8.8 38.6
## [3,] 8.8 43.4
## [4,] 8.8 48.2
## [5,] 8.8 53.0
## [6,] 13.6 24.2
## [7,] 13.6 29.0
## [8,] 13.6 33.8
## [9,] 13.6 38.6
## [10,] 13.6 43.4
## node1 node2 node3
## [1,] 1 186 188
## [2,] 1 188 189
## [3,] 1 189 191
## [4,] 1 191 194
## [5,] 1 194 7
## [6,] 1 7 8
## [7,] 1 8 2
## [8,] 1 2 186
## [9,] 2 459 460
## [10,] 2 460 184
VT25=TriMesh(coord[[1]],25,list(as.matrix(coord[[2]]),as.matrix(coord[[3]])))
head(VT25$V,10);head(VT25$Tr,10)
## X Y
## [1,] 6.88 36.68
## [2,] 6.88 39.56
## [3,] 6.88 42.44
## [4,] 6.88 45.32
## [5,] 6.88 48.20
## [6,] 6.88 51.08
## [7,] 9.76 28.04
## [8,] 9.76 30.92
## [9,] 9.76 33.80
## [10,] 9.76 36.68
## node1 node2 node3
## [1,] 1 528 529
## [2,] 1 529 530
## [3,] 1 530 531
## [4,] 1 531 532
## [5,] 1 532 9
## [6,] 1 9 10
## [7,] 1 10 2
## [8,] 1 2 528
## [9,] 2 803 804
## [10,] 2 804 527
VT100=TriMesh(coord[[1]],100,list(as.matrix(coord[[2]]),as.matrix(coord[[3]])))
head(VT100$V,10);head(VT100$Tr,10)
## X Y
## [1,] 4.72 38.84
## [2,] 4.72 39.56
## [3,] 4.72 40.28
## [4,] 4.72 41.00
## [5,] 4.72 41.72
## [6,] 4.72 42.44
## [7,] 4.72 43.16
## [8,] 5.44 35.24
## [9,] 5.44 35.96
## [10,] 5.44 36.68
## node1 node2 node3
## [1,] 1 8898 8899
## [2,] 1 8899 8900
## [3,] 1 8900 13
## [4,] 1 13 2
## [5,] 1 2 8898
## [6,] 2 9452 9453
## [7,] 2 9453 8898
## [8,] 2 13 14
## [9,] 2 14 3
## [10,] 2 3 9452
IMPORTANTE: ¿Cuál es el óptimo de ‘finura’ de la triangulación? Todavía por definir.
Brain Application:
Ahora toca utilizar estas coordenadas que hemos obtenido para usarlas en nuestros datos. En primer lugar voy a exportar las coordenadas como un objeto a parte para no tener que correr todo esto en un futuro para tenerlas. Claro que estas son las coordenadas en Z=30 pero, de todos modos, está bien tenerlas ya exportadas.
Vale, estas son las coordenadas para la triangulación, aunque no es la triangulación per se ya que esta requiere una serie de coordenadas de los vértices y luego qué vértices se requieren para cada triángulo en cuestión. Esto realmente ya lo hemos hecho de forma un poco indirecta, es lo que se muestra en:
head(VT15$V,10);head(VT15$Tr,10)
## X Y
## [1,] 8.8 33.8
## [2,] 8.8 38.6
## [3,] 8.8 43.4
## [4,] 8.8 48.2
## [5,] 8.8 53.0
## [6,] 13.6 24.2
## [7,] 13.6 29.0
## [8,] 13.6 33.8
## [9,] 13.6 38.6
## [10,] 13.6 43.4
## node1 node2 node3
## [1,] 1 186 188
## [2,] 1 188 189
## [3,] 1 189 191
## [4,] 1 191 194
## [5,] 1 194 7
## [6,] 1 7 8
## [7,] 1 8 2
## [8,] 1 2 186
## [9,] 2 459 460
## [10,] 2 460 184
Estos data.frames lo que nos estarían diciendo es: VT15_V son las coordenadas en el plano cartesiano que van a ser los vértices para los triángulos, son numerosos (aunque no se aprecia demasiado en el ‘head’ que muestro en código) y de ellos sale VT15_Tr que es una combinación de 3 vértices que forman cada uno de los triángulos para la regresión.
Procedemos a un ejemplo básico y aquí voy a dejar este post, ya que el objetivo primario de este era el obtener las nuevas triangulaciones a medida para mis cerebros. En un post posterior lo aplicaré a comparaciones pero esto ya es demasiado peso para RStudio de modo que voy a comprobar que está todo OK y pasamos a otro sub-tema.
Cargamos la muestra sample.csv que ya tenemos usado previamente, viene siendo la muestra de un participante en un Z=30, datos estándar como los que vamos a utilizar a lo largo de la tesis. Estos datos tienen que convertirse en una matriz, asegurarse de que las filas están sin nombrar para evitar conflictos, y finalmente separamos dos columnas que serán las coordenadas y una tercera columna que son los valores de PET.
sample <- read_csv("sample.csv", col_types = cols(pet = col_number(), x = col_integer(), y = col_integer())) # Datos 'sample'
sample <- as.matrix(sample) # Como matriz
rownames(sample)<-NULL # filas sin nombre
colnames(sample) <- c("X","Y","pet") # Nombres de las columnas
head(sample,10) # Vemos una muestra
## X Y pet
## [1,] 1 1 NA
## [2,] 1 2 NA
## [3,] 1 3 NA
## [4,] 1 4 NA
## [5,] 1 5 NA
## [6,] 1 6 NA
## [7,] 1 7 NA
## [8,] 1 8 NA
## [9,] 1 9 NA
## [10,] 1 10 NA
Z <- sample[,1:2] # Esto serán las coordenadas
head(Z)
## X Y
## [1,] 1 1
## [2,] 1 2
## [3,] 1 3
## [4,] 1 4
## [5,] 1 5
## [6,] 1 6
Ya=sample[,3] # Los datos de respuesta PET serán 'Ya' (no Y sin más porque después queremos comparar con Yb
Ya <- replace(Ya,is.na(Ya),0) # SCC no trabaja con NaN, hay que asegurar que está todo con 0's
Ahora ya tenemos definidos los nombres de coordenadas y valores de PET, lo siguiente que necesitamos son los vértices y triángulos que hemos definido para nuestros cerebros a partir de las coordenadas frontera que teníamos definidas.
Esto es tremendamente demandante de recursos del ordenador. He probado con VT100 (finura/n = 100) y mi ordenador no es capaz de trabajar con tal cantidad de datos, también he probado con VT25 y el ordenador es capaz de trabajar, pero requiere horas para sacar resultados y finalmente no es capaz de mostrarlo ya que no puede escribir estos datos por el excesivo tamaño del vector (mas de 30Gb). Al final, de forma temporal por lo menos, he utilizado los de VT8 que son los recomendados por el artículo de Wang et al. aunque a mí, a ojo de buen cubero, me parece que deja que desear.
Para mantener la nomenclatura del paquete más o menos similar hago las siguientes asignaciones:
Brain.V <- VT8[[1]]
Brain.Tr <- VT8[[2]]
head(Brain.V);head(Brain.Tr)
## X Y
## [1,] 13 32
## [2,] 13 41
## [3,] 13 50
## [4,] 13 59
## [5,] 22 14
## [6,] 22 23
## node1 node2 node3
## [1,] 1 54 57
## [2,] 1 57 60
## [3,] 1 60 62
## [4,] 1 62 63
## [5,] 1 63 64
## [6,] 1 64 66
plot(Brain.V,main="New Triangulation Points for my data")
Ahora hay toda una serie de parámetros que deben definirse para la triangulación. En primer lugar hay que trasponer la matriz de datos PET para que el número de filas de coordenadas corresponda con el número de columnas de la variable respuesta.
Ya <- t(Ya) # Trasposición de la variable respuesta
V.band=Brain.V # Los triángulos en cuestión que pueden ser distintos para la imagen A y la B, en este caso solo estamos haciendo el SCC de una sola imagen de modo que V.est,V.band,Tr.est,Tr.band serán iguales dos a dos.
V.est=V.band
Tr.band=Brain.Tr
Tr.est=Tr.band
d.est=5 # degree of spline for mean function
d.band=2 # degree of spline for SCC
r=1 # smoothing parameter
lambda=10^{seq(-6,3,0.5)} # penalty parameters
alpha.grid=c(0.10,0.05,0.01) # vector of confidence levels
Y finalmente este es el código para sacar las imágenes SCC, primero la triangulación en cuestión (el dibujo con triángulos) y luego plotear los SCC’s en función de dicha triangulación y con los márgenes superiores e inferiores para los diferentes alfas. Estas imágenes se van a mejorar en el siguiente post/sub-sección pero para el caso que nos atañe ahora mismo, es suficiente.
Lo que no sé es si se verán correctamente ya que hace falta darle a Enter en linea de comando así que quizás markdown no es capaz.
# Run one sample SCC construction one time:
SCC_sample=scc.image(Ya=Ya,Z=Z,d.est=d.est,d.band=d.band,r=r,
V.est.a=V.est,Tr.est.a=Tr.est,
V.band.a=V.band,Tr.band.a=Tr.band,
penalty=TRUE,lambda=lambda,alpha.grid=alpha.grid,adjust.sigma=TRUE)
plot(SCC_sample)
En efecto, el RMarkdown no puede correr los SCC porque pide que se presione Enter en linea de comando de modo que directamente muestro aquí las imágenes resultantes (solo para 0’05, también muestra para 0’1 y 0’01).
library(magick)
library(cowplot)
library(ggplot2)
pc1 <- ggdraw() + draw_image("mean2.png")
pc2 <- ggdraw() + draw_image("lower0052.png")
pc3 <- ggdraw() + draw_image("upper0052.png")
plot_grid(
pc1, pc2, pc3,
ncol = 1)
Juan A. Arias
Investigador Doctoral
Mis intereses incluyen Biología, Neurociencias, Bioestadística, Neuroimagen, Aprendizaje Estadístico y Comunicación Científica.