Source code for flimage.series

import h5py

import lmfit
import numpy as np
from skimage.restoration import denoise_tv_chambolle

from .core import FLImage
from .meta import FLMetaDict

[docs]class FLSeries(object): _instances = 0 def __init__(self, flimage_list=[], meta_data={}, h5file=None, h5mode="a", identifier=None): """Fluorescence microscopy series data Parameters ---------- flimage_list: list A list of instances of :class:`flimage.FLImage`. meta_data: dict Meta data associated with the input data (see :const:`flimage.META_KEYS`). This overrides the meta data of the FLImages in `flimage_list` and, if `h5file` is given and `h5mode` is not "r", overrides the meta data in `h5file`. h5file: str, h5py.Group, h5py.File, or None A path to an hdf5 data file where all data is cached. If set to `None` (default), all data will be handled in memory using the "core" driver of the :mod:`h5py`'s :class:`h5py:File` class. If the file does not exist, it is created. If the file already exists, it is opened with the file mode defined by `hdf5_mode`. If this is an instance of h5py.Group or h5py.File, then this will be used to internally store all data. If `h5file` is given and `flimage_list` is not empty, all FLImages in `flimage_list` are appended to `h5file` in the given order. h5mode: str Valid file modes are (only applies if `h5file` is a path): - "r": Readonly, file must exist - "r+": Read/write, file must exist - "w": Create file, truncate if exists - "w-" or "x": Create file, fail if exists - "a": Read/write if exists, create otherwise (default) """ if flimage_list and not isinstance(flimage_list, list): msg = "`flimage_list` must be a list!" if isinstance(flimage_list, str): msg += " Did you mean `h5file={}`?".format(flimage_list) raise ValueError(msg) if isinstance(h5file, h5py.Group): self.h5 = h5file self._do_h5_cleanup = False else: if h5file is None: h5kwargs = {"name": "flseries{}.h5".format( FLSeries._instances), "driver": "core", "backing_store": False, "mode": "w"} else: h5kwargs = {"name": h5file, "mode": h5mode} self.h5 = h5py.File(**h5kwargs) self._do_h5_cleanup = True FLSeries._instances += 1 if meta_data and h5mode == "r": msg = "`h5mode` must not be 'r' if `meta_data` is given!" raise ValueError(msg) # make sure self.h5 is not itself a FLImage file if "fluorescence" in self.h5: raise ValueError( "`h5file` is an FLImage file, not an FLSeries file!") # Write QPimage data to h5 file for fli in flimage_list: self.add_flimage(fli) # Update meta data if meta_data: meta = FLMetaDict(meta_data) for ii in range(len(self)): flii = self.get_flimage(index=ii) for mk in meta: flii.h5.attrs[mk] = meta[mk] # Set identifier if identifier: self.h5.attrs["identifier"] = identifier def __contains__(self, flid): """test whether an FLImage with the given identifier exists""" for ii in range(len(self)): fli = self[ii] if "identifier" in fli and fli["identifier"] == flid: exists = True break else: exists = False return exists def __enter__(self): return self def __exit__(self, exc_type, exc_val, exc_tb): if self._do_h5_cleanup: self.h5.flush() self.h5.close() def __getitem__(self, index): return self.get_flimage(index) def __iter__(self): for ii in range(len(self)): yield self[ii] def __len__(self): keys = list(self.h5.keys()) keys = [kk for kk in keys if kk.startswith("fli_")] return len(keys) @property def identifier(self): """unique identifier of the series""" if "identifier" in self.h5.attrs: return self.h5.attrs["identifier"] else: return None
[docs] def add_flimage(self, fli, identifier=None): """Add a FLImage instance to the FLSeries Parameters ---------- fli: flimage.FLImage The FLImage that is added to the series identifier: str Identifier key for `fli` """ if not isinstance(fli, FLImage): raise ValueError("`fli` must be instance of FLImage!") if "identifier" in fli and identifier is None: identifier = fli["identifier"] if identifier and identifier in self: msg = "The identifier '{}' already ".format(identifier) \ + "exists! You can either change the identifier of " \ + " '{}' or remove it.".format(fli) raise ValueError(msg) # determine number of flimages num = len(self) # indices start at zero; do not add 1 name = "fli_{}".format(num) group = self.h5.create_group(name) fli.copy(h5file=group) if identifier: # set identifier group.attrs["identifier"] = identifier
[docs] def bleach_correction(self, h5out=None, border_px=20, flscorr=None, count=None, max_count=None): """Perform a correction for photobleaching Bleaching is modeled with an exponential and an offset. Parameters ---------- h5out: path or h5py.Group A new FLSeries will be written to this HDF5 file or h5py group. border_px: int Number of border pixels to include for background estimation. flscorr: flimage.series.FLSeries Apply the background correction to this FLSeries instead of the current instance. Use this in combination with filtered versions of the same series. Returns ------- bg: float Background value subtracted from each image flint: 1d ndarray Fluorescence intensity trace extracted from the series decay: 1d ndarray Exponential fit to `flint` times: 1d ndarray Recording times corresponding to the indices in `flint` and `decay` Notes ----- It is recommended to first denoise the fluorescence data with `FLSeries.denoise` do avoid an amplification of the background noise. """ if max_count is not None: max_count.value += 3*len(self) if flscorr is None: flscorr = self # get border data border = np.ones(self[0].shape, dtype=bool) border[border_px:-border_px, border_px:-border_px] = False bgavg = np.zeros(len(self)) for ii, fli in enumerate(self): bgavg[ii] = np.mean(fli.fl[border]) if count is not None: count.value += 1 bg = np.mean(bgavg) # intensity values flint = np.zeros(len(self)) for ii, fli in enumerate(self): fl = fli.fl - bg flint[ii] = np.sum(fl[fl > 0]) if count is not None: count.value += 1 # time values times = np.array([fli["time"] for fli in self]) # fit exponential decay = fit_exponential(times, flint) corr = decay[0] / decay if h5out: # perform bleach correction with FLSeries(h5file=h5out) as fls: for ii, fli in enumerate(flscorr): fld = (fli.fl - bg) * corr[ii] fln = FLImage(data=fld, meta_data=fli.meta) fls.add_flimage(fln) if count is not None: count.value += 1 return bg, flint, decay, times
[docs] def denoise(self, h5file, count=None, max_count=None): """Denoise fluorescence data `h5file` specifies a path or HDF5 group for the FLSeries to which the denoised fluorescence data is written.""" if max_count is not None: max_count.value += len(self) with FLSeries(h5file=h5file) as fls: for fli in self: fldn = denoise_tv_chambolle(fli.fl, weight=2) fln = FLImage(data=fldn, meta_data=fli.meta) fls.add_flimage(fln) if count is not None: count.value += 1
[docs] def get_flimage(self, index): """Return a single FLImage of the series Parameters ---------- index: int or str Index or identifier of the FLImage Notes ----- Instead of ``fls.get_flimage(index)``, it is possible to use the short-hand ``fls[index]``. """ if isinstance(index, str): # search for the identifier for ii in range(len(self)): fli = self[ii] if "identifier" in fli and fli["identifier"] == index: group = self.h5["fli_{}".format(ii)] break else: msg = "FLImage identifier '{}' not found!".format(index) raise KeyError(msg) else: # integer index if index < -len(self): msg = "Index {} out of bounds for FLSeries of size {}!".format( index, len(self)) raise ValueError(msg) elif index < 0: index += len(self) name = "fli_{}".format(index) if name in self.h5: group = self.h5[name] else: msg = "Index {} not found for FLSeries of length {}".format( index, len(self)) raise KeyError(msg) return FLImage(h5file=group)
def fit_exponential(times, flint): def model_fct(params, t): ampl = params["amplitude"].value decay = params["bleach_decay"].value offset = params["offset"].value model = ampl*np.exp(-decay*t) + offset return model def residual(params, t, data): return model_fct(params, t)-data params = lmfit.Parameters() params.add("amplitude", value=flint.max()-flint.min(), min=flint.min(), max=2*flint.max()) dt = times[-1] - times[0] df = flint[-1] - flint[0] params.add("bleach_decay", value=-dt/df*flint[0]/20, min=0, brute_step=dt/20) params.add("offset", value=flint.min(), min=0) out = lmfit.minimize(residual, params, args=(times, flint)) model = model_fct(out.params, times) return model