metrics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df7898a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install hausdorff\n",
    "%pip install numba\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import glob\n",
    "import numpy as np\n",
    "import nibabel as nib\n",
    "from ACDCUNet import build_dict_images, build_dict_images_pred\n",
    "\n",
    "import statistics\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "from hausdorff import hausdorff_distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d744133d",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "author: Clément Zotti (clement.zotti@usherbrooke.ca)\n",
    "date: April 2017\n",
    "\n",
    "DESCRIPTION :\n",
    "The script provide helpers functions to handle nifti image format:\n",
    "    - load_nii()\n",
    "    - save_nii()\n",
    "\n",
    "to generate metrics for two images:\n",
    "    - metrics()\n",
    "\n",
    "And it is callable from the command line (see below).\n",
    "Each function provided in this script has comments to understand\n",
    "how they works.\n",
    "\n",
    "HOW-TO:\n",
    "\n",
    "This script was tested for python 3.4.\n",
    "\n",
    "First, you need to install the required packages with\n",
    "    pip install -r requirements.txt\n",
    "\n",
    "After the installation, you have two ways of running this script:\n",
    "    1) python metrics.py ground_truth/patient001_ED.nii.gz prediction/patient001_ED.nii.gz\n",
    "    2) python metrics.py ground_truth/ prediction/\n",
    "\n",
    "The first option will print in the console the dice and volume of each class for the given image.\n",
    "The second option wiil ouput a csv file where each images will have the dice and volume of each class.\n",
    "\n",
    "\n",
    "Link: http://acdc.creatis.insa-lyon.fr\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "HEADER = [\n",
    "    \"Name\",\n",
    "    \"Dice LV\",\n",
    "    \"Volume LV pred\",\n",
    "    \"Volume LV GT\",\n",
    "    \"Err LV(ml)\",\n",
    "    \"Dice RV\",\n",
    "    \"Volume RV pred\",\n",
    "    \"Volume RV GT\",\n",
    "    \"Err RV(ml)\",\n",
    "    \"Dice MYO\",\n",
    "    \"Volume MYO pred\",\n",
    "    \"Volume MYO GT\",\n",
    "    \"Err MYO(ml)\",\n",
    "]\n",
    "\n",
    "\n",
    "#\n",
    "# Utils functions used to sort strings into a natural order\n",
    "#\n",
    "def conv_int(i):\n",
    "    return int(i) if i.isdigit() else i\n",
    "\n",
    "\n",
    "def natural_order(sord):\n",
    "    \"\"\"\n",
    "    Sort a (list,tuple) of strings into natural order.\n",
    "\n",
    "    Ex:\n",
    "\n",
    "    ['1','10','2'] -> ['1','2','10']\n",
    "\n",
    "    ['abc1def','ab10d','b2c','ab1d'] -> ['ab1d','ab10d', 'abc1def', 'b2c']\n",
    "\n",
    "    \"\"\"\n",
    "    if isinstance(sord, tuple):\n",
    "        sord = sord[0]\n",
    "    return [conv_int(c) for c in re.split(r\"(\\d+)\", sord)]\n",
    "\n",
    "\n",
    "#\n",
    "# Utils function to load and save nifti files with the nibabel package\n",
    "#\n",
    "\n",
    "img_path = \"ACDC\\database\"\n",
    "\n",
    "\n",
    "def load_nii(img_path):\n",
    "    \"\"\"\n",
    "    Function to load a 'nii' or 'nii.gz' file, The function returns\n",
    "    everyting needed to save another 'nii' or 'nii.gz'\n",
    "    in the same dimensional space, i.e. the affine matrix and the header\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "\n",
    "    img_path: string\n",
    "    String with the path of the 'nii' or 'nii.gz' image file name.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Three element, the first is a numpy array of the image values,\n",
    "    the second is the affine transformation of the image, and the\n",
    "    last one is the header of the image.\n",
    "    \"\"\"\n",
    "    nimg = nib.load(img_path)\n",
    "    return nimg.get_fdata(), nimg.affine, nimg.header\n",
    "\n",
    "\n",
    "def save_nii(img_path, data, affine, header):\n",
    "    \"\"\"\n",
    "    Function to save a 'nii' or 'nii.gz' file.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "\n",
    "    img_path: string\n",
    "    Path to save the image should be ending with '.nii' or '.nii.gz'.\n",
    "\n",
    "    data: np.array\n",
    "    Numpy array of the image data.\n",
    "\n",
    "    affine: list of list or np.array\n",
    "    The affine transformation to save with the image.\n",
    "\n",
    "    header: nib.Nifti1Header\n",
    "    The header that define everything about the data\n",
    "    (pleasecheck nibabel documentation).\n",
    "    \"\"\"\n",
    "    nimg = nib.Nifti1Image(data, affine=affine, header=header)\n",
    "    nimg.to_filename(img_path)\n",
    "\n",
    "\n",
    "#\n",
    "# Functions to process files, directories and metrics\n",
    "#\n",
    "def metrics(img_gt, img_pred, voxel_size):\n",
    "    \"\"\"\n",
    "    Function to compute the metrics between two segmentation maps given as input.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    img_gt: np.array\n",
    "    Array of the ground truth segmentation map.\n",
    "\n",
    "    img_pred: np.array\n",
    "    Array of the predicted segmentation map.\n",
    "\n",
    "    voxel_size: list, tuple or np.array\n",
    "    The size of a voxel of the images used to compute the volumes.\n",
    "\n",
    "    Return\n",
    "    ------\n",
    "    A list of metrics in this order, [Dice LV, Volume LV, Volume GT, Err LV(ml),\n",
    "    Dice RV, Volume RV, Volume GT, Err RV(ml), Dice MYO, Volume MYO, Volume GT, Err MYO(ml)]\n",
    "    \"\"\"\n",
    "\n",
    "    if img_gt.ndim != img_pred.ndim:\n",
    "        raise ValueError(\n",
    "            \"The arrays 'img_gt' and 'img_pred' should have the \"\n",
    "            \"same dimension, {} against {}\".format(img_gt.ndim, img_pred.ndim)\n",
    "        )\n",
    "\n",
    "    res = []\n",
    "    # Loop on each classes of the input images\n",
    "    for c in [3, 1, 2]:\n",
    "        # Copy the gt image to not alterate the input\n",
    "        gt_c_i = np.copy(img_gt)\n",
    "        gt_c_i[gt_c_i != c] = 0\n",
    "\n",
    "        # Copy the pred image to not alterate the input\n",
    "        pred_c_i = np.copy(img_pred)\n",
    "        pred_c_i[pred_c_i != c] = 0\n",
    "\n",
    "        # Clip the value to compute the volumes\n",
    "        gt_c_i = np.clip(gt_c_i, 0, 1)\n",
    "        pred_c_i = np.clip(pred_c_i, 0, 1)\n",
    "\n",
    "        # Compute the Dice\n",
    "        # dice = dc(gt_c_i, pred_c_i)\n",
    "        dice = 1\n",
    "\n",
    "        # Eventueel alternatief\n",
    "        gt_volume = np.sum(gt_c_i)\n",
    "        pred_volume = np.sum(pred_c_i)\n",
    "        intersect = np.sum(gt_c_i * pred_c_i)\n",
    "        dice = (2 * intersect) / (gt_volume + pred_volume)\n",
    "\n",
    "        # Compute volume\n",
    "        volpred = pred_c_i.sum() * np.prod(voxel_size) / 1000.0\n",
    "        volgt = gt_c_i.sum() * np.prod(voxel_size) / 1000.0\n",
    "\n",
    "        res += [dice, volpred, volgt, volpred - volgt]\n",
    "\n",
    "    return res\n",
    "\n",
    "\n",
    "def compute_metrics_on_files(path_gt, path_pred):\n",
    "    \"\"\"\n",
    "    Function to give the metrics for two files\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "\n",
    "    path_gt: string\n",
    "    Path of the ground truth image.\n",
    "\n",
    "    path_pred: string\n",
    "    Path of the predicted image.\n",
    "    \"\"\"\n",
    "    gt, _, header = load_nii(path_gt)\n",
    "    pred, _, _ = load_nii(path_pred)\n",
    "    zooms = header.get_zooms()\n",
    "\n",
    "    name = os.path.basename(path_gt)\n",
    "    name = name.split(\".\")[0]\n",
    "    res = metrics(gt, pred, zooms)\n",
    "    res = [\"{:.3f}\".format(r) for r in res]\n",
    "\n",
    "    formatting = \"{:<20}\" + \"{:>12}\" * len(res)\n",
    "    output = formatting.format(name, *res)\n",
    "\n",
    "    print(formatting.format(*HEADER))\n",
    "    print(output)\n",
    "\n",
    "    # formatting = \"{:>14}, {:>7}, {:>9}, {:>10}, {:>7}, {:>9}, {:>10}, {:>8}, {:>10}, {:>11}\"\n",
    "    # print(formatting.format(*HEADER))\n",
    "    # print(formatting.format(name, *res))\n",
    "\n",
    "    # return [name, *res]\n",
    "    return res\n",
    "\n",
    "\n",
    "def compute_metrics_on_directories(dir_gt, dir_pred):\n",
    "    \"\"\"\n",
    "    Function to generate a csv file for each images of two directories.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "\n",
    "    path_gt: string\n",
    "    Directory of the ground truth segmentation maps.\n",
    "\n",
    "    path_pred: string\n",
    "    Directory of the predicted segmentation maps.\n",
    "    \"\"\"\n",
    "    lst_gt = sorted(glob(os.path.join(dir_gt, \"*\")), key=natural_order)\n",
    "    lst_pred = sorted(glob(os.path.join(dir_pred, \"*\")), key=natural_order)\n",
    "\n",
    "    res = []\n",
    "    for p_gt, p_pred in zip(lst_gt, lst_pred):\n",
    "        if os.path.basename(p_gt) != os.path.basename(p_pred):\n",
    "            raise ValueError(\n",
    "                \"The two files don't have the same name\"\n",
    "                \" {}, {}.\".format(os.path.basename(p_gt), os.path.basename(p_pred))\n",
    "            )\n",
    "\n",
    "        gt, _, header = load_nii(p_gt)\n",
    "        pred, _, _ = load_nii(p_pred)\n",
    "        zooms = header.get_zooms()\n",
    "        res.append(metrics(gt, pred, zooms))\n",
    "\n",
    "    lst_name_gt = [os.path.basename(gt).split(\".\")[0] for gt in lst_gt]\n",
    "    res = [\n",
    "        [\n",
    "            n,\n",
    "        ]\n",
    "        + r\n",
    "        for r, n in zip(res, lst_name_gt)\n",
    "    ]\n",
    "    df = pd.DataFrame(res, columns=HEADER)\n",
    "    df.to_csv(\"results_{}.csv\".format(time.strftime(\"%Y%m%d_%H%M%S\")), index=False)\n",
    "\n",
    "\n",
    "def main(path_gt, path_pred):\n",
    "    \"\"\"\n",
    "    Main function to select which method to apply on the input parameters.\n",
    "    \"\"\"\n",
    "    if os.path.isfile(path_gt) and os.path.isfile(path_pred):\n",
    "        compute_metrics_on_files(path_gt, path_pred)\n",
    "    elif os.path.isdir(path_gt) and os.path.isdir(path_pred):\n",
    "        compute_metrics_on_directories(path_gt, path_pred)\n",
    "    else:\n",
    "        raise ValueError(\"The paths given needs to be two directories or two files.\")\n",
    "\n",
    "\n",
    "# if __name__ == \"__main__\":\n",
    "#     parser = argparse.ArgumentParser(\n",
    "#         description=\"Script to compute ACDC challenge metrics.\")\n",
    "#     parser.add_argument(\"GT_IMG\", type=str, help=\"Ground Truth image\")\n",
    "#     parser.add_argument(\"PRED_IMG\", type=str, help=\"Predicted image\")\n",
    "#     args = parser.parse_args()\n",
    "#     main(args.GT_IMG, args.PRED_IMG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "142caccb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Testing\n",
    "path_gt = os.path.join('ACDC','database','testing','patient101','patient101_frame01.nii.gz')\n",
    "dir_gt = os.path.join('ACDC','database','testing','patient101','patient101_frame01_gt.nii.gz')\n",
    "\n",
    "path_pred = os.path.join('ACDC','database','testing','patient101','patient101_frame01.nii.gz')\n",
    "dir_pred = os.path.join('ACDC','database','testing','patient101','patient101_frame01_ml_pred.nii.gz')\n",
    "\n",
    "test_files = compute_metrics_on_files(dir_gt, dir_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0451a847",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = \"ACDC/database\"\n",
    "test_dict_gt = build_dict_images(data_path, mode='testing').ravel()\n",
    "test_dict_pred = build_dict_images_pred(data_path, mode='testing')\n",
    "\n",
    "# Ground truth\n",
    "# First file\n",
    "print(len(test_dict_gt))\n",
    "image_paths_gt_ED = [d['image'] for d in test_dict_gt if d[\"first_file\"]]\n",
    "label_paths_gt_ED = [d['label'] for d in test_dict_gt if d[\"first_file\"]]\n",
    "\n",
    "# Last file\n",
    "image_paths_gt_ES = [d['image'] for d in test_dict_gt if not d[\"first_file\"]]\n",
    "label_paths_gt_ES = [d['label'] for d in test_dict_gt if not d[\"first_file\"]]\n",
    "\n",
    "print('Image path gt ED is ', image_paths_gt_ED) \n",
    "print('Label path gt ED is ', label_paths_gt_ED)\n",
    "print('Image path gt ES is ', image_paths_gt_ES)\n",
    "print('Label path gt ES is ', label_paths_gt_ES)\n",
    "\n",
    "# Prediction\n",
    "image_paths_pred_ED = [d['image'] for d in test_dict_pred if d[\"first_file\"]]\n",
    "label_paths_pred_ED = [d['label'] for d in test_dict_pred if d[\"first_file\"]]\n",
    "\n",
    "image_paths_pred_ES = [d['image'] for d in test_dict_pred if not d[\"first_file\"]]\n",
    "label_paths_pred_ES = [d['label'] for d in test_dict_pred if not d[\"first_file\"]]\n",
    "\n",
    "print('Image path prediction ED is ', image_paths_pred_ED)\n",
    "print('Label path prediction ED is ', label_paths_pred_ED)\n",
    "print('Image path prediction ES is ', image_paths_pred_ES)\n",
    "print('Label path prediction ES is ', label_paths_pred_ES)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c7be742",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_files_ED = [compute_metrics_on_files(label_paths_gt_ED[i], label_paths_pred_ED[i]) for i in range(len(label_paths_gt_ED))]\n",
    "test_files_ES = [compute_metrics_on_files(label_paths_gt_ES[i], label_paths_pred_ES[i]) for i in range(len(label_paths_gt_ES))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "049a0262",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ED\n",
    "LV_GT_ED = np.array([float(metric[2]) for metric in test_files_ED])\n",
    "LV_pred_ED = np.array([float(metric[1]) for metric in test_files_ED])\n",
    "LV_dice_ED = np.array([float(metric[0]) for metric in test_files_ED])\n",
    "LV_err_ED = np.array([float(metric[3]) for metric in test_files_ED])\n",
    "\n",
    "RV_GT_ED = np.array([float(metric[6]) for metric in test_files_ED])\n",
    "RV_pred_ED = np.array([float(metric[5]) for metric in test_files_ED])\n",
    "RV_dice_ED = np.array([float(metric[4]) for metric in test_files_ED])\n",
    "RV_err_ED = np.array([float(metric[7]) for metric in test_files_ED])\n",
    "\n",
    "MYO_GT_ED = np.array([float(metric[10]) for metric in test_files_ED])\n",
    "MYO_pred_ED = np.array([float(metric[9]) for metric in test_files_ED])\n",
    "MYO_dice_ED = np.array([float(metric[8]) for metric in test_files_ED])\n",
    "MYO_err_ED = np.array([float(metric[11]) for metric in test_files_ED])\n",
    "\n",
    "# ES\n",
    "LV_GT_ES = np.array([float(metric[2]) for metric in test_files_ES])\n",
    "LV_pred_ES = np.array([float(metric[1]) for metric in test_files_ES])\n",
    "LV_dice_ES = np.array([float(metric[0]) for metric in test_files_ES])\n",
    "LV_err_ES = np.array([float(metric[3]) for metric in test_files_ES])\n",
    "\n",
    "RV_GT_ES = np.array([float(metric[6]) for metric in test_files_ES])\n",
    "RV_pred_ES = np.array([float(metric[5]) for metric in test_files_ES])\n",
    "RV_dice_ES = np.array([float(metric[4]) for metric in test_files_ES])\n",
    "RV_err_ES = np.array([float(metric[7]) for metric in test_files_ES])\n",
    "\n",
    "MYO_GT_ES = np.array([float(metric[10]) for metric in test_files_ES])\n",
    "MYO_pred_ES = np.array([float(metric[9]) for metric in test_files_ES])\n",
    "MYO_dice_ES = np.array([float(metric[8]) for metric in test_files_ES])\n",
    "MYO_err_ES = np.array([float(metric[11]) for metric in test_files_ES])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7c36b2e0",
   "metadata": {},
   "source": [
    "Other ways for visualization\n",
    "> Comparison graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fd62eb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_comparison_graph(gt, pred, dice, err, color_boundaries=20, part=\"LV\", mode=\"ED\"):\n",
    "    plt.clf()\n",
    "    vmin = min(err)\n",
    "    vmax = max(err)\n",
    "\n",
    "    colors = ['red', 'yellow', 'green', 'yellow', 'red']\n",
    "    color_positions = [vmin, -1 * color_boundaries, 0, color_boundaries, vmax]\n",
    "    norm = plt.Normalize(vmin, vmax)\n",
    "    print(norm(color_positions))\n",
    "\n",
    "    colormap = list(zip(norm(color_positions), colors))\n",
    "    cmap = LinearSegmentedColormap.from_list(\"custom_cmap\", colormap)\n",
    "\n",
    "    plt.scatter(gt, pred, c=err, cmap=cmap, vmin=vmin, vmax=vmax, s=dice * 100)\n",
    "    plt.plot([min(gt), max(gt)], [min(gt), max(gt)], 'k--')\n",
    "    plt.xlabel(f'Ground Truth {part} Volume [ml]')\n",
    "    plt.ylabel(f'Predicted {part} Volume [ml]')\n",
    "    plt.title(f'Ground Truth vs. Predicted Volume of {part} {mode}')\n",
    "\n",
    "    cbar = plt.colorbar()\n",
    "    cbar.set_label('Error Value')\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1541570e",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_comparison_graph(LV_GT_ED, LV_pred_ED, LV_dice_ED, LV_err_ED)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93c71590",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_comparison_graph(RV_GT_ED, RV_pred_ED, RV_dice_ED, RV_err_ED, color_boundaries=18, part=\"RV\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa1aa317",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_comparison_graph(MYO_GT_ED, MYO_pred_ED, MYO_dice_ED, MYO_err_ED, color_boundaries=18, part=\"MYO\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "5b91e382",
   "metadata": {},
   "source": [
    "> Correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a21ffc52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlation coefficients\n",
    "def calculate_correlation(gt, pred, err):\n",
    "    gt_mean = statistics.mean(gt)\n",
    "    gt_std = statistics.stdev(gt)\n",
    "    pred_mean = statistics.mean(pred)\n",
    "    pred_std = statistics.stdev(pred)\n",
    "    covariance = sum((x - gt_mean) * (y - pred_mean) for x, y in zip(gt, pred)) / len(gt)\n",
    "    corr = covariance / (gt_std * pred_std)\n",
    "    bias = statistics.mean(err)\n",
    "    err_std = statistics.stdev(err)\n",
    "    loa = 1.96*err_std\n",
    "    return corr, bias, loa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4ecf214",
   "metadata": {},
   "outputs": [],
   "source": [
    "LV_corr, LV_bias, LV_LOA = calculate_correlation(LV_GT_ED, LV_pred_ED, LV_err_ED)\n",
    "print('The correlation coefficient for left ventricle EDV is ', LV_corr)\n",
    "print('The bias for the left ventricle EDV is ', LV_bias)\n",
    "print('The LOA of the left ventricle EDV is ', LV_LOA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12c061d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "RV_corr, RV_bias, RV_LOA = calculate_correlation(RV_GT_ED, RV_pred_ED, RV_err_ED)\n",
    "print('The correlation coefficient for right ventricle EDV is ', RV_corr)\n",
    "print('The bias for the right ventricle EDV is ', RV_bias)\n",
    "print('The LOA of the right ventricle EDV is ', RV_LOA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccb53139",
   "metadata": {},
   "outputs": [],
   "source": [
    "MYO_corr, MYO_bias, MYO_LOA = calculate_correlation(MYO_GT_ED, MYO_pred_ED, MYO_err_ED)\n",
    "print('The correlation coefficient for myocardium EDV is ', MYO_corr)\n",
    "print('The bias for the left myocardium EDV is ', MYO_bias)\n",
    "print('The LOA of the left myocardium EDV is ', MYO_LOA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d85b01dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Average dice\n",
    "\n",
    "av_LV_dice = statistics.mean(LV_dice_ED)\n",
    "av_RV_dice = statistics.mean(RV_dice_ED)\n",
    "av_MYO_dice = statistics.mean(MYO_dice_ED)\n",
    "print('Average dice left ventricle is ', av_LV_dice)\n",
    "print('Average dice right ventricle is ', av_RV_dice)\n",
    "print('Average dice myocardium is ', av_MYO_dice)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ce64f7e5",
   "metadata": {},
   "source": [
    "ES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "437c2b32",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_comparison_graph(LV_GT_ES, LV_pred_ES, LV_dice_ES, LV_err_ES, color_boundaries=15, part=\"LV\", mode=\"ES\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ff66a64",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_comparison_graph(RV_GT_ES, RV_pred_ES, RV_dice_ES, RV_err_ES, color_boundaries=10, part=\"RV\", mode=\"ES\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "931d0abb",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_comparison_graph(MYO_GT_ES, MYO_pred_ES, MYO_dice_ES, MYO_err_ES, color_boundaries=30, part=\"MYO\", mode=\"ES\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d30933c",
   "metadata": {},
   "outputs": [],
   "source": [
    "LV_corr, LV_bias, LV_LOA = calculate_correlation(LV_GT_ES, LV_pred_ES, LV_err_ES)\n",
    "print('The correlation coefficient for left ventricle ESV is ', LV_corr)\n",
    "print('The bias for the left ventricle ESV is ', LV_bias)\n",
    "print('The LOA of the left ventricle ESV is ', LV_LOA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "823112fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "RV_corr, RV_bias, RV_LOA = calculate_correlation(RV_GT_ES, RV_pred_ES, RV_err_ES)\n",
    "print('The correlation coefficient for right ventricle ESV is ', RV_corr)\n",
    "print('The bias for the right ventricle ESV is ', RV_bias)\n",
    "print('The LOA of the right ventricle ESV is ', RV_LOA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50408bf1",
   "metadata": {},
   "outputs": [],
   "source": [
    "MYO_corr, MYO_bias, MYO_LOA = calculate_correlation(MYO_GT_ES, MYO_pred_ES, MYO_err_ES)\n",
    "print('The correlation coefficient for myocardium ESV is ', MYO_corr)\n",
    "print('The bias for the left myocardium ESV is ', MYO_bias)\n",
    "print('The LOA of the left myocardium ESV is ', MYO_LOA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf0f0d91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Average dice\n",
    "av_LV_dice = statistics.mean(LV_dice_ES)\n",
    "av_RV_dice = statistics.mean(RV_dice_ES)\n",
    "av_MYO_dice = statistics.mean(MYO_dice_ES)\n",
    "print('Average dice left ventricle ESV is ', av_LV_dice)\n",
    "print('Average dice right ventricle ESV is ', av_RV_dice)\n",
    "print('Average dice myocardium ESV is ', av_MYO_dice)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "ecf8346f",
   "metadata": {},
   "source": [
    "EJ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1671753",
   "metadata": {},
   "outputs": [],
   "source": [
    "EJ_LV_GT = [((LV_GT_ED[i] - LV_GT_ES[i])/ LV_GT_ED[i]) * 100 for i in range(len(LV_GT_ED))]\n",
    "EJ_RV_GT = [((RV_GT_ED[i] - RV_GT_ES[i])/ RV_GT_ED[i]) * 100 for i in range(len(RV_GT_ED))]\n",
    "\n",
    "EJ_LV_pred = [((LV_pred_ED[i] - LV_pred_ES[i])/ (LV_pred_ED[i]+0.1)) * 100 for i in range(len(LV_pred_ED))]\n",
    "EJ_RV_pred = [((RV_pred_ED[i] - RV_pred_ES[i])/ (RV_pred_ED[i]+0.1)) * 100 for i in range(len(RV_pred_ED))]\n",
    "\n",
    "EJ_LV_err = [EJ_LV_GT[i] - EJ_LV_pred[i] for i in range(len(EJ_LV_GT))]\n",
    "EJ_RV_err = [EJ_RV_GT[i] - EJ_RV_pred[i] for i in range(len(EJ_RV_GT))]\n",
    "\n",
    "EJ_LV_corr, EJ_LV_bias, LV_LOA = calculate_correlation(EJ_LV_GT, EJ_LV_pred, EJ_LV_err)\n",
    "print('The correlation of the ejection fraction for left ventricle is ', EJ_LV_corr)\n",
    "print('The mean bias of the ejection fraction of left ventricle is ', EJ_LV_bias)\n",
    "print('The LOA of the ejection fraction of the left ventricle is ', LV_LOA)\n",
    "\n",
    "EJ_RV_corr, EJ_RV_bias, RV_LOA = calculate_correlation(EJ_RV_GT, EJ_RV_pred, EJ_RV_err)\n",
    "print('The correlation of the ejection fraction for right ventricle is ', EJ_RV_corr)\n",
    "print('The mean bias of the ejection fraction of right ventricle is ', EJ_RV_bias)\n",
    "print('The LOA of the ejection fraction of the right ventricle is ', RV_LOA)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "de9c181c",
   "metadata": {},
   "source": [
    "Best and worst results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9a463de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This determines the best and worst performing predictions based on the dice score of all the segments and both time frames\n",
    "# Uses the LV_dice_ES list to find back the correct patient, but could be any of the lists\n",
    "\n",
    "\n",
    "score_sum = {LV_dice_ES[x]: sum([LV_dice_ES[x], RV_dice_ES[x], MYO_dice_ES[x], LV_dice_ED[x], RV_dice_ED[x], MYO_dice_ED[x]]) for x in range(len(LV_dice_ES))}\n",
    "max_value = max(score_sum, key=score_sum.get)\n",
    "min_value = min(score_sum, key=score_sum.get)\n",
    "print(min_value)\n",
    "del score_sum[min_value]\n",
    "\n",
    "min_value = min(score_sum, key=score_sum.get)\n",
    "print(min_value)\n",
    "\n",
    "max_index = np.where(LV_dice_ES == max_value)[0][0]\n",
    "min_index = np.where(LV_dice_ES == min_value)[0][0]\n",
    "\n",
    "print(\"Max value:\", max_value)\n",
    "print(\"Max index:\", max_index)\n",
    "\n",
    "print(\"Min value:\", min_value)\n",
    "print(\"Min index:\", min_index)\n",
    "\n",
    "# Show slices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ef0ed36",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Best image path gt ES is ', image_paths_gt_ES[max_index])\n",
    "print('Best label path gt ES is ', label_paths_gt_ES[max_index])\n",
    "\n",
    "print('Best image path prediction ES is ', image_paths_pred_ES[max_index])\n",
    "print('Best label path prediction ES is ', label_paths_pred_ES[max_index])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35b72c7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_slice(index, image_gt, label_gt, image_pred, label_pred,  channel_index=0, frame='ES'):\n",
    "    image_path_gt = image_gt[index]\n",
    "    label_path_gt = label_gt[index]\n",
    "    image_path_pred = image_pred[index]\n",
    "    label_path_pred = label_pred[index]\n",
    "\n",
    "    # Load image and label data\n",
    "    image_gt = nib.load(image_path_gt).get_fdata()\n",
    "    label_gt_ES = nib.load(label_path_gt).get_fdata()\n",
    "    image_pred = nib.load(image_path_pred).get_fdata()\n",
    "    label_pred = nib.load(label_path_pred).get_fdata()\n",
    "    return (\n",
    "        [\n",
    "        image_gt[..., channel_index],\n",
    "        label_gt_ES[..., channel_index],\n",
    "        image_pred[..., channel_index],\n",
    "        label_pred[..., channel_index]\n",
    "        ]\n",
    "    )\n",
    "def display_slices(min_index, max_index,  image_gt, label_gt, image_pred, label_pred,  channel_index=0, frame='ES'):\n",
    "    best_image_gt_ES_channel, best_label_gt_ES_channel, best_image_pred_ES_channel, best_label_pred_ES_channel = get_slice(max_index, image_gt, label_gt, image_pred, label_pred, channel_index)\n",
    "    worst_image_gt_ES_channel, worst_label_gt_ES_channel, worst_image_pred_ES_channel, worst_label_pred_ES_channel = get_slice(min_index, image_gt, label_gt, image_pred, label_pred, channel_index)\n",
    "\n",
    "    plt.subplot(2, 2, 1)\n",
    "    plt.imshow(best_image_gt_ES_channel, cmap='gray')\n",
    "    plt.title(f\"Ground truth best prediction {frame}\")\n",
    "    plt.imshow(best_label_gt_ES_channel, alpha=0.5, cmap='Reds')\n",
    "\n",
    "\n",
    "    plt.subplot(2, 2, 2)\n",
    "    plt.imshow(best_image_pred_ES_channel, cmap='gray')\n",
    "    plt.title(f\"Best prediction {frame}\")\n",
    "    plt.imshow(best_label_pred_ES_channel, cmap='Reds', alpha=0.5)\n",
    "\n",
    "    plt.subplot(2, 2, 3)\n",
    "    plt.imshow(worst_image_gt_ES_channel, cmap='gray')\n",
    "    plt.title(f\"Ground truth worst prediction {frame}\")\n",
    "    plt.imshow(worst_label_gt_ES_channel, alpha=0.5, cmap='Reds')\n",
    "\n",
    "    plt.subplot(2, 2, 4)\n",
    "    plt.imshow(worst_image_pred_ES_channel, cmap='gray')\n",
    "    plt.title(f\"Worst prediction {frame}\")\n",
    "    plt.imshow(worst_label_pred_ES_channel, cmap='Reds', alpha=0.5)\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8074a3d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "display_slices(min_index, max_index, image_paths_gt_ES, label_paths_gt_ES, image_paths_pred_ES, label_paths_pred_ES, channel_index=4, frame='ES')\n",
    "compute_metrics_on_files(label_paths_gt_ED[max_index], label_paths_pred_ED[max_index])\n",
    "compute_metrics_on_files(label_paths_gt_ES[max_index], label_paths_pred_ES[max_index])\n",
    "print()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2d5604a",
   "metadata": {},
   "outputs": [],
   "source": [
    "display_slices(min_index, max_index, image_paths_gt_ED, label_paths_gt_ED, image_paths_pred_ED, label_paths_pred_ED, channel_index=4, frame='ED')\n",
    "compute_metrics_on_files(label_paths_gt_ED[min_index], label_paths_pred_ED[min_index])\n",
    "compute_metrics_on_files(label_paths_gt_ES[min_index], label_paths_pred_ES[min_index])\n",
    "print()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "c4ee1d4a",
   "metadata": {},
   "source": [
    "Hausdorff distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a633891f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First file\n",
    "label_paths_gt_ED = [d['label'] for d in test_dict_gt if d[\"first_file\"]]\n",
    "print(len(label_paths_gt_ED))\n",
    "\n",
    "# Last file\n",
    "label_paths_gt_ES = [d['label'] for d in test_dict_gt if not d[\"first_file\"]]\n",
    "\n",
    "# Prediction\n",
    "label_paths_pred_ED = [d['label'] for d in test_dict_pred if d[\"first_file\"]]\n",
    "\n",
    "label_paths_pred_ES = [d['label'] for d in test_dict_pred if not d[\"first_file\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3a8a4e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# GT and pred of ED\n",
    "dir_gt_ED = []\n",
    "dir_pred_ED = []\n",
    "for label_path_gt_ED, label_path_pred_ED in zip(label_paths_gt_ED, label_paths_pred_ED):\n",
    "    label_GT_ED = nib.load(label_path_gt_ED).get_fdata()\n",
    "    label_pred_ED = nib.load(label_path_pred_ED).get_fdata()\n",
    "    dir_gt_ED.append(label_GT_ED)\n",
    "    dir_pred_ED.append(label_pred_ED)\n",
    "\n",
    "# GT and pred of ES\n",
    "dir_gt_ES = []\n",
    "dir_pred_ES = []\n",
    "for label_path_gt, label_path_pred in zip(label_paths_gt_ES, label_paths_pred_ES):\n",
    "    label_GT_ES = nib.load(label_path_gt).get_fdata()\n",
    "    label_pred = nib.load(label_path_pred).get_fdata()\n",
    "    dir_gt_ES.append(label_GT_ES)\n",
    "    dir_pred_ES.append(label_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "106f8cb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "label_gt_ED = []\n",
    "label_gt_ED_LV = []\n",
    "label_gt_ED_wall = []\n",
    "label_gt_ED_RV = []\n",
    "\n",
    "for i in range(50):\n",
    "    label_gt_ED_single = np.array(dir_gt_ED[i], dtype=np.float32)\n",
    "    label_gt_ED.append(label_gt_ED_single)\n",
    "    label_gt_ED_LV.append(np.where(label_gt_ED_single == 1,1,0))\n",
    "    label_gt_ED_wall.append(np.where(label_gt_ED_single == 2,1,0))\n",
    "    label_gt_ED_RV.append(np.where(label_gt_ED_single == 3,1,0))\n",
    "\n",
    "\n",
    "label_pred_ED = []\n",
    "label_pred_ED_LV = []\n",
    "label_pred_ED_wall = []\n",
    "label_pred_ED_RV = []\n",
    "\n",
    "for i in range(50):\n",
    "    label_pred_ED_single = np.array(dir_pred_ED[i], dtype=np.float32)\n",
    "    label_pred_ED.append(label_pred_ED_single)\n",
    "    label_pred_ED_LV.append(np.where(label_pred_ED_single == 1,1,0))\n",
    "    label_pred_ED_wall.append(np.where(label_pred_ED_single == 2,1,0))\n",
    "    label_pred_ED_RV.append(np.where(label_pred_ED_single == 3,1,0))\n",
    "\n",
    "label_gt_ES = []\n",
    "label_gt_ES_LV = []\n",
    "label_gt_ES_wall = []\n",
    "label_gt_ES_RV = []\n",
    "\n",
    "for i in range(50):\n",
    "    label_gt_ES_single = np.array(dir_gt_ES[i], dtype=np.float32)\n",
    "    label_gt_ES.append(label_gt_ES_single)\n",
    "    label_gt_ES_LV.append(np.where(label_gt_ES_single == 1,1,0))\n",
    "    label_gt_ES_wall.append(np.where(label_gt_ES_single == 2,1,0))\n",
    "    label_gt_ES_RV.append(np.where(label_gt_ES_single == 3,1,0))\n",
    "\n",
    "\n",
    "label_pred = []\n",
    "label_pred_ES_LV = []\n",
    "label_pred_ES_wall = []\n",
    "label_pred_ES_RV = []\n",
    "\n",
    "for i in range(50):\n",
    "    label_pred_ES_single = np.array(dir_pred_ES[i], dtype=np.float32)\n",
    "    label_pred.append(label_pred_ES_single)\n",
    "    label_pred_ES_LV.append(np.where(label_pred_ES_single == 1,1,0))\n",
    "    label_pred_ES_wall.append(np.where(label_pred_ES_single == 2,1,0))\n",
    "    label_pred_ES_RV.append(np.where(label_pred_ES_single == 3,1,0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "385e0f03",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(label_pred_ES_RV))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a682e84",
   "metadata": {},
   "outputs": [],
   "source": [
    "hd_LV_ED = []\n",
    "hd_wall_ED = []\n",
    "hd_RV_ED = []\n",
    "hd_LV_ES = []\n",
    "hd_wall_ES = []\n",
    "hd_RV_ES = []\n",
    "\n",
    "for i in range(50):\n",
    "    # Calculate Hausdorff distance for ED\n",
    "    hd_lv_ed = max([hausdorff_distance(label_pred_ED_LV[i][:,:,y], label_gt_ED_LV[i][:,:,y], distance='euclidean') for y in range(label_pred_ED_LV[i].shape[2])])\n",
    "    hd_wall_ed = max([hausdorff_distance(label_pred_ED_wall[i][:,:,y], label_gt_ED_wall[i][:,:,y], distance='euclidean') for y in range(label_pred_ED_wall[i].shape[2])])\n",
    "    hd_rv_ed = max([hausdorff_distance(label_pred_ED_RV[i][:,:,y], label_gt_ED_RV[i][:,:,y], distance='euclidean') for y in range(label_pred_ED_RV[i].shape[2])])\n",
    "    hd_LV_ED.append(hd_lv_ed)\n",
    "    hd_wall_ED.append(hd_wall_ed)\n",
    "    hd_RV_ED.append(hd_rv_ed)\n",