diff --git a/Modules/Biophotonics/python/iMC/mc/batches.py b/Modules/Biophotonics/python/iMC/mc/batches.py index 58693ae906..1058c6d906 100644 --- a/Modules/Biophotonics/python/iMC/mc/batches.py +++ b/Modules/Biophotonics/python/iMC/mc/batches.py @@ -1,245 +1,262 @@ ''' Created on Oct 15, 2015 @author: wirkert ''' 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._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.df.shape[0] class GenericBatch(AbstractBatch): """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 # 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 = 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) # "normalize" d to 2mm # first extract all layers from df self.df 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 # 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 = 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) 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 # 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 = 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) +class GenericMeanScatteringBatch(GenericBatch): + """three layer batch simulating colonic tissue""" + + def __init__(self): + super(GenericMeanScatteringBatch, self).__init__() + + def append_one_layer(self, saO2, nr_samples): + """helper function to create parameters for one layer""" + + # create as generic batch + super(GenericMeanScatteringBatch, self).append_one_layer(saO2, + nr_samples) + self._nr_layers -= 1 # we're not finished + + # restrict exponential scattering to mean value for soft tissue. + self.df["layer" + str(self._nr_layers), "b_mie"] = 1.286 + + self._nr_layers += 1 + + class ColonMuscleMeanScatteringBatch(ColonMuscleBatch): """three layer batch simulating colonic tissue""" def __init__(self): super(ColonMuscleMeanScatteringBatch, self).__init__() def append_one_layer(self, saO2, n, d_ranges, nr_samples): """helper function to create parameters for one layer""" # create as generic batch super(ColonMuscleMeanScatteringBatch, self).append_one_layer(saO2, n, d_ranges, nr_samples) self._nr_layers -= 1 # we're not finished # restrict exponential scattering to mean value for soft tissue. self.df["layer" + str(self._nr_layers), "b_mie"] = 1.286 self._nr_layers += 1 -# 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) + +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__() + + def append_one_layer(self, vhb, sao2, a_mie, b_mie, d, n, g, nr_samples): + """helper function to create parameters for one layer""" + + # create layer elements + self.df["layer" + str(self._nr_layers), "vhb"] = vhb + self.df["layer" + str(self._nr_layers), "sao2"] = sao2 + self.df["layer" + str(self._nr_layers), "a_mie"] = a_mie + self.df["layer" + str(self._nr_layers), "b_mie"] = b_mie + self.df["layer" + str(self._nr_layers), "d"] = d + self.df["layer" + str(self._nr_layers), "n"] = n + self.df["layer" + str(self._nr_layers), "g"] = g + self._nr_layers += 1 + + 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.1, 0.02], [0.7, 0.1], 18.9*100., 1.286, + 500 * 10 ** -6, 1.38, 0.9, + nr_samples) + self.append_one_layer(0.04, 0.7, 18.9*100., 1.286, 500 * 10 ** -6, + 1.36, 0.9, + nr_samples) + self.append_one_layer(0.04, 0.7, 18.9*100., 1.286, 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 640e08e8fb..418f141a21 100644 --- a/Modules/Biophotonics/python/iMC/mc/factories.py +++ b/Modules/Biophotonics/python/iMC/mc/factories.py @@ -1,100 +1,101 @@ ''' Created on Oct 15, 2015 @author: wirkert ''' from mc.tissuemodels import AbstractTissue, GenericTissue from mc.batches import AbstractBatch from mc.batches import GenericBatch, LessGenericBatch, GenericMeanScatteringBatch from mc.batches import ColonMuscleBatch, ColonMuscleMeanScatteringBatch +from mc.batches import VisualizationBatch 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 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 GenericMeanScatteringFactory(GenericMcFactory): def create_batch_to_simulate(self): return GenericMeanScatteringBatch() def __init__(self): ''' Constructor ''' class ColonMuscleMeanScatteringFactory(GenericMcFactory): def create_batch_to_simulate(self): return ColonMuscleMeanScatteringBatch() 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/tissuemodels.py b/Modules/Biophotonics/python/iMC/mc/tissuemodels.py index 504d310db9..d645131bcc 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_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 + \ + " - vhb: " + "%.1f" % self.uas[i].bvf * 100. + \ + "%; sao2: " + "%.1f" % self.uas[i].saO2 * 100. + \ "%; 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/linear.py b/Modules/Biophotonics/python/iMC/regression/linear.py index 0905b513cc..2833005c14 100644 --- a/Modules/Biophotonics/python/iMC/regression/linear.py +++ b/Modules/Biophotonics/python/iMC/regression/linear.py @@ -1,125 +1,80 @@ ''' Created on Oct 19, 2015 @author: wirkert ''' import numpy as np from scipy.interpolate import interp1d from mc.usuag import get_haemoglobin_extinction_coefficients class LinearSaO2Unmixing(object): ''' classdocs ''' def __init__(self, use_LCTF=False): - eHb02 = 0 - eHb = 0 - if use_LCTF: - # oxygenated haemoglobin extinction coefficients - eHbO2 = np.array([34772.8, - 27840.93333, - 23748.8 , - 21550.8 , - 21723.46667, - 28064.8 , - 39131.73333, - 45402.93333, - 42955.06667, - 40041.73333, - 42404.4 , - 36333.6 , - 22568.26667, - 6368.933333, - 1882.666667, - 1019.333333, - 664.6666667, - 473.3333333, - 376.5333333, - 327.2 , - 297.0666667],) - # deoxygenated haemoglobin extinction coefficients - eHb = [18031.73333 , - 15796.8 , - 17365.33333 , - 21106.53333 , - 26075.06667 , - 32133.2 , - 39072.66667 , - 46346.8 , - 51264 , - 50757.33333 , - 45293.33333 , - 36805.46667 , - 26673.86667 , - 17481.73333 , - 10210.13333 , - 7034 , - 5334.533333 , - 4414.706667 , - 3773.96 , - 3257.266667 , - 2809.866667] - else: - window_size = 6 - eHbO2_map, eHb_map = get_haemoglobin_extinction_coefficients() - parse_wavelengths = np.arange(450, 720, 2) * 10 ** -9 - eHbO2_parsed = eHbO2_map(parse_wavelengths) - eHb_parsed = eHb_map(parse_wavelengths) - eHbO2 = np.convolve(eHbO2_parsed, - np.ones(window_size, dtype=float) / - float(window_size), mode="same") - eHb = np.convolve(eHb_parsed, - np.ones(window_size, dtype=float) / - float(window_size), mode="same") - - eHbO2_map2 = interp1d(parse_wavelengths, eHbO2) - eHb_map2 = interp1d(parse_wavelengths, eHb) - wavelengths = np.arange(470, 680, 10) * 10 ** -9 - - eHbO2 = eHbO2_map2(wavelengths) - eHb = eHb_map2(wavelengths) - + # oxygenated haemoglobin extinction coefficients + eHbO2 = np.array([ \ + 50104 , + 33209.2, + 319.6 , + 32613.2, + 26629.2, + 20035.2, + 3200 + , 290 + ]) + # deoxygenated haemoglobin extinction coefficients + eHb = np.array([ + 37020 , + 16156.4 , + 3226.56 , + 53788 , + 14550 , + 25773.6 , + 14677.2 + , 1794.28 + ]) nr_total_wavelengths = len(eHb) # to account for scattering losses we allow a constant offset scattering = np.ones(nr_total_wavelengths) # put eHbO2, eHb and scattering term in one measurement matrix self.H = np.vstack((eHbO2, eHb, scattering)).T - self.selected_features = np.arange(0, nr_total_wavelengths, 1) + self.lsq_solution_matrix = np.dot(np.linalg.inv(np.dot(self.H.T, + self.H)), + self.H.T) + self.selected_features = np.arange(0, nr_total_wavelengths, 1) def set_selected_features(self, selected_features): """ Parameters: selected_features: index array with indices of selected features. Example np.array([0, 3, 19]) would select the 1st, 4th, and 19th feature. Will be initialized by default to all features. """ self.selected_features = selected_features def fit(self, X, y, weights=None): """only implemented to fit to the standard sklearn framework.""" pass def predict(self, X): """predict like in sklearn: Parameters: X: nrsamples x nr_features matrix of samples to predict for regression Returns: y: array of shape [nr_samples] with values for predicted oxygenation """ - # sortout not selected features - H_sel = self.H[self.selected_features, :] - X_sel = X[:, self.selected_features] # do least squares estimation - oxy_test, deoxy, s = np.linalg.lstsq(H_sel, X_sel.T)[0] + oxy_test, deoxy, s = np.dot(self.lsq_solution_matrix, X.T) # calculate oxygenation = oxygenated blood / total blood saO2 = oxy_test / (oxy_test + deoxy) return np.clip(saO2, 0., 1.) diff --git a/Modules/Biophotonics/python/iMC/regression/preprocessing.py b/Modules/Biophotonics/python/iMC/regression/preprocessing.py index 1c8c462e3e..d4f224ccd4 100644 --- a/Modules/Biophotonics/python/iMC/regression/preprocessing.py +++ b/Modules/Biophotonics/python/iMC/regression/preprocessing.py @@ -1,94 +1,93 @@ ''' Created on Oct 26, 2015 @author: wirkert ''' import numpy as np from sklearn.preprocessing import Normalizer def preprocess2(df, nr_samples=None, snr=None, movement_noise_sigma=None, magnification=None, bands_to_sortout=None): - # first set 0 reflectances to nan df["reflectances"] = df["reflectances"].replace(to_replace=0., value=np.nan) # remove nan df.dropna(inplace=True) # extract nr_samples samples from data if nr_samples is not None: df = df.sample(nr_samples) # get reflectance and oxygenation X = df.reflectances - if bands_to_sortout.size > 0: + if bands_to_sortout is not None and bands_to_sortout.size > 0: X.drop(X.columns[bands_to_sortout], axis=1, inplace=True) snr = np.delete(snr, bands_to_sortout) X = X.values 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 camera_noise = 0. if snr is not None: sigmas = X / snr noises = np.random.normal(loc=0., scale=1, size=X.shape) camera_noise = sigmas*noises movement_noise = 0. if movement_noise_sigma is not None: nr_bands = X.shape[1] nr_samples = X.shape[0] - # we assume pretty high correlation between neighboring bands + # we assume no correlation between neighboring bands CORRELATION_COEFFICIENT = 0.0 movement_variance = movement_noise_sigma ** 2 movement_variances = np.ones(nr_bands) * movement_variance movement_covariances = np.ones(nr_bands-1) * CORRELATION_COEFFICIENT * \ movement_variance movement_covariance_matrix = np.diag(movement_variances) + \ np.diag(movement_covariances, -1) + \ np.diag(movement_covariances, 1) # percentual sample errors sample_errors_p = np.random.multivariate_normal(mean=np.zeros(nr_bands), cov=movement_covariance_matrix, size=nr_samples) # errors w.r.t. the curve height. movement_noise = X * sample_errors_p X += camera_noise + movement_noise X = np.clip(X, 0.00001, 1.) # do normalizations X = normalize(X) return X, y def preprocess(batch, nr_samples=None, snr=None, movement_noise_sigma=None, magnification=None, bands_to_sortout=None): X, y = preprocess2(batch, nr_samples, snr, movement_noise_sigma, magnification, bands_to_sortout) 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_analyze_da_in_silico.py b/Modules/Biophotonics/python/iMC/script_analyze_da_in_silico.py index d2f2031f7f..346a2a40a8 100644 --- a/Modules/Biophotonics/python/iMC/script_analyze_da_in_silico.py +++ b/Modules/Biophotonics/python/iMC/script_analyze_da_in_silico.py @@ -1,346 +1,348 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ +import seaborn as sns + # import everything as done in IPCAI in silico evaluation from script_analyze_ipcai_in_silico import * # additionally we need the weights estimation functionality from regression.domain_adaptation import estimate_weights_random_forests # also the folder has changed there we store the data. sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2016_02_22_DA_in_silico" class WeightedBatch(luigi.Task): which_source = luigi.Parameter() which_target = luigi.Parameter() def requires(self): return tasks_mc.CameraBatch(self.which_source), \ tasks_mc.CameraBatch(self.which_target) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "adapted_" + self.which_source + "_with_" + self.which_target + ".txt")) def run(self): # get data df_source = pd.read_csv(self.input()[0].path, header=[0, 1]) df_target = pd.read_csv(self.input()[1].path, header=[0, 1]) # first extract X_source and X_target, preprocessed at standard noise # level X_source, y_source = preprocess(df_source, w_percent=w_standard) X_target, y_target = preprocess(df_target, w_percent=w_standard) # train a classifier to determine probability for specific class weights = estimate_weights_random_forests(X_source, X_target, X_source) # add weight to dataframe df_source["weights"] = weights # finally save the dataframe with the added weights df_source.to_csv(self.output().path, index=False) class DaNoisePlot(luigi.Task): """ Very similar to NoisePlot in IPCAI in silico evaluation but with weighted data coming in. """ which_train = luigi.Parameter() which_test = luigi.Parameter() def requires(self): return WeightedBatch(self.which_train, self.which_test), \ tasks_mc.CameraBatch(self.which_test) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "noise_da_plot_train_" + self.which_train + "_test_" + self.which_test + ".png")) def run(self): # get data df_train = pd.read_csv(self.input()[0].path, header=[0, 1]) df_test = pd.read_csv(self.input()[1].path, header=[0, 1]) df = evaluate_data(df_train, noise_levels, df_test, noise_levels) standard_plotting(df) # finally save the figure plt.savefig(self.output().path, dpi=500, bbox_inches='tight') class GeneratingDistributionPlot(luigi.Task): which_source = luigi.Parameter() which_target = luigi.Parameter() def requires(self): return WeightedBatch(self.which_source, self.which_target), \ tasks_mc.CameraBatch(self.which_target) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "generating_distribution_" + self.which_source + "_adapted_to_" + self.which_target + ".png")) def run(self): # get data df_source = pd.read_csv(self.input()[0].path, header=[0, 1]) df_target = pd.read_csv(self.input()[1].path, header=[0, 1]) # create dataframe suited for plotting # we do a weighted sampling with replacement to be able to create some # plots there the data distribution is visible. nr_samples = 100 # first data weighted by domain adaptation df_source_adapted = df_source["layer0"].copy() df_source_adapted["weights"] = df_source["weights"] df_source_adapted["data"] = "adapted" df_source_adapted = df_source_adapted.sample(n=nr_samples, replace=True, weights="weights") # now original source data df_source = df_source["layer0"].copy() df_source["weights"] = 1 # we set the weights here to 1 df_source["data"] = "source" df_source = df_source.sample(n=nr_samples, replace=True, weights="weights") # now the target data df_target = df_target["layer0"] df_target["weights"] = 1 df_target["data"] = "target" df_target = df_target.sample(n=nr_samples, replace=True, weights="weights") # now merge all three dataframes to the dataframe used for the plotting df = pd.concat([df_source, df_source_adapted, df_target]) # since we already sampled we can get rid of weights df.drop("weights", axis=1, inplace=True) # d to um df["d"] *= 10**6 # vhb and sao2 to % df["vhb"] *= 100 df["sao2"] *= 100 # do some serious plotting g = sns.pairplot(df, vars=["vhb", "sao2", "d"], hue="data", markers=["o", "s", "D"]) # tidy up plot g.add_legend() # finally save the figure plt.savefig(self.output().path, dpi=500, bbox_inches='tight') class WeightDistributionPlot(luigi.Task): which_source = luigi.Parameter() which_target = luigi.Parameter() def requires(self): return WeightedBatch(self.which_source, self.which_target) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "weight_distribution_" + self.which_source + "_adapted_to_" + self.which_target + ".png")) def run(self): # get data df_source = pd.read_csv(self.input().path, header=[0, 1]) df_source["weights"].plot.hist(bins=100) plt.axvline(x=1, ymin=0, ymax=df_source.shape[0]) # TODO: add cumsum on top # finally save the figure plt.savefig(self.output().path, dpi=500, bbox_inches='tight') class FeatureDistributionPlot(luigi.Task): which_source = luigi.Parameter() which_target = luigi.Parameter() def requires(self): return WeightedBatch(self.which_source, self.which_target), \ tasks_mc.CameraBatch(self.which_target) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "feature_distribution_" + self.which_source + "_adapted_to_" + self.which_target + ".png")) def run(self): # get data df_source = pd.read_csv(self.input()[0].path, header=[0, 1]) df_target = pd.read_csv(self.input()[1].path, header=[0, 1]) df_f_source = format_dataframe_for_distribution_plotting(df_source) df_f_target = format_dataframe_for_distribution_plotting(df_target) df_f_adapted = format_dataframe_for_distribution_plotting(df_source, weights=df_source["weights"].values.squeeze()) # build a combined df df_f_source["data"] = "source" df_f_target["data"] = "target" df_f_adapted["data"] = "adapted" df = pd.concat([df_f_source, df_f_target, df_f_adapted]) # do the plotting grid = sns.FacetGrid(df, col="w", hue="data", col_wrap=3, size=1.5) grid.map(plt.plot, "bins", "frequency") # tidy up plot grid.fig.tight_layout(w_pad=1) grid.add_legend() grid.set(xticks=(0, 1)) # finally save the figure plt.savefig(self.output().path, dpi=500) class DAvNormalPlot(luigi.Task): which_train = luigi.Parameter() which_test = luigi.Parameter() which_train_no_covariance_shift = luigi.Parameter() def requires(self): return WeightedBatch(self.which_train, self.which_test), \ tasks_mc.CameraBatch(self.which_test), \ tasks_mc.CameraBatch(self.which_train_no_covariance_shift) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "da_v_normal_train_" + self.which_train + "_test_" + self.which_test + ".png")) def run(self): # get data df_train = pd.read_csv(self.input()[0].path, header=[0, 1]) df_test = pd.read_csv(self.input()[1].path, header=[0, 1]) df_train_no_covariance_shift = pd.read_csv(self.input()[2].path, header=[0, 1]) evaluation_setups = [EvaluationStruct("Proposed", rf)] # evaluate the different methods df_adapted = evaluate_data(df_train, noise_levels, df_test, noise_levels, evaluation_setups=evaluation_setups) df_adapted["data"] = "adapted" df_no_adaptation = evaluate_data( df_train.drop("weights", axis=1), noise_levels, df_test, noise_levels, evaluation_setups=evaluation_setups) df_no_adaptation["data"] = "source" df_no_covariance_shift = evaluate_data( df_train_no_covariance_shift, noise_levels, df_test, noise_levels, evaluation_setups=evaluation_setups) df_no_covariance_shift["data"] = "target" df = pd.concat([df_adapted, df_no_adaptation, df_no_covariance_shift]) # plot it sns.boxplot(data=df, x="noise added [sigma %]", y="Errors", hue="data", hue_order=["source", "adapted", "target"], fliersize=0) # tidy up plot plt.ylim((0, 40)) plt.legend(loc='upper left') # finally save the figure plt.savefig(self.output().path, dpi=500) def format_dataframe_for_distribution_plotting(df, weights=None): if weights is None: weights = np.ones(df.shape[0]) bins = np.arange(0, 1, 0.01) # we're only interested in reflectance information df_formatted = df.loc[:, "reflectances"] # to [nm] for plotting df_formatted.rename(columns=lambda x: float(x)*10**9, inplace=True) # transform data to a histogram df_formatted = df_formatted.apply(lambda x: pd.Series(np.histogram(x, bins=bins, weights=weights, normed=True)[0]), axis=0) # convert to long form using bins as identifier df_formatted["bins"] = bins[1:] df_formatted = pd.melt(df_formatted, id_vars=["bins"], var_name="w", value_name="frequency") return df_formatted if __name__ == '__main__': logging.basicConfig(filename='log/da_in_silico_plots_' + 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() source_domain = "ipcai_revision_generic" target_domain = "ipcai_revision_colon_test" sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) # check how the graph looks with same domains for training and testing w.add(DaNoisePlot(which_train="ipcai_revision_colon_train", which_test="ipcai_revision_colon_test")) # check how the graph looks without domain adaptation w.add(NoisePlot(which_train=source_domain, which_test=target_domain)) # Set a different testing domain to evaluate domain sensitivity w.add(DaNoisePlot(which_train=source_domain, which_test=target_domain)) w.add(WeightDistributionPlot(which_source=source_domain, which_target=target_domain)) w.add(FeatureDistributionPlot(which_source=source_domain, which_target=target_domain)) # also plot distributions for equal domains to check for errors in data w.add(FeatureDistributionPlot(which_source="ipcai_revision_colon_mean_scattering_train", which_target="ipcai_revision_colon_mean_scattering_test")) # plot how the generating model data (e.g. sao2 and vhb) is distributed w.add(GeneratingDistributionPlot(which_source=source_domain, which_target=target_domain)) w.add(DAvNormalPlot(which_train=source_domain, which_test=target_domain, which_train_no_covariance_shift="ipcai_revision_colon_train")) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_silico.py b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_silico.py index 41766af8fb..230b3262b5 100644 --- a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_silico.py +++ b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_silico.py @@ -1,334 +1,373 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import logging import datetime from collections import namedtuple +import matplotlib import numpy as np import pandas as pd from pandas import DataFrame import luigi import matplotlib.pyplot as plt from sklearn.ensemble.forest import RandomForestRegressor -import seaborn as sns import tasks_mc import scriptpaths as sp from regression.preprocessing import preprocess, preprocess2 from regression.linear import LinearSaO2Unmixing sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2016_02_22_IPCAI_revision_in_silico" sp.FINALS_FOLDER = "finals" sp.RECORDED_WAVELENGTHS = np.arange(470, 680, 10) * 10 ** -9 sp.MC_DATA_FOLDER = "mc_data" -w_standard = 0.1 # for this evaluation we add 10% noise -sns.set_style("whitegrid") -# sns.set_context("poster") +w_standard = 10. # for this evaluation we add 10% noise + +font = {'family' : 'normal', + 'size' : 20} + +matplotlib.rc('font', **font) + # setup standard random forest rf = RandomForestRegressor(10, min_samples_leaf=10, max_depth=9, n_jobs=-1) EvaluationStruct = namedtuple("EvaluationStruct", "name regressor") # standard evaluation setup standard_evaluation_setups = [EvaluationStruct("Linear Beer-Lambert", LinearSaO2Unmixing()) , EvaluationStruct("Proposed", rf)] # color palette -my_colors = [sns.xkcd_rgb["pale red"], sns.xkcd_rgb["medium green"]] +my_colors = ["red", "green"] # standard noise levels -noise_levels = np.arange(0.00, 0.18, 0.02) +noise_levels = np.array([1,2,3,4,5,6,7,8,9,10, + 15,20,30,40,50,100,150,200]).astype("float") class TrainingSamplePlot(luigi.Task): which_train = luigi.Parameter() which_test = luigi.Parameter() def requires(self): return tasks_mc.CameraBatch(self.which_train), \ tasks_mc.CameraBatch(self.which_test) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "sample_plot_train_" + self.which_train + "_test_" + self.which_test + - ".png")) + ".pdf")) def run(self): # get data df_train = pd.read_csv(self.input()[0].path, header=[0, 1]) df_test = pd.read_csv(self.input()[1].path, header=[0, 1]) # for this plot we write a custom evaluation function as it is built # a little different # create a new dataframe which will hold all the generated errors df = pd.DataFrame() - nr_training_samples = np.arange(10, 10010, 50).astype(int) + nr_training_samples = np.arange(10, 15010, 50).astype(int) # not very pythonic, don't care for n in nr_training_samples: - X_test, y_test = preprocess(df_test, w_percent=w_standard) + X_test, y_test = preprocess(df_test, snr=w_standard) # only take n samples for training X_train, y_train = preprocess(df_train, nr_samples=n, - w_percent=w_standard) + snr=w_standard) regressor = rf regressor.fit(X_train, y_train) y_pred = regressor.predict(X_test) # save results to a dataframe errors = np.abs(y_pred - y_test) errors = errors.reshape(len(errors), 1) current_df = DataFrame(errors * 100, columns=["Errors"]) current_df["Method"] = "Proposed" - current_df["Number Samples"] = n + current_df["Number Samples"] = n / 10**3. df = pd.concat([df, current_df], ignore_index=True) logging.info( "Finished training classifier with {0} samples".format( str(n))) df = df.groupby("Number Samples").describe() # get the error description in the rows: df = df.unstack(-1) # get rid of multiindex by dropping "Error" level df.columns = df.columns.droplevel(0) plt.figure() - plt.plot(df.index, df["50%"], color=sns.xkcd_rgb["medium green"]) + plt.plot(df.index, df["50%"], color="green") # tidy up the plot - plt.xlabel("number of training samples") + plt.xlabel("number of training samples / 1000") plt.ylabel("absolute error [%]") plt.ylim((0, 20)) - plt.xlim((0, 15000)) + plt.xlim((0, 15)) + plt.grid() # finally save the figure - plt.savefig(self.output().path, dpi=500, + plt.savefig(self.output().path, mode="pdf", dpi=500, bbox_inches='tight') class VhbPlot(luigi.Task): which_train = luigi.Parameter() which_test = luigi.Parameter() def requires(self): return tasks_mc.CameraBatch(self.which_train), \ tasks_mc.CameraBatch(self.which_test) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "vhb_noise_plot_train_" + self.which_train + "_test_" + self.which_test + - ".png")) + ".pdf")) @staticmethod - def preprocess_vhb(batch, nr_samples=None, w_percent=None, + def preprocess_vhb(batch, nr_samples=None, snr=None, magnification=None, bands_to_sortout=None): """ For evaluating vhb we extract labels for vhb instead of sao2""" - X, y = preprocess2(batch, nr_samples, w_percent, + X, y = preprocess2(batch, nr_samples, snr, magnification, bands_to_sortout) - return X, y["vhb"] + return X, y["vhb"].values def run(self): # get data df_train = pd.read_csv(self.input()[0].path, header=[0, 1]) df_test = pd.read_csv(self.input()[1].path, header=[0, 1]) # for vhb we only evaluate the proposed method since the linear # beer-lambert is not applicable evaluation_setups = [EvaluationStruct("Proposed", rf)] df = evaluate_data(df_train, noise_levels, df_test, noise_levels, evaluation_setups=evaluation_setups, - preprocessing=self.preprocess_vhb, ) - standard_plotting(df, color_palette=[sns.xkcd_rgb["medium green"]], + preprocessing=self.preprocess_vhb) + standard_plotting(df, color_palette=["green"], xytext_position=(2, 3)) plt.ylim((0, 4)) # finally save the figure plt.savefig(self.output().path, dpi=500, bbox_inches='tight') class NoisePlot(luigi.Task): which_train = luigi.Parameter() which_test = luigi.Parameter() def requires(self): return tasks_mc.CameraBatch(self.which_train), \ tasks_mc.CameraBatch(self.which_test) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "noise_plot_train_" + self.which_train + "_test_" + self.which_test + - ".png")) + ".pdf")) def run(self): # get data df_train = pd.read_csv(self.input()[0].path, header=[0, 1]) df_test = pd.read_csv(self.input()[1].path, header=[0, 1]) df = evaluate_data(df_train, noise_levels, df_test, noise_levels) standard_plotting(df) # finally save the figure - plt.savefig(self.output().path, dpi=500, + plt.savefig(self.output().path, mode="pdf", dpi=500, bbox_inches='tight') class WrongNoisePlot(luigi.Task): which_train = luigi.Parameter() which_test = luigi.Parameter() + train_snr = luigi.FloatParameter() def requires(self): return tasks_mc.CameraBatch(self.which_train), \ tasks_mc.CameraBatch(self.which_test) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, - "wrong_noise_plot_train_" + + str(self.train_snr) + + "_wrong_noise_plot_train_" + self.which_train + "_test_" + self.which_test + - ".png")) + ".pdf")) def run(self): # get data df_train = pd.read_csv(self.input()[0].path, header=[0, 1]) df_test = pd.read_csv(self.input()[1].path, header=[0, 1]) # do same as in NoisePlot but with standard noise input - df = evaluate_data(df_train, np.ones_like(noise_levels) * w_standard, + df = evaluate_data(df_train, + np.ones_like(noise_levels) * self.train_snr, df_test, noise_levels) standard_plotting(df) # finally save the figure - plt.savefig(self.output().path, dpi=500, + plt.savefig(self.output().path, mode="pdf", dpi=500, bbox_inches='tight') +class VisualizationPlot(luigi.Task): + def requires(self): + return tasks_mc.CameraBatch("ipcai_visualization_batch_0") + + def output(self): + return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, + sp.RESULTS_FOLDER, + sp.FINALS_FOLDER, + "visualization1.pdf")) + + def run(self): + # get data + df = pd.read_csv(self.input().path, header=[0, 1]) + + wavelengths = np.array(df["reflectances"].columns.astype(float)*10**9) + plt.figure() + plt.grid() + plt.plot(wavelengths, df["reflectances"].loc[1, :], "*") + plt.xlabel("\lambda [nm]") + plt.ylabel("r") + plt.savefig(self.output().path+"3.pdf", mode="pdf", dpi=500, + bbox_inches='tight') + + + def evaluate_data(df_train, w_train, df_test, w_test, evaluation_setups=None, preprocessing=None): """ Our standard method to evaluate the data. It will fill a DataFrame df which saves the errors for each evaluated setup""" if evaluation_setups is None: evaluation_setups = standard_evaluation_setups if preprocessing is None: preprocessing = preprocess if ("weights" in df_train) and df_train["weights"].size > 0: weights = df_train["weights"].as_matrix().squeeze() else: weights = np.ones(df_train.shape[0]) # create a new dataframe which will hold all the generated errors df = pd.DataFrame() for one_w_train, one_w_test in zip(w_train, w_test): # setup testing function - X_test, y_test = preprocessing(df_test, w_percent=one_w_test) + X_test, y_test = preprocessing(df_test, snr=one_w_test) # extract noisy data - X_train, y_train = preprocessing(df_train, w_percent=one_w_train) + X_train, y_train = preprocessing(df_train, snr=one_w_train) for e in evaluation_setups: regressor = e.regressor regressor.fit(X_train, y_train, weights) y_pred = regressor.predict(X_test) # save results to a dataframe errors = np.abs(y_pred - y_test) errors = errors.reshape(len(errors), 1) current_df = DataFrame(errors * 100, columns=["Errors"]) current_df["Method"] = e.name - current_df["noise added [sigma %]"] = int(one_w_test * 100) + current_df["SNR"] = int(one_w_test) df = pd.concat([df, current_df], ignore_index=True) return df def standard_plotting(df, color_palette=None, xytext_position=None): if color_palette is None: color_palette = my_colors if xytext_position is None: xytext_position = (2, 15) plt.figure() # group it by method and noise level and get description on the errors - df_statistics = df.groupby(['Method', 'noise added [sigma %]']).describe() + df_statistics = df.groupby(['Method', 'SNR']).describe() # get the error description in the rows: df_statistics = df_statistics.unstack(-1) # get rid of multiindex by dropping "Error" level df_statistics.columns = df_statistics.columns.droplevel(0) # iterate over methods to plot linegraphs with error tube # probably this can be done nicer, but no idea how exactly for color, method in zip( color_palette, df_statistics.index.get_level_values("Method").unique()): df_method = df_statistics.loc[method] plt.plot(df_method.index, df_method["50%"], color=color, label=method) plt.fill_between(df_method.index, df_method["25%"], df_method["75%"], facecolor=color, edgecolor=color, alpha=0.5) - plt.xlabel("noise added [sigma %]") - plt.ylabel("absolute error [%]") # add annotation: - median_proposed = df_statistics.T["Proposed", 0].loc["50%"] - plt.gca().annotate('error: ' + - "{0:.1f}".format(median_proposed) + "%", - xy=(0, median_proposed), xytext=xytext_position, - va="center", ha="center", - bbox=dict(boxstyle="round4", - fc=sns.xkcd_rgb["medium green"], - alpha=0.5), - arrowprops=dict(arrowstyle="-|>", - connectionstyle="arc3,rad=-0.2", - fc="w")) + # median_proposed = df_statistics.T["Proposed", 0].loc["50%"] + # plt.gca().annotate('error: ' + + # "{0:.1f}".format(median_proposed) + "%", + # xy=(0, median_proposed), xytext=xytext_position, + # va="center", ha="center", + # bbox=dict(boxstyle="round4", + # fc=sns.xkcd_rgb["medium green"], + # alpha=0.5), + # arrowprops=dict(arrowstyle="-|>", + # connectionstyle="arc3,rad=-0.2", + # fc="w")) # tidy up the plot - plt.ylim((0, 20)) - plt.xlim((0, 15)) + plt.ylim((0, 40)) + plt.gca().set_xticks(np.arange(0, 200, 10), minor=True) + plt.xlabel("SNR") + plt.ylabel("absolute error [%]") + plt.grid() plt.legend() if __name__ == '__main__': logging.basicConfig(filename='log/in_silico_plots_' + 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() train = "ipcai_revision_colon_mean_scattering_train" test = "ipcai_revision_colon_mean_scattering_test" sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) w.add(TrainingSamplePlot(which_train=train, which_test=test)) w.add(NoisePlot(which_train=train, which_test=test)) - w.add(WrongNoisePlot(which_train=train, which_test=test)) + w.add(WrongNoisePlot(which_train=train, which_test=test, train_snr=10.)) + w.add(WrongNoisePlot(which_train=train, which_test=test, train_snr=50.)) + w.add(WrongNoisePlot(which_train=train, which_test=test, train_snr=200.)) # Set a different testing domain to evaluate domain sensitivity - w.add(NoisePlot(which_train=train, which_test="ipcai_revision_generic")) + w.add(NoisePlot(which_train=train, + which_test="ipcai_revision_generic_mean_scattering_test")) w.add(VhbPlot(which_train=train, which_test=test)) + w.add(VisualizationPlot()) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo.py b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo.py index d0647b7fb0..5f9f210c71 100644 --- a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo.py +++ b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo.py @@ -1,424 +1,424 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import Image import ImageEnhance import pickle import copy 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_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 from msi.msi import Msi 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, normalize sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2015_11_12_IPCAI_Neil_InVivo" sp.DATA_FOLDER = "colon_images/not_registered" sp.MC_DATA_FOLDER = "mc_data2" sp.FINALS_FOLDER = "OxySecondCamera" sp.RECORDED_WAVELENGTHS = np.arange(400, 730, 10) * 10 ** -9 def resort_image_wavelengths(collapsed_image): return collapsed_image def resort_wavelengths(msi): """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 sp.resort_wavelengths = resort_wavelengths sp.bands_to_sortout = np.array([0, 1, 2, 3, 4, 5, 6 , 7, 8, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32]) class IPCAICreateOxyImageTask(luigi.Task): folder_name = luigi.Parameter() df_prefix = luigi.Parameter() def requires(self): return IPCAIEstimateTissueParametersTask(folder_name=self.folder_name, batch_prefix= self.df_prefix), \ MultiSpectralImageFromPNGFiles(folder_name=self.folder_name) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, self.folder_name + "_" + self.df_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) # segmentation window: seg_windows = np.array([ [135, 332, 227, 154 ], [434, 280, 161, 161 ], [428, 381, 99 , 99 ], [701, 368, 161, 148 ], [756, 363, 151, 139 ], [630, 340, 150, 129 ], [468, 505, 97 , 108 ], [679, 499, 114, 114 ], [493, 513, 89 , 89 ], [336, 72 , 177 , 152 ], [313, 122, 132, 139 ], [697, 311, 113, 103 ], [881, 400, 128, 112 ], [634, 452, 124, 103 ], [632, 452, 120, 108 ], [675, 478, 122, 126 ], [774, 577, 128, 126 ], [490, 501, 130, 99 ], [495, 597, 112, 91 ], [499, 482, 103, 95 ], [681, 408, 122, 85 ], [747, 439, 120, 99 ], [634, 509, 99 , 103 ], [571, 499, 108, 85 ], [760, 495, 97 , 85 ], [636, 415, 93 , 77 ], [872, 445, 112, 103 ], [659, 462, 118, 103 ], [610, 452, 89 , 73 ], [840, 616, 89 , 75 ], [620, 419, 97 , 81 ], [597, 404, 79 , 83 ], [657, 400, 97 , 81 ], [501, 341, 95 , 77 ], [612, 398, 89 , 83 ], [655, 437, 79 , 68 ], [776, 394, 89 , 83 ], [597, 462, 93 , 71 ], [681, 505, 81 , 81 ], [566, 427, 83 , 71 ], [727, 427, 81 , 69 ], [772, 486, 103, 89 ], [515, 511, 83 , 73 ], [673, 412, 108, 91 ], [895, 505, 101, 93 ], [474, 365, 120, 101 ]]) # plot f, axarr = plt.subplots(1, 2) # create artificial rgb rgb_image = msi.get_image()[:, :, [-1, 6, 0]] 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(2.) plot_image = np.array(im) top_left_axis = axarr[0] top_left_axis.imshow(plot_image, interpolation='nearest') top_left_axis.set_title("synthetic rgb", fontsize=5) top_left_axis.xaxis.set_visible(False) top_left_axis.yaxis.set_visible(False) plt.set_cmap("jet") # postprocess parametric maps image[np.isnan(image)] = 0. image[np.isinf(image)] = 0. image[image > 1.] = 1. image[image < 0.] = 0. image[0, 0] = 0. image[0, 1] = 1. try: folder_nr = int(self.folder_name[-2:]) except ValueError: folder_nr = int(self.folder_name[-1]) seg_window = seg_windows[folder_nr - 1, :] evaluation_image = copy.deepcopy(image) # set out-segmented pixels to -1 evaluation_image[(msi.get_image() >= np.max(msi.get_image())).any(axis=2)] = -1 evaluation_image[(msi.get_image() <= 0).any(axis=2)] = -1 # only take elements in window evaluation_image = evaluation_image[ seg_window[1]:seg_window[1] + seg_window[3], seg_window[0]:seg_window[0] + seg_window[2]] # throw away elements not in window result = np.mean(evaluation_image[evaluation_image >= 0]) logging.info("mean of image is: " + str(result)) fd = open(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, 'results_summary.csv'), 'a') fd.write(str(folder_nr) + ", " + str(result) + "\n") fd.close() # write the windowed evaluation image sitk_img = sitk.GetImageFromArray(evaluation_image, isVector=True) sitk.WriteImage(sitk_img, self.output().path + "window.nrrd") # plot parametric maps at.plot_axis(axarr[1], image[:, :], "oxygenation [%]") plt.savefig(self.output().path, dpi=1000, bbox_inches='tight') class IPCAIEstimateTissueParametersTask(luigi.Task): folder_name = luigi.Parameter() df_prefix = luigi.Parameter() def requires(self): return IPCAIPreprocessMSI(folder_name=self.folder_name), \ IPCAITrainRegressor(batch_prefix= self.df_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, "processed", sp.FINALS_FOLDER, self.folder_name + "_" + self.df_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 sitk_image = estimate_image(msi, e) # save data outFile = self.output().open('w') outFile.close() sitk.WriteImage(sitk_image, self.output().path) class IPCAITrainRegressor(luigi.Task): df_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "reg_" + self.df_prefix)) def requires(self): - return tasks_mc.CameraBatch(self.df_prefix) + return tasks_mc.SpectroCamBatch(self.df_prefix) def run(self): # extract data from the batch f = open(self.input().path, "r") batch_train = pickle.load(f) X, y = preprocess(batch_train, w_percent=0.1, bands_to_sortout=sp.bands_to_sortout) # train regressor reg = RandomForestRegressor(max_depth=10, n_estimators=10, min_samples_leaf=1, 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): folder_name = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.folder_name + "_preprocessed.nrrd")) def requires(self): return IPCAICorrectImagingSetupTask(folder_name=self.folder_name) 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""" folder_name = luigi.Parameter() def requires(self): return MultiSpectralImageFromPNGFiles(folder_name=self.folder_name), \ FlatfieldFromMeasurement() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.folder_name + "_image_setup_corrected.nrrd")) def run(self): """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) # msi.set_image(gaussian_filter(msi.get_image().astype(np.float), (1, 1, 0))) # correct image by flatfield msimani.flatfield_correction(msi, flatfield) sp.resort_wavelengths(msi) pre.touch_and_save_msi(msi, self.output()) class FlatfieldFromMeasurement(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "flatfield" + ".nrrd")) def run(self): msi = Msi() msi.set_image(np.array([ 0.90733563 , 3.6583764 , 27.30637258, 52.35073931, 103.0408631, 176.812168 , 392.361853 , 620.0531964, 712.2926615, 1052.643461, 1202.266783, 1483.629833, 1730.168541, 2041.730527, 2296.256574, 2399.004343, 2536.580738, 2736.328832, 2760.046341, 2927.307208, 3539.747558, 3420.653065, 3229.74619 , 3435.280705, 3108.530329, 2817.607701, 2612.795573, 2053.155475, 1343.305582, 721.058596 , 340.6692005, 189.2180019, 77.61268264 ])) writer = NrrdWriter(msi) writer.write(self.output().path) class MultiSpectralImageFromPNGFiles(luigi.Task): folder_name = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.folder_name + ".nrrd")) def run(self): reader = PngReader() msi = reader.read(os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, self.folder_name, "image sample")) writer = NrrdWriter(msi) writer.write(self.output().path) if __name__ == '__main__': # configure logging logging.basicConfig(filename='invivo_neil_colon' + str(datetime.datetime.now()) + '.log', level=logging.INFO) ch = logging.StreamHandler() ch.setLevel(logging.INFO) logger = logging.getLogger() logger.addHandler(ch) # configure luigi 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) all_folders = os.listdir(image_file_folder) # each complete stack has elements F0 - F7 only_colon_folders = [ f for f in all_folders if f.startswith("SaO2") ] for f in only_colon_folders: main_task = IPCAICreateOxyImageTask( folder_name=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_liver.py b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_liver.py index a86459b78e..0f555ad8e0 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,478 +1,479 @@ # -*- 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 pandas as pd import luigi import matplotlib.pyplot as plt import SimpleITK as sitk from sklearn.ensemble.forest import RandomForestRegressor import seaborn as sns import tasks_analyze as at import tasks_mc import scriptpaths as sp from msi.msi import Msi 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.preprocessing import preprocess, preprocess2 from msi.io.tiffringreader import TiffRingReader sp.ROOT_FOLDER = "/media/wirkert/data/Data/2015_11_12_IPCAI_in_vivo" sp.DATA_FOLDER = "all_liver_images" -sp.FINALS_FOLDER = "liver_all" +sp.FINALS_FOLDER = "liver_all_SNR10" sp.FLAT_FOLDER = "/media/wirkert/data/Data/" + \ "2016_01_19_Flatfield_and_Dark_Laparoscope/Flatfields" sp.DARK_FOLDER = "/media/wirkert/data/Data/" + \ "2016_01_19_Flatfield_and_Dark_Laparoscope/Dark" AVERAGE_FOLDER = "average_intensity" AVERAGE_FINALS_FOLDER = "finals" sns.set_style("whitegrid") font = {'family' : 'normal', 'size' : 30} # sp.RECORDED_WAVELENGTHS = np.arange(470, 680, 10) * 10 ** -9 new_order = [0, 1, 2, 3, 4, 5, 6, 7] def get_image_files_from_folder(folder): # small helper function to get all the image files in a folder image_files = [f for f in os.listdir(folder) if os.path.isfile(os.path.join(folder, f))] image_files.sort() image_files = [f for f in image_files if f.endswith(".tiff")] return image_files def resort_image_wavelengths(collapsed_image): return collapsed_image[:, new_order] def resort_wavelengths(msi): """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 sp.resort_wavelengths = resort_wavelengths sp.bands_to_sortout = np.array([]) class ResultsFile(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "results.csv")) class OxyAndVhbOverTimeTask(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "liver_oxy_over_time.pdf")) def requires(self): return ResultsFile() def run(self): df = pd.read_csv(self.input().path, index_col=0) # determine times from start: image_name_strings = df["image name"].values time_strings = map(lambda s: s[ s.find("2014-08-03_")+11:s.find("2014-08-03_")+19], image_name_strings) time_in_s = map(lambda s: int(s[0:2]) * 3600 + int(s[3:5]) * 60 + int(s[6:]), time_strings) df["time since drug delivery [s]"] = np.array(time_in_s) - time_in_s[0] # print oxy over time as scatterplot. ax = df.plot.scatter(x="time since drug delivery [s]", y="oxygenation mean [%]", - s=100, + s=100, alpha=0.5, fontsize=30) ax.set_xlim((-1, 70)) plt.axvline(x=0, ymin=0, ymax=1, linewidth=2) plt.axvline(x=56, ymin=0, ymax=1, linewidth=2) ax.annotate('drug delivery', xy=(0, ax.get_ylim()[1]), xycoords='data', xytext=(0, 0), fontsize=30, textcoords='offset points') ax.annotate('porcine death', xy=(56, ax.get_ylim()[1]), xycoords='data', xytext=(-100, 0), fontsize=30, textcoords='offset points') ax.yaxis.label.set_size(30) ax.xaxis.label.set_size(30) df.to_csv(self.input().path) # create and save vhb plot plt.savefig(self.output().path, dpi=250, bbox_inches='tight', mode="pdf") # print vhb over time as scatterplot. ax = df.plot.scatter(x="time since drug delivery [s]", y="blood volume fraction mean [%]", - s=100, + s=100, alpha=0.5, fontsize=30) ax.set_xlim((-1, 70)) plt.axvline(x=0, ymin=0, ymax=1, linewidth=2) plt.axvline(x=56, ymin=0, ymax=1, linewidth=2) ax.annotate('drug delivery', xy=(0, ax.get_ylim()[1]), xycoords='data', xytext=(0, 0), fontsize=30, textcoords='offset points') ax.annotate('porcine death', xy=(56, ax.get_ylim()[1]), xycoords='data', xytext=(-100, 0), fontsize=30, textcoords='offset points') ax.yaxis.label.set_size(30) ax.xaxis.label.set_size(30) plt.savefig(self.output().path + "_vhb_mean.pdf", dpi=250, bbox_inches='tight', mode="pdf") class IPCAICreateOxyImageTask(luigi.Task): image_name = luigi.Parameter() df_prefix = luigi.Parameter() def requires(self): return IPCAITrainRegressor(df_prefix=self.df_prefix), \ Flatfield(flatfield_folder=sp.FLAT_FOLDER), \ - SingleMultispectralImage(image=self.image_name) + SingleMultispectralImage(image=self.image_name), \ + Dark(dark_folder=sp.DARK_FOLDER) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, self.image_name + "_" + self.df_prefix + "_oxy_summary" + ".png")) def run(self): nrrd_reader = NrrdReader() tiff_ring_reader = TiffRingReader() # read the flatfield flat = nrrd_reader.read(self.input()[1].path) + dark = nrrd_reader.read(self.input()[3].path) # read the msi nr_filters = len(sp.RECORDED_WAVELENGTHS) msi, segmentation = tiff_ring_reader.read(self.input()[2].path, nr_filters) # only take into account not saturated pixels. segmentation = np.max(msi.get_image(), axis=-1) < 2000. # read the regressor e_file = open(self.input()[0].path, 'r') e = pickle.load(e_file) # correct image setup position_filter_nr_in_string = self.image_name.find(" 2014") - 1 filter_nr = int(self.image_name[ position_filter_nr_in_string:position_filter_nr_in_string+1]) 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) + msimani.image_correction(msi, flat, dark) # 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, time = 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 # first oxygenation plt.figure() oxy_image = image[:, :, 0] segmentation[0, 0] = 1 segmentation[0, 1] = 1 oxy_image = np.ma.masked_array(oxy_image, ~segmentation) 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.0 oxy_image[0, 1] = 1. at.plot_axis(plt.gca(), oxy_image, "oxygenation [%]. Mean " + "{0:.1f}".format((oxy_mean * 100.)) + "%") plt.savefig(self.output().path, dpi=250, bbox_inches='tight') # second blood volume fraction plt.figure() vhb_image = image[:, :, 1] vhb_image = np.ma.masked_array(vhb_image, ~segmentation) vhb_image[np.isnan(vhb_image)] = 0. vhb_image[np.isinf(vhb_image)] = 0. vhb_image[0, 0] = 0.0 - vhb_image[0, 1] = 0.05 - vhb_image = np.clip(vhb_image, 0.0, 0.05) + vhb_image[0, 1] = 0.1 + vhb_image = np.clip(vhb_image, 0.0, 0.1) vhb_mean = np.mean(vhb_image) at.plot_axis(plt.gca(), vhb_image, "vhb [%]. Mean " + "{0:.1f}".format((vhb_mean * 100.)) + "%") plt.savefig(self.output().path + "vhb.png", dpi=250, bbox_inches='tight') # store results summary in dataframe df_image_results = pd.DataFrame(data=np.expand_dims([self.image_name, oxy_mean * 100., vhb_mean * 100., time], 0), columns=["image name", "oxygenation mean [%]", "blood volume fraction mean [%]", "time to estimate"]) results_file = os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "results.csv") if os.path.isfile(results_file): df_results = pd.read_csv(results_file, index_col=0) df_results = pd.concat((df_results, df_image_results)).reset_index( drop=True ) else: df_results = df_image_results df_results.to_csv(results_file) # 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): df_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "reg_" + self.df_prefix)) def requires(self): return tasks_mc.CameraBatch(self.df_prefix) def run(self): # extract data from the batch df_train = pd.read_csv(self.input().path, header=[0, 1]) df_snrs = pd.read_csv("/media/wirkert/data/Data/" + "2016_01_08_Calibration_Measurements/processed/" + "Carets_Laparoscope/finals/" + "snr_for_bands.txt", index_col=0) - X, y = preprocess2(df_train, snr=df_snrs["SNR"].values, - movement_noise_sigma=0.1, - bands_to_sortout=sp.bands_to_sortout) + X, y = preprocess2(df_train, snr=10., + 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) # reg = LinearSaO2Unmixing() # save regressor regressor_file = self.output().open('w') pickle.dump(reg, regressor_file) regressor_file.close() class SingleMultispectralImage(luigi.Task): image = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, \ self.image)) class Flatfield(luigi.Task): flatfield_folder = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "flatfield" + ".nrrd")) def run(self): tiff_ring_reader = TiffRingReader() nr_filters = len(sp.RECORDED_WAVELENGTHS) # analyze all the first image files image_files = get_image_files_from_folder(self.flatfield_folder) image_files = filter(lambda image_name: "F0" in image_name, image_files) # helper function to take maximum of two images def maximum_of_two_images(image_1, image_name_2): image_2 = tiff_ring_reader.read(os.path.join(self.flatfield_folder, image_name_2), nr_filters)[0].get_image() return np.maximum(image_1, image_2) # now reduce to maximum of all the single images flat_maximum = reduce(lambda x, y: maximum_of_two_images(x, y), image_files, 0) msi = Msi(image=flat_maximum) msi.set_wavelengths(sp.RECORDED_WAVELENGTHS) # write flatfield as nrrd writer = NrrdWriter(msi) writer.write(self.output().path) class Dark(luigi.Task): dark_folder = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "dark" + ".nrrd")) def run(self): tiff_ring_reader = TiffRingReader() nr_filters = len(sp.RECORDED_WAVELENGTHS) # analyze all the first image files image_files = get_image_files_from_folder(self.dark_folder) image_files = filter(lambda image_name: "F0" in image_name, image_files) # returns the mean dark image vector of all inputted dark image # overly complicated TODO SW: make this simple code readable. dark_means = map(lambda image_name: msimani.calculate_mean_spectrum( tiff_ring_reader.read(os.path.join(self.dark_folder, image_name), nr_filters)[0]), image_files) dark_means_sum = reduce(lambda x, y: x+y.get_image(), dark_means, 0) final_dark_mean = dark_means_sum / len(dark_means) msi = Msi(image=final_dark_mean) msi.set_wavelengths(sp.RECORDED_WAVELENGTHS) # write flatfield as nrrd writer = NrrdWriter(msi) writer.write(self.output().path) if __name__ == '__main__': # root folder there the data lies logging.basicConfig(filename='liver' + 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) image_file_folder = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER) onlyfiles = get_image_files_from_folder(image_file_folder) first_invivo_image_files = filter(lambda image_name: "0 2014" in image_name, onlyfiles) for f in first_invivo_image_files: main_task = IPCAICreateOxyImageTask(image_name=f, df_prefix="ipcai_revision_colon_mean_scattering_train") w.add(main_task) w.run() oxygenation_over_time_task = OxyAndVhbOverTimeTask() w.add(oxygenation_over_time_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 ca224804dd..ae580696cb 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,450 +1,458 @@ # -*- 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 pandas as pd import luigi import matplotlib.pyplot as plt import SimpleITK as sitk from sklearn.ensemble.forest import RandomForestRegressor -import seaborn as sns import matplotlib import tasks_analyze as at import tasks_mc import scriptpaths as sp from msi.msi import Msi 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 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_jet_new_scattering" +sp.FINALS_FOLDER = "small_bowel_images_jet_new_scattering_SN10_more_brighness" sp.FLAT_FOLDER = "/media/wirkert/data/Data/" + \ "2016_01_19_Flatfield_and_Dark_Laparoscope/Flatfields" sp.DARK_FOLDER = "/media/wirkert/data/Data/" + \ "2016_01_19_Flatfield_and_Dark_Laparoscope/Dark" AVERAGE_FOLDER = "average_intensity" AVERAGE_FINALS_FOLDER = "finals" -sns.set_style("whitegrid") -sns.set_context("talk") font = {'family' : 'normal', - 'size' : 18} + 'size' : 25} matplotlib.rc('font', **font) # sp.RECORDED_WAVELENGTHS = np.arange(470, 680, 10) * 10 ** -9 new_order = [0, 1, 2, 3, 4, 5, 6, 7] def get_image_files_from_folder(folder): # small helper function to get all the image files in a folder image_files = [f for f in os.listdir(folder) if os.path.isfile(os.path.join(folder, f))] image_files.sort() image_files = [f for f in image_files if f.endswith(".tiff")] return image_files def resort_image_wavelengths(collapsed_image): return collapsed_image[:, new_order] def resort_wavelengths(msi): """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 sp.resort_wavelengths = resort_wavelengths sp.bands_to_sortout = np.array([]) class ResultsFile(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "results.csv")) class OxyOverTimeTask(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "colon_oxy_over_time.pdf")) def requires(self): return ResultsFile() def run(self): df = pd.read_csv(self.input().path, index_col=0) # determine times from start: image_name_strings = df["image name"].values - time_strings = map(lambda s: s[27:35], image_name_strings) + time_strings = map(lambda s: s[ + s.find("2015-10-08_")+11:s.find("2015-10-08_")+19], + image_name_strings) time_in_s = map(lambda s: int(s[0:2]) * 3600 + int(s[3:5]) * 60 + int(s[6:]), time_strings) df["time since first frame [s]"] = np.array(time_in_s) - time_in_s[0] - + df["time since applying first clip [s]"] = df["time since first frame [s]"] - 4 # print oxy over time as scatterplot. plt.figure() - ax = df.plot.scatter(x="time since first frame [s]", - y="oxygenation mean [%]", fontsize=16, - s=50) - ax.set_xlim((0, 600)) - - plt.axvline(x=44, ymin=0, ymax=1, linewidth=2) - plt.axvline(x=104, ymin=0, ymax=1, linewidth=2) - plt.axvline(x=204, ymin=0, ymax=1, linewidth=2) - ax.annotate('first clip', xy=(44, ax.get_ylim()[1]), + ax = df.plot.scatter(x="time since applying first clip [s]", + y="oxygenation mean [%]", fontsize=30, + s=50, alpha=0.5) + ax.set_xlim((-5, 600)) + + plt.axvline(x=0, ymin=0, ymax=1, linewidth=2) + plt.axvline(x=39, ymin=0, ymax=1, linewidth=2) + plt.axvline(x=100, ymin=0, ymax=1, linewidth=2) + plt.axvline(x=124, ymin=0, ymax=1, linewidth=2) + ax.annotate('1', xy=(0, ax.get_ylim()[1]), + fontsize=18, + color="blue", + xycoords='data', xytext=(-5, 0), + textcoords='offset points') + ax.annotate('2', xy=(39, ax.get_ylim()[1]), fontsize=18, - xycoords='data', xytext=(-40, 0), + color="blue", + xycoords='data', xytext=(-5, 0), textcoords='offset points') - ax.annotate('second clip', xy=(104, ax.get_ylim()[1]), + ax.annotate('3', xy=(100, ax.get_ylim()[1]), fontsize=18, - xycoords='data', xytext=(-15, 0), + color="blue", + xycoords='data', xytext=(-5, 0), textcoords='offset points') - ax.annotate('third clip', xy=(204, ax.get_ylim()[1]), + ax.annotate('4', xy=(124, ax.get_ylim()[1]), fontsize=18, - xycoords='data', xytext=(0, 0), + color="blue", + xycoords='data', xytext=(-5, 0), textcoords='offset points') - ax.yaxis.label.set_size(22) - ax.xaxis.label.set_size(22) + + plt.grid() df.to_csv(self.input().path) # save plt.savefig(self.output().path, dpi=250, bbox_inches='tight', mode="pdf") class IPCAICreateOxyImageTask(luigi.Task): image_name = luigi.Parameter() df_prefix = luigi.Parameter() def requires(self): return IPCAITrainRegressor(df_prefix=self.df_prefix), \ Flatfield(flatfield_folder=sp.FLAT_FOLDER), \ SingleMultispectralImage(image=self.image_name), \ Dark(dark_folder=sp.DARK_FOLDER) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, self.image_name + "_" + self.df_prefix + "_summary" + ".png")) def run(self): df_results = pd.DataFrame() nrrd_reader = NrrdReader() tiff_ring_reader = TiffRingReader() # read the flatfield flat = nrrd_reader.read(self.input()[1].path) dark = nrrd_reader.read(self.input()[3].path) # read the msi nr_filters = len(sp.RECORDED_WAVELENGTHS) msi, segmentation = tiff_ring_reader.read(self.input()[2].path, nr_filters) # only take into account not saturated pixels. segmentation = np.logical_and(segmentation, (np.max(msi.get_image(), axis=-1) < 1000.)) # 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.image_correction(msi, flat, dark) # 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, time = 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.) + im = enh_brightness.enhance(10.) 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) 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.0 oxy_image[0, 1] = 1. # 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.)) + "%") df_image_results = pd.DataFrame(data=np.expand_dims([self.image_name, oxy_mean * 100., time], 0), columns=["image name", "oxygenation mean [%]", "time to estimate"]) results_file = os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "results.csv") if os.path.isfile(results_file): df_results = pd.read_csv(results_file, index_col=0) df_results = pd.concat((df_results, df_image_results)).reset_index( drop=True ) else: df_results = df_image_results df_results.to_csv(results_file) plt.savefig(self.output().path, dpi=250, 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): df_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "reg_" + self.df_prefix)) def requires(self): - return tasks_mc.CameraBatch(self.df_prefix) + return tasks_mc.SpectroCamBatch(self.df_prefix) def run(self): # extract data from the batch df_train = pd.read_csv(self.input().path, header=[0, 1]) df_snrs = pd.read_csv("/media/wirkert/data/Data/" + "2016_01_08_Calibration_Measurements/processed/" + "Carets_Laparoscope/finals/" + "snr_for_bands.txt", index_col=0) - X, y = preprocess(df_train, snr=df_snrs["SNR"].values, - movement_noise_sigma=0.1, + X, y = preprocess(df_train, snr=10., 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) # reg = LinearSaO2Unmixing() # save regressor regressor_file = self.output().open('w') pickle.dump(reg, regressor_file) regressor_file.close() class SingleMultispectralImage(luigi.Task): image = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, \ self.image)) class Flatfield(luigi.Task): flatfield_folder = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "flatfield" + ".nrrd")) def run(self): tiff_ring_reader = TiffRingReader() nr_filters = len(sp.RECORDED_WAVELENGTHS) # analyze all the first image files image_files = get_image_files_from_folder(self.flatfield_folder) image_files = filter(lambda image_name: "F0" in image_name, image_files) # helper function to take maximum of two images def maximum_of_two_images(image_1, image_name_2): image_2 = tiff_ring_reader.read(os.path.join(self.flatfield_folder, image_name_2), nr_filters)[0].get_image() return np.maximum(image_1, image_2) # now reduce to maximum of all the single images flat_maximum = reduce(lambda x, y: maximum_of_two_images(x, y), image_files, 0) msi = Msi(image=flat_maximum) msi.set_wavelengths(sp.RECORDED_WAVELENGTHS) # write flatfield as nrrd writer = NrrdWriter(msi) writer.write(self.output().path) class Dark(luigi.Task): dark_folder = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "dark" + ".nrrd")) def run(self): tiff_ring_reader = TiffRingReader() nr_filters = len(sp.RECORDED_WAVELENGTHS) # analyze all the first image files image_files = get_image_files_from_folder(self.dark_folder) image_files = filter(lambda image_name: "F0" in image_name, image_files) # returns the mean dark image vector of all inputted dark image # overly complicated TODO SW: make this simple code readable. dark_means = map(lambda image_name: msimani.calculate_mean_spectrum( tiff_ring_reader.read(os.path.join(self.dark_folder, image_name), nr_filters)[0]), image_files) dark_means_sum = reduce(lambda x, y: x+y.get_image(), dark_means, 0) final_dark_mean = dark_means_sum / len(dark_means) msi = Msi(image=final_dark_mean) msi.set_wavelengths(sp.RECORDED_WAVELENGTHS) # write flatfield as nrrd writer = NrrdWriter(msi) 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) image_file_folder = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER) onlyfiles = get_image_files_from_folder(image_file_folder) first_invivo_image_files = filter(lambda image_name: "F0" in image_name, onlyfiles) for f in first_invivo_image_files: main_task = IPCAICreateOxyImageTask(image_name=f, df_prefix="ipcai_revision_colon_mean_scattering_train") w.add(main_task) w.run() oxygenation_over_time_task = OxyOverTimeTask() w.add(oxygenation_over_time_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_analyze_neil_data.py b/Modules/Biophotonics/python/iMC/script_analyze_neil_data.py index b6ac569e5f..f90d793a2e 100644 --- a/Modules/Biophotonics/python/iMC/script_analyze_neil_data.py +++ b/Modules/Biophotonics/python/iMC/script_analyze_neil_data.py @@ -1,280 +1,280 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import Image import ImageEnhance import numpy as np import luigi import matplotlib.pyplot as plt import SimpleITK as sitk import tasks_analyze as at import tasks_preprocessing as pre import tasks_regression as rt 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 sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2014_12_08_Neil_Pig_Oxygenation_Experiments/organized_recordings" sp.DATA_FOLDER = "colon_images3" sp.FINALS_FOLDER = "Oxy3" def resort_wavelengths(msi): """Neil organized his wavelengths differently. This function fixes this.""" collapsed_image = imgmani.collapse_image(msi.get_image()) collapsed_image = collapsed_image[:, [1, 5, 3, 0, 6, 7, 2, 4]] msi.set_image(np.reshape(collapsed_image, msi.get_image().shape)) # rebind this method since in this recoding the wavelengths got messed up sp.resort_wavelengths = resort_wavelengths # the 700nm band got messed up in Neils data (way to dark). thus sort # it out from further processing. sp.bands_to_sortout = np.array([7]) class NeilCreateOxyImageTask(luigi.Task): image_prefix = luigi.Parameter() - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() def requires(self): return NeilEstimateTissueParametersTask(image_prefix=self.image_prefix, batch_prefix= - self.batch_prefix), \ + self.df_prefix), \ NeilCorrectImagingSetupTask(image_prefix=self.image_prefix), \ NeilReprojectReflectancesTask(image_prefix=self.image_prefix, 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.image_prefix + "_" + - self.batch_prefix + + self.df_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) f, axarr = plt.subplots(1, 3) # create artificial rgb rgb_image = msi.get_image()[:, :, [2, 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(2.) plot_image = np.array(im) top_left_axis = axarr[0] top_left_axis.imshow(plot_image, interpolation='nearest') top_left_axis.set_title("synthetic rgb", fontsize=5) top_left_axis.xaxis.set_visible(False) top_left_axis.yaxis.set_visible(False) plt.set_cmap("jet") # plot parametric maps image[image[:, :, 0] > 0.2] = 0.2 # set maximum bvf to 0.2 at.plot_axis(axarr[1], image[:, :, 0], "blood volume fraction [%]") image[0, 0, 1] = 0.; image[1, 0, 1] = 1. at.plot_axis(axarr[2], image[:, :, 1], "oxygenation [%]") plt.savefig(self.output().path, dpi=1000, bbox_inches='tight') class NeilEstimateTissueParametersTask(luigi.Task): image_prefix = luigi.Parameter() - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() def requires(self): return NeilPreprocessMSI(image_prefix=self.image_prefix), \ 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.image_prefix + "_" + - self.batch_prefix + + self.df_prefix + "estimate.nrrd")) def run(self): sitk_image = at.estimate_image(self.input()[0].path, self.input()[1].path) outFile = self.output().open('w') outFile.close() sitk.WriteImage(sitk_image, self.output().path) class NeilReprojectReflectancesTask(luigi.Task): image_prefix = luigi.Parameter() - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() def requires(self): return NeilEstimateTissueParametersTask(image_prefix=self.image_prefix, 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.image_prefix + "_backprojection_" + - self.batch_prefix + + self.df_prefix + "estimate.nrrd")) def run(self): sitk_image = at.estimate_image(self.input()[0].path, self.input()[1].path) outFile = self.output().open('w') outFile.close() sitk.WriteImage(sitk_image, self.output().path) class NeilPreprocessMSI(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 NeilCorrectImagingSetupTask(image_prefix=self.image_prefix) def run(self): reader = NrrdReader() image = reader.read(self.input().path) image.set_image(imgmani.sortout_bands(image.get_image(), sp.bands_to_sortout)) norm.standard_normalizer.normalize(image) pre.touch_and_save_msi(image, self.output()) class NeilCorrectImagingSetupTask(luigi.Task): """unfortunately filter were ordered differently in neils. this task is to do all the fiddeling to get it right again and corrects for dark image and flatfield""" image_prefix = luigi.Parameter() def requires(self): return MultiSpectralImageFromPNGFiles(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""" # 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) # correct for differently sorted filters if necessary resort_wavelengths(msi) 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, "TestImage_141211_10.19.22_")) writer = NrrdWriter(msi) 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) writer.write(self.output().path) class MultiSpectralImageFromPNGFiles(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 = PngReader() 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_F7_files = [ f for f in onlyfiles if f.endswith("_F7.png") ] files_prefixes = [ f[:-5] for f in only_F7_files] for f in files_prefixes: main_task = NeilCreateOxyImageTask( image_prefix=f, batch_prefix= "jacques_no_billi_generic_scattering_") 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 a32e77545c..19525fbd69 100644 --- a/Modules/Biophotonics/python/iMC/script_calculate_spectra.py +++ b/Modules/Biophotonics/python/iMC/script_calculate_spectra.py @@ -1,133 +1,133 @@ ''' 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_BATCHES = 1 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/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): 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.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_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='log/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, NR_BATCHES, 1) for i in BATCH_NUMBERS: - colon_task = CreateSpectraTask("ipcai_revision_colon_mean_scattering", + colon_task = CreateSpectraTask("ipcai_visualization_batch", i, NR_ELEMENTS_IN_BATCH, - mcfac.ColonMuscleMeanScatteringFactory()) + mcfac.VisualizationMcFactory()) w.add(colon_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_invertabilitytest.py b/Modules/Biophotonics/python/iMC/script_invertabilitytest.py index 357e99506f..f34e4d8bca 100644 --- a/Modules/Biophotonics/python/iMC/script_invertabilitytest.py +++ b/Modules/Biophotonics/python/iMC/script_invertabilitytest.py @@ -1,132 +1,132 @@ ''' Created on Aug 24, 2015 @author: wirkert ''' import os import pickle import matplotlib.pyplot as plt import matplotlib import numpy as np import luigi import tasks_regression as rt import scriptpaths as sp sp.FINALS_FOLDER = "invertabilitytest" class TestInvertability(luigi.Task): - batch_prefix = luigi.Parameter() + df_prefix = luigi.Parameter() file_name_prefix_testing = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "invertability_" + self.file_name_prefix_testing + ".png")) def requires(self): return rt.ParameterFile(self.file_name_prefix_testing), \ rt.ReflectanceFile(self.file_name_prefix_testing), \ rt.TrainForest(file_name_prefix_training= - self.batch_prefix), \ + self.df_prefix), \ rt.TrainForestAugmentedData(file_name_prefix_training= - self.batch_prefix) + self.df_prefix) def run(self): # get the parameters to test test_params, test_refl = rt.extract_mc_data(self.input()[0].path, self.input()[1].path) test_params = test_params[:, 0:4] # use regressors to calculate the errors # first random forest rf_regressor_file = self.input()[2].open('r') rf_regressor = pickle.load(rf_regressor_file) rf_absErrors = np.abs(rf_regressor.predict(test_refl) - test_params) rf_r2Score = rf_regressor.score(test_refl, test_params) # now augmented rf rf_aug_regressor_file = self.input()[3].open('r') rf_aug_regressor = pickle.load(rf_aug_regressor_file) rf_aug_absErrors = np.abs(rf_aug_regressor.predict(test_refl) - test_params) rf_aug_r2Score = rf_aug_regressor.score(test_refl, test_params) rf_r2Scores = [rf_r2Score, rf_aug_r2Score] # now for the plotting part colors = ['cyan', 'magenta'] matplotlib.rcParams.update({'font.size': 5}) f, axarr = plt.subplots(1, 4) bvf_axis = axarr[0] dsp_axis = axarr[1] r_axis = axarr[2] sao2_axis = axarr[3] boxplots = [] data = [rf_absErrors[:, 0] , rf_aug_absErrors[:, 0]] boxplots.append(bvf_axis.boxplot(data, 0, "", patch_artist=True)) bvf_axis.set_title("blood volume fraction") bvf_axis.set_ylabel('absolute error') bvf_axis.get_xaxis().set_ticks([]) data = [rf_absErrors[:, 1], rf_aug_absErrors[:, 1]] boxplots.append(dsp_axis.boxplot(data, 0, "", patch_artist=True)) dsp_axis.set_title("density scattering particles") dsp_axis.get_xaxis().set_ticks([]) data = [rf_absErrors[:, 2], rf_aug_absErrors[:, 2]] boxplots.append(r_axis.boxplot(data, 0, "", patch_artist=True)) r_axis.set_title("mucosa thickness") r_axis.get_xaxis().set_ticks([]) data = [rf_absErrors[:, 3], rf_aug_absErrors[:, 3]] boxplots.append(sao2_axis.boxplot(data, 0, "", patch_artist=True)) sao2_axis.set_title("oxygenation") sao2_axis.get_xaxis().set_ticks([]) for boxplot in boxplots: for i, (patch, color) in enumerate(zip(boxplot['boxes'], colors)): patch.set_facecolor(color) head, filename = os.path.split(self.input()[2 + i].path) plt.figtext(0.30, 0.08 - 0.02 * i, filename + ". rf_r2Score=" + str(rf_r2Scores[i]), backgroundcolor=colors[i], color='black', weight='roman', size='x-small') print("absolute error distribution BVF, Volume fraction") print("median: " + str(np.median(rf_absErrors, axis=0))) print("lower quartile: " + str(np.percentile(rf_absErrors, 25, axis=0))) print("higher quartile: " + str(np.percentile(rf_absErrors, 75, axis=0))) print("rf_r2Score", str(rf_r2Score)) plt.savefig(self.output().path, dpi=500) if __name__ == '__main__': # root folder there the data lies luigi.interface.setup_interface_logging() sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) main_task = TestInvertability( file_name_prefix_training= "2015June0809:51PMNoisyRandomTraining", file_name_prefix_testing= "2015June0902:43AMNoisyRandomTesting") w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_visualize_spectra.py b/Modules/Biophotonics/python/iMC/script_visualize_spectra.py index ce209f65c3..7b3b8974ff 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() + df_prefix = luigi.Parameter() batch_nr = luigi.IntParameter() nr_samples = luigi.IntParameter() def requires(self): - return tasks_mc.CreateSpectraTask(self.batch_prefix, self.batch_nr, + return tasks_mc.CreateSpectraTask(self.df_prefix, self.batch_nr, self.nr_samples, mcfac.VisualizationMcFactory()), \ - tasks_mc.CameraBatch(self.batch_prefix) + tasks_mc.CameraBatch(self.df_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")) + self.df_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 # 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/spielwiese.py b/Modules/Biophotonics/python/iMC/spielwiese.py index 85c9e64585..9db0b1a652 100644 --- a/Modules/Biophotonics/python/iMC/spielwiese.py +++ b/Modules/Biophotonics/python/iMC/spielwiese.py @@ -1,75 +1,43 @@ ''' Created on Oct 19, 2015 @author: wirkert ''' import pickle +import os import skimage.restoration as res import numpy as np import matplotlib.pylab as plt import SimpleITK as sitk from msi.io.nrrdreader import NrrdReader from msi.io.nrrdwriter import NrrdWriter import msi.msimanipulations as msimani import msi.plot as msiplot import msi.imgmani as imgmani from scipy.ndimage.filters import gaussian_filter def save_as_np(filename, prefix): f = open(filename) batch = pickle.load(f) np.save("/media/wirkert/data/Data/temp/" + prefix + "_reflectances", batch.reflectances) np.save("/media/wirkert/data/Data/temp/" + prefix + "_saO2", batch.layers[0][:, 1]) if __name__ == '__main__': - r = NrrdReader() + image_file_folder = "/media/wirkert/data/Data/2015_11_12_IPCAI_in_silico/mc_data2" - msi = r.read("/media/wirkert/data/Data/2015_11_12_IPCAI_in_vivo/processed/" + - "image sample .nrrd") - seg = r.read("/media/wirkert/data/Data/2015_11_12_IPCAI_in_vivo/" + - "noise_segmentation.nrrd") + onlyfiles = [ os.path.join(image_file_folder, f) for f in os.listdir(image_file_folder) if + os.path.isfile(os.path.join(image_file_folder, f)) ] + onlyfiles.sort() - image = msi.get_image() - -# sitk_reader = sitk.ImageFileReader() -# sitk_reader.SetFileName("/media/wirkert/data/Data/2015_11_12_IPCAI_in_vivo/processed/" + -# "image sample .nrrd") -# sitk_image = sitk_reader.Execute() -# -# -# gaussian = sitk.SmoothingRecursiveGaussianImageFilter() -# gaussian.SetSigma(100.) -# smooth_image = gaussian.Execute(sitk_image) - - filtered_image = gaussian_filter(image.astype(np.float), 2) - msi.set_image(filtered_image) - msimani.apply_segmentation(msi, seg) - msi.set_wavelengths(np.arange(470, 680, 10)) - msiplot.plotMeanError(msi) -# plt.show() -# -# image_res = res.denoise_tv_chambolle(image.astype(np.float), -# weight=1000, n_iter_max=100, -# multichannel=True) - - w = NrrdWriter(msi) - w.write("/media/wirkert/data/Data/2015_11_12_IPCAI_in_vivo/test.nrrd") - -# img = sitk.GetImageFromArray(smooth_image, isVector=True) -# sitk.WriteImage(img, "/media/wirkert/data/Data/2015_11_12_IPCAI_in_vivo/test.nrrd") - -# save_as_np("/media/wirkert/data/Data/" + -# "2015_11_12_IPCAI_in_silico/processed/" + -# "generic_tissue_train_all_camera.imc", -# "hard_train_2") -# -# save_as_np("/media/wirkert/data/Data/" + -# "2015_11_12_IPCAI_in_silico/processed/" + -# "generic_tissue_gaussian_all_camera.imc", -# "gaussian_train_2") + for f in onlyfiles: + file = open(f) + b = pickle.load(file) + print f + print np.sum(b.reflectances <= 0) + print np.sum(np.isnan(b.reflectances)) diff --git a/Modules/Biophotonics/python/iMC/tasks_preprocessing.py b/Modules/Biophotonics/python/iMC/tasks_preprocessing.py index bbad2d120f..317346e34a 100644 --- a/Modules/Biophotonics/python/iMC/tasks_preprocessing.py +++ b/Modules/Biophotonics/python/iMC/tasks_preprocessing.py @@ -1,270 +1,267 @@ ''' Created on Aug 26, 2015 @author: wirkert ''' import os import numpy as np import copy import luigi import SimpleITK as sitk 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 import scriptpaths as sp def calc_max_flatfield(flatfield1, flatfield2): """ calculate maximum of two flatfields """ flatfield = copy.copy(flatfield1) maximum_of_image_array = np.maximum(flatfield1.get_image(), flatfield2.get_image()) flatfield.set_image(maximum_of_image_array) return flatfield def smooth(msi): """ helper function to gaussian smooth msi channel by channel. """ img = sitk.GetImageFromArray(msi.get_image(), isVector=True) smoothFilter = sitk.SmoothingRecursiveGaussianImageFilter() smoothFilter.SetSigma(4) img_smoothed = smoothFilter.Execute(img) img_array = sitk.GetArrayFromImage(img_smoothed) msi.set_image(img_array) - - - def touch_and_save_msi(msi, outfile): """ saves msi as a nrrd to outfile. if the directory / file does not exist it will be created """ # touch file so path definately exists _outFile = outfile.open('w') _outFile.close() # use nrrd writer to write file _out = outfile writer = NrrdWriter(msi) writer.write(_out.path) class MultiSpectralImageFile(luigi.Task): """ the unaltered file c.f. hard disk """ imageName = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, self.imageName)) class ShortExposiorMSI(luigi.Task): """ the unaltered file c.f. hard disk """ imageName = luigi.Parameter() def output(self): assert(self.imageName.endswith("_long.nrrd")) short_exposior_imagename = self.imageName.replace("_long.nrrd", ".nrrd") return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, short_exposior_imagename)) class WhiteSegmentation(luigi.Task): """ the unaltered file c.f. hard disk """ imageName = luigi.Parameter() def output(self): assert(self.imageName.endswith("_long.nrrd")) white_seg_imagename = \ self.imageName.replace("_long.nrrd", "_white_segmentation.nrrd") return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, white_seg_imagename)) def get_ambient_factor(flatfield, short_exposior, segmentation): """calculate the changes to the flatfield introduced by the organ.""" flatfield_copy = copy.deepcopy(flatfield) se_copy = copy.deepcopy(short_exposior) msimani.apply_segmentation(flatfield_copy, segmentation) msimani.apply_segmentation(se_copy, segmentation) # calculate Fwb(x1,y1) msimani.calculate_mean_spectrum(se_copy) # calculate Fspace(x1,y1) msimani.calculate_mean_spectrum(flatfield_copy) # calculate Famb ambient_factor = se_copy.get_image() / flatfield_copy.get_image() return ambient_factor class CorrectImagingSetupTask(luigi.Task): """unfortunately filter were ordered weirdly. this task is to do all the fiddeling to get it right again""" imageName = luigi.Parameter() def requires(self): return MultiSpectralImageFile(imageName=self.imageName), \ PreprocessFlatfields(), \ PreprocessDark(), \ ShortExposiorMSI(imageName=self.imageName), \ WhiteSegmentation(imageName=self.imageName) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.imageName + "_image_setup_corrected.nrrd")) def run(self): """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) dark = reader.read(self.input()[2].path) se_msi = reader.read(self.input()[3].path) white_seg = reader.read(self.input()[4].path) # calculate and apply ambient factor f_ab = get_ambient_factor(flatfield, se_msi, white_seg) flatfield.set_image(flatfield.get_image() * f_ab) # add integration times for proper image correction if "long" in self.input()[0].path: msi.add_property({'integration times': np.array([150., 250., 117., 160., 150., 175., 82., 70.])}) else: msi.add_property( {'integration times': np.array([30., 50., 47., 32., 30., 35., 35., 60.])}) flatfield.add_property({'integration times': np.array([30., 50., 47., 32., 30., 35., 35., 60.])}) # correct image by flatfield and dark image msimani.image_correction(msi, flatfield, dark) # correct for differently sorted filters if necessary sp.resort_wavelengths(msi) touch_and_save_msi(msi, self.output()) class PreprocessFlatfields(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "flatfield.nrrd")) def requires(self): return luigi.Task.requires(self) def run(self): reader = NrrdReader() flatfield_folder = os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER) flatfield_nrrds = os.listdir(flatfield_folder) # map to full file path flatfield_nrrds = \ [os.path.join(flatfield_folder, f) for f in flatfield_nrrds] flatfields = map(reader.read, flatfield_nrrds) max_flatfield = reduce(calc_max_flatfield, flatfields) # apply sitk gaussian smoothing to flatfield result smooth(max_flatfield) touch_and_save_msi(max_flatfield, self.output()) class PreprocessDark(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "dark.nrrd")) def requires(self): return luigi.Task.requires(self) def run(self): reader = NrrdReader() dark_folder = os.path.join(sp.ROOT_FOLDER, sp.DARK_FOLDER) dark_nrrds = os.listdir(dark_folder) # map to full file path dark_nrrds = [os.path.join(dark_folder, d) for d in dark_nrrds] dark = reader.read(dark_nrrds[0]) # just take the first dark image # alternatively multiple dark images could be averaged to one. # apply sitk gaussian smoothing to dark result smooth(dark) touch_and_save_msi(dark, self.output()) class PreprocessMSI(luigi.Task): imageName = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.imageName + "_preprocessed.nrrd")) def requires(self): return CorrectImagingSetupTask(imageName=self.imageName) def run(self): reader = NrrdReader() image = reader.read(self.input().path) image.set_image(imgmani.sortout_bands(image.get_image(), sp.bands_to_sortout)) norm.standard_normalizer.normalize(image) touch_and_save_msi(image, self.output()) class SegmentMSI(luigi.Task): """ in this class we segment the msi. We filter out saturated (dark and bright) pixels. The remaining pixels are taken for domain adaptation """ imageName = luigi.Parameter() def requires(self): return MultiSpectralImageFile(self.imageName), \ PreprocessMSI(self.imageName) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.imageName + "_segmentation.nrrd")) def run(self): reader = NrrdReader() msi_image = reader.read(self.input()[0].path).get_image() preprocessed_msi_image = reader.read(self.input()[1].path).get_image() max_low_wavelengths = \ np.max(preprocessed_msi_image[:, :, [0, 1, 3, 4, 5]], axis=-1) min_high_wavelengths = \ np.min(preprocessed_msi_image[:, :, [2, 6, 7]], axis=-1) # do "blood test" segmentation = max_low_wavelengths < min_high_wavelengths # filter dark spots segmentation = np.logical_and(segmentation, (np.max(msi_image, axis=-1) > 400.)) # filter bright spots segmentation = np.logical_and(segmentation, (np.max(msi_image, axis=-1) < 4000.)) img = sitk.GetImageFromArray(np.uint8(segmentation)) outFile = self.output().open('w') outFile.close() sitk.WriteImage(img, self.output().path)