Tags: dataformat ricoh meg

Getting started with Ricoh data

Ricoh took over the MEG branch from Yokogawa in 2016 and has started manufacturing and installing the MEG system. The Ricoh system has some updates over the Yokogawa system.

  • This page is under construction [Last updated on June 12th by KK].
  • Ricoh MEG Reader Toolbox (ver. 1.0) was merged into FieldTrip on June 11th, 2018
  • Yokogawa MEG Reader Toolbox has been upgraded from ver. 1.4 to 1.5. The reader ver. 1.5 allows you to analyze EEG data.

Introduction

The RICOH MEG is an MEG system that is developed by Ricoh Company, Ltd. Ricoh took over the MEG branch from Yokogawa in 2016 and has started manufacturing and installing the MEG system, continuing to support Yokogawa systems. They released their first MEG system at the end of 2017. A current RICOH MEG system is a successor to Yokogawa MEG systems and has most of the same characteristics as those of Yokogawa systems, especially those of the 160-channel Yokogawa MEG system.

The following data files can be read and used in FieldTrip: the files with the extensions of ‘.con’, ‘.ave’, and ‘.mrk’. All required reading low-level functions are located in the external/ricoh_meg_reader directory as a set of pre-compiled p-files and will be called by the appropriate FieldTrip functions. The low-level functions are officially supplied by Ricoh Company, Ltd. for using their data in open source software environments.

The data files, .con, .ave, and .mrk, are pre-selected standard files of the RICOH MEG system. Each contains

  • .con: continuous data,
  • .ave: data averaged by segments of interest,
  • .mrk: Head Position Indicator coils (HPI) magnetic marker coil positions,

respectively. All these files also contain information of channel types and gradiometer positions. These files can be read in FieldTrip and allow you to conduct various analyses. To fully utilize FieldTrip functions, it is however recommended that you should use yet another .con file, which can be called as an “exported .con file”. The exported .con file contains all different data objects related to your recording including simultaneously measured EEG data, timings of events, and digitized points in addition to MEG data, HPIs, and gradiometer labels and positions. This .con file can be generated by using ‘Export menu’ on the RICOH acquisition software.

In general, the first step of analyzing MEG data is composed of MEG data read-in, trial selection, and MEG-MRI co-registration. This page therefore describes how to get started reading and using these file types in FieldTrip, mainly focusing on how to use the exported .con file in such pre-processing and co-registration. It is worth noting that almost all the contents in this page can be utilized for Yokogawa data (Yokogawa’s .con, .ave, and .mrk files) because the extensions of ‘.con’, ‘.ave’, and ‘.mrk’ are common between Ricoh and Yokogawa systems and the corresponding files have the same data structure.

The functions in FieldTrip that allows you to execute the pre-processing and co-registration are: ft_read_header, ft_preprocessing, ft_read_event, ft_read_headshape, and ft_read_sens. This page will describe how these FieldTrip functions perform for Ricoh and Yokogawa data by showing several examples on pre-processing and co-registration. As for the co-registration, you need to prepare MRI file with formats supported in FieldTrip. An MRI file with NIfTI (.nii) or DICOM format is assumed to be used in this page.

Set the path

To get started, you should add the FieldTrip main directory to your path, and execute the ft_defaults function, which sets the defaults and configures up the minimal required path settings. See also this frequently asked question. You also need to set the path to your data files.

addpath <path_to_fieldtrip>
ft_defaults

%% MEG data
meg_path = <path_to_meg_data_file>;
meg_file = 'continuous_data_export.con';  % your exported con-filename
dataset = fullfile(meg_path, meg_file);

%% MRI data
mri_path = <path_to_meg_data>;
mri_file = 'mri_data.nii';  % your NIfTI filename or the first DICOM file

Read MEG data

To check if you can read in the data, you try the FieldTrip functions, ft_read_header, ft_preprocessing, ft_read_event, in the command window.

Read header

ft_read_header provides the header information embedded in an exported .con file:

hdr = ft_read_header(dataset)

>> hdr =
struct array with fields
         orig: [1x1 struct]
           Fs: 2000
       nChans: 251
     nSamples: 664000
  nSamplesPre: 0
      nTrials: 1
        label: {251x1 cell}
         grad: [1x1 struct]
     chantype: {251x1 cell}
     chanunit: {251x1 cell}

