From f331c8fbcd2527f842b668762abb30509e9e94ef Mon Sep 17 00:00:00 2001 From: Lorenzo ZOLFANELLI Date: Mon, 25 Jan 2021 18:42:45 +0100 Subject: [PATCH] Add functions to read only rectangular regions. --- nd2reader/parser.py | 106 ++++++++++++++++++++++++++++++++++++++++++++ nd2reader/reader.py | 15 +++++++ 2 files changed, 121 insertions(+) diff --git a/nd2reader/parser.py b/nd2reader/parser.py index 8534e79..8147128 100644 --- a/nd2reader/parser.py +++ b/nd2reader/parser.py @@ -6,6 +6,7 @@ import six import warnings from pims.base_frames import Frame import numpy as np +from tqdm import tqdm from nd2reader.common import get_version, read_chunk from nd2reader.label_map import LabelMap @@ -78,6 +79,38 @@ class Parser(object): else: return Frame(image, frame_no=frame_number, metadata=self._get_frame_metadata()) + def get_slice_by_attributes(self, xywh, frame_number, field_of_view, channel, z_level, height, width): + """Gets a rectangular slice of an image based on its attributes alone + + Args: + xywh: tuples containing (x, y, w, h) values of the + rectangular region to load + frame_number: the frame number + field_of_view: the field of view + channel_name: the color channel name + z_level: the z level + height: the height of the image + width: the width of the image + + Returns: + Frame: the requested image + + """ + frame_number = 0 if frame_number is None else frame_number + field_of_view = 0 if field_of_view is None else field_of_view + channel = 0 if channel is None else channel + z_level = 0 if z_level is None else z_level + + image_group_number = self._calculate_image_group_number(frame_number, field_of_view, z_level) + try: + timestamp, raw_image_data = self._get_raw_slice_data( + xywh, image_group_number, channel, height, width + ) + except (TypeError): + return Frame([], frame_no=frame_number, metadata=self._get_frame_metadata()) + else: + return Frame(raw_image_data, frame_no=frame_number, metadata=self._get_frame_metadata()) + def get_image_by_attributes(self, frame_number, field_of_view, channel, z_level, height, width): """Gets an image based on its attributes alone @@ -246,6 +279,79 @@ class Parser(object): """ return {channel: n for n, channel in enumerate(self.metadata["channels"])} + def _get_raw_slice_data(self, xywh, image_group_number, channel, height, width): + """Reads the raw bytes and the timestamp of a rectangular slice + of an image. + + Args: + xywh: tuples containing (x, y, w, h) values of the + rectangular region to load + image_group_number: the image group number (see _calculate_image_group_number) + channel: the position (int) of the channel to load + height: the height of the image + width: the width of the image + + Returns: + + """ + size_c = len(self.metadata["channels"]) + + x0, y0, w, h = xywh + chunk_location = self._label_map.get_image_data_location(image_group_number) + fh = self._fh + + if chunk_location is None or fh is None: + return None + fh.seek(chunk_location) + # The chunk metadata is always 16 bytes long + chunk_metadata = fh.read(16) + header, relative_offset, data_length = struct.unpack("IIQ", chunk_metadata) + if header != 0xabeceda: + raise ValueError("The ND2 file seems to be corrupted.") + # We start at the location of the chunk metadata, skip over the metadata, and then proceed to the + # start of the actual data field, which is at some arbitrary place after the metadata. + fh.seek(chunk_location + 16 + relative_offset) + + # Read timestamp (8 bytes) + timestamp = struct.unpack("d", fh.read(8))[0] + + # Stitched Images: evaluate number of bytes to strip + # (with stitched images sometimes after each row we have a regular number of extra bytes) + n_unwanted_bytes = (data_length-8) % (height*width) + assert 0 == n_unwanted_bytes % height + rowskip = n_unwanted_bytes // height + + # Read ROI: row-by-row + image_start_pos = chunk_location + 16 + relative_offset + 8 + + line_bytemask = np.zeros(size_c, dtype=np.bool) + line_bytemask[channel] = True + line_bytemask = np.tile(line_bytemask.repeat(2),w) + + def get_line(y): + fh.seek(image_start_pos + size_c*2*((width)*y+x0) + y*rowskip) + return np.frombuffer(fh.read(size_c*2*w), np.byte)[line_bytemask] + + data = [get_line(y) for y in tqdm(range(y0, y0+h))] + data = bytes().join(data) + + image_group_data = array.array("H", data) + true_channels_no = int(len(image_group_data) / (h * w)) + + image_data = np.reshape(image_group_data, (h, w, true_channels_no)) + + missing_channels = ~np.any(image_data, axis=(0, 1)) + image_data[..., missing_channels] = np.full( + (h, w, missing_channels.sum()), np.nan) + + if np.any(missing_channels): + warnings.warn( + "ND2 file contains gap frames which are represented by " + + "np.nan-filled arrays; to convert to zeros use e.g. " + + "np.nan_to_num(array)") + return timestamp, image_data[...,0] + + def _get_raw_image_data(self, image_group_number, channel_offset, height, width): """Reads the raw bytes and the timestamp of an image. diff --git a/nd2reader/reader.py b/nd2reader/reader.py index 4e14f42..295bded 100644 --- a/nd2reader/reader.py +++ b/nd2reader/reader.py @@ -69,6 +69,21 @@ class ND2Reader(FramesSequenceND): except KeyError: return 0 + def get_roi(self, roi, c=0, t=0, z=0, x=0, y=0, v=0): + height = self.metadata['height'] + width = self.metadata['width'] + ylim = roi[0].indices(height) + xlim = roi[1].indices(width) + + y = ylim[0] + x = xlim[0] + w = xlim[1]-xlim[0] + h = ylim[1]-ylim[0] + + return self._parser.get_slice_by_attributes( + (x, y, w, h), t, v, c, z, height, width + ) + def get_frame_2D(self, c=0, t=0, z=0, x=0, y=0, v=0): """Gets a given frame using the parser Args: