Source code for allensdk.brain_observatory.roi_masks

# Allen Institute Software License - This software license is the 2-clause BSD
# license plus a third clause that prohibits redistribution for commercial
# purposes without further permission.
#
# Copyright 2017. Allen Institute. All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright notice,
# this list of conditions and the following disclaimer in the documentation
# and/or other materials provided with the distribution.
#
# 3. Redistributions for commercial purposes are not permitted without the
# Allen Institute's written permission.
# For purposes of this license, commercial purposes is the incorporation of the
# Allen Institute's software into anything for which you will charge fees or
# other compensation. Contact terms@alleninstitute.org for commercial licensing
# opportunities.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
#
import numpy as np
import math
import scipy.ndimage.morphology as morphology
import logging
import h5py

# constants used for accessing border array
RIGHT_SHIFT = 0
LEFT_SHIFT = 1
DOWN_SHIFT = 2
UP_SHIFT = 3


[docs]class Mask(object): ''' Abstract class to represent image segmentation mask. Its two main subclasses are RoiMask and NeuropilMask. The former represents the mask of a region of interest (ROI), such as a cell observed in 2-photon imaging. The latter represents the neuropil around that cell, and is useful when subtracting the neuropil signal from the measured ROI signal. This class should not be instantiated directly. Parameters ---------- image_w: integer Width of image that ROI resides in image_h: integer Height of image that ROI resides in label: text User-defined text label to identify mask mask_group: integer User-defined number to help put masks into different categories ''' def __init__(self, image_w, image_h, label, mask_group): ''' Mask class constructor. The Mask class is designed to be abstract and it should not be instantiated directly. ''' self.img_rows = image_h self.img_cols = image_w # initialize to invalid state. Mask must be manually initialized # by pixel list or mask array self.x = 0 self.width = 0 self.y = 0 self.height = 0 self.mask = None self.overlaps_motion_border = False # label is for distinguishing neuropil from ROI, in case # these masks are mixed together self.label = label # auxiliary metadata. if a particula mask is part of an group, # that data can be stored here self.mask_group = mask_group def __str__(self): return "%s: TL=%d,%d w,h=%d,%d\n%s" % (self.label, self.x, self.y, self.width, self.height, str(self.mask))
[docs] def init_by_pixels(self, border, pix_list): ''' Initialize mask using a list of mask pixels Parameters ---------- border: float[4] Coordinates defining useable area of image. See create_roi_mask() pix_list: integer[][2] List of pixel coordinates (x,y) that define the mask ''' assert pix_list.shape[1] == 2, "Pixel list not properly formed" array = np.zeros((self.img_rows, self.img_cols), dtype=bool) # pix_list stores array of [x,y] coordinates array[pix_list[:, 1], pix_list[:, 0]] = 1 self.init_by_mask(border, array)
[docs] def get_mask_plane(self): ''' Returns mask content on full-size image plane Returns ------- numpy 2D array [img_rows][img_cols] ''' mask = np.zeros((self.img_rows, self.img_cols)) mask[self.y:self.y + self.height, self.x:self.x + self.width] = self.mask return mask
[docs]def create_roi_mask(image_w, image_h, border, pix_list=None, roi_mask=None, label=None, mask_group=-1): ''' Conveninece function to create and initializes an RoiMask Parameters ---------- image_w: integer Width of image that ROI resides in image_h: integer Height of image that ROI resides in border: float[4] Coordinates defining useable area of image. If the entire image is usable, and masks are valid anywhere in the image, this should be [(image_w-1), 0, (image_h-1), 0]. The following constants help describe the array order: RIGHT_SHIFT = 0 LEFT_SHIFT = 1 DOWN_SHIFT = 2 UP_SHIFT = 3 When parts of the image are unusable, for example due motion correction shifting of different image frames, the border array should store the usable image area pix_list: integer[][2] List of pixel coordinates (x,y) that define the mask roi_mask: integer[image_h][image_w] Image-sized array that describes the mask. Active parts of the mask should have values >0. Background pixels must be zero label: text User-defined text label to identify mask mask_group: integer User-defined number to help put masks into different categories Returns ------- RoiMask object ''' m = RoiMask(image_w, image_h, label, mask_group) if pix_list is not None: m.init_by_pixels(border, pix_list) elif roi_mask is not None: m.init_by_mask(border, roi_mask) else: assert False, "Must specify either roi_mask or pix_list" return m
[docs]class RoiMask(Mask): def __init__(self, image_w, image_h, label, mask_group): ''' RoiMask class constructor Parameters ---------- image_w: integer Width of image that ROI resides in image_h: integer Height of image that ROI resides in label: text User-defined text label to identify mask mask_group: integer User-defined number to help put masks into different categories ''' super(RoiMask, self).__init__(image_w, image_h, label, mask_group)
[docs] def init_by_mask(self, border, array): ''' Initialize mask using spatial mask Parameters ---------- border: float[4] Coordinates defining useable area of image. See create_roi_mask(). roi_mask: integer[image height][image width] Image-sized array that describes the mask. Active parts of the mask should have values >0. Background pixels must be zero ''' # find lowest and highest non-zero indices on each axis px = np.argwhere(array) (top, left), (bottom, right) = px.min(0), px.max(0) # left and right border insets l_inset = math.ceil(border[RIGHT_SHIFT]) r_inset = math.floor(self.img_cols - border[LEFT_SHIFT]) - 1 # top and bottom border insets t_inset = math.ceil(border[DOWN_SHIFT]) b_inset = math.floor(self.img_rows - border[UP_SHIFT]) - 1 # if ROI crosses border, it's considered invalid if left < l_inset or right > r_inset: self.overlaps_motion_border = True if top < t_inset or bottom > b_inset: self.overlaps_motion_border = True # self.x = left self.width = right - left + 1 self.y = top self.height = bottom - top + 1 # make copy of mask self.mask = array[top:bottom + 1, left:right + 1]
[docs]def create_neuropil_mask(roi, border, combined_binary_mask, label=None): ''' Conveninece function to create and initializes a Neuropil mask. Neuropil masks are defined as the region around an ROI, up to 13 pixels out, that does not include other ROIs Parameters ---------- roi: RoiMask object The ROI that the neuropil masks will be based on border: float[4] Coordinates defining useable area of image. See create_roi_mask(). combined_binary_mask List of pixel coordinates (x,y) that define the mask combined_binary_mask: integer[image_h][image_w] Image-sized array that shows the position of all ROIs in the image. ROI masks should have a value of one. Background pixels must be zero. In other words, ithe combined_binary_mask is a bitmap union of all ROI masks label: text User-defined text label to identify the mask Returns ------- NeuropilMask object ''' # combined_binary_mask is a bitmap union of ALL ROI masks # create a binary mask of the ROI binary_mask = np.zeros((roi.img_rows, roi.img_cols)) binary_mask[roi.y:roi.y + roi.height, roi.x:roi.x + roi.width] = roi.mask binary_mask = binary_mask > 0 # dilate the mask binary_mask_dilated = morphology.binary_dilation( binary_mask, structure=np.ones((3, 3)), iterations=13) # T/F # eliminate ROIs from the dilation binary_mask_dilated = binary_mask_dilated > combined_binary_mask # create mask from binary dilation m = NeuropilMask(w=roi.img_cols, h=roi.img_rows, label=label, mask_group=roi.mask_group) m.init_by_mask(border, binary_mask_dilated) return m
[docs]class NeuropilMask(Mask): def __init__(self, w, h, label, mask_group): ''' NeuropilMask class constructor. This class should be created by calling create_neuropil_mask() Parameters ---------- label: text User-defined text label to identify mask mask_group: integer User-defined number to help put masks into different categories ''' super(NeuropilMask, self).__init__(w, h, label, mask_group)
[docs] def init_by_mask(self, border, array): ''' Initialize mask using spatial mask Parameters ---------- border: float[4] Coordinates defining useable area of image. See create_roi_mask(). array: integer[image height][image width] Image-sized array that describes the mask. Active parts of the mask should have values >0. Background pixels must be zero ''' # find lowest and highest non-zero indices on each axis px = np.argwhere(array) (top, left), (bottom, right) = px.min(0), px.max(0) # left and right border insets l_inset = math.ceil(border[RIGHT_SHIFT]) r_inset = math.floor(self.img_cols - border[LEFT_SHIFT]) - 1 # top and bottom border insets t_inset = math.ceil(border[DOWN_SHIFT]) b_inset = math.floor(self.img_rows - border[UP_SHIFT]) - 1 # restrict neuropil masks to center area of frame (ie, exclude # areas that overlap with movement correction buffer) if left < l_inset: left = l_inset if right < l_inset: right = l_inset if right > r_inset: right = r_inset if left > r_inset: left = r_inset if top < t_inset: top = t_inset if bottom < t_inset: bottom = t_inset if bottom > b_inset: bottom = b_inset if top > b_inset: top = b_inset # self.x = left self.width = right - left + 1 self.y = top self.height = bottom - top + 1 # make copy of mask self.mask = array[top:bottom + 1, left:right + 1]
[docs]def calculate_traces(stack, mask_list, block_size=100): ''' Calculates the average response of the specified masks in the image stack Parameters ---------- stack: float[image height][image width] Image stack that masks are applied to mask_list: list<Mask> List of masks Returns ------- float[number masks][number frames] This is the average response for each Mask in each image frame ''' traces = np.zeros((len(mask_list), stack.shape[0]), dtype=float) num_frames = stack.shape[0] # make sure masks are numpy objects mask_areas = np.zeros(len(mask_list), dtype=float) valid_masks = np.ones(len(mask_list), dtype=bool) for i,mask in enumerate(mask_list): if not isinstance(mask.mask, np.ndarray): mask.mask = np.array(mask.mask) # compute mask areas mask_areas[i] = mask.mask.sum() # if the mask is empty, the trace is nan if mask_areas[i] == 0: logging.warning("mask '%d/%s' is empty", i, mask.label) traces[i,:] = np.nan valid_masks[i] = False # if the mask overlaps the motion border, the trace is nan if mask.overlaps_motion_border: logging.warning("mask '%d/%s' overlaps with motion border", i, mask.label) traces[i,:] = np.nan valid_masks[i] = False # calculate traces for frame_num in range(0, num_frames, block_size): if frame_num % 1000 == 0: logging.debug("frame " + str(frame_num) + " of " + str(num_frames)) frames = stack[frame_num:frame_num+block_size] for i in range(len(mask_list)): if valid_masks[i] is False: continue mask = mask_list[i] subframe = frames[:,mask.y:mask.y + mask.height, mask.x:mask.x + mask.width] total = subframe[:, mask.mask].sum(axis=1) traces[i, frame_num:frame_num+block_size] = total / mask_areas[i] return traces
[docs]def calculate_roi_and_neuropil_traces(movie_h5, roi_mask_list, motion_border): """ get roi and neuropil masks """ # a combined binary mask for all ROIs (this is used to # subtracted ROIs from annuli mask_array = create_roi_mask_array(roi_mask_list) combined_mask = mask_array.max(axis=0) logging.info("%d total ROIs" % len(roi_mask_list)) # create neuropil masks for the central ROIs neuropil_masks = [] for m in roi_mask_list: nmask = create_neuropil_mask(m, motion_border, combined_mask, "neuropil for " + m.label) neuropil_masks.append(nmask) # calculate fluorescence traces for valid ROI and neuropil masks # create a combined list and calculate these together (this lets us # read the large image stack only once) combined_list = [] for m in roi_mask_list: combined_list.append(m) for n in neuropil_masks: combined_list.append(n) with h5py.File(movie_h5, "r") as movie_f: stack_frames = movie_f["data"] logging.info("Calculating %d traces (neuropil + ROI) over %d frames" % (len(combined_list), len(stack_frames))) traces = calculate_traces(stack_frames, combined_list) roi_traces = traces[:len(roi_mask_list)] neuropil_traces = traces[len(roi_mask_list):] return roi_traces, neuropil_traces
[docs]def create_roi_mask_array(rois): '''Create full image mask array from list of RoiMasks. Parameters ---------- rois: list<RoiMask> List of roi masks. Returns ------- np.ndarray: NxWxH array Boolean array of of len(rois) image masks. ''' if rois: height = rois[0].img_rows width = rois[0].img_cols masks = np.zeros((len(rois), height, width), dtype=np.uint8) for i, roi in enumerate(rois): masks[i, :, :] = roi.get_mask_plane() else: masks = None return masks