BIDS import: Added support for ACPC and CapTrak coordinate systems

Author: ftadel

doc/updates.txt

@@ -1,5 +1,6 @@
 November 2022
 - Forward: Display leadfield sensitivity (surface, MRI, isosurface)
+- BIDS import: Added support for ACPC and CapTrak coordinate systems
 --------------------------------------------------------------
 October 2022
 - EEG: Project channel files between subjects

toolbox/anatomy/cs_compute.m

@@ -5,7 +5,7 @@
 %
 % INPUT:
 %     - sMri   : Brainstorm MRI structure
-%     - csname : Coordinate system for which we need to evaluate the transformation {'scs','mni','acpc','tal'}
+%     - csname : Coordinate system for which we need to evaluate the transformation {'scs','mni','acpc','tal','captrak'}
 
 % @=============================================================================
 % This function is part of the Brainstorm software:
@@ -71,7 +71,6 @@
             sMri.SCS.T      = Transf.T;
         end
 
-        
     % ===== MRI => MNI =====
     case 'mni'
         error('To estimate the MNI coordinates: right-click on the MRI > MNI normalization.');
@@ -80,7 +79,7 @@
     case 'acpc'
         % The necessary points are not defined
         if isempty(sMri) || ~isfield(sMri, 'NCS') || ~isfield(sMri.NCS, 'AC') || ~isfield(sMri.NCS, 'PC') || ~isfield(sMri.NCS, 'IH') || (length(sMri.NCS.AC)~=3) || (length(sMri.NCS.PC)~=3) || (length(sMri.NCS.IH)~=3)
-            disp('BST> Cannot compute MRI=>TAL transformation: Missing fiducial points.');
+            disp('BST> Cannot compute MRI=>ACPC transformation: Missing fiducial points.');
             return;
         end
         % Get coordinates in meters
@@ -108,6 +107,38 @@
         Transf.R = transform(1:3,1:3);
         Transf.T = transform(1:3,4);
 
+    % ===== MRI => CapTrak =====
+    case 'captrak'
+        % The necessary points are not defined
+        if isempty(sMri) || ~isfield(sMri, 'SCS') || ~isfield(sMri.SCS, 'NAS') || ~isfield(sMri.SCS, 'LPA') || ~isfield(sMri.SCS, 'RPA') || (length(sMri.SCS.NAS)~=3) || (length(sMri.SCS.LPA)~=3) || (length(sMri.SCS.RPA)~=3)
+            disp('BST> Cannot compute MRI=>SCS transformation: Missing fiducial points.');
+            return;
+        end
+        % Get coordinates in meters
+        NAS = double(sMri.SCS.NAS(:))' ./ 1000;
+        LPA = double(sMri.SCS.LPA(:))' ./ 1000;
+        RPA = double(sMri.SCS.RPA(:))' ./ 1000;
+        % X axis: From LPA through RPA exactly
+        dirx   = RPA - LPA;
+        dirx   = dirx/norm(dirx);
+        % Y axis: Orthogonal to the X-axis through the nasion (NAS)
+        origin = LPA + dirx * sum((NAS - LPA) .* dirx);
+        diry = NAS - origin;
+        diry   = diry/norm(diry);
+        % Z axis: Orthogonal to the XY-plane through the vertex of the head
+        dirz = cross(dirx,diry);
+        % Compute the rotation matrix
+        rot = eye(4);
+        rot(1:3,1:3) = inv(eye(3) / [dirx; diry; dirz]);
+        % compute the translation matrix
+        tra = eye(4);
+        tra(1:4,4) = [-origin(:); 1];
+        % Combine these to compute the full homogeneous transformation matrix
+        transform = rot * tra;
+        % Return in split format
+        Transf.R = transform(1:3,1:3);
+        Transf.T = transform(1:3,4);
+
     % ===== SCS => TAL =====
     % Not a real TALAIRACH system, used only from tess_envelope.m
     case 'tal'

toolbox/anatomy/cs_convert.m

@@ -6,8 +6,8 @@
 %
 % INPUT: 
 %     - sMri         : Brainstorm MRI structure
-%     - src          : Current coordinate system {'voxel','mri','scs','mni','acpc','world'}
-%     - dest         : Target coordinate system {'voxel','mri','scs','mni','acpc','world'}
+%     - src          : Current coordinate system {'voxel','mri','scs','mni','acpc','world','captrak'}
+%     - dest         : Target coordinate system {'voxel','mri','scs','mni','acpc','world','captrak'}
 %     - P            : a Nx3 matrix of point coordinates to convert
 %     - isNanAllowed : If 0, throw an error and if the conversion of some points P leads to NaN,
 %                      and if there is no linear MNI registration to be used instead
@@ -29,6 +29,10 @@
 %               Axis Y: Negative y-axis is passing from AC through PC
 %               Axis Z: Passing through a mid-hemispheric point in the superior direction
 %     - world : Transformation available in the initial file loaded as the default MRI (vox2ras/qform/world transformation)
+%     - captrak: RAS orientation and the origin approximately between LPA and RPA
+%               Axis X: From LPA through RPA exactly
+%               Axis Y: Orthogonal to the X-axis through the nasion (NAS)
+%               Axis Z: Orthogonal to the XY-plane through the vertex of the head
 
 % @=============================================================================
 % This function is part of the Brainstorm software:
@@ -108,6 +112,16 @@
     scs2acpc = [tACPC.R, tACPC.T; 0 0 0 1];
 end
 
+% ===== COMPUTE CAPTRAK TRANSFORMATION =====
+if strcmpi(src, 'captrak') || strcmpi(dest, 'captrak')
+    tCapTrak = cs_compute(sMri, 'captrak');
+    if isempty(tCapTrak) || isempty(tCapTrak.R)
+        P = [];
+        return;
+    end
+    scs2captrak = [tCapTrak.R, tCapTrak.T; 0 0 0 1];
+end
+
 % ===== CONVERT SRC => MRI =====
 % Evaluate the transformation to apply
 switch lower(src)
@@ -151,6 +165,9 @@
     case 'acpc'
         % ACPC => SCS => MRI
         RT1 = inv(scs2acpc);
+    case 'captrak'
+        % CapTrak => SCS => MRI
+        RT1 = inv(scs2captrak);
     case 'world'
         RT1 = world2mri;
     otherwise
@@ -196,6 +213,9 @@
     case 'acpc'
         % MRI => SCS => ACPC
         RT2 = scs2acpc;
+    case 'captrak'
+        % MRI => SCS => CapTrak
+        RT2 = scs2captrak;
     case 'world'
         RT2 = mri2world;
     otherwise

toolbox/io/import_channel.m

@@ -136,13 +136,13 @@
         FileUnits = 'm';
 
     % ===== EEG ONLY =====
-    case {'BIDS-SCANRAS-MM', 'BIDS-MNI-MM', 'BIDS-ACPC-MM', 'BIDS-ALS-MM'}
+    case {'BIDS-SCANRAS-MM', 'BIDS-MNI-MM', 'BIDS-ACPC-MM', 'BIDS-ALS-MM', 'BIDS-CAPTRAK-MM'}
         ChannelMat = in_channel_bids(ChannelFile, 0.001);
         FileUnits = 'mm';
-    case {'BIDS-SCANRAS-CM', 'BIDS-MNI-CM', 'BIDS-ACPC-CM', 'BIDS-ALS-CM'}
+    case {'BIDS-SCANRAS-CM', 'BIDS-MNI-CM', 'BIDS-ACPC-CM', 'BIDS-ALS-CM', 'BIDS-CAPTRAK-CM'}
         ChannelMat = in_channel_bids(ChannelFile, 0.01);
         FileUnits = 'cm';
-    case {'BIDS-SCANRAS-M', 'BIDS-MNI-M', 'BIDS-ACPC-M', 'BIDS-ALS-M'}
+    case {'BIDS-SCANRAS-M', 'BIDS-MNI-M', 'BIDS-ACPC-M', 'BIDS-ALS-M', 'BIDS-CAPTRAK-M'}
         ChannelMat = in_channel_bids(ChannelFile, 1);
         FileUnits = 'm';
 
@@ -410,7 +410,7 @@
         warning(['WARNING: When importing ACPC positions for multiple subjects: the ACPC transformation from the first subject is used for all of them.' 10 ...
                  'Please consider importing your subjects seprately.']);
     end
-    % If we know the destination study: convert from MNI to SCS coordinates
+    % If we know the destination study: convert from ACPC to SCS coordinates
     if ~isempty(iStudies)
         % Get the subject for the first study
         sStudy = bst_get('Study', iStudies(1));
@@ -433,6 +433,36 @@
     % Do not convert the positions to SCS
     isAlignScs = 0;
 
