{

 "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

}