TDA(topological data analysis)で、一次元の輪、2次元の空洞を数える
パッケージ
library(TDA)
library(foreach)
library(ggplot2)
library(tidyverse)
library(readxl)
SDの計算
下の画像のようにROIをとってImageJで各ROIのStDvを取得
全ROIの平均値をSDとする
> read_excel("C:/Users/owner/Desktop/sd.xlsx")
# A tibble: 55 x 3
Area StdDev mA
<dbl> <dbl> <dbl>
1 746 25.2 50.0
2 746 21.9 50.0
3 746 21.8 50.0
4 746 21.1 50.0
5 746 24.1 50.0
6 746 17.6 100
7 746 14.6 100
8 746 15.1 100
9 746 15.0 100
10 746 17.2 100
# ... with 45 more rows
read_excel("C:/Users/owner/Desktop/sd.xlsx") %>%
select(-Area) %>%
group_by(mA) %>%
summarise(SD = mean(StdDev)) -> SD_tbl
> SD_tbl
# A tibble: 11 x 2
mA SD
<dbl> <dbl>
1 50.0 22.8
2 100 15.9
3 150 13.2
4 200 11.2
5 250 10.1
6 300 9.58
7 350 8.46
8 400 7.99
9 450 7.68
10 500 7.32
11 550 7.13
ImageJで画像の座標を取得
Save XY Coordinatesで(x,y,z,Value)を取得
今回は管電流が50から550mAまで50mA刻みのCT水ファントム画像を用いました。だから11枚の画像
read_delim("20170930.txt",
"\t",escape_double = FALSE,
col_names = FALSE,
col_types = cols(X4 = col_double()),
trim_ws = TRUE) %>%
rename(x = V1, y = V2, z = V3, HU = V4) %>%
as.tbl() %>%
dplyr::filter(HU >= 100) %>%
group_by(z) %>%
nest() -> dat
slice_Number <- seq_along(dat$z) ## foreachで使う
> dat$data[[1]]
# A tibble: 57,275 x 3
x y HU
<int> <int> <dbl>
1 243 123 293
2 244 123 302
3 245 123 297
4 246 123 302
5 247 123 313
6 248 123 268
7 249 123 296
8 250 123 289
9 251 123 319
10 252 123 287
# ... with 57,265 more rows
TDA
TDA::alphashapeDiag()を使って計算
shapeDiag <- foreach(i = slice_Number) %do% alphaShapeDiag(dat$data[[i]])
パーシステント図(PD)をデータフレームに変換
mutate()でDeathとBirthの差を付け加える
shapeDiag_df <- foreach(i = slice_Number) %do% {
tibble::tibble(
Dimension = as.factor(shapeDiag[[i]]$diagram[,1]),
Birth = shapeDiag[[i]]$diagram[,2],
Death = shapeDiag[[i]]$diagram[,3]
) %>%
mutate(SuvivalInterval = Death - Birth)
}
PDの座標からDimension毎に個数(holes)、生存半径区間のMax,Min,Median,Meanの集計
summarise_df <- foreach(i = slice_Number, .combine = rbind) %do% {
shapeDiag_df[[i]] %>%
group_by(Dimension) %>%
summarise(.,
holes = n(),
Max = max(SuvivalInterval),
Min = min(SuvivalInterval),
Median = median(SuvivalInterval),
Mean = mean(SuvivalInterval)
) %>%
mutate(.,mAs = mAs[i])
}
生成される半径が各画像のSD以上またはSD以下で区切って集計
Birth ≦ SD
summarise_sdun_df <- foreach(i = slice_Number, .combine = rbind) %do% {
shapeDiag_df[[i]] %>%
dplyr::group_by(.,Dimension) %>%
dplyr::filter(.,Birth <= SD_tbl$SD[i]) %>%
dplyr::summarise(.,
holes = n(),
Max = max(SuvivalInterval),
Min = min(SuvivalInterval),
Median = median(SuvivalInterval),
Mean = mean(SuvivalInterval)
) %>%
mutate(.,mAs = mAs[i])
}
Birth ≧ SD
summarise_sdov_df <- foreach(i = slice_Number, .combine = rbind) %do% {
shapeDiag_df[[i]] %>%
dplyr::group_by(.,Dimension) %>%
dplyr::filter(.,Birth >= SD_tbl$SD[i]) %>%
dplyr::summarise(.,
holes = n(),
Max = max(SuvivalInterval),
Min = min(SuvivalInterval),
Median = median(SuvivalInterval),
Mean = mean(SuvivalInterval)
) %>%
mutate(.,mAs = mAs[i])
}
Dimension毎にネストさせてggplotでmAと穴or空洞の数(holes)の関係
## dimensionでネストさせたtbl_df
Interval_tbl <- summarise_df %>%
group_by(Dimension) %>%
nest()
## Birth <= SD
Interval_sdun_tbl <- summarise_sdun_df %>%
group_by(Dimension) %>%
nest()
## Birth >= SD
Interval_sdov_tbl <- summarise_sdov_df %>%
group_by(Dimension) %>%
nest()
dim <- seq_along(Interval_tbl$Dimension)
## ggplot_穴の数
gg_dimk_holes <- foreach(k = dim) %do% {
ggplot(Interval_tbl$data[[k]],aes(x = mAs,y = holes))+
xlim(0,550)+
geom_point()+
geom_smooth(se = FALSE,method = "loess") +
ggtitle(str_c("dim", sep = k-1, "_holes"))
}
gg_dimk_holes_sdun <- foreach(k = dim) %do% {
ggplot(Interval_sdun_tbl$data[[k]],aes(x = mAs,y = holes))+
xlim(0,550)+
geom_point()+
geom_smooth(se = FALSE,method = "loess") +
ggtitle(str_c("dim", sep = k-1, "_holes"),subtitle = "<= under_sd")
}
gg_dimk_holes_sdov <- foreach(k = seq_along(Interval_sdov_tbl$Dimension)) %do% {
ggplot(Interval_sdov_tbl$data[[k]],aes(x = mAs,y = holes))+
xlim(0,550)+
geom_point()+
geom_smooth(se = FALSE,method = "loess") +
ggtitle(str_c("dim", sep = k, "_holes"),subtitle = "over_sd")
}
こんな感じ
各次元(Dimension)でbirth,Death が一致した行の個数(holes)を算出することでPD作成
PD_df <- foreach(i = slice_Number) %do% {
shapeDiag_df[[i]] %>%
group_by(Dimension,Birth,Death) %>%
summarise(holes = n(),
log = log(n()),
loglog = log(log(n()))
) %>%
ungroup()
}
PD_df_nest <- foreach(i = slice_Number) %do% {
shapeDiag_df[[i]] %>%
group_by(Dimension,Birth,Death) %>%
summarise(holes = n(),
log = log(n()),
loglog = log(log(n()))
) %>%
group_by(Dimension) %>%
nest()
}
重複した点があると、ggplotではPDを書くのが遅くなる。重複が7000を超える座標が多く存在するため
PDのプロット
## PD
gg_PD <- foreach(i = slice_Number) %do% {
PD_df[[i]] %$%
ggplot(.,aes(x = Birth, y = Death, shape = Dimension, colour = Dimension)) +
geom_point() +
geom_abline(slope = 1) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) +
ggtitle(str_c("mAs",sep = i*50,"_PD"))
}
## holesは重複度が線形スケール
gg_PD_holes <- foreach(i = slice_Number) %do% {
foreach(k = dim) %do% {
PD_df_nest[[i]]$data[[k]] %$%
ggplot(.,aes(x = Birth, y = Death, colour = holes)) +
geom_point(size = rel(0.2)) +
geom_abline(slope = 1) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) +
viridis::scale_colour_viridis() +
ggtitle(str_c("mAs",sep = i*50,"_PD_holes"),subtitle = str_c("dim",sep = k))
}
}
## logは重複度がlog()スケール
gg_PD_log <- foreach(i = slice_Number) %do% {
foreach(k = dim) %do% {
PD_df_nest[[i]]$data[[k]] %$%
ggplot(.,aes(x = Birth, y = Death, colour = log)) +
geom_point(size = rel(0.1)) +
geom_abline(slope = 1) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) +
viridis::scale_colour_viridis() +
ggtitle(str_c("mAs",sep = i*50,"_PD_log"),subtitle = str_c("dim",sep = k))
}
}
## loglogは重複度loglogスケール
gg_PD_loglog <- foreach(i = slice_Number) %do% {
foreach(k = dim) %do% {
PD_df_nest[[i]]$data[[k]] %$%
ggplot(.,aes(x = Birth, y = Death, colour = loglog)) +
geom_point(size = rel(0.1)) +
geom_abline(slope = 1) +
theme_bw() +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) +
viridis::scale_colour_viridis() +
ggtitle(str_c("mAs",sep = i*50,"_PD_loglog"),subtitle = str_c("dim",sep = k))
}
}
PD
重複度holes
重複度log
重複度loglog
loglogスケールでは、log1 = 0 なので log(log(1)) = Inf になりInfの部分、すなわち重複していないPDの点はグレーで表示されます