The header contains a lot of information about the measurement parameters. In this example 251 channels were recorded, and the sampling frequency was a 2000 Hz. The field ‘hdr.orig’ contains all the original header information.

The header information can be also extracted from the raw data files (the raw .con, .ave, and .mrk files). As mentioned above, the extensions of ‘.con’, ‘.ave’, and ‘.mrk,’ are common between Ricoh and Yokogawa systems. In the read-in process, the Ricoh’s extensions are distinguished from the Yokogawa’s ones through a system-version check (see the function isricohmegfile.m located in the fileio/private directory). The system version larger than and equal to 3 (i.e., ver >= 3) represents that the file is generated by Ricoh system. A file generated by Yokogawa system has a system version less than 3. You can check a version of your file by trying hdr.orig.ver.

The channel labels can be found in ‘hdr.label’. In this example, the labels are listed as follow

   'AG001' 'AG002' ... 'AG160' 'TRIG161' ... 'TRIG174' 'RM175' ... 'RM186' 'EEG187' ... 'EEG251'

‘AG*’ denotes the axial gradiometer, ‘TRIG*’ the trigger channel, ‘RM*’ the reference magnetometer, and ‘EEG*’ the EEG channel. The definitions of the gradiometer (‘AG*’) such as the positions and directions of each sensor, which are needed in source reconstruction, are contained in hdr.grad. Current Ricoh system has the same type of sensor as that of Yokogawa system, and thus you will find a name of ‘Yokogawa160’ in hdr.grad.type.

Read data

The signal time course of the axial gradiometers can be checked by using ft_preprocessing:

cfg = [];
cfg.dataset = dataset;
cfg.channel = 'AG*';
cfg.hpfilter = 'yes'; % to remove slow drift
cfg.hpfreq = 1;
cfg.hpfiltord = 5;
data = ft_preprocessing(cfg);

%% Display waveforms using ft_databrowser
cfg = [];
cfg.blocksize = 10;
cfg = ft_databrowser(cfg, data);

If you use the option cfg.channel = 'TRIG\*', you can check the waveform of the trigger signals. You can also check the data by calling a low-level function, ft_read_data. For example,

dat = ft_read_data(dataset, 'begsample', 100000, 'endsample', 120000, 'chanindx', 12);

The options such as chanindx should be specified in key-value pairs, see ft_read_data. When only the filename is specified, all time-course data of every channels with large amount of memory will be read. To avoid to load such large amount of data, you should specify channel index when you directly use the low-level read function.

Read events

In order to select pieces of data around the events in which you are interest, you can call the FieldTrip function, ft_definetrial, either using a generic definition or using your own trialfun. The trial function calls the low-level reading function, ft_read_event. The function reads event information and represents it in a common data-independent format. This event-reading function works for an exported .con file as follows:

event = ft_read_event(dataset,'threshold', 1.6)

>> event =
204x1 struct array with fields
  sample
  value
  type
  offset
  duration

ft_read_event detects all the events embedded in the exported .con file. You can deal with three types of events:

  1. User-defined event [triginfo]
  2. Annotation [annotations]
  3. Analog trigger [analogtrig (trial for Yokogawa data)]

The event types ‘triginfo’ and ‘annotations’, are available only for a Ricoh’s exported .con file. Unlike the exported .con file, the two event types are not defined in original .con and Yokogawa’s data files. Only an event type, ‘analogtrig’ (or ‘trial’ for Yokogawa data), is available for them. The option of threshold in the above example denotes a threshold value in the unit of volt (1.6 volts) and works as a threshold for analog trigger events.

Multiple events involved in the user-defined event (event.type = 'triginfo') can be categorized by event.value. The labels read in ‘event.value’ can be provided by setting your favorite names when you conduct recordings. For example, if your recording is related to so-called odd-ball task, you can name events as (event.value = ) ‘standard’ or ‘deviant’ at the trigger name boxes of the Ricoh acquisition software. You can also input a simple task name, such as ‘AEF’ and ‘SEF’, at the name boxes. The (event.type = ) ‘annotations’ denotes the events that the user put annotations on recording data. The Ricoh software allows you to put annotations on your recording data by using five predefined category labels. The categories defined are: ‘undefined’, ‘epileptic’, ‘others’, ‘noise’, and ‘text’. These annotation categories are allocated in ‘event.value’ as 0, 10, 20, 30, and 40, respectively. For example, you can easily review epileptic activities by selecting the events with event.type = 'annotations'; and event.value = 10;, identifying when such an epileptic activity happens by referring to event.sample. The (event.type = ) ‘analogtrig’ literally represents the analog triggers. The trigger channels used in your recording are listed in event.value (e.g., ‘TRIG162’).

