Datasets
Models
Downloads: 1
Diff of
/.ipynb_checkpoints/Untitled-Copy1-checkpoint.ipynb
[000000]
..
[4cda31]
Switch to side-by-side view
--- a +++ b/.ipynb_checkpoints/Untitled-Copy1-checkpoint.ipynb @@ -0,0 +1,505 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import h5py\n", + "import imageio\n", + "import numpy as np\n", + "import pandas as pd\n", + "import nibabel as nib\n", + "import seaborn as sns\n", + "from utils import metrics\n", + "from data import base_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "from utils import visualizer\n", + "import matplotlib.pylab as plt\n", + "from scipy.ndimage.measurements import center_of_mass" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas\n", + "from matplotlib.figure import Figure" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DM 0.27433953948483814 -0.17179240013074187\n" + ] + } + ], + "source": [ + "Clinical_params = pd.DataFrame({'Group':[], 'SubjectID':[], 'rESS':[], 'cESS':[],\n", + " 'RV_EDV_ml':[], 'RV_ESV_ml':[], 'RV_EF':[], \n", + " 'LV_EDV_ml':[], 'LV_ESV_ml':[], 'LV_EF':[], 'LV_mass_g':[], 'Frame':[]})\n", + "\n", + "\n", + "df = pd.read_csv('../private_data/MARTINOS.csv', index_col=0)\n", + "\n", + "Err = []\n", + "Ecc = []\n", + "for subject_index in df.index: \n", + " group = df.iloc[subject_index].Group\n", + " sid = df.iloc[subject_index].SubjectID\n", + " \n", + " if group == 'DM':\n", + " M_nifti = nib.load('../private_data_results/MARTINOS/MSTAT_DM_VOL%d_V1_segmentation.nii'%(sid))\n", + " u_HF = h5py.File('../private_data_results/MARTINOS/MSTAT_DM_VOL%d_V1_motion.h5'%(sid), 'r')\n", + " elif sid > 100:\n", + " M_nifti = nib.load('../private_data_results/MARTINOS/MSTAT_%d_V2_segmentation.nii'%(sid))\n", + " u_HF = h5py.File('../private_data_results/MARTINOS/MSTAT_%d_V2_motion.h5'%(sid), 'r')\n", + " else:\n", + " try:\n", + " M_nifti = nib.load('../private_data_results/MARTINOS/MSTAT_VOL%d_V1_segmentation.nii'%(sid))\n", + " u_HF = h5py.File('../private_data_results/MARTINOS/MSTAT_VOL%d_V1_motion.h5'%(sid), 'r')\n", + " except:\n", + " M_nifti = nib.load('../private_data_results/MARTINOS/MSTAT_VOL%d_V2_segmentation.nii'%(sid))\n", + " u_HF = h5py.File('../private_data_results/MARTINOS/MSTAT_VOL%d_V2_motion.h5'%(sid), 'r')\n", + " \n", + " \n", + " M_nifti = nifti_dataset.resample_nifti(M_nifti, in_plane_resolution_mm=1.25, number_of_slices=16)\n", + " \n", + " m = M_nifti.get_fdata()\n", + " u = load_HF(u_HF)\n", + " \n", + " center = center_of_mass(m[:,:,:,0]==3)\n", + " \n", + " u=base_dataset._roll2center_crop(u,center)\n", + " m=base_dataset._roll2center_crop(m,center)\n", + " \n", + " esid = (m[:,:,3:-3]==3).sum(axis=(0,1,2)).argmin()\n", + " \n", + " \n", + " for t in [esid]:\n", + " # CALCULATE PARAMS AT END-SYSTOLE ONLY\n", + " \n", + " E = strain.MyocardialStrain(m[:,:,:,0],u[:,:,:,:,t])\n", + " E.calculate_strain()\n", + "\n", + " rESS = E.Err[E.mask_rot==2].mean()\n", + " cESS = E.Ecc[E.mask_rot==2].mean()\n", + " \n", + " print(group, rESS, cESS)\n", + "\n", + " \n", + " break" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "from scipy.ndimage import rotate\n", + "from scipy.ndimage import gaussian_filter\n", + "from scipy.interpolate import interp1d, interp2d\n", + "from scipy.ndimage.measurements import center_of_mass\n", + "\n", + "class MyocardialStrain():\n", + " \n", + " def __init__(self, mask, flow):\n", + " \n", + " self.mask = mask\n", + " self.flow = flow\n", + " \n", + " assert len(mask.shape) == 3\n", + " assert len(flow.shape) == 4\n", + " assert mask.shape == flow.shape[:3]\n", + " assert flow.shape[-1] == 3\n", + " \n", + " def calculate_strain(self, dx=1, dy=1, dz=1, lv_label=3):\n", + " \n", + " cx, cy, cz = center_of_mass(self.mask==lv_label)\n", + " nx, ny, nz = self.mask.shape\n", + " \n", + " self.flow_rot = _roll_to_center(self.flow, cx, cy)\n", + " self.mask_rot = _roll_to_center(self.mask, cx, cy)\n", + "\n", + " ux, uy, uz = np.array_split(self.flow_rot, 3, -1)\n", + " Uxx,Uxy,Uxz = np.gradient(np.squeeze(ux),dx,dy,dz)\n", + " Uyx,Uyy,Uyz = np.gradient(np.squeeze(uy),dx,dy,dz)\n", + " Uzx,Uzy,Uzz = np.gradient(np.squeeze(uz),dx,dy,dz)\n", + "\n", + " self.E_cart = np.zeros((nx,ny,nz,3,3))\n", + " for i in range(nx):\n", + " for j in range(ny):\n", + " for k in range(nz):\n", + " Ugrad = [[Uxx[i,j,k], Uxy[i,j,k], Uxz[i,j,k]], \n", + " [Uyx[i,j,k], Uyy[i,j,k], Uyz[i,j,k]],\n", + " [Uzx[i,j,k], Uzy[i,j,k], Uzz[i,j,k]]]\n", + " F = np.array(Ugrad) + np.identity(3)\n", + " e = 0.5*(np.matmul(F.T, F) - np.identity(3))\n", + " self.E_cart[i,j,k] += e\n", + "\n", + " self.Ezz = self.E_cart[:,:,:,2,2]\n", + " self.Err, self.Ecc = self._convert_to_polar(self.E_cart[:,:,:,:2,:2])\n", + "\n", + " \n", + " \n", + " \n", + " def _convert_to_polar(self, E):\n", + "\n", + " phi = _polar_grid(*E.shape[:2])[0]\n", + " Err = np.zeros(self.mask.shape)\n", + " Ecc = np.zeros(self.mask.shape)\n", + " for k in range(self.mask.shape[-1]):\n", + " cos = np.cos(np.deg2rad(phi))\n", + " sin = np.sin(np.deg2rad(phi))\n", + " \n", + " Exx, Exy, Eyx, Eyy = E[:,:,k,0,0],E[:,:,k,0,1],E[:,:,k,1,0],E[:,:,k,1,1]\n", + " Err[:,:,k] += cos*( cos*Exx+sin*Exy) +sin*( cos*Eyx+sin*Eyy)\n", + " Ecc[:,:,k] += -sin*(-sin*Exx+cos*Exy) +cos*(-sin*Eyx+cos*Eyy)\n", + "\n", + " return Err, Ecc\n", + " \n", + " \n", + "\n", + "def _roll(x, rx, ry):\n", + " x = np.roll(x, rx, axis=0)\n", + " return np.roll(x, ry, axis=1)\n", + "\n", + "def _roll_to_center(x, cx, cy):\n", + " nx, ny = x.shape[:2]\n", + " return _roll(x, int(nx//2-cx), int(ny//2-cy))\n", + "\n", + "def _polar_grid(nx=128, ny=128):\n", + " x, y = np.meshgrid(np.linspace(-nx//2, nx//2, nx), np.linspace(-ny//2, ny//2, ny))\n", + " phi = (np.rad2deg(np.arctan2(y, x)) + 180).T\n", + " r = np.sqrt(x**2+y**2+1e-8)\n", + " return phi, r" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DM 0.27433953948483814 -0.17179240013074187\n" + ] + } + ], + "source": [ + "E = MyocardialStrain(m[:,:,:,0],u[:,:,:,:,t])\n", + "E.calculate_strain()\n", + "\n", + "rESS = E.Err[E.mask_rot==2].mean()\n", + "cESS = E.Ecc[E.mask_rot==2].mean()\n", + "\n", + "print(group, rESS, cESS)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "mask = m[:,:,:,0]\n", + "flow = u[:,:,:,:,t]" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "lv_label = 3\n", + "\n", + "cx, cy, cz = center_of_mass(mask==lv_label)\n", + "nx, ny, nz = mask.shape\n", + "\n", + "cx, cy, cz\n", + "\n", + "flow_rot = _roll_to_center(flow, cx, cy)\n", + "mask_rot = _roll_to_center(mask, cx, cy)\n", + "\n", + "dx=1; dy=1; dz=1\n", + "\n", + "ux, uy, uz = np.array_split(flow_rot, 3, -1)\n", + "Uxx,Uxy,Uxz = np.gradient(np.squeeze(ux),dx,dy,dz)\n", + "Uyx,Uyy,Uyz = np.gradient(np.squeeze(uy),dx,dy,dz)\n", + "Uzx,Uzy,Uzz = np.gradient(np.squeeze(uz),dx,dy,dz)" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": {}, + "outputs": [], + "source": [ + "def_grad = np.array([[Uxx,Uxy,Uxz],\n", + " [Uyx,Uyy,Uyz],\n", + " [Uzx,Uzy,Uzz]])" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "metadata": {}, + "outputs": [], + "source": [ + "I = np.identity(3)[:,:,None,None,None]\n", + "I = np.repeat(I,repeats=128, axis=2)\n", + "I = np.repeat(I,repeats=128, axis=3)\n", + "I = np.repeat(I,repeats=16, axis=4)" + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "metadata": {}, + "outputs": [], + "source": [ + "F = def_grad + I" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 64, + "metadata": {}, + "outputs": [], + "source": [ + "for i in range(nx):\n", + " for j in range(ny):\n", + " for k in range(nz):\n", + " \n", + " Ugrad = [[Uxx[i,j,k], Uxy[i,j,k], Uxz[i,j,k]], \n", + " [Uyx[i,j,k], Uyy[i,j,k], Uyz[i,j,k]],\n", + " [Uzx[i,j,k], Uzy[i,j,k], Uzz[i,j,k]]]\n", + " \n", + " Fijk = np.array(Ugrad) + np.identity(3)\n", + " \n", + " assert (Fijk == F[:,:,i,j,k]).all()" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[1., 0., 0.],\n", + " [0., 1., 0.],\n", + " [0., 0., 1.]])" + ] + }, + "execution_count": 58, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Fijk" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[0., 0., 0.],\n", + " [0., 0., 0.],\n", + " [0., 0., 0.]])" + ] + }, + "execution_count": 59, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "F[:,:,i,j,k]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "F = np.array(Ugrad) + np.identity(3)\n", + " e = 0.5*(np.matmul(F.T, F) - np.identity(3))\n", + " self.E_cart[i,j,k] += e\n" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "metadata": {}, + "outputs": [], + "source": [ + "np.repeat?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# Manuel A. Morales (moralesq@mit.edu)\n", + "# Harvard-MIT Department of Health Sciences & Technology \n", + "# Athinoula A. Martinos Center for Biomedical Imaging\n", + "\n", + "import os\n", + "import h5py\n", + "import glob\n", + "import warnings\n", + "import numpy as np\n", + "import nibabel as nib\n", + "\n", + "from dipy.align.reslice import reslice\n", + "from data.base_dataset import BaseDataset, Transforms\n", + "from data.image_folder import make_dataset\n", + "\n", + "\n", + "class H5PYDataset(BaseDataset):\n", + "\n", + " def __init__(self, opt):\n", + " BaseDataset.__init__(self, opt)\n", + " self.filenames = sorted(make_dataset(opt.dataroot, opt.max_dataset_size, 'H5PY'))\n", + " \n", + " def __len__(self):\n", + " return len(self.filenames)\n", + " \n", + " def __getitem__(self, idx): \n", + "\n", + " HF = h5py.File(self.filenames[idx], 'r')\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "import h5py\n", + "from utils import strain\n", + "from data import nifti_dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def load_HF(HF):\n", + " output = []\n", + " for frame_id in range(len(HF.keys())):\n", + " key = 'frame_%d'%(frame_id)\n", + " for subkey in HF[key].keys():\n", + " output += [np.array(HF[key][subkey])]\n", + "\n", + " HF.close()\n", + " return np.stack(output,-1)" + ] + }, + { + "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.6.9" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}