+%% ===== CAPTRAK TRANSFORMATION =====
+elseif ismember(FileFormat, {'BIDS-CAPTRAK-MM', 'BIDS-CAPTRAK-CM', 'BIDS-CAPTRAK-M'})
+    % Warning for multiple studies
+    if (length(iStudies) > 1)
+        warning(['WARNING: When importing CapTrak positions for multiple subjects: the CapTrak transformation from the first subject is used for all of them.' 10 ...
+                 'Please consider importing your subjects seprately.']);
+    end
+    % If we know the destination study: convert from CapTrak to SCS coordinates
+    if ~isempty(iStudies)
+        % Get the subject for the first study
+        sStudy = bst_get('Study', iStudies(1));
+        sSubject = bst_get('Subject', sStudy.BrainStormSubject);
+        % Get the subject's MRI
+        if isempty(sSubject.Anatomy) || isempty(sSubject.Anatomy(1).FileName)
+            error('You need the subject anatomy in order to load sensor positions in CapTrak coordinates.');
+        end
+        % Load the MRI
+        MriFile = file_fullpath(sSubject.Anatomy(1).FileName);
+        sMri = in_mri_bst(MriFile);
+        if ~isfield(sMri, 'SCS') || ~isfield(sMri.SCS, 'R') || isempty(sMri.SCS.R) || ~isfield(sMri.SCS, 'NAS') || isempty(sMri.SCS.NAS) || ~isfield(sMri.SCS, 'LPA') || isempty(sMri.SCS.LPA) || ~isfield(sMri.SCS, 'RPA') || isempty(sMri.SCS.RPA)
+            error(['All fiducials must be defined for this subject (NAS,LPA,RPA)' 10 'in order to load sensor positions in CapTrak coordinates.']);
+        end
+        % Convert all the coordinates: CapTrak => SCS
+        fcnTransf = @(Loc)cs_convert(sMri, 'captrak', 'scs', Loc')';
+        AllChannelMats = channel_apply_transf(ChannelMat, fcnTransf, [], 1);
+        ChannelMat = AllChannelMats{1};
+    end
+    % Do not convert the positions to SCS
+    isAlignScs = 0;
+
 %% ===== MRI/NII TRANSFORMATION =====
 % If the SCS coordinates are not defined (NAS/LPA/RPA fiducials), try to use the MRI=>subject transformation available in the MRI (eg. NIfTI sform/qform)
 % Only available if there is one study in output

toolbox/process/functions/process_import_bids.m

@@ -20,12 +20,15 @@
 %    - EEG  : https://openneuro.org/datasets/ds004024 : OK
 %    - iEEG : https://openneuro.org/datasets/ds003688 : ERROR: Wrong interpretation of ACPC coordinates (easier to see in ECOG for sub-02)
 %    - iEEG : https://openneuro.org/datasets/ds003848 : WARNING: Impossible to interepret correctly the coordinates in electrodes.tsv
-%    - iEEG : https://openneuro.org/datasets/ds004085 : 
+%    - iEEG : https://openneuro.org/datasets/ds004085 : ERROR: FreeSurfer not imported (not in a "freesurfer" pipeline subfolder)
+%                                                       ERROR: Some subjects (1,2,10,11) have incorrect channel names in the electrodes.tsv => Not imported
+%                                                       ERROR: Inaccurate electrode positioning because AC/PC landmarks are not defined in the MRI
 %    - iEEG : https://openneuro.org/datasets/ds004126 : OK (ACPC OK)
 %    - STIM : https://openneuro.org/datasets/ds002799 : WARNING: No channel file imported because there are no SEEG recordings
 %    - MEG  : https://openneuro.org/datasets/ds000117 : 
+%    - MEG  : https://openneuro.org/datasets/ds000246 : 
 %    - MEG  : https://openneuro.org/datasets/ds000247 : 
-%    - MEG  : https://openneuro.org/datasets/ds004107 : 
+%    - MEG  : https://openneuro.org/datasets/ds004107 : WARNING: Multiple NAS/LPA/RPA in T1w.json, one for each MEG session => Used the average for both sessions
 %    - Tutorial FEM  : https://neuroimage.usc.edu/brainstorm/Tutorials/FemMedianNerve   :
 %    - Tutorial ECOG : https://neuroimage.usc.edu/brainstorm/Tutorials/ECoG             :
 %    - Tutorial SEEG : https://neuroimage.usc.edu/brainstorm/Tutorials/Epileptogenicity : 
@@ -762,6 +765,8 @@
                     if ~isempty(electrodesCoordSystem)
                         if strcmpi(electrodesCoordSystem, 'ACPC')
                             electrodesSpace = 'ACPC';
+                        elseif strcmpi(electrodesCoordSystem, 'CapTrak')
+                            electrodesSpace = 'CapTrak';
                         elseif ~isempty(strfind(electrodesCoordSystem, 'MNI')) || ~isempty(strfind(electrodesCoordSystem, 'IXI')) || ~isempty(strfind(electrodesCoordSystem, 'ICBM'))  || ~isempty(strfind(electrodesCoordSystem, 'fs')) 
                             electrodesSpace = 'MNI';
                         elseif ismember(upper(electrodesCoordSystem), {'CTF', 'EEGLAB', 'EEGLAB-HJ', 'ElektaNeuromag', '4DBti', 'KitYokogawa', 'ChietiItab'})