SampleInterpolateValues3D.m

%% SampleInterpolateValues3D
% Demonstrates the capabilities of pf2.probe.plot.interpolateValues3D
% for 3D brain surface visualization with interpolated fNIRS data.
%
% Sections:
%   1-2.   Basic device plotting (fNIR1000, Hitachi 3x5)
%   3.     Interpolation types (quadratic vs linear)
%   4.     Multiple probes on one brain
%   5.     10-20 reference labels
%   6.     Custom label styling
%   7.     Brain styling (wireframe, low-res)
%   8.     Split colormap with custom range
%   9.     Colorbar title
%  10.     Log scale
%  11.     Buffer distance comparison
%  12.     Camera positions
%  13.     NaN channel masking
%  14.     Dark background styling
%  15-16.  EEG montage / fNIRS over EEG
%  17-20.  Optical scattering pathways
%  21-25.  Brodmann areas
%  26-27.  Voxel brain basics
%  28-31.  Voxel brain lighting styles
%  32-35.  Advanced visualization
%  36-37.  ROI definition and plotting
%  38-42.  Misc features
%  43-45.  Coordinate transformation (CapTrak/subject-space to MNI)
%
% See also: SampleInterpolateValues3D_Animation

%% 1. Basic device plotting — fNIR1000
% Pass config name as second arg to specify device without globals
figure;
pf2.probe.plot.interpolateValues3D([], {'fNIR_Devices_fNIR1000.cfg'});

%% 2. Basic device plotting — Hitachi 3x5
figure;
pf2.probe.plot.interpolateValues3D([], {'Hitachi_ETG4000_3x5.cfg'});

%% 3. Interpolation types
% Compare quadratic (default) vs linear interpolation
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];

figure;
subplot(1, 2, 1);
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'ChannelLabels', false, 'SDLabels', false);
title('Quadratic (default)');

subplot(1, 2, 2);
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'ChannelLabels', false, 'SDLabels', false, ...
    'InterpolateType', 'linear');
title('Linear');

%% 4. Multiple probes on one brain
% Combine two different devices on the same 3D brain
hitachi35 = pf2.import.sampleData.Hitachi_ETG4000_3x5();
fnir1200 = pf2.import.sampleData.fNIR1200();

figure;
pf2.probe.plot.interpolateValues3D( ...
    {1:16, 1:22}, {fnir1200, hitachi35}, ...
    'SDLabels', true, 'showReference', false, ...
    'bufferDistance', 20, 'cmap', 'cool');

%% 5. 10-20 reference labels
% Overlay international 10-20 landmark labels
figure;
pf2.probe.plot.interpolateValues3D(1:22, {'Hitachi_ETG4000_3x5.cfg'}, ...
    'I1020_labels', {'T7', 'CPz'});

%% 6. Custom label styling
% Change source/detector sphere colors and font color
figure;
pf2.probe.plot.interpolateValues3D(1:16, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'labelspherecolors', ["k", "b"], 'labelfontcolor', 'g');

%% 7. Brain styling
% White brain with black wireframe, low-res mesh
figure;
pf2.probe.plot.interpolateValues3D(1:16, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'brainColor', [1 1 1], 'brainLineColor', 'k', 'useHighRes', false);

%% 8. Split colormap with custom range
% Use different colormaps above and below a threshold
figure;
pf2.probe.plot.interpolateValues3D(-7:8, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'minVal', [-2, 2], 'maxVal', 6, ...
    'cmap', 'autumn', 'cmap_lower', 'winter', ...
    'ChannelLabels', false, 'SDLabels', false);

%% 9. Colorbar title
figure;
pf2.probe.plot.interpolateValues3D(1:22, {'Hitachi_ETG4000_3x5.cfg'}, ...
    'interpolateType', 'quadratic', 'colorbarStr', "HbO");

%% 10. Log scale
% Useful for p-value maps or data spanning orders of magnitude
figure;
data = [ones(1, 9) 1e-6 1e-5 1e-2 ones(1, 10)];
pf2.probe.plot.interpolateValues3D(data, {'Hitachi_ETG4000_3x5.cfg'}, ...
    'logScale', true, 'interpolateType', 'quadratic', ...
    'minVal', 1e-8, 'maxVal', 0.9);

%% 11. Buffer distance comparison
% Larger buffer = more brain surface colored around each optode
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];

