from __future__ import division
import logging
import os
import sys
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 NUMPY_SITK_TYPE_LOOKUP.items()}
# 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