377 lines
11 KiB
Rust
377 lines
11 KiB
Rust
//! BTi/4D-Neuroimaging MEG Dataset Reader
|
|
//!
|
|
//! Main reader for BTi directory-based datasets.
|
|
|
|
use crate::{IoError, IoResult, NeuroReader};
|
|
use byteorder::{BigEndian, ReadBytesExt};
|
|
use std::fs::{self, File};
|
|
use std::io::{BufReader, Read, Seek, SeekFrom};
|
|
use std::path::{Path, PathBuf};
|
|
|
|
use super::config::BtiConfig;
|
|
use super::constants::*;
|
|
|
|
/// BTi MEG dataset reader
|
|
///
|
|
/// Reads data from 4D-Neuroimaging (BTi) MEG systems.
|
|
#[derive(Debug)]
|
|
pub struct BtiReader {
|
|
/// Path to the data directory
|
|
data_path: PathBuf,
|
|
/// Parsed configuration
|
|
config: BtiConfig,
|
|
/// Path to the main data file (PDF)
|
|
pdf_path: PathBuf,
|
|
/// Data type in the PDF file
|
|
data_type: BtiDataType,
|
|
/// Channel names (cached)
|
|
channel_names: Vec<String>,
|
|
/// Data offset in bytes (after header)
|
|
data_offset: usize,
|
|
}
|
|
|
|
impl BtiReader {
|
|
/// Open a BTi data directory or PDF file
|
|
///
|
|
/// # Arguments
|
|
///
|
|
/// * `path` - Path to either:
|
|
/// - A directory containing `config` and data files
|
|
/// - A direct path to a PDF data file (e.g., `c,rfDC`)
|
|
pub fn open(path: impl AsRef<Path>) -> IoResult<Self> {
|
|
let path = path.as_ref();
|
|
|
|
// Determine if path is directory or file
|
|
let (data_dir, pdf_path) = if path.is_dir() {
|
|
// Find the data file in the directory
|
|
let pdf_path = Self::find_data_file(path)?;
|
|
(path.to_path_buf(), pdf_path)
|
|
} else {
|
|
// Use the file directly, parent as data dir
|
|
let data_dir = path
|
|
.parent()
|
|
.map(|p| p.to_path_buf())
|
|
.unwrap_or_else(|| PathBuf::from("."));
|
|
(data_dir, path.to_path_buf())
|
|
};
|
|
|
|
// Try to load config file
|
|
let config_path = data_dir.join("config");
|
|
let mut config = if config_path.exists() {
|
|
BtiConfig::from_file(&config_path)?
|
|
} else {
|
|
// Create minimal config from PDF header
|
|
Self::config_from_pdf(&pdf_path)?
|
|
};
|
|
|
|
// Read PDF header for additional info
|
|
let (data_type, data_offset) = Self::read_pdf_header(&pdf_path)?;
|
|
|
|
// Update config from PDF if needed
|
|
if config.epoch_size == 0 {
|
|
let file_size = fs::metadata(&pdf_path)?.len();
|
|
let data_size = file_size - data_offset as u64;
|
|
let n_samples = data_size as usize / (config.n_channels * data_type.size());
|
|
config.epoch_size = n_samples / config.n_epochs.max(1);
|
|
}
|
|
|
|
// Cache channel names
|
|
let channel_names = if config.channels.is_empty() {
|
|
(0..config.n_channels)
|
|
.map(|i| format!("MEG{:03}", i + 1))
|
|
.collect()
|
|
} else {
|
|
config.channels.iter().map(|c| c.name.clone()).collect()
|
|
};
|
|
|
|
Ok(Self {
|
|
data_path: data_dir,
|
|
config,
|
|
pdf_path,
|
|
data_type,
|
|
channel_names,
|
|
data_offset,
|
|
})
|
|
}
|
|
|
|
/// Find the main data file in a BTi directory
|
|
fn find_data_file(dir: &Path) -> IoResult<PathBuf> {
|
|
// Common BTi data file patterns
|
|
let patterns = [
|
|
"c,rfDC", // Most common continuous data file
|
|
"c,rfhp", // High-pass filtered continuous
|
|
"e,rfDC", // Event-related data
|
|
"c,rf", // Generic continuous
|
|
"pdf", // Processed data file
|
|
];
|
|
|
|
for pattern in &patterns {
|
|
let path = dir.join(pattern);
|
|
if path.exists() {
|
|
return Ok(path);
|
|
}
|
|
}
|
|
|
|
// Try to find any file starting with 'c,' or 'e,'
|
|
for entry in fs::read_dir(dir)? {
|
|
let entry = entry?;
|
|
let name = entry.file_name();
|
|
let name_str = name.to_string_lossy();
|
|
if name_str.starts_with("c,") || name_str.starts_with("e,") {
|
|
return Ok(entry.path());
|
|
}
|
|
}
|
|
|
|
Err(IoError::FileNotFound(format!(
|
|
"No BTi data file found in {}",
|
|
dir.display()
|
|
)))
|
|
}
|
|
|
|
/// Create config from PDF header
|
|
fn config_from_pdf(pdf_path: &Path) -> IoResult<BtiConfig> {
|
|
let file = File::open(pdf_path)?;
|
|
let mut reader = BufReader::new(file);
|
|
|
|
// Read magic header
|
|
let mut magic = [0u8; 8];
|
|
reader.read_exact(&mut magic)?;
|
|
|
|
let _version = if &magic == PDF_MAGIC_V1 {
|
|
1
|
|
} else if &magic == PDF_MAGIC_V2 {
|
|
2
|
|
} else {
|
|
// Try to continue anyway with default assumptions
|
|
1
|
|
};
|
|
|
|
// Read header fields (big-endian)
|
|
reader.seek(SeekFrom::Start(PDF_NCHAN_OFFSET as u64))?;
|
|
let n_channels = reader.read_i16::<BigEndian>()? as usize;
|
|
|
|
reader.seek(SeekFrom::Start(PDF_NEPOCH_OFFSET as u64))?;
|
|
let n_epochs = reader.read_i32::<BigEndian>()? as usize;
|
|
|
|
reader.seek(SeekFrom::Start(PDF_EPOCH_SIZE_OFFSET as u64))?;
|
|
let epoch_size = reader.read_i32::<BigEndian>()? as usize;
|
|
|
|
reader.seek(SeekFrom::Start(PDF_SFREQ_OFFSET as u64))?;
|
|
let sfreq = reader.read_f64::<BigEndian>()?;
|
|
|
|
Ok(BtiConfig::from_pdf_header(
|
|
sfreq.max(1.0),
|
|
n_channels.max(1),
|
|
n_epochs.max(1),
|
|
epoch_size,
|
|
))
|
|
}
|
|
|
|
/// Read PDF header to determine data type and offset
|
|
fn read_pdf_header(pdf_path: &Path) -> IoResult<(BtiDataType, usize)> {
|
|
let file = File::open(pdf_path)?;
|
|
let mut reader = BufReader::new(file);
|
|
|
|
// Read magic header
|
|
let mut magic = [0u8; 8];
|
|
reader.read_exact(&mut magic)?;
|
|
|
|
// Read data type
|
|
reader.seek(SeekFrom::Start(PDF_DTYPE_OFFSET as u64))?;
|
|
let dtype_code = reader.read_i16::<BigEndian>()?;
|
|
|
|
let data_type = BtiDataType::try_from(dtype_code).unwrap_or(BtiDataType::Short);
|
|
|
|
Ok((data_type, PDF_HEADER_SIZE))
|
|
}
|
|
|
|
/// Get configuration info
|
|
pub fn config(&self) -> &BtiConfig {
|
|
&self.config
|
|
}
|
|
|
|
/// Get path to the data directory
|
|
pub fn path(&self) -> &Path {
|
|
&self.data_path
|
|
}
|
|
|
|
/// Get path to the PDF data file
|
|
pub fn pdf_path(&self) -> &Path {
|
|
&self.pdf_path
|
|
}
|
|
|
|
/// Read raw data from the PDF file
|
|
///
|
|
/// Returns data in channel-major format: [ch0_s0, ch0_s1, ..., ch1_s0, ...]
|
|
pub fn read_data(&mut self, tmin: f64, tmax: f64) -> IoResult<Vec<f64>> {
|
|
let sfreq = self.config.sfreq;
|
|
let n_channels = self.config.n_channels;
|
|
let total_samples = self.config.n_samples();
|
|
|
|
// Convert time to sample indices
|
|
let start_sample = ((tmin * sfreq).floor() as usize).min(total_samples);
|
|
let end_sample = ((tmax * sfreq).ceil() as usize).min(total_samples);
|
|
|
|
if start_sample >= end_sample {
|
|
return Ok(Vec::new());
|
|
}
|
|
|
|
let n_samples = end_sample - start_sample;
|
|
|
|
// Allocate output buffer (channel-major format)
|
|
let mut data = vec![0.0f64; n_channels * n_samples];
|
|
|
|
// Get calibration factors
|
|
let cals: Vec<f64> = if self.config.channels.is_empty() {
|
|
vec![1.0; n_channels]
|
|
} else {
|
|
self.config.channels.iter().map(|c| c.cal).collect()
|
|
};
|
|
|
|
// Open data file
|
|
let file = File::open(&self.pdf_path)?;
|
|
let mut reader = BufReader::new(file);
|
|
|
|
// Seek to start of data range
|
|
// BTi data is typically stored as: all channels for sample 0, all channels for sample 1, etc.
|
|
let sample_size = n_channels * self.data_type.size();
|
|
let data_offset = self.data_offset + start_sample * sample_size;
|
|
reader.seek(SeekFrom::Start(data_offset as u64))?;
|
|
|
|
// Read samples based on data type
|
|
match self.data_type {
|
|
BtiDataType::Short => {
|
|
for s in 0..n_samples {
|
|
for ch in 0..n_channels {
|
|
let raw_value = reader.read_i16::<BigEndian>()?;
|
|
data[ch * n_samples + s] = raw_value as f64 * cals[ch];
|
|
}
|
|
}
|
|
}
|
|
BtiDataType::Long => {
|
|
for s in 0..n_samples {
|
|
for ch in 0..n_channels {
|
|
let raw_value = reader.read_i32::<BigEndian>()?;
|
|
data[ch * n_samples + s] = raw_value as f64 * cals[ch];
|
|
}
|
|
}
|
|
}
|
|
BtiDataType::Float => {
|
|
for s in 0..n_samples {
|
|
for ch in 0..n_channels {
|
|
let raw_value = reader.read_f32::<BigEndian>()?;
|
|
data[ch * n_samples + s] = raw_value as f64 * cals[ch];
|
|
}
|
|
}
|
|
}
|
|
BtiDataType::Double => {
|
|
for s in 0..n_samples {
|
|
for ch in 0..n_channels {
|
|
let raw_value = reader.read_f64::<BigEndian>()?;
|
|
data[ch * n_samples + s] = raw_value * cals[ch];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
Ok(data)
|
|
}
|
|
|
|
/// Get number of MEG channels
|
|
pub fn n_meg_channels(&self) -> usize {
|
|
self.config
|
|
.channels
|
|
.iter()
|
|
.filter(|c| matches!(c.kind, super::config::BtiChannelKind::Meg))
|
|
.count()
|
|
}
|
|
|
|
/// Get number of EEG channels
|
|
pub fn n_eeg_channels(&self) -> usize {
|
|
self.config
|
|
.channels
|
|
.iter()
|
|
.filter(|c| matches!(c.kind, super::config::BtiChannelKind::Eeg))
|
|
.count()
|
|
}
|
|
}
|
|
|
|
impl NeuroReader for BtiReader {
|
|
fn read_header(&mut self) -> IoResult<()> {
|
|
// Header is already parsed in open()
|
|
Ok(())
|
|
}
|
|
|
|
fn sfreq(&self) -> f64 {
|
|
self.config.sfreq
|
|
}
|
|
|
|
fn n_channels(&self) -> usize {
|
|
self.config.n_channels
|
|
}
|
|
|
|
fn n_samples(&self) -> usize {
|
|
self.config.n_samples()
|
|
}
|
|
|
|
fn channel_names(&self) -> Vec<String> {
|
|
self.channel_names.clone()
|
|
}
|
|
|
|
fn read_raw_data(&mut self, tmin: f64, tmax: f64) -> IoResult<Vec<f64>> {
|
|
self.read_data(tmin, tmax)
|
|
}
|
|
}
|
|
|
|
#[cfg(test)]
|
|
mod tests {
|
|
use super::*;
|
|
|
|
#[test]
|
|
fn test_data_file_patterns() {
|
|
// Test the expected file patterns
|
|
let patterns = ["c,rfDC", "c,rfhp", "e,rfDC", "c,rf", "pdf"];
|
|
assert_eq!(patterns[0], "c,rfDC");
|
|
assert_eq!(patterns[1], "c,rfhp");
|
|
assert!(patterns[0].starts_with("c,"));
|
|
}
|
|
|
|
#[test]
|
|
fn test_time_to_sample_conversion() {
|
|
let sfreq: f64 = 1017.25;
|
|
let tmin: f64 = 0.5;
|
|
let tmax: f64 = 1.5;
|
|
let total_samples: usize = 5000;
|
|
|
|
let start_sample = ((tmin * sfreq).floor() as usize).min(total_samples);
|
|
let end_sample = ((tmax * sfreq).ceil() as usize).min(total_samples);
|
|
|
|
assert_eq!(start_sample, 508);
|
|
assert_eq!(end_sample, 1526);
|
|
}
|
|
|
|
#[test]
|
|
fn test_sample_size_calculation() {
|
|
let n_channels = 148;
|
|
|
|
let short_size = n_channels * BtiDataType::Short.size();
|
|
assert_eq!(short_size, 148 * 2);
|
|
|
|
let float_size = n_channels * BtiDataType::Float.size();
|
|
assert_eq!(float_size, 148 * 4);
|
|
}
|
|
|
|
#[test]
|
|
fn test_channel_major_indexing() {
|
|
// Test channel-major indexing pattern
|
|
let n_channels = 3;
|
|
let n_samples = 10;
|
|
|
|
let ch = 1;
|
|
let s = 5;
|
|
let idx = ch * n_samples + s;
|
|
|
|
assert_eq!(idx, 15);
|
|
}
|
|
}
|