figure;
subplot(1, 2, 1);
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'ChannelLabels', false, 'SDLabels', false, ...
    'InterpolateType', 'linear', 'bufferDistance', 25, ...
    'showColorbar', false);
title('Buffer 25');

subplot(1, 2, 2);
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'ChannelLabels', false, 'SDLabels', false, ...
    'InterpolateType', 'linear', 'bufferDistance', 50, ...
    'showColorbar', false);
title('Buffer 50');

%% 12. Camera positions
% Compare front and left views on multi-probe setup
hitachi35 = pf2.import.sampleData.Hitachi_ETG4000_3x5();
fnir1200 = pf2.import.sampleData.fNIR1200();

figure;
subplot(2, 1, 1);
pf2.probe.plot.interpolateValues3D( ...
    {1:16, 1:22}, {fnir1200, hitachi35}, ...
    'initCamPosition', 'front', ...
    'ChannelLabels', false, 'SDLabels', false);
title('Front');

subplot(2, 1, 2);
pf2.probe.plot.interpolateValues3D( ...
    {1:16, 1:22}, {fnir1200, hitachi35}, ...
    'initCamPosition', 'left', ...
    'ChannelLabels', false, 'SDLabels', false);
title('Left');

%% 13. NaN channel masking
% Channels set to NaN are excluded from interpolation.
% Useful for highlighting only significant channels (e.g. from stats).

% All channels have data
allData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];

% Only channels 3, 7, 10 are "significant" — rest are masked out
maskedData = nan(1, 16);
maskedData([3, 7, 10]) = allData([3, 7, 10]);

figure;
subplot(1, 2, 1);
pf2.probe.plot.interpolateValues3D(allData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'ChannelLabels', true, 'SDLabels', false, ...
    'showColorbar', false);
title('All channels');

subplot(1, 2, 2);
pf2.probe.plot.interpolateValues3D(maskedData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'ChannelLabels', true, 'SDLabels', false, ...
    'showColorbar', false);
title('Only significant channels');

%% 14. Dark background styling
% Dark figure background with light brain wireframe
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];

figure('Color', [0.1 0.1 0.1]);
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'brainColor', [0.3 0.3 0.3], ...
    'brainLineColor', [0.4 0.4 0.4], ...
    'backgroundColor', [0.1 0.1 0.1], ...
    'labelfontcolor', 'w', ...
    'labelspherecolors', ["w", "c"], ...
    'cmap', 'hot', ...
    'ChannelLabels', true, 'SDLabels', false);

%% 15. Full EEG montage
% Plot data at standard 10-20 positions across the whole scalp
figure;
data = [-2.7, -2.8, -2.6,...
        -3.5, -2.4, -1.8, -2.2, -1.7, -1.6, -1.1, -1.8, -1.9,...
        -1.7, -1.0, -0.8, -0.7, -0.5, -0.4, -1.2,...
        -3.0, -1.0,  0.3,  0.8,  0.8,  1.2,  1.0,  0.4, -1.0,...
        -1.4,  0.2,  1.3,  2.5,  3.0,  2.8,  2.4,  2.0,  0.0,...
         0.0,  0.9,  2.3,  2.5,  3.0,  2.8,  3.4,  3.0,  1.6,...
         2.4,  1.8,  2.2,  1.5, -1.0,...
         1.0,  0.0,  1.0];
labels = {'Fp1', 'Fpz', 'Fp2',...
    'F7', 'F5', 'F3', 'F1', 'Fz', 'F2', 'F4', 'F6', 'F8',...
    'FC5', 'FC3', 'FC1', 'FCz', 'FC2', 'FC4', 'FC6', ...
    'T7', 'C5', 'C3', 'C1', 'Cz', 'C2', 'C4', 'C6', 'T8',...
    'TP7', 'CP5', 'CP3', 'CP1', 'CPz', 'CP2', 'CP4', 'CP6', 'TP8',...
    'P7', 'P5', 'P3', 'P1', 'Pz', 'P2', 'P4', 'P6', 'P8', ...
    'PO7', 'PO3', 'POz', 'PO4', 'PO8', ...
    'O1', 'Oz', 'O2'};
pf2.probe.plot.interpolateValues3D(data, ...
    'useEEG', true, 'I1020_labels', labels, ...
    'cmap', 'jet', 'interpolateType', 'linear', 'bufferDistance', 25);

