Source code for ClearMap.ImageProcessing.Thresholding.Thresholding

"""
Thresholding
============

This module contains vairous thresholding routines, including 
hysteresis thresholding.

Module towards smart thresholding routines.
"""
__author__    = 'Christoph Kirst <ckirst@rockefeller.edu>'
__license__   = 'MIT License <http://www.opensource.org/licenses/mit-license.php>'
__copyright__ = 'Copyright (c) 2017 by Christoph Kirst, The Rockefeller University, New York City'

import os
import numpy as np

import ClearMap.IO.IO as io;

import pyximport;
#pyximport.install(setup_args={"include_dirs":np.get_include()}, reload_support=True)

old_get_distutils_extension = pyximport.pyximport.get_distutils_extension

[docs] def new_get_distutils_extension(modname, pyxfilename, language_level=None): extension_mod, setup_args = old_get_distutils_extension(modname, pyxfilename, language_level) extension_mod.language='c++' return extension_mod,setup_args
pyximport.pyximport.get_distutils_extension = new_get_distutils_extension pyximport.install(setup_args = {"include_dirs" : [np.get_include(), os.path.dirname(os.path.abspath(__file__))]}, reload_support=True) from . import ThresholdingCode as code ############################################################################### ### Hysteresis thresholding ###############################################################################
[docs] def threshold(source, sink = None, threshold = None, hysteresis_threshold = None, seeds = None, background = None): """Hysteresis thresholding. Arguments --------- source : array Input source. sink : array or None If None, a new array is allocated. threshold : float The threshold for the initial seeds. hysteresis_threshold : float or None The hysteresis threshold to extend the initial seeds. If None, no hysteresis thresholding is performed. seeds : array or None The seeds from which to start the hysteresis thresholds background : array or None Exclude this area from the hysteresis thresholding. Returns ------- sink : array Thresholded output. """ if threshold is None and seeds is None: raise ValueError('The threshold and seeds cannot both be None!'); if source is sink: raise NotImplementedError("Cannot perform operation in place.") if sink is None: sink = np.zeros(source.shape, dtype = 'int8', order = io.order(source)); source_flat = source.reshape(-1, order=io.order(source)); sink_flat = sink.reshape(-1, order=io.order(sink)); strides = np.array(io.element_strides(source)); if seeds is None: seeds = np.where(source_flat >= threshold)[0]; else: seeds_flat = seeds.reshape(-1, order=io.order(seeds)); seeds = np.where(seeds_flat)[0]; if hysteresis_threshold is not None: parameter_index = np.zeros(0, dtype = int); parameter_double = np.array([hysteresis_threshold], dtype = float); if background is not None: background_flat = background.reshape(-1, order=io.order(background)); background_flat = background_flat.view(dtype='uint8') code.threshold_to_background(source_flat, sink_flat, background_flat, strides, seeds, parameter_index, parameter_double); else: code.threshold(source_flat, sink_flat, strides, seeds, parameter_index, parameter_double); else: sink_flat[seeds] = 1; return sink;
############################################################################### ### Tests ############################################################################### def _test(): import numpy as np import ClearMap.Visualization.Plot3d as p3d import ClearMap.ImageProcessing.Thresholding.Thresholding as th r = np.arange(50); x,y,z = np.meshgrid(r,r,r); x,y,z = [i - 25 for i in (x,y,z)]; d = np.exp(-(x*x + y*y +z*z)/10.0**2) t = th.threshold(d, sink = None, threshold = 0.9, hysteresis_threshold=0.5); p3d.plot([[d,t]])