Cropped_image_mne_bids.py

##Modified Version  from https://mne.tools/mne-bids/stable/auto_examples/convert_mri_and_trans.html#sphx-glr-auto-examples-convert-mri-and-trans-py
##mne-bids version 0.7

# Authors: Stefan Appelhoff <stefan.appelhoff@mailbox.org>
#          Alex Rockhill <aprockhill206@gmail.com>
#          Alex Gramfort <alexandre.gramfort@inria.fr>
# License: BSD (3-clause)

###############################################################################
# We are importing everything we need for this example:

import os.path as op
import shutil

import numpy as np
import matplotlib.pyplot as plt

import nibabel as nib
from nilearn.plotting import plot_anat

import mne
from mne.datasets import sample
from mne.source_space import head_to_mri

from mne_bids import (write_raw_bids, BIDSPath, write_anat,
                      get_head_mri_trans, print_dir_tree)

###############################################################################
# We will be using the `MNE sample data <mne_sample_data_>`_ and write a basic
# BIDS dataset. For more information, you can checkout the respective
# :ref:`example <ex-convert-mne-sample>`.

data_path = sample.data_path()
event_id = {'Auditory/Left': 1, 'Auditory/Right': 2, 'Visual/Left': 3,
            'Visual/Right': 4, 'Smiley': 5, 'Button': 32}
raw_fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw.fif')
events_data = op.join(data_path, 'MEG', 'sample', 'sample_audvis_raw-eve.fif')
output_path = op.abspath(op.join(data_path, '..', 'MNE-sample-data-bids'))

###############################################################################
# To ensure the output path doesn't contain any leftover files from previous
# tests and example runs, we simply delete it.
#
# .. warning:: Do not delete directories that may contain important data!
#

if op.exists(output_path):
    shutil.rmtree(output_path)

###############################################################################
# Read the input data and store it as BIDS data.

raw = mne.io.read_raw_fif(raw_fname)
raw.info['line_freq'] = 60  # specify power line frequency as required by BIDS

sub = '01'
ses = '01'
task = 'audiovisual'
run = '01'
bids_path = BIDSPath(subject=sub, session=ses, task=task,
                     run=run, root=output_path)
write_raw_bids(raw, bids_path, events_data=events_data,
               event_id=event_id, overwrite=True)

###############################################################################
# Print the directory tree
print_dir_tree(output_path)

###############################################################################
# Writing T1 image
# ----------------
#
# Now let's assume that we have also collected some T1 weighted MRI data for
# our subject. And furthermore, that we have already aligned our coordinate
# frames (using e.g., the `coregistration GUI`_) and obtained a transformation
# matrix :code:`trans`.

# Get the path to our MRI scan
t1_mgh_fname = op.join(data_path, 'subjects', 'sample', 'mri', 'T1.mgz')

# Load the transformation matrix and show what it looks like
trans_fname = op.join(data_path, 'MEG', 'sample',
                      'sample_audvis_raw-trans.fif')
trans = mne.read_trans(trans_fname)
print(trans)


#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
########## CHANGES TO CODE ################
import nibabel as nb
mri = nb.load(t1_mgh_fname)
dat = mri.get_fdata().astype(np.int16)
nifti_out = nb.Nifti1Image(dat, mri.affine)
nifti_out.to_filename('./t1w.nii')

dat2 = dat[:-40,:-80,:-20]
nifti_cropped_out = nb.Nifti1Image(dat2, mri.affine)
nifti_cropped_out.to_filename('./croppedT1.nii')
###################################



###############################################################################
# We can save the MRI to our existing BIDS directory and at the same time
# create a JSON sidecar file that contains metadata, we will later use to
# retrieve our transformation matrix :code:`trans`. The metadata will here
# consist of the coordinates of three anatomical landmarks (LPA, Nasion and
# RPA (=left and right preauricular points) expressed in voxel coordinates
# w.r.t. the T1 image.

# First create the BIDSPath object.
t1w_bids_path = \
    BIDSPath(subject=sub, session=ses, root=output_path, suffix='T1w')

# We use the write_anat function
t1w_bids_path = write_anat(
    image='./croppedT1.nii', #t1_mgh_fname,  # path to the MRI scan
    bids_path=t1w_bids_path,
    raw=raw,  # the raw MEG data file connected to the MRI
    trans=trans,  # our transformation matrix
    verbose=True, # this will print out the sidecar file
    overwrite=True
)
anat_dir = t1w_bids_path.directory

###############################################################################
# Let's have another look at our BIDS directory
print_dir_tree(output_path)

###############################################################################
# Our BIDS dataset is now ready to be shared. We can easily estimate the
# transformation matrix using ``MNE-BIDS`` and the BIDS dataset.
estim_trans = get_head_mri_trans(bids_path=bids_path)

###############################################################################
# Finally, let's use the T1 weighted MRI image and plot the anatomical
# landmarks Nasion, LPA, and RPA onto the brain image. For that, we can
# extract the location of Nasion, LPA, and RPA from the MEG file, apply our
# transformation matrix :code:`trans`, and plot the results.

# Get Landmarks from MEG file, 0, 1, and 2 correspond to LPA, NAS, RPA
# and the 'r' key will provide us with the xyz coordinates. The coordinates
# are expressed here in MEG Head coordinate system.
pos = np.asarray((raw.info['dig'][0]['r'],
                  raw.info['dig'][1]['r'],
                  raw.info['dig'][2]['r']))

# We now use the ``head_to_mri`` function from MNE-Python to convert MEG
# coordinates to MRI scanner RAS space. For the conversion we use our
# estimated transformation matrix and the MEG coordinates extracted from the
# raw file. `subjects` and `subjects_dir` are used internally, to point to
# the T1-weighted MRI file: `t1_mgh_fname`. Coordinates are is mm.
mri_pos = head_to_mri(pos=pos,
                      subject='sample',
                      mri_head_t=estim_trans,
                      subjects_dir=op.join(data_path, 'subjects'))

# Our MRI written to BIDS, we got `anat_dir` from our `write_anat` function
t1_nii_fname = op.join(anat_dir, 'sub-01_ses-01_T1w.nii.gz')

# Plot it
fig, axs = plt.subplots(3, 1, figsize=(7, 7), facecolor='k')
for point_idx, label in enumerate(('LPA', 'NAS', 'RPA')):
    plot_anat(t1_nii_fname, axes=axs[point_idx],
              cut_coords=mri_pos[point_idx, :],
              title=label, vmax=160)
plt.show()

###############################################################################
# Writing defaced and anonymized T1 image
# ---------------------------------------
#
# We can deface the MRI for anonymization by passing ``deface=True``.
t1w_bids_path = write_anat(
    image='./croppedT1.nii', # path to the MRI scan
    bids_path=bids_path,
    raw=raw,  # the raw MEG data file connected to the MRI
    trans=trans,  # our transformation matrix
    deface=True,
    overwrite=True,
    verbose=True  # this will print out the sidecar file
)
anat_dir = t1w_bids_path.directory

# Our MRI written to BIDS, we got `anat_dir` from our `write_anat` function
t1_nii_fname = op.join(anat_dir, 'sub-01_ses-01_T1w.nii.gz')

# Plot it
fig, ax = plt.subplots()
plot_anat(t1_nii_fname, axes=ax, title='Defaced', vmax=160)
plt.show()