diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..95602c0 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,24 @@ +FROM ubuntu +MAINTAINER Jim Rybarski + +RUN apt-get update && apt-get install -y \ + gcc \ + gfortran \ + libblas-dev \ + liblapack-dev \ + libatlas-dev \ + tk \ + tk-dev \ + libpng12-dev \ + python \ + python-dev \ + python-pip \ + libfreetype6-dev \ + python-skimage + +RUN pip install numpy +RUN pip install --upgrade scikit-image + +COPY . /opt/nd2reader +WORKDIR /opt/nd2reader +RUN python setup.py install diff --git a/nd2reader/__init__.py b/nd2reader/__init__.py index 264eb55..4f47559 100644 --- a/nd2reader/__init__.py +++ b/nd2reader/__init__.py @@ -1,26 +1,30 @@ +# -*- coding: utf-8 -*- + +import array +from datetime import datetime import logging -from nd2reader.service import BaseNd2 from nd2reader.model import Image, ImageSet +from nd2reader.parser import Nd2Parser +import re +import struct log = logging.getLogger(__name__) -log.setLevel(logging.DEBUG) +log.addHandler(logging.StreamHandler()) +log.setLevel(logging.WARN) -class Nd2(BaseNd2): - def __init__(self, filename): +class Nd2(Nd2Parser): + def __init__(self, filename, image_sets=False): super(Nd2, self).__init__(filename) - - def get_image(self, time_index, fov, channel_name, z_level): - image_set_number = self._calculate_image_set_number(time_index, fov, z_level) - timestamp, raw_image_data = self._reader.get_raw_image_data(image_set_number, self.channel_offset[channel_name]) - return Image(timestamp, raw_image_data, fov, channel_name, z_level, self.height, self.width) + self._use_image_sets = image_sets def __iter__(self): - """ - Just return every image in order (might not be exactly the order that the images were physically taken, but it will - be within a few seconds). A better explanation is probably needed here. + if self._use_image_sets: + return self.image_sets() + else: + return self.images() - """ + def images(self): for i in range(self._image_count): for fov in range(self.field_of_view_count): for z_level in range(self.z_level_count): @@ -29,19 +33,152 @@ class Nd2(BaseNd2): if image.is_valid: yield image - def image_sets(self, field_of_view, time_indices=None, channels=None, z_levels=None): + def image_sets(self): + for time_index in xrange(self.time_index_count): + image_set = ImageSet() + for fov in range(self.field_of_view_count): + for channel_name in self.channels: + for z_level in xrange(self.z_level_count): + image = self.get_image(time_index, fov, channel_name, z_level) + if image.is_valid: + image_set.add(image) + yield image_set + + def get_image(self, time_index, fov, channel_name, z_level): + image_set_number = self._calculate_image_set_number(time_index, fov, z_level) + timestamp, raw_image_data = self._get_raw_image_data(image_set_number, self._channel_offset[channel_name]) + return Image(timestamp, raw_image_data, fov, channel_name, z_level, self.height, self.width) + + @property + def channels(self): + metadata = self.metadata['ImageMetadataSeq']['SLxPictureMetadata']['sPicturePlanes'] + try: + validity = self.metadata['ImageMetadata']['SLxExperiment']['ppNextLevelEx'][''][0]['ppNextLevelEx'][''][0]['pItemValid'] + except KeyError: + # If none of the channels have been deleted, there is no validity list, so we just make one + validity = [True for _ in metadata] + # Channel information is contained in dictionaries with the keys a0, a1...an where the number + # indicates the order in which the channel is stored. So by sorting the dicts alphabetically + # we get the correct order. + for (label, chan), valid in zip(sorted(metadata['sPlaneNew'].items()), validity): + if not valid: + continue + yield chan['sDescription'] + + @property + def height(self): """ - Gets all the images for a given field of view and + :return: height of each image, in pixels + """ - timepoint_set = xrange(self.time_index_count) if time_indices is None else time_indices - channel_set = [channel.name for channel in self.channels] if channels is None else channels - z_level_set = xrange(self.z_level_count) if z_levels is None else z_levels + return self.metadata['ImageAttributes']['SLxImageAttributes']['uiHeight'] - for timepoint in timepoint_set: - image_set = ImageSet() - for channel_name in channel_set: - for z_level in z_level_set: - image = self.get_image(timepoint, field_of_view, channel_name, z_level) - if image.is_valid: - image_set.add(image) - yield image_set \ No newline at end of file + @property + def width(self): + """ + :return: width of each image, in pixels + + """ + return self.metadata['ImageAttributes']['SLxImageAttributes']['uiWidth'] + + @property + def absolute_start(self): + if self._absolute_start is None: + for line in self.metadata['ImageTextInfo']['SLxImageTextInfo'].values(): + absolute_start_12 = None + absolute_start_24 = None + # ND2s seem to randomly switch between 12- and 24-hour representations. + try: + absolute_start_24 = datetime.strptime(line, "%m/%d/%Y %H:%M:%S") + except ValueError: + pass + try: + absolute_start_12 = datetime.strptime(line, "%m/%d/%Y %I:%M:%S %p") + except ValueError: + pass + if not absolute_start_12 and not absolute_start_24: + continue + self._absolute_start = absolute_start_12 if absolute_start_12 else absolute_start_24 + return self._absolute_start + + @property + def channel_count(self): + pattern = r""".*?λ\((\d+)\).*?""" + try: + count = int(re.match(pattern, self._dimensions).group(1)) + except AttributeError: + return 1 + else: + return count + + @property + def field_of_view_count(self): + """ + The metadata contains information about fields of view, but it contains it even if some fields + of view were cropped. We can't find anything that states which fields of view are actually + in the image data, so we have to calculate it. There probably is something somewhere, since + NIS Elements can figure it out, but we haven't found it yet. + + """ + pattern = r""".*?XY\((\d+)\).*?""" + try: + count = int(re.match(pattern, self._dimensions).group(1)) + except AttributeError: + return 1 + else: + return count + + @property + def time_index_count(self): + """ + The number of image sets. If images were acquired using some kind of cycle, all images at each step in the + program will have the same timestamp (even though they may have varied by a few seconds in reality). For example, + if you have four fields of view that you're constantly monitoring, and you take a bright field and GFP image of + each, and you repeat that process 100 times, you'll have 800 individual images. But there will only be 400 + time indexes. + + :rtype: int + + """ + pattern = r""".*?T'\((\d+)\).*?""" + try: + count = int(re.match(pattern, self._dimensions).group(1)) + except AttributeError: + return 1 + else: + return count + + @property + def z_level_count(self): + pattern = r""".*?Z\((\d+)\).*?""" + try: + count = int(re.match(pattern, self._dimensions).group(1)) + except AttributeError: + return 1 + else: + return count + + @property + def _channel_offset(self): + """ + Image data is interleaved for each image set. That is, if there are four images in a set, the first image + will consist of pixels 1, 5, 9, etc, the second will be pixels 2, 6, 10, and so forth. Why this would be the + case is beyond me, but that's how it works. + + """ + channel_offset = {} + for n, channel in enumerate(self.channels): + self._channel_offset[channel.name] = n + return channel_offset + + def _get_raw_image_data(self, image_set_number, channel_offset): + chunk = self._label_map["ImageDataSeq|%d!" % image_set_number] + data = self._read_chunk(chunk) + timestamp = struct.unpack("d", data[:8])[0] + # The images for the various channels are interleaved within each other. + image_data = array.array("H", data) + image_data_start = 4 + channel_offset + return timestamp, image_data[image_data_start::self.channel_count] + + def _calculate_image_set_number(self, time_index, fov, z_level): + return time_index * self.field_of_view_count * self.z_level_count + (fov * self.z_level_count + z_level) diff --git a/nd2reader/model/__init__.py b/nd2reader/model/__init__.py index 1ae84e4..800cb7b 100644 --- a/nd2reader/model/__init__.py +++ b/nd2reader/model/__init__.py @@ -1,33 +1,10 @@ import numpy as np import skimage.io import logging -from io import BytesIO -import array -import struct - log = logging.getLogger(__name__) -class Channel(object): - def __init__(self, name, camera, exposure_time): - self._name = name - self._camera = camera - self._exposure_time = exposure_time - - @property - def name(self): - return self._name - - @property - def camera(self): - return self._camera - - @property - def exposure_time(self): - return self._exposure_time - - class ImageSet(object): """ A group of images that share the same timestamp. NIS Elements doesn't store a unique timestamp for every @@ -96,172 +73,4 @@ class Image(object): def show(self): skimage.io.imshow(self.data) - skimage.io.show() - - -class MetadataItem(object): - def __init__(self, start, data): - self._datatype = ord(data[start]) - self._label_length = 2 * ord(data[start + 1]) - self._data = data - - @property - def is_valid(self): - return self._datatype > 0 - - @property - def key(self): - return self._data[2:self._label_length].decode("utf16").encode("utf8") - - @property - def length(self): - return self._length - - @property - def data_start(self): - return self._label_length + 2 - - @property - def _body(self): - """ - All data after the header. - - """ - return self._data[self.data_start:] - - def _get_bytes(self, count): - return self._data[self.data_start: self.data_start + count] - - @property - def value(self): - parser = {1: self._parse_unsigned_char, - 2: self._parse_unsigned_int, - 3: self._parse_unsigned_int, - 5: self._parse_unsigned_long, - 6: self._parse_double, - 8: self._parse_string, - 9: self._parse_char_array, - 11: self._parse_metadata_item - } - return parser[self._datatype]() - - def _parse_unsigned_char(self): - self._length = 1 - return self._unpack("B", self._get_bytes(self._length)) - - def _parse_unsigned_int(self): - self._length = 4 - return self._unpack("I", self._get_bytes(self._length)) - - def _parse_unsigned_long(self): - self._length = 8 - return self._unpack("Q", self._get_bytes(self._length)) - - def _parse_double(self): - self._length = 8 - return self._unpack("d", self._get_bytes(self._length)) - - def _parse_string(self): - # the string is of unknown length but ends at the first instance of \x00\x00 - stop = self._body.index("\x00\x00") - self._length = stop - return self._body[:stop - 1].decode("utf16").encode("utf8") - - def _parse_char_array(self): - array_length = self._unpack("Q", self._get_bytes(8)) - self._length = array_length + 8 - return array.array("B", self._body[8:array_length]) - - def _parse_metadata_item(self): - count, length = struct.unpack("