## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  eval = FALSE,
  comment = "#>"
)
#> library(rgl)
#> knitr::knit_hooks$set(webgl = hook_webgl)

## ----setup--------------------------------------------------------------------
library(diegr)
library(rlang)

## -----------------------------------------------------------------------------
# # development version:
# # install.packages("devtools")
# devtools::install_github("gerslovaz/diegr")
# # install.packages("diegr") # after deploying to CRAN

## -----------------------------------------------------------------------------
# # in interactive boxplot of signal we choose subject,
# # channel and concrete time range,
# # in which the boxplots should be plotted
# # example for subject 1 and channel E34 in time points from 10 to 20
# boxplot_epoch(epochdata, amplitude = "signal", subject = 1, channel = "E34",
#               time_lim = 10:20)
# 
# # customizing the output title and axis labels using |> operator
# p1 <- boxplot_epoch(epochdata, amplitude = "signal", subject = 1, channel = "E34",
#                     time_lim = 10:20)
# p1 |> plotly::layout(xaxis = list(title = ""),
#                      yaxis = list(title = "Amplitude"),
#                      title = list(text = "Outlier detection for subject 1 and channel E34",
#                                   font = list(size = 11, color = "darkblue"))
#                      )
# # note: the output is not plotted to reduce the vignette file size

## -----------------------------------------------------------------------------
# # outlier detection for subject 1 and channel E34 using IQR method
# outliers_epoch(epochdata, subject = 1, sensor = "E34", method = "iqr")

## -----------------------------------------------------------------------------
# preparing data for later visualisation 
# filtering signal without outlier epochs 14 and 15
# and computing epoch level jackknife average from baseline corrected data
edata <- epochdata |>
  dplyr::filter(subject == 2, epoch %in% 1:13)
edata <- baseline_correction(edata, baseline_range = 1:10)

s1 <- compute_mean(edata, amplitude = "signal_base", time = 10,
                   level = "epoch", group = "space",
                   type = "jack")

## ----fig.alt="An average amplitude time-series plot showing the brain's electrical activity (in microvolts) over time (in milliseconds), time-locked to a stimulus event at 0 ms. The red line represents the average amplitude, and the shaded red area represents the confidence interval."----
# computing jackknife average from channel E65 
s2 <- compute_mean(edata, amplitude = "signal_base", channel = "E65",
                   level = "epoch", group = "time",
                   type = "jack")

# plotting the average, the zero time point (stimulus) is set to 10
plot_time_mean(s2, t0 = 10)

## ----fig.show='hold', fig.alt="A schematic representations of an HD-EEG sensor net, showing the arrangement of 204 labeled channels distributed across the entire scalp area together with denser point mesh created for interpolation."----
# creating a mesh with default settings
sensors204 <- unique(epochdata$sensor)
M1 <- point_mesh(template = "HCGSN256", sensor_select = sensors204)
# plot output in 2D
plot_point_mesh(M1$D2, sensor_select = sensors204)

# creating a circular mesh, only 2D coordinates
M2 <- point_mesh(dimension = 2, n = 3000, template = "HCGSN256",
                 sensor_select = sensors204, type = "circle")
# plotting a mesh - function allows different options of the result plot
plot_point_mesh(M2$D2, sensor_select = sensors204, col_sensors = "purple",
                label_sensors = TRUE, cex = 0.1)

## ----fig.show='hold', fig.alt="Fig 1: A top-down topographic map of a high-density EEG amplitude in topographic colour scale with contours and black points on sensor locations. The amplitude legend is on the left side of the scalp projection. Fig 2: The same plot without contours and legend, but with sensor labels as text displayed together with black point on sensor positions. Fig 3: A ggplot divided into three panels with topographic maps of average signal (middle) and its lower (left) and upper (right) confidence interval bounds. The colour legend position is bottom of the facets."----
# Plot average topographic map of baseline corrected signal for subject 1 
# from the time point 10 (the time of the stimulus)
# the outlier epoch 14 is extracted before computing average

# preparing data
edata <- epochdata |>
  dplyr::filter(subject == 1 & time %in% 1:10)
data_base <- baseline_correction(edata, baseline_range = 1:10)
data_mean <- compute_mean(data_base, amplitude = "signal_base", time = 10,
                          type = "point", ex_epoch = 14, group = "space")

