---
title: "EEG Source Localization with PhysioEEG"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{EEG Source Localization}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

# EEG Source Localization with PhysioEEG

This vignette demonstrates EEG source localization methods available in
PhysioEEG. Source localization (also called source imaging) addresses the
inverse problem: estimating brain source activity from scalp-recorded
potentials. The EEG inverse problem is ill-posed (infinitely many source
configurations can produce the same scalp pattern), so additional constraints
are needed to obtain a unique solution.

PhysioEEG provides:

- **Forward modeling** with spherical head and simplified BEM models
- **Distributed source imaging** via eLORETA, sLORETA, and MNE
- **Beamforming** via LCMV and DICS spatial filters
- **Source-space frequency analysis** for band-specific source maps
- **Visualization** of source estimates on a 2D brain projection

## Setup

```{r setup, eval=FALSE}
library(PhysioEEG)

# Create 19-channel EEG data (10-20 system)
pe <- make_eeg(n_time = 2500, n_channels = 19, sr = 500)

# Assign standard electrode positions (required for source localization)
pe <- eegMontage(pe, system = "10-20")

# Verify positions are set
head(SummarizedExperiment::colData(pe))
```

For best results, preprocess the data before source localization:

```{r preprocess, eval=FALSE}
# Bandpass filter to remove drift and high-frequency noise
pe <- eegFilter(pe, lowcut = 1, highcut = 40)

# Re-reference to average (recommended for source imaging)
pe <- eegRereference(pe, ref_type = "average", assay_name = "filtered")
```

## Forward Model

The forward model describes how neural source activity maps to scalp
electrode potentials. It is represented by a leadfield matrix of dimensions
(n_electrodes x n_sources * 3), where each source has three orientation
columns (x, y, z dipole directions).

### Spherical Head Model

The simplest approach models the head as an infinite homogeneous conducting
medium with conductivity 0.33 S/m:

```{r forward-spherical, eval=FALSE}
# Build forward model with 100 source locations
fm_sphere <- eegForwardModel(pe, method = "spherical", n_sources = 100,
                              assay_name = "rereferenced")

# Inspect the leadfield
dim(fm_sphere$leadfield)
# [1]  19 300  (19 electrodes x 100 sources x 3 orientations)

# Source positions are distributed inside a sphere of radius 0.8
head(fm_sphere$source_positions)

# Electrode positions
print(fm_sphere$electrode_positions)
```

### Simplified BEM Model

The Berg & Scherg (1994) 3-shell approximation accounts for the different
conductivities of brain, skull, and scalp using virtual dipoles at scaled
eccentricities:

```{r forward-bem, eval=FALSE}
fm_bem <- eegForwardModel(pe, method = "bem_simplified", n_sources = 100,
                           assay_name = "rereferenced")
dim(fm_bem$leadfield)
# Same dimensions, but different values due to 3-shell correction
```

### Leadfield Interpretation

The leadfield encodes the sensitivity pattern of each electrode to each
source:

```{r forward-interpret, eval=FALSE}
# For source 1, orientation x: how much does each electrode "see"?
source1_x <- fm_sphere$leadfield[, 1]
names(source1_x) <- fm_sphere$electrode_positions$label
print(round(source1_x, 4))

# Electrodes closest to the source have the largest values
```

## Source Estimation

Source estimation uses the forward model and data to reconstruct brain
activity. PhysioEEG supports three distributed source imaging methods:

### MNE (Minimum Norm Estimate)

MNE minimizes the L2 norm of the source distribution. It provides the
minimum-energy solution but tends to favor superficial sources:

```{r source-mne, eval=FALSE}
pe_mne <- eegSourceEstimate(pe, fm_sphere, method = "mne",
                             lambda = 0.05,
                             assay_name = "rereferenced")

# Source estimates are stored as a time x (n_sources*3) matrix
src_data <- SummarizedExperiment::assay(pe_mne, "source")
dim(src_data)
# [1] 2500  300  (2500 time points x 100 sources x 3 orientations)
```

### sLORETA (Standardized Low-Resolution Brain Electromagnetic Tomography)

sLORETA normalizes MNE by the resolution matrix diagonal, providing
zero-localization error for point sources:

```{r source-sloreta, eval=FALSE}
pe_sloreta <- eegSourceEstimate(pe, fm_sphere, method = "sloreta",
                                 lambda = 0.05,
                                 assay_name = "rereferenced")
```

### eLORETA (Exact Low-Resolution Electromagnetic Tomography)

eLORETA uses iterative weight estimation to achieve exact localization for
point sources, even with noise:

```{r source-eloreta, eval=FALSE}
pe_eloreta <- eegSourceEstimate(pe, fm_sphere, method = "eloreta",
                                 lambda = 0.05,
                                 assay_name = "rereferenced")
```

### Comparing Methods

The three methods differ in their spatial resolution, depth bias, and noise
sensitivity:

```{r source-compare, eval=FALSE}
# Compute source power (RMS across orientations) for each method
.rms_power <- function(pe_src, n_sources) {
  src <- SummarizedExperiment::assay(pe_src, "source")
  power <- numeric(n_sources)
  for (j in seq_len(n_sources)) {
    cols <- ((j - 1) * 3 + 1):(j * 3)
    power[j] <- mean(rowSums(src[, cols]^2))
  }
  power
}

mne_power <- .rms_power(pe_mne, 100)
sloreta_power <- .rms_power(pe_sloreta, 100)
eloreta_power <- .rms_power(pe_eloreta, 100)

# Find peak source for each method
cat("MNE peak source:", which.max(mne_power), "\n")
cat("sLORETA peak source:", which.max(sloreta_power), "\n")
cat("eLORETA peak source:", which.max(eloreta_power), "\n")
```

## Beamformer

Beamformers are spatial filters that estimate the power at each source
location by constructing a filter that passes activity from the target
location while suppressing contributions from elsewhere.

### LCMV Beamformer

The Linearly Constrained Minimum Variance (LCMV) beamformer operates in the
time domain using the data covariance matrix:

```{r beamformer-lcmv, eval=FALSE}
pe_lcmv <- eegBeamformer(pe, fm_sphere, method = "lcmv",
                          assay_name = "rereferenced")

# Source power is stored as a n_sources x 1 matrix
bf_power <- SummarizedExperiment::assay(pe_lcmv, "beamformer")
dim(bf_power)
# [1] 100   1

# Find the source with maximum power
peak_source <- which.max(bf_power[, 1])
cat("LCMV peak source:", peak_source, "\n")
cat("Peak power:", bf_power[peak_source, 1], "\n")
```

### DICS Beamformer

Dynamic Imaging of Coherent Sources (DICS) operates in the frequency domain,
using the cross-spectral density matrix in a specified frequency band:

```{r beamformer-dics, eval=FALSE}
# Alpha-band (8-13 Hz) source localization
pe_dics <- eegBeamformer(pe, fm_sphere, method = "dics",
                          freq_range = c(8, 13),
                          assay_name = "rereferenced")

dics_power <- SummarizedExperiment::assay(pe_dics, "beamformer")
peak_dics <- which.max(dics_power[, 1])
cat("DICS alpha peak source:", peak_dics, "\n")

# Compare with beta-band source
pe_dics_beta <- eegBeamformer(pe, fm_sphere, method = "dics",
                               freq_range = c(13, 30),
                               assay_name = "rereferenced")
dics_beta <- SummarizedExperiment::assay(pe_dics_beta, "beamformer")
cat("DICS beta peak source:", which.max(dics_beta[, 1]), "\n")
```

## Source Power

`eegSourcePower()` computes frequency-band source maps from previously
estimated source time series. It requires prior source estimation via
`eegSourceEstimate()`:

```{r source-power, eval=FALSE}
# Compute band power at each source location
bp <- eegSourcePower(pe_mne, source_assay = "source")
head(bp)
# source_id, band, power

# Reshape to wide format for comparison
bp_wide <- reshape(bp, direction = "wide", idvar = "source_id",
                   timevar = "band", v.names = "power")
head(bp_wide)
```

### Custom Frequency Bands

```{r source-power-custom, eval=FALSE}
# Narrow-band source analysis
custom_bands <- list(
  mu     = c(8, 12),
  beta   = c(15, 25),
  gamma1 = c(30, 45)
)

bp_custom <- eegSourcePower(pe_mne, bands = custom_bands,
                             source_assay = "source")
head(bp_custom)
```

## Visualization

PhysioEEG provides `eegPlotSource()` for visualizing source localization
results on a 2D brain projection.

### Scatter Plot

Source locations are shown as points, sized and colored by amplitude. Only
sources above a percentile threshold are displayed:

```{r vis-scatter, eval=FALSE}
# Create source data for visualization
src_positions <- fm_sphere$source_positions
src_amplitudes <- bf_power[, 1]

src_df <- data.frame(
  x = src_positions$x,
  y = src_positions$y,
  amplitude = src_amplitudes
)

# Plot with 80th percentile threshold (show top 20% sources)
eegPlotSource(pe, source_data = src_df, method = "scatter",
              threshold_pct = 80)
```

### Flat Map

The flat map mode uses interpolation to create a continuous image of source
activity:

