CoCalc -- seminar2_1

CoCalc provides the best real-time collaborative environment for Jupyter Notebooks, LaTeX documents, and SageMath, scalable from individual users to large groups and classes!

GitHub Repository: adasegroup/NEUROML2022
Path: blob/main/seminar2/seminar2_1_mri.ipynb
Views: ⁶³

Kernel: Python 3 (ipykernel)

MRI data analysis, sources databases, tools and QC

PLAN:

MRI data formats \
accessing data voxel and meta data \
transformations \
preprocessing \
visualization

In [ ]:

from pydicom.data import get_testdata_files
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
from ipywidgets import interact, interactive, fixed, interact_manual
from IPython.display import Image
import ipywidgets as widgets
import nibabel.freesurfer.mghformat as mgh
import nibabel as nib
import pydicom
from glob import glob
import os
import numpy as np
from nilearn.input_data import NiftiMasker
from nilearn import plotting
%matplotlib inline

1) MRI File Formats

All volume-based formats store 3D or 4D arrays of voxels in some fashion with a variety of additional meta-data. Anatomical images are typically 3D while EPIs are typically 4D (x,y,z, and time).

In [ ]:

Image(filename = "/workspace/assets/voxel.png", width=1000, height=800)

voxel: A three-dimensional pixel and the basic unit of spatial measurement in MRI.

volume: The three dimensional array covering the brain. Volumes are composed of voxels All volume files contain both meta-data and voxels. The meta-data is just a set of information about the file’s contents while the voxels are a 3D or 4D array of values.

slice: A two-dimensional 'view' of the three-dimensional volume obtained by taking all of the elements in two of the dimensions for a fixed location in the third dimension.

Timecourse or Timeseries: a set of numbers representing a measurement (like BOLD activation) taken over time.

Here are the most typical file formats of neuroimaging data:

DICOM (.dcm, .ima) \
NIFTI (.nii/nii.gz, gii.gz) \
MGH (.mgh, .mgz) \
Custom (.dtseries, .label, .surf)

DICOM file format

In [ ]:

# Single slice of Brain MRI

dcm_path = '/workspace/data/Brats_kaggle/DICOM/00008/T1w/'
g = glob(os.path.join(dcm_path, '*'))
!tree ./data/Brats_kaggle/DICOM/ --filelimit 2

# Print out the first 5 file names to verify we're in the right folder.

print('Total of %d DICOM images.\nFirst 5 filenames:' % len(g))
print('Standart DICOM directory tree: \n')
print('\n'.join(g[:]))

Common MRI DICOM sctructure

In [ ]:

dcm_file = pydicom.dcmread(os.path.join(dcm_path, 'Image-20.dcm'))
plt.imshow(dcm_file.pixel_array, cmap='gray', origin='upper')
rows = int(dcm_file.Rows)
cols = int(dcm_file.Columns)
print('Dimension of one dcm file:', dcm_file.pixel_array.shape)
plt.show()

Lets stack all DICOM's and plot em aall:

In [ ]:

def load_scan(path):
    slices = [pydicom.dcmread(dcm_path + '/' + s) for s in os.listdir(dcm_path)]
    slices.sort(key = lambda x: int(x.InstanceNumber))
    try:
        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] - slices[1].ImagePositionPatient[2])
    except:
        slice_thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation)
        
    for s in slices:
        s.SliceThickness = slice_thickness
    
    image = np.stack([s.pixel_array for s in slices])
    return np.array(image)

In [ ]:

dcm_stack = load_scan(dcm_path)
rows=3
cols=3
start_with=5
show_every=3
fig,ax = plt.subplots(rows,cols,figsize=[15,15])
for i in range(rows*cols):
    ind = start_with + i*show_every
    ax[int(i/rows),int(i % rows)].imshow(dcm_stack[ind],cmap='gray', origin="upper")
    ax[int(i/rows),int(i % rows)].set_title('slice %d' % ind)
    ax[int(i/rows),int(i % rows)].set_xlabel('xlabel')
    ax[int(i/rows),int(i % rows)].set_ylabel('ylabel')
plt.show()

DICOM Headers

Another good way to look at the meta-data in a volume file is to load it with the relevant programming environment and examine the data-structures there. Here are a few examples. Accordingly, we need to be able to, at a minimum, store some amount of information about the coordinate system employed in any MRI volume file, and ideally some amount of information about how to precisely align the brain to some standard orientation.

In [ ]:

pat_name = dcm_file.PatientName
display_name = pat_name.family_name + ", " + pat_name.given_name
print("Patient id.......:", dcm_file.PatientID)
print("Modality.........:", dcm_file.Modality)
print("Rows.............:", dcm_file.Rows)
print("Columns..........:", dcm_file.Columns)
print("Pixel Spacing....:", dcm_file.PixelSpacing)
print("Slide Thickness.:", dcm_file.SliceThickness)
print("Patient Position.:", dcm_file.PatientPosition)

In [ ]:

#dcm_header = pydicom.dcmread(dcm_file)
dcm_file

Data manipulations

Assume that you recieved new dicom from hospital or downloaded challenge data and you need to convert to convinient NIFTI, let's review the tool distinctive from simple one represented above:

dcm2niix Robust .dcm to nii/nii.gz converter. Params %i %n %p stands for patient ID, subject name, protocol name DICOM header keys. Try to run with params %p_%t_%s \

https://www.nitrc.org/plugins/mwiki/index.php/dcm2nii:MainPage#General_Usage

In [ ]:

!dcm2niix -o /workspace/ -f %i_%n_%p -z y ./data/Brats_kaggle/DICOM/00008/T1w

Plotting results of convertation:

In [ ]:

plotting.plot_anat('/workspace/00008_00008_T1w.nii.gz', bg_img=None)

Nifti File Format

The Neuroimaging Informatics Technology Initiative (nifti) is format to store radiological information, the first three dimensions are reserved to define the three spatial dimensions — x, y and z —, while the fourth dimension is reserved to define the time points — t. The remaining dimensions, from fifth to seventh, are for other uses. The fifth dimension, however, can still have some predefined uses, such as to store voxel-specific distributional parameters or to hold vector-based data.

In [ ]:

nii_file = nib.load('/workspace/00008_00008_T1w.nii.gz')

Nifti Header

In [ ]:

nii_header = nii_file.header
print("Image dimensions...........................:", nii_header['dim'])
print("Voxel size and time interval...............:", nii_header['pixdim'])
print("Slice Order................................:", nii_header['slice_code'])
print("Human readable text........................:", nii_header['descrip'])
print("Three rows of sform affine transformation..:"'\n', nii_header['srow_x'],'\n', nii_header['srow_y'],'\n', nii_header['srow_z'])

MGH File Format

MGH File Format is a default format of freesurfer preprocessing software. That can do the following

In [ ]:

Image(filename = "/workspace/assets/freesurfer.png", width=1000, height=800)

In [ ]:

mgh_file = mgh.load('./data/freesurfer_preproc/100206/mri/T1.mgz')
mgh_file.get_fdata().shape

In [ ]:

# Define plot function
def show_slices(image, axis1="x", axis2="y", axis3="z"):
    slice_0 = image[round(len(image[0])/2), :, :]
    slice_1 = image[:, round(len(image[1])/2), :]
    slice_2 = image[:, :, round(len(image[2])/2)]
    image = ([slice_0, slice_1, slice_2])
    fig, axes = plt.subplots(1, len(image), figsize=[15,15])
    for i, slice in enumerate(image):
        axes[i].imshow(slice.T, cmap="gray", origin="upper")
        axes[0].set(xlabel=axis2, ylabel=axis3)
        axes[1].set(xlabel=axis1, ylabel=axis3)
        axes[2].set(xlabel=axis1, ylabel=axis2)
    plt.show()

In [ ]:

mgh_image = np.array(mgh_file.dataobj).astype(np.float64)
show_slices(mgh_image)

MGH Header

In [ ]:

mgh_header = mgh_file.header
print("Image dimensions........:", mgh_header['dims'])
print("Orientation matrix...: "'\n', mgh_header.get_affine())

In [ ]:

plotting.plot_anat(mgh_file, bg_img=None)

In [ ]:

Image(filename = "/workspace/assets/Orientation.png", width=600, height=800)

DICOM patient coordinate and Nifti coodrinate system when Patient is Head First Supine. The arrows of axes indicate the positive directions, LPS for DICOM and RAS for Nifti. For Dicom images, the origin (0,0,0) refers to magnet isocenter.

Voxel size and Pixel spacing are important as we for preprocessing we can use only isomorphic voxel

BIDS Brain Imaging Data Structure. Simplify representation of DICOM neuroimaging data and grants huge advantages for experiments and version control of dataset. Moreover, reduces dataset storage range.

In [ ]:

Image(filename = "/workspace/assets/BIDS.png", width=800, height=800)

There is a few robust DICOM converters:

HeudiConv

dcm2bids

They both use DICOM headers, to distinquish sequence protocols. So fully anonomized DICOMs can't be converted using this software. Let's look how heudiconv runs: First you should implement heudiconv.py. The only required function for a heuristic, infotodict is used to both define the conversion outputs and specify the criteria for scan to output association. Conversion outputs are defined as keys, a tuple consisting of a template path used for the basis of outputs, as well as a tuple of output types. Valid types include nii, nii.gz, and dicom.

In [ ]:

!/opt/miniconda-latest/envs/neuro/bin/heudiconv -d /workspace/data/Brats_kaggle/DICOM/{subject}/*/* \
    -f /workspace/heuristic.py \
    -s 00008 \
    -c dcm2niix \
    -b --overwrite \
    -o /workspace/BIDS > /workspace/heudiconv_log 2>&1

Heudiconv output BIDS directory tree:

In [ ]:

!tree /workspace/BIDS

In [ ]:

!rm -r /workspace/BIDS/*
!rm -r /workspace/BIDS/.heudiconv

In [ ]: