-
Notifications
You must be signed in to change notification settings - Fork 7
/
test_cifti.R
128 lines (99 loc) · 2.92 KB
/
test_cifti.R
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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
# rm(list = ls())
library(xml2)
library(gifti)
library(rgl)
library(RColorBrewer)
library(matrixStats)
library(cifti)
onefile=TRUE
gzipped=FALSE
verbose=TRUE
warn=-1
reorient=FALSE
call=NULL
fname = "~/Downloads/cifti-2_test_data/Conte69.MyelinAndCorrThickness.32k_fs_LR.dscalar.nii"
# fname = "~/Downloads/cifti-2_test_data/hcp_atlas_clean.dtseries.nii"
surf_fname = "~/Downloads/cifti-2_test_data/Conte69.L.inflated.32k_fs_LR.surf.gii"
surf = readgii(surf_fname)
make_attr = function(x, which, value) {
mapply(function(w, v) {
attr(x, which = w) = v
}, which, value, SIMPLIFY = FALSE)
return(x)
}
# xpath = switch(cifti_version,
# # "1.0" = '/CIFTI/Matrix/Volume',
# # "1" = '/CIFTI/Matrix/Volume',
# "2.0" = './BrainModel',
# "2" = './BrainModel'
# )
nim = nifti_2_hdr(fname)
# xmldata = cifti_xml_txt(fname)
data = read_cifti(fname)
n = sapply(data$BrainModel, attr, "BrainStructure")
names(data$BrainModel) = n
ind = data$BrainModel$CIFTI_STRUCTURE_CORTEX_LEFT + 1
data = data$data
doc = cifti_xml(fname)
doc = xml_ns_strip(doc)
# xml_attrs(doc)
# L = list(nim)
# return(nim)
doc_attrs = xml_attrs(doc)
cifti_version = doc_attrs["Version"]
stopifnot(cifti_version %in% c("2.0", "2"))
nodes = xml_find_all(
doc,
"/CIFTI/Matrix/MatrixIndicesMap")
run_attr = xml_attrs(nodes)
types = c("Volume", "Surface",
"Parcel", "NamedMap",
"BrainModel")
nodeset = xml_find_all(nodes, "./Volume")
vols = parse_volume(nodeset)
nodeset = xml_find_all(nodes, "./BrainModel")
verts = parse_brain_model(nodeset)
nodeset = xml_find_all(nodes, "./NamedMap")
lt = parse_named_map(nodeset)
nodeset = xml_find_all(nodes, "./Surface")
surfs = parse_surface(nodeset)
nodeset = xml_find_all(nodes, "./Parcel")
pars = parse_parcel(nodeset)
xL = L = surf$data
reg_verts = L$pointset
faces = as.vector(t(L$triangle) + 1)
surf_verts = L$pointset[faces,]
# norms = L$normals[faces,]
# rgl.triangles(x = surf_verts)
# ind = verts[[1]] + 1
# faces = L$faces
# run_verts = run_verts[faces, ]
good_data = rep(NA, length = nrow(L$pointset))
good_data[ind] = data[seq(length(ind)), 1]
cdata = good_data[faces]
# d = data[ind, 1]
# nfaces = length(faces)
# cdata = vector(length = nfaces)
# cdata[ faces %in% ind ] = d[ faces[ faces %in% ind ] ]
# dim(surf_verts)
# keep_ind = faces %in% ind
cols = brewer.pal(10, "Spectral")
# cols = brewer.pal(10, "Reds")
mypal = colorRampPalette(colors = cols)
# mypal = colorRampPalette(colors = c("red", "black"))
n = 20
# breaks = seq(min(cdata), max(cdata), length.out = n)
grid = 0.4
breaks = seq(min(cdata, na.rm = TRUE),
max(cdata, na.rm = TRUE) + grid, by = grid)
ints = cut(cdata, include.lowest = TRUE,
breaks = breaks)
ints = as.integer(ints)
# stopifnot(!any(is.na(ints)))
cols = mypal(n)[ints]
cols = scales::alpha(cols, 1)
# cols = scales::alpha(cols, 1)
rgl.open()
rgl.triangles(x = surf_verts, color = cols)
# rgl.triangles(x = surf_verts, coloo)
# rgl.triangles(x = nii_verts)