diff --git a/Modules/Biophotonics/python/iMC/mc/mcmanipulations.py b/Modules/Biophotonics/python/iMC/mc/mcmanipulations.py index c3e1376ed4..4e5a9346d0 100644 --- a/Modules/Biophotonics/python/iMC/mc/mcmanipulations.py +++ b/Modules/Biophotonics/python/iMC/mc/mcmanipulations.py @@ -1,51 +1,52 @@ ''' Created on Oct 19, 2015 @author: wirkert ''' import numpy as np from scipy.interpolate import interp1d def fold_by_sliding_average(batch, window_size): """take a batch and apply a sliding average with given window size to the reflectances. window_size is elements to the left and to the right. There will be some boundary effect on the edges.""" sliding_average = lambda x: np.convolve(x, np.ones(window_size, dtype=float) / float(window_size), mode="same") batch.reflectances = np.apply_along_axis(sliding_average, axis=1, arr=batch.reflectances) def interpolate_wavelengths(batch, new_wavelengths): """ interpolate image data to fit new_wavelengths. Current implementation performs simple linear interpolation. Neither existing nor new wavelengths need to be sorted. """ interpolator = interp1d(batch.wavelengths, - batch.reflectances, assume_sorted=False) + batch.reflectances, assume_sorted=False, + bounds_error=False) batch.reflectances = interpolator(new_wavelengths) batch.wavelengths = new_wavelengths def sortout_bands(batch, bands_to_sortout=None): """ TODO: Documentation and test """ if bands_to_sortout is not None: # get rid of sorted out bands batch.reflectances = np.delete(batch.reflectances, bands_to_sortout, axis=1) batch.wavelengths = np.delete(batch.wavelengths, bands_to_sortout) def select_n(batch, n): """ randomly select n elements from batch. TODO: Test """ perms = np.random.permutation(batch.reflectances.shape[0]) first_n_perms = perms[0:n] batch.reflectances = batch.reflectances[first_n_perms, :] for i, l in enumerate(batch.layers): batch.layers[i] = np.atleast_2d(l)[first_n_perms, :] - return batch \ No newline at end of file + return batch diff --git a/Modules/Biophotonics/python/iMC/msi/plot.py b/Modules/Biophotonics/python/iMC/msi/plot.py index fdb4d880a1..41d12fb2d0 100644 --- a/Modules/Biophotonics/python/iMC/msi/plot.py +++ b/Modules/Biophotonics/python/iMC/msi/plot.py @@ -1,84 +1,84 @@ # -*- coding: utf-8 -*- """ Created on Thu Aug 13 11:13:31 2015 @author: wirkert """ import copy import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.axes_grid1 import make_axes_locatable import imgmani import msimanipulations as msimani def plot(msi, axes=None, color=None): """ create a plot for the Msi with x axes being the wavelengths and y-axes being the corresponding image values (e.g. reflectances, absorptions) Takes image masks into account: doesn't plot a spectrum containing masked elements """ if axes is None: axes = plt.gca() sortedIndices = sorted(range(len(msi.get_wavelengths())), key=lambda k: msi.get_wavelengths()[k]) sortedWavelenghts = msi.get_wavelengths()[sortedIndices] # reshape to collapse all but last dimension (which contains reflectances) collapsedImage = imgmani.collapse_image(msi.get_image()) # todo: simply use np.ma.compress_rows # print "filtered ", filteredImage.shape i = 0 for i in range(collapsedImage.shape[0]): if (collapsedImage[i, 0] is not np.ma.masked): axes.plot(sortedWavelenghts, collapsedImage[i, :][sortedIndices], "-o", color=color) def plot_images(msi): """plot the images as a 2d image array, one image for each wavelength.""" nr_wavelengths = msi.get_image().shape[-1] f, axarr = plt.subplots(1, nr_wavelengths) for i, a in enumerate(axarr): one_band_image = imgmani.get_bands(msi.get_image(), i) im = a.imshow(np.squeeze(one_band_image)) a.set_title("band nr " + str(i), fontsize=5) divider_dsp = make_axes_locatable(a) cax_dsp = divider_dsp.append_axes("right", size="10%", pad=0.1) cbar = plt.colorbar(im, cax=cax_dsp) cbar.ax.tick_params(labelsize=5) a.xaxis.set_visible(False) a.yaxis.set_visible(False) def plotMeanError(msi, axes=None): """ create a plot for the Msi with x axes being the wavelengths and y-axes being the corresponding mean image values (e.g. reflectances, absorptions). Plots also standard deviation bands Takes image masks into account: doesn't plot a spectrum containing masked elements """ if axes is None: axes = plt.gca() # sort the wavelengths sortedIndices = sorted(range(len(msi.get_wavelengths())), key=lambda k: msi.get_wavelengths()[k]) sortedWavelenghts = msi.get_wavelengths()[sortedIndices] # copy the msi, since it will be altered (mean will be built) msi_copy = copy.deepcopy(msi) image = msi_copy.get_image() image = imgmani.collapse_image(image) std_curve = np.ma.std(image, axis=0) msimani.calculate_mean_spectrum(msi_copy) # calculate std - print "percentual std: " + str(std_curve / msi_copy.get_image()) + print "percentual std: " + str(std_curve / msi_copy.get_image() * 100.) # plot as errorbar axes.errorbar(sortedWavelenghts, msi_copy.get_image()[sortedIndices], yerr=std_curve, fmt='-o') diff --git a/Modules/Biophotonics/python/iMC/regression/preprocessing.py b/Modules/Biophotonics/python/iMC/regression/preprocessing.py index fa6bf4aabd..6a27afedb4 100644 --- a/Modules/Biophotonics/python/iMC/regression/preprocessing.py +++ b/Modules/Biophotonics/python/iMC/regression/preprocessing.py @@ -1,65 +1,74 @@ ''' Created on Oct 26, 2015 @author: wirkert ''' import copy import numpy as np from sklearn.preprocessing import Normalizer import mc.mcmanipulations as mcmani -def preprocess(batch, nr_samples=None, w_percent=None, magnification=None, +def preprocess2(batch, nr_samples=None, w_percent=None, magnification=None, bands_to_sortout=None): + working_batch = copy.deepcopy(batch) mcmani.sortout_bands(working_batch, bands_to_sortout) # get reflectance and oxygenation X = working_batch.reflectances - y = working_batch.layers[0][:, 1] - y = y[:, np.newaxis] + y = working_batch.layers[0][:, [0, 1]] # remove nan no_nan = ~np.isnan(X).any(axis=1) X = X[no_nan, :] y = y[no_nan, :] # remove zero no_zero = ~(X <= 0).any(axis=1) X = X[no_zero, :] y = y[no_zero, :] # extract nr_samples samples from data if nr_samples is not None: X = X[0:nr_samples] y = y[0:nr_samples] # do data magnification if magnification is not None: X_temp = X y_temp = y for i in range(magnification - 1): X = np.vstack((X, X_temp)) y = np.vstack((y, y_temp)) # add noise to reflectances if w_percent is not None and (w_percent > 0.): noises = np.random.normal(loc=0., scale=w_percent, size=X.shape) X += noises * X X = np.clip(X, 0.00001, 1.) # do normalizations X = normalize(X) return X, np.squeeze(y) + +def preprocess(batch, nr_samples=None, w_percent=None, magnification=None, + bands_to_sortout=None): + X, y = preprocess2(batch, nr_samples, w_percent, magnification, + bands_to_sortout) + + return X, y[:, 1] + + 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_ipcai_in_silico.py b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_silico.py index 36cdf82027..e0f0a2d9a0 100644 --- a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_silico.py +++ b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_silico.py @@ -1,543 +1,637 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import pickle from collections import namedtuple -import itertools -import time - import luigi import matplotlib.pyplot as plt -from sklearn.neighbors import KNeighborsRegressor +import matplotlib from sklearn.ensemble.forest import RandomForestRegressor -import sklearn.svm -import sklearn.grid_search -import sklearn.cross_validation -from sklearn.linear_model import * -from sklearn.decomposition import PCA -import scipy.io import tasks_mc import scriptpaths as sp -from regression.preprocessing import preprocess +from regression.preprocessing import preprocess, preprocess2 from regression.linear import LinearSaO2Unmixing from regression.domain_adaptation import * from regression.estimation import standard_score sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2015_11_12_IPCAI_in_silico" sp.FINALS_FOLDER = "Images" sp.RECORDED_WAVELENGTHS = np.arange(470, 680, 10) * 10 ** -9 sp.MC_DATA_FOLDER = "mc_data2" w_standard = 0.1 # for this evaluation we add 5% noise +matplotlib.rcParams.update({'font.size': 22}) EvaluationStruct = namedtuple("EvaluationStruct", "name regressor data color alpha") +# setup standard random forest +rf = RandomForestRegressor(10, min_samples_leaf=10, max_depth=9, n_jobs=-1) + class TrainingSamplePlots(luigi.Task): which = luigi.Parameter() def requires(self): return tasks_mc.CameraBatch("ipcai_" + self.which + "_train"), \ tasks_mc.CameraBatch("ipcai_" + self.which + "_test") def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "samples_plot_" + self.which + ".png")) def run(self): # get data f = open(self.input()[0].path, "r") batch_train = pickle.load(f) f = open(self.input()[1].path, "r") batch_test = pickle.load(f) # setup testing function and data X_test, y_test = preprocess(batch_test, w_percent=w_standard) X_train, y_train = preprocess(batch_train, w_percent=w_standard) # setup structures to save test data evaluation_setups = [ - EvaluationStruct("Linear Beer-Lambert", LinearSaO2Unmixing(), - [], "red", 0.5), - EvaluationStruct("Proposed", RandomForestRegressor(10, - max_depth=10, - min_samples_leaf=1, - n_jobs=-1), [], "green", 0.5)] +# EvaluationStruct("Linear Beer-Lambert", +# LinearSaO2Unmixing(), [], "red", 0.5), + EvaluationStruct("Proposed", rf, "green", 0.5)] # do validation - nr_training_samples = np.arange(10, 10010, 50).astype(int) + nr_training_samples = np.arange(10, 15060, 50).astype(int) # not very pythonic, don't care for n in nr_training_samples: print "number samples ", str(n) for e in evaluation_setups: X_j = X_train[0:n] y_j = y_train[0:n] lr = e.regressor lr.fit(X_j, y_j) # save results e.data.append(standard_score(lr, X_test, y_test) * 100.) # make a nice plot for results plt.figure() legend_handles = [] for e in evaluation_setups: h, = plt.plot(nr_training_samples, e.data, color=e.color, label=e.name) legend_handles.append(h) # minor_ticks_x = np.arange(0, 1000, 50) # minor_ticks_y = np.arange(5, 16, 0.5) # plt.gca().set_xticks(minor_ticks_x, minor=True) # plt.gca().set_yticks(minor_ticks_y, minor=True) plt.grid() - plt.legend(handles=legend_handles) - plt.title("dependency on number of training samples") + # plt.legend(handles=legend_handles) + # plt.title("dependency on number of training samples") plt.xlabel("number of training samples") plt.ylabel("absolute error [%]") - plt.ylim((4, 18)) + + for tick in plt.gca().xaxis.get_major_ticks(): + tick.label.set_fontsize(14) + plt.ylim((2, 10)) + plt.xlim((0, 15000)) plt.savefig(self.output().path, dpi=500, bbox_inches='tight') class WrongNoisePlots(luigi.Task): which = luigi.Parameter() def requires(self): return tasks_mc.CameraBatch("ipcai_" + self.which + "_train"), \ tasks_mc.CameraBatch("ipcai_" + self.which + "_test") def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "noise_wrong_plot_" + self.which + ".png")) def run(self): # get data f = open(self.input()[0].path, "r") batch_train = pickle.load(f) f = open(self.input()[1].path, "r") batch_test = pickle.load(f) - # setup structures to save test data - rf = RandomForestRegressor(10, max_depth=10, n_jobs=-1) - evaluation_setups = [ - EvaluationStruct("Linear Beer-Lambert", LinearSaO2Unmixing(), [], - "red", 0.25) + EvaluationStruct("Linear Beer-Lambert", + LinearSaO2Unmixing(), [], "red", 0.25) , EvaluationStruct("Proposed", rf, [], "green", 0.5) ] # iterate over different levels of noise - noise_levels = np.arange(0.00, 1.00, 0.02) + noise_levels = np.arange(0.00, 0.17, 0.02) for w in noise_levels: # setup testing function X_test, y_test = preprocess(batch_test, w_percent=w) # extract noisy data X_train, y_train = preprocess(batch_train, w_percent=0.1) for e in evaluation_setups: lr = e.regressor lr.fit(X_train, y_train) y_pred = lr.predict(X_test) # save results errors = np.abs(y_pred - y_test) e.data.append(errors) # make a nice plot for results plt.figure() legend_handles = [] median0 = 0 for e in evaluation_setups: p25 = lambda x: np.percentile(x, 25) * 100. p75 = lambda x: np.percentile(x, 75) * 100. m = lambda x: np.median(x) * 100. plt.fill_between(np.array(noise_levels) * 100., map(p25, e.data), map(p75, e.data), edgecolor=e.color, facecolor=e.color, alpha=e.alpha) h, = plt.plot(np.array(noise_levels) * 100., map(m, e.data), color=e.color, label=e.name) legend_handles.append(h) if e.name is "Proposed": median0 = m(e.data[0]) - minor_ticks = np.arange(0, 80, 5) + minor_ticks = np.arange(0, 15, 5) plt.gca().set_xticks(minor_ticks, minor=True) plt.gca().set_yticks(minor_ticks, minor=True) plt.gca().annotate('error: ' + "{0:.2f}".format( np.median(evaluation_setups[-1].data[0]) * 100) + "%", - xy=(0, median0), xytext=(10, 35), size=10, + xy=(0, median0), xytext=(5, 10), size=15, va="center", ha="center", bbox=dict(boxstyle="round4", fc="g", alpha=0.5), arrowprops=dict(arrowstyle="-|>", connectionstyle="arc3,rad=-0.2", fc="w"), ) plt.grid() plt.legend(handles=legend_handles) - plt.title("dependency on wrong training noise [w_training = 10%]") + # plt.title("dependency on wrong training noise [w_training = 10%]") plt.xlabel("noise added for testing [sigma %]") plt.ylabel("absolute error [%]") - plt.ylim((0, 80)) + plt.ylim((0, 20)) + plt.xlim((0, 15)) + plt.savefig(self.output().path, dpi=500, + bbox_inches='tight') + + + +class NoisePlotsBVF(luigi.Task): + + which = luigi.Parameter() + + def requires(self): + return tasks_mc.CameraBatch("ipcai_" + self.which + "_train"), \ + tasks_mc.CameraBatch("ipcai_" + self.which + "_test") + + def output(self): + return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, + sp.RESULTS_FOLDER, + sp.FINALS_FOLDER, + "noise_plot_bvf_" + self.which + ".png")) + + def run(self): + # get data + f = open(self.input()[0].path, "r") + batch_train = pickle.load(f) + f = open(self.input()[1].path, "r") + batch_test = pickle.load(f) + + evaluation_setups = [ + EvaluationStruct("Proposed", rf, [], "green", 0.5) + ] + # iterate over different levels of noise + noise_levels = np.arange(0.00, 0.17, 0.02) + for w in noise_levels: + # setup testing function + X_test, y_test = preprocess2(batch_test, w_percent=w) + # extract noisy data + X_train, y_train = preprocess2(batch_train, w_percent=w) + for e in evaluation_setups: + lr = e.regressor + lr.fit(X_train, y_train) + y_pred = lr.predict(X_test) + # save results + errors = np.abs(y_pred[:, 0] - y_test[:, 0]) + e.data.append(errors) + # make a nice plot for results + plt.figure() + legend_handles = [] + for e in evaluation_setups: + p25 = lambda x: np.percentile(x, 25) * 100. + p75 = lambda x: np.percentile(x, 75) * 100. + m = lambda x: np.median(x) * 100. + plt.fill_between(np.array(noise_levels) * 100., + map(p25, e.data), + map(p75, e.data), + edgecolor=e.color, facecolor=e.color, + alpha=e.alpha) + h, = plt.plot(np.array(noise_levels) * 100., + map(m, e.data), color=e.color, + label=e.name) + legend_handles.append(h) + minor_ticks = np.arange(0, 15, 5) + plt.gca().set_xticks(minor_ticks, minor=True) + + plt.grid() + # plt.legend(handles=legend_handles) + # plt.title("dependency on noise") + plt.xlabel("noise added [sigma %]") + plt.ylabel("absolute error [%]") + plt.ylim((0, 4)) + plt.xlim((0, 15)) + plt.savefig(self.output().path, dpi=500, + bbox_inches='tight') + + + +class DomainNoisePlot(luigi.Task): + + def requires(self): + return tasks_mc.CameraBatch("ipcai_colon_muscle_train"), \ + tasks_mc.CameraBatch("ipcai_generic_test") + + def output(self): + return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, + sp.RESULTS_FOLDER, + sp.FINALS_FOLDER, + "noise_da_plot_.png")) + + def run(self): + # get data + f = open(self.input()[0].path, "r") + batch_train = pickle.load(f) + f = open(self.input()[1].path, "r") + batch_test = pickle.load(f) + + evaluation_setups = [ + EvaluationStruct("Linear Beer-Lambert", + LinearSaO2Unmixing(), [], "red", 0.25) + , EvaluationStruct("Proposed", rf, [], "green", 0.5) + ] + # iterate over different levels of noise + noise_levels = np.arange(0.00, 0.17, 0.02) + for w in noise_levels: + # setup testing function + X_test, y_test = preprocess(batch_test, w_percent=w) + # extract noisy data + X_train, y_train = preprocess(batch_train, w_percent=w) + for e in evaluation_setups: + lr = e.regressor + lr.fit(X_train, y_train) + y_pred = lr.predict(X_test) + # save results + errors = np.abs(y_pred - y_test) + e.data.append(errors) + # make a nice plot for results + plt.figure() + legend_handles = [] + for e in evaluation_setups: + p25 = lambda x: np.percentile(x, 25) * 100. + p75 = lambda x: np.percentile(x, 75) * 100. + m = lambda x: np.median(x) * 100. + plt.fill_between(np.array(noise_levels) * 100., + map(p25, e.data), + map(p75, e.data), + edgecolor=e.color, facecolor=e.color, + alpha=e.alpha) + h, = plt.plot(np.array(noise_levels) * 100., + map(m, e.data), color=e.color, + label=e.name) + legend_handles.append(h) + if e.name is "Proposed": + median0 = m(e.data[0]) + minor_ticks = np.arange(0, 15, 5) + plt.gca().set_xticks(minor_ticks, minor=True) + plt.gca().set_yticks(minor_ticks, minor=True) + + plt.gca().annotate('error: ' + + "{0:.2f}".format( + np.median(evaluation_setups[-1].data[0]) + * 100) + "%", + xy=(0, median0), xytext=(5, 12), size=15, + va="center", ha="center", bbox=dict( + boxstyle="round4", fc="g", alpha=0.5), + arrowprops=dict(arrowstyle="-|>", + connectionstyle="arc3,rad=-0.2", + fc="w"), + ) + plt.grid() + plt.legend(handles=legend_handles) + + # plt.title("dependency on noise") + plt.xlabel("noise added [sigma %]") + plt.ylabel("absolute error [%]") + plt.ylim((0, 20)) + plt.xlim((0, 15)) plt.savefig(self.output().path, dpi=500, bbox_inches='tight') class NoisePlots(luigi.Task): which = luigi.Parameter() def requires(self): return tasks_mc.CameraBatch("ipcai_" + self.which + "_train"), \ tasks_mc.CameraBatch("ipcai_" + self.which + "_test") def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "noise_plot_" + self.which + ".png")) def run(self): # get data f = open(self.input()[0].path, "r") batch_train = pickle.load(f) f = open(self.input()[1].path, "r") batch_test = pickle.load(f) - # setup structures to save test data - rf = RandomForestRegressor(10, max_depth=10, n_jobs=-1) - evaluation_setups = [ EvaluationStruct("Linear Beer-Lambert", LinearSaO2Unmixing(), [], "red", 0.25) , EvaluationStruct("Proposed", rf, [], "green", 0.5) ] # iterate over different levels of noise - noise_levels = np.arange(0.00, 1.02, 0.02) + noise_levels = np.arange(0.00, 0.17, 0.02) for w in noise_levels: # setup testing function X_test, y_test = preprocess(batch_test, w_percent=w) # extract noisy data X_train, y_train = preprocess(batch_train, w_percent=w) for e in evaluation_setups: lr = e.regressor lr.fit(X_train, y_train) y_pred = lr.predict(X_test) # save results errors = np.abs(y_pred - y_test) e.data.append(errors) # make a nice plot for results plt.figure() legend_handles = [] for e in evaluation_setups: p25 = lambda x: np.percentile(x, 25) * 100. p75 = lambda x: np.percentile(x, 75) * 100. m = lambda x: np.median(x) * 100. plt.fill_between(np.array(noise_levels) * 100., map(p25, e.data), map(p75, e.data), edgecolor=e.color, facecolor=e.color, alpha=e.alpha) h, = plt.plot(np.array(noise_levels) * 100., map(m, e.data), color=e.color, label=e.name) legend_handles.append(h) if e.name is "Proposed": median0 = m(e.data[0]) - minor_ticks = np.arange(0, 80, 5) + minor_ticks = np.arange(0, 15, 5) plt.gca().set_xticks(minor_ticks, minor=True) plt.gca().set_yticks(minor_ticks, minor=True) plt.gca().annotate('error: ' + "{0:.2f}".format( np.median(evaluation_setups[-1].data[0]) * 100) + "%", - xy=(0, median0), xytext=(10, 35), size=10, va="center", + xy=(0, median0), xytext=(5, 10), size=15, va="center", ha="center", bbox=dict(boxstyle="round4", fc="g", alpha=0.5), arrowprops=dict(arrowstyle="-|>", connectionstyle="arc3,rad=-0.2", fc="w"), ) plt.grid() plt.legend(handles=legend_handles) - plt.title("dependency on noise") + # plt.title("dependency on noise") plt.xlabel("noise added [sigma %]") plt.ylabel("absolute error [%]") - plt.ylim((0, 80)) + plt.ylim((0, 20)) + plt.xlim((0, 15)) plt.savefig(self.output().path, dpi=500, bbox_inches='tight') class WeightsCalculation(luigi.Task): + source_domain = luigi.Parameter() + target_domain = luigi.Parameter() def requires(self): - return tasks_mc.CameraBatch("generic_tissue_no_aray_train"), \ - tasks_mc.CameraBatch("colon_muscle_tissue_d_ranges_train") + return tasks_mc.CameraBatch(self.source_domain), \ + tasks_mc.CameraBatch(self.target_domain) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "weights.npy")) def run(self): # get data f = open(self.input()[0].path, "r") batch_train = pickle.load(f) f = open(self.input()[1].path, "r") batch_da_train = pickle.load(f) X_train, y_train = preprocess(batch_train, w_percent=w_standard, - magnification=10) + magnification=None) X_da_train, y_da_train = preprocess(batch_da_train, - w_percent=w_standard, magnification=10) + w_percent=w_standard, + magnification=None) X_ws, y_ws = preprocess(batch_train, w_percent=w_standard) weights = estimate_weights_random_forests(X_train, X_da_train, X_ws) - bvfs = batch_train.layers[0][:, 0] - saO2s = batch_train.layers[0][:, 1] - a_mies = batch_train.layers[0][:, 2] - ds = batch_train.layers[0][:, 4] + + no_zero = ~(batch_train.reflectances <= 0).any(axis=1) + bvfs = batch_train.layers[0][no_zero, 0] + saO2s = batch_train.layers[0][no_zero, 1] + a_mies = batch_train.layers[0][no_zero, 2] + ds = batch_train.layers[0][no_zero, 4] correct_high_weights = np.logical_and(bvfs < 0.10, ds < 1000 * 10 ** -6, ds > 600 * 10 ** -6) plt.figure() plt.hist(bvfs, weights=weights, alpha=0.5, label="found_bvfs") plt.hist(bvfs[correct_high_weights], alpha=0.5, label="true_bvfs") plt.title("weight distr bvf") plt.legend(loc='upper right') plt.figure() plt.hist(saO2s, weights=weights, alpha=0.5, label="found_saO2s") plt.hist(saO2s[correct_high_weights], alpha=0.5, label="true_saO2s") plt.title("weight distr saO2") plt.legend(loc='upper right') plt.figure() plt.hist(a_mies, weights=weights, alpha=0.5, label="found_a_mies") plt.hist(a_mies[correct_high_weights], alpha=0.5, label="true_a_mies") plt.title("weight distr a_mie") plt.legend(loc='upper right') plt.figure() plt.hist(ds, weights=weights, alpha=0.5, label="found_ds") plt.hist(ds[correct_high_weights], alpha=0.5, label="true_ds") plt.title("weight distr d") plt.legend(loc='upper right') plt.show() np.save(self.output().path, weights) -class DAPlots(luigi.Task): - +class NewDAPlot(luigi.Task): def requires(self): - return tasks_mc.CameraBatch("generic_tissue_no_aray_train"), \ - tasks_mc.CameraBatch("colon_muscle_tissue_d_ranges_train"), \ - tasks_mc.CameraBatch("colon_muscle_tissue_d_ranges_train"), \ - WeightsCalculation() + return tasks_mc.CameraBatch("ipcai_generic"), \ + tasks_mc.CameraBatch("ipcai_colon_muscle_test"), \ + WeightsCalculation("ipcai_generic", "ipcai_colon_muscle_test") def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, - "covariance_shift_plot.png")) + "new_da.png")) def run(self): # get data f = open(self.input()[0].path, "r") batch_train = pickle.load(f) f = open(self.input()[1].path, "r") - batch_da_train = pickle.load(f) + batch_test = pickle.load(f) f = open(self.input()[2].path, "r") - batch_da_test = pickle.load(f) - weights = np.load(self.input()[3].path) - X_train, y_train = preprocess(batch_train, w_percent=w_standard) - X_da_train, y_da_train = preprocess(batch_da_train, - w_percent=w_standard) - X_da_test, y_da_test = preprocess(batch_da_test, - w_percent=w_standard) - scipy.io.savemat('/home/wirkert/workspace/matlab/domain_adaptation/GMKLIEP_demo/x_nu', mdict={'arr': X_da_train}) - scipy.io.savemat('/home/wirkert/workspace/matlab/domain_adaptation/GMKLIEP_demo/x_de', mdict={'arr': X_train}) - scipy.io.savemat('/home/wirkert/workspace/matlab/domain_adaptation/GMKLIEP_demo/bvf_de', mdict={'arr': batch_train.layers[0][:, 0]}) - # setup structures to save test data - kf = sklearn.cross_validation.KFold(X_train.shape[0], 10, shuffle=True) - param_grid_svm = [ - {'C': np.logspace(1, 8, 8), - 'kernel':['rbf'], - 'gamma': np.logspace(-3, 3, 7), - 'nu':np.array([0.2, 0.5, 0.7])} ] - svm = sklearn.grid_search.GridSearchCV(sklearn.svm.NuSVR(), - param_grid_svm, cv=kf, n_jobs=-1) - svm = sklearn.svm.SVR(kernel="rbf", degree=3, C=100., gamma=10.) - # use the random forest regressor - param_grid_rf = [ - {"n_estimators": np.array([50]), - "max_depth": np.array([30]), # np.arange(5, 40, 2), - "max_features": np.array([8]), # np.array(range(X_train.shape[1] - 1)) + 1, - "min_weight_fraction_leaf": np.array([0.01, 0.03, 0.05, 0.07, 0.1, 0.2])}] # , 7, 10, 15])}] - rf = sklearn.grid_search.GridSearchCV(RandomForestRegressor(50, - max_depth=15, n_jobs=-1), - param_grid_rf, cv=kf, n_jobs=-1) - # do not do cv for now to save some time. - rf = RandomForestRegressor(10, max_depth=10, # min_weight_fraction_leaf=0.05, - n_jobs=-1) - + weights = np.load(f) evaluation_setups = [ - EvaluationStruct("classic", LinearSaO2Unmixing(), [], - "red", 0.5) - , EvaluationStruct("lr", LinearRegression(), [], - "yellow", 0.5) - , EvaluationStruct("knn", - KNeighborsRegressor(), - [], "yellow", 0.5) - , EvaluationStruct("rf", rf, [], "green", 0.5) - , EvaluationStruct("svr", svm, [], "yellow", 0.5) -# , EvaluationStruct("weighted rf", rf, [], "green", 0.5) - , EvaluationStruct("weighted svm", svm, [], "yellow", 0.5) - , EvaluationStruct("svm no cov shift", svm, [], "blue", 0.5) + EvaluationStruct("Linear Beer-Lambert", + LinearSaO2Unmixing(), [], "red", 0.25) + , EvaluationStruct("Proposed", rf, [], "yellow", 0.5) + , EvaluationStruct("Proposed DA", rf, [], "green", 0.5) ] - - bvfs = batch_train.layers[0][:, 0] - saO2s = batch_train.layers[0][:, 1] - a_mies = batch_train.layers[0][:, 2] - ds = batch_train.layers[0][:, -1] - -# weights = scipy.io.loadmat("/media/wirkert/data/Data/2015_11_12_IPCAI_in_silico/weights_LHSS.mat") -# weights = weights['r'] -# weights = scipy.io.loadmat("/media/wirkert/data/Data/2015_11_12_IPCAI_in_silico/weights.mat") -# weights = weights['wh_x_de'] -# weights = np.squeeze(weights) - + X_train, y_train = preprocess(batch_train, w_percent=w_standard) + X_da_test, y_da_test = preprocess(batch_test, w_percent=w_standard) + X_train_sampled, y_train_sampled, weights_sampled = \ + resample(X_train, y_train, weights, 8000) f, axarr = plt.subplots(X_train.shape[1], 1) - -# X_train, y_train, weights = resample(X_train, y_train, weights, 5000) -# pca = PCA(copy=True, n_components=5, whiten=True) -# pca.fit(X_train) -# X_train = pca.transform(X_train) -# X_da_train = pca.transform(X_da_train) -# X_da_test = pca.transform(X_da_test) for i in range(X_train.shape[1]): y, binEdges = np.histogram(X_train[:, i], weights=weights) bincenters = 0.5 * (binEdges[1:] + binEdges[:-1]) axarr[i].plot(bincenters, y, '-+', label="adapted") y, binEdges = np.histogram(X_train[:, i]) bincenters = 0.5 * (binEdges[1:] + binEdges[:-1]) axarr[i].plot(bincenters, y, '-o', label="original") - y, binEdges = np.histogram(X_da_train[:, i]) + y, binEdges = np.histogram(X_da_test[:, i], + weights=np.ones(X_da_test.shape[0]) * + X_train.shape[0] / X_da_test.shape[0]) bincenters = 0.5 * (binEdges[1:] + binEdges[:-1]) axarr[i].plot(bincenters, y, '-*', label="searched") + y, binEdges = np.histogram(X_train_sampled[:, i], + weights=np.ones(X_train_sampled.shape[0]) * + X_train.shape[0] / X_train_sampled.shape[0]) + bincenters = 0.5 * (binEdges[1:] + binEdges[:-1]) + axarr[i].plot(bincenters, y, '-^', label="sampled") + + axarr[i].xaxis.set_visible(False) + axarr[i].yaxis.set_visible(False) axarr[i].set_title("feature " + str(i) + " distribution") + plt.legend(loc="best") # plt.tight_layout(h_pad=.5) -# plt.show() - -# for i in range(10): -# relevant_elements = np.logical_and(bvfs > 0.01 * i, -# bvfs < 0.01 * (i + 1)) -# relevant_saO2s = saO2s[relevant_elements] -# relevant_weights = weights[relevant_elements] -# plt.figure() -# plt.hist(relevant_saO2s) # , weights=relevant_weights) -# plt.savefig(self.output().path + "saO2_distribution_" + str(i) + -# ".png") - + plt.show() - for e in evaluation_setups: - # train - start = time.time() - if "no cov shift" in e.name: - e.regressor.fit(X_da_train[0:-1], y_da_train[0:-1]) - elif "weighted" in e.name: - # print "number of used samples", np.sum(weights > 0.4) - # rw = weights > 1. - # e.regressor.fit(X_train[rw, :], y_train[rw]) - e.regressor.fit(X_train, y_train, weights) - # tree.export_graphviz(rf.estimators_[0], out_file='tree.dot') - else: - e.regressor.fit(X_train, y_train) - # print some interesting statistics - end = time.time() - print "time to train", e.name, ":", str(end - start), "s" - if isinstance(e.regressor, sklearn.grid_search.GridSearchCV): - print e.regressor.best_estimator_ - - # evaluate - start = time.time() - y_pred = e.regressor.predict(X_da_test) - end = time.time() - print "time to evaluate ", str(X_da_test.shape[0]), "samples: ", \ - str(end - start), "s" - - # save results - e.data.append(np.abs(y_pred - y_da_test)) - # collect data for boxplot - boxplot_data = [] - for e in evaluation_setups: - boxplot_data.append(100. * e.data[0]) + # iterate over different levels of noise + noise_levels = np.arange(0.00, 0.17, 0.02) + for w in noise_levels: + # setup testing function + X_test, y_test = preprocess(batch_test, w_percent=w) + # extract noisy data + X_train, y_train = preprocess(batch_train, w_percent=w) + for e in evaluation_setups: + lr = e.regressor + if "DA" in e.name: + lr.fit(X_train, y_train, weights) + else: + lr.fit(X_train, y_train) + y_pred = lr.predict(X_test) + # save results + errors = np.abs(y_pred - y_test) + e.data.append(errors) # make a nice plot for results - fig = plt.figure() - ax = fig.add_subplot(111) - bp = ax.boxplot(boxplot_data, showfliers=False, - showmeans=True, patch_artist=True) - # nicen it up - for e, box in zip(evaluation_setups, bp['boxes']): - # change outline color - box.set(color="black", linewidth=2) - # change fill color - box.set(facecolor=e.color, alpha=e.alpha) - # change color and linewidth of the whiskers - x_tick_labels = [] - for e, whisker, cap, median in itertools.izip( - evaluation_setups, - bp['whiskers'], bp['caps'], bp['medians']): - whisker.set(color="black", linewidth=2) - cap.set(color="black", linewidth=2) - median.set(color="black", linewidth=2) - x_tick_labels.append(e.name) - ax.set_xticklabels(x_tick_labels) - ax.get_xaxis().tick_bottom() - ax.get_yaxis().tick_left() - ax.axvline(x=len(evaluation_setups) - 1.5, - linestyle="--", color="black", linewidth=2) + plt.figure() + legend_handles = [] + for e in evaluation_setups: + p25 = lambda x: np.percentile(x, 25) * 100. + p75 = lambda x: np.percentile(x, 75) * 100. + m = lambda x: np.median(x) * 100. + plt.fill_between(np.array(noise_levels) * 100., + map(p25, e.data), + map(p75, e.data), + edgecolor=e.color, facecolor=e.color, + alpha=e.alpha) + h, = plt.plot(np.array(noise_levels) * 100., + map(m, e.data), color=e.color, + label=e.name) + legend_handles.append(h) + if e.name is "Proposed": + median0 = m(e.data[0]) + minor_ticks = np.arange(0, 15, 5) + plt.gca().set_xticks(minor_ticks, minor=True) + plt.gca().set_yticks(minor_ticks, minor=True) + + plt.gca().annotate('error: ' + + "{0:.2f}".format( + np.median(evaluation_setups[-1].data[0]) + * 100) + "%", + xy=(0, median0), xytext=(5, 10), size=15, va="center", + ha="center", bbox=dict(boxstyle="round4", + fc="g", alpha=0.5), + arrowprops=dict(arrowstyle="-|>", + connectionstyle="arc3,rad=-0.2", + fc="w"), + ) plt.grid() - plt.title("performance under covariance shift") - plt.xlabel("method") + plt.legend(handles=legend_handles) + + # plt.title("dependency on noise") + plt.xlabel("noise added [sigma %]") plt.ylabel("absolute error [%]") - # plt.ylim((0, 40)) - plt.savefig(self.output().path + "_temp.png", dpi=500, + plt.ylim((0, 20)) + plt.xlim((0, 15)) + plt.savefig(self.output().path + "temp.png", dpi=500, bbox_inches='tight') + + if __name__ == '__main__': # root folder there the data lies luigi.interface.setup_interface_logging() sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) w.add(WrongNoisePlots(which="colon_muscle")) - w.add(WrongNoisePlots(which="generic")) + # w.add(WrongNoisePlots(which="generic")) w.add(TrainingSamplePlots(which="colon_muscle")) w.add(TrainingSamplePlots(which="generic")) w.add(NoisePlots(which="colon_muscle")) - w.add(NoisePlots(which="generic")) + w.add(DomainNoisePlot()) + w.add(NoisePlotsBVF(which="colon_muscle")) + # w.add(NoisePlots(which="generic")) + w.add(NewDAPlot()) 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 9c791b6cef..8fcc647551 100644 --- a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo.py +++ b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo.py @@ -1,303 +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/" + \ - "2014_12_08_Neil_Pig_Oxygenation_Experiments/organized_recordings" -sp.DATA_FOLDER = "colon_images_selection" -sp.FINALS_FOLDER = "Oxy_selection" + "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(470, 680, 10) * 10 ** -9 - -new_order = [1, 5, 3, 0, 6, 7, 2, 4] +sp.RECORDED_WAVELENGTHS = np.arange(400, 730, 10) * 10 ** -9 def resort_image_wavelengths(collapsed_image): - return collapsed_image[:, new_order] + 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([7]) # [0, 1, 2, 20, 19, 18, 17, 16]) +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): - image_prefix = luigi.Parameter() + folder_name = luigi.Parameter() batch_prefix = luigi.Parameter() def requires(self): - return IPCAIEstimateTissueParametersTask(image_prefix=self.image_prefix, + return IPCAIEstimateTissueParametersTask(folder_name=self.folder_name, batch_prefix= self.batch_prefix), \ - IPCAICorrectImagingSetupTask(image_prefix=self.image_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.image_prefix + "_" + + self.folder_name + "_" + self.batch_prefix + "_summary.png")) def run(self): nrrdReader = NrrdReader() # read in the parametric image param_image = nrrdReader.read(self.input()[0].path) image = param_image.get_image() + # read the multispectral image msi = nrrdReader.read(self.input()[1].path) + # 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()[:, :, [2, 3, 1]] + 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") - # plot parametric maps + # 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): - image_prefix = luigi.Parameter() + folder_name = luigi.Parameter() batch_prefix = luigi.Parameter() def requires(self): - return IPCAIPreprocessMSI(image_prefix=self.image_prefix), \ + return IPCAIPreprocessMSI(folder_name=self.folder_name), \ IPCAITrainRegressor(batch_prefix= self.batch_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, "processed", sp.FINALS_FOLDER, - self.image_prefix + "_" + + self.folder_name + "_" + self.batch_prefix + "estimate.nrrd")) def run(self): # load data reader = NrrdReader() msi = reader.read(self.input()[0].path) e_file = open(self.input()[1].path, 'r') e = pickle.load(e_file) # estimate 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): batch_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "reg_" + self.batch_prefix)) def requires(self): return tasks_mc.CameraBatch(self.batch_prefix) def run(self): # extract data from the batch f = open(self.input().path, "r") batch_train = pickle.load(f) X, y = preprocess(batch_train, w_percent=0.1, bands_to_sortout=sp.bands_to_sortout) # train regressor - reg = RandomForestRegressor(max_depth=30, n_estimators=10, - min_samples_leaf=5, n_jobs=-1) - reg = LinearSaO2Unmixing() + 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): - image_prefix = luigi.Parameter() + folder_name = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, - self.image_prefix + "_preprocessed.nrrd")) + self.folder_name + "_preprocessed.nrrd")) def requires(self): - return IPCAICorrectImagingSetupTask(image_prefix=self.image_prefix) + 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""" - image_prefix = luigi.Parameter() + folder_name = luigi.Parameter() def requires(self): - return MultiSpectralImageFromPNGFiles(image_prefix=self.image_prefix), \ - FlatfieldFromPNGFiles() + return MultiSpectralImageFromPNGFiles(folder_name=self.folder_name), \ + FlatfieldFromMeasurement() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, - self.image_prefix + + 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 FlatfieldFromPNGFiles(luigi.Task): +class FlatfieldFromMeasurement(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_")) + 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) - 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() + folder_name = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, - self.image_prefix + + self.folder_name + ".nrrd")) def run(self): reader = PngReader() msi = reader.read(os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, - self.image_prefix)) + self.folder_name, + "image sample")) writer = NrrdWriter(msi) writer.write(self.output().path) if __name__ == '__main__': - # root folder there the data lies + # 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) - onlyfiles = [ f for f in os.listdir(image_file_folder) if - os.path.isfile(os.path.join(image_file_folder, f)) ] + all_folders = os.listdir(image_file_folder) # each complete stack has elements F0 - F7 - only_colon_images = [ f for f in onlyfiles if - f.endswith("_F7.png") ] - files_prefixes = [ f[:-5] for f in only_colon_images] + only_colon_folders = [ f for f in all_folders if + f.startswith("SaO2") ] - for f in files_prefixes: + for f in only_colon_folders: main_task = IPCAICreateOxyImageTask( - image_prefix=f, + folder_name=f, batch_prefix= - "colon_muscle_tissue_d_ranges_test") + "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 e61bf1999f..50b926ce63 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,317 +1,355 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import Image import ImageEnhance import pickle import time import numpy as np import luigi import matplotlib.pyplot as plt import SimpleITK as sitk from sklearn.ensemble.forest import RandomForestRegressor from scipy.ndimage.filters import gaussian_filter import tasks_analyze as at import tasks_preprocessing as pre import tasks_mc import scriptpaths as sp from msi.io.pngreader import PngReader from msi.io.nrrdreader import NrrdReader from msi.io.nrrdwriter import NrrdWriter import msi.msimanipulations as msimani import msi.imgmani as imgmani import msi.normalize as norm from regression.estimation import estimate_image from regression.linear import LinearSaO2Unmixing -from regression.preprocessing import preprocess, normalize +from regression.preprocessing import preprocess2, normalize from msi.io.tiffreader import TiffReader sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2015_11_12_IPCAI_in_vivo" sp.DATA_FOLDER = "liver_images" sp.FINALS_FOLDER = "liver_oxy" sp.FLAT_FOLDER = "flatfields_liver" +sp.MC_DATA_FOLDER = "mc_data2" # sp.RECORDED_WAVELENGTHS = np.arange(470, 680, 10) * 10 ** -9 new_order = [0, 1, 2, 3, 4, 5, 6, 7] def resort_image_wavelengths(collapsed_image): return collapsed_image[:, new_order] def resort_wavelengths(msi): """Neil organized his wavelengths differently. 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, 20, 19, 18, 17, 16]) class IPCAICreateOxyImageTask(luigi.Task): image_prefix = luigi.Parameter() batch_prefix = luigi.Parameter() def requires(self): return IPCAIEstimateTissueParametersTask(image_prefix=self.image_prefix, batch_prefix= self.batch_prefix), \ + MultiSpectralImageFromTiffFiles(image_prefix=self.image_prefix), \ IPCAICorrectImagingSetupTask(image_prefix=self.image_prefix) + def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, self.image_prefix + "_" + self.batch_prefix + "_summary.png")) def run(self): nrrdReader = NrrdReader() # read in the parametric image param_image = nrrdReader.read(self.input()[0].path) image = param_image.get_image() # read the multispectral image msi = nrrdReader.read(self.input()[1].path) + msi_original = nrrdReader.read(self.input()[2].path) + + plt.figure() + + + # filter dark spots + segmentation = np.max(msi.get_image(), axis=-1) < 100. + # filter bright spots + segmentation = np.logical_or(segmentation, + (np.max(msi.get_image(), axis=-1) > 2000.)) - f, axarr = plt.subplots(1, 2) # create artificial rgb - rgb_image = msi.get_image()[:, :, [2, 3, 1]] + msimani.apply_segmentation(msi_original, ~segmentation) + rgb_image = msi_original.get_image()[:, :, [2, 3, 1]] rgb_image /= np.max(rgb_image) rgb_image *= 255. rgb_image = rgb_image.astype(np.uint8) im = Image.fromarray(rgb_image, 'RGB') enh_brightness = ImageEnhance.Brightness(im) - im = enh_brightness.enhance(2.) + im = enh_brightness.enhance(1.) plot_image = np.array(im) - top_left_axis = axarr[0] + top_left_axis = plt.gca() top_left_axis.imshow(plot_image, interpolation='nearest') - top_left_axis.set_title("synthetic rgb", + top_left_axis.set_title("synthetic rgb (specular pixels black)", fontsize=5) top_left_axis.xaxis.set_visible(False) top_left_axis.yaxis.set_visible(False) + plt.savefig(self.output().path, dpi=1000, bbox_inches='tight') plt.set_cmap("jet") + plt.figure() # plot parametric maps - 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. - at.plot_axis(axarr[1], image[:, :], - "oxygenation [%]") - plt.savefig(self.output().path, dpi=1000, bbox_inches='tight') + # oxy_image = image[:, :, 1] # image[500:1000, 1500:2500, 1] + oxy_image = np.ma.masked_array(image[:, :, 1], segmentation) + oxy_image[np.isnan(oxy_image)] = 0. + oxy_image[np.isinf(oxy_image)] = 0. + oxy_image[oxy_image > 1.] = 1. + oxy_image[oxy_image < 0.] = 0. + oxy_mean = np.mean(oxy_image) + oxy_image[0, 0] = 0. + oxy_image[0, 1] = 1. + at.plot_axis(plt.gca(), oxy_image[:, :] * 100, + "oxygenation [%]. Mean " + + "{0:.1f}".format((oxy_mean * 100.)) + "%") + plt.savefig(self.output().path + "_oxy.png", dpi=1000, bbox_inches='tight') + + bvf_image = np.ma.masked_array(image[:, :, 0], segmentation) + bvf_image[np.isnan(bvf_image)] = 0. + bvf_image[np.isinf(bvf_image)] = 0. + bvf_image[bvf_image > 1.] = 1. + bvf_image[bvf_image < 0.] = 0. + bvf_mean = np.mean(bvf_image) + bvf_image[0, 0] = 0. + bvf_image[1, 1] = 0.1 + at.plot_axis(plt.gca(), bvf_image[:, :] * 100, + "blood volume fraction [%]. Mean " + "{0:.1f}".format((bvf_mean * 100.)) + "%") + plt.savefig(self.output().path + "_bvf.png", dpi=1000, bbox_inches='tight') + + fd = open(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, + sp.FINALS_FOLDER, + 'results_summary.csv'), 'a') + fd.write(self.image_prefix + ", " + str(int(oxy_mean * 100)) + + "," + "{0:.1f}".format((bvf_mean * 100.)) + "\n") + fd.close() class IPCAIEstimateTissueParametersTask(luigi.Task): image_prefix = luigi.Parameter() batch_prefix = luigi.Parameter() def requires(self): return IPCAIPreprocessMSI(image_prefix=self.image_prefix), \ IPCAITrainRegressor(batch_prefix= self.batch_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, "processed", sp.FINALS_FOLDER, self.image_prefix + "_" + self.batch_prefix + "estimate.nrrd")) def run(self): # load data reader = NrrdReader() msi = reader.read(self.input()[0].path) e_file = open(self.input()[1].path, 'r') e = pickle.load(e_file) # estimate start = time.time() sitk_image = estimate_image(msi, e) end = time.time() print "time necessary for estimating image parameters: ", \ str(end - start), "s" # save data outFile = self.output().open('w') outFile.close() sitk.WriteImage(sitk_image, self.output().path) class IPCAITrainRegressor(luigi.Task): batch_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "reg_" + self.batch_prefix)) def requires(self): return tasks_mc.CameraBatch(self.batch_prefix) def run(self): # extract data from the batch f = open(self.input().path, "r") batch_train = pickle.load(f) # TODO: fix hack by sorting always to ascending wavelengths batch_train.reflectances = resort_image_wavelengths( batch_train.reflectances) batch_train.wavelengths = batch_train.wavelengths[new_order] - X, y = preprocess(batch_train, - w_percent=0.05, bands_to_sortout=sp.bands_to_sortout) + X, y = preprocess2(batch_train, + w_percent=0.1, bands_to_sortout=sp.bands_to_sortout) # train regressor reg = RandomForestRegressor(max_depth=10, n_estimators=10, n_jobs=-1) # reg = LinearSaO2Unmixing() reg.fit(X, y) # reg = LinearSaO2Unmixing() # save regressor f = self.output().open('w') pickle.dump(reg, f) f.close() class IPCAIPreprocessMSI(luigi.Task): image_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.image_prefix + "_preprocessed.nrrd")) def requires(self): return IPCAICorrectImagingSetupTask(image_prefix=self.image_prefix) def run(self): reader = NrrdReader() msi = reader.read(self.input().path) # sortout unwanted bands msi.set_image(imgmani.sortout_bands(msi.get_image(), sp.bands_to_sortout)) # zero values would lead to infinity logarithm, thus clip. msi.set_image(np.clip(msi.get_image(), 0.00001, 2. ** 64)) # normalize to get rid of lighting intensity norm.standard_normalizer.normalize(msi) # transform to absorption msi.set_image(-np.log(msi.get_image())) # normalize by l2 for stability norm.standard_normalizer.normalize(msi, "l2") pre.touch_and_save_msi(msi, self.output()) class IPCAICorrectImagingSetupTask(luigi.Task): """corrects for dark image and flatfield""" image_prefix = luigi.Parameter() def requires(self): return MultiSpectralImageFromTiffFiles(image_prefix=self.image_prefix), \ FlatfieldFromPNGFiles() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.image_prefix + "_image_setup_corrected.nrrd")) def run(self): """sort wavelengths and normalize by respective integration times""" # read inputs reader = NrrdReader() msi = reader.read(self.input()[0].path) flatfield = reader.read(self.input()[1].path) # correct image by flatfield msimani.flatfield_correction(msi, flatfield) pre.touch_and_save_msi(msi, self.output()) class FlatfieldFromPNGFiles(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "flatfield" + ".nrrd")) def run(self): reader = PngReader() msi = reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "Flatfield9cm_I000x_F")) writer = NrrdWriter(msi) msi.set_image(msi.get_image().astype(float)) seg_reader = NrrdReader() segmentation = seg_reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "segmentation.nrrd")) # msimani.apply_segmentation(msi, segmentation) # msimani.calculate_mean_spectrum(msi) # apply gaussian smoothing for error reduction img = gaussian_filter(msi.get_image(), sigma=(5, 5, 0), order=0) msi.set_image(img) writer.write(self.output().path) class MultiSpectralImageFromTiffFiles(luigi.Task): image_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.image_prefix + ".nrrd")) def run(self): reader = TiffReader() msi = reader.read(os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, self.image_prefix)) writer = NrrdWriter(msi) writer.write(self.output().path) if __name__ == '__main__': # root folder there the data lies luigi.interface.setup_interface_logging() sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) # run over all subfolders (non-recursively) # to collect the data and generate the results image_file_folder = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER) onlyfiles = [ f for f in os.listdir(image_file_folder) if os.path.isfile(os.path.join(image_file_folder, f)) ] # each complete stack has elements F0 - F7 only_colon_images = [ f for f in onlyfiles if f.startswith("1_") ] files_prefixes = [ f[:-28] for f in only_colon_images] for f in files_prefixes: main_task = IPCAICreateOxyImageTask( image_prefix=f, batch_prefix= - "colon_muscle_tissue_d_ranges_test") + "ipcai_colon_muscle_train") w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_small_bowel.py b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_vivo_small_bowel.py index 91aecd75be..4524bbd118 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,319 +1,292 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import Image import ImageEnhance import pickle -import time +import 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 sklearn.neighbors import KNeighborsRegressor import sklearn.svm from scipy.ndimage.filters import gaussian_filter import tasks_analyze as at -import tasks_preprocessing as pre import tasks_mc import scriptpaths as sp from msi.io.pngreader import PngReader from msi.io.nrrdreader import NrrdReader from msi.io.nrrdwriter import NrrdWriter import msi.msimanipulations as msimani import msi.imgmani as imgmani import msi.normalize as norm from regression.estimation import estimate_image from regression.linear import LinearSaO2Unmixing -from regression.preprocessing import preprocess, normalize -from msi.io.tiffreader import TiffReader - +from regression.preprocessing import preprocess, preprocess2 +from msi.io.tiffringreader import TiffRingReader +from sklearn.linear_model.base import LinearRegression sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2015_11_12_IPCAI_in_vivo" -sp.DATA_FOLDER = "small_bowel_selection" -sp.FINALS_FOLDER = "small_bowel_oxy" +sp.DATA_FOLDER = "small_bowel_selection2" +sp.FINALS_FOLDER = "small_bowel_selection2" sp.FLAT_FOLDER = "flatfields_liver" # sp.RECORDED_WAVELENGTHS = np.arange(470, 680, 10) * 10 ** -9 new_order = [0, 1, 2, 3, 4, 5, 6, 7] def resort_image_wavelengths(collapsed_image): return collapsed_image[:, new_order] def resort_wavelengths(msi): """Neil organized his wavelengths differently. 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, 20, 19, 18, 17, 16]) class IPCAICreateOxyImageTask(luigi.Task): - image_prefix = luigi.Parameter() + image_name = luigi.Parameter() batch_prefix = luigi.Parameter() def requires(self): - return IPCAIEstimateTissueParametersTask(image_prefix=self.image_prefix, - batch_prefix= - self.batch_prefix), \ - IPCAICorrectImagingSetupTask(image_prefix=self.image_prefix) - + return IPCAITrainRegressor(batch_prefix=self.batch_prefix), \ + FlatfieldFromPNGFiles() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, - self.image_prefix + "_" + + self.image_name + "_" + self.batch_prefix + "_summary.png")) def run(self): - nrrdReader = NrrdReader() - # read in the parametric image - param_image = nrrdReader.read(self.input()[0].path) - image = param_image.get_image() - # read the multispectral image - msi = nrrdReader.read(self.input()[1].path) + nrrd_reader = NrrdReader() + tiff_ring_reader = TiffRingReader() + # read the flatfield + flat = nrrd_reader.read(self.input()[1].path) + # read the msi + nr_filters = len(sp.RECORDED_WAVELENGTHS) + full_image_name = os.path.join(sp.ROOT_FOLDER, + sp.DATA_FOLDER, + self.image_name) + msi = tiff_ring_reader.read(full_image_name, nr_filters) + # read the segmentation + seg_path = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, + "segmentation_" + self.image_name + ".nrrd") + segmentation = nrrd_reader.read(seg_path) + segmentation = segmentation.get_image() + # read the regressor + e_file = open(self.input()[0].path, 'r') + e = pickle.load(e_file) - f, axarr = plt.subplots(1, 2) + # correct image setup + filter_nr = int(self.image_name[-6:-5]) + original_order = np.arange(nr_filters) + new_image_order = np.concatenate(( + original_order[nr_filters - filter_nr:], + original_order[:nr_filters - filter_nr])) + # resort msi to restore original order + msimani.get_bands(msi, new_image_order) + # correct by flatfield + msimani.flatfield_correction(msi, flat) # create artificial rgb rgb_image = msi.get_image()[:, :, [2, 3, 1]] rgb_image /= np.max(rgb_image) rgb_image *= 255. + + # preprocess the image + # sortout unwanted bands + print "1" + msi.set_image(imgmani.sortout_bands(msi.get_image(), + sp.bands_to_sortout)) + # zero values would lead to infinity logarithm, thus clip. + msi.set_image(np.clip(msi.get_image(), 0.00001, 2. ** 64)) + # normalize to get rid of lighting intensity + norm.standard_normalizer.normalize(msi) + # transform to absorption + msi.set_image(-np.log(msi.get_image())) + # normalize by l2 for stability + norm.standard_normalizer.normalize(msi, "l2") + print "2" + # estimate + sitk_image = estimate_image(msi, e) + image = sitk.GetArrayFromImage(sitk_image) + + plt.figure() + print "3" rgb_image = rgb_image.astype(np.uint8) im = Image.fromarray(rgb_image, 'RGB') enh_brightness = ImageEnhance.Brightness(im) - im = enh_brightness.enhance(2.) + im = enh_brightness.enhance(5.) plot_image = np.array(im) - top_left_axis = axarr[0] + top_left_axis = plt.gca() 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.savefig(self.output().path, dpi=500, bbox_inches='tight') + plt.figure() plt.set_cmap("jet") - + print "4" # plot 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. - at.plot_axis(axarr[1], image[:, :], - "oxygenation [%]") - plt.savefig(self.output().path, dpi=1000, bbox_inches='tight') - - -class IPCAIEstimateTissueParametersTask(luigi.Task): - image_prefix = luigi.Parameter() - batch_prefix = luigi.Parameter() + oxy_image = image[:, :] + oxy_image = np.ma.masked_array(image[:, :, 1], (segmentation < 1)) + oxy_image[np.isnan(oxy_image)] = 0. + oxy_image[np.isinf(oxy_image)] = 0. + oxy_image[oxy_image > 1.] = 1. + oxy_image[oxy_image < 0.] = 0. + oxy_mean = np.mean(oxy_image) + oxy_image[0, 0] = 0. + oxy_image[0, 1] = 1. + at.plot_axis(plt.gca(), oxy_image[:, :], + "oxygenation [%]. Mean " + + "{0:.1f}".format((oxy_mean * 100.)) + "%") + + plt.savefig(self.output().path + "_oxy.png", dpi=500, bbox_inches='tight') +# plt.figure() +# bvf_image = np.ma.masked_array(image[:, :, 0], (segmentation < 1)) +# # bvf_image = image[:, :, 0] +# bvf_image[np.isnan(bvf_image)] = 0. +# bvf_image[np.isinf(bvf_image)] = 0. +# bvf_image[bvf_image > 1.] = 1. +# bvf_image[bvf_image < 0.] = 0. +# bvf_mean = np.mean(bvf_image) +# bvf_image[0, 0] = 0. +# bvf_image[1, 1] = 0.1 +# at.plot_axis(plt.gca(), bvf_image[:, :] * 100, +# "blood volume fraction [%]. Mean " + "{0:.1f}".format((bvf_mean * 100.)) + "%") +# +# plt.savefig(self.output().path + "_bvf.png", dpi=500, bbox_inches='tight') + print "4.5" + +# fd = open(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, +# sp.FINALS_FOLDER, +# 'results_summary.csv'), 'a') +# fd.write(self.image_name + ", " + str(int(oxy_mean * 100)) + +# ", " + "{0:.1f}".format((bvf_mean * 100.)) + "\n") +# fd.close() + print "5" - def requires(self): - return IPCAIPreprocessMSI(image_prefix=self.image_prefix), \ - IPCAITrainRegressor(batch_prefix= - self.batch_prefix) - def output(self): - return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, "processed", - sp.FINALS_FOLDER, - self.image_prefix + "_" + - self.batch_prefix + - "estimate.nrrd")) - - def run(self): - # load data - reader = NrrdReader() - msi = reader.read(self.input()[0].path) - e_file = open(self.input()[1].path, 'r') - e = pickle.load(e_file) - # estimate - start = time.time() - sitk_image = estimate_image(msi, e) - end = time.time() - print "time necessary for estimating image parameters: ", \ - str(end - start), "s" - # save data - outFile = self.output().open('w') - outFile.close() - sitk.WriteImage(sitk_image, self.output().path) class IPCAITrainRegressor(luigi.Task): batch_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "reg_" + self.batch_prefix)) def requires(self): return tasks_mc.CameraBatch(self.batch_prefix) def run(self): # extract data from the batch f = open(self.input().path, "r") batch_train = pickle.load(f) # TODO: fix hack by sorting always to ascending wavelengths batch_train.reflectances = resort_image_wavelengths( batch_train.reflectances) batch_train.wavelengths = batch_train.wavelengths[new_order] - X, y = preprocess(batch_train, - w_percent=0.10, bands_to_sortout=sp.bands_to_sortout) + X, y = preprocess2(batch_train, + w_percent=0.1, bands_to_sortout=sp.bands_to_sortout) # train regressor reg = RandomForestRegressor(max_depth=10, n_estimators=10, n_jobs=-1) + # reg = KNeighborsRegressor(algorithm="auto") + # reg = LinearRegression() # reg = sklearn.svm.SVR(kernel="rbf", degree=3, C=100., gamma=10.) # reg = LinearSaO2Unmixing() - reg.fit(X, y) + reg.fit(X, y[:, 1]) # reg = LinearSaO2Unmixing() # save regressor f = self.output().open('w') pickle.dump(reg, f) f.close() -class IPCAIPreprocessMSI(luigi.Task): - image_prefix = luigi.Parameter() - - def output(self): - return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, - sp.RESULTS_FOLDER, - self.image_prefix + "_preprocessed.nrrd")) - - def requires(self): - return IPCAICorrectImagingSetupTask(image_prefix=self.image_prefix) - - def run(self): - reader = NrrdReader() - msi = reader.read(self.input().path) - # sortout unwanted bands - msi.set_image(imgmani.sortout_bands(msi.get_image(), - sp.bands_to_sortout)) - # zero values would lead to infinity logarithm, thus clip. - msi.set_image(np.clip(msi.get_image(), 0.00001, 2. ** 64)) - # normalize to get rid of lighting intensity - norm.standard_normalizer.normalize(msi) - # transform to absorption - msi.set_image(-np.log(msi.get_image())) - # normalize by l2 for stability - norm.standard_normalizer.normalize(msi, "l2") - pre.touch_and_save_msi(msi, self.output()) - - -class IPCAICorrectImagingSetupTask(luigi.Task): - """corrects for dark image - and flatfield""" - image_prefix = luigi.Parameter() - - def requires(self): - return MultiSpectralImageFromTiffFiles(image_prefix=self.image_prefix), \ - FlatfieldFromPNGFiles() - - def output(self): - return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, - sp.RESULTS_FOLDER, - self.image_prefix + - "_image_setup_corrected.nrrd")) - - def run(self): - """sort wavelengths and normalize by respective integration times""" - # read inputs - reader = NrrdReader() - msi = reader.read(self.input()[0].path) - flatfield = reader.read(self.input()[1].path) - # correct image by flatfield - msimani.flatfield_correction(msi, flatfield) - pre.touch_and_save_msi(msi, self.output()) - - class FlatfieldFromPNGFiles(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "flatfield" + ".nrrd")) def run(self): reader = PngReader() msi = reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "Flatfield9cm_I000x_F")) writer = NrrdWriter(msi) msi.set_image(msi.get_image().astype(float)) seg_reader = NrrdReader() segmentation = seg_reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "segmentation.nrrd")) # msimani.apply_segmentation(msi, segmentation) # msimani.calculate_mean_spectrum(msi) # apply gaussian smoothing for error reduction img = gaussian_filter(msi.get_image(), sigma=(5, 5, 0), order=0) msi.set_image(img) writer.write(self.output().path) -class MultiSpectralImageFromTiffFiles(luigi.Task): - image_prefix = luigi.Parameter() - - def output(self): - return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, - sp.RESULTS_FOLDER, - self.image_prefix + - ".nrrd")) - - def run(self): - reader = TiffReader() - msi = reader.read(os.path.join(sp.ROOT_FOLDER, - sp.DATA_FOLDER, - self.image_prefix)) - writer = NrrdWriter(msi) - writer.write(self.output().path) - if __name__ == '__main__': # root folder there the data lies + logging.basicConfig(filename='small_bowel' + str(datetime.datetime.now()) + + '.log', level=logging.INFO) luigi.interface.setup_interface_logging() + ch = logging.StreamHandler() + ch.setLevel(logging.INFO) + logger = logging.getLogger() + logger.addHandler(ch) + sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) # run over all subfolders (non-recursively) # to collect the data and generate the results image_file_folder = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER) onlyfiles = [ f for f in os.listdir(image_file_folder) if os.path.isfile(os.path.join(image_file_folder, f)) ] - # each complete stack has elements F0 - F7 - only_colon_images = [ f for f in onlyfiles if - f.endswith("F0.tiff") ] - files_prefixes = [ f[:-6] for f in only_colon_images] + onlyfiles.sort() + only_first_files = [ f for f in onlyfiles if + f.endswith("F0.tiff") ] - for f in files_prefixes: + for f in only_first_files: main_task = IPCAICreateOxyImageTask( - image_prefix=f, + image_name=f, batch_prefix= - "colon_muscle_tissue_d_ranges_train") + "ipcai_colon_muscle_train") w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_calculate_spectra.py b/Modules/Biophotonics/python/iMC/script_calculate_spectra.py index 824905404b..50716c6f63 100644 --- a/Modules/Biophotonics/python/iMC/script_calculate_spectra.py +++ b/Modules/Biophotonics/python/iMC/script_calculate_spectra.py @@ -1,49 +1,40 @@ ''' Created on Sep 9, 2015 @author: wirkert ''' import logging import datetime import numpy as np import luigi import scriptpaths as sp import tasks_mc import mc.factories as mcfac # general output path config sp.ROOT_FOLDER = \ "/media/wirkert/data/Data/2015_11_12_IPCAI_in_silico" sp.RESULTS_FOLDER = "mc_data" if __name__ == '__main__': logging.basicConfig(filename='in_silico' + 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, 10, 1) + BATCH_NUMBERS = np.arange(0, 100, 1) for i in BATCH_NUMBERS: - main_task = tasks_mc.CreateSpectraTask("ipcai_generic", i, 1000, - mcfac.GenericMcFactory()) - w.add(main_task) - w.run() - main_task = tasks_mc.CreateSpectraTask("ipcai_colon_muscle", i, 1000, - mcfac.ColonMuscleMcFactory()) - w.add(main_task) - w.run() - main_task = tasks_mc.CreateSpectraTask("ipcai_generic_", i + - np.max(BATCH_NUMBERS), 1000, - mcfac.GenericMcFactory()) + main_task = tasks_mc.CreateSpectraTask("ipcai_less_generic", i, 1000, + mcfac.LessGenericMcFactory()) w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_plot_one_spectrum.py b/Modules/Biophotonics/python/iMC/script_plot_one_spectrum.py index f9b44a8844..0558b261fb 100644 --- a/Modules/Biophotonics/python/iMC/script_plot_one_spectrum.py +++ b/Modules/Biophotonics/python/iMC/script_plot_one_spectrum.py @@ -1,70 +1,70 @@ ''' Created on Oct 12, 2015 @author: wirkert ''' import pickle import logging import numpy as np import matplotlib.pyplot as plt import luigi import tasks_regression as rt from msi.plot import plot from msi.msi import Msi import msi.normalize as norm import scriptpaths as sp sp.ROOT_FOLDER = "/media/wirkert/data/Data/2015_xxxx_plot_one_spectrum" # the wavelengths recorded by our camera RECORDED_WAVELENGTHS = \ - np.array([580, 470, 660, 560, 480, 511, 600, 700]) * 10 ** -9 -PARAMS = np.array([0.05, # bvf - 0.0, # SaO2 + np.arange(450, 720, 2) * 10 ** -9 +PARAMS = np.array([0.1, # bvf + 0.7, # SaO2 0.0, # billirubin 500., # a_mie 0.0, # a_ray 1.091, # b (for scattering 500. * 10 ** -6]) # d_muc class PlotOneSpectrum(luigi.Task): batch_prefix = luigi.Parameter() def requires(self): return rt.TrainForestForwardModel(self.batch_prefix) def run(self): f = file(self.input().path, "r") rf = pickle.load(f) f.close() refl = rf.predict(PARAMS) msi = Msi(refl) msi.set_wavelengths(RECORDED_WAVELENGTHS) - norm.standard_normalizer.normalize(msi) + # norm.standard_normalizer.normalize(msi) plot(msi) plt.gca().set_xlabel("wavelength") plt.gca().set_ylabel("normalized reflectance") plt.grid() plt.ylim([0.0, 0.4]) plt.title("bvf: " + str(PARAMS[0]) + "; saO2: " + str(PARAMS[1]) + "; bili: " + str(PARAMS[2]) + "; a_mie: " + str(PARAMS[3]) + "; a_ray: " + str(PARAMS[4]) + "; d_muc: " + str(PARAMS[6])) plt.show() if __name__ == '__main__': logging.basicConfig(level=logging.INFO) luigi.interface.setup_interface_logging() sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) main_task = PlotOneSpectrum(batch_prefix= "jacques_no_billi_generic_scattering_") w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/script_visualize_spectra.py b/Modules/Biophotonics/python/iMC/script_visualize_spectra.py index 1dce34c368..20a8ac8e78 100644 --- a/Modules/Biophotonics/python/iMC/script_visualize_spectra.py +++ b/Modules/Biophotonics/python/iMC/script_visualize_spectra.py @@ -1,89 +1,92 @@ ''' Created on Sep 10, 2015 @author: wirkert ''' import logging import os import pickle import numpy as np import matplotlib.pyplot as plt +import matplotlib import luigi +from regression.preprocessing import preprocess2, normalize import scriptpaths as sp import tasks_mc import mc.plot import mc.factories as mcfac # general output path config sp.ROOT_FOLDER = \ "/media/wirkert/data/Data/2015_06_01_Filtertransmittance_Spectrometer" sp.DATA_FOLDER = "spectrometer_measurements" sp.MC_DATA_FOLDER = "processed" SPECTRAL_PLOTS = "spectral_plots" +sp.RECORDED_WAVELENGTHS = np.arange(470, 690, 10) * 10 ** -9 +# matplotlib.rcParams.update({'font.size': 18}) class VisualizeSpectraTask(luigi.Task): batch_prefix = luigi.Parameter() batch_nr = luigi.IntParameter() nr_samples = luigi.IntParameter() def requires(self): return tasks_mc.CreateSpectraTask(self.batch_prefix, self.batch_nr, self.nr_samples, mcfac.VisualizationMcFactory()), \ tasks_mc.CameraBatch(self.batch_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, SPECTRAL_PLOTS, "plotted_spectra_" + self.batch_prefix + "_" + str(self.batch_nr) + ".png")) def run(self): f = file(self.input()[0].path, "r") batch = pickle.load(f) f = file(self.input()[1].path, "r") camera_batch = pickle.load(f) f.close() - f, axarr = plt.subplots(2, 1) + 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 - org_ax = axarr[0] + # axis for plotting the folded spectrum + normalization - fold_ax = axarr[1] - mc.plot.plot(batch, org_ax) - mc.plot.plot(camera_batch, fold_ax) + + mc.plot.plot(batch) # tidy up and save plot - major_ticks = np.arange(450, 720, 50) * 10 ** -9 - minor_ticks = np.arange(450, 720, 10) * 10 ** -9 - org_ax.set_title(self.batch_prefix) - org_ax.set_xticks(major_ticks) - org_ax.set_xticks(minor_ticks, minor=True) - fold_ax.set_title("response to camera filter sensitivities") - fold_ax.set_xticks(major_ticks) - fold_ax.set_xticks(minor_ticks, minor=True) - org_ax.grid() - fold_ax.grid() - plt.savefig(self.output().path, dpi=500) + 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("vizualization_generic", - 0, 10) + main_task = VisualizeSpectraTask("paper2", + 0, 1) w.add(main_task) w.run() diff --git a/Modules/Biophotonics/python/iMC/tasks_mc.py b/Modules/Biophotonics/python/iMC/tasks_mc.py index d162d99449..40948f5631 100644 --- a/Modules/Biophotonics/python/iMC/tasks_mc.py +++ b/Modules/Biophotonics/python/iMC/tasks_mc.py @@ -1,218 +1,218 @@ ''' Created on Sep 10, 2015 @author: wirkert ''' ''' Created on Aug 31, 2015 @author: wirkert ''' import os import numpy as np import pickle import time import logging import luigi import scriptpaths as sp from mc.sim import SimWrapper, get_diffuse_reflectance from mc.batches import join_batches import mc.mcmanipulations as mcmani from msi.io.spectrometerreader import SpectrometerReader from msi.io.msiwriter import MsiWriter from msi.io.msireader import MsiReader from msi.normalize import NormalizeMean import msi.msimanipulations as msimani # experiment configuration MCI_FILENAME = "./temp.mci" MCO_FILENAME = "temp.mco" # this path definitly needs to be adapted by you PATH_TO_MCML = "/home/wirkert/workspace/monteCarlo/gpumcml/fast-gpumcml/" EXEC_MCML = "gpumcml.sm_20" class SpectrometerFile(luigi.Task): input_file = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.SPECTROMETER_FOLDER, self.input_file)) class FilterTransmission(luigi.Task): input_file = luigi.Parameter() def requires(self): return SpectrometerFile(self.input_file) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "processed_transmission" + self.input_file + ".msi")) def run(self): reader = SpectrometerReader() filter_transmission = reader.read(self.input().path) # filter high and low wavelengths wavelengths = filter_transmission.get_wavelengths() fi_image = filter_transmission.get_image() fi_image[wavelengths < 450 * 10 ** -9] = 0.0 fi_image[wavelengths > 720 * 10 ** -9] = 0.0 # filter elements farther away than +- 30nm name_to_float = float(os.path.splitext(self.input_file)[0]) fi_image[wavelengths < (name_to_float - 30) * 10 ** -9] = 0.0 fi_image[wavelengths > (name_to_float + 30) * 10 ** -9] = 0.0 # elements < 0 are set to 0. fi_image[fi_image < 0.0] = 0.0 # write it writer = MsiWriter(filter_transmission) writer.write(self.output().path) class JoinBatches(luigi.Task): batch_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.batch_prefix + "_" + "all" + ".imc")) def run(self): path = os.path.join(sp.ROOT_FOLDER, sp.MC_DATA_FOLDER) # get all files in the search path files = [ f for f in os.listdir(path) \ if os.path.isfile(os.path.join(path, f)) ] # from these get only those who start with correct batch prefix batch_file_names = \ [ f for f in files if f.startswith(self.batch_prefix)] batch_files = \ [open(os.path.join(path, f), 'r') for f in batch_file_names] # load these files batches = [pickle.load(f) for f in batch_files] # now join them to one batch joined_batch = reduce(join_batches, batches) # write it joined_batch_file = open(self.output().path, 'w') # there seems to be a bug in pickle. thus, unfortunately the generator # has to be removed before saving joined_batch.generator = None pickle.dump(joined_batch, joined_batch_file) class CameraBatch(luigi.Task): """takes a batch of reference data and converts it to the spectra processed by the camera""" batch_prefix = luigi.Parameter() def requires(self): return JoinBatches(self.batch_prefix) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.batch_prefix + "_all_camera.imc")) def run(self): f = file(self.input().path, "r") batch = pickle.load(f) f.close() # camera batch creation: camera_batch = batch - mcmani.fold_by_sliding_average(camera_batch, 5) + mcmani.fold_by_sliding_average(camera_batch, 6) mcmani.interpolate_wavelengths(camera_batch, sp.RECORDED_WAVELENGTHS) # write it joined_batch_file = open(self.output().path, 'w') pickle.dump(camera_batch, joined_batch_file) class CreateSpectraTask(luigi.Task): batch_prefix = luigi.Parameter() batch_nr = luigi.IntParameter() nr_samples = luigi.IntParameter() factory = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, self.batch_prefix + "_" + str(self.batch_nr) + ".imc")) def run(self): start = time.time() # setup simulation wrapper self.sim_wrapper = SimWrapper() self.sim_wrapper.set_mci_filename(MCI_FILENAME) self.sim_wrapper.set_mcml_executable(PATH_TO_MCML + EXEC_MCML) # setup model self.tissue_model = self.factory.create_tissue_model() self.tissue_model.set_mci_filename(self.sim_wrapper.mci_filename) self.tissue_model.set_mco_filename(MCO_FILENAME) # setup array in which data shall be stored batch = self.factory.create_batch_to_simulate() batch.create_parameters(self.nr_samples) for i in range(batch.nr_elements()): self.tissue_model.set_batch_element(batch, i) logging.info("running simulation " + str(i) + " for\n" + str(self.tissue_model)) start = time.time() # map the wavelengths array to reflectance list reflectances = map(self.wavelength_to_reflectance, batch.wavelengths) end = time.time() # store in batch_nr batch.reflectances[i, :] = np.asarray(reflectances) # success! logging.info("successfully ran simulation in " + "{:.2f}".format(end - start) + " seconds") # convert the lists in batch to np arrays batch.reflectances = np.array(batch.reflectances) # clean up temporarily created files os.remove(MCI_FILENAME) created_mco_file = PATH_TO_MCML + "/" + MCO_FILENAME if os.path.isfile(created_mco_file): os.remove(created_mco_file) # save the created output f = open(self.output().path, 'w') # there seems to be a bug in pickle. thus, unfortunately the generator # has to be removed before saving batch.generator = None pickle.dump(batch, f) f.close() end = time.time() logging.info("time for creating batch of mc data: %.f s" % (end - start)) def wavelength_to_reflectance(self, wavelength): """helper function to determine the reflectance for a given wavelength using the current model and simulation """ self.tissue_model.set_wavelength(wavelength) self.tissue_model.create_mci_file() if os.path.isfile(PATH_TO_MCML + EXEC_MCML): self.sim_wrapper.run_simulation() return get_diffuse_reflectance(PATH_TO_MCML + MCO_FILENAME) else: raise IOError("path to gpumcml not valid") def get_transmission_data(input_path, desired_wavelengths): """small helper method to get filter transmission from file at input_path. The wavelengths are interpolated to the desired ones""" reader = MsiReader() filter_transmission = reader.read(input_path) msimani.interpolate_wavelengths(filter_transmission, desired_wavelengths) # normalize to one normalizer = NormalizeMean() normalizer.normalize(filter_transmission) return filter_transmission