from __future__ import division
import logging
import os
import sys
from six import iteritems
import numpy as np
import SimpleITK as sitk
from skimage.draw import polygon
from allensdk.config.manifest import Manifest
if sys.version_info[0] > 2:
failed_import = (ImportError, ModuleNotFoundError)
else:
failed_import = (ImportError,)
# use np_sitk_convert or sitk_np_convert to access
# TODO: check if this already exists. If not: add more dtypes
# it does not
NUMPY_SITK_TYPE_LOOKUP = {np.dtype(np.float32): sitk.sitkFloat32}
SITK_NUMPY_TYPE_LOOKUP = {v: k for k, v in iteritems(NUMPY_SITK_TYPE_LOOKUP)}
# ITK/Numpy
[docs]def set_image_spacing(image, spacing, origin=True):
'''
'''
spacing = np.array(spacing)
image.SetSpacing(spacing.tolist())
if origin:
image.SetOrigin((0.5 * spacing).tolist())
[docs]def new_image(dims, spacing, dtype, origin=True):
'''
'''
if len(dims) == 2:
image = sitk.Image(dims[0], dims[1], dtype)
elif len(dims) == 3:
image = sitk.Image(dims[0], dims[1], dims[2], dtype)
set_image_spacing(image, spacing, origin)
return image
[docs]def image_from_array(array, spacing, origin=True):
'''
'''
image = sitk.GetImageFromArray(array)
set_image_spacing(image, spacing, origin)
return image
[docs]def np_sitk_convert(np_type):
'''
'''
return NUMPY_SITK_TYPE_LOOKUP[np_type]
[docs]def sitk_np_convert(sitk_type):
'''
'''
return SITK_NUMPY_TYPE_LOOKUP[sitk_type]
# Math
[docs]def compute_coarse_parameters(in_dims, in_spacing, out_spacing, reduce_level):
'''
'''
reduce_factor = pow(2, reduce_level)
fradius = np.divide(out_spacing, in_spacing) / 2.0 / reduce_factor
coarse_grid_radius = np.round(fradius)
coarse_grid_size = (coarse_grid_radius * 2 + 1) * reduce_factor
coarse_grid_spacing = np.multiply(in_spacing, coarse_grid_size)
coarse_grid_dims = np.ceil(
np.divide(in_dims, coarse_grid_size)
).astype(int)
return coarse_grid_dims, coarse_grid_spacing, coarse_grid_radius
[docs]def block_apply(in_image, out_shape, dtype, blocks, fn):
'''
'''
out_image = np.zeros(out_shape, dtype=dtype)
for ii, row_block in enumerate(blocks[0]):
for jj, col_block in enumerate(blocks[1]):
out_image[ii, jj] = fn(in_image[row_block[0]:row_block[1],
col_block[0]:col_block[1]])
return out_image
[docs]def grid_image_blocks(in_shape, in_spacing, out_spacing):
'''
'''
blocks = []
out_shape = []
for dim in range(len(in_shape)):
in_px_centers = np.arange(in_spacing[dim]*0.5,
in_shape[dim]*in_spacing[dim],
in_spacing[dim])
out_px_edges = np.arange(out_spacing[dim],
(in_shape[dim]-0.5)*in_spacing[dim],
out_spacing[dim])
dig = np.digitize(in_px_centers, out_px_edges)
inds = np.where(np.diff(dig) > 0)[0] + 1
inds = [0] + inds.tolist() + [in_shape[dim]]
dim_blocks = [
(int(inds[i]), int(inds[i+1])) for i in range(len(inds)-1)
]
out_shape.append(len(dim_blocks))
blocks.append(dim_blocks)
return blocks, out_shape
# Polygons
[docs]def rasterize_polygons(shape, scale, polys):
canvas = np.zeros(shape, dtype=np.uint8)
for points in polys:
rpts = np.array([
int(np.around(item[1] * scale[1])) for item in points
])
cpts = np.array([
int(np.around(item[0] * scale[0])) for item in points
])
poly = polygon(rpts, cpts)
canvas[poly] = 1
return canvas
# Transforms
[docs]def resample_into_volume(image, transform, z, vol, dtype=sitk.sitkFloat32):
'''
'''
if transform is None:
transform = sitk.Transform()
timage = sitk.Resample(image, transform, sitk.sitkLinear, 0.0, dtype)
tvol = sitk.JoinSeries(timage)
return sitk.Paste(vol, tvol, tvol.GetSize(), destinationIndex=[0, 0, z])
[docs]def resample_volume(volume, dims, spacing, interpolator=None, transform=None):
'''
'''
if transform is None:
transform = sitk.Transform()
if interpolator is None:
interpolator = sitk.sitkLinear
ref = new_image(dims, spacing, sitk.sitkFloat32, False)
return sitk.Resample(volume, ref, transform, interpolator)
[docs]def write_volume(volume,
name,
prefix=None,
specify_resolution=None,
extension='.nrrd',
paths=None):
if prefix is None:
path = name
else:
path = os.path.join(prefix, name)
if specify_resolution is not None:
if isinstance(specify_resolution, (float, np.floating)) and \
specify_resolution % 1.0 == 0:
specify_resolution = int(specify_resolution)
path = path + '_{0}'.format(specify_resolution)
path = path + extension
logging.info('writing {0} volume to {1}'.format(name, path))
Manifest.safe_make_parent_dirs(path)
volume.SetOrigin([0, 0, 0])
sitk.WriteImage(volume, str(path), True)
if paths is not None:
paths.append(path)
def __read_segmentation_image_with_kakadu(path):
if not os.path.exists(path):
raise OSError('file not found at {}'.format(path))
return jpeg_twok.read(path).T
def __read_intensity_image_with_kakadu(path, reduce_level, channel):
if not os.path.exists(path):
raise OSError('file not found at {}'.format(path))
return jpeg_twok.read(path, reduce_level, channel).T
def __read_segmentation_image_with_glymur(path):
return glymur.Jp2k(path)[:]
def __read_intensity_image_with_glymur(path):
return glymur.Jp2k(path)[:]
try:
# we use a proprietary library called kakadu internally
# (jpeg_twok is a python interface around that library)
# kakadu offers really good performance as well as support for
# advanced jp2 features
# however, since it is proprietary, we can't share it
# alongside the allensdk,
# so we default to glymur (a python openjpeg) for external users.
sys.path.append('/shared/bioapps/itk/itk_shared/jp2/build')
import jpeg_twok
read_segmentation_image = __read_segmentation_image_with_kakadu
read_intensity_image = __read_intensity_image_with_kakadu
except failed_import:
import glymur
read_segmentation_image = __read_segmentation_image_with_glymur
read_intensity_image = __read_intensity_image_with_glymur