%% 16. fNIRS device with 10-20 reference landmarks
% Show fNIRS probe with nearby 10-20 positions labeled for reference
figure;
fnir1200 = pf2.import.sampleData.fNIR1200();
pf2.probe.plot.interpolateValues3D(1:16, fnir1200, ...
    'cmap', 'jet', ...
    'I1020_labels', {'C5', 'C3', 'C1', 'Cz', 'C2', 'C4', 'C6'}, ...
    'SDLabels', true);

%% 17. Optical scattering pathways — single device
figure;
pf2.probe.plot.interpolateValues3D(1:16, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'showScattering', true, 'brainAlpha', 0.25);

%% 18. Optical scattering pathways — Hitachi
figure;
pf2.probe.plot.interpolateValues3D(1:22, {'Hitachi_ETG4000_3x5.cfg'}, ...
    'showScattering', true, 'brainAlpha', 0.7, 'scatteringFactor', 0.5);

%% 19. Optical scattering — multi-device (NIRX + fNIR)
figure;
pf2.probe.plot.interpolateValues3D( ...
    {([1:16]*5/16+3)*-1, [1:48]*5/48+3}, ...
    {'fNIR_Devices_fNIR1000', 'NIRX_Sport_16x16_parietal'}, ...
    'initCamPosition', 'left', ...
    'ChannelLabels', true, 'SDLabels', true, ...
    'showScattering', true, ...
    'cmap', 'autumn', 'cmap_lower', 'winter', 'minVal', [-1, 1], ...
    'brainAlpha', 0.85, 'brainColor', [1 1 1]*0.8);

%% 20. Optical scattering — NIRX only
figure;
pf2.probe.plot.interpolateValues3D( ...
    1:48, {'NIRX_Sport_16x16_parietal'}, ...
    'initCamPosition', 'left', ...
    'ChannelLabels', true, 'SDLabels', true, ...
    'showScattering', true);

%% ========================================================================
%  BRODMANN AREAS
%  ========================================================================

%% 21. Brodmann areas — all regions on surface mesh
% Pass true to show all 55 Brodmann areas color-coded on the brain.
% Each area gets a unique color from the BA_cmp colormap, with a legend.
figure;
pf2.probe.plot.interpolateValues3D([], {'fNIR_Devices_fNIR1000.cfg'}, ...
    'BrodmannAreas', true, ...
    'ChannelLabels', false, 'SDLabels', false);

%% 22. Brodmann areas — selected regions (dorsolateral PFC)
% Pass a numeric array to highlight only specific areas.
% BA 9, 10, 46 = dorsolateral prefrontal cortex
figure;
pf2.probe.plot.interpolateValues3D([], {'fNIR_Devices_fNIR1000.cfg'}, ...
    'BrodmannAreas', [9, 10, 46], ...
    'ChannelLabels', true, 'SDLabels', false, ...
    'initCamPosition', 'front');

%% 23. Brodmann areas — with fNIRS data overlay
% Combine Brodmann area outlines with actual channel data.
% Data is interpolated on top while BA regions show underneath.
figure;
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'BrodmannAreas', [9, 10, 46], ...
    'brainAlpha', 0.8, ...
    'ChannelLabels', true, 'SDLabels', false, ...
    'initCamPosition', 'front');

%% 24. Brodmann areas — voxel mode
% Volumetric voxel rendering of Brodmann areas from the MNI atlas.
% Much heavier but shows 3D structure of each area.
figure;
pf2.probe.plot.interpolateValues3D([], {'fNIR_Devices_fNIR1000.cfg'}, ...
    'BrodmannAreas', [9, 10, 46], ...
    'useVoxelBrodmannAreas', true, ...
    'showVoxelBrain', true, ...
    'ChannelLabels', false, 'SDLabels', false);

%% 25. Brodmann areas — custom colormap
% Change the color scheme used for Brodmann area regions
figure;
pf2.probe.plot.interpolateValues3D([], {'Hitachi_ETG4000_3x5.cfg'}, ...
    'BrodmannAreas', [7, 39, 40, 22], ...
    'BA_cmp', @parula, ...
    'ChannelLabels', true, 'SDLabels', false, ...
    'initCamPosition', 'left');

