diff --git a/Modules/Biophotonics/python/iMC/mc/mcmanipulations.py b/Modules/Biophotonics/python/iMC/mc/dfmanipulations.py similarity index 100% rename from Modules/Biophotonics/python/iMC/mc/mcmanipulations.py rename to Modules/Biophotonics/python/iMC/mc/dfmanipulations.py diff --git a/Modules/Biophotonics/python/iMC/mc/test/test_dfmanipulations.py b/Modules/Biophotonics/python/iMC/mc/test/test_dfmanipulations.py new file mode 100644 index 0000000000..7e62886813 --- /dev/null +++ b/Modules/Biophotonics/python/iMC/mc/test/test_dfmanipulations.py @@ -0,0 +1,78 @@ +''' +Created on Oct 19, 2015 + +@author: wirkert +''' +import unittest + +import numpy as np +from pandas.util.testing import assert_frame_equal + +from mc.batches import ColonMuscleBatch +import mc.dfmanipulations as dfmani + +class Test(unittest.TestCase): + + def setUp(self): + # create a colon batch with 2 samples + self.test_batch = ColonMuscleBatch() + self.test_batch.create_parameters(2) + + # artificially add 10 fake "reflectances" to this batch + # at 10 fake "wavelengths" + WAVELENGHTS = np.linspace(450, 720, 10) + reflectance1 = np.arange(0, 30, 3) + reflectance2 = np.arange(30, 60, 3) + for w in WAVELENGHTS: + self.test_batch.df["reflectances", w] = np.NAN + for r1, r2, w in zip(reflectance1, reflectance2, WAVELENGHTS): + self.test_batch.df["reflectances", w][0] = r1 + self.test_batch.df["reflectances", w][1] = r2 + + # shortcut to dataframe that we are interested in: + self.df = self.test_batch.df + + def test_sliding_average(self): + # by test design folding should not alter elements (only at boundaries, + # which are excluded by array slicing: + expected_elements = self.df.reflectances.iloc[:, 1:-1].copy() + dfmani.fold_by_sliding_average(self.df, 3) + + assert_frame_equal(self.df.reflectances, expected_elements) + + def test_interpolation(self): + new_wavelengths = [465, 615, 555] + + dfmani.interpolate_wavelengths(self.df, new_wavelengths) + + expected = np.array([[1.5, 16.5, 10.5], [31.5, 46.5, 40.5]]) + np.testing.assert_almost_equal(self.df.reflectances.as_matrix(), + expected, + err_msg="test if interpolation " + + "works fine on batches") + + def test_select_n(self): + """ this is less a test and more a showing of how to select n elements + from a dataframe.""" + # draw one sample. Look into documentation for sample to see all the + # options. Sample is quite powerfull. + self.df = self.df.sample(1) + self.assertEqual(self.df.shape[0], 1, + "one sample selected") + + def test_sortout_bands(self): + """ this is less a test and more a showing of how to sortout specific + bands from a dataframe """ + # drop the 510 and 720 nm band + band_names_to_sortout = [510, 720] + self.df.drop(band_names_to_sortout, axis=1, level=1, inplace=True) + + df_r = self.df["reflectances"] + self.assertTrue(not (510 in df_r.columns)) + self.assertTrue(not 720 in df_r.columns) + self.assertTrue(690 in df_r.columns) + + +if __name__ == "__main__": + # import sys;sys.argv = ['', 'Test.testName'] + unittest.main() diff --git a/Modules/Biophotonics/python/iMC/mc/test/test_mcmanipulations.py b/Modules/Biophotonics/python/iMC/mc/test/test_mcmanipulations.py deleted file mode 100644 index 2311302501..0000000000 --- a/Modules/Biophotonics/python/iMC/mc/test/test_mcmanipulations.py +++ /dev/null @@ -1,58 +0,0 @@ -''' -Created on Oct 19, 2015 - -@author: wirkert -''' -import unittest - -import numpy as np - -from mc.batches import AbstractBatch -import mc.mcmanipulations as mcmani - -class Test(unittest.TestCase): - - def setUp(self): - self.test_batch = AbstractBatch() - a = np.arange(0, 30, 3) - b = np.arange(30, 60, 3) - self.test_batch.reflectances = np.vstack((a, b)) - layer1 = np.arange(0, 10, 1) - layer1 = layer1[:, np.newaxis] - layer2 = np.arange(0, 10, 1) - layer2 = layer2[:, np.newaxis] - layer2 = np.vstack((layer2, layer2)) - self.test_batch.layers = [layer1, layer2] - self.test_batch._wavelengths = np.arange(0, 10, 1) - - def test_sliding_average(self): - # by design folding should not alter elements (only at boundaries, - # which are excluded by array slicing: - expected_elements = self.test_batch.reflectances[:, 1:-2] - mcmani.fold_by_sliding_average(self.test_batch, 3) - np.testing.assert_almost_equal(self.test_batch.reflectances[:, 1:-2], - expected_elements, - err_msg="test if sliding average " + - "function works on batches") - - def test_interpolation(self): - mcmani.interpolate_wavelengths(self.test_batch, - np.array([0.5, 5.5, 3.5])) - expected = np.array([[1.5, 16.5, 10.5], [31.5, 46.5, 40.5]]) - np.testing.assert_almost_equal(self.test_batch.reflectances, expected, - err_msg="test if interpolation " + - "works fine on batches") - - def test_select_n(self): - mcmani.select_n(self.test_batch, 1) - self.assertEqual(self.test_batch.reflectances.shape[0], 1, - "1 reflectances selected") - self.assertEqual(self.test_batch.layers[0].shape[0], 1, - "1 elements selected from layer 1") - self.assertEqual(self.test_batch.layers[1].shape[0], 1, - "1 elements selected from layer 2") - - -if __name__ == "__main__": - # import sys;sys.argv = ['', 'Test.testName'] - unittest.main() diff --git a/Modules/Biophotonics/python/iMC/mc/tissuemodels.py b/Modules/Biophotonics/python/iMC/mc/tissuemodels.py index 0a6714d5d9..504d310db9 100644 --- a/Modules/Biophotonics/python/iMC/mc/tissuemodels.py +++ b/Modules/Biophotonics/python/iMC/mc/tissuemodels.py @@ -1,132 +1,132 @@ ''' Created on Sep 9, 2015 @author: wirkert ''' import numpy as np from mc.sim import MciWrapper from mc.usuag import Ua, UsgJacques class AbstractTissue(object): ''' Initializes a abstract tissue model" ''' def set_nr_photons(self, nr_photons): self._mci_wrapper.set_nr_photons(nr_photons) def set_mci_filename(self, mci_filename): self._mci_wrapper.set_mci_filename(mci_filename) def set_mco_filename(self, mco_filename): self._mci_wrapper.set_mco_filename(mco_filename) def set_wavelength(self, wavelength): self.wavelength = wavelength def create_mci_file(self): # set layers for i, ua in enumerate(self.uas): self._mci_wrapper.set_layer(i, # layer nr self.ns[i], # refraction index self.uas[i](self.wavelength), # ua self.usgs[i](self.wavelength)[0], # us self.usgs[i](self.wavelength)[1], # g self.ds[i]) # d # now that the layers have been updated: create file self._mci_wrapper.create_mci_file() def __str__(self): """ Overwrite this method! print the current model""" model_string = \ "" return model_string def __init__(self, ns, uas, usgs, ds): self._mci_wrapper = MciWrapper() self.uas = uas self.usgs = usgs self.ds = ds self.ns = ns # initially create layers. these will be overwritten as soon # as create_mci_file is called. for i in enumerate(uas): self._mci_wrapper.add_layer() class GenericTissue(AbstractTissue): ''' Initializes a n-layer generic tissue model ''' - def set_batch_element(self, df, element): + def set_dataframe_element(self, df, element): """take the element element of the batch and set the tissue to resemble the structure specified by this""" layers = [l for l in df.columns.levels[0] if "layer" in l] for i, l in enumerate(layers): self.set_layer(i, df[l, "vhb"][element], df[l, "sao2"][element], df[l, "a_mie"][element], df[l, "b_mie"][element], df[l, "d"][element], df[l, "n"][element], df[l, "g"][element]) def set_layer(self, layer_nr=0, bvf=None, saO2=None, a_mie=None, b_mie=None, d=None, n=None, g=None): """Helper function to set one layer.""" if bvf is None: bvf = 0.02 if saO2 is None: saO2 = 0.7 if a_mie is None: a_mie = 10. * 100 if d is None: d = 500. * 10 ** -6 if b_mie is None: b_mie = 1.286 if n is None: n = 1.38 if g is None: g = 0. # build obejct for absorption coefficient determination self.uas[layer_nr].bvf = bvf self.uas[layer_nr].saO2 = saO2 # and one for scattering coefficient self.usgs[layer_nr].a_mie = a_mie self.usgs[layer_nr].a_ray = 0. self.usgs[layer_nr].b_mie = b_mie self.usgs[layer_nr].g = g self.ds[layer_nr] = d self.ns[layer_nr] = n def __str__(self): """print the current model""" model_string = "" for i, ua in enumerate(self.uas): layer_string = "layer " + str(i) + \ " - vhb: " + "%.2f" % self.uas[i].bvf + \ "%; sao2: " + "%.2f" % self.uas[i].saO2 + \ "%; a_mie: " + "%.2f" % (self.usgs[i].a_mie / 100.) + \ "cm^-1; a_ray: " + "%.2f" % (self.usgs[i].a_ray / 100.) + \ "cm^-1; b_mie: " + "%.3f" % self.usgs[i].b_mie + \ "; d: " + "%.0f" % (self.ds[i] * 10 ** 6) + "um" + \ "; n: " + "%.2f" % (self.ns[i]) + \ "; g: " + "%.2f" % (self.usgs[i].g) + "\n" model_string += layer_string return model_string def __init__(self): uas = [Ua(), Ua(), Ua()] usgs = [UsgJacques(), UsgJacques(), UsgJacques()] ds = np.ones_like(uas, dtype=float) * 500.*10 ** -6 ns = np.ones_like(uas, dtype=float) * 1.38 super(GenericTissue, self).__init__(ns, uas, usgs, ds) diff --git a/Modules/Biophotonics/python/iMC/regression/preprocessing.py b/Modules/Biophotonics/python/iMC/regression/preprocessing.py index 6a27afedb4..6b3a01901d 100644 --- a/Modules/Biophotonics/python/iMC/regression/preprocessing.py +++ b/Modules/Biophotonics/python/iMC/regression/preprocessing.py @@ -1,74 +1,67 @@ ''' Created on Oct 26, 2015 @author: wirkert ''' -import copy - import numpy as np from sklearn.preprocessing import Normalizer -import mc.mcmanipulations as mcmani -def preprocess2(batch, nr_samples=None, w_percent=None, magnification=None, +def preprocess2(df, nr_samples=None, w_percent=None, magnification=None, bands_to_sortout=None): - working_batch = copy.deepcopy(batch) - mcmani.sortout_bands(working_batch, bands_to_sortout) - # get reflectance and oxygenation - X = working_batch.reflectances - y = working_batch.layers[0][:, [0, 1]] + df.drop(bands_to_sortout, axis=1, level=1, inplace=True) + # first set 0 reflectances to nan + df["reflectances"] = df["reflectances"].replace(to_replace=0., + value=np.nan) # remove nan - no_nan = ~np.isnan(X).any(axis=1) - X = X[no_nan, :] - y = y[no_nan, :] - # remove zero - no_zero = ~(X <= 0).any(axis=1) - X = X[no_zero, :] - y = y[no_zero, :] + df.dropna(inplace=True) # extract nr_samples samples from data if nr_samples is not None: - X = X[0:nr_samples] - y = y[0:nr_samples] + df = df.sample(nr_samples) + + # get reflectance and oxygenation + X = df.reflectances + y = df.layer0[["sao2", "vhb"]] # do data magnification if magnification is not None: X_temp = X y_temp = y for i in range(magnification - 1): X = np.vstack((X, X_temp)) y = np.vstack((y, y_temp)) # add noise to reflectances if w_percent is not None and (w_percent > 0.): noises = np.random.normal(loc=0., scale=w_percent, size=X.shape) X += noises * X X = np.clip(X, 0.00001, 1.) # do normalizations X = normalize(X) - return X, np.squeeze(y) + return X, y def preprocess(batch, nr_samples=None, w_percent=None, magnification=None, bands_to_sortout=None): X, y = preprocess2(batch, nr_samples, w_percent, magnification, bands_to_sortout) - return X, y[:, 1] + return X, y["sao2"] def normalize(X): # normalize reflectances normalizer = Normalizer(norm='l1') X = normalizer.transform(X) # reflectances to absorption absorptions = -np.log(X) X = absorptions # get rid of sorted out bands normalizer = Normalizer(norm='l2') X = normalizer.transform(X) return X diff --git a/Modules/Biophotonics/python/iMC/script_calculate_spectra.py b/Modules/Biophotonics/python/iMC/script_calculate_spectra.py index a4f161fce6..e8a9c310c5 100644 --- a/Modules/Biophotonics/python/iMC/script_calculate_spectra.py +++ b/Modules/Biophotonics/python/iMC/script_calculate_spectra.py @@ -1,137 +1,136 @@ ''' Created on Sep 9, 2015 @author: wirkert ''' import logging import datetime import os import time import numpy as np import luigi import scriptpaths as sp import mc.factories as mcfac from mc.sim import SimWrapper, get_diffuse_reflectance # parameter setting NR_BATCHES = 100 NR_ELEMENTS_IN_BATCH = 1000 # the wavelengths to be simulated WAVELENGHTS = np.arange(450, 720, 2) * 10 ** -9 NR_PHOTONS = 10 ** 6 # general output path config sp.ROOT_FOLDER = \ - "/media/wirkert/data/Data/2015_11_12_IPCAI_in_silico" -sp.RESULTS_FOLDER = "mc_data_after_revision" + "/media/wirkert/data/Data/2016_02_22_IPCAI_revision_in_silico" +sp.RESULTS_FOLDER = "mc_data" # experiment configuration MCI_FILENAME = "./temp.mci" MCO_FILENAME = "temp.mco" # this path definitly needs to be adapted by you PATH_TO_MCML = "/home/wirkert/workspace/monteCarlo/gpumcml/fast-gpumcml/" EXEC_MCML = "gpumcml.sm_20" class CreateSpectraTask(luigi.Task): - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() batch_nr = luigi.IntParameter() nr_samples = luigi.IntParameter() factory = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, - self.batch_prefix + "_" + + self.df_prefix + "_" + str(self.batch_nr) + ".txt")) def run(self): start = time.time() # setup simulation wrapper self.sim_wrapper = SimWrapper() self.sim_wrapper.set_mci_filename(MCI_FILENAME) self.sim_wrapper.set_mcml_executable(PATH_TO_MCML + EXEC_MCML) # setup model self.tissue_model = self.factory.create_tissue_model() self.tissue_model.set_mci_filename(self.sim_wrapper.mci_filename) self.tissue_model.set_mco_filename(MCO_FILENAME) self.tissue_model.set_nr_photons(NR_PHOTONS) # setup array in which data shall be stored batch = self.factory.create_batch_to_simulate() batch.create_parameters(self.nr_samples) # dataframe created by batch: df = batch.df # add reflectance column to dataframe for w in WAVELENGHTS: df["reflectances", w] = np.NAN # for each instance of our tissue model for i in range(df.shape[0]): - self.tissue_model.set_batch_element(df, i) + self.tissue_model.set_dataframe_element(df, i) logging.info("running simulation " + str(i) + " for\n" + str(self.tissue_model)) start = time.time() # map the _wavelengths array to reflectance list reflectances = map(self.wavelength_to_reflectance, WAVELENGHTS) end = time.time() # store in dataframe for r, w in zip(reflectances, WAVELENGHTS): df["reflectances", w][i] = r # success! logging.info("successfully ran simulation in " + "{:.2f}".format(end - start) + " seconds") # clean up temporarily created files os.remove(MCI_FILENAME) created_mco_file = PATH_TO_MCML + "/" + MCO_FILENAME if os.path.isfile(created_mco_file): os.remove(created_mco_file) # save the created output f = open(self.output().path, 'w') df.to_csv(f) end = time.time() logging.info("time for creating batch of mc data: %.f s" % (end - start)) def wavelength_to_reflectance(self, wavelength): """helper function to determine the reflectance for a given wavelength using the current model and simulation """ self.tissue_model.set_wavelength(wavelength) self.tissue_model.create_mci_file() if os.path.isfile(PATH_TO_MCML + EXEC_MCML): self.sim_wrapper.run_simulation() return get_diffuse_reflectance(PATH_TO_MCML + MCO_FILENAME) else: raise IOError("path to gpumcml not valid") - if __name__ == '__main__': - logging.basicConfig(filename='calculate_spectra' + + logging.basicConfig(filename='log/calculate_spectra_' + str(datetime.datetime.now()) + - '.log', level=logging.INFO) + '.log', level=logging.INFO) ch = logging.StreamHandler() ch.setLevel(logging.INFO) logger = logging.getLogger() logger.addHandler(ch) luigi.interface.setup_interface_logging() sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) BATCH_NUMBERS = np.arange(0, NR_BATCHES, 1) for i in BATCH_NUMBERS: - colon_task = CreateSpectraTask("ipcai_revision_colon", i, + colon_task = CreateSpectraTask("ipcai_revision_colon_train", i, NR_ELEMENTS_IN_BATCH, mcfac.ColonMuscleMcFactory()) generic_task = CreateSpectraTask("ipcai_revision_generic", i, NR_ELEMENTS_IN_BATCH, mcfac.GenericMcFactory()) w.add(colon_task) w.add(generic_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_plot_one_spectrum.py b/Modules/Biophotonics/python/iMC/script_plot_one_spectrum.py index 0558b261fb..6ff6481992 100644 --- a/Modules/Biophotonics/python/iMC/script_plot_one_spectrum.py +++ b/Modules/Biophotonics/python/iMC/script_plot_one_spectrum.py @@ -1,70 +1,70 @@ ''' Created on Oct 12, 2015 @author: wirkert ''' import pickle import logging import numpy as np import matplotlib.pyplot as plt import luigi import tasks_regression as rt from msi.plot import plot from msi.msi import Msi import msi.normalize as norm import scriptpaths as sp sp.ROOT_FOLDER = "/media/wirkert/data/Data/2015_xxxx_plot_one_spectrum" # the wavelengths recorded by our camera RECORDED_WAVELENGTHS = \ np.arange(450, 720, 2) * 10 ** -9 PARAMS = np.array([0.1, # bvf 0.7, # SaO2 0.0, # billirubin 500., # a_mie 0.0, # a_ray 1.091, # b (for scattering 500. * 10 ** -6]) # d_muc class PlotOneSpectrum(luigi.Task): - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() def requires(self): - return rt.TrainForestForwardModel(self.batch_prefix) + return rt.TrainForestForwardModel(self.df_prefix) def run(self): f = file(self.input().path, "r") rf = pickle.load(f) f.close() refl = rf.predict(PARAMS) msi = Msi(refl) msi.set_wavelengths(RECORDED_WAVELENGTHS) # norm.standard_normalizer.normalize(msi) plot(msi) plt.gca().set_xlabel("wavelength") plt.gca().set_ylabel("normalized reflectance") plt.grid() plt.ylim([0.0, 0.4]) plt.title("bvf: " + str(PARAMS[0]) + "; saO2: " + str(PARAMS[1]) + "; bili: " + str(PARAMS[2]) + "; a_mie: " + str(PARAMS[3]) + "; a_ray: " + str(PARAMS[4]) + "; d_muc: " + str(PARAMS[6])) plt.show() if __name__ == '__main__': logging.basicConfig(level=logging.INFO) luigi.interface.setup_interface_logging() sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) main_task = PlotOneSpectrum(batch_prefix= "jacques_no_billi_generic_scattering_") w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/setup.py b/Modules/Biophotonics/python/iMC/setup.py index 1468a1eff5..d4e414a683 100644 --- a/Modules/Biophotonics/python/iMC/setup.py +++ b/Modules/Biophotonics/python/iMC/setup.py @@ -1,25 +1,25 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 7 18:41:50 2015 @author: wirkert """ from setuptools import setup, find_packages setup(name='MultiSpectralImaging', version='0.1', description='Multi spectral imaging utilities', author='Sebastian Wirkert', author_email='s.wirkert@dkfz-heidelberg.de', license='BSD', packages=find_packages(exclude=['test*']), package_dir={}, package_data={'data':['*.txt', '*.mci', '*.nrrd']}, install_requires=['numpy>=1.8.2', 'scikit-learn>0.16.1', 'scipy', - 'SimpleITK>=0.9.0', + 'SimpleITK>=0.9.0', 'subprocess32', 'scikit-image', 'pypng', 'pandas>0.17'], entry_points={} # for scripts, add later ) diff --git a/Modules/Biophotonics/python/iMC/tasks_analyze.py b/Modules/Biophotonics/python/iMC/tasks_analyze.py index a967088a5c..ab79f6baec 100644 --- a/Modules/Biophotonics/python/iMC/tasks_analyze.py +++ b/Modules/Biophotonics/python/iMC/tasks_analyze.py @@ -1,205 +1,205 @@ ''' Created on Aug 27, 2015 @author: wirkert ''' import os import numpy as np import Image import ImageEnhance import pickle import luigi import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import make_axes_locatable import SimpleITK as sitk from msi.io.nrrdreader import NrrdReader import msi.msimanipulations as msimani import tasks_preprocessing as ppt import tasks_regression as rt import scriptpaths as sp from regression.estimation import estimate_image def plot_axis(axis, image, title): im2 = axis.imshow(image, interpolation='nearest', alpha=1.0) - axis.set_title(title, fontsize=5) - divider = make_axes_locatable(axis) - cax = divider.append_axes("right", size="10%", pad=0.05) - cbar = plt.colorbar(im2, cax=cax) - cbar.ax.tick_params(labelsize=5) + # axis.set_title(title, fontsize=5) + # divider = make_axes_locatable(axis) + # cax = divider.append_axes("right", size="10%", pad=0.05) + # cbar = plt.colorbar(im2, cax=cax) + # cbar.ax.tick_params(labelsize=5) axis.xaxis.set_visible(False) axis.yaxis.set_visible(False) - return im2, cbar + return im2 # , cbar class EstimateTissueParametersTask(luigi.Task): imageName = luigi.Parameter() - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() def requires(self): return ppt.PreprocessMSI(imageName=self.imageName), \ rt.TrainForest(batch_prefix= - self.batch_prefix) + self.df_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, "processed", sp.FINALS_FOLDER, self.imageName + "_" + - self.batch_prefix + + self.df_prefix + "estimate.nrrd")) def run(self): reader = NrrdReader() msi = reader.read(self.input()[0].path) # load random forest e_file = open(self.input()[1].path, 'r') e = pickle.load(e_file) sitk_image = estimate_image(msi, e) outFile = self.output().open('w') outFile.close() sitk.WriteImage(sitk_image, self.output().path) class ReprojectReflectancesTask(luigi.Task): imageName = luigi.Parameter() - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() def requires(self): return EstimateTissueParametersTask(imageName=self.imageName, batch_prefix= - self.batch_prefix), \ + self.df_prefix), \ rt.TrainForestForwardModel(batch_prefix= - self.batch_prefix) + self.df_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, "processed", sp.FINALS_FOLDER, self.imageName + "_backprojection_" + - self.batch_prefix + + self.df_prefix + "estimate.nrrd")) def run(self): reader = NrrdReader() msi = reader.read(self.input()[0].path) # load random forest e_file = open(self.input()[1].path, 'r') e = pickle.load(e_file) sitk_image = estimate_image(msi, e) outFile = self.output().open('w') outFile.close() sitk.WriteImage(sitk_image, self.output().path) def r2_evaluation(reflectances, backprojections): # our own version of r2_score. Cannot use the sklearn # one since this calculates one mean r2_score. # we want the r2_score for each element in the image. mean = np.mean(reflectances, axis=2) # unfortunately keepdims does not work for np.mean, thus: new_shape_hack = mean.shape + (1,) mean = np.reshape(mean, new_shape_hack) ss_res = np.sum(np.square(reflectances - backprojections), axis=-1) ss_tot = np.sum(np.square(reflectances - mean), axis=-1) r2_image = 1 - ss_res / ss_tot r2_median = np.median(r2_image) # santiy check: one element equals the "professional" version. # assert(r2_score(reflectances[0, 0, :], backprojections[0, 0, :]) # == r2_image[0, 0]) return r2_median, r2_image class CreateNiceParametricImagesTask(luigi.Task): imageName = luigi.Parameter() - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() def requires(self): return EstimateTissueParametersTask(imageName=self.imageName, batch_prefix= - self.batch_prefix), \ + self.df_prefix), \ ppt.CorrectImagingSetupTask(imageName=self.imageName), \ ppt.SegmentMSI(imageName=self.imageName), \ ppt.PreprocessMSI(imageName=self.imageName), \ ReprojectReflectancesTask(imageName=self.imageName, batch_prefix= - self.batch_prefix) + self.df_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, self.imageName + "_" + - self.batch_prefix + + self.df_prefix + "_summary.png")) def run(self): # read the segmentation reader = sitk.ImageFileReader() reader.SetFileName(self.input()[2].path) segmentation = reader.Execute() segmentation = sitk.GetArrayFromImage(segmentation) # read the multispectral image nrrdReader = NrrdReader() msi = nrrdReader.read(self.input()[1].path) msimani.apply_segmentation(msi, segmentation) # and the normalized multispectra image norm_msi = nrrdReader.read(self.input()[3].path) msimani.apply_segmentation(norm_msi, segmentation) # read in the parametric image param_image = nrrdReader.read(self.input()[0].path) msimani.apply_segmentation(param_image, segmentation) image = param_image.get_image() # and the reprojection reprojection = nrrdReader.read(self.input()[4].path) msimani.apply_segmentation(reprojection, segmentation) f, axarr = plt.subplots(2, 3) # create artificial rgb rgb_image = msi.get_image()[:, :, [7, 5, 1]] rgb_image /= np.max(rgb_image) rgb_image *= 255. rgb_image = rgb_image.astype(np.uint8) im = Image.fromarray(rgb_image, 'RGB') enh_brightness = ImageEnhance.Brightness(im) im = enh_brightness.enhance(5.) plot_image = np.array(im) top_left_axis = axarr[0, 0] top_left_axis.imshow(plot_image, interpolation='nearest') top_left_axis.set_title("synthetic rgb (segmented pixels black)", fontsize=5) top_left_axis.xaxis.set_visible(False) top_left_axis.yaxis.set_visible(False) plt.set_cmap("jet") # plot error map r2_median, r2_image = r2_evaluation(norm_msi.get_image(), reprojection.get_image()) bottom_left_axis = axarr[1, 0] r2_image_smaller_zero = r2_image < 0 r2_image[r2_image_smaller_zero] = 0 plot_axis(bottom_left_axis, r2_image, "r2 score (cut off at 0).\n" + "median r2_score = " + str(r2_median)) # plot parametric maps plot_axis(axarr[0, 1], image[:, :, 0], "blood volume fraction [%]") image[0, 0, 1] = 0.; image[1, 0, 1] = 1. plot_axis(axarr[0, 2], image[:, :, 1], "oxygenation [%]") # image[0, 0, 2] = 175 * 10 ** -6; image[1, 0, 0] = 750 * 10 ** -6 plot_axis(axarr[1, 1], image[:, :, 2] * 100, "bilirubin concentration [mg/dl]") plot_axis(axarr[1, 2], image[:, :, 3] / 100., "mie scatter parameter [1/cm]") plt.savefig(self.output().path, dpi=1000, bbox_inches='tight') diff --git a/Modules/Biophotonics/python/iMC/tasks_mc.py b/Modules/Biophotonics/python/iMC/tasks_mc.py index e54d60bfb5..97f51f41dd 100644 --- a/Modules/Biophotonics/python/iMC/tasks_mc.py +++ b/Modules/Biophotonics/python/iMC/tasks_mc.py @@ -1,139 +1,119 @@ ''' Created on Sep 10, 2015 @author: wirkert ''' -''' -Created on Aug 31, 2015 - -@author: wirkert -''' - import os -import pickle + +import pandas as pd import luigi import scriptpaths as sp -import mc.mcmanipulations as mcmani +import mc.dfmanipulations as dfmani from msi.io.spectrometerreader import SpectrometerReader from msi.io.msiwriter import MsiWriter from msi.io.msireader import MsiReader from msi.normalize import NormalizeMean import msi.msimanipulations as msimani - - - class SpectrometerFile(luigi.Task): input_file = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.SPECTROMETER_FOLDER, self.input_file)) class FilterTransmission(luigi.Task): input_file = luigi.Parameter() def requires(self): return SpectrometerFile(self.input_file) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "processed_transmission" + self.input_file + ".msi")) def run(self): reader = SpectrometerReader() filter_transmission = reader.read(self.input().path) # filter high and low _wavelengths wavelengths = filter_transmission.get_wavelengths() fi_image = filter_transmission.get_image() fi_image[wavelengths < 450 * 10 ** -9] = 0.0 fi_image[wavelengths > 720 * 10 ** -9] = 0.0 # filter elements farther away than +- 30nm name_to_float = float(os.path.splitext(self.input_file)[0]) fi_image[wavelengths < (name_to_float - 30) * 10 ** -9] = 0.0 fi_image[wavelengths > (name_to_float + 30) * 10 ** -9] = 0.0 # elements < 0 are set to 0. fi_image[fi_image < 0.0] = 0.0 # write it writer = MsiWriter(filter_transmission) writer.write(self.output().path) class JoinBatches(luigi.Task): - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, - self.batch_prefix + "_" + - "all" + ".imc")) + self.df_prefix + "_" + + "all" + ".txt")) def run(self): path = os.path.join(sp.ROOT_FOLDER, sp.MC_DATA_FOLDER) # get all files in the search path files = [ f for f in os.listdir(path) \ if os.path.isfile(os.path.join(path, f)) ] # from these get only those who start with correct batch prefix - batch_file_names = \ - [ f for f in files if f.startswith(self.batch_prefix)] - batch_files = \ - [open(os.path.join(path, f), 'r') for f in batch_file_names] + df_file_names = \ + [ os.path.join(path, f) for f in files if f.startswith(self.df_prefix)] # load these files - batches = [pickle.load(f) for f in batch_files] + dfs = [pd.read_csv(f, header=[0, 1]) for f in df_file_names] # now join them to one batch - joined_batch = reduce(join_batches, batches) + joined_df = pd.concat(dfs, ignore_index=True) # write it - joined_batch_file = open(self.output().path, 'w') - # there seems to be a bug in pickle. thus, unfortunately the _generator - # has to be removed before saving - joined_batch._generator = None - pickle.dump(joined_batch, joined_batch_file) + joined_df.to_csv(self.output().path, index=False) class CameraBatch(luigi.Task): """takes a batch of reference data and converts it to the spectra processed by the camera""" - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() def requires(self): - return JoinBatches(self.batch_prefix) + return JoinBatches(self.df_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, - self.batch_prefix + "_all_camera.imc")) + self.df_prefix + "_all_camera.txt")) def run(self): - f = file(self.input().path, "r") - batch = pickle.load(f) - f.close() + # load dataframe + df = pd.read_csv(self.input().path, header=[0, 1]) # camera batch creation: - camera_batch = batch - mcmani.fold_by_sliding_average(camera_batch, 6) - mcmani.interpolate_wavelengths(camera_batch, - sp.RECORDED_WAVELENGTHS) + dfmani.fold_by_sliding_average(df, 6) + dfmani.interpolate_wavelengths(df, sp.RECORDED_WAVELENGTHS) # write it - joined_batch_file = open(self.output().path, 'w') - pickle.dump(camera_batch, joined_batch_file) - - + df.to_csv(self.output().path, index=False) def get_transmission_data(input_path, desired_wavelengths): """small helper method to get filter transmission from file at input_path. The _wavelengths are interpolated to the desired ones""" reader = MsiReader() filter_transmission = reader.read(input_path) msimani.interpolate_wavelengths(filter_transmission, desired_wavelengths) # normalize to one normalizer = NormalizeMean() normalizer.normalize(filter_transmission) return filter_transmission