506 lines (505 with data), 13.9 kB

{
 "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
}