Source code for colicoords.data_models

import mahotas as mh
import numpy as np
import math
from scipy.ndimage.interpolation import rotate as scipy_rotate


# https://stackoverflow.com/questions/26598109/preserve-custom-attributes-when-pickling-subclass-of-numpy-array
[docs]class BinaryImage(np.ndarray): """ Binary image data class. Attributes ---------- name : :obj:`str` Name identifying the data element. metadata : :obj:`dict` Optional dict for metadata, load/save not implemented. """ def __new__(cls, input_array, name=None, metadata=None): if input_array is None: return None obj = np.asarray(input_array).view(cls) obj.name = name obj.metadata = metadata obj.dclass = 'binary' return obj def __array_finalize__(self, obj): if obj is None: return self.name = getattr(obj, 'name', None) self.metadata = getattr(obj, 'metadata', None) self.dclass = getattr(obj, 'dclass', 'binary') def __array_wrap__(self, out_arr, context=None): return np.ndarray.__array_wrap__(self, out_arr, context) def __reduce__(self): # Get the parent's __reduce__ tuple pickled_state = super(BinaryImage, self).__reduce__() # Create our own tuple to pass to __setstate__ new_state = pickled_state[2] + (self.name, self.metadata, self.dclass) # Return a tuple that replaces the parent's __setstate__ tuple with our own return pickled_state[0], pickled_state[1], new_state def __setstate__(self, state): self.name, self.metadata, self.dclass = state[-3:] super(BinaryImage, self).__setstate__(state[0:-3]) @property def orientation(self): """:obj:`float`: The main image axis orientation in degrees""" return _calc_orientation(self)
[docs]class BrightFieldImage(np.ndarray): """ Brightfield image data class. Attributes ---------- name : :obj:`str` Name identifying the data element. metadata : :obj:`dict` Optional dict for metadata, load/save not implemented. """ def __new__(cls, input_array, name=None, metadata=None): if input_array is None: return None obj = np.asarray(input_array).view(cls) obj.name = name obj.metadata = metadata obj.dclass = 'brightfield' return obj def __array_finalize__(self, obj): if obj is None: return self.name = getattr(obj, 'name', None) self.metadata = getattr(obj, 'metadata', None) self.dclass = getattr(obj, 'dclass', 'brightfield') def __array_wrap__(self, out_arr, context=None): return np.ndarray.__array_wrap__(self, out_arr, context) def __reduce__(self): # Get the parent's __reduce__ tuple pickled_state = super(BrightFieldImage, self).__reduce__() # Create our own tuple to pass to __setstate__ new_state = pickled_state[2] + (self.name, self.metadata, self.dclass) # Return a tuple that replaces the parent's __setstate__ tuple with our own return (pickled_state[0], pickled_state[1], new_state) def __setstate__(self, state): self.name, self.metadata, self.dclass = state[-3:] super(BrightFieldImage, self).__setstate__(state[0:-3]) @property def orientation(self): """:obj:`float`: The main image axis orientation in degrees""" return _calc_orientation(self)
[docs]class FluorescenceImage(np.ndarray): """ Fluorescence image data class. Attributes ---------- name : :obj:`str` Name identifying the data element. metadata : :obj:`dict` Optional dict for metadata, load/save not implemented. """ def __new__(cls, input_array, name=None, metadata=None): if input_array is None: return None obj = np.asarray(input_array).view(cls) obj.name = name obj.metadata = metadata obj.dclass = 'fluorescence' return obj def __array_finalize__(self, obj): if obj is None: return self.name = getattr(obj, 'name', None) self.metadata = getattr(obj, 'metadata', None) self.dclass = getattr(obj, 'dclass', 'fluorescence') def __array_wrap__(self, out_arr, context=None): return np.ndarray.__array_wrap__(self, out_arr, context) def __reduce__(self): # Get the parent's __reduce__ tuple pickled_state = super(FluorescenceImage, self).__reduce__() # Create our own tuple to pass to __setstate__ new_state = pickled_state[2] + (self.name, self.metadata, self.dclass) # Return a tuple that replaces the parent's __setstate__ tuple with our own return (pickled_state[0], pickled_state[1], new_state) def __setstate__(self, state): self.name, self.metadata, self.dclass = state[-3:] super(FluorescenceImage, self).__setstate__(state[0:-3]) @property def orientation(self): """:obj:`float`: The main image axis orientation in degrees""" return _calc_orientation(self)
[docs]class STORMTable(np.ndarray): """ STORM table data class. Attributes ---------- name : :obj:`str` Name identifying the data element. metadata : :obj:`dict` Optional dict for metadata, load/save not implemented. """ def __new__(cls, input_array, name=None, metadata=None): if input_array is None: return None obj = np.asarray(input_array).view(cls) obj.name = name obj.metadata = metadata obj.dclass = 'storm' return obj def __array_finalize__(self, obj): if obj is None: return self.name = getattr(obj, 'name', None) self.metadata = getattr(obj, 'metadata', None) self.dclass = getattr(obj, 'dclass', 'storm') def __array_wrap__(self, out_arr, context=None): return np.ndarray.__array_wrap__(self, out_arr, context) def __reduce__(self): # Get the parent's __reduce__ tuple pickled_state = super(STORMTable, self).__reduce__() # Create our own tuple to pass to __setstate__ new_state = pickled_state[2] + (self.name, self.metadata, self.dclass) # Return a tuple that replaces the parent's __setstate__ tuple with our own return (pickled_state[0], pickled_state[1], new_state) def __setstate__(self, state): self.name, self.metadata, self.dclass = state[-3:] super(STORMTable, self).__setstate__(state[0:-3])
#todo this shoud be a dict? (use open microscopy format?) (XML)
[docs]class MetaData(dict): pass
[docs]class Data(object): """ Main data class holding data from different input channels. The data class is designed to combine and organize all different channels (brightfield, binary, fluorescence, storm) into one object. The class provides basic functionality such as rotation and slicing. Data elements can be accessed from `data_dict` or by attribute '<class>_<name>', where class can be either 'flu', 'storm'. Binary and brightfield can bre accessed as properties. Attributes ---------- data_dict : :obj:`dict` Dictionary with all data elements by their name. flu_dict : :obj:`dict` Subset of `data_dict` with all Fluorescence data elements. storm_dict : :obj:`dict` Subset of `data_dict` with all STORM data elements. """ def __init__(self): self.data_dict = {} self.flu_dict = {} self.bf_dict = {} self.storm_dict = {} self.shape = None self.ndim = None self._idx = 0
[docs] def add_data(self, data, dclass, name=None, metadata=None): """ Add data to form a new data element. Parameters ---------- data : array_like Input data. Either np.ndarray with ndim 2 or 3 (images / movies) or numpy structured array for STORM data. dclass : :obj:`str` Data class. Must be either 'binary', 'brightfield', 'fluorescence' or 'storm'. name : :obj:`str`, optional The name to identify the data element. Default is equal to the data class. metadata : :obj:`dict` Associated metadata (load/save metadata currently not supported) """ if name in ['', u'', r'', None]: name = dclass else: name = str(name) if name in self.names: raise ValueError('Data element name {} is already used'.format(name)) if dclass == 'binary': if self.binary_img is not None: raise ValueError('Binary image has to be unique and is already given') if not np.issubdtype(data.dtype, np.integer): raise TypeError('Invalid data type {} for data class binary'.format(data.dtype)) self._check_shape(data.shape, data.ndim) self.data_dict[name] = BinaryImage(data.astype(np.int32), name=name, metadata=metadata) elif dclass == 'brightfield': self._check_shape(data.shape, data.ndim) b = BrightFieldImage(data, name=name, metadata=metadata) self.bf_dict[name] = b elif dclass == 'fluorescence': assert name assert name not in self.flu_dict self._check_shape(data.shape, data.ndim) f = FluorescenceImage(data, name=name, metadata=metadata) setattr(self, 'flu_' + name, f) self.flu_dict[name] = f elif dclass == 'storm': assert name assert name not in self.storm_dict for field in ['x', 'y', 'frame']: assert field in data.dtype.names s = STORMTable(data, name=name, metadata=metadata) self.storm_dict[name] = s setattr(self, 'storm_' + name, s) else: raise ValueError('Invalid data class {}'.format(dclass)) self.data_dict.update(self.bf_dict) self.data_dict.update(self.flu_dict) self.data_dict.update(self.storm_dict)
[docs] def prune(self, data_name): """ Removes localizations from the STORM-dataset with name `data_name` which lie outside of the associated image. Parameters ---------- data_name : :obj:`str` Name of the data element to prune. Returns ------- None """ #todo np.clip? storm = self.data_dict.pop(data_name) self.storm_dict.pop(data_name) assert isinstance(storm, STORMTable) xmax, ymax = self.shape[1], self.shape[0] bools = (storm['x'] < 0) + (storm['x'] > xmax) + (storm['y'] < 0) + (storm['y'] > ymax) storm_out = storm[bools].copy() self.add_data(storm_out, storm.dclass, name=data_name)
[docs] def copy(self): """ Copy the data object. Returns ------- data : :class:`~colicoords.data_models.Data` Copied data object. """ data = Data() for v in self.data_dict.values(): data.add_data(np.copy(v), v.dclass, name=v.name, metadata=v.metadata) return data
@property def dclasses(self): """:obj:`list`: List of all data classes in the ``Data`` object.""" return [d.dclass for d in self.data_dict.values()] @property def names(self): """:obj:`list`: List of all data names in the ``Data`` object.""" return [d.name for d in self.data_dict.values()] @property def binary_img(self): """:class:`~np.ndarray`: Returns the binary image if present, else ``None``""" try: return self.data_dict['binary'] except KeyError: return None @property def bf_img(self): """:class:`~np.ndarray`: Returns the brightfield image if present, else ``None``""" try: return self.data_dict['brightfield'] except KeyError: return None
[docs] def rotate(self, theta): """ Rotate all data elements and return a new ``Data`` object with rotated data elements. Parameters ---------- theta : :obj:`float` Rotation angle in degrees. Returns ------- data : :class:`colicoords.data_models.Data` Rotated ``Data`` """ data = Data() for v in self.data_dict.values(): if v.dclass == 'storm': rotated = _rotate_storm(v, -theta, shape=self.shape) else: rotated = scipy_rotate(v, -theta, mode='nearest', axes=(-2, -1)) #todo check dis data.add_data(rotated, v.dclass, name=v.name, metadata=v.metadata) return data
def _check_shape(self, shape, ndim): if self.shape: if not ((shape == self.shape) and (ndim == self.ndim)): if not ((ndim == self.ndim + 1) and (shape[1:] == self.shape)): raise ValueError("Invalid shape") else: self.shape = shape self.ndim = ndim def __len__(self): if not hasattr(self, 'ndim'): return 0 elif self.ndim == 3: return self.shape[0] elif self.ndim == 2: return 1 else: raise ValueError def __iter__(self): self._idx = 0 return self def __next__(self): if not hasattr(self, 'ndim'): raise StopIteration if self.ndim == 2: if self._idx == 0: self._idx += 1 return self else: self._idx = 0 raise StopIteration if self._idx >= len(self): self._idx = 0 raise StopIteration else: data = self[self._idx] self._idx += 1 return data def __getitem__(self, key): #todo len (key) means 2d slicing! data = Data() for v in self.data_dict.values(): if v.dclass == 'storm': #todo needs testing if type(key) == int: #Select the appropriate frames, STORM frame numbers start counting at 1 xmin = 0 ymin = 0 b_z = v['frame'] == key + 1 b_xy = True elif len(key) == 3: #3d slicing, slices the frames? #todo 3d slicing by frame! raise NotImplementedError() elif len(key) == 2: ymin, ymax, ystep = key[0].start, key[0].stop, key[0].step xmin, xmax, xstep = key[1].start, key[1].stop, key[1].step xmin = xmin if xmin else 0 ymin = ymin if ymin else 0 xmax = xmax if xmax else self.shape[1] ymax = ymax if ymax else self.shape[0] if ystep is not None or xstep is not None: raise ValueError('Cannot specify slice steps for slicing images in x and y dimenions') #Create boolean array to mask entries withing the chosen range b_xy = (v['x'] > xmin) * (v['x'] < xmax) * (v['y'] > ymin) * (v['y'] < ymax) b_z = True # Choose selected data and copy, rezero x and y b_overall = b_z * b_xy table_out = v[b_overall].copy() table_out['x'] -= xmin table_out['y'] -= ymin data.add_data(table_out, v.dclass, name=v.name, metadata=v.metadata) else: #key = (slice(None), *key) if v.ndim == 3 and type(key) == tuple else key #print(key)#todo this needs maaaasive testing data.add_data(v[key], v.dclass, name=v.name, metadata=v.metadata) return data next = __next__
[docs]class CellListData(object): """ Data class for CellList with common attributes for all cells. Individual data elements are accessed per cell. Parameters ---------- cell_list : :obj:`list` or :class:`numpy.ndarray` List of array of :class:`~colicoords.cell.Cell` objects. """ def __init__(self, cell_list): self.cell_list = cell_list @property def shape(self): """:obj:`tuple`: Tuple of cell's data element's shape if they are all equal, else `None`""" sh0, sh1 = np.array([c.data.shape for c in self.cell_list]).T if np.all(sh0[0] == sh0) and np.all(sh1[0] == sh1): return sh0[0], sh1[1] else: return None @property def names(self): """:obj:`list`: List of all data names in the ``Data`` objects of the cells, if all are equal, else `None`.""" names = np.array([c.data.names for c in self.cell_list]) if np.all(names[0] == names): return list(names[0]) else: return None @property def dclasses(self): """:obj:`list`: List of all data classes in the ``Data`` objects of the cells, if all are equal, else `None`.""" dclasses = np.array([c.data.dclasses for c in self.cell_list]) if np.all(dclasses[0] == dclasses): return list(dclasses[0]) else: return None
def _rotate_storm(storm_data, theta, shape=None): theta *= np.pi / 180 # to radians x = storm_data['x'].copy() y = storm_data['y'].copy() if shape: ymax = shape[0] xmax = shape[1] ynew = np.abs(xmax * np.sin(-theta)) + np.abs(ymax * np.cos(-theta)) xnew = np.abs(xmax * np.cos(-theta)) + np.abs(ymax * np.sin(-theta)) else: xmax = int(storm_data['x'].max()) + 2 ymax = int(storm_data['y'].max()) + 2 ynew = 0 xnew = 0 x -= xmax / 2 y -= ymax / 2 xr = x * np.cos(theta) + y * np.sin(theta) yr = y * np.cos(theta) - x * np.sin(theta) xr += xnew / 2 yr += ynew / 2 storm_out = np.copy(storm_data) storm_out['x'] = xr storm_out['y'] = yr return storm_out def _calc_orientation(img): com = mh.center_of_mass(img) mu00 = mh.moments(img, 0, 0, com) mu11 = mh.moments(img, 1, 1, com) mu20 = mh.moments(img, 2, 0, com) mu02 = mh.moments(img, 0, 2, com) mup_20 = mu20 / mu00 mup_02 = mu02 / mu00 mup_11 = mu11 / mu00 theta_rad = 0.5 * math.atan(2 * mup_11 / (mup_20 - mup_02)) # todo math -> numpy theta = theta_rad * (180 / math.pi) if (mup_20 - mup_02) > 0: theta += 90 return theta