In this tutorial we will learn how to
The MITK-MSI software provides a wrapper to the popular MCML approach to simulate how light travels through tissue. This wrapper can be found in mc/sim.py. In this tutorial we will utilize our tissue model which uses this wrapper to create the reflectance spectra.
As a prerequisit, you need a MCML monte carlo simulation which uses the format specified here. I tested this software with the GPU accelerated version which can be found here.
In [1]:
# 1.1 create spectra - setup simulation environment
# some necessary imports
import logging
import numpy as np
import os
# everything related to the simulation wrapper
from mc import sim
# the factories create batches (tissue samples) and suited tissue models
from mc import factories
# the tissueparser reads a tissue init file
# examples can be found in mc/data/tissues
from mc import tissueparser
# Where does your monte carlo simulation executable resides in?
MCML_EXECUTABLE = "/home/wirkert/workspace/monteCarlo/gpumcml/fast-gpumcml/gpumcml.sm_20"
# The wavelengths for which we want to run our simulation
WAVELENGTHS = np.arange(450, 978, 2) * 10 ** -9
# how many spectra do you want to simulate:
NR_ELEMENTS_IN_BATCH = 10
# number of batch to simulate
BATCH_NR = 0
# batch_prefix
BATCH_PREFIX = "tutorial_batch"
# the file where the tissue is specified.
# switch out this file with your tissue specification to
# create samples from a different tissue spec.
TISSUE_CONFIG_FILE = "../mc/data/tissues/laparoscopic_ipcai_colon_2016_08_23.ini"
# The MCML needs a simulation input file, where shall it be created?
MCI_FILENAME = "./" + BATCH_PREFIX + "_" + str(BATCH_NR) + ".mci"
# we want to create standard colonic tissue as specified in the IPCAI 2016 publication
# "Robust Near Real-Time Estimation of Physiological Parameters from Megapixel
# Multispectral Images with Inverse Monte Carlo and Random Forest Regression"
factory = factories.GenericMcFactory()
# create a simulation wrapper
sim_wrapper = sim.SimWrapper()
# our simulation needs to know where the input file for the simulation
# shall resign (will be automatically created)
sim_wrapper.set_mci_filename(MCI_FILENAME)
# also it needs to know where the simulation executable shall lie in
sim_wrapper.set_mcml_executable(MCML_EXECUTABLE)
# create the tissue model
# it is responsible for writing the simulation input file
tissue_model = factory.create_tissue_model()
# tell it where the input file shall lie in
tissue_model.set_mci_filename(sim_wrapper.mci_filename)
# tell it how much photons shall be simulated. Will be set to 10**6 by standard,
# this is just an example
tissue_model.set_nr_photons(10**6)
# enter total number of simulations for mci file generation
total_nr_mc_simulations = NR_ELEMENTS_IN_BATCH * len(WAVELENGTHS)
tissue_model._mci_wrapper.set_nr_runs(total_nr_mc_simulations)
# create the header for the mci file
tissue_model.create_mci_file()
In [3]:
# 1.2 create spectra - create tissue samples for simulation
# create a tissueparser which reads the tissue in the ini file
tissue_instance = tissueparser.read_tissue_config(TISSUE_CONFIG_FILE)
# setup batch with tissue instances which should be simulated
batch = factory.create_batch_to_simulate()
# tell the batch whcih tissue it should create samples from
batch.set_tissue_instance(tissue_instance)
# we want to simulate ten tissue instances in this example
df = batch.create_tissue_samples(NR_ELEMENTS_IN_BATCH)
# lets have a look at the dataframe. Each row corresponds to one tissue instance,
# each tissue instance is defined by various layers, which all have certain parameters
# like e.g. oxygenation (here sao2)
df
Out[3]:
In [4]:
# 1.3 create spectra - run simulation
# add reflectance column to dataframe
for w in WAVELENGTHS:
df["reflectances", w] = np.NAN # the reflectances have not been calculated yet, thus set no nan
# small helper function to create mco filenames:
def _create_mco_filename_for(prefix, batch, simulation):
return str(prefix) + "_Bat_" + str(batch) + "_Sim_" + str(simulation) + "_"
for i in range(df.shape[0]):
# set the desired element in the dataframe to be simulated
base_mco_filename = _create_mco_filename_for(BATCH_PREFIX,
BATCH_NR,
i)
tissue_model.set_base_mco_filename(base_mco_filename)
tissue_model.set_tissue_instance(df.loc[i, :])
tissue_model.update_mci_file(WAVELENGTHS)
# Run simulations for the created mci file
sim_wrapper.run_simulation()
# get information from created mco files
for i in range(df.shape[0]):
for wavelength in WAVELENGTHS:
# for simulation get which mco file was created
simulation_path = os.path.split(sim_wrapper.mcml_executable)[0]
base_mco_filename = _create_mco_filename_for(BATCH_PREFIX,
BATCH_NR,
i)
mco_filename = base_mco_filename + str(wavelength) + '.mco'
# get diffuse reflectance from simulation
df["reflectances", wavelength][i] = \
sim.get_diffuse_reflectance(os.path.join(simulation_path, mco_filename))
# delete created mco file
os.remove(os.path.join(simulation_path, mco_filename))
# Hooray, finished,
# now our dataframe also contains reflectances for each tissue instance:
df["reflectances"]
Out[4]:
In [5]:
# 2.1 analyse spectra - plot reflectances
# the usual settings for plotting in ipython notebooks
import matplotlib.pylab as plt
%matplotlib inline
# let's have a look at our reflectances
df["reflectances"].T.plot(kind="line")
plt.ylabel("reflectance")
plt.xlabel("wavelengths [m]")
# put legend outside of plot
plt.gca().legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid()
In [6]:
# 2.1 analyse spectra - show distribution of blood volume fraction (vhb) and sao2
# now we need some special pandas functions
import pandas as pd
# we're interested in the distribution of vhb and sao2 in the first layer (layer0)
df_vhb_sao2 = df["layer0"][["vhb", "sao2"]]
# plot a scatter matrix showing the distribution of vhb and sao2.
# of course, with this little data this does not really make sense,
# however it is a useful tool for analysis if much data is available
pd.tools.plotting.scatter_matrix(df_vhb_sao2, alpha=0.75, figsize=(6, 6))
plt.show()
In [7]:
# 3.1 manipulate spectra - apply sliding average
# in 3.1 and 3.2 we will adapt the generated spectra to an imaginary imaging system
# This system has filters with 20nm bandwith (taken care of in 3.1)
# and takes multispectral images in 10nm steps (taken care of in 3.2)
# the module mc.dfmanipulations was written to provide some basic,
# often needed manipulations of the calculated spectra
# all dmfmanipulations are performed inplace, however, the df is also returned.
import mc.dfmanipulations as dfmani
# first copy to not loose our original data
df2 = df.copy()
# e.g. we cann apply a sliding average to our data. This is usefull if
# we want to see e.g. how the reflectance was recorded by bands with a certain width
# a sliding average of 11 will take the five left and five right of the current reflectance
# and average. Because we take 2nm steps of reflectance in our simulation this means
# a 20nm window.
dfmani.fold_by_sliding_average(df2, 11)
# lets again plot the reflectances
df2["reflectances"].T.plot(kind="line")
plt.ylabel("reflectance")
plt.xlabel("wavelengths [m]")
# put legend outside of plot
plt.gca().legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid()
# we can see that the bump at 560nm is "smoother"
In [8]:
# 3.2 manipulate spectra - select certain wavelenghts
# our imaginary imaging system takes images in 10nm steps from 470 to 660nm
imaging_system_wavelengths = np.arange(470, 670, 10) * 10**-9
df3 = df2.copy()
dfmani.interpolate_wavelengths(df3, imaging_system_wavelengths)
# let's look at the newly created reflectances
df3["reflectances"].T.plot(kind="line", marker='o')
plt.ylabel("reflectance")
plt.xlabel("wavelengths [m]")
# put legend outside of plot
plt.gca().legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid()
In [8]:
# that's it, folks, if you want, you can save the created dataframe easily to csv:
df.to_csv("results.csv", index=False)
In [ ]: