diff --git a/Modules/Biophotonics/python/iMC/msi/io/msiwriter.py b/Modules/Biophotonics/python/iMC/msi/io/msiwriter.py index 2541d2ebd8..815c465a54 100644 --- a/Modules/Biophotonics/python/iMC/msi/io/msiwriter.py +++ b/Modules/Biophotonics/python/iMC/msi/io/msiwriter.py @@ -1,26 +1,25 @@ ''' Created on Aug 25, 2015 @author: wirkert ''' import pickle from writer import Writer class MsiWriter(Writer): ''' The MsiReader writing the Msi as a serialized python object. - This is the prefered way of writing an Msi. ''' def __init__(self, msiToWrite): """ initialize the write with a specific multi spectral image (class Msi) """ self.msiToWrite = msiToWrite def write (self, uri): f = open(uri, 'w') pickle.dump(self.msiToWrite, f) f.close() diff --git a/Modules/Biophotonics/python/iMC/msi/io/reader.py b/Modules/Biophotonics/python/iMC/msi/io/reader.py index fa46016d0b..99d6b678c3 100644 --- a/Modules/Biophotonics/python/iMC/msi/io/reader.py +++ b/Modules/Biophotonics/python/iMC/msi/io/reader.py @@ -1,20 +1,19 @@ # -*- coding: utf-8 -*- """ Created on Fri Aug 7 13:39:16 2015 @author: wirkert """ + class Reader: """ - This should be redone. Reader could take care - of everything but reading in one image wavelength. - This should be specialized for tiffs, pngs, ... + Abstract reader base class """ def __init__(self): pass def read(self, fileToRead): return None diff --git a/Modules/Biophotonics/python/iMC/msi/io/writer.py b/Modules/Biophotonics/python/iMC/msi/io/writer.py index f735141b54..9342ff00fe 100644 --- a/Modules/Biophotonics/python/iMC/msi/io/writer.py +++ b/Modules/Biophotonics/python/iMC/msi/io/writer.py @@ -1,14 +1,18 @@ # -*- coding: utf-8 -*- """ Created on Thu Aug 13 09:47:31 2015 @author: wirkert """ + class Writer(): + """ + Abstract writer base class + """ def __init__(self): pass def write(self, fileToWrite): return None \ No newline at end of file diff --git a/Modules/Biophotonics/python/iMC/msi/msimanipulations.py b/Modules/Biophotonics/python/iMC/msi/msimanipulations.py index 19f406c723..9169fa35c1 100644 --- a/Modules/Biophotonics/python/iMC/msi/msimanipulations.py +++ b/Modules/Biophotonics/python/iMC/msi/msimanipulations.py @@ -1,120 +1,142 @@ # -*- coding: utf-8 -*- """ Created on Thu Aug 13 13:42:00 2015 @author: wirkert """ import logging import copy import numpy as np from scipy.interpolate import interp1d from imgmani import collapse_image import imgmani from msi import Msi +''' +The msi manipulations module includes usefull convenience operations on msis. + +E.g. calculate_mean_spectrum to calculate the average spectrum on all the image +data or interpolate_wavelengths to change to a different wavelenght set by +simple interpolation. + +All methods take a msi and change it in place. They also return the same msi +object for convenience (can e.g. be used to chain msi manipulations). +''' + + def apply_segmentation(msi, segmentation): - """ TODO """ + """ applies a segmentation to an msi. + + If the msi has imaging data of n x m x nr_wavelengths the segmentation + has to be a numpy array of n x m size. pixelvalues with values different + than zero will be included in the segmentation. + By applying the segmentation, not segmented elements will be np.ma.masked. + + Alternatively, one can input a msi with the mentioned n x m numpy array + in from of a msi as segmentation (for convenience to be able to use the + same reader for msis and segmentations) + """ if (isinstance(segmentation, Msi)): # expects just an image, but if a Msi is passed it's also ok segmentation = segmentation.get_image() segmentation = np.squeeze(segmentation) mask = (0 == segmentation) # mask needs to be expanded to cover all wavlengths wholeMask = np.zeros_like(msi.get_image(), dtype="bool") # doesn't seem elegant for i in range(wholeMask.shape[-1]): wholeMask[:, :, i] = mask msi.set_mask(wholeMask) return msi def calculate_mean_spectrum(msi): """ reduce this image to only its mean spectrum. If the msi.get_image() is a masked array these values will be ignored when calculating the mean spectrum """ # reshape to collapse all but last dimension (which contains reflectances) collapsedImage = collapse_image(msi.get_image()) msi.set_image(np.mean(collapsedImage, axis=0)) # returns the same msi. return msi def interpolate_wavelengths(msi, newWavelengths): """ interpolate image data to fit newWavelengths. Current implementation performs simple linear interpolation. Neither existing nor new wavelengths need to be sorted. """ interpolator = interp1d(msi.get_wavelengths(), msi.get_image(), assume_sorted=False) msi.set_image(interpolator(newWavelengths), wavelengths=newWavelengths) return msi def normalize_integration_times(msi): """ divides by integration times """ if ('integration times' not in msi.get_properties()): logging.warn("warning: trying to normalize integration times for " "image without the integration time property") return msi original_shape = msi.get_image().shape collapsed_image = collapse_image(msi.get_image()) collapsed_image = collapsed_image / msi.get_properties()['integration times'] msi.set_image(collapsed_image.reshape(original_shape)) msi.add_property({'integration times': np.ones_like(msi.get_properties()['integration times'])}) return msi def dark_correction(msi, dark): """" subtract dark current from multi spectral image. The dark msi should either be of the same shape as msi or 1xnr_wavelengths (see tests).""" msi.set_image(msi.get_image() - dark.get_image()) return msi def flatfield_correction(msi, flatfield): """ divide by flatfield to remove dependencies on light source form and imaging system. The flatfield msi should either be of the same shape as msi or 1xnr_wavelengths (see tests).""" # we copy the flatfield to ensure it is unchanged by the normalization flatfield_copy = copy.copy(flatfield) normalize_integration_times(flatfield_copy) normalize_integration_times(msi) msi.set_image(msi.get_image() / flatfield_copy.get_image()) return msi def image_correction(msi, flatfield, dark): """ do the usual image correction: msi = ((msi - dark) / integration_time) / ((flatfield - dark) / integration_time) this function changes only the msi, flatfield and dark image are left unchanged. """ # we need a copy of flatfield since otherwise the dark correction # changes the flatfield dark_correction(msi, dark) flatfield_copy = copy.copy(flatfield) dark_correction(flatfield_copy, dark) flatfield_correction(msi, flatfield_copy) return msi def get_bands(msi, bands): """ TODO SW: document and test """ msi.set_image(imgmani.get_bands(msi.get_image(), bands)) if msi.get_wavelengths() is not None: msi.set_wavelengths(msi.get_wavelengths()[bands]) return msi diff --git a/Modules/Biophotonics/python/iMC/msi/plot.py b/Modules/Biophotonics/python/iMC/msi/plot.py index 41d12fb2d0..0ae653b2ce 100644 --- a/Modules/Biophotonics/python/iMC/msi/plot.py +++ b/Modules/Biophotonics/python/iMC/msi/plot.py @@ -1,84 +1,85 @@ # -*- coding: utf-8 -*- """ Created on Thu Aug 13 11:13:31 2015 @author: wirkert """ import copy +import logging 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() * 100.) + logging.info("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')