The first 14 components of the event structure in the above example are:

table = struct2table(event);

>> table(1:14,:)
ans =
14x5 table
  sample        type           value      offset    duration
  ______    _____________    _________    ______    ________
   2380     'triginfo'       'AEF'        []        []
   2380     'analogtrig'     'TRIG162'    []        []
   9046     'triginfo'       'AEF'        []        []
   9046     'analogtrig'     'TRIG162'    []        []
   9283     'annotations'    [     10]    []        []
  13829     'annotations'    [     10]    []        []
  15712     'triginfo'       'AEF'        []        []
  15712     'analogtrig'     'TRIG162'    []        []
  18603     'annotations'    [     40]    []        []
  22378     'triginfo'       'AEF'        []        []
  22378     'analogtrig'     'TRIG162'    []        []
  23796     'annotations'    [     30]    []        []
  29044     'triginfo'       'AEF'        []        []
  29044     'analogtrig'     'TRIG162'    []        []

Trial Selection

In the same way as in the Preprocessing - Segmenting and reading trial-based EEG and MEG data tutorial, you can define segments of epochs-of-interest (trials) in your recorded MEG data using the FieldTrip function, ft_definetrial. For example,

%% Define trials
cfg = [];
cfg.dataset             = dataset;
cfg.trialfun            = 'ft_trialfun_general';
cfg.trialdef.eventtype  = 'triginfo';
cfg.trialdef.eventvalue = 'AEF';
cfg.trialdef.prestim    = 0.1; % sec
cfg.trialdef.poststim   = 0.5; % sec
cfg = ft_definetrial(cfg);
trl = cfg.trl;

%% Segment data into pieces
cfg = [];
cfg.trl = trl;
data = ft_redefinetrial(cfg, data);

If you change the options for eventtype and eventvalue as

cfg.trialdef.eventtype  = 'annotations';
cfg.trialdef.eventvalue = [30 40];

you can select the events that were annotated as ‘noise’ and ‘text’. Alternatively, if you use the options

cfg.trialdef.eventtype  = 'analogtrig';
cfg.trialdef.eventvalue = 'TRIG162';

you can select the trials based on the triggers signals of the channel ‘TRIG162’. This trigger-based trial selection can also be achieved by using your own trial function

%% Define trials
cfg                      = [];
cfg.dataset              = dataset;
cfg.trialfun             = 'yourowntrialfun';
cfg.trialdef.eventtype   = 'analogtrig'; % use 'trial' for Yokogawa data
cfg.trialdef.eventvalue  = 'TRIG162';
cfg.trialdef.trigchannel = cfg.trialdef.eventvalue;
cfg.trialdef.threshold   = 1.6; %default
cfg.trialdef.detectflank = 'up'; %default
cfg.trialdef.prestim     = 0.1; %sec
cfg.trialdef.poststim    = 0.5; %sec
cfg = ft_definetrial(cfg);
trl = cfg.trl;

An example of your own trial function, yourowntrialfun is

function [trl, event] = yourowntrialfun(cfg)

% read the header information and the events from the data
hdr = ft_read_header(cfg.dataset);
nSamples = hdr.nSamples * hdr.nTrials;
event = ft_read_event(cfg.dataset, 'threshold', cfg.trialdef.threshold, 'detectflank', cfg.trialdef.detectflank);

% trigger events detected by the designated 'eventvalue' (i.e., trigger channel)
sample = [event(find(strcmp(cfg.trialdef.eventvalue, {event.value}))).sample]';

% For yokogawa data, use
% sample = [event(find(strcmp(cfg.trialdef.eventvalue, {event.type}))).sample]';

% creating your own definition of trials
trl = [];
pretrig = cfg.trialdef.prestim*hdr.Fs;
posttrig = cfg.trialdef.poststim*hdr.Fs;
for j = 1:length(sample);
  trlbegin = sample(j) - pretrig;
  trlend   = sample(j) + posttrig - 1;
  offset   = -pretrig;
  if ( 0 < trlbegin ) && ( trlend < nSamples )
    newtrl   = [trlbegin trlend offset];
    trl      = [trl; newtrl];
  end
end

MRI-MEG Co-registration

