Visualization/dMRI Reconstruction Project.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# dMRI Data Reconstruction\n",
    "\n",
    "**Contributor:** Erjun Zhang\n",
    "\n",
    "**Update time:** June 12, 2020\n",
    "\n",
    "**Status:** Complete first phase, and still working on this project untill the end of 2020 year\n",
    "\n",
    "Welcome anyone, who has questions about this project, to contact the contributor through [email](erjun.zhang@polymtl.ca).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This project is about diffusion magnetic resonance (MR) data processing and analysis. It mainly consists of three parts: brain diffusion MR data preprocessing, diffusion MRI images reconstruction, data visualization and left and right hemispherical preprocessed MR images classification. The whole procedures can be found in [this Jupyter Notebook file](https://github.com/brainhack-school2020/BHS_Project_dMRI/blob/master/dMRI%20Reconstruction%20Project.ipynb). Explanations about procedures results and other details are given in it.\n",
    "\n",
    "With reproducibility being a primary concern, this project was completed by using open-source softwares/tools (Python, FSL, DIPYPE...) and dataset (dHCP and PRIME). It can be as a simple tutorial/example for new students in neuroscience (diffusion MR brain imaging) to familiar/start with their further neuroscience study. \n",
    "\n",
    "Goals:\n",
    "\n",
    "* **Data Preprocesssing:** Get preprocessed diffusion MR images from raw MRI data;\n",
    "* **Image Reconstruction:** Reconstruct diffusion tensor images from the preprocessed data;\n",
    "* **Machine Learning Classsification:** By using machine learning, randomly choose left or right hemispheres on each slice of MRI images preprocessed above and find method to give results that which part it is belonging to (right or left?). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Preprocessing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data preprocessing is quit important for dMRI reconstruction. Different data preprocessing may lead to different reconstruction image qualities, which will make the comparation of different reconstruct methods unreliable. Thus, here we first preprocessing MRI by following same steps: denosing, topup (susceptibility-induced distortion correction) and eddy current-induced distortion and motion correction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### __Import python libraries__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os #TO control directories\n",
    "import nibabel as nib # read and save medical images\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "%matplotlib inline\n",
    "import plotly\n",
    "import plotly.express as px\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "from skimage import io # image/video read and show\n",
    "\n",
    "import timeit #compute time, useage: timeit.timeit()\n",
    "import math\n",
    "import time\n",
    "import warnings\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import nipype.interfaces.fsl as fsl #topup\n",
    "from nipype.interfaces.fsl import TOPUP\n",
    "from nipype.interfaces.fsl import ApplyTOPUP\n",
    "from nipype.interfaces.fsl import Eddy\n",
    "from nipype.testing import anatfile\n",
    "\n",
    "from dipy.denoise.localpca import mppca #denoising\n",
    "from dipy.io import read_bvals_bvecs\n",
    "from dipy.core.gradients import gradient_table\n",
    "from dipy.reconst.dti import TensorModel\n",
    "from dipy.reconst.dti import fractional_anisotropy\n",
    "from dipy.reconst.dti import color_fa\n",
    "import dipy.reconst.dki as dki\n",
    "\n",
    "\n",
    "# Self-defined functions\n",
    "from Extract_b0_Image import Extract_b0_Image\n",
    "from nib_rdshow_img import nib_rdshow_img\n",
    "from nib_read_img import nib_read_img\n",
    "from nib_show_img import nib_show_img\n",
    "from vol_plot import vol_plot\n",
    "from interact_vol_plot import interact_vol_plot\n",
    "\n",
    "# Import libararies used to data analysis\n",
    "import seaborn as sns\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import classification_report,confusion_matrix\n",
    "from sklearn.neighbors import KNeighborsClassifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set data path and load data\n",
    "\n",
    "This project used data from online dataset offered by:\n",
    "1. [The Developing Human Connectome Project](http://www.developingconnectome.org/second-data-release/). It consists of over 800 neonatal scans and over 250 fetal scans and can be used for data analysis after image reconstruction. \n",
    "2. [PRIME](http://fcon_1000.projects.nitrc.org/indi/PRIMEdownloads.html): used this dataset to reconstruct diffusion Images. This also can be downloaded [here](https://drive.google.com/file/d/1zgxynxjUCETBC6MAl4rfh0sL0WhFtKA9/view?usp=sharing) directly. \n",
    "3. Since during preprocessing, we used epi data with two opposite phase-encoding directions to correct distortions, other data can also be used as the source data if it meets this requirement.\n",
    "4. Data used for analysis is generated from image data after preprocessing. The final data can be found in this [github project folder](https://github.com/brainhack-school2020/BHS_Project_dMRI/blob/master/brain_water.csv)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data type (NII.GZ)\n",
    "\n",
    "- NII (NifTI, format by Neuroimaging Informatics Technology Initiative ) file is commonly used format for multi-dimensional (can be up to 7-dimensional) neuroimaging data. Fisrt four dimension: spatial dimensions and time. \n",
    "\n",
    "- GZ means gzip-compressed NII files.\n",
    "* nib.nifil.Nifti1Image: three parts included in, namely, image data array, an affine array and image metadata.\n",
    "* image metadata: machine info., voxel size and slices\n",
    "\n",
    "Thus, in order to know the exact position of each voxel, we have to combine image data array and affine array. For more information, please check [fMRI Processing based on python](https://ff120.github.io/2016/06/12/%E8%AE%A4%E7%9F%A5%E7%A5%9E%E7%BB%8F%E7%A7%91%E5%AD%A6%E4%B8%93%E9%A2%98/%E4%BD%BF%E7%94%A8Python%E5%A4%84%E7%90%86fMRI%E6%95%B0%E6%8D%AE/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### __Denoising__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we use [Marcenko-Pastur PCA algorithm](https://dipy.org/documentation/1.0.0./examples_built/denoise_mppca/) to denoise images. This algorithm has been shown to provide an optimal compromise between noise suppression and loss of anatomical information for different techniques such as DTI.\n",
    "\n",
    "During the denoising, mppca use a 3D sliding window (decised by denoising radius, pathc_radius) to denoise. Basicaly, this 3D sliding window voxles should be larger than DTI volumes.\n",
    "\n",
    "* Input: AP_file and PA_file, patch radius\n",
    "* Output: 'Denoised_AP_dwi.nii.gz' and 'Denoised_AP_dwi.nii.gz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extract b0 images\n",
    "\n",
    "* Input file: ap_file, or pa_file, or denoised_AP_file, or denoised_PA_file\n",
    "* Corresponding extract 3D file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Merge images\n",
    "\n",
    "* Input files: 'Extract_denoised_AP_b0.nii.gz' +'Extract_denoised_PA_b0.nii.gz'\n",
    "* Output files: 'Extract_denoised_AP_b0_merged.nii.gz', 'topup_encoding.txt'\n",
    "\n",
    "While edit topup encoding file, phase encoded direciton can be found in .json file. The fourth parameter is the [time duration](file:///home/erjun/Downloads/topup(2f)ExampleTopupFollowedByApplytopup.html) between the readout of the centre of the first echo and the centre of the last echo. It can also be found in .json file or use parameters in .json file to calculate it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### TOPUP\n",
    "* Input files: 'Extract_denoised_AP_b0_merged.nii.gz' +'topup_encoding.txt'\n",
    "* Output files: 'Extract_denoised_AP_b0_merged_base_fieldcoef.nii.gz', 'Extract_denoised_AP_b0_merged_base_movpar.txt', 'Extract_denoised_AP_b0_merged_corrected.nii.gz',Extract_denoised_AP_b0_merged_field.nii.gz, Extract_denoised_AP_b0_merged_topup.log, jac_01.nii.gz, jac_02.nii.gz, warpfield_01.nii.gz, warpfield_02.nii.gz, xfm_01.mat, xfm_02.mat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ApplyTOPUP\n",
    "* Input files: 'Extract_denoised_AP_b0.nii.gz','Extract_denoised_PA_b0.nii.gz'\n",
    "* Output files: Extract_denoised_AP_b0_merged_corrected.nii.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### bet\n",
    "\n",
    "* Input files: 'Extract_denoised_AP_b0_corrected.nii.gz'\n",
    "* Output files: bet_brain.nii.gz, bet_brain_mask.nii.gz, index.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EDDY\n",
    "\n",
    "During MRI scannning, subject movements and eddy current-induced distortions may occur. These distortion can be corrected by using FSL.Eddy. \n",
    "* Input files: denoised_AP_file, bet_brain_mask.nii.gz, index.txt\n",
    "* Output files: eddy_corrected_AP.nii.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DTI Reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, I first show the basic images generate during preprocessing and final image reconstruction process. Then I will go the data visualization part. And before that, to create images more eassily, I will definie several image showing function first.\n",
    "\n",
    "All the images are generated from this project and base on the data preprocessing and iamge reconstruction part.\n",
    "\n",
    "Let's go to check what basic images we aready have!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Basic Images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# DTI Images\n",
    "# set plot background\n",
    "#plt.style.use('seaborn-dark')\n",
    "plt.style.use('grayscale')\n",
    "# plot paramter maps  \n",
    "axf1=np.fliplr(DTI_RGB[:,:,28,:])\n",
    "#axf1=axf1.T\n",
    "axf2=np.fliplr(DTI_MD[:,:,28])\n",
    "axf2=axf2.T\n",
    "axf3=np.fliplr(DTI_RD[:,:,28])\n",
    "axf3=axf3.T\n",
    "axf4=np.fliplr(DTI_AD[:,:,28])\n",
    "axf4=axf4.T\n",
    "fig, [ax0,ax2, ax3, ax4] = plt.subplots(1,4,figsize=(12,10),subplot_kw={'xticks': [], 'yticks': []})\n",
    "ax0.imshow(axf1); ax0.set_title('Color coded FA',fontweight='bold',size=13)\n",
    "#ax1.imshow(DTI_FA[:,:,28]); ax1.set_title('Fractional anisotropy',fontweight='bold',size=13)\n",
    "ax2.imshow(axf2); ax2.set_title('Mean diffusivity',fontweight='bold',size=13)\n",
    "ax3.imshow(axf3); ax3.set_title('Radial diffusivity',fontweight='bold',size=13)\n",
    "ax4.imshow(axf4); ax4.set_title('Axial diffusivity',fontweight='bold',size=13)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# DKI images  \n",
    "axt0=np.fliplr(DKI_AD[:,:,28])\n",
    "axt0=axt0.T\n",
    "\n",
    "axt1=np.fliplr(DKI_RD[:,:,28])\n",
    "axt1=axt1.T\n",
    "\n",
    "axt2=np.fliplr(DKI_MD[:,:,28])\n",
    "axt2=axt2.T\n",
    "\n",
    "axt3=np.fliplr(DKI_AK[:,:,28])\n",
    "axt3=axt3.T\n",
    "\n",
    "axt4=np.fliplr(DKI_RK[:,:,28])\n",
    "axt4=axt4.T\n",
    "\n",
    "axt5=np.fliplr(DKI_MK[:,:,28])\n",
    "axt5=axt5.T\n",
    "\n",
    "fig, ([ax0, ax1, ax2],[ax3, ax4, ax5]) = plt.subplots(2,3,figsize=(10,8),subplot_kw={'xticks': [], 'yticks': []})\n",
    "ax0.imshow(axt0); ax0.set_title('Axial diffusivity',fontweight='bold',size=10)\n",
    "ax1.imshow(axt0); ax1.set_title('Radial diffusivity',fontweight='bold',size=10)\n",
    "ax2.imshow(axt0); ax2.set_title('Mean diffusivity',fontweight='bold',size=10)\n",
    "ax3.imshow(axt0); ax3.set_title('Axial kurtosis',fontweight='bold',size=10)\n",
    "ax4.imshow(axt0); ax4.set_title('Radial kurtosis',fontweight='bold',size=10)\n",
    "ax5.imshow(axt0); ax5.set_title('Mean kurtosis',fontweight='bold',size=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function, used to show images, definition\n",
    "# To read images data from nii.gz file: nib_read_img(path)\n",
    "# To show images from nii.gz file: nib_show_img(img0,slices,intenseScale)\n",
    "# To read and show images from nii.gz file: nib_rdshow_img(Images,Slices,IntenseScale,TitleImg)\n",
    "# Function to show volume slices images: vol_plot(x) and interact_vol_plot(x,IntenseScale):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3D volume slice image\n",
    "\n",
    "This images will quickly check the general apperance of you images during the slice-by-slice animation. One can also rotate it and zoom in to see clearly.\n",
    "\n",
    "Here, as an example, I just show FA map generated from DTI model. Once can change it to other maps, such as DTI_MD or DKI_FA.\n",
    "\n",
    "  \t1. Show different slices of images;\n",
    "  \t2. Change observe view point;\n",
    "  \t3. Play animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#vol_plot(DTI_FA[:,:,:]) # or can use interact_vol_plot(DTI_MD,90) to show images\n",
    "interact_vol_plot(DTI_FA,0.90)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data visualization (Preprocessing part)\n",
    "\n",
    "I am curious about what changed of the iamges during our data preprocessing. This interactive images can be used to reach this purpose.\n",
    "\n",
    "One can control Slices to check different slice changes among them. You can control the level of high color scale to get higher image contrast. You can put cursor on images to get the exact color value if you want.\n",
    "\n",
    "  \t1. Select different maps;\n",
    "  \t2. Select slices;\n",
    "  \t3. Select maximum diffusion sthrength;\n",
    "  \t4. Display title"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set image data for visualization (preprocessing part)\n",
    "# This will be used for visualization of preprocessing part 1 and part 2\n",
    "F11='Extract_AP_b0.nii.gz'\n",
    "F12='Extract_PA_b0.nii.gz'\n",
    "F13='Extract_denoised_AP_b0.nii.gz'\n",
    "F14='Extract_denoised_PA_b0.nii.gz'\n",
    "F15='Extract_denoised_AP_b0_merged1.nii.gz'\n",
    "F16='Extract_denoised_AP_b0_merged_corrected1.nii.gz'\n",
    "F17='bet_brain.nii.gz'\n",
    "F18='bet_brain_mask.nii.gz'\n",
    "F19='Extract_eddy_correct_AP.nii.gz'\n",
    "\n",
    "\n",
    "PreproImage= dict()\n",
    "PreproImage['Original_AP_Image']=F11\n",
    "PreproImage['Original_PA_Image']=F12\n",
    "PreproImage['Denoised_AP_Image']=F13\n",
    "PreproImage['Denoised_PA_Image']=F14\n",
    "PreproImage['Merged_Image']=F15\n",
    "PreproImage['TOPUP_AP_Image']=F16\n",
    "PreproImage['Brain_After_bet']=F17\n",
    "\n",
    "PreproImage['Brain_mask']=F18\n",
    "PreproImage['Eddy_Correct_AP_Image']=F19\n",
    "\n",
    "# Load DTI and DKI results\n",
    "# \n",
    "Images= dict()\n",
    "Images['DTI_FA_RGB']=(data_path+'/DTI/FA_RGB.nii.gz')\n",
    "Images['DTI_AD']=(data_path+'/DTI/AD.nii.gz')\n",
    "Images['DTI_FA']=data_path+'/DTI/FA.nii.gz'\n",
    "Images['DTI_MD']=data_path+'/DTI/MD.nii.gz'\n",
    "Images['DTI_RD']=data_path+'/DTI/RD.nii.gz'\n",
    "# DKI\n",
    "Images['DKI_AD']=data_path+'/DKI/dki_AD.nii.gz'\n",
    "Images['DKI_FA']=data_path+'/DKI/dki_FA.nii.gz'\n",
    "Images['DKI_MD']=data_path+'/DKI/dki_MD.nii.gz'\n",
    "Images['DKI_RD']=data_path+'/DKI/dki_RD.nii.gz'\n",
    "Images['DKI_AK']=data_path+'/DKI/AK.nii.gz'\n",
    "Images['DKI_MK']=data_path+'/DKI/MK.nii.gz'\n",
    "Images['DKI_RK']=data_path+'/DKI/RK.nii.gz'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preprocessing image visualization 1\n",
    "interactive(nib_rdshow_img,Images=PreproImage,Slices=widgets.IntSlider(min=0,max=43,step=1,value=28),\\\n",
    "         IntenseScale=widgets.IntSlider(min=0,max=100,step=1,value=90),\\\n",
    "           TitleImg='Reconstructed Images of Different Models')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Preprocessing image visualization 2\n",
    "\n",
    "1. play different maps;\n",
    "2. Select slices;\n",
    "3. Select maximum diffusion sthrength;\n",
    "4. Display title"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preprocessing image visualization 2\n",
    "# Function below is used to show correct figure title while animaiton playing\n",
    "def nib_rdshow_play(Slices,IntenseScale,NoX):\n",
    "    warnings.filterwarnings('ignore')\n",
    "    if(NoX==0):\n",
    "        Image=F11\n",
    "        TitleImg='Original_AP_Image'\n",
    "    elif(NoX==1):\n",
    "        Image=F12 #'epi_b0_merged.nii.gz'\n",
    "        TitleImg='Original_PA_Image'\n",
    "    elif(NoX==2):\n",
    "        Image=F13\n",
    "        TitleImg='Denoised_AP_Image'\n",
    "    elif(NoX==3):\n",
    "        Image=F14\n",
    "        TitleImg='Denoised_PA_Image'\n",
    "    elif(NoX==4):\n",
    "        Image=F15 #'epi_b0_merged.nii.gz'\n",
    "        TitleImg='Merged_Image'\n",
    "    elif(NoX==5):\n",
    "        Image=F16\n",
    "        TitleImg='TOPUP_AP_Image'\n",
    "    elif(NoX==6):\n",
    "        Image=F17\n",
    "        TitleImg='Brain_After_bet'\n",
    "    elif(NoX==7):\n",
    "        Image=F18\n",
    "        TitleImg='Brain_mask'\n",
    "    else:\n",
    "        Image=F15\n",
    "        TitleImg='Eddy_Correct_AP_Image'\n",
    "    image_data = nib.load(Image).get_data()\n",
    "    img00=image_data[:,:,Slices]\n",
    "    zmax0=img00.max()\n",
    "    fig=px.imshow(image_data[:,:,Slices],color_continuous_scale=\"Viridis\",\\\n",
    "                  zmin=0,zmax=zmax0*IntenseScale/100,\\\n",
    "                  labels={},template=\"plotly_white\")\n",
    "    fig.update_xaxes(showticklabels=False)\n",
    "    fig.update_yaxes(showticklabels=False)\n",
    "    fig.update_layout(coloraxis_showscale=False)\n",
    "    fig.update_layout(title=TitleImg)\n",
    "    fig.show()\n",
    "\n",
    "interactive(nib_rdshow_play,Slices=widgets.IntSlider(min=0,max=43,step=1,value=28),\\\n",
    "         IntenseScale=widgets.IntSlider(min=0,max=100,step=1,value=90),\\\n",
    "           NoX = widgets.Play(value=0,min=0,max=8,step=1,interval=2000,description=\"Press play\",\\\n",
    "                                   disabled=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data visualization (Reconstruction part)\n",
    "\n",
    "This part is used to compare images generated from DTI and DKI models. The founction is almost the same as virsualizaiotn above. But to make us better choose iamges, we use dropdown bar to replace slider animation.\n",
    "\n",
    "1. Play different maps;\n",
    "2. Select slices;\n",
    "3. Slect maximum diffusion sthrength;\n",
    "4. Display give title"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reconstruction image visualization 1\n",
    "interactive(nib_rdshow_img,Images=Images,Slices=widgets.IntSlider(min=0,max=43,step=1,value=22),\\\n",
    "         IntenseScale=widgets.IntSlider(min=0,max=100,step=1,value=90),\\\n",
    "           TitleImg='Reconstructed Images of Different Models')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Analysis\n",
    "\n",
    "1. 80 Slices of brain images (40 left hemisphere slice and 40 right hemisphere) were extracted and transfered to brain parts classification;\n",
    "2. Left hemisphere brain image and right hemisphere brain iamge were taged with 0 and 1 specifically;\n",
    "3. Slice number was treated as one of features used for classification;\n",
    "4. Diffusino strength of each slice hemisphere was calculated  as the second feature of classification;\n",
    "5. Total effictive voxel of eahc slice hemisphere was calculated as the third feature for classification;\n",
    "6. Dataset was splitted into train dataset (70%) and test dataset (30%);\n",
    "7. Use KNN to classify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import libararies used to data analysis\n",
    "import seaborn as sns\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import classification_report,confusion_matrix\n",
    "from sklearn.neighbors import KNeighborsClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shapeT18=np.shape(T18)\n",
    "data_MR=np.zeros([80,4])\n",
    "print('MR data shape:',shapeT18)\n",
    "print('dMRI data shape:',np.shape(data_MR))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate data of right heimesphere and left heimesphere\n",
    "for i in range(0,39):\n",
    "    data_MR[i,0]=i\n",
    "    data_MR[i+39,0]=i\n",
    "    data_MR[i,3]=0\n",
    "    data_MR[i+39,3]=1\n",
    "    \n",
    "    EDDY_valueA=T19[0:47,:,i]\n",
    "    EDDY_valueA=np.sum(EDDY_valueA)\n",
    "    data_MR[i,2]=5*np.log2(1+EDDY_valueA)\n",
    "    EDDY_valueB=T19[48:95,:,i]\n",
    "    EDDY_valueB=5*np.sum(EDDY_valueB)\n",
    "    data_MR[i+39,2]=np.log2(1+EDDY_valueB)\n",
    "    \n",
    "    Num_voxA=T18[0:47,:,i]\n",
    "    Num_voxA=np.sum(Num_voxA)\n",
    "    data_MR[i,1]=Num_voxA\n",
    "    Num_voxB=T18[48:95,:,i]\n",
    "    Num_voxB=np.sum(Num_voxB)\n",
    "    data_MR[i+39,1]=Num_voxB\n",
    "np.savetxt('brain_water.csv',data_MR,delimiter=',')\n",
    "brain_water_csv = pd.read_csv('brain_water.csv',names=['Num_slice','Num_vox','Semi_AB','Position'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Details of data waiting for analyzing\n",
    "brain_water_csv=brain_water_csv.dropna() # Delete NAN data line\n",
    "brain_water_csv.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Details of data waiting for analyzing\n",
    "brain = brain_water_csv.describe()\n",
    "brain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#sns.pairplot(brain_water_csv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Checek frequency of right and left hemispheres\n",
    "brain_water_csv['Position'].plot.hist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check brain diffusion water data location\n",
    "brain_water_csv['Semi_AB'].plot.density()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set train dataset (catergoricalfetures and target class)\n",
    "x1=(brain_water_csv['Num_slice'])\n",
    "x2=(brain_water_csv['Num_vox'])\n",
    "x3=(brain_water_csv['Semi_AB'])\n",
    "x=np.vstack((x1,x2,x3))\n",
    "x=x.T\n",
    "y=5*(brain_water_csv['Position'])\n",
    "y=y.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split dataset into train dataset and test dataset\n",
    "x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.3, random_state=123)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use KNN with different K value to train and test\n",
    "error_rate = []\n",
    "for i in range(1,40):    \n",
    "    knn = KNeighborsClassifier(n_neighbors=i)\n",
    "    knn.fit(x_train,y_train)\n",
    "    pred_i = knn.predict(x_test)\n",
    "    error_rate.append(np.mean(pred_i != y_test))\n",
    "\n",
    "plt.figure(figsize=(12,9))\n",
    "plt.plot(range(1,40),error_rate,color='gray', linestyle='dashed', marker='o',\n",
    "         markerfacecolor='black', markersize=5)\n",
    "plt.title('Error Rate vs. K Value', fontsize=20)\n",
    "plt.xlabel('K value',fontsize=14)\n",
    "plt.ylabel('Error Rate',fontsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use KNN to train and test\n",
    "knn = KNeighborsClassifier(n_neighbors=1)\n",
    "knn.fit(x_train,y_train)\n",
    "pred = knn.predict(x_test)\n",
    "conf_mat=confusion_matrix(y_test,pred)\n",
    "print('Conf_mat',conf_mat)\n",
    "\n",
    "# Show clasification report\n",
    "print(classification_report(y_test,pred))\n",
    "print(\"Misclassification error rate:\",round(np.mean(pred!=y_test),3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "This project started with medical imaging data format, then preprocessed the diffusion weighted images, which includes MP-PCA denoising, FSL TOPUP distortion correction and head movement correction. After this, DTI images were reconstructed from these preprocessed images successfully. Additionally, a dataset for left and right hemispheres classification was generated from these preprocessed images. By using the KNN method, classification accuracy of 92% was reached. As a starter, I will continue this project in the future. Soon, self-made models, instead of DTI model, will be used to reconstruct diffusion images. After that, I would like to use machine learning methods to classify potential changes of brain microstructures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}