%% ========================================================================
%  ADVANCED VISUALIZATION
%  ========================================================================

%% 26. Voxel brain rendering (without Brodmann)
% Show the MNI T1 volumetric brain as voxels instead of the surface mesh
figure;
pf2.probe.plot.interpolateValues3D(1:16, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'showVoxelBrain', true, ...
    'ChannelLabels', true, 'SDLabels', false);

%% 27. Voxel brain with data overlay and custom tint
% Channel data projected onto the volumetric voxel brain
figure;
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'showVoxelBrain', true, ...
    'voxelColor', [0.92, 0.68, 0.68], ...
    'ChannelLabels', true, 'SDLabels', false);

%% 28. Voxel brain — realistic wet lighting
% Smooth lighting with subtle specular sheen, like a wet brain specimen
figure('Color', 'k');
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'showVoxelBrain', true, ...
    'voxelLighting', 'realistic', ...
    'voxelColor', [0.85, 0.75, 0.75], ...
    'ChannelLabels', true, 'SDLabels', false);

%% 29. Voxel brain — dramatic lighting
% High-contrast poster-style with deep sulci shadows and sharp highlights
figure('Color', 'k');
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'showVoxelBrain', true, ...
    'voxelLighting', 'dramatic', ...
    'voxelColor', [0.8, 0.72, 0.72], ...
    'ChannelLabels', false, 'SDLabels', false);

%% 30. Voxel brain — clinical/radiologist style
% Grayscale brain on dark background with sharp nearest-neighbor projection
figure('Color', 'k');
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'showVoxelBrain', true, ...
    'voxelLighting', 'clinical', ...
    'voxelColor', [0.85, 0.85, 0.85], ...
    'backgroundColor', [0 0 0], ...
    'interpolateType', 'nearest', ...
    'ChannelLabels', false, 'SDLabels', false);

%% 31. Voxel brain — dramatic with Brodmann areas
% Dramatic lighting on Brodmann voxel areas
figure('Color', 'k');
pf2.probe.plot.interpolateValues3D(1:16, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'showVoxelBrain', true, ...
    'useVoxelBrodmannAreas', true, ...
    'BrodmannAreas', [9, 10, 46], ...
    'voxelLighting', 'dramatic', ...
    'ChannelLabels', true, 'SDLabels', false);

%% 32. Glass brain — transparent with scattering paths
% Semi-transparent brain showing light propagation through tissue
figure;
pf2.probe.plot.interpolateValues3D(1:16, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'showScattering', true, ...
    'brainAlpha', 0.15, ...
    'brainColor', [0.85 0.85 0.9], ...
    'ChannelLabels', false, 'SDLabels', true, ...
    'initCamPosition', 'left');

%% 33. All four interpolation types compared
% nearest (sharp), linear, quadratic (default), cubic (smoothest)
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];
interpTypes = {'nearest', 'linear', 'quadratic', 'cubic'};

figure;
for idx = 1:4
    subplot(2, 2, idx);
    pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
        'ChannelLabels', false, 'SDLabels', false, ...
        'showColorbar', false, ...
        'InterpolateType', interpTypes{idx});
    title(interpTypes{idx});
end

%% 34. Custom camera angle — top-down and tilted
figure;
subplot(1, 2, 1);
pf2.probe.plot.interpolateValues3D(1:22, {'Hitachi_ETG4000_3x5.cfg'}, ...
    'initCamPosition', 'top', ...
    'ChannelLabels', true, 'SDLabels', false, ...
    'showColorbar', false);
title('Top');

subplot(1, 2, 2);
pf2.probe.plot.interpolateValues3D(1:22, {'Hitachi_ETG4000_3x5.cfg'}, ...
    'initCamPosition', [500, 800, 600], ...
    'camTarget', [0, -20, 0], ...
    'ChannelLabels', true, 'SDLabels', false, ...
    'showColorbar', false);
title('Custom [500, 800, 600]');

%% 35. Diverging colormap — activation vs deactivation
% Split colormap centered at zero: warm = increase, cool = decrease
figure;
activationData = randn(1, 22) * 2;  % mix of positive and negative
pf2.probe.plot.interpolateValues3D(activationData, {'Hitachi_ETG4000_3x5.cfg'}, ...
    'cmap', 'autumn', 'cmap_lower', 'winter', ...
    'minVal', [-2, 2], 'maxVal', 4, ...
    'colorbarStr', 'HbO (\muM)', ...
    'interpolateType', 'quadratic', ...
    'ChannelLabels', true, 'SDLabels', false);

%% 36. Publication figure — white background, saved to file
% Clean white styling suitable for journal figures
figure('Color', 'w');
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'brainColor', [0.95 0.92 0.90], ...
    'backgroundColor', [1 1 1], ...
    'cmap', 'hot', ...
    'colorbarStr', 'F-stat', ...
    'ChannelLabels', true, 'SDLabels', false, ...
    'initCamPosition', 'front', ...
    'savePath', 'sample_3d_brain.png', ...
    'saveWidth', 800, 'saveHeight', 600, 'saveDPI', 300);

%% ========================================================================
%  ROI MODE
%  ========================================================================

%% 37. ROI definition and 2D ROI plotting
% Define regions of interest, then visualize with imageROIvalues
% ROIs group channels together for region-level analysis
figure;
fnir1200 = pf2.import.sampleData.fNIR1200();
processed = processFNIRS2(fnir1200);

% Define 3 ROIs: left, center, right prefrontal
processed = pf2.probe.roi.defineROI(processed, ...
    {[1,2,3,4,5], [6,7,8,9,10,11], [12,13,14,15,16]}, ...
    {'Left_PFC', 'Center_PFC', 'Right_PFC'});

% Plot per-ROI mean values as a 2D heatmap
roiMeans = [0.8, 1.5, 0.3];
pf2.probe.plot.imageROIvalues(processed, roiMeans, ...
    'title', 'ROI Mean HbO', 'colorbarTitle', '\muM');

%% 38. ROI interpolated 2D topography
% Smooth interpolated 2D surface from ROI values
figure;
pf2.probe.plot.interpolateROIvalues(roiMeans, processed);

%% ========================================================================
%  MISC FEATURES
%  ========================================================================

%% 39. Optode direction lines
% Show direction vectors emanating from each source/detector
figure;
pf2.probe.plot.interpolateValues3D(1:16, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'optodeLines', true, ...
    'SDLabels', true, 'ChannelLabels', false);

%% 40. Projected optode locations
% Projects optode positions onto the brain surface mesh
figure;
pf2.probe.plot.interpolateValues3D(1:16, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'useProjectedOptodeLocations', true, ...
    'SDLabels', true, 'ChannelLabels', true);

%% 41. showProbe3D — quick probe visualization
% Convenience wrapper: shows probe geometry without data coloring
figure;
fnir1200 = pf2.import.sampleData.fNIR1200();
processed = processFNIRS2(fnir1200);
pf2.probe.plot.showProbe3D(processed, 'initCamPosition', 'front');

%% 42. High-res vs low-res mesh comparison
fnirData = [1, 2, 1, 2, 2, 3, 3, 3, 4, 3, 4, 2, 3, 3, 1, 1];

figure;
subplot(1, 2, 1);
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'useHighRes', true, ...
    'ChannelLabels', false, 'SDLabels', false, ...
    'showColorbar', false);
title('High-res mesh');

subplot(1, 2, 2);
pf2.probe.plot.interpolateValues3D(fnirData, {'fNIR_Devices_fNIR1000.cfg'}, ...
    'useHighRes', false, ...
    'ChannelLabels', false, 'SDLabels', false, ...
    'showColorbar', false);
title('Low-res mesh');

%% ========================================================================
%  COORDINATE TRANSFORMATION (CapTrak/Subject-Space to MNI)
%  ========================================================================
% SNIRF files may contain probe coordinates in subject-specific coordinate
% systems (e.g., CapTrak based on LPA/RPA/Nasion landmarks) rather than
% standard MNI space. When landmarks (10-20 electrode positions) are
% available in the SNIRF file, interpolateValues3D can automatically
% transform coordinates to MNI for proper alignment with the brain template.
%
% The 'transformToMNI' parameter controls this behavior:
%   'auto' (default) - Transform if coordinate system is not MNI and
%                      landmarks are available in the data
%   true             - Always transform (error if no landmarks available)
%   false            - Never transform, use original coordinates as-is

%% 43. Automatic MNI transformation — SNIRF with CapTrak coordinates
% When importing SNIRF files with subject-specific coordinates (e.g., CapTrak),
% the transformation happens automatically if landmarks are present.
% This example uses a SNIRF file from the FRESH dataset (OpenNeuro).