The registration between MRI and MEG is essential for source-space analysis on MEG data. The goal of the co-registration is to transform the positions of the MRI voxel and the sensor array into a common head coordinate. This part describes how to co-register MRI and MEG data recorded by Ricoh system, showing its examples. Although FieldTrip supports various anatomical MRI data formats as presented here, anatomical MRI data are assumed to be saved as NIfTI (.nii) or DICOM files in this page.

The co-registration is usually done in two steps; the first step is an approximate alignment that is based on fiducial points, and the second step is a refinement with additional head-shape digitized points. You practically have two possible cases regarding to the first step:

  1. To use HPIs (marker coils) as the fiducial points
  2. To use three anatomical landmarks (nasion, left, and right pre-auricular points) as the fiducial points

For the case 1, the positions corresponding to the HPIs are needed to be displayed on the MRI (by vitamin-drug dot images in most cases). For the case 2, it is necessary that the three landmarks, NAS, LPA, and RPA, were picked by a 3D scanner (digitizer) and already have been coregistered with HPIs (i.e., MEG). For Ricoh and Yokogawa data (.con, .ave, .mrk), the HPI positions, the anatomical landmarks, and the additional digitized points can be elicited by the FieldTrip function, ‘ft_read_headshape’.

In the Ricoh and Yokogawa system, the coordinate system is defined as:

  • Origin is at the center of the helmet
  • X-axis goes towards the nose
  • Y-axis goes towards the left
  • Z-axis goes towards the top of the head

This coordinate system is same as the CTF coordinate system (see the faq). It is therefore convenient to introduce the CTF coordinate system for defining a head coordinate system. In the co-registration procedure described below, there is a step of defining a head coordinate system with the CTF convention.

Read MRI

Load MRI and rescale it:

%% Load MRI
mri = ft_read_mri(fullfile(mri_path, mri_file));
mri = ft_convert_units(mri, 'mm');

%% Rescale MRI
cfg = [];
mri = ft_volumereslice(cfg, mri)

If your MRI has HPI positions, e.g., represented by vitamin E capsules, load their positions [otherwise skip the code below)

%% Load HPIs from .txt file (vitamin E capsule positions in MRI)
hpi_mri = load(fullfile(fid_path, fid_file)) % The unit is needed to be in mm

%% Derive the voxel locations
hpi_vox = ft_warp_apply(inv(mri.transform), hpi_mri);

The “fid_hpi.txt” file contains vitamin-E dot positions (with the units of mm) in the MRI with the order of ‘HPI_3’, ‘HPI_1’, ‘HPI_2’, ‘HPI_4’, and ‘HPI_5’. Including the above five HPIs, several special points are defined in the Ricoh system (this is the same as Yokogawa system.); the eight special points are defined in the Ricoh system in the convention with the following order:

  1. fidt9: Left pre-auricular point [lpa]
  2. HPI_1: Left pre-auricular magnetic coil position [LPA]
  3. HPI_4: Forehead left magnetic coil position [LPF]
  4. HPI_3: Forehead center magnetic coil position [CPF]
  5. HPI_5: Forehead right magnetic coil position [RPF]
  6. HPI_2: Right pre-auricular magnetic coil position [RPA]
  7. fidt10: Right pre-auricular point [rpa]
  8. fidnz: Nasion [nas]

In accordance with this convention, the user of the Ricoh system is required to point these HPIs in this sequence when using a 3D scanner (digitizer). After selecting these 8 positions, the user should digitize the positions of EEG electrodes and head surface points.

Define a head coordinate system for MRI

Define a head coordinate system by manually taking up the anatomical landmarks, NAS, LPA, and RPA. A GUI window will pop up. The way to use the GUI to specify the fiducial points will be instructed in your command window.

%% Define manually the fiducials and define a head coordinate system
cfg = [];
cfg.method = 'interactive';
cfg.coordsys = 'ctf';
mri = ft_volumerealign(cfg, mri);

%% Fiducials on MRI
fid_pos = mri.cfg.fiducial;
fid_vox = [fid_pos.nas; fid_pos.lpa; fid_pos.rpa];

%% Represent fiducial points and HPIs in the head coordinate system for MRI
fid_mri = ft_warp_apply(mri.transform, fid_vox, 'homogeneous');

Read HPIs and anatomical landmarks that are represented in MEG coordinate

Read the HPIs and anatomical landmarks by calling ft_read_headshap

