# -*- coding: utf-8 -*-
"""
LocalStatistics
===============
Module provides functions to calculate local data and statistics of an image
and apply a function to those. It is useful for local and adaptive image
processing of large 3d images.
Note
----
The module provides ways to speed up the local statistics by only sampling
on a sub-grid of the image and resample the result to the full image shape.
"""
__author__ = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__ = 'GPLv3 - GNU General Pulic License v3 (see LICENSE.txt)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__ = 'http://idisco.info'
__download__ = 'http://www.github.com/ChristophKirst/ClearMap2'
import numpy as np
import scipy.ndimage as ndi
###############################################################################
### Local image processing
###############################################################################
[docs]
def apply_local_function(source, function, selem=(50,50), spacing=None, step=None, interpolate=2,
mask=None, fshape=None, dtype=None, return_centers=False):
"""Calculate local histograms on a sub-grid, apply a scalar valued function and resample to original image shape.
Arguments
---------
source : array
The source to process.
function : function
Function to apply to the linear array of the local source data.
If the function does not return a scalar, fshape has to be given.
selem : tuple or array or None
The structural element to use to extract the local image data.
If tuple, use a rectangular region of this shape. If array, the array
is assumed to be bool and acts as a local mask around the center point.
spacing : tuple or None
The spacing between sample points. If None, use shape of structural_element.
step : tuple of int or None
If tuple, subsample the local region by these step. Note that the
structural_element is applied after this subsampling.
interpolate : int or None
If int, resample the result back to the original source shape using this
order of interpolation. If None, return the results on the sub-grid.
mask : array or None
Optional mask to use.
fshape : tuple or None
If tuple, this is the shape of the function output.
If None assumed to be (1,).
dtype : dtype or None
Optional data type for the result.
return_centers : bool
If True, additionaly return the centers of the sampling.
Returns
-------
local : array
The reuslt of applying the function to the local samples.
cetners : array
Optional cttners of the sampling.
"""
if spacing is None:
spacing = selem;
shape = source.shape;
ndim = len(shape);
if step is None:
step = (None,) * ndim
if len(spacing) != ndim or len(step) != ndim:
raise ValueError('Dimension mismatch in the parameters!')
#histogram centers
n_centers = tuple(s//h for s,h in zip(shape, spacing))
left = tuple((s - (n-1) * h)//2 for s,n,h in zip(shape, n_centers, spacing));
#center points
centers = np.array(np.meshgrid(*[range(l, s, h) for l,s,h in zip(left, shape, spacing)], indexing = 'ij'));
#centers = np.reshape(np.moveaxis(centers, 0, -1),(-1,len(shape)));
centers = np.moveaxis(centers, 0, -1)
#create result
rshape = (1,) if fshape is None else fshape;
rdtype = source.dtype if dtype is None else dtype;
results = np.zeros(n_centers + rshape, dtype = rdtype);
#calculate function
centers_flat = np.reshape(centers, (-1,ndim));
results_flat = np.reshape(results, (-1,) + rshape);
#structuring element
if isinstance(selem, np.ndarray):
selem_shape = selem.shape
else:
selem_shape = selem;
selem = None
hshape_left = tuple(h//2 for h in selem_shape);
hshape_right = tuple(h - l for h,l in zip(selem_shape, hshape_left));
for result, center in zip(results_flat,centers_flat):
sl = tuple(slice(max(0,c-l), min(c+r,s), d) for c,l,r,s,d in zip(center, hshape_left, hshape_right, shape, step));
if selem is None:
if mask is not None:
data = source[sl][mask[sl]];
else:
data = source[sl].flatten();
else:
slm = tuple(slice(None if c-l >= 0 else min(l-c,m), None if c+r <= s else min(m - (c + r - s), m), d) for c,l,r,s,d,m in zip(center, hshape_left, hshape_right, shape, step, selem_shape));
data = source[sl];
if mask is not None:
data = data[np.logical_and(mask[sl], selem[slm])]
else:
data = data[selem[slm]];
#print result.shape, data.shape, function(data)
result[:] = function(data);
#resample
if interpolate:
res_shape = results.shape[:len(shape)];
zoom = tuple(float(s) / float(r) for s,r in zip(shape, res_shape));
results_flat = np.reshape(results, res_shape + (-1,));
results_flat = np.moveaxis(results_flat, -1, 0);
full = np.zeros(shape + rshape, dtype = results.dtype);
full_flat = np.reshape(full, shape + (-1,));
full_flat = np.moveaxis(full_flat, -1, 0);
#print results_flat.shape, full_flat.shape
for r,f in zip(results_flat, full_flat):
f[:] = ndi.zoom(r, zoom=zoom, order = interpolate);
results = full;
if fshape is None:
results.shape = results.shape[:-1];
if return_centers:
return results, centers
else:
return results
[docs]
def local_histogram(source, max_bin = 2**12, selem = (50,50), spacing = None, step = None, interpolate = None, mask = None, dtype = None, return_centers = False):
"""Calculate local histograms on a sub-grid.
Arguments
---------
source : array
The source to process.
selem : tuple or array or None
The structural element to use to extract the local image data.
If tuple, use a rectangular region of this shape. If array, the array
is assumed to be bool and acts as a local mask around the center point.
spacing : tuple or None
The spacing between sample points. If None, use shape of structural_element.
step : tuple of int or None
If tuple, subsample the local region by these step. Note that the
structural_element is applied after this subsampling.
interpolate : int or None
If int, resample the result back to the original source shape using this
order of interpolation. If None, return the results on the sub-grid.
mask : array or None
Optional mask to use.
max_bin : int
Maximal bin value to account for.
return_centers : bool
If True, additionaly return the centers of the sampling.
Returns
-------
histograms : array
The local histograms.
cetners : array
Optional centers of the sampling.
Note
----
For speed, this function works only for uint sources as the histogram is
calculated directly via the source values. The source values should be
smaller than max_bin.
"""
def _hist(data):
data, counts = np.unique(data, return_counts=True);
histogram = np.zeros(max_bin, dtype=int);
histogram[data] = counts;
return histogram;
return apply_local_function(source, selem=selem, spacing=spacing, step=step, interpolate=interpolate, mask=mask, dtype=dtype, return_centers=return_centers,
function=_hist, fshape = (max_bin,));
[docs]
def local_percentile(source, percentile, selem = (50,50), spacing = None, step = None, interpolate = 1, mask = None, dtype = None, return_centers = False):
"""Calculate local percentile.
Arguments
---------
source : array
The source to process.
percentile : float or array
The percentile(s) to estimate locally.
selem : tuple or array or None
The structural element to use to extract the local image data.
If tuple, use a rectangular region of this shape. If array, the array
is assumed to be bool and acts as a local mask around the center point.
spacing : tuple or None
The spacing between sample points. If None, use shape of structural_element.
step : tuple or int or None
If tuple, subsample the local region by these steps. Note that the
structural_element is applied after this subsampling.
interpolate : int or None
If int, resample the result back to the original source shape using this
order of interpolation. If None, return the results on the sub-grid.
mask : array or None
Optional mask to use.
return_centers : bool
If True, additionally return the centers of the sampling.
Returns
-------
percentiles : array
The local percentiles.
centers : array
Optional centers of the sampling.
"""
if isinstance(percentile, (tuple, list)):
percentile = np.array([100*p for p in percentile]);
fshape = (len(percentile),)
def _percentile(data):
if len(data) == 0:
return np.array((0,) * len(percentile));
return np.percentile(data, percentile, axis = None);
else:
percentile = 100 * percentile;
fshape = None;
def _percentile(data):
if len(data) == 0:
return 0;
return np.percentile(data, percentile, axis = None);
return apply_local_function(source, selem=selem, spacing=spacing, step=step, interpolate=interpolate, mask=mask, dtype=dtype, return_centers=return_centers,
function=_percentile, fshape=fshape);
###############################################################################
### Tests
###############################################################################
def _test():
"""Tests."""
import numpy as np
import ClearMap.Visualization.Plot3d as p3d
import ClearMap.ImageProcessing.LocalStatistics as ls
from importlib import reload
reload(ls)
source = np.random.rand(100,200,150) + np.arange(100)[:,None,None];
p = ls.local_percentile(source, percentile=0.5, selem=(30,30,30), interpolate=1);
p3d.plot([source, p])
#def apply_histogram_function(histograms, function, shape = None, interpolation = 2):
#
# hist_shape = histograms.shape[:-1];
# max_bin = histograms.shape[-1];
# result = np.zeros(np.prod(hist_shape), dtype = float);
# histograms_flat = np.reshape(histograms, (-1,max_bin));
# for i,h in enumerate(histograms_flat):
# result[i] = function(h);
# result.shape = hist_shape;
#
# if shape is not None:
# zoom = tuple(float(s) / float(h) for s,h in zip(shape, hist_shape))
# result = ndi.zoom(result, zoom=zoom, order = interpolation);
#
# return result;