```{r vis-flatmap, eval=FALSE}
eegPlotSource(pe, source_data = src_df, method = "flatmap",
              threshold_pct = 70)
```

### Comparing Methods Visually

```{r vis-compare, eval=FALSE}
# LCMV result
lcmv_df <- data.frame(
  x = src_positions$x,
  y = src_positions$y,
  amplitude = bf_power[, 1]
)
p_lcmv <- eegPlotSource(pe, source_data = lcmv_df, method = "scatter",
                          threshold_pct = 80)

# DICS alpha result
dics_df <- data.frame(
  x = src_positions$x,
  y = src_positions$y,
  amplitude = dics_power[, 1]
)
p_dics <- eegPlotSource(pe, source_data = dics_df, method = "scatter",
                          threshold_pct = 80)

# Display side by side (requires patchwork or gridExtra)
print(p_lcmv)
print(p_dics)
```

## Method Comparison

### Distributed Source Imaging vs. Beamforming

| Feature | Distributed (MNE/LORETA) | Beamformer (LCMV/DICS) |
|---------|--------------------------|------------------------|
| Output | Full time course per source | Power per source |
| Multiple sources | Well-suited | Can suppress correlated sources |
| Depth bias | MNE biased to surface | Less depth bias |
| Regularization | Lambda parameter | Data covariance |
| Frequency specificity | Post-hoc via eegSourcePower | DICS: built-in |

### Choosing a Method

- **MNE**: Good starting point, interpretable, but surface-biased.
- **sLORETA**: Best for localization accuracy of focal sources; zero
  localization error for point sources in noiseless data.
- **eLORETA**: Exact localization even with noise, but more computationally
  expensive due to iterative estimation.
- **LCMV**: Effective for broadband source power estimation; robust to noise.
- **DICS**: Best for frequency-specific source analysis (e.g., alpha
  generators, beta oscillations in motor cortex).

## Best Practices

### Electrode Density

Source localization accuracy improves with more electrodes. The 19-channel
10-20 system provides limited spatial resolution. For research applications,
64 or more channels are recommended:

```{r best-density, eval=FALSE}
# For demonstration with 19 channels
pe_19 <- make_eeg(n_time = 2500, n_channels = 19, sr = 500)
pe_19 <- eegMontage(pe_19, system = "10-20")
fm_19 <- eegForwardModel(pe_19, method = "spherical", n_sources = 50)
# Limited spatial resolution with 19 channels
```

### SNR Considerations

Higher signal-to-noise ratio (SNR) improves source estimation. Strategies
to improve SNR include:

- Averaging across epochs (ERP paradigms)
- Narrowband filtering before source estimation
- Artifact rejection and ICA cleaning
- Increasing regularization (lambda) for noisy data

```{r best-snr, eval=FALSE}
# Higher regularization for noisy data
pe_noisy <- eegSourceEstimate(pe, fm_sphere, method = "sloreta",
                               lambda = 0.2,
                               assay_name = "rereferenced")

# Lower regularization for clean data
pe_clean <- eegSourceEstimate(pe, fm_sphere, method = "sloreta",
                               lambda = 0.01,
                               assay_name = "rereferenced")
```

### Regularization

The regularization parameter `lambda` controls the trade-off between data
fit and solution smoothness:

- **Small lambda** (e.g., 0.01): Better data fit, but sensitive to noise
  and may produce spurious peaks.
- **Large lambda** (e.g., 0.5): Smoother solution, more robust to noise,
  but reduced spatial resolution.
- **Typical range**: 0.01 to 0.1 for well-preprocessed data.

```{r best-regularization, eval=FALSE}
# Sweep regularization values
lambdas <- c(0.001, 0.01, 0.05, 0.1, 0.5)
for (lam in lambdas) {
  pe_src <- eegSourceEstimate(pe, fm_sphere, method = "mne",
                               lambda = lam,
                               assay_name = "rereferenced")
  src <- SummarizedExperiment::assay(pe_src, "source")
  max_power <- max(rowSums(src^2))
  cat(sprintf("lambda = %.3f -> max source power = %.2e\n", lam, max_power))
}
```

### Preprocessing Recommendations

1. **Filter**: Bandpass 1-40 Hz (or narrower for specific questions)
2. **Re-reference**: Average reference is standard for source imaging
3. **Artifact rejection**: Remove or correct eye blinks, muscle, and
   cardiac artifacts
4. **Electrode positions**: Use digitized positions when available;
   standard montage positions are an approximation
5. **Forward model**: Use BEM for more realistic head modeling when
   individual MRIs are available (not implemented in this package;
   consider MNE-Python for advanced forward modeling)