%% Read HPIs & anatomical landmarks in MEG coordinate system from .con, .ave, or .mrk files
headshape = ft_read_headshape(fullfile(meg_path, meg_file));
headshape = ft_convert_units(headshape, 'mm');

%% HPIs (marker coils)
coil_pos = headshape.fid.pos(end-4:end,:);
coil_label = headshape.fid.label(end-4:end,:);

%% Anatomical landmarks
fid_pos = headshape.fid.pos(1:3,:);
fid_label = headshape.fid.label(1:3,:);

For an exported .con file that contains the information of digitized points, headshape provides the full of the digitized-point information. Furthermore, headshape.fid provides the list of the anatomical fiducial points and HPIs measured by using the digitizer. An example is as follows:

headshape = ft_read_headshape(fullfile(meg_path, meg_file));

>> headshape =
struct with field
    pos: [106x3 double]
  label: {106x1 cell}
    fid: [1x1 struct]
   unit: 'cm'

>> headshape.fid.label
ans =
13x1 cell-array
  'nas'
  'lpa'
  'rpa'
  'HPI_3'
  'HPI_1'
  'HPI_2'
  'HPI_4'
  'HPI_5'
  'CPF'
  'LPA'
  'RPA'
  'LPF'
  'RPF'

The first three components of headshape.fid denote the digitized anatomical landmarks, NAS, LPA, and RPA. The middle five components, ‘HPI_3’, ‘HPI_1’, ‘HPI_2’, ‘HPI_4’, and ‘HPI_5’, denote the “digitized” HPI points that corresponds to the HPIs (marker coils), ‘CPF’, ‘LPA’, ‘RPA’, ‘LPF’, and ‘RPF’, respectively. The values of the positions of these digitized points are given by those in the MEG (i.e., dewar) coordinate. The last five components of ‘headshape.fid’ are secured for reading the HPI (marker coils) information and accordingly provides the information of the HPI positions with the label of ‘CPF’, ‘LPA’, ‘RPA’, ‘LPF’, and ‘RPF’.

On contrary to an exported .con file, an original (not-exported) data file (.con, .ave, .mrk) contains only HPI information without digitized points. The headshape.fid of an original data file provides only marker-coil information

headshape = ft_read_headshape(fullfile(meg_path, meg_file));

>> headshape =
struct with field
   pos: []
   fid: [1x1 struct]
  unit: 'cm'

>> headshape.fid.label
ans =
5x1 cell-array
  'nas'
  'lpa'
  'rpa'
  'Marker4'
  'Marker5'

The labels ‘nas’, ‘lpa’, ‘rpa’, ‘Marker4’, and ‘Marker5’, represent the magnetic-marker coils at ‘CPF’, ‘LPA’, ‘RPA’, ‘LPF’, and ‘RPF’ in order. The name of ‘nas’ is just a label, not indicating the anatomical nasion. The values of headshape.fid.pos represent the marker coil positions estimated through MEG measurement. It should be noted that the case 2 as the first step of co-registration is not applicable to an original data file, and only the case 1 is applicable.

Derive the transformation matrix from the MEG to Head coordinate systems

Now you conduct the alignment to derive the transformation matrix from the MEG to the defined head coordinate systems. There are two routes in accordance with the cases 1 and 2.

Case 1: Using HPI (marker coil) position

%% Get the transformation matrix from MEG to Head
[TR, TT, ER, t, info] = icp(hpi_mri',coil_pos');
meg2ctf = [[TR TT]; 0 0 0 1];

Case 2: Using anatomical landmarks

%% Define a head coordinate for MEG
meg2head = ft_headcoordinates(fid_pos(1,:), fid_pos(2,:), fid_pos(3,:), 'ctf');
fid_pos = ft_warp_apply(meg2head, fid_pos, 'homogeneous');

%% Derive the transformation matrix between head coordinates for MRI and MRI
[TR, TT, ER, t, info] = icp(fid_mri',fid_pos');
head2head = [[TR TT]; 0 0 0 1];

%% Get the transformation matrix from MEG to Head
meg2ctf = meg2head*head2head;

For both cases, the iterative closest point (ICP) algorithm, in which [TR, TT] = icp(q,p) returns the rotation matrix TR and translation vector TT that minimizes the distances from (TR * p + TT) to q, can be employed to derive the transformation matrix.

Co-register the sensor array with MRI

This is the last part of the first step. Load the gradiometer positions and transform them using the transformation matrix ‘meg2ctf’

% Read gradiometer definition in MEG coordinate system
grad = ft_read_sens(fullfile(meg_path, meg_file), 'senstype', 'meg');
grad = ft_convert_units(grad, 'mm');

% Co-registration
mri_coreg = mri;
realignedgrad = ft_transform_geometry(meg2ctf, grad);

% Save data
save(fullfile(output_path, 'realignedgrad'), '-struct', 'realignedgrad');

% Represent realigned HPIs & anatomical landmarks using the head coordinat
headshape = ft_transform_geometry(meg2ctf, headshape);
coil_pos = headshape.fid.pos(end-4:end,:);
fid_pos = headshape.fid.pos(1:3,:);

Now the MEG, MRI, and digitizer are all living in the common head coordinate system.

Realign the scalp surface with the digitized cloud [2nd Step]

To refine the above fiducial-points based registration, it is recommended to utilize the full of digitized head points by matching them with the scalp surface that is extracted from MRI. This procedure is applicable only to an exported .con file. This refinement procedure can be done by employing the FieldTrip function, ft_volumerealign, with the option cfg.headshape.icp = 'yes'.

%% volumerealign with IC
cfg = [];
cfg.method = 'headshape';
cfg.headshape.interactive = 'yes';
cfg.headshape.icp = 'yes';
cfg.headshape.headshape = headshape;
cfg.coordsys = 'ctf';
mri_coreg = ft_volumerealign(cfg, mri_coreg);

%% For representing anatomical landmarks and HPIs in the head coordinate system
fid_mri = ft_warp_apply(mri_coreg.transform, fid_vox, 'homogeneous');
hpi_mri = ft_warp_apply(mri_coreg.transform, hpi_vox, 'homogeneous');

%% Save data
save(fullfile(output_path, 'mri_coreg'), '-struct', 'mri_coreg')

Check the result

%% Plot the head model and sensor array
figure
ft_plot_sens(realignedgrad, 'style','*b')
hold on
ft_plot_mesh(mesh_scalp, 'edgecolor','none','facealpha',0.6,'facecolor',[0.6 0.6 0.8]);
ft_plot_headmodel(headmodel,'edgecolor','b')
ft_plot_headshape(headshape)
ft_plot_axes([], 'unit', 'mm');
% plot the fiducial points
plot3(coil_pos(1,1), coil_pos(1,2), coil_pos(1,3), 'r.', 'MarkerSize', 25);
plot3(coil_pos(2,1), coil_pos(2,2), coil_pos(2,3), 'r.', 'MarkerSize', 25);
plot3(coil_pos(3,1), coil_pos(3,2), coil_pos(3,3), 'r.', 'MarkerSize', 25);
plot3(coil_pos(4,1), coil_pos(4,2), coil_pos(4,3), 'r.', 'MarkerSize', 25);
plot3(coil_pos(5,1), coil_pos(5,2), coil_pos(5,3), 'r.', 'MarkerSize', 25);
plot3(fid_pos(1,1), fid_pos(1,2), fid_pos(1,3), 'k.', 'MarkerSize', 25);
plot3(fid_pos(2,1), fid_pos(2,2), fid_pos(2,3), 'k.', 'MarkerSize', 25);
plot3(fid_pos(3,1), fid_pos(3,2), fid_pos(3,3), 'k.', 'MarkerSize', 25);
plot3(fid_mri(1,1), fid_mri(1,2), fid_mri(1,3), 'y.', 'MarkerSize', 30);
plot3(fid_mri(2,1), fid_mri(2,2), fid_mri(2,3), 'y.', 'MarkerSize', 30);
plot3(fid_mri(3,1), fid_mri(3,2), fid_mri(3,3), 'y.', 'MarkerSize', 30);
plot3(hpi_mri(1,1), hpi_mri(1,2), hpi_mri(1,3), 'g.', 'MarkerSize', 25);
plot3(hpi_mri(2,1), hpi_mri(2,2), hpi_mri(2,3), 'g.', 'MarkerSize', 25);
plot3(hpi_mri(3,1), hpi_mri(3,2), hpi_mri(3,3), 'g.', 'MarkerSize', 25);
plot3(hpi_mri(4,1), hpi_mri(4,2), hpi_mri(4,3), 'g.', 'MarkerSize', 25);
plot3(hpi_mri(5,1), hpi_mri(5,2), hpi_mri(5,3), 'g.', 'MarkerSize', 25);