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 e0f0a2d9a0..8f4ee0040d 100644 --- a/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_silico.py +++ b/Modules/Biophotonics/python/iMC/script_analyze_ipcai_in_silico.py @@ -1,637 +1,663 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import pickle from collections import namedtuple import luigi import matplotlib.pyplot as plt import matplotlib from sklearn.ensemble.forest import RandomForestRegressor import tasks_mc import scriptpaths as sp 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", rf, "green", 0.5)] # do validation 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.xlabel("number of training samples") plt.ylabel("absolute error [%]") 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) 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=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, 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.legend(handles=legend_handles) # plt.title("dependency on wrong training noise [w_training = 10%]") plt.xlabel("noise added for testing [sigma %]") plt.ylabel("absolute error [%]") 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") + tasks_mc.CameraBatch("ipcai_generic") 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) 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, 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.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 WeightsCalculation(luigi.Task): source_domain = luigi.Parameter() target_domain = luigi.Parameter() def requires(self): 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")) + "weights_" + self.source_domain + "to" + + self.target_domain + ".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=None) X_da_train, y_da_train = preprocess(batch_da_train, 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) 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 NewDAPlot(luigi.Task): + source_train_set = luigi.Parameter() + target_train_set = luigi.Parameter() + target_test_set = luigi.Parameter() + def requires(self): - return tasks_mc.CameraBatch("ipcai_generic"), \ - tasks_mc.CameraBatch("ipcai_colon_muscle_test"), \ - WeightsCalculation("ipcai_generic", "ipcai_colon_muscle_test") + return tasks_mc.CameraBatch(self.source_train_set), \ + tasks_mc.CameraBatch(self.target_test_set), \ + tasks_mc.CameraBatch(self.target_train_set), \ + WeightsCalculation(self.source_train_set, self.target_train_set) def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "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_test = pickle.load(f) f = open(self.input()[2].path, "r") + batch_target_train = pickle.load(f) + f = open(self.input()[3].path, "r") weights = np.load(f) evaluation_setups = [ - EvaluationStruct("Linear Beer-Lambert", - LinearSaO2Unmixing(), [], "red", 0.25) + EvaluationStruct("Proposed no DS", + rf, [], "red", 0.25) , EvaluationStruct("Proposed", rf, [], "yellow", 0.5) , EvaluationStruct("Proposed DA", rf, [], "green", 0.5) + , EvaluationStruct("Proposed Fake", rf, [], "black", 0.1) ] 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) + X_target_train, y_target_train = preprocess(batch_target_train, + 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) 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_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") +# +# 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() # 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) + + X_fake_train, y_fake_train = preprocess2(batch_train, + w_percent=w) + correct_bvf_train = y_fake_train[:, 0] < 0.1 + X_fake_train = X_fake_train[correct_bvf_train, :] + y_fake_train = y_fake_train[correct_bvf_train, 1] + for e in evaluation_setups: lr = e.regressor if "DA" in e.name: lr.fit(X_train, y_train, weights) + elif "DS" in e.name: + lr.fit(X_target_train, y_target_train) + elif "Fake" in e.name: + lr.fit(X_fake_train, y_fake_train) 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 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]) + 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.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 + "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(TrainingSamplePlots(which="colon_muscle")) w.add(TrainingSamplePlots(which="generic")) w.add(NoisePlots(which="colon_muscle")) w.add(DomainNoisePlot()) w.add(NoisePlotsBVF(which="colon_muscle")) # w.add(NoisePlots(which="generic")) - w.add(NewDAPlot()) + w.add(NewDAPlot(source_train_set="ipcai_less_generic", + target_train_set="ipcai_colon_muscle_train", + target_test_set="ipcai_colon_muscle_test")) 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 071dee0395..7309121468 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,294 +1,297 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 14 11:09:18 2015 @author: wirkert """ import os import Image import ImageEnhance import pickle import logging import datetime import numpy as np import luigi import matplotlib.pyplot as plt import SimpleITK as sitk from sklearn.ensemble.forest import RandomForestRegressor -from sklearn.neighbors import KNeighborsRegressor -import sklearn.svm from scipy.ndimage.filters import gaussian_filter import tasks_analyze as at import tasks_mc import scriptpaths as sp from msi.io.pngreader import PngReader from msi.io.nrrdreader import NrrdReader from msi.io.nrrdwriter import NrrdWriter import msi.msimanipulations as msimani import msi.imgmani as imgmani import msi.normalize as norm from regression.estimation import estimate_image from regression.linear import LinearSaO2Unmixing from regression.preprocessing import preprocess, preprocess2 from msi.io.tiffringreader import TiffRingReader from sklearn.linear_model.base import LinearRegression sp.ROOT_FOLDER = "/media/wirkert/data/Data/" + \ "2015_11_12_IPCAI_in_vivo" sp.DATA_FOLDER = "small_bowel_images" -sp.FINALS_FOLDER = "small_bowel_images" +sp.FINALS_FOLDER = "small_bowel_images_colorbar20_80" sp.FLAT_FOLDER = "flatfields_liver" # sp.RECORDED_WAVELENGTHS = np.arange(470, 680, 10) * 10 ** -9 new_order = [0, 1, 2, 3, 4, 5, 6, 7] def resort_image_wavelengths(collapsed_image): return collapsed_image[:, new_order] def resort_wavelengths(msi): """Neil organized his wavelengths differently. 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_name = luigi.Parameter() batch_prefix = luigi.Parameter() def requires(self): return IPCAITrainRegressor(batch_prefix=self.batch_prefix), \ FlatfieldFromPNGFiles() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, self.image_name + "_" + self.batch_prefix + - "_summary.png")) + "_summary.png.png")) def run(self): nrrd_reader = NrrdReader() tiff_ring_reader = TiffRingReader() # read the flatfield flat = nrrd_reader.read(self.input()[1].path) # read the msi nr_filters = len(sp.RECORDED_WAVELENGTHS) full_image_name = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, self.image_name) msi, segmentation = tiff_ring_reader.read(full_image_name, nr_filters) + segmentation = np.logical_and(segmentation, + (np.max(msi.get_image(), + axis=-1) < 1000.)) # read the segmentation # # seg_path = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER, # "segmentation_" + self.image_name + ".nrrd") # segmentation = nrrd_reader.read(seg_path) # segmentation = segmentation.get_image() # read the regressor e_file = open(self.input()[0].path, 'r') e = pickle.load(e_file) # correct image setup filter_nr = int(self.image_name[-6:-5]) original_order = np.arange(nr_filters) new_image_order = np.concatenate(( original_order[nr_filters - filter_nr:], original_order[:nr_filters - filter_nr])) # resort msi to restore original order msimani.get_bands(msi, new_image_order) # correct by flatfield msimani.flatfield_correction(msi, flat) # create artificial rgb rgb_image = msi.get_image()[:, :, [2, 3, 1]] rgb_image /= np.max(rgb_image) rgb_image *= 255. # preprocess the image # sortout unwanted bands print "1" msi.set_image(imgmani.sortout_bands(msi.get_image(), sp.bands_to_sortout)) # zero values would lead to infinity logarithm, thus clip. msi.set_image(np.clip(msi.get_image(), 0.00001, 2. ** 64)) # normalize to get rid of lighting intensity norm.standard_normalizer.normalize(msi) # transform to absorption msi.set_image(-np.log(msi.get_image())) # normalize by l2 for stability norm.standard_normalizer.normalize(msi, "l2") print "2" # estimate sitk_image = estimate_image(msi, e) image = sitk.GetArrayFromImage(sitk_image) plt.figure() print "3" rgb_image = rgb_image.astype(np.uint8) im = Image.fromarray(rgb_image, 'RGB') enh_brightness = ImageEnhance.Brightness(im) im = enh_brightness.enhance(5.) plot_image = np.array(im) top_left_axis = plt.gca() top_left_axis.imshow(plot_image, interpolation='nearest') top_left_axis.xaxis.set_visible(False) top_left_axis.yaxis.set_visible(False) - plt.savefig(self.output().path, dpi=500, bbox_inches='tight') - plt.figure() plt.set_cmap("jet") print "4" # plot parametric maps oxy_image = image[:, :] + segmentation[0, 0] = 1 + segmentation[0, 1] = 1 oxy_image = np.ma.masked_array(image[:, :], (segmentation < 1)) oxy_image[np.isnan(oxy_image)] = 0. oxy_image[np.isinf(oxy_image)] = 0. - oxy_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. + oxy_image[0, 0] = 0.2 + oxy_image[0, 1] = 0.8 + oxy_image = np.clip(oxy_image, 0.2, 0.8) +# oxy_image[oxy_image > 0.8] = 0.8 +# oxy_image[oxy_image < 0.2] = 0.2 + at.plot_axis(plt.gca(), oxy_image[:, :], "oxygenation [%]. Mean " + "{0:.1f}".format((oxy_mean * 100.)) + "%") - plt.savefig(self.output().path + "_oxy.png", dpi=500, bbox_inches='tight') + plt.savefig(self.output().path, dpi=500, bbox_inches='tight') # plt.figure() # bvf_image = np.ma.masked_array(image[:, :, 0], (segmentation < 1)) # # bvf_image = image[:, :, 0] # bvf_image[np.isnan(bvf_image)] = 0. # bvf_image[np.isinf(bvf_image)] = 0. # bvf_image[bvf_image > 1.] = 1. # bvf_image[bvf_image < 0.] = 0. # bvf_mean = np.mean(bvf_image) # bvf_image[0, 0] = 0. # bvf_image[1, 1] = 0.1 # at.plot_axis(plt.gca(), bvf_image[:, :] * 100, # "blood volume fraction [%]. Mean " + "{0:.1f}".format((bvf_mean * 100.)) + "%") # # plt.savefig(self.output().path + "_bvf.png", dpi=500, bbox_inches='tight') print "4.5" # fd = open(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, # sp.FINALS_FOLDER, # 'results_summary.csv'), 'a') # fd.write(self.image_name + ", " + str(int(oxy_mean * 100)) + # ", " + "{0:.1f}".format((bvf_mean * 100.)) + "\n") # fd.close() print "5" plt.close("all") class IPCAITrainRegressor(luigi.Task): batch_prefix = luigi.Parameter() def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, sp.FINALS_FOLDER, "reg_" + self.batch_prefix)) def requires(self): return tasks_mc.CameraBatch(self.batch_prefix) def run(self): # extract data from the batch f = open(self.input().path, "r") batch_train = pickle.load(f) # TODO: fix hack by sorting always to ascending wavelengths batch_train.reflectances = resort_image_wavelengths( batch_train.reflectances) batch_train.wavelengths = batch_train.wavelengths[new_order] X, y = preprocess2(batch_train, w_percent=0.1, bands_to_sortout=sp.bands_to_sortout) # train regressor reg = RandomForestRegressor(10, min_samples_leaf=10, max_depth=9, n_jobs=-1) # reg = KNeighborsRegressor(algorithm="auto") # reg = LinearRegression() # reg = sklearn.svm.SVR(kernel="rbf", degree=3, C=100., gamma=10.) # reg = LinearSaO2Unmixing() reg.fit(X, y[:, 1]) # reg = LinearSaO2Unmixing() # save regressor f = self.output().open('w') pickle.dump(reg, f) f.close() class FlatfieldFromPNGFiles(luigi.Task): def output(self): return luigi.LocalTarget(os.path.join(sp.ROOT_FOLDER, sp.RESULTS_FOLDER, "flatfield" + ".nrrd")) def run(self): reader = PngReader() msi = reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "Flatfield9cm_I000x_F")) writer = NrrdWriter(msi) msi.set_image(msi.get_image().astype(float)) seg_reader = NrrdReader() segmentation = seg_reader.read(os.path.join(sp.ROOT_FOLDER, sp.FLAT_FOLDER, "segmentation.nrrd")) # msimani.apply_segmentation(msi, segmentation) # msimani.calculate_mean_spectrum(msi) # apply gaussian smoothing for error reduction img = gaussian_filter(msi.get_image(), sigma=(5, 5, 0), order=0) msi.set_image(img) writer.write(self.output().path) if __name__ == '__main__': # root folder there the data lies logging.basicConfig(filename='small_bowel' + str(datetime.datetime.now()) + '.log', level=logging.INFO) luigi.interface.setup_interface_logging() ch = logging.StreamHandler() ch.setLevel(logging.INFO) logger = logging.getLogger() logger.addHandler(ch) sch = luigi.scheduler.CentralPlannerScheduler() w = luigi.worker.Worker(scheduler=sch) # run over all subfolders (non-recursively) # to collect the data and generate the results image_file_folder = os.path.join(sp.ROOT_FOLDER, sp.DATA_FOLDER) onlyfiles = [ f for f in os.listdir(image_file_folder) if os.path.isfile(os.path.join(image_file_folder, f)) ] onlyfiles.sort() onlyfiles = [ f for f in onlyfiles if f.endswith(".tiff") ] for f in onlyfiles: main_task = IPCAICreateOxyImageTask( image_name=f, batch_prefix= "ipcai_colon_muscle_train") w.add(main_task) w.run()