V5_run_sma_preprocess.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load libraries and directories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display, HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from IPython import get_ipython\n",
    "from tqdm.notebook import tqdm\n",
    "import pickle\n",
    "import os\n",
    "import pprint\n",
    "pp = pprint.PrettyPrinter(indent=1)\n",
    "\n",
    "# Custom modules for debugging\n",
    "from SliceViewer import ImageSliceViewer3D, ImageSliceViewer3D_1view,ImageSliceViewer3D_2views\n",
    "from investigate import *\n",
    "\n",
    "#pd.set_option(\"display.max_rows\", 10)\n",
    "      \n",
    "import json\n",
    "from run_sma_experiment import find_l3_images,output_images\n",
    "import pprint\n",
    "from L3_finder import *\n",
    "\n",
    "# Custom functions\n",
    "def save_object(obj, filename):\n",
    "    with open(filename, 'wb') as output:  # Overwrites any existing file.\n",
    "        pickle.dump(obj, output, pickle.HIGHEST_PROTOCOL)\n",
    "\n",
    "def load_object(filename):        \n",
    "    with open(filename, 'rb') as input:\n",
    "        return pickle.load(input)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_ipython().run_line_magic('tb', '')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In CCHMC's workflow, there were two dicom dumps, with different folder structures and naming convention:\n",
    "<br>\n",
    "Folder structure: <br>\n",
    "Dump-1: Project_folder/Patient_folder/Series_Folder/dicom_files <br>\n",
    "Dump-2: Project_folder/Patient_folder/Study_folder/Series_Folder/dicom_files<br> \n",
    "<br>\n",
    "Naming Convention for patient folder: <br>\n",
    "Dump-1: PATID-GMRN-PATID-STUDYNAME <br>\n",
    "Dump-2: PT-PATID-PATID "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Select which dump you are processing here [dump1: 1, dump2: 2]\n",
    "dump = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cwd = os.getcwd()\n",
    "data = '/tf/data'\n",
    "output = '/tf/pickles/v5_8pts'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 1 - Load list of normal patients filtered from Epic data and select those patients from the DICOM dump of all patients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load normal patient list\n",
    "infile  = 'patlist_with_validBMI_corrected_v5.csv'\n",
    "df_P = pd.read_csv(infile, index_col=False)\n",
    "df_P = df_P.loc[:, ~df_P.columns.str.contains('^Unnamed')]\n",
    "df_P = df_P[['GIVEN_MRN','PAT_ID','ACC']]\n",
    "print('Columns of df_P: ', list(df_P))\n",
    "print('Length of df_P: ', len(df_P))\n",
    "display(df_P.head(10))\n",
    "#print('# of Unique patients: ', len(df_P.subject_id.unique()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pats = next(os.walk(data))[1]\n",
    "print('Total patient folders in data dir: ',len(pats))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if dump == 1:\n",
    "    patids = [pat.split('-')[0] for pat in pats]\n",
    "elif dump == 2:\n",
    "    patids = [pat.split('-')[-1] for pat in pats]\n",
    "\n",
    "valid_ids = [valid_id for valid_id,valid_dir in zip(patids,pats) if valid_id in df_P.PAT_ID.values]\n",
    "\n",
    "valid_ids = set(valid_ids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('valid ids: ',len(valid_ids))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2 - Load each study into subject object\n",
    "<br>\n",
    "Subject object defined in L3finder.ingest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Import modules and config file\n",
    "configfile = os.path.join(cwd,'config/debug_ES/series_filter_ds1.json')\n",
    "with open(configfile, \"r\") as f:\n",
    "        config = json.load(f)\n",
    "\n",
    "config = config[\"series_filter\"]        \n",
    "print('Current config dict: ')\n",
    "pp.pprint(config)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if dump==2:\n",
    "    config[\"new_tim_dicom_dir_structure\"] = False\n",
    "elif dump==1:\n",
    "    config[\"new_tim_dicom_dir_structure\"] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Debug\n",
    "print(\"Finding subjects\")\n",
    "\n",
    "subjects = list(\n",
    "    find_subjects(\n",
    "        config[\"dicom_dir\"],\n",
    "        new_tim_dir_structure=config[\"new_tim_dicom_dir_structure\"]\n",
    "    )\n",
    ")\n",
    "\n",
    "print('Subjects found: ', len(subjects))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section-3 - check if there are subjects with multiple folders (studies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subjects = [subject for subject in subjects if subject.id_ in valid_ids]\n",
    "print('Subjects found: ', len(subjects))\n",
    "print('Valid Subjects: ', len(valid_ids))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find Duplicate Subjects\n",
    "unique_subjects = []\n",
    "duplicate_subjects = []\n",
    "for subject in subjects:\n",
    "    if subject.id_ not in unique_subjects:\n",
    "        unique_subjects.append(subject.id_)\n",
    "    else:\n",
    "        duplicate_subjects.append(subject.id_)\n",
    "\n",
    "print('Duplicates: ',len(duplicate_subjects)           )\n",
    "print('Uniques: ',len(unique_subjects))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Save subjects without duplicates\n",
    "save_object(subjects, os.path.join(output,'subjects_noduplicates.pkl'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 4 - Load each series into series object and keep only axials and sagittals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "subjects = load_object(os.path.join(output,'subjects_noduplicates.pkl'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(subjects)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Find series images\n",
    "print(\"Finding series\")\n",
    "series = list(flatten(tqdm((s.find_series() for s in subjects),total=len(subjects))))\n",
    "print(\"Total number of series found: \", len(series))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Separate axial and sagittal series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# Debug\n",
    "from L3_finder import *\n",
    "from l3finder.ingest import *\n",
    "from multiprocessing import get_context\n",
    "from multiprocessing import set_start_method\n",
    "#set_start_method(\"spawn\")\n",
    "\n",
    "if __name__=='__main__':\n",
    "    # Find series images\n",
    "    print(\"Finding series\")\n",
    "    series = list(flatten(s.find_series() for s in subjects))\n",
    "\n",
    "    # Separate series\n",
    "    print(\"Separating series\")\n",
    "    #sagittal_series, axial_series, excluded_series = separate_series(series)\n",
    "    \n",
    "    excluded_series = []\n",
    "\n",
    "    sag_filter = functools.partial(\n",
    "        same_orientation,\n",
    "        orientation='sagittal',\n",
    "        excluded_series=excluded_series\n",
    "    )\n",
    "    \n",
    "    axial_filter = functools.partial(\n",
    "        same_orientation,\n",
    "        orientation='axial',\n",
    "        excluded_series=excluded_series\n",
    "    )\n",
    "\n",
    "    def pool_filter(pool, func, candidates):\n",
    "        return [\n",
    "            c for c, keep\n",
    "            in zip(candidates, tqdm(pool.imap(func, candidates),total=len(candidates)))\n",
    "            if keep]\n",
    "    \n",
    "    print('Filtering series using ', multiprocessing.cpu_count(), ' cores:')\n",
    "    with get_context(\"spawn\").Pool() as p:\n",
    "        sagittal_series = pool_filter(p, sag_filter, series)\n",
    "        print(\"Processed Sagittals\")\n",
    "        axial_series = pool_filter(p, axial_filter, series)\n",
    "        print(\"Processed Axials\")\n",
    "        p.close()\n",
    "        p.join()\n",
    "\n",
    "    \n",
    "    \n",
    "    print(\"Series seperated\")\n",
    "\n",
    "#remove_start_method(\"spawn\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Length of valid pats: \", len(subjects))\n",
    "print(\"Length of sagittal series\", len(sagittal_series))\n",
    "print(\"Length of axial series\", len(axial_series))\n",
    "#print(\"Length of excluded series\", len(excluded_series))\n",
    "#print(\"Length of all series in dataset\", len(series))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save required objects\n",
    "save_object(axial_series, os.path.join(output,'axial_series.pkl'))\n",
    "save_object(sagittal_series, os.path.join(output,'sagittal_series.pkl'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 5 - Investigate subjects and series using pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axial_series = load_object(os.path.join(output,'axial_series.pkl'))\n",
    "sagittal_series = load_object(os.path.join(output,'sagittal_series.pkl'))\n",
    "subjects = load_object(os.path.join(output,'subjects_noduplicates.pkl'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_a = get_summary_dfs(axial_series,sagittal_series,subjects)\n",
    "save_object(df_a, os.path.join(output,'df_a.pkl'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(df_a.head(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_a_axials = get_summary_by_serieslength(axial_series)\n",
    "df_a_sags = get_summary_by_serieslength(sagittal_series)\n",
    "save_object(df_a_axials, os.path.join(output,'df_a_axials.pkl'))\n",
    "save_object(df_a_sags, os.path.join(output,'df_a_sags.pkl'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Length of subjects with atleast 1 axial or sagittal series: \", len(df_a))\n",
    "print(\"Length of subjects with atleast 1 axial series: \", len(df_a_axials['ID'].unique()))\n",
    "print(\"Length of subjects with atleast 1 sagittal series: \", len(df_a_sags['ID'].unique()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Patients without Axial\n",
    "pats = [pat for pat in df_a['ID'].values if pat not in df_a_axials['ID'].values]\n",
    "print(len(pats))\n",
    "print(pats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Patients without Sagittal\n",
    "pats = [pat for pat in df_a['ID'].values if pat not in df_a_sags['ID'].values]\n",
    "print(len(pats))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imseries = get_subject_series('Z837620','Z837620-SE-6-Vol_Body_Vol._0.5',subjects)\n",
    "print(imseries.orientation,' ' , imseries.slice_thickness)\n",
    "imdata = imseries.pixel_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "print(imdata.shape)\n",
    "ImageSliceViewer3D(imdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_summary_by_serieslength(df_a_axials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print_summary_by_serieslength(df_a_sags)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 6 - Create dataframe of optimal axial sagittal pairs\n",
    "\n",
    "The function filter_finalpairs in investigate.py is used"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axial_series = load_object(os.path.join(output,'axial_series.pkl'))\n",
    "sagittal_series = load_object(os.path.join(output,'sagittal_series.pkl'))\n",
    "subjects = load_object(os.path.join(output,'subjects_noduplicates.pkl'))\n",
    "\n",
    "df_a_axials = load_object(os.path.join(output,'df_a_axials.pkl'))\n",
    "df_a_sags = load_object(os.path.join(output,'df_a_sags.pkl'))\n",
    "df_a = load_object(os.path.join(output,'df_a.pkl'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "from L3_finder import *\n",
    "from l3finder.ingest import *\n",
    "from multiprocessing import get_context\n",
    "from multiprocessing import set_start_method\n",
    "#set_start_method(\"spawn\")\n",
    "df_filt = None\n",
    "\n",
    "if __name__=='__main__':\n",
    "    # Find series images\n",
    "    print(\"Finding IDs\")\n",
    "    \n",
    "    IDs = [s.id_ for s in subjects]\n",
    "    pair_filter = functools.partial(\n",
    "        filter_finalpairs,\n",
    "        df_ax=df_a_axials,\n",
    "        df_sag=df_a_sags,\n",
    "        subjects=subjects\n",
    "    )\n",
    "    \n",
    "    def pool_filter(pool, func, candidates):\n",
    "        return [a for a in tqdm(pool.imap(func, candidates),total=len(candidates))]\n",
    "        \n",
    "    print('Filtering series using ', multiprocessing.cpu_count(), ' cores:')\n",
    "        \n",
    "    with get_context(\"spawn\").Pool(4) as p:\n",
    "        result_list = pool_filter(p, pair_filter, IDs)\n",
    "        p.close()\n",
    "        p.join()\n",
    "    \n",
    "    print('parallel processing over')\n",
    "     # Start from here\n",
    "    df_filt  = pd.DataFrame(columns=['ID','Axial','Sagittal','Overlap','MissingScore','PairValidity', \n",
    "                                'AxSlices','SagSlices','AxThick','SagThick'])\n",
    "    for i,op in enumerate(result_list):\n",
    "        df_filt.loc[i] = op\n",
    "    \n",
    "\n",
    "print(\"Processed\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_object(df_filt, os.path.join(output,'df_filteredpairs.pkl'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 7: Investigate the dataframe for missing and low quality pairs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Load all params\n",
    "df_filt = load_object(os.path.join(output,'df_filteredpairs.pkl'))\n",
    "axial_series = load_object(os.path.join(output,'axial_series.pkl'))\n",
    "sagittal_series = load_object(os.path.join(output,'sagittal_series.pkl'))\n",
    "subjects = load_object(os.path.join(output,'subjects_noduplicates.pkl'))\n",
    "\n",
    "df_a_axials = load_object(os.path.join(output,'df_a_axials.pkl'))\n",
    "df_a_sags = load_object(os.path.join(output,'df_a_sags.pkl'))\n",
    "df_a = load_object(os.path.join(output,'df_a.pkl'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make sure filtered df and subjects are equal length\n",
    "print('Length of filtered df: ',len(df_filt))\n",
    "print('Length of subjects: ',len(subjects))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# View and Remove subjects without axial series\n",
    "df_noaxials =  df_filt[df_filt['Axial'].isnull()]\n",
    "print('Number of subjects without Axials: ', len(df_noaxials))\n",
    "display(df_noaxials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove subjects without axials from subjects list:\n",
    "subjects = [s for s in subjects if s.id_ not in df_noaxials.ID.values]\n",
    "print('Length of subjects with axials: ',len(subjects))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print cases that don't have sagittals\n",
    "df_nosags = df_filt[df_filt['Sagittal'].isnull()]\n",
    "print(\"Number of cases without sagittals: \", len(df_nosags))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Investigate cases with less than 0.7 overlap and  < 0.9 Missing Score [Tracks Slices missing from stack]\n",
    "df_pooroverlap = df_filt[(df_filt['Overlap'] < 0.7) | (df_filt['MissingScore'] < 0.9)]\n",
    "print('Cases with overlap < 0.7: ', len(df_pooroverlap))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "display(df_pooroverlap.sort_values(by=['MissingScore'],ascending=[True]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Handy Functions to investigate the poor pairs\n",
    "def get_ax_sag(df,ind):\n",
    "    global subjects\n",
    "    subid = df.loc[ind,'ID']\n",
    "    axid = df.loc[ind,'Axial']\n",
    "    sagid = df.loc[ind,'Sagittal']\n",
    "    ax = get_subject_series(subid,axid,subjects)\n",
    "    sag = get_subject_series(subid,sagid,subjects)\n",
    "    return ax,sag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calculate_missing_slices_sagittals(get_ax_sag(df_pooroverlap,85)[1]),verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calculate_series_overlap(*get_ax_sag(df_pooroverlap,85),verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Based on investigation, eliminate series and subjects not eligible and create final df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Keep Sagittals only when overlap > 0.7\n",
    "df_final = df_filt.copy()\n",
    "for ind,row in df_final.iterrows():\n",
    "    if (not row['Overlap']) or (row['Overlap'] < 0.7):\n",
    "            df_final.loc[ind,'Sagittal'] = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print cases that don't have sagittals\n",
    "print(\"Number of cases without sagittals in filter df: \", len(df_nosags))\n",
    "df_nosags2 = df_final[df_final['Sagittal'].isnull()]\n",
    "print(\"Number of cases without sagittals in final df: \", len(df_nosags2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get final \n",
    "\n",
    "final_df_file = 'df_final_dump_'+str(dump)+'_8pats.pkl'\n",
    "final_subs_file = 'subjects_final_dump_'+str(dump)+'_8pats.pkl'\n",
    "\n",
    "save_object(df_final, os.path.join(output,final_df_file))\n",
    "save_object(subjects, os.path.join(output,final_subs_file))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "file_extension": ".py",
  "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"
  },
  "mimetype": "text/x-python",
  "name": "python",
  "npconvert_exporter": "python",
  "pygments_lexer": "ipython3",
  "version": 3
 },
 "nbformat": 4,
 "nbformat_minor": 4
}