% Note: Update path to your SNIRF file location
snirfPath = fullfile(pf2_base.pf2_defaultRootPath(), ...
    'tests', 'benchmarks', 'data', 'dataset_I_auditory', ...
    'sub-01', 'ses-01', 'nirs', 'sub-01_ses-01_task-FRESHAUDIO_nirs.snirf');

if isfile(snirfPath)
    % Import SNIRF with landmarks (CapTrak coordinate system)
    data = pf2.import.importSNIRF(snirfPath, false);
    processed = processFNIRS2(data);

    % Check what coordinate system was detected
    fprintf('Coordinate system: %s\n', data.device.CoordinateSystem);
    fprintf('Landmarks available: %d\n', height(data.device.Landmarks));

    % Get HbO at a single timepoint
    t_idx = round(size(processed.HbO, 1) / 2);
    hboVals = processed.HbO(t_idx, :);

    % Plot with auto transformation (default) — coordinates transformed to MNI
    figure;
    pf2.probe.plot.interpolateValues3D(hboVals, processed, ...
        'initCamPosition', 'front', ...
        'colorbarStr', 'HbO (\muM\cdotmm)', ...
        'ChannelLabels', true, 'SDLabels', false);
    title(sprintf('Auto MNI Transform (t=%.1fs)', processed.time(t_idx)));
else
    fprintf('SNIRF file not found at: %s\n', snirfPath);
    fprintf('Download FRESH dataset from OpenNeuro or update the path.\n');
end

%% 44. Compare transformed vs original coordinates
% Side-by-side comparison showing the effect of MNI transformation.
% Without transformation, CapTrak coordinates may not align with the brain.

if isfile(snirfPath)
    figure;

    % With transformation (properly aligned)
    subplot(1, 2, 1);
    pf2.probe.plot.interpolateValues3D(hboVals, processed, ...
        'transformToMNI', true, ...
        'initCamPosition', 'front', ...
        'ChannelLabels', false, 'SDLabels', false, ...
        'showColorbar', false);
    title('transformToMNI = true');

    % Without transformation (may be misaligned)
    subplot(1, 2, 2);
    pf2.probe.plot.interpolateValues3D(hboVals, processed, ...
        'transformToMNI', false, ...
        'initCamPosition', 'front', ...
        'ChannelLabels', false, 'SDLabels', false, ...
        'showColorbar', false);
    title('transformToMNI = false (original coords)');
end

%% 45. Direct use of pf2.probe.transformToMNI
% You can also use the transformation function directly to transform
% arbitrary coordinates using the landmarks from an fNIRS struct.

if isfile(snirfPath)
    % Get original probe coordinates
    origCoords = [processed.device.Probe{1}.OptPos.x, ...
                  processed.device.Probe{1}.OptPos.y, ...
                  processed.device.Probe{1}.OptPos.z];

    % Transform to MNI space
    [mniCoords, T, matchInfo] = pf2.probe.transformToMNI(origCoords, processed, ...
        'Verbose', true);

    % Display transformation results
    fprintf('\nTransformation summary:\n');
    fprintf('  Matched landmarks: %d\n', matchInfo.nMatched);
    fprintf('  Scale factor: %.4f\n', T.s);
    fprintf('  RMS error: %.2f mm\n', T.error);

    % Compare coordinate ranges
    fprintf('\nOriginal coordinate ranges (mm):\n');
    fprintf('  X: [%.1f, %.1f]\n', min(origCoords(:,1)), max(origCoords(:,1)));
    fprintf('  Y: [%.1f, %.1f]\n', min(origCoords(:,2)), max(origCoords(:,2)));
    fprintf('  Z: [%.1f, %.1f]\n', min(origCoords(:,3)), max(origCoords(:,3)));

    fprintf('\nMNI coordinate ranges (mm):\n');
    fprintf('  X: [%.1f, %.1f]\n', min(mniCoords(:,1)), max(mniCoords(:,1)));
    fprintf('  Y: [%.1f, %.1f]\n', min(mniCoords(:,2)), max(mniCoords(:,2)));
    fprintf('  Z: [%.1f, %.1f]\n', min(mniCoords(:,3)), max(mniCoords(:,3)));
end