code/gamos-geom-mpt-helper.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# GAMOS Geometry file creation\n",
    "\n",
    "GAMOS requires a geometry file to define the structure of world used in each simulations. When optical properties are needed for the simulation, they must be defined for each material at every wavelength of interest. This is accomplished through material property tables. Writing these tables manually is time-consuming and prone to input errors. The following code blocks are used to read an Excel file and create a multi-layer geometry with a variety of properties. The code is currently limited to BOX geometries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas\n",
    "import numpy\n",
    "# Define input and output directories\n",
    "input_dir = '../data/'\n",
    "output_dir = '../results/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def print_section_break(title, output_str=''):\n",
    "    output_str += '\\n//{}//\\n'.format('-'*76)\n",
    "    output_str += '// {}\\n'.format(title.upper())\n",
    "    return output_str\n",
    "def print_materail_init(materials, output_str=''):\n",
    "    for each_material in sorted(materials):\n",
    "        \n",
    "        all_keys = list(materials[each_material]['key'])\n",
    "        name_index = all_keys.index('name')\n",
    "        if not materials[each_material]['value'][name_index].startswith('G4_'):\n",
    "            density_index = all_keys.index('density')\n",
    "            mixc_index = all_keys.index('mixture_content')\n",
    "            mixr_index = all_keys.index('mixture_ratio')\n",
    "\n",
    "\n",
    "            material_names = materials[each_material]['value'][mixc_index].split(',')\n",
    "            num_materials = len(material_names)\n",
    "            material_ratios = materials[each_material]['value'][mixr_index].split(',')\n",
    "\n",
    "\n",
    "            output_str += '\\n:MIXT\\t{}\\t{}\\t{}\\n'.format(materials[each_material]['value'][name_index],\n",
    "                                                       materials[each_material]['value'][density_index],\n",
    "                                                       num_materials)\n",
    "            i=0\n",
    "            for each_name in material_names:\n",
    "                output_str += '{}\\t{}\\n'.format(each_name,material_ratios[i])\n",
    "                i += 1\n",
    "    output_str += '\\n'\n",
    "    return output_str\n",
    "    \n",
    "def get_property_list(input_series, delim):\n",
    "    cur_list = list(input_series.round(7).values.astype(str))\n",
    "    output_str = delim.join(cur_list)\n",
    "    return output_str\n",
    "    \n",
    "def print_material_properties(materials, material_properties, search='tissue', output_str=''):\n",
    "    \n",
    "    for each_material in sorted(materials):\n",
    "        if search in each_material:\n",
    "            all_keys = list(materials[each_material]['key'])\n",
    "            mpt_index = all_keys.index('mpt_name')\n",
    "            cur_mpt_name = materials[each_material]['value'][mpt_index]\n",
    "\n",
    "            name_index = all_keys.index('name')\n",
    "            cur_material = materials[each_material]['value'][name_index]\n",
    "\n",
    "            all_wavelengths_str = get_property_list(material_properties[each_material]['mpt_wavelength'], '*wvnm ')\n",
    "            all_rindex_str = get_property_list(material_properties[each_material]['mpt_rindex'], ' ')\n",
    "\n",
    "            output_str += '\\n:MATE_PROPERTIES_TABLE\\t{}\\n'.format(cur_mpt_name)\n",
    "            output_str += ':MATEPT_ADD_WAVELENGTHS\\t{}\\t{}*wvnm\\n'.format(cur_mpt_name, all_wavelengths_str)\n",
    "            output_str += ':MATEPT_ADD_PROPERTY\\t{}\\tRINDEX\\t{}\\n'.format(cur_mpt_name,all_rindex_str )\n",
    "            if 'mpt_miescat' in material_properties[each_material].columns:\n",
    "                all_mus_str = get_property_list(material_properties[each_material]['mpt_miescat'], '*/mm ')\n",
    "                output_str += ':MATEPT_ADD_PROPERTY\\t{}\\tMIE_SCATCOEF\\t{}*/mm\\n'.format(cur_mpt_name,all_mus_str )\n",
    "            if 'mpt_abs' in material_properties[each_material].columns:\n",
    "                all_mua_str = get_property_list(material_properties[each_material]['mpt_abs'], '*/mm ')\n",
    "                output_str += ':MATEPT_ADD_PROPERTY\\t{}\\tABSCOEF\\t{}*/mm\\n'.format(cur_mpt_name,all_mua_str )\n",
    "            if 'mpt_g' in material_properties[each_material].columns:\n",
    "                all_gval_str = get_property_list(material_properties[each_material]['mpt_g'], ' ')\n",
    "                output_str += ':MATEPT_ADD_PROPERTY\\t{}\\tMIE_GVALUE\\t{}\\n'.format(cur_mpt_name,all_gval_str )\n",
    "            if 'mpt_fl_abs' in material_properties[each_material].columns:\n",
    "                all_fl_absval_str = get_property_list(material_properties[each_material]['mpt_fl_abs'], '*/mm ' )\n",
    "                output_str += ':MATEPT_ADD_PROPERTY\\t{}\\tFLUOR_ABSCOEF\\t{}*/mm\\n'.format(cur_mpt_name,all_fl_absval_str )\n",
    "            if 'mpt_fl_ems' in material_properties[each_material].columns:\n",
    "                all_fl_emsval_str = get_property_list(material_properties[each_material]['mpt_fl_ems'], ' ')\n",
    "                output_str += ':MATEPT_ADD_PROPERTY\\t{}\\tFLUOR_EMISSION\\t{}\\n'.format(cur_mpt_name,all_fl_emsval_str )\n",
    "                output_str += ':MATEPT_ADD_CONST_PROPERTY\\t{}\\tFLUOR_QUANTUMYIELD\\t{:1.1f}\\n'.format(cur_mpt_name, float(qys[each_material]))\n",
    "                output_str += ':MATEPT_ADD_CONST_PROPERTY\\t{}\\tFLUOR_LIFETIME\\t{}\\n'.format(cur_mpt_name, lifetimes[each_material]) \n",
    "                                                                                                   \n",
    "                                                                                                   \n",
    "            output_str += ':MATEPT_ATTACH_TO_MATERIAL\\t{}\\t{}\\n'.format(cur_mpt_name, cur_material)\n",
    "    output_str += '\\n'\n",
    "    return output_str\n",
    "    \n",
    "def print_build_and_place_volumes(volumes, output_str=''):\n",
    "    rm_name = 'RM00'\n",
    "    output_str += ':ROTM\\t{}\\t0.\\t0.\\t0.\\n'.format(rm_name)\n",
    "    for each_name in volumes['name'].values:\n",
    "        #print(each_volume)\n",
    "        cur_volume = all_volumes[all_volumes['name']==each_name]\n",
    "        geom = cur_volume['geometry'].values[0]\n",
    "        x = cur_volume['x'].values[0]\n",
    "        x_unit = cur_volume['x_unit'].values[0]\n",
    "        y = cur_volume['y'].values[0]\n",
    "        y_unit = cur_volume['y_unit'].values[0]\n",
    "        z = cur_volume['z'].values[0]\n",
    "        z_unit = cur_volume['z_unit'].values[0]\n",
    "        material = cur_volume['material'].values[0]\n",
    "        vis = cur_volume['vis'].values[0]\n",
    "        output_str += '\\n:VOLU\\t{}\\t{}\\t{}*{}\\t{}*{}\\t{}*{}\\t{}\\n'.format(each_name,\n",
    "                                                                          geom,\n",
    "                                                                          x,x_unit,\n",
    "                                                                          y,y_unit,\n",
    "                                                                          z,z_unit,\n",
    "                                                                          material)\n",
    "        output_str += ':VIS\\t{}\\t{}\\n'.format(each_name, vis.upper())\n",
    "        \n",
    "        if each_name != 'world':\n",
    "            place_x = cur_volume['place_x'].astype(float).values[0]\n",
    "            place_x_unit = cur_volume['place_x_unit'].values[0]\n",
    "            place_y = cur_volume['place_y'].astype(float).values[0]\n",
    "            place_y_unit = cur_volume['place_y_unit'].values[0]\n",
    "            place_z = cur_volume['place_z'].astype(float).values[0]\n",
    "            place_z_unit = cur_volume['place_z_unit'].values[0]\n",
    "            mother = cur_volume['mother'].values[0]\n",
    "            \n",
    "            \n",
    "            output_str += ':PLACE\\t{}\\t1\\t{}\\t{}\\t{:0.2f}*{}\\t{:0.2f}*{}\\t{:0.2f}*{}\\n'.format(each_name,\n",
    "                                                                                            mother,\n",
    "                                                                                            rm_name,\n",
    "                                                                                            place_x, place_x_unit,\n",
    "                                                                                            place_y, place_y_unit,\n",
    "                                                                                            place_z, place_z_unit)\n",
    "    return output_str\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Define input Excel file\n",
    "raw_file = os.path.join(input_dir, 'PtG4_MC_geometry_LaRochelle.xlsx')\n",
    "# First few rows have generic material properties \n",
    "#which are not wavelength-dependent\n",
    "header_row = 7\n",
    "xls_raw = pandas.ExcelFile(raw_file)\n",
    "all_materials = {}\n",
    "all_materials_meta = {}\n",
    "active_sheets =[]\n",
    "lifetimes = {}\n",
    "qys = {}\n",
    "\n",
    "# Each material has it's own sheet in the Excel which define the properties\n",
    "# The keywords 'skin', 'fat', or 'muscle' must be in the sheet name\n",
    "# If the keyword 'tumor' is in the sheet name, additional fluoresence properties will be parsed\n",
    "# If 'air' is in the sheet name, less properties are needed\n",
    "# The keyword 'volumes' should only be in one sheet name. This defines the size and placement of the BOX geometries\n",
    "for each_sheet_name in xls_raw.sheet_names:\n",
    "    if ('skin' in each_sheet_name) or ('fat' in each_sheet_name) or ('muscle' in each_sheet_name):\n",
    "        active_sheets.append(each_sheet_name)\n",
    "        all_materials[each_sheet_name] = xls_raw.parse(each_sheet_name, header=header_row-1)\n",
    "        all_materials_meta[each_sheet_name] = xls_raw.parse(each_sheet_name, usecols=[0,1], skip_footer=all_materials[each_sheet_name].shape[0]+1)\n",
    "    elif 'tumor' in each_sheet_name:\n",
    "        all_materials[each_sheet_name] = xls_raw.parse(each_sheet_name, header=header_row+2-1)\n",
    "        all_materials_meta[each_sheet_name] = xls_raw.parse(each_sheet_name, usecols=[0,1], skip_footer=all_materials[each_sheet_name].shape[0]+1)\n",
    "        lifetimes[each_sheet_name] = all_materials_meta['tumor'][all_materials_meta['tumor']['key']=='mpt_lifetime']['value'].values[0]\n",
    "        qys[each_sheet_name] = all_materials_meta['tumor'][all_materials_meta['tumor']['key']=='mpt_qy']['value'].values[0]\n",
    "    elif 'air' in each_sheet_name:\n",
    "        all_materials[each_sheet_name] = xls_raw.parse(each_sheet_name, header=header_row-2-1)\n",
    "        all_materials_meta[each_sheet_name] = xls_raw.parse(each_sheet_name, usecols=[0,1], skip_footer=all_materials[each_sheet_name].shape[0]+1)\n",
    "    elif 'volumes' in each_sheet_name:\n",
    "        all_volumes = xls_raw.parse(each_sheet_name)\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_volumes[all_volumes['name']=='world']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_volumes['name'].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_materials_meta.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "all_materials['skin1'].columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_materials['tumor'].columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Create a single string defining the geometry file\n",
    "file_str = ''\n",
    "# Create custom materials of custom density by mixing predefined matirals at specific ratios\n",
    "file_str = print_section_break('Initialize the materials', output_str=file_str)\n",
    "file_str = print_materail_init(all_materials_meta, output_str=file_str)\n",
    "# Define the optical properties of all skin layes\n",
    "file_str = print_section_break('Define skin properties (Layers 1-7)', output_str=file_str)\n",
    "file_str = print_material_properties(all_materials_meta, all_materials, search='skin', output_str=file_str)\n",
    "# Define optical properties of fat and muscle layers\n",
    "file_str = print_section_break('Define fat and muscle', output_str=file_str)\n",
    "file_str = print_material_properties(all_materials_meta, all_materials, search='fat', output_str=file_str)\n",
    "file_str = print_material_properties(all_materials_meta, all_materials, search='muscle', output_str=file_str)\n",
    "# Define tumor optical properties, including fluoresence\n",
    "file_str = print_section_break('Define tumor inclusion properties\\n', output_str=file_str)\n",
    "file_str = print_material_properties(all_materials_meta, all_materials, search='tumor', output_str=file_str)\n",
    "# Add air optical properties\n",
    "file_str = print_section_break('Define general world optical properties\\n', output_str=file_str)\n",
    "file_str = print_material_properties(all_materials_meta, all_materials, search='air', output_str=file_str)\n",
    "# Define material shape and placements\n",
    "file_str = print_section_break('Build and place volumes', output_str=file_str)\n",
    "file_str = print_build_and_place_volumes(all_volumes, output_str=file_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "print(file_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Write the above sting to an output file\n",
    "# Note: For the simulations in the corresponding manuscript the line\n",
    "#  :VOLU\ttumor_inclusion\tBOX\t0.5*cm\t0.5*cm\t5.0*mm\ttissue_tumor\n",
    "# was replaced with\n",
    "#  :VOLU\ttumor_inclusion\tORB\t5.0*mm\ttissue_tumor\n",
    "with open(os.path.join(output_dir, 'world.geom'), 'w') as f:\n",
    "    f.write(file_str)\n",
    "          "
   ]
  },
  {
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}