# plotting the base topographic polygon map with contours and legend
# interval (-30,15) is selected in consideration of the signal progress
topo_plot(data = data_mean, amplitude = "average", 
          col_range = c(-30, 15), contour = TRUE)

# plotting the same map without contours and legend 
# but with sensor labels and adding the title
g1 <- topo_plot(data = data_mean, amplitude = "average", col_range = c(-30, 15),
                label_sensors = TRUE, show_legend = FALSE)
g1 + ggplot2::ggtitle("Subject 1, time of the stimulus")

# plotting the average together with CI bounds using plot_topo_mean
plot_topo_mean(data = data_mean, template = "HCGSN256", col_range = c(-30, 15))

## -----------------------------------------------------------------------------
# # plotting the scalp polygon map of previously prepared signal s1
# #open3d()
# scalp_plot(s1, amplitude = "average", col_range = c(-30, 15))
# 
# # note: the output is not plotted to reduce the vignette file size

## -----------------------------------------------------------------------------
# # filtering the example data
# s1e05 <- epochdata |> dplyr::filter(subject == 1 & epoch == 5 & time %in% 10:20)
# # Plot animation with setting the range of time:
# animate_topo(s1e05, t_lim = c(0,50))
# # Export the gif animation - this code will create the animation inside your working directory.
# # If you want to export it elsewhere, set the whole path in output_path parameter.
# animate_topo(s1e05, t_lim = c(0,50), output_path = "example_animation.gif")

## -----------------------------------------------------------------------------
# # Plot animation with default mesh and triangulation:
# animate_scalp(s1e05, amplitude = "signal")
# # export the video - the .mp4 extension is required for correct functioning
# rgl::par3d(windowRect = c(100, 100, 800, 800))
# animate_scalp(s1e05, amplitude = "signal",
#               frames_dir = "your_frames_dir_path",
#               output_path = "your_created_video_path.mp4")

## -----------------------------------------------------------------------------
# #open3d()
# scalp_plot(s1, amplitude = "average", col_range = c(-30,15))
# rgl::play3d(rgl::spin3d(axis = c(0, 0, 1), rpm = 10), duration = 10)
# 
# # note: the output is not plotted to reduce the vignette file size

## -----------------------------------------------------------------------------
# # using IQR method for selected sensor from central region
# outliers_epoch(epochdata, subject = 1, sensor = "E59", method = "iqr")
# outliers_epoch(epochdata, subject = 2, sensor = "E59", method = "iqr")
# 
# # verify results by visual inspection
# interactive_waveforms(epochdata, amplitude = "signal", subject = 1,
#                       channel= "E59", level = "epoch", t0 = 10)
# interactive_waveforms(epochdata, amplitude = "signal", subject = 2,
#                       channel= "E59", level = "epoch", t0 = 10)

## -----------------------------------------------------------------------------
# sub1clear <- epochdata |>
#   filter(subject == 1, !epoch == 14)
# sub2clear <- epochdata |>
#   filter(subject == 2, !epoch %in% c(14, 15))
# sub1base <- baseline_correction(sub1clear, baseline_range = 1:10)
# sub2base <- baseline_correction(sub2clear, baseline_range = 1:10)

## -----------------------------------------------------------------------------
# sub1mean <- compute_mean(sub1base, amplitude = "signal_base", time = 35,
#                          group = "space", level = "epoch", R = 1000)
# sub2mean <- compute_mean(sub2base, amplitude = "signal_base", time = 35,
#                          group = "space", level = "epoch", R = 1000)

## -----------------------------------------------------------------------------
# # prepare a mesh for plotting
# M <- point_mesh(dimension = 2, n = 4000, template = "HCGSN256",
#                 sensor_select = unique(sub1mean$sensor), type = "polygon")
# # compute consistent col_range across subjects
# # a) find the range of average amplitude
# c(min(sub1mean$average, sub2mean$average), max(sub1mean$average, sub2mean$average))
# # -15.634288   8.609518
# # b) make the range symmetric
# cr_subjects <- c(-16,16)
# # plot the average topographic maps with the same color range for both subjects
# plot_topo_mean(sub1mean, M, template = "HCGSN256", col_range = cr_subjects,
#                contour = TRUE)
# plot_topo_mean(sub2mean, M, template = "HCGSN256", col_range = cr_subjects,
#                contour = TRUE)

