Calculate_SSIM_PSNR.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "11ee4019-f711-41cd-84a2-2e1a0950f390",
   "metadata": {},
   "source": [
    "# Calculate SSIM/PSNR values for ROI"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e98e9052-196d-42f8-901c-2f9cb9cedef5",
   "metadata": {},
   "source": [
    "## Generate Masks\n",
    "Clone https://github.com/jasonccai/HeadCTSegmentation and download the model weights. \n",
    "Copy z_unet.py and weights.hdf5 into this folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81b93582-374b-4d79-9c00-ffd49afdfa46",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"CUDA_DEVICE_ORDER\"]=\"PCI_BUS_ID\"\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"2\"\n",
    "import sys\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from skimage.metrics import structural_similarity as ssim\n",
    "from skimage.metrics import peak_signal_noise_ratio as psnr\n",
    "import pandas as pd\n",
    "from scipy.ndimage import binary_fill_holes\n",
    "from scipy.stats import shapiro\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f188b613-b1b4-4b04-83c5-92495363bd68",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load segmentation model\n",
    "from z_unet import unet\n",
    "\n",
    "U_Net = unet(nb_classes=17, savefolder=\"./\", flag=\"a\")\n",
    "U_Net.load_weights(\"./weights.hdf5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7013e845-59c2-4490-96a1-0f9c2b96a374",
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(\"./Masks/\"):\n",
    "    os.mkdir(\"./Masks/\")\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae5ab885-6172-4ed5-a657-939d4f167b75",
   "metadata": {
    "tags": []
   },
   "source": [
    "### RSNA test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a45ff16a-bec9-42f6-9f67-59108ce2c34d",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_list = pd.read_csv(\"./test_data.csv\")[\"filename\"].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea11b288-580c-4a7e-9ee7-e006243b1511",
   "metadata": {},
   "outputs": [],
   "source": [
    "#test segmentation model\n",
    "for filename in test_list[0:3]:\n",
    "    filename = filename.replace(\".dcm\", \".npy\")\n",
    "    img = np.load(\"./Data/4096/\" + filename)*4095 - 1024\n",
    "\n",
    "    pred = U_Net.predict(img[None, :, :, None])\n",
    "    pred_re = pred.reshape(512, 512, 17)\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.imshow(img.squeeze(), cmap='gray', vmin=0, vmax=80)\n",
    "    plt.imshow(1-pred_re[:, :, 0], alpha=0.4, cmap='Reds', vmin=0, vmax=1)\n",
    "    segmentation = img.copy()\n",
    "    segmentation[pred_re[:, :, 0]>0.3] = 0 \n",
    "    segmentation[pred_re[:, :, 0]<=0.3] = 1\n",
    "    mask = binary_fill_holes(segmentation)\n",
    "\n",
    "    plt.figure()\n",
    "    plt.imshow(mask, cmap='gray')\n",
    "    #np.save(f\"./Masks/mask_{filename}\", segmentation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ff437a5-92e3-4d20-9b3e-0e7f86b45de4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#loop trough entire testset and save ROI masks for each file \n",
    "for i, filename in enumerate(test_list):\n",
    "    sys.stdout.write('\\r'+f'{i}/{len(test_list)}')    \n",
    "    filename = filename.replace(\".dcm\", \".npy\")\n",
    "    img = np.load(\"./Data/4096/\" + filename) - 1024\n",
    "\n",
    "    pred = U_Net.predict(img[None, :, :, None], verbose=False)\n",
    "    pred_re = pred.reshape(512, 512, 17)\n",
    "    \n",
    "    segmentation = img.copy()\n",
    "    segmentation[pred_re[:, :, 0]>0.3] = 0 \n",
    "    segmentation[pred_re[:, :, 0]<=0.3] = 1\n",
    "    \n",
    "    np.save(f\"./Masks/mask_{filename}\", segmentation)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "219f77e8-2356-4756-8aa8-c608022ed393",
   "metadata": {},
   "source": [
    "### Repeat for CQ500 dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a992334e-d641-4b73-bc14-f4f09926b03d",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_list = os.listdir(\"/data-pool/data_no_backup/ga63cun/CQ500/4096/\")\n",
    "#test segmentation model\n",
    "for filename in test_list[1440:1443]:\n",
    "    filename = filename.replace(\".dcm\", \".npy\")\n",
    "    print(filename)\n",
    "    img = np.load(\"./Data/CQ500/4096/\" + filename)*4095 - 1024\n",
    "\n",
    "    pred = U_Net.predict(img[None, :, :, None])\n",
    "    pred_re = pred.reshape(512, 512, 17)\n",
    "    \n",
    "    plt.figure()\n",
    "    plt.imshow(img.squeeze(), cmap='gray', vmin=0, vmax=80)\n",
    "    plt.imshow(1-pred_re[:, :, 0], alpha=0.4, cmap='Reds', vmin=0, vmax=1)\n",
    "    segmentation = img.copy()\n",
    "    segmentation[pred_re[:, :, 0]>0.3] = 0 \n",
    "    segmentation[pred_re[:, :, 0]<=0.3] = 1\n",
    "    mask = binary_fill_holes(segmentation)\n",
    "\n",
    "    plt.figure()\n",
    "    plt.imshow(mask, cmap='gray')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c94eddd-2f08-482b-a7b4-6041141fcd7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#loop trough entire testset and save ROI masks for each file \n",
    "for i, filename in enumerate(test_list[::-1]):\n",
    "    sys.stdout.write('\\r'+f'{i}/{len(test_list)}')    \n",
    "    filename = filename.replace(\".dcm\", \".npy\")\n",
    "    img = np.load(\"./Data/CQ500/4096/\" + filename)*4095 - 1024\n",
    "\n",
    "    pred = U_Net.predict(img[None, :, :, None], verbose=False)\n",
    "    pred_re = pred.reshape(512, 512, 17)\n",
    "    \n",
    "    segmentation = img.copy()\n",
    "    segmentation[pred_re[:, :, 0]>0.3] = 0 \n",
    "    segmentation[pred_re[:, :, 0]<=0.3] = 1\n",
    "    \n",
    "    np.save(f\"./Masks/mask_{filename}\", segmentation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9402594e-5d87-4fa9-b2bf-f78c01fb28a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "del U_Net"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9572558c-30fa-47f6-95f9-257c015c428e",
   "metadata": {},
   "source": [
    "## Calculate PSNR/SSIM"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e141606f-79d5-4346-b7e6-f5333adc42d6",
   "metadata": {},
   "source": [
    "### RSNA test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25365e51-a89b-461e-9521-fbe770125a8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.util.arraycrop import crop\n",
    "\n",
    "def ROI_ssim(img1, img2, mask, data_range=1):\n",
    "    \n",
    "    \"\"\"\n",
    "    Calculate the SSIM of two images for a region of interest (ROI). \n",
    "    \n",
    "    Parameters:\n",
    "    img1, img2: array\n",
    "    The two images to be compared\n",
    "    mask: array\n",
    "    binary mask of ROI\n",
    "    data_range: int/float\n",
    "    data range of input data\n",
    "    \n",
    "    returns:\n",
    "    ssim_roi: float\n",
    "    ssim in region of interest\n",
    "    \"\"\"\n",
    "    ssim_value, ssim_map = ssim(img1, img2, data_range=1, full=True)\n",
    "    \n",
    "    #crop to ignore edge effects. Only relevant if mask = 1 for all pixels. \n",
    " \n",
    "    pad = (7 - 1) // 2\n",
    "    ssim_map, mask = crop(ssim_map, pad), crop(mask, pad)\n",
    "    \n",
    "    ssim_roi = np.sum(ssim_map*mask)/np.sum(mask)\n",
    "    return ssim_roi\n",
    "\n",
    "\n",
    "def ROI_psnr(img1, img2, mask, data_range=1):\n",
    "    \n",
    "    \"\"\"\n",
    "    Calculate the PSNR of two images for a region of interest (ROI).\n",
    "    \n",
    "    Parameters:\n",
    "    img1, img2: array\n",
    "    The two images to be compared\n",
    "    mask: array\n",
    "    binary mask of ROI\n",
    "    data_range: int/float\n",
    "    data range of input data\n",
    "    \n",
    "    \n",
    "    returns:\n",
    "    psnr_roi: float\n",
    "    psnr in region of interest\n",
    "    \"\"\"\n",
    "    err_roi = np.sum((img1*mask - img2*mask) ** 2, dtype=np.float64)/np.sum(mask)\n",
    "    psnr_roi = 10 * np.log10((data_range ** 2) / err_roi)\n",
    "    return psnr_roi\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dabd9259-9914-4890-90e7-c7f79811df6e",
   "metadata": {},
   "source": [
    "## check if everything works as expected"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "144985b9-2713-4f6d-9f15-ab536735182a",
   "metadata": {},
   "outputs": [],
   "source": [
    "test_list = pd.read_csv(\"./test_data.csv\")[\"filename\"].values\n",
    "\n",
    "for i, filename in enumerate(test_list[0:3]):\n",
    "    #sys.stdout.write('\\r'+f'{i}/{len(test_list)}')    \n",
    "    filename = filename.replace(\".dcm\", \".npy\")\n",
    "    gt = np.load(\"./Data/4096/\" + filename)/4095\n",
    "    sparse = np.load(\"./Data/256/\" + filename)/4095\n",
    "    mask = np.load(\"./Masks/mask_\" + filename)\n",
    "    mask = binary_fill_holes(mask)\n",
    "    \n",
    "    #uncomment to test if ROI_ functions\n",
    "    #mask = np.ones_like(mask)\n",
    "\n",
    "    #skip iterations where mask is zero\n",
    "    if 1 not in mask:\n",
    "        continue\n",
    "    \n",
    "    fig, (ax1, ax2, ax3) = plt.subplots(1, 3)\n",
    "    ax1.imshow(gt, cmap='gray', vmin=(0+1024)/4095, vmax=(80+1024)/4095)\n",
    "    ax2.imshow(sparse, cmap='gray', vmin=(0+1024)/4095, vmax=(80+1024)/4095)\n",
    "    ax3.imshow(mask, cmap='gray', vmin=0, vmax=1)\n",
    "    \n",
    "    psnr_roi = ROI_psnr(gt, sparse, mask)\n",
    "    ssim_roi = ROI_ssim(gt, sparse, mask)\n",
    "    \n",
    "    psnr_full = psnr(gt, sparse, data_range=1)\n",
    "    ssim_full = ssim(gt, sparse, data_range=1)\n",
    "    \n",
    "    print(psnr_roi, psnr_full, ssim_roi, ssim_full)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "696d99be-44dd-4da0-b9e8-da0a0561a727",
   "metadata": {},
   "source": [
    "## loop trough testset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3b7ef69-9bfb-4b38-a065-7ec096e2404d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.restoration import denoise_tv_chambolle\n",
    "from models import *\n",
    "\n",
    "test_list = pd.read_csv(\"./Data/test_data.csv\")[\"filename\"].values\n",
    "\n",
    "unet = U_Net()\n",
    "\n",
    "\n",
    "weights = {2048:0.001, 1024:0.002, 512:0.009, 256:0.045, 128:0.153, 64:0.570}\n",
    "\n",
    "results = pd.DataFrame()\n",
    "angles = [2048, 1024, 512, 256, 128, 64]\n",
    "#angles = [2048]\n",
    "for angle in angles:\n",
    "    print('\\n', angle)\n",
    "    \n",
    "    #load U-Net for specific angle\n",
    "    unet = U_Net()\n",
    "    unet.load_weights(f'./model_weights/U-Net/{angle}/model_75.h5')\n",
    "    ssim_list_sparse = []\n",
    "    psnr_list_sparse = []\n",
    "\n",
    "    ssim_list_U_Net = []\n",
    "    psnr_list_U_Net = []\n",
    "    \n",
    "    ssim_list_TV = []\n",
    "    psnr_list_TV = []\n",
    "    \n",
    "    ssim_full_list_U_Net = []\n",
    "\n",
    "    for i, filename in enumerate(test_list):\n",
    "        sys.stdout.write('\\r'+f'{i}/{len(test_list)}')    \n",
    "        filename = filename.replace(\".dcm\", \".npy\")\n",
    "        gt = np.load(\"./Data/4096/\" + filename)/4095\n",
    "        sparse = np.load(f\"./Data/{angle}/\" + filename)/4095\n",
    "        mask = np.load(\"./Masks/mask_\" + filename)\n",
    "        mask = binary_fill_holes(mask)\n",
    "\n",
    "        #skip iterations where mask is zero\n",
    "        if 1 not in mask:\n",
    "            continue\n",
    "        \n",
    "        psnr_list_sparse.append(ROI_psnr(gt, sparse, mask))\n",
    "        ssim_list_sparse.append(ROI_ssim(gt, sparse, mask))\n",
    "        \n",
    "        #calculate psnr + ssim for predicted data\n",
    "        img_U_Net = unet.predict(sparse[np.newaxis, :, :, np.newaxis], verbose=False)\n",
    "        img_U_Net = np.squeeze(img_U_Net)\n",
    "        psnr_list_U_Net.append(ROI_psnr(gt, img_U_Net, mask))\n",
    "        ssim_list_U_Net.append(ROI_ssim(gt, img_U_Net, mask))\n",
    "        ssim_full_list_U_Net.append(ssim(gt, img_U_Net, data_range=1))\n",
    "        \n",
    "        #calculate psnr + ssim for tv data\n",
    "        img_tv = denoise_tv_chambolle(sparse, weight=weights[angle])\n",
    "        psnr_list_TV.append(ROI_psnr(gt, img_tv, mask))\n",
    "        ssim_list_TV.append(ROI_ssim(gt, img_tv, mask))\n",
    "    \n",
    "    results[f'psnr_TV_{angle}'] = psnr_list_TV\n",
    "    results[f'ssim_TV_{angle}'] = ssim_list_TV\n",
    "    \n",
    "    results[f'psnr_sparse_{angle}'] = psnr_list_sparse\n",
    "    results[f'ssim_sparse_{angle}'] = ssim_list_sparse\n",
    "\n",
    "    results[f'psnr_U_Net_{angle}'] = psnr_list_U_Net\n",
    "    results[f'ssim_U_Net_{angle}'] = ssim_list_U_Net\n",
    "    results[f'ssim_full_U_Net_{angle}'] = ssim_full_list_U_Net\n",
    "\n",
    "results.to_csv('./ROI_psnr_ssim_RSNA.csv')\n",
    "results.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2feb1b08-f156-4a35-8fb3-e5b29b8cf286",
   "metadata": {},
   "source": [
    "### CQ500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64860fc3-8058-4511-b363-138e6b042619",
   "metadata": {},
   "outputs": [],
   "source": [
    "from skimage.restoration import denoise_tv_chambolle\n",
    "from models import *\n",
    "\n",
    "test_list = os.listdir(\"/data-pool/data_no_backup/ga63cun/CQ500/4096/\")\n",
    "test_list = [element for element in test_list if \".npy\" in element]\n",
    "\n",
    "unet = U_Net()\n",
    "\n",
    "\n",
    "weights = {2048:0.001, 1024:0.002, 512:0.009, 256:0.045, 128:0.153, 64:0.570}\n",
    "\n",
    "results = pd.DataFrame()\n",
    "angles = [2048, 1024, 512, 256, 128, 64]\n",
    "#angles = [2048]\n",
    "for angle in angles:\n",
    "    print('\\n', angle)\n",
    "    \n",
    "    #load U-Net for specific angle\n",
    "    unet = U_Net()\n",
    "    unet.load_weights(f'./model_weights/U-Net/{angle}/model_75.h5')\n",
    "    ssim_list_sparse = []\n",
    "    psnr_list_sparse = []\n",
    "\n",
    "    ssim_list_U_Net = []\n",
    "    psnr_list_U_Net = []\n",
    "    \n",
    "    ssim_list_TV = []\n",
    "    psnr_list_TV = []\n",
    "    \n",
    "    ssim_full_list_U_Net = []\n",
    "\n",
    "    for i, filename in enumerate(test_list):\n",
    "        sys.stdout.write('\\r'+f'{i}/{len(test_list)}')    \n",
    "        filename = filename.replace(\".dcm\", \".npy\")\n",
    "        gt = np.load(\"./CQ500/Data/4096/\" + filename)/4095\n",
    "        sparse = np.load(f\"./CQ500/Data/{angle}/\" + filename)/4095\n",
    "        mask = np.load(\"./Masks/mask_\" + filename)\n",
    "        mask = binary_fill_holes(mask)\n",
    "\n",
    "        #skip iterations where mask is zero\n",
    "        if 1 not in mask:\n",
    "            continue\n",
    "        \n",
    "        psnr_list_sparse.append(ROI_psnr(gt, sparse, mask))\n",
    "        ssim_list_sparse.append(ROI_ssim(gt, sparse, mask))\n",
    "        \n",
    "        #calculate psnr + ssim for predicted data\n",
    "        img_U_Net = unet.predict(sparse[np.newaxis, :, :, np.newaxis], verbose=False)\n",
    "        img_U_Net = np.squeeze(img_U_Net)\n",
    "        psnr_list_U_Net.append(ROI_psnr(gt, img_U_Net, mask))\n",
    "        ssim_list_U_Net.append(ROI_ssim(gt, img_U_Net, mask))\n",
    "        ssim_full_list_U_Net.append(ssim(gt, img_U_Net, data_range=1))\n",
    "        \n",
    "        #calculate psnr + ssim for tv dataf\n",
    "        img_tv = denoise_tv_chambolle(sparse, weight=weights[angle])\n",
    "        psnr_list_TV.append(ROI_psnr(gt, img_tv, mask))\n",
    "        ssim_list_TV.append(ROI_ssim(gt, img_tv, mask))\n",
    "    \n",
    "    results[f'psnr_TV_{angle}'] = psnr_list_TV\n",
    "    results[f'ssim_TV_{angle}'] = ssim_list_TV\n",
    "    \n",
    "    results[f'psnr_sparse_{angle}'] = psnr_list_sparse\n",
    "    results[f'ssim_sparse_{angle}'] = ssim_list_sparse\n",
    "\n",
    "    results[f'psnr_U_Net_{angle}'] = psnr_list_U_Net\n",
    "    results[f'ssim_U_Net_{angle}'] = ssim_list_U_Net\n",
    "    results[f'ssim_full_U_Net_{angle}'] = ssim_full_list_U_Net\n",
    "\n",
    "    results.to_csv('./ROI_psnr_ssim_CQ500.csv')\n",
    "results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2addbf6-cc94-4499-a634-bc7c6ca6498f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e11837d7-98cb-489a-b72d-dbc78d070176",
   "metadata": {},
   "source": [
    "## Calculate mean values and CIs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fbd8df52-f03f-43da-8061-dd16cc3a22cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats as st\n",
    "\n",
    "def wilco_func(x, y, alternative='two-sided', method='auto'):\n",
    "    out = st.wilcoxon(x, y, alternative=alternative, method=method)\n",
    "    return round(out[0], 5), round(out[1], 5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3271024-6f39-4626-834b-4f978aa0767e",
   "metadata": {},
   "source": [
    "### RSNA test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7068eea9-a3d5-4c95-b8a2-bccb1c8b28f0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "df_mean = pd.DataFrame()\n",
    "df_rsna = pd.read_csv(\"./ROI_psnr_ssim_RSNA.csv\")\n",
    "df_rsna = df_rsna.replace([np.inf, -np.inf], np.nan).dropna()\n",
    "angles = [2048, 1024, 512, 256, 128, 64]\n",
    "\n",
    "for angle in angles:\n",
    "    column = []\n",
    "    \n",
    "    for typ in [\"ssim\", \"psnr\"]:\n",
    "        for proc in [\"sparse\", \"U_Net\", \"TV\"]:\n",
    "\n",
    "            print(angle, proc)\n",
    "            metric_list = df_rsna[f'{typ}_{proc}_{angle}']\n",
    "            res = st.bootstrap((metric_list,), np.mean, n_resamples=1000, confidence_level=0.95)\n",
    "            metric_interval = res.confidence_interval\n",
    "            \n",
    "            column.append(f'{round(metric_list.mean(), 3)}')\n",
    "            column.append(f'({round(metric_interval[0], 3)}-{round(metric_interval[1], 3)})')\n",
    "            \n",
    "    df_mean[f\"{angle}\"] = column\n",
    "\n",
    "            #check for normal distribution\n",
    "            #print(f\"{typ}\", shapiro(metric_list))\n",
    "df_mean.to_csv(\"./table_ROI_SSIM_PSNR_raw.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b614b527-d9e8-4629-a2f7-7aa1cf6d816c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check if distributions are significantly different\n",
    "angles = [2048, 1024, 512, 256, 128, 64]\n",
    "typ=\"ssim\"\n",
    "for angle in angles:\n",
    "    stat_sparse_tv = wilco_func(df_rsna[f'{typ}_sparse_{angle}'], df_rsna[f'{typ}_TV_{angle}'])\n",
    "    stat_sparse_U_Net = wilco_func(df_rsna[f'{typ}_sparse_{angle}'], df_rsna[f'{typ}_U_Net_{angle}'])\n",
    "    stat_TV_U_Net = wilco_func(df_rsna[f'{typ}_TV_{angle}'], df_rsna[f'{typ}_U_Net_{angle}'])\n",
    "\n",
    "    print(f\"{angle}, FBP-TV: {stat_sparse_tv}, FBP-U-Net: {stat_sparse_U_Net}, TV-U-Net: {stat_TV_U_Net}\")\n",
    "    \n",
    "typ=\"psnr\"\n",
    "for angle in angles:\n",
    "    stat_sparse_tv = wilco_func(df_rsna[f'{typ}_sparse_{angle}'], df_rsna[f'{typ}_TV_{angle}'])\n",
    "    stat_sparse_U_Net = wilco_func(df_rsna[f'{typ}_sparse_{angle}'], df_rsna[f'{typ}_U_Net_{angle}'])\n",
    "    stat_TV_U_Net = wilco_func(df_rsna[f'{typ}_TV_{angle}'], df_rsna[f'{typ}_U_Net_{angle}'])\n",
    "\n",
    "    print(f\"{angle}, FBP-TV: {stat_sparse_tv}, FBP-U-Net: {stat_sparse_U_Net}, TV-U-Net: {stat_TV_U_Net}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "748498ab-48e8-4638-b98a-e595c0a2d4c2",
   "metadata": {},
   "source": [
    "### CQ500 test split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f320b0ae-6638-4167-acdd-e111ac32db9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_mean = pd.DataFrame()\n",
    "df_cq = pd.read_csv(\"./ROI_psnr_ssim_CQ500.csv\")\n",
    "df_cq = df_cq.replace([np.inf, -np.inf], np.nan).dropna()\n",
    "angles = [2048, 1024, 512, 256, 128, 64]\n",
    "\n",
    "for angle in angles:\n",
    "    column = []\n",
    "    \n",
    "    for typ in [\"ssim\", \"psnr\"]:\n",
    "        for proc in [\"sparse\", \"U_Net\", \"TV\"]:\n",
    "\n",
    "            print(angle, proc)\n",
    "            metric_list = df_cq[f'{typ}_{proc}_{angle}']\n",
    "            res = st.bootstrap((metric_list,), np.mean, n_resamples=1000, confidence_level=0.95)\n",
    "            metric_interval = res.confidence_interval\n",
    "            \n",
    "            column.append(f'{round(metric_list.mean(), 3)}')\n",
    "            column.append(f'({round(metric_interval[0], 3)}-{round(metric_interval[1], 3)})')\n",
    "            \n",
    "    df_mean[f\"{angle}\"] = column\n",
    "\n",
    "            #check for normal distribution\n",
    "            #print(f\"{typ}\", shapiro(metric_list))\n",
    "df_mean.to_csv(\"./tableCQ500_ROI_SSIM_PSNR_raw.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acaa2f15-b07c-4741-ab39-0d86376dd5b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#check if distributions are significantly different\n",
    "angles = [2048, 1024, 512, 256, 128, 64]\n",
    "typ=\"ssim\"\n",
    "for angle in angles:\n",
    "    stat_sparse_tv = wilco_func(df_cq[f'{typ}_sparse_{angle}'], df_cq[f'{typ}_TV_{angle}'])\n",
    "    stat_sparse_U_Net = wilco_func(df_cq[f'{typ}_sparse_{angle}'], df_cq[f'{typ}_U_Net_{angle}'])\n",
    "    stat_TV_U_Net = wilco_func(df_cq[f'{typ}_TV_{angle}'], df_cq[f'{typ}_U_Net_{angle}'])\n",
    "\n",
    "    print(f\"{angle}, FBP-TV: {stat_sparse_tv}, FBP-U-Net: {stat_sparse_U_Net}, TV-U-Net: {stat_TV_U_Net}\")\n",
    "    \n",
    "typ=\"psnr\"\n",
    "for angle in angles:\n",
    "    stat_sparse_tv = wilco_func(df_cq[f'{typ}_sparse_{angle}'], df_cq[f'{typ}_TV_{angle}'])\n",
    "    stat_sparse_U_Net = wilco_func(df_cq[f'{typ}_sparse_{angle}'], df_cq[f'{typ}_U_Net_{angle}'])\n",
    "    stat_TV_U_Net = wilco_func(df_cq[f'{typ}_TV_{angle}'], df_cq[f'{typ}_U_Net_{angle}'])\n",
    "\n",
    "    print(f\"{angle}, FBP-TV: {stat_sparse_tv}, FBP-U-Net: {stat_sparse_U_Net}, TV-U-Net: {stat_TV_U_Net}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "192899f8-d106-42c1-9193-65039d08da02",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}