---
title: "MRI physiological signal processing"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{MRI physiological signal processing}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  fig.width  = 7,
  fig.height = 3.5
)
```

Philips MRI scanners record cardiac, respiratory, and motion signals
alongside the BOLD acquisition via the Physiological Monitoring Unit (PMU).
mrpheus can read these `.log` files, detect QRS complexes with the
Pan-Tompkins algorithm, and compute instantaneous heart rate — producing
cardiac regressors ready for RETROICOR or HRV analysis.

```{r libs, message = FALSE}
library(mrpheus)
library(ggplot2)
library(gsignal)
```

```{r theme, include = FALSE}
theme_mri <- function() {
  theme_minimal(base_size = 13) +
    theme(
      panel.grid        = element_blank(),
      axis.line         = element_line(colour = "black"),
      axis.line.x.top   = element_blank(),
      axis.line.y.right = element_blank()
    )
}
```

---

## Reading a Philips PMU log

`read_philips_physlog()` parses any Philips PMU `.log` file and returns a
`mrpheus_physlog` object containing the signal matrix, event markers, and
recording metadata.

```{r read}
path <- system.file("extdata", "example_physlog.log", package = "mrpheus")
rec  <- read_philips_physlog(path)
rec
```

The `system` argument sets the hardware preset:

```{r system, eval = FALSE}
# Wireless VCG (Achieva/Ingenia dStream) — default
rec <- read_philips_physlog(path, system = "wBTU")

# Older wired ECG system (500 Hz)
rec <- read_philips_physlog(path, system = "wired")

# Site-specific sampling rate
rec <- read_philips_physlog(path, system = "custom", sfreq = 400)
```

The resolved sampling frequency is always at `rec$HDR$sfreq` so nothing
downstream needs to hardcode it.

### Signal channels

```{r channels}
ecg  <- as.double(rec$C[, "v1raw"])
ppu  <- as.double(rec$C[, "ppu"])
resp <- as.double(rec$C[, "resp"])
```

The wBTU layout provides ten channels:

| Channel | Contents |
|---|---|
| `v1raw`, `v2raw` | Raw VCG electrode voltages |
| `v1`, `v2` | Processed VCG signals |
| `ppu` | Peripheral pulse unit (PPG, foot) |
| `resp` | Respiratory belt |
| `gx`, `gy`, `gz` | 3-axis accelerometer |
| `mark` | Bit-encoded event marker |

```{r plot-signals, fig.height = 5, echo = FALSE}
sfreq <- rec$HDR$sfreq
t     <- seq_along(ecg) / sfreq

df_signals <- data.frame(
  t       = rep(t, 3),
  value   = c(ecg, ppu, resp),
  channel = rep(c("ECG (v1raw)", "PPU", "RESP"), each = length(t))
)
df_signals$channel <- factor(df_signals$channel,
                              levels = c("ECG (v1raw)", "PPU", "RESP"))

ggplot(df_signals, aes(t, value)) +
  geom_line(linewidth = 0.3, colour = "#014370") +
  facet_wrap(~channel, ncol = 1, scales = "free_y") +
  labs(x = "Time (s)", y = "Amplitude (a.u.)") +
  theme_mri() +
  theme(strip.text = element_text(face = "bold", size = 11))
```

### Event markers

```{r markers}
rec$I$ScannerStart
rec$I$ScannerStop
head(rec$I$VcgOnset)
```

---

## QRS detection

`detect_qrs()` implements the Pan-Tompkins (1985) algorithm and returns a
`mrpheus_qrs` object.

```{r qrs}
qrs <- detect_qrs(ecg, fs = rec$HDR$sfreq)
qrs
```

R-peak indices are in `qrs$qrs_i`. Account for `qrs$delay` if you need
indices aligned to the raw (pre-filter) ECG.

```{r plot-peaks, echo = FALSE}
bp     <- gsignal::butter(3L, c(5, 15) * 2 / sfreq, type = "pass")
ecg_bp <- as.double(gsignal::filtfilt(bp, ecg))

df_bp <- data.frame(t = t, ecg = ecg_bp)
df_rpeaks <- data.frame(t = qrs$qrs_i / sfreq, amp = qrs$qrs_amp)

ggplot(df_bp, aes(t, ecg)) +
  geom_line(linewidth = 0.3, colour = "#1B6799") +
  geom_point(data = df_rpeaks, aes(x = t, y = amp),
             colour = "#FC544A", size = 2, shape = 21,
             fill = "#FC544A", alpha = 0.85) +
  labs(x = "Time (s)", y = "Amplitude (normalised)",
       title = "Bandpass ECG with detected R-peaks") +
  theme_mri()
```

```{r rr}
rr_ms <- diff(qrs$qrs_i) / rec$HDR$sfreq * 1000
cat("Mean RR :", round(mean(rr_ms)), "ms\n")
cat("SDNN    :", round(sd(rr_ms)), "ms\n")
```

### PPU as an alternative

In neonates, PPU often gives a cleaner signal than the ECG electrodes.
`detect_qrs()` accepts any pulsatile channel:

```{r ppu, eval = FALSE}
qrs_ppu <- detect_qrs(as.double(rec$C[, "ppu"]), fs = rec$HDR$sfreq)
```

---

## Heart rate signal

`compute_hr_signal()` converts R-peak indices into a sample-by-sample
instantaneous HR trace (bpm). It reads `fs` from the `mrpheus_qrs` object
directly:

```{r hr}
hr <- compute_hr_signal(qrs)
```

```{r plot-hr, echo = FALSE}
df_hr <- data.frame(t = seq_along(hr) / sfreq, hr = hr)

ggplot(subset(df_hr, hr > 0), aes(t, hr)) +
  geom_line(linewidth = 0.4, colour = "#FFA75D") +
  labs(x = "Time (s)", y = "Heart rate (bpm)",
       title = "Instantaneous heart rate") +
  theme_mri()
```

---

## Epoch extraction

To extract the physiology aligned to the fMRI scan, use the scanner-stop
marker. This is more reliable than the start marker because the exact
recording onset varies but the stop is always clean:

```{r epoch}
sfreq    <- rec$HDR$sfreq
enddelay <- 42L
tr       <- 0.392
vols     <- 47L

stopidx  <- tail(rec$I$ScannerStop, 1L) - enddelay
startidx <- max(1L, stopidx - round(vols * tr * sfreq))
epoch    <- rec$C[startidx:stopidx, ]
cat("Epoch:", nrow(epoch), "samples /", round(nrow(epoch) / sfreq, 1), "s\n")
```

### RR intervals per TR

```{r rr-tr}
qrs_epoch <- detect_qrs(as.double(epoch[, "v1raw"]), fs = sfreq)

rr_t <- qrs_epoch$qrs_i[-1] / sfreq
rr_s <- diff(qrs_epoch$qrs_i) / sfreq

tr_times <- seq(0, vols * tr, by = tr)
rr_tr    <- approx(rr_t, rr_s, xout = tr_times, method = "constant")$y

cat("RR per TR (first 5):", round(rr_tr[1:5] * 1000), "ms\n")
```

`rr_tr` is a length-`vols` vector ready for RETROICOR phase computation or
inclusion as a nuisance regressor in a GLM.

---

## Full pipeline

```{r pipeline, eval = FALSE}
library(mrpheus)

rec      <- read_philips_physlog("sub-01_ses-01_physlog.log")
sfreq    <- rec$HDR$sfreq
stopidx  <- tail(rec$I$ScannerStop, 1L) - 42L
startidx <- stopidx - round(2300L * 0.392 * sfreq)
epoch    <- rec$C[startidx:stopidx, ]
qrs      <- detect_qrs(as.double(epoch[, "v1raw"]), fs = sfreq)
hr       <- compute_hr_signal(qrs)
rr_tr    <- approx(
  x      = qrs$qrs_i[-1] / sfreq,
  y      = diff(qrs$qrs_i) / sfreq,
  xout   = seq(0, 2300 * 0.392, by = 0.392),
  method = "constant"
)$y
```

---

## References

Pan J, Tompkins WJ. A real-time QRS detection algorithm.
*IEEE Trans Biomed Eng.* 1985;32(3):230–236.
<https://doi.org/10.1109/TBME.1985.325532>

Groot PFC. *ReadPhilipsScanPhysLog* (MATLAB). AMC Amsterdam. 2010.

Glover GH, Li TQ, Ress D. Image-based method for retrospective correction of
physiological motion effects in fMRI: RETROICOR. *Magn Reson Med.*
2000;44(1):162–167.
