Source code for allensdk.core.reference_space

# 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.
#
from __future__ import division, print_function, absolute_import
from collections import defaultdict
import operator as op
import functools
import os

from scipy.misc import imresize
from scipy.ndimage.interpolation import zoom
import numpy as np
import nrrd

from allensdk.core.structure_tree import StructureTree


[docs]class ReferenceSpace(object): @property def direct_voxel_map(self): if not hasattr(self, '_direct_voxel_map'): self.direct_voxel_counts() return self._direct_voxel_map @direct_voxel_map.setter def direct_voxel_map(self, data): self._direct_voxel_map = data @property def total_voxel_map(self): if not hasattr(self, '_total_voxel_map'): self.total_voxel_counts() return self._total_voxel_map @total_voxel_map.setter def total_voxel_map(self, data): self._total_voxel_map = data def __init__(self, structure_tree, annotation, resolution): '''Handles brain structures in a 3d reference space Parameters ---------- structure_tree : StructureTree Defines the heirarchy and properties of the brain structures. annotation : numpy ndarray 3d volume whose elements are structure ids. resolution : length-3 tuple of numeric Resolution of annotation voxels along each dimension. ''' self.structure_tree = structure_tree self.resolution = resolution self.annotation = np.ascontiguousarray(annotation)
[docs] def direct_voxel_counts(self): '''Determines the number of voxels directly assigned to one or more structures. Returns ------- dict : Keys are structure ids, values are the number of voxels directly assigned to those structures. ''' uniques = np.unique(self.annotation, return_counts=True) found = {k: v for k, v in zip(*uniques) if k != 0} self._direct_voxel_map = {k: (found[k] if k in found else 0) for k in self.structure_tree.node_ids()}
[docs] def total_voxel_counts(self): '''Determines the number of voxels assigned to a structure or its descendants Returns ------- dict : Keys are structure ids, values are the number of voxels assigned to structures' descendants. ''' self._total_voxel_map = {} for stid in self.structure_tree.node_ids(): desc_ids = self.structure_tree.descendant_ids([stid])[0] self._total_voxel_map[stid] = sum([self.direct_voxel_map[dscid] for dscid in desc_ids])
[docs] def remove_unassigned(self, update_self=True): '''Obtains a structure tree consisting only of structures that have at least one voxel in the annotation. Parameters ---------- update_self : bool, optional If True, the contained structure tree will be replaced, Returns ------- list of dict : elements are filtered structures ''' structures = self.structure_tree.filter_nodes( lambda x: self.total_voxel_map[x['id']] > 0) if update_self: self.structure_tree = StructureTree(structures) return structures
[docs] def make_structure_mask(self, structure_ids, direct_only=False): '''Return an indicator array for one or more structures Parameters ---------- structure_ids : list of int Make a mask that indicates the union of these structures' voxels direct_only : bool, optional If True, only include voxels directly assigned to a structure in the mask. Otherwise include voxels assigned to descendants. Returns ------- numpy ndarray : Same shape as annotation. 1 inside mask, 0 outside. ''' if direct_only: mask = np.zeros(self.annotation.shape, dtype=np.uint8, order='C') for stid in structure_ids: if self.direct_voxel_map[stid] == 0: continue mask[self.annotation == stid] = True return mask else: structure_ids = self.structure_tree.descendant_ids(structure_ids) structure_ids = set(functools.reduce(op.add, structure_ids)) return self.make_structure_mask(structure_ids, direct_only=True)
[docs] def many_structure_masks(self, structure_ids, output_cb=None, direct_only=False): '''Build one or more structure masks and do something with them Parameters ---------- structure_ids : list of int Specify structures to be masked output_cb : function, optional Must have the following signature: output_cb(structure_id, fn). On each requested id, fn will be curried to make a mask for that id. Defaults to returning the structure id and mask. direct_only : bool, optional If True, only include voxels directly assigned to a structure in the mask. Otherwise include voxels assigned to descendants. Yields ------- Return values of output_cb called on each structure_id, structure_mask pair. Notes ----- output_cb is called on every yield, so any side-effects (such as writing to a file) will be carried out regardless of what you do with the return values. You do actually have to iterate through the output, though. ''' if output_cb is None: output_cb = ReferenceSpace.return_mask_cb for stid in structure_ids: yield output_cb(stid, functools.partial(self.make_structure_mask, [stid], direct_only))
[docs] def check_coverage(self, structure_ids, domain_mask): '''Determines whether a spatial domain is completely covered by structures in a set. Parameters ---------- structure_ids : list of int Specifies the set of structures to check. domain_mask : numpy ndarray Same shape as annotation. 1 inside the mask, 0 out. Specifies spatial domain. Returns ------- numpy ndarray : 1 where voxels are missing from the candidate, 0 where the candidate exceeds the domain ''' candidate_mask = self.make_structure_mask(structure_ids) return domain_mask - candidate_mask
[docs] def validate_structures(self, structure_ids, domain_mask): '''Determines whether a set of structures produces an exact and nonoverlapping tiling of a spatial domain Parameters ---------- structure_ids : list of int Specifies the set of structures to check. domain_mask : numpy ndarray Same shape as annotation. 1 inside the mask, 0 out. Specifies spatial domain. Returns ------- set : Ids of structures that are the ancestors of other structures in the supplied set. numpy ndarray : Indicator for missing voxels. ''' return [self.structure_tree.has_overlaps(structure_ids), self.check_coverage(structure_ids, domain_mask)]
[docs] def downsample(self, target_resolution): '''Obtain a smaller reference space by downsampling Parameters ---------- target_resolution : tuple of numeric Resolution in microns of the output space. interpolator : string Method used to interpolate the volume. Currently only 'nearest' is supported Returns ------- ReferenceSpace : A new ReferenceSpace with the same structure tree and a downsampled annotation. ''' factors = [ float(ii / jj) for ii, jj in zip(self.resolution, target_resolution)] target = zoom(self.annotation, factors, order=0) return ReferenceSpace(self.structure_tree, target, target_resolution)
[docs] def get_slice_image(self, axis, position, cmap=None): '''Produce a AxBx3 RGB image from a slice in the annotation Parameters ---------- axis : int Along which to slice the annotation volume. 0 is coronal, 1 is horizontal, and 2 is sagittal. position : int In microns. Take the slice from this far along the specified axis. cmap : dict, optional Keys are structure ids, values are rgb triplets. Defaults to structure rgb_triplets. Returns ------- np.ndarray : RGB image array. Notes ----- If you assign a custom colormap, make sure that you take care of the background in addition to the structures. ''' if cmap is None: cmap = self.structure_tree.get_colormap() cmap[0] = [0, 0, 0] position = int(np.around(position / self.resolution[axis])) image = np.squeeze(self.annotation.take([position], axis=axis)) return np.reshape([cmap[point] for point in image.flat], list(image.shape) + [3]).astype(np.uint8)
[docs] @staticmethod def return_mask_cb(structure_id, fn): '''A basic callback for many_structure_masks ''' return structure_id, fn()
[docs] @staticmethod def check_and_write(base_dir, structure_id, fn): '''A many_structure_masks callback that writes the mask to a nrrd file if the file does not already exist. ''' mask_path = os.path.join(base_dir, 'structure_{0}.nrrd'.format(structure_id)) if not os.path.exists(mask_path): nrrd.write(mask_path, fn()) return structure_id