diff --git a/Modules/Biophotonics/python/iMC/mc/batches.py b/Modules/Biophotonics/python/iMC/mc/batches.py index 8d50fd3e79..0a8b26a0f0 100644 --- a/Modules/Biophotonics/python/iMC/mc/batches.py +++ b/Modules/Biophotonics/python/iMC/mc/batches.py @@ -1,267 +1,227 @@ ''' Created on Oct 15, 2015 @author: wirkert ''' -import copy - import numpy as np - +from pandas import DataFrame +import pandas as pd class AbstractBatch(object): """summarizes a batch of simulated mc spectra""" def __init__(self): - self.reflectances = np.array([]) - self.layers = [] # list of np arrays, one for each layer - self.nr_photons = 10 ** 6 - self.wavelengths = np.arange(450, 720, 2) * 10 ** -9 - # by default random uniform samples will be generated - self.generator = np.random.random_sample + self._nr_layers = 0 # internally keeps track of number of layers + my_index = pd.MultiIndex(levels=[[], []], + labels=[[], []]) + self.df = DataFrame(columns=my_index) def create_parameters(self, nr_samples): """create the parameters for the batch, the simulation has to create the resulting reflectances""" pass def nr_elements(self): - return self.layers[0].shape[0] - - def join(self, other_batch): - """helper method to join two batches""" - joined_batch = copy.deepcopy(self) + return self.df.shape[0] - f = lambda x, y: np.append(x, y, axis=0) - joined_batch.layers = map(f, self.layers, other_batch.layers) - joined_batch.reflectances = f(joined_batch.reflectances, - other_batch.reflectances) - # next two should be the same for two batches from the same experiment - np.testing.assert_almost_equal(joined_batch.wavelengths, - other_batch.wavelengths) - np.testing.assert_almost_equal(joined_batch.nr_photons, - other_batch.nr_photons) - return joined_batch class GenericBatch(AbstractBatch): - """generic three layer batch """ + """generic 3-layer batch with each layer having the same oxygenation """ def __init__(self): super(GenericBatch, self).__init__() def append_one_layer(self, saO2, nr_samples): """helper function to create parameters for one layer""" # scales data to lie between maxi and mini instead of 0 and 1 scale = lambda x, mini, maxi: x * (maxi - mini) + mini - - samples = self.generator(size=(nr_samples, 7)) - # scale samples to correct ranges - samples[:, 0] = scale(samples[:, 0], 0., 1.00) - # saO2 is the same for each layer, since the blood "diffuses" in tissue - samples[:, 1] = saO2 # saO2 - samples[:, 2] = scale(samples[:, 2], 5.* 100., 30.* 100.) # a_mie - samples[:, 3] = scale(samples[:, 3], 0.* 100., 60.* 100.) * 0 # a_ray - # d will be normalized later to 2mm total depth - samples[:, 4] = scale(samples[:, 4], 0.0, 1.) # d - samples[:, 5] = scale(samples[:, 5], 1.33, 1.54) # n - samples[:, 6] = scale(samples[:, 6], 0.80, 0.95) # g - # append as last layer - self.layers.append(samples) + # shortcut to random generator + gen = np.random.random_sample + gen_n = np.random.normal + + # create layer elements + self.df["layer" + str(self._nr_layers), "vhb"] = \ + scale(gen(nr_samples), 0, 1.) + self.df["layer" + str(self._nr_layers), "sao2"] = \ + saO2 + self.df["layer" + str(self._nr_layers), "a_mie"] = \ + np.clip(gen_n(loc=18.9, scale=10.2, size=nr_samples), + 0.1, np.inf) * 100 # to 1/m + self.df["layer" + str(self._nr_layers), "b_mie"] = \ + np.clip(gen_n(loc=1.286, scale=0.521, size=nr_samples), 0, np.inf) + self.df["layer" + str(self._nr_layers), "d"] = \ + scale(gen(nr_samples), 0, 1.) + self.df["layer" + str(self._nr_layers), "n"] = \ + scale(gen(nr_samples), 1.33, 1.54) + self.df["layer" + str(self._nr_layers), "g"] = \ + scale(gen(nr_samples), 0.8, 0.95) + self._nr_layers += 1 def create_parameters(self, nr_samples): """Create generic three layer batch with a total diameter of 2mm. saO2 is the same in all layers, but all other parameters vary randomly within each layer""" - saO2 = self.generator(size=nr_samples) + saO2 = np.random.random_sample(size=nr_samples) + # create three layers with random samples self.append_one_layer(saO2, nr_samples) self.append_one_layer(saO2, nr_samples) self.append_one_layer(saO2, nr_samples) - # create empty reflectances matrix - self.reflectances = np.zeros((nr_samples, len(self.wavelengths))) - # set all weights to 1 - self.weights = np.ones(nr_samples, dtype=float) - # "normalize" d to 2mm - ds = map(lambda x: x[:, 4], self.layers) - total_d = reduce(lambda x, y: x + y, ds) - desired_d = 2000. * 10 ** -6 - for l in self.layers: - l[:, 4] = l[:, 4] / total_d * desired_d - l[l[:, 4] < 25 * 10 ** -6, :] = 25 * 10 ** -6 + # "normalize" d to 2mm + # first extract all layers from df + self.df -class GaussianBatch(GenericBatch): - """three layer batch with gaussian sampled tissue parameters""" - - def __init__(self): - super(GaussianBatch, self).__init__() - - def append_one_layer(self, saO2, nr_samples): - """helper function to create parameters for one layer""" - # scales data to lie between maxi and mini instead of 0 and 1 - scale = lambda x, mini, maxi: x * (maxi - mini) + mini - # rescale gaussian - scale_gaussian = lambda x, u, sigma: x * sigma + u - samples = self.generator(size=(nr_samples, 7)) - # scale samples to correct ranges - # bvf - samples[:, 0] = np.random.normal(size=nr_samples) - samples[:, 0] = np.clip(scale_gaussian(samples[:, 0], 0.04, 0.02), - 0., - 1.) - # saO2 is the same for each layer, since the blood "diffuses" in tissue - samples[:, 1] = scale(saO2, 0., 1.) - # a_mie - samples[:, 2] = scale(samples[:, 2], 5. *100, 30.* 100) - # a_ray - samples[:, 3] = scale(samples[:, 3], 0.*100., 60.*100.) * 0. - # d will be normalized later to 2mm total depth - samples[:, 4] = scale(samples[:, 4], 0., 1.) - samples[:, 5] = scale(samples[:, 5], 1.33, 1.54) # n - samples[:, 6] = scale(samples[:, 6], 0.80, 0.95) # g - # append as last layer - self.layers.append(samples) + layers = [l for l in self.df.columns.levels[0] if "layer" in l] + # summarize all ds + sum_d = 0 + for l in layers: + sum_d += self.df[l, "d"] + for l in layers: + self.df[l, "d"] = self.df[l, "d"] / sum_d * 2000. * 10 ** -6 + self.df[l, "d"] = np.clip(self.df[l, "d"], 25 * 10 ** -6, np.inf) class LessGenericBatch(AbstractBatch): """less generic three layer batch. This only varies blood volume fraction w.r.t. the ColonMuscleBatch. Let's see if DA works in this case.""" def __init__(self): super(LessGenericBatch, self).__init__() + def append_one_layer(self, saO2, n, d_ranges, nr_samples): """helper function to create parameters for one layer""" + # scales data to lie between maxi and mini instead of 0 and 1 scale = lambda x, mini, maxi: x * (maxi - mini) + mini - samples = self.generator(size=(nr_samples, 7)) - # scale samples to correct ranges - # bvf - samples[:, 0] = scale(samples[:, 0], 0., 1.) - # saO2 is the same for each layer, since the blood "diffuses" in tissue - samples[:, 1] = saO2 # saO2 - samples[:, 2] = scale(samples[:, 2], 5.* 100., 30.* 100.) # a_mie - samples[:, 3] = scale(samples[:, 3], 0.* 100., 60.* 100.) * 0. # a_ray - # d - samples[:, 4] = scale(samples[:, 4], d_ranges[0], d_ranges[1]) # d - samples[:, 5] = n - samples[:, 6] = scale(samples[:, 6], 0.80, 0.95) # g - # append as last layer - self.layers.append(samples) + # shortcut to random generator + gen = np.random.random_sample + + # create as generic batch + super(LessGenericBatch, self).append_one_layer(saO2, nr_samples) + self._nr_layers -= 1 # we're not finished + + # but some changes in specific layer elements + # more specific layer thicknesses + self.df["layer" + str(self._nr_layers), "d"] = \ + scale(gen(nr_samples), d_ranges[0], d_ranges[1]) + # more specific n + self.df["layer" + str(self._nr_layers), "n"] = \ + n + + self._nr_layers += 1 def create_parameters(self, nr_samples): """Create generic three layer batch with a total diameter of 2mm. saO2 is the same in all layers, but all other parameters vary randomly within each layer""" - saO2 = self.generator(size=nr_samples) + saO2 = np.random.random_sample(size=nr_samples) n = np.ones_like(saO2) # create three layers with random samples # muscle - self.append_one_layer(saO2, n * 1.36, (600.*10 ** -6, 1010.*10 ** -6), nr_samples) + self.append_one_layer(saO2, n * 1.36, (600.*10 ** -6, 1010.*10 ** -6), + nr_samples) # submucosa - self.append_one_layer(saO2, n * 1.36, (415.*10 ** -6, 847.*10 ** -6), nr_samples) + self.append_one_layer(saO2, n * 1.36, (415.*10 ** -6, 847.*10 ** -6), + nr_samples) # mucosa - self.append_one_layer(saO2, n * 1.38, (395.*10 ** -6, 603.*10 ** -6), nr_samples) - # create empty reflectances matrix - self.reflectances = np.zeros((nr_samples, len(self.wavelengths))) - # set all weights to 1 - self.weights = np.ones(nr_samples, dtype=float) + self.append_one_layer(saO2, n * 1.38, (395.*10 ** -6, 603.*10 ** -6), + nr_samples) -class ColonMuscleBatch(AbstractBatch): + +class ColonMuscleBatch(GenericBatch): """three layer batch simulating colonic tissue""" def __init__(self): super(ColonMuscleBatch, self).__init__() def append_one_layer(self, saO2, n, d_ranges, nr_samples): """helper function to create parameters for one layer""" + # scales data to lie between maxi and mini instead of 0 and 1 scale = lambda x, mini, maxi: x * (maxi - mini) + mini - samples = self.generator(size=(nr_samples, 7)) - # scale samples to correct ranges - # bvf - samples[:, 0] = scale(samples[:, 0], 0., 0.1) - # saO2 is the same for each layer, since the blood "diffuses" in tissue - samples[:, 1] = saO2 # saO2 - samples[:, 2] = scale(samples[:, 2], 5.* 100., 30.* 100.) # a_mie - samples[:, 3] = scale(samples[:, 3], 0.* 100., 60.* 100.) * 0. # a_ray - # d - samples[:, 4] = scale(samples[:, 4], d_ranges[0], d_ranges[1]) # d - samples[:, 5] = n - samples[:, 6] = scale(samples[:, 6], 0.80, 0.95) # g - # append as last layer - self.layers.append(samples) + # shortcut to random generator + gen = np.random.random_sample + + # create as generic batch + super(ColonMuscleBatch, self).append_one_layer(saO2, nr_samples) + self._nr_layers -= 1 # we're not finished + + # but some changes in specific layer elements + # less blood + self.df["layer" + str(self._nr_layers), "vhb"] = \ + scale(gen(nr_samples), 0, 0.1) + # more specific layer thicknesses + self.df["layer" + str(self._nr_layers), "d"] = \ + scale(gen(nr_samples), d_ranges[0], d_ranges[1]) + # more specific n + self.df["layer" + str(self._nr_layers), "n"] = \ + n + + self._nr_layers += 1 def create_parameters(self, nr_samples): """Create generic three layer batch with a total diameter of 2mm. saO2 is the same in all layers, but all other parameters vary randomly within each layer""" - saO2 = self.generator(size=nr_samples) + saO2 = np.random.random_sample(size=nr_samples) n = np.ones_like(saO2) # create three layers with random samples # muscle - self.append_one_layer(saO2, n * 1.36, (600.*10 ** -6, 1010.*10 ** -6), nr_samples) - # submucosa - self.append_one_layer(saO2, n * 1.36, (415.*10 ** -6, 847.*10 ** -6), nr_samples) - # mucosa - self.append_one_layer(saO2, n * 1.38, (395.*10 ** -6, 603.*10 ** -6), nr_samples) - # create empty reflectances matrix - self.reflectances = np.zeros((nr_samples, len(self.wavelengths))) - # set all weights to 1 - self.weights = np.ones(nr_samples, dtype=float) - - -class VisualizationBatch(AbstractBatch): - """batch used for visualization of different spectra. Feel free to adapt - for your visualization purposes.""" - - def __init__(self): - super(VisualizationBatch, self).__init__() - # self.wavelengths = np.arange(470, 680, 10) * 10 ** -9 - - def append_one_layer(self, bvf, saO2, a_mie, a_ray, d, n, g, nr_samples): - """helper function to create parameters for one layer""" - samples = np.zeros((nr_samples, 7)) - # scale samples to correct ranges - samples[:, 0] = bvf - samples[:, 1] = saO2 - samples[:, 2] = a_mie - samples[:, 3] = a_ray - # d will be normalized later to 2mm total depth - samples[:, 4] = d - samples[:, 5] = n - samples[:, 6] = g - # append as last layer - self.layers.append(samples) - - def create_parameters(self, nr_samples): - # bvf = np.linspace(0.0, .1, nr_samples) - # saO2 = np.linspace(0., 1., nr_samples) - # d = np.linspace(175, 735, nr_samples) * 10 ** -6 - # a_mie = np.linspace(5., 30., nr_samples) * 100 - # a_ray = np.linspace(0., 60., nr_samples) * 100 - # n = np.linspace(1.33, 1.54, nr_samples) - # g = np.linspace(0, 0.95, nr_samples) - # create three layers with random samples - self.append_one_layer(0.02, 0.1, 30.*100., 0., 500 * 10 ** -6, - 1.38, 0.9, + self.append_one_layer(saO2, n * 1.36, (600.*10 ** -6, 1010.*10 ** -6), nr_samples) - self.append_one_layer(0.04, 0.7, 5.0 * 100, 0.*100, 500 * 10 ** -6, - 1.36, 0.9, + # submucosa + self.append_one_layer(saO2, n * 1.36, (415.*10 ** -6, 847.*10 ** -6), nr_samples) - self.append_one_layer(0.04, 0.7, 5.0 * 100, 0.*100, 500 * 10 ** -6, - 1.36, 0.9, + # mucosa + self.append_one_layer(saO2, n * 1.38, (395.*10 ** -6, 603.*10 ** -6), nr_samples) - # create empty reflectances matrix - self.reflectances = np.zeros((nr_samples, len(self.wavelengths))) - # set all weights to 1 - self.weights = np.ones(nr_samples, dtype=float) - -def join_batches(batch_1, batch_2): - return batch_1.join(batch_2) +# class VisualizationBatch(AbstractBatch): +# """batch used for visualization of different spectra. Feel free to adapt +# for your visualization purposes.""" +# +# def __init__(self): +# super(VisualizationBatch, self).__init__() +# # self._wavelengths = np.arange(470, 680, 10) * 10 ** -9 +# +# def append_one_layer(self, bvf, saO2, a_mie, a_ray, d, n, g, nr_samples): +# """helper function to create parameters for one layer""" +# samples = np.zeros((nr_samples, 7)) +# # scale samples to correct ranges +# samples[:, 0] = bvf +# samples[:, 1] = saO2 +# samples[:, 2] = a_mie +# samples[:, 3] = a_ray +# # d will be normalized later to 2mm total depth +# samples[:, 4] = d +# samples[:, 5] = n +# samples[:, 6] = g +# # append as last layer +# self.layers.append(samples) +# +# def create_parameters(self, nr_samples): +# # bvf = np.linspace(0.0, .1, nr_samples) +# # saO2 = np.linspace(0., 1., nr_samples) +# # d = np.linspace(175, 735, nr_samples) * 10 ** -6 +# # a_mie = np.linspace(5., 30., nr_samples) * 100 +# # a_ray = np.linspace(0., 60., nr_samples) * 100 +# # n = np.linspace(1.33, 1.54, nr_samples) +# # g = np.linspace(0, 0.95, nr_samples) +# # create three layers with random samples +# self.append_one_layer(0.02, 0.1, 30.*100., 0., 500 * 10 ** -6, +# 1.38, 0.9, +# nr_samples) +# self.append_one_layer(0.04, 0.7, 5.0 * 100, 0.*100, 500 * 10 ** -6, +# 1.36, 0.9, +# nr_samples) +# self.append_one_layer(0.04, 0.7, 5.0 * 100, 0.*100, 500 * 10 ** -6, +# 1.36, 0.9, +# nr_samples) diff --git a/Modules/Biophotonics/python/iMC/mc/factories.py b/Modules/Biophotonics/python/iMC/mc/factories.py index ee914879cc..3e91ca64d6 100644 --- a/Modules/Biophotonics/python/iMC/mc/factories.py +++ b/Modules/Biophotonics/python/iMC/mc/factories.py @@ -1,85 +1,74 @@ ''' Created on Oct 15, 2015 @author: wirkert ''' from mc.tissuemodels import AbstractTissue, GenericTissue -from mc.batches import AbstractBatch, GenericBatch, VisualizationBatch, \ - GaussianBatch, ColonMuscleBatch, LessGenericBatch +from mc.batches import AbstractBatch, GenericBatch, ColonMuscleBatch class AbstractMcFactory(object): ''' Monte Carlo Factory. Will create fitting models and batches, dependent on your task ''' def create_tissue_model(self): return AbstractTissue() def create_batch_to_simulate(self): return AbstractBatch() def __init__(self): ''' Constructor ''' class GenericMcFactory(AbstractMcFactory): def create_tissue_model(self): return GenericTissue() def create_batch_to_simulate(self): return GenericBatch() def __init__(self): ''' Constructor ''' -class GaussianMcFactory(GenericMcFactory): - - def create_batch_to_simulate(self): - return GaussianBatch() - - def __init__(self): - ''' - Constructor - ''' - class LessGenericMcFactory(GenericMcFactory): def create_batch_to_simulate(self): return LessGenericBatch() def __init__(self): ''' Constructor ''' class ColonMuscleMcFactory(GenericMcFactory): def create_batch_to_simulate(self): return ColonMuscleBatch() def __init__(self): ''' Constructor ''' class VisualizationMcFactory(AbstractMcFactory): def create_tissue_model(self): return GenericTissue() def create_batch_to_simulate(self): return VisualizationBatch() def __init__(self): ''' Constructor ''' diff --git a/Modules/Biophotonics/python/iMC/mc/mcmanipulations.py b/Modules/Biophotonics/python/iMC/mc/mcmanipulations.py index 4e5a9346d0..a97696e8ab 100644 --- a/Modules/Biophotonics/python/iMC/mc/mcmanipulations.py +++ b/Modules/Biophotonics/python/iMC/mc/mcmanipulations.py @@ -1,52 +1,52 @@ ''' Created on Oct 19, 2015 @author: wirkert ''' import numpy as np from scipy.interpolate import interp1d def fold_by_sliding_average(batch, window_size): """take a batch and apply a sliding average with given window size to the reflectances. window_size is elements to the left and to the right. There will be some boundary effect on the edges.""" sliding_average = lambda x: np.convolve(x, np.ones(window_size, dtype=float) / float(window_size), mode="same") batch.reflectances = np.apply_along_axis(sliding_average, axis=1, arr=batch.reflectances) def interpolate_wavelengths(batch, new_wavelengths): """ interpolate image data to fit new_wavelengths. Current implementation - performs simple linear interpolation. Neither existing nor new wavelengths + performs simple linear interpolation. Neither existing nor new _wavelengths need to be sorted. """ - interpolator = interp1d(batch.wavelengths, + interpolator = interp1d(batch._wavelengths, batch.reflectances, assume_sorted=False, bounds_error=False) batch.reflectances = interpolator(new_wavelengths) - batch.wavelengths = new_wavelengths + batch._wavelengths = new_wavelengths def sortout_bands(batch, bands_to_sortout=None): """ TODO: Documentation and test """ if bands_to_sortout is not None: # get rid of sorted out bands batch.reflectances = np.delete(batch.reflectances, bands_to_sortout, axis=1) - batch.wavelengths = np.delete(batch.wavelengths, + batch._wavelengths = np.delete(batch._wavelengths, bands_to_sortout) def select_n(batch, n): """ randomly select n elements from batch. TODO: Test """ perms = np.random.permutation(batch.reflectances.shape[0]) first_n_perms = perms[0:n] batch.reflectances = batch.reflectances[first_n_perms, :] for i, l in enumerate(batch.layers): batch.layers[i] = np.atleast_2d(l)[first_n_perms, :] return batch diff --git a/Modules/Biophotonics/python/iMC/mc/plot.py b/Modules/Biophotonics/python/iMC/mc/plot.py index a5f3236e32..2706a15c2b 100644 --- a/Modules/Biophotonics/python/iMC/mc/plot.py +++ b/Modules/Biophotonics/python/iMC/mc/plot.py @@ -1,43 +1,43 @@ ''' Created on Oct 16, 2015 @author: wirkert ''' import numpy as np import matplotlib.pyplot as plt class PlotFunctor(object): """helping functor necessary because we need to save color for plotting""" def __init__(self, axes, wavelengths, nr_plot_elements): self.axes = axes self.sortedIndices = sorted(range(len(wavelengths)), key=lambda k: wavelengths[k]) self.sortedWavelenghts = wavelengths[self.sortedIndices] self.nr_plot_elements = nr_plot_elements self.i = 0 def __call__(self, r): pass # set color so it slowly moves from blue to red plt_color = (1. / float(self.nr_plot_elements) * self.i, 0., 1. - (1. / float(self.nr_plot_elements) * self.i)) self.axes.plot(self.sortedWavelenghts, r[self.sortedIndices], "-o", color=plt_color) self.i += 1 return self.i def plot(batch, axes=None): if axes is None: axes = plt.gca() - f = PlotFunctor(axes, batch.wavelengths, batch.reflectances.shape[0]) + f = PlotFunctor(axes, batch._wavelengths, batch.reflectances.shape[0]) np.apply_along_axis(f, axis=1, arr=batch.reflectances) diff --git a/Modules/Biophotonics/python/iMC/mc/test/test_mcmanipulations.py b/Modules/Biophotonics/python/iMC/mc/test/test_mcmanipulations.py index 0222a8436c..2311302501 100644 --- a/Modules/Biophotonics/python/iMC/mc/test/test_mcmanipulations.py +++ b/Modules/Biophotonics/python/iMC/mc/test/test_mcmanipulations.py @@ -1,58 +1,58 @@ ''' 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) + 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 09f1f36fd2..0a6714d5d9 100644 --- a/Modules/Biophotonics/python/iMC/mc/tissuemodels.py +++ b/Modules/Biophotonics/python/iMC/mc/tissuemodels.py @@ -1,128 +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 three layer generic tissue model + Initializes a n-layer generic tissue model ''' - def set_batch_element(self, batch, element): + def set_batch_element(self, df, element): """take the element element of the batch and set the tissue to resemble the structure specified by this""" - self._mci_wrapper.set_nr_photons(batch.nr_photons) - for i, l in enumerate(batch.layers): + layers = [l for l in df.columns.levels[0] if "layer" in l] + for i, l in enumerate(layers): self.set_layer(i, - l[element, 0], # bvf - l[element, 1], # sao2 - l[element, 2], # a_mie - l[element, 3], # a_ray - l[element, 4], # d - l[element, 5], # n - l[element, 6]) # g + 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, a_ray=None, d=None, + 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 a_ray is None: - a_ray = 0. + 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 = a_ray + 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) + \ - " - bvf: " + "%.2f" % self.uas[i].bvf + \ - "%; SaO2: " + "%.2f" % self.uas[i].saO2 + \ + " - 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/script_analyze_ipcai_in_vivo_liver.py b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_liver.py index 50b926ce63..768542708c 100644 --- a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_liver.py +++ b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_liver.py @@ -1,355 +1,355 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import Image import ImageEnhance import pickle import time import numpy as np import luigi import matplotlib.pyplot as plt import SimpleITK as sitk from sklearn.ensemble.forest import RandomForestRegressor from scipy.ndimage.filters import gaussian_filter import tasks_analyze as at import tasks_preprocessing as pre import tasks_mc import scriptpaths as sp from msi.io.pngreader import PngReader from msi.io.nrrdreader import NrrdReader from msi.io.nrrdwriter import NrrdWriter import msi.msimanipulations as msimani import msi.imgmani as imgmani import msi.normalize as norm from regression.estimation import estimate_image from regression.linear import LinearSaO2Unmixing from regression.preprocessing import preprocess2, normalize from msi.io.tiffreader import TiffReader sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2015_11_12_IPCAI_in_vivo" sp.DATA_FOLDER = "liver_images" sp.FINALS_FOLDER = "liver_oxy" sp.FLAT_FOLDER = "flatfields_liver" sp.MC_DATA_FOLDER = "mc_data2" # sp.RECORDED_WAVELENGTHS = np.arange(470, 680, 10) * 10 ** -9 new_order = [0, 1, 2, 3, 4, 5, 6, 7] def resort_image_wavelengths(collapsed_image): return collapsed_image[:, new_order] def resort_wavelengths(msi): - """Neil organized his wavelengths differently. + """Neil organized his _wavelengths differently. This function fixes this.""" collapsed_image = imgmani.collapse_image(msi.get_image()) collapsed_image = resort_image_wavelengths(collapsed_image) msi.set_image(np.reshape(collapsed_image, msi.get_image().shape)) -# rebind this method since in this recoding the wavelengths got messed up +# rebind this method since in this recoding the _wavelengths got messed up sp.resort_wavelengths = resort_wavelengths sp.bands_to_sortout = np.array([]) # [0, 1, 2, 20, 19, 18, 17, 16]) class IPCAICreateOxyImageTask(luigi.Task): image_prefix = luigi.Parameter() batch_prefix = luigi.Parameter() def requires(self): return IPCAIEstimateTissueParametersTask(image_prefix=self.image_prefix, batch_prefix= self.batch_prefix), \ MultiSpectralImageFromTiffFiles(image_prefix=self.image_prefix), \ IPCAICorrectImagingSetupTask(image_prefix=self.image_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, self.image_prefix + "_" + self.batch_prefix + "_summary.png")) def run(self): nrrdReader = NrrdReader() # read in the parametric image param_image = nrrdReader.read(self.input()[0].path) image = param_image.get_image() # read the multispectral image msi = nrrdReader.read(self.input()[1].path) msi_original = nrrdReader.read(self.input()[2].path) plt.figure() # filter dark spots segmentation = np.max(msi.get_image(), axis=-1) < 100. # filter bright spots segmentation = np.logical_or(segmentation, (np.max(msi.get_image(), axis=-1) > 2000.)) # create artificial rgb msimani.apply_segmentation(msi_original, ~segmentation) rgb_image = msi_original.get_image()[:, :, [2, 3, 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(1.) plot_image = np.array(im) top_left_axis = plt.gca() top_left_axis.imshow(plot_image, interpolation='nearest') top_left_axis.set_title("synthetic rgb (specular pixels black)", fontsize=5) top_left_axis.xaxis.set_visible(False) top_left_axis.yaxis.set_visible(False) plt.savefig(self.output().path, dpi=1000, bbox_inches='tight') plt.set_cmap("jet") plt.figure() # plot parametric maps # oxy_image = image[:, :, 1] # image[500:1000, 1500:2500, 1] oxy_image = np.ma.masked_array(image[:, :, 1], segmentation) oxy_image[np.isnan(oxy_image)] = 0. oxy_image[np.isinf(oxy_image)] = 0. oxy_image[oxy_image > 1.] = 1. oxy_image[oxy_image < 0.] = 0. oxy_mean = np.mean(oxy_image) oxy_image[0, 0] = 0. oxy_image[0, 1] = 1. at.plot_axis(plt.gca(), oxy_image[:, :] * 100, "oxygenation [%]. Mean " + "{0:.1f}".format((oxy_mean * 100.)) + "%") plt.savefig(self.output().path + "_oxy.png", dpi=1000, bbox_inches='tight') bvf_image = np.ma.masked_array(image[:, :, 0], segmentation) bvf_image[np.isnan(bvf_image)] = 0. bvf_image[np.isinf(bvf_image)] = 0. bvf_image[bvf_image > 1.] = 1. bvf_image[bvf_image < 0.] = 0. bvf_mean = np.mean(bvf_image) bvf_image[0, 0] = 0. bvf_image[1, 1] = 0.1 at.plot_axis(plt.gca(), bvf_image[:, :] * 100, "blood volume fraction [%]. Mean " + "{0:.1f}".format((bvf_mean * 100.)) + "%") plt.savefig(self.output().path + "_bvf.png", dpi=1000, bbox_inches='tight') fd = open(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, 'results_summary.csv'), 'a') fd.write(self.image_prefix + ", " + str(int(oxy_mean * 100)) + "," + "{0:.1f}".format((bvf_mean * 100.)) + "\n") fd.close() class IPCAIEstimateTissueParametersTask(luigi.Task): image_prefix = luigi.Parameter() batch_prefix = luigi.Parameter() def requires(self): return IPCAIPreprocessMSI(image_prefix=self.image_prefix), \ IPCAITrainRegressor(batch_prefix= self.batch_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, "processed", sp.FINALS_FOLDER, self.image_prefix + "_" + self.batch_prefix + "estimate.nrrd")) def run(self): # load data reader = NrrdReader() msi = reader.read(self.input()[0].path) e_file = open(self.input()[1].path, 'r') e = pickle.load(e_file) # estimate start = time.time() sitk_image = estimate_image(msi, e) end = time.time() print "time necessary for estimating image parameters: ", \ str(end - start), "s" # save data outFile = self.output().open('w') outFile.close() sitk.WriteImage(sitk_image, self.output().path) class IPCAITrainRegressor(luigi.Task): batch_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "reg_" + self.batch_prefix)) def requires(self): return tasks_mc.CameraBatch(self.batch_prefix) def run(self): # extract data from the batch f = open(self.input().path, "r") batch_train = pickle.load(f) - # TODO: fix hack by sorting always to ascending wavelengths + # TODO: fix hack by sorting always to ascending _wavelengths batch_train.reflectances = resort_image_wavelengths( batch_train.reflectances) - batch_train.wavelengths = batch_train.wavelengths[new_order] + batch_train._wavelengths = batch_train._wavelengths[new_order] X, y = preprocess2(batch_train, w_percent=0.1, bands_to_sortout=sp.bands_to_sortout) # train regressor reg = RandomForestRegressor(max_depth=10, n_estimators=10, n_jobs=-1) # reg = LinearSaO2Unmixing() reg.fit(X, y) # reg = LinearSaO2Unmixing() # save regressor f = self.output().open('w') pickle.dump(reg, f) f.close() class IPCAIPreprocessMSI(luigi.Task): image_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.image_prefix + "_preprocessed.nrrd")) def requires(self): return IPCAICorrectImagingSetupTask(image_prefix=self.image_prefix) def run(self): reader = NrrdReader() msi = reader.read(self.input().path) # sortout unwanted bands msi.set_image(imgmani.sortout_bands(msi.get_image(), sp.bands_to_sortout)) # zero values would lead to infinity logarithm, thus clip. msi.set_image(np.clip(msi.get_image(), 0.00001, 2. ** 64)) # normalize to get rid of lighting intensity norm.standard_normalizer.normalize(msi) # transform to absorption msi.set_image(-np.log(msi.get_image())) # normalize by l2 for stability norm.standard_normalizer.normalize(msi, "l2") pre.touch_and_save_msi(msi, self.output()) class IPCAICorrectImagingSetupTask(luigi.Task): """corrects for dark image and flatfield""" image_prefix = luigi.Parameter() def requires(self): return MultiSpectralImageFromTiffFiles(image_prefix=self.image_prefix), \ FlatfieldFromPNGFiles() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.image_prefix + "_image_setup_corrected.nrrd")) def run(self): - """sort wavelengths and normalize by respective integration times""" + """sort _wavelengths and normalize by respective integration times""" # read inputs reader = NrrdReader() msi = reader.read(self.input()[0].path) flatfield = reader.read(self.input()[1].path) # correct image by flatfield msimani.flatfield_correction(msi, flatfield) pre.touch_and_save_msi(msi, self.output()) class FlatfieldFromPNGFiles(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "flatfield" + ".nrrd")) def run(self): reader = PngReader() msi = reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "Flatfield9cm_I000x_F")) writer = NrrdWriter(msi) msi.set_image(msi.get_image().astype(float)) seg_reader = NrrdReader() segmentation = seg_reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "segmentation.nrrd")) # msimani.apply_segmentation(msi, segmentation) # msimani.calculate_mean_spectrum(msi) # apply gaussian smoothing for error reduction img = gaussian_filter(msi.get_image(), sigma=(5, 5, 0), order=0) msi.set_image(img) writer.write(self.output().path) class MultiSpectralImageFromTiffFiles(luigi.Task): image_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.image_prefix + ".nrrd")) def run(self): reader = TiffReader() msi = reader.read(os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, self.image_prefix)) writer = NrrdWriter(msi) writer.write(self.output().path) if __name__ == '__main__': # root folder there the data lies luigi.interface.setup_interface_logging() sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) # run over all subfolders (non-recursively) # to collect the data and generate the results image_file_folder = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER) onlyfiles = [ f for f in os.listdir(image_file_folder) if os.path.isfile(os.path.join(image_file_folder, f)) ] # each complete stack has elements F0 - F7 only_colon_images = [ f for f in onlyfiles if f.startswith("1_") ] files_prefixes = [ f[:-28] for f in only_colon_images] for f in files_prefixes: main_task = IPCAICreateOxyImageTask( image_prefix=f, batch_prefix= "ipcai_colon_muscle_train") w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_small_bowel.py b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_small_bowel.py index 7309121468..4e894cae80 100644 --- a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_small_bowel.py +++ b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_small_bowel.py @@ -1,297 +1,297 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import Image import ImageEnhance import pickle import logging import datetime import numpy as np import luigi import matplotlib.pyplot as plt import SimpleITK as sitk from sklearn.ensemble.forest import RandomForestRegressor from scipy.ndimage.filters import gaussian_filter import tasks_analyze as at import tasks_mc import scriptpaths as sp from msi.io.pngreader import PngReader from msi.io.nrrdreader import NrrdReader from msi.io.nrrdwriter import NrrdWriter import msi.msimanipulations as msimani import msi.imgmani as imgmani import msi.normalize as norm from regression.estimation import estimate_image from regression.linear import LinearSaO2Unmixing from regression.preprocessing import preprocess, preprocess2 from msi.io.tiffringreader import TiffRingReader from sklearn.linear_model.base import LinearRegression sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2015_11_12_IPCAI_in_vivo" sp.DATA_FOLDER = "small_bowel_images" sp.FINALS_FOLDER = "small_bowel_images_colorbar20_80" sp.FLAT_FOLDER = "flatfields_liver" # sp.RECORDED_WAVELENGTHS = np.arange(470, 680, 10) * 10 ** -9 new_order = [0, 1, 2, 3, 4, 5, 6, 7] def resort_image_wavelengths(collapsed_image): return collapsed_image[:, new_order] def resort_wavelengths(msi): - """Neil organized his wavelengths differently. + """Neil organized his _wavelengths differently. This function fixes this.""" collapsed_image = imgmani.collapse_image(msi.get_image()) collapsed_image = resort_image_wavelengths(collapsed_image) msi.set_image(np.reshape(collapsed_image, msi.get_image().shape)) -# rebind this method since in this recoding the wavelengths got messed up +# rebind this method since in this recoding the _wavelengths got messed up sp.resort_wavelengths = resort_wavelengths sp.bands_to_sortout = np.array([]) # [0, 1, 2, 20, 19, 18, 17, 16]) class IPCAICreateOxyImageTask(luigi.Task): image_name = luigi.Parameter() batch_prefix = luigi.Parameter() def requires(self): return IPCAITrainRegressor(batch_prefix=self.batch_prefix), \ FlatfieldFromPNGFiles() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, self.image_name + "_" + self.batch_prefix + "_summary.png.png")) def run(self): nrrd_reader = NrrdReader() tiff_ring_reader = TiffRingReader() # read the flatfield flat = nrrd_reader.read(self.input()[1].path) # read the msi nr_filters = len(sp.RECORDED_WAVELENGTHS) full_image_name = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, self.image_name) msi, segmentation = tiff_ring_reader.read(full_image_name, nr_filters) segmentation = np.logical_and(segmentation, (np.max(msi.get_image(), axis=-1) < 1000.)) # read the segmentation # # seg_path = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, # "segmentation_" + self.image_name + ".nrrd") # segmentation = nrrd_reader.read(seg_path) # segmentation = segmentation.get_image() # read the regressor e_file = open(self.input()[0].path, 'r') e = pickle.load(e_file) # correct image setup filter_nr = int(self.image_name[-6:-5]) original_order = np.arange(nr_filters) new_image_order = np.concatenate(( original_order[nr_filters - filter_nr:], original_order[:nr_filters - filter_nr])) # resort msi to restore original order msimani.get_bands(msi, new_image_order) # correct by flatfield msimani.flatfield_correction(msi, flat) # create artificial rgb rgb_image = msi.get_image()[:, :, [2, 3, 1]] rgb_image /= np.max(rgb_image) rgb_image *= 255. # preprocess the image # sortout unwanted bands print "1" msi.set_image(imgmani.sortout_bands(msi.get_image(), sp.bands_to_sortout)) # zero values would lead to infinity logarithm, thus clip. msi.set_image(np.clip(msi.get_image(), 0.00001, 2. ** 64)) # normalize to get rid of lighting intensity norm.standard_normalizer.normalize(msi) # transform to absorption msi.set_image(-np.log(msi.get_image())) # normalize by l2 for stability norm.standard_normalizer.normalize(msi, "l2") print "2" # estimate sitk_image = estimate_image(msi, e) image = sitk.GetArrayFromImage(sitk_image) plt.figure() print "3" 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 = plt.gca() top_left_axis.imshow(plot_image, interpolation='nearest') top_left_axis.xaxis.set_visible(False) top_left_axis.yaxis.set_visible(False) plt.set_cmap("jet") print "4" # plot parametric maps oxy_image = image[:, :] segmentation[0, 0] = 1 segmentation[0, 1] = 1 oxy_image = np.ma.masked_array(image[:, :], (segmentation < 1)) oxy_image[np.isnan(oxy_image)] = 0. oxy_image[np.isinf(oxy_image)] = 0. oxy_mean = np.mean(oxy_image) oxy_image[0, 0] = 0.2 oxy_image[0, 1] = 0.8 oxy_image = np.clip(oxy_image, 0.2, 0.8) # oxy_image[oxy_image > 0.8] = 0.8 # oxy_image[oxy_image < 0.2] = 0.2 at.plot_axis(plt.gca(), oxy_image[:, :], "oxygenation [%]. Mean " + "{0:.1f}".format((oxy_mean * 100.)) + "%") plt.savefig(self.output().path, dpi=500, bbox_inches='tight') # plt.figure() # bvf_image = np.ma.masked_array(image[:, :, 0], (segmentation < 1)) # # bvf_image = image[:, :, 0] # bvf_image[np.isnan(bvf_image)] = 0. # bvf_image[np.isinf(bvf_image)] = 0. # bvf_image[bvf_image > 1.] = 1. # bvf_image[bvf_image < 0.] = 0. # bvf_mean = np.mean(bvf_image) # bvf_image[0, 0] = 0. # bvf_image[1, 1] = 0.1 # at.plot_axis(plt.gca(), bvf_image[:, :] * 100, # "blood volume fraction [%]. Mean " + "{0:.1f}".format((bvf_mean * 100.)) + "%") # # plt.savefig(self.output().path + "_bvf.png", dpi=500, bbox_inches='tight') print "4.5" # fd = open(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, # sp.FINALS_FOLDER, # 'results_summary.csv'), 'a') # fd.write(self.image_name + ", " + str(int(oxy_mean * 100)) + # ", " + "{0:.1f}".format((bvf_mean * 100.)) + "\n") # fd.close() print "5" plt.close("all") class IPCAITrainRegressor(luigi.Task): batch_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "reg_" + self.batch_prefix)) def requires(self): return tasks_mc.CameraBatch(self.batch_prefix) def run(self): # extract data from the batch f = open(self.input().path, "r") batch_train = pickle.load(f) - # TODO: fix hack by sorting always to ascending wavelengths + # TODO: fix hack by sorting always to ascending _wavelengths batch_train.reflectances = resort_image_wavelengths( batch_train.reflectances) - batch_train.wavelengths = batch_train.wavelengths[new_order] + batch_train._wavelengths = batch_train._wavelengths[new_order] X, y = preprocess2(batch_train, w_percent=0.1, bands_to_sortout=sp.bands_to_sortout) # train regressor reg = RandomForestRegressor(10, min_samples_leaf=10, max_depth=9, n_jobs=-1) # reg = KNeighborsRegressor(algorithm="auto") # reg = LinearRegression() # reg = sklearn.svm.SVR(kernel="rbf", degree=3, C=100., gamma=10.) # reg = LinearSaO2Unmixing() reg.fit(X, y[:, 1]) # reg = LinearSaO2Unmixing() # save regressor f = self.output().open('w') pickle.dump(reg, f) f.close() class FlatfieldFromPNGFiles(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "flatfield" + ".nrrd")) def run(self): reader = PngReader() msi = reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "Flatfield9cm_I000x_F")) writer = NrrdWriter(msi) msi.set_image(msi.get_image().astype(float)) seg_reader = NrrdReader() segmentation = seg_reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "segmentation.nrrd")) # msimani.apply_segmentation(msi, segmentation) # msimani.calculate_mean_spectrum(msi) # apply gaussian smoothing for error reduction img = gaussian_filter(msi.get_image(), sigma=(5, 5, 0), order=0) msi.set_image(img) writer.write(self.output().path) if __name__ == '__main__': # root folder there the data lies logging.basicConfig(filename='small_bowel' + str(datetime.datetime.now()) + '.log', level=logging.INFO) luigi.interface.setup_interface_logging() ch = logging.StreamHandler() ch.setLevel(logging.INFO) logger = logging.getLogger() logger.addHandler(ch) sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) # run over all subfolders (non-recursively) # to collect the data and generate the results image_file_folder = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER) onlyfiles = [ f for f in os.listdir(image_file_folder) if os.path.isfile(os.path.join(image_file_folder, f)) ] onlyfiles.sort() onlyfiles = [ f for f in onlyfiles if f.endswith(".tiff") ] for f in onlyfiles: main_task = IPCAICreateOxyImageTask( image_name=f, batch_prefix= "ipcai_colon_muscle_train") w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_calculate_spectra.py b/Modules/Biophotonics/python/iMC/script_calculate_spectra.py index 50716c6f63..a4f161fce6 100644 --- a/Modules/Biophotonics/python/iMC/script_calculate_spectra.py +++ b/Modules/Biophotonics/python/iMC/script_calculate_spectra.py @@ -1,40 +1,137 @@ ''' 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 tasks_mc 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" +sp.RESULTS_FOLDER = "mc_data_after_revision" + +# 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() + 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 + "_" + + 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) + 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='in_silico' + str(datetime.datetime.now()) + + logging.basicConfig(filename='calculate_spectra' + + str(datetime.datetime.now()) + '.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, 100, 1) + BATCH_NUMBERS = np.arange(0, NR_BATCHES, 1) for i in BATCH_NUMBERS: - main_task = tasks_mc.CreateSpectraTask("ipcai_less_generic", i, 1000, - mcfac.LessGenericMcFactory()) - w.add(main_task) + colon_task = CreateSpectraTask("ipcai_revision_colon", 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_visualize_spectra.py b/Modules/Biophotonics/python/iMC/script_visualize_spectra.py index 20a8ac8e78..ce209f65c3 100644 --- a/Modules/Biophotonics/python/iMC/script_visualize_spectra.py +++ b/Modules/Biophotonics/python/iMC/script_visualize_spectra.py @@ -1,92 +1,92 @@ ''' Created on Sep 10, 2015 @author: wirkert ''' import logging import os import pickle import numpy as np import matplotlib.pyplot as plt import matplotlib import luigi from regression.preprocessing import preprocess2, normalize import scriptpaths as sp import tasks_mc import mc.plot import mc.factories as mcfac # general output path config sp.ROOT_FOLDER = \ "/media/wirkert/data/Data/2015_06_01_Filtertransmittance_Spectrometer" sp.DATA_FOLDER = "spectrometer_measurements" sp.MC_DATA_FOLDER = "processed" SPECTRAL_PLOTS = "spectral_plots" sp.RECORDED_WAVELENGTHS = np.arange(470, 690, 10) * 10 ** -9 # matplotlib.rcParams.update({'font.size': 18}) class VisualizeSpectraTask(luigi.Task): batch_prefix = luigi.Parameter() batch_nr = luigi.IntParameter() nr_samples = luigi.IntParameter() def requires(self): return tasks_mc.CreateSpectraTask(self.batch_prefix, self.batch_nr, self.nr_samples, mcfac.VisualizationMcFactory()), \ tasks_mc.CameraBatch(self.batch_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, SPECTRAL_PLOTS, "plotted_spectra_" + self.batch_prefix + "_" + str(self.batch_nr) + ".png")) def run(self): f = file(self.input()[0].path, "r") batch = pickle.load(f) f = file(self.input()[1].path, "r") camera_batch = pickle.load(f) f.close() - batch.wavelengths *= 10 ** 9 + batch._wavelengths *= 10 ** 9 # X, y = preprocess2(camera_batch, w_percent=0.1) # camera_batch.reflectances = X # f, axarr = plt.subplots(2, 1) # axis to plot the original data c.f. the mcml simulation # axis for plotting the folded spectrum + normalization mc.plot.plot(batch) # tidy up and save plot plt.grid() plt.ylabel("r", fontsize=30) plt.xlabel("$\lambda$ [nm]", fontsize=30) # org_ax.axes.xaxis.set_ticklabels([]) # org_ax.xaxis.set_visible(False) plt.show() # plt.savefig(self.output().path, dpi=250) 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 = VisualizeSpectraTask("paper2", 0, 1) w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/setup.py b/Modules/Biophotonics/python/iMC/setup.py index 2a8ab46196..1468a1eff5 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', - 'scikit-image', 'pypng'], + 'scikit-image', 'pypng', 'pandas>0.17'], entry_points={} # for scripts, add later ) diff --git a/Modules/Biophotonics/python/iMC/tasks_mc.py b/Modules/Biophotonics/python/iMC/tasks_mc.py index 40948f5631..e54d60bfb5 100644 --- a/Modules/Biophotonics/python/iMC/tasks_mc.py +++ b/Modules/Biophotonics/python/iMC/tasks_mc.py @@ -1,218 +1,139 @@ ''' Created on Sep 10, 2015 @author: wirkert ''' ''' Created on Aug 31, 2015 @author: wirkert ''' import os -import numpy as np import pickle -import time -import logging import luigi import scriptpaths as sp -from mc.sim import SimWrapper, get_diffuse_reflectance -from mc.batches import join_batches import mc.mcmanipulations as mcmani 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 -# 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 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 + # 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() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.batch_prefix + "_" + "all" + ".imc")) 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] # load these files batches = [pickle.load(f) for f in batch_files] # now join them to one batch joined_batch = reduce(join_batches, batches) # write it joined_batch_file = open(self.output().path, 'w') - # there seems to be a bug in pickle. thus, unfortunately the generator + # there seems to be a bug in pickle. thus, unfortunately the _generator # has to be removed before saving - joined_batch.generator = None + joined_batch._generator = None pickle.dump(joined_batch, joined_batch_file) class CameraBatch(luigi.Task): """takes a batch of reference data and converts it to the spectra processed by the camera""" batch_prefix = luigi.Parameter() def requires(self): return JoinBatches(self.batch_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.batch_prefix + "_all_camera.imc")) def run(self): f = file(self.input().path, "r") batch = pickle.load(f) f.close() # camera batch creation: camera_batch = batch mcmani.fold_by_sliding_average(camera_batch, 6) mcmani.interpolate_wavelengths(camera_batch, sp.RECORDED_WAVELENGTHS) # write it joined_batch_file = open(self.output().path, 'w') pickle.dump(camera_batch, joined_batch_file) -class CreateSpectraTask(luigi.Task): - batch_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 + "_" + - str(self.batch_nr) + ".imc")) - - 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) - # setup array in which data shall be stored - batch = self.factory.create_batch_to_simulate() - batch.create_parameters(self.nr_samples) - - for i in range(batch.nr_elements()): - self.tissue_model.set_batch_element(batch, 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, - batch.wavelengths) - end = time.time() - # store in batch_nr - batch.reflectances[i, :] = np.asarray(reflectances) - # success! - logging.info("successfully ran simulation in " + - "{:.2f}".format(end - start) + " seconds") - # convert the lists in batch to np arrays - batch.reflectances = np.array(batch.reflectances) - # 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') - # there seems to be a bug in pickle. thus, unfortunately the generator - # has to be removed before saving - batch.generator = None - pickle.dump(batch, f) - f.close() - - 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") 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""" + 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