Source code for ClearMap.ParallelProcessing.DataProcessing.MeasurePointList

# -*- coding: utf-8 -*-
"""
MeasurePointList
================

Measurements on a subset of points in large arrays.

Paralllel measuremnts at specified points of the data only.
Useful to speed up processing in large arrays and only a smaller number 
of measurement points.

See also
--------
:mod:`ClearMap.Analysis.Measurements`.
"""
__author__    = 'Christoph Kirst <christoph.kirst.ck@gmail.com>'
__license__   = 'GPLv3 - GNU General Pulic License v3 (see LICENSE)'
__copyright__ = 'Copyright © 2020 by Christoph Kirst'
__webpage__   = 'http://idisco.info'
__download__  = 'http://www.github.com/ChristophKirst/ClearMap2'


import numpy as np;

import pyximport; 
pyximport.install(setup_args={"include_dirs": [np.get_include()]}, reload_support=True)
 
import ClearMap.ParallelProcessing.DataProcessing.ArrayProcessing as ap

import ClearMap.ParallelProcessing.DataProcessing.MeasurePointListCode as code


###############################################################################
### Measure extrema
###############################################################################

[docs] def measure_max(source, points, search, max_search_indices, sink = None, processes = None, verbose = False): """Find local maximum in a large array for a list of center points. Arguments --------- source : array Data source. points : array List of linear indices of center points. search : array List of linear indices to add to the center index defining the local search area. max_search_indices : array The maximal index in the search array for each point to use. sink : array or None Optional sink for result indices. processes : int or None Number of processes to use. verbose : bool If True, print progress info. Returns ------- sink : array Linear array with length of points containing the local maxima. """ processes, timer = ap.initialize_processing(processes=processes, verbose=verbose, function='measure_max'); source, source_buffer, source_shape, source_strides = ap.initialize_source(source, as_1d=True, return_shape=True, return_strides=True); sink, sink_buffer = ap.initialize_sink(sink=sink, shape=points.shape[0], dtype=source.dtype); if sink_buffer.shape[0] != points.shape[0]: raise RuntimeError('Sink has invalid size %r not %r' % (sink.shape, points.shape)); points, points_buffer = ap.initialize_source(points); if points_buffer.ndim == 1: points_buffer = points_buffer[:,None]; #print(source1d.shape, shape, strides, points.shape, search.shape, max_search_indices.shape, sink.shape) code.measure_max(source_buffer, source_shape, source_strides, points_buffer, search, max_search_indices, sink_buffer, processes); ap.finalize_processing(verbose=verbose, function='measure_max', timer=timer); return sink;
[docs] def measure_min(source, points, search, max_search_indices, sink = None, processes = None, verbose = False): """Find local minimum in a large array for a list of center points. Arguments --------- source : array Data source. points : array List of linear indices of center points. search : array List of linear indices to add to the center index defining the local search area. max_search_indices : array The maximal index in the search array for each point to use. sink : array or None Optional sink for result indices. processes : int or None Number of processes to use. verbose : bool If True, print progress info. Returns ------- sink : array Linear array with length of points containing the local minima. """ processes, timer = ap.initialize_processing(processes=processes, verbose=verbose, function='measure_max'); source, source_buffer, source_shape, source_strides = ap.initialize_source(source, as_1d=True, return_shape=True, return_strides=True); sink, sink_buffer = ap.initialize_sink(sink=sink, shape=points.shape[0], dtype=source.dtype); if sink.shape[0] != points.shape[0]: raise RuntimeError('Sink has invalid size %r not %r' % (sink.shape, points.shape)); points, points_buffer = ap.initialize_source(points); if points_buffer.ndim == 1: points_buffer = points_buffer[:,None]; code.measure_min(source_buffer, source_shape, source_strides, points_buffer, search, max_search_indices, sink_buffer, processes); ap.finalize_processing(verbose=verbose, function='measure_max', timer=timer); return sink;
[docs] def measure_mean(source, points, search, max_search_indices, sink = None, processes = None, verbose = False): """Find local mean in a large array for a list of center points. Arguments --------- source : array Data source. points : array List of linear indices of center points. search : array List of linear indices to add to the center index defining the local search area. max_search_indices : array The maximal index in the search array for each point to use. sink : array or None Optional sink for result indices. processes : int or None Number of processes to use. verbose : bool If True, print progress info. Returns ------- sink : array Linear array with length of points containing the local mean. """ processes, timer = ap.initialize_processing(processes=processes, verbose=verbose, function='measure_mean'); source, source_buffer, source_shape, source_strides = ap.initialize_source(source, as_1d=True, return_shape=True, return_strides=True); sink, sink_buffer = ap.initialize_sink(sink=sink, shape=points.shape[0], dtype=float); points_buffer = ap.initialize_source(points); if points_buffer.ndim == 1: points_buffer = points_buffer[:,None]; if sink.shape[0] != points_buffer.shape[0]: raise RuntimeError('Sink has invalid size %r not %r' % (sink.shape, points_buffer.shape)); code.measure_mean(source_buffer, source_shape, source_strides, points_buffer, search, max_search_indices, sink_buffer, processes); ap.finalize_processing(verbose=verbose, function='measure_mean', timer=timer); return sink;
[docs] def measure_sum(source, points, search, max_search_indices, sink = None, processes = None, verbose = False): """Find local mean in a large array for a list of center points. Arguments --------- source : array Data source. points : array List of linear indices of center points. search : array List of linear indices to add to the center index defining the local search area. max_search_indices : array The maximal index in the search array for each point to use. sink : array or None Optional sink for result indices. processes : int or None Number of processes to use. verbose : bool If True, print progress info. Returns ------- sink : array Linear array with length of points containing the local mean. """ processes, timer = ap.initialize_processing(processes=processes, verbose=verbose, function='measure_mean'); source, source_buffer, source_shape, source_strides = ap.initialize_source(source, as_1d=True, return_shape=True, return_strides=True); sink, sink_buffer = ap.initialize_sink(sink=sink, shape=points.shape[0], dtype=float); points_buffer = ap.initialize_source(points); if points_buffer.ndim == 1: points_buffer = points_buffer[:,None]; if sink.shape[0] != points_buffer.shape[0]: raise RuntimeError('Sink has invalid size %r not %r' % (sink.shape, points_buffer.shape)); code.measure_sum(source_buffer, source_shape, source_strides, points_buffer, search, max_search_indices, sink_buffer, processes); ap.finalize_processing(verbose=verbose, function='measure_mean', timer=timer); return sink;
############################################################################### ### Find in local neighbourhood ###############################################################################
[docs] def find_smaller_than_value(source, points, search, value, sink = None, processes = None, verbose = False): """Find index in local search indices with a voxel with value smaller than a specified value for a list of points. Arguments --------- source : array Data source. points : array List of linear indices of center points. search : array List of linear indices to add to the center index defining the local search area. value : float Search for first voxel in local area with value smaller than this value. sink : array or None Optional sink for result indices. processes : int or None Number of processes to use. verbose : bool If True, print progress info. Returns ------- sink : array Linear array with length of points containing the first search index with voxel below value. """ processes, timer = ap.initialize_processing(processes=processes, verbose=verbose, function='find_smaller_than_value'); source, source_buffer, source_shape, source_strides = ap.initialize_source(source, as_1d=True, return_shape=True, return_strides=True); sink, sink_buffer = ap.initialize_sink(sink=sink, shape=points.shape[0], dtype=int); if sink.shape[0] != points.shape[0]: raise RuntimeError('Sink has invalid size %r not %r' % (sink.shape, points.shape)); points, points_buffer = ap.initialize_source(points); if points_buffer.ndim == 1: points_buffer = points_buffer[:,None]; code.find_smaller_than_value(source_buffer, source_shape, source_strides, points_buffer, search, value, sink_buffer, processes); ap.finalize_processing(verbose=verbose, function='find_smaller_than_value', timer=timer); return sink;
[docs] def find_smaller_than_fraction(source, points, search, fraction, sink = None, processes = None, verbose = False): """Find index in local search indices with a voxel with value smaller than a fraction of the value of the center voxel for a list of points. Arguments --------- source : array Data source. points : array List of linear indices of center points. search : array List of linear indices to add to the center index defining the local search area. fraction : float Search for first voxel in local area with value smaller than this fraction of the center value. sink : array or None Optional sink for result indices. processes : int or None Number of processes to use. verbose : bool If True, print progress info. Returns ------- sink : array Linear array with length of points containing the first search index with voxel below the fraction of the center value. """ processes, timer = ap.initialize_processing(processes=processes, verbose=verbose, function='find_smaller_than_fraction'); source, source_buffer, source_shape, source_strides = ap.initialize_source(source, as_1d=True, return_shape=True, return_strides=True); sink, sink_buffer = ap.initialize_sink(sink=sink, shape=points.shape[0], dtype=int); if sink.shape[0] != points.shape[0]: raise RuntimeError('Sink has invalid size %r not %r' % (sink.shape, points.shape)); points, points_buffer = ap.initialize_source(points); if points.ndim == 1: points = points[:,None]; code.find_smaller_than_fraction(source_buffer, source_shape, source_strides, points_buffer, search, fraction, sink_buffer, processes); ap.finalize_processing(verbose=verbose, function='find_smaller_than_fraction', timer=timer); return sink;
[docs] def find_smaller_than_values(source, points, search, values, sink = None, processes = None, verbose = False): """Find index in local search indices with a voxel with value smaller than a fraction of the value of the center voxel for a list of points. Arguments --------- source : array Data source. points : array List of linear indices of center points. search : array List of linear indices to add to the center index defining the local search area. fraction : float Search for first voxel in local area with value smaller than this fraction of the center value. sink : array or None Optional sink for result indices. processes : int or None Number of processes to use. verbose : bool If True, print progress info. Returns ------- sink : array Linear array with length of points containing the first search index with voxel below the fraction of the center value. """ processes, timer = ap.initialize_processing(processes=processes, verbose=verbose, function='find_smaller_than_values'); source, source_buffer, source_shape, source_strides = ap.initialize_source(source, as_1d=True, return_shape=True, return_strides=True); sink, sink_buffer = ap.initialize_sink(sink=sink, shape=points.shape, dtype=int); if sink.shape[0] != points.shape[0]: raise RuntimeError('Sink has invalid size %r not %r' % (sink.shape, points.shape)); points, points_buffer = ap.initialize_source(points); if points_buffer.ndim == 1: points_buffer = points_buffer[:,None]; code.find_smaller_than_fraction(source_buffer, source_shape, source_strides, points_buffer, search, values, sink_buffer, processes); ap.finalize_processing(verbose=verbose, function='find_smaller_than_values', timer=timer); return sink;
############################################################################### ### Tests ###############################################################################
[docs] def test(): import numpy as np import ClearMap.Analysis.Measurements.MeasureRadius as mr; import ClearMap.ParallelProcessing.DataProcessing.MeasurePointList as mpl from importlib import reload reload(mpl); reload(mpl.code); # find_smaller_than_value goal = (50,55,42); d = np.ones((100,100,100)); d[goal] = 0.3; strides = np.array(d.strides) / np.array(d.itemsize); search, dist = mr.indices_from_center(radius = 15, strides = strides); indices = np.array([np.ravel_multi_index((50,50,50), d.shape)]); result = mpl.find_smaller_than_value(d, indices, search, max_value = 0.5, out = None, processes = None); coords = np.unravel_index(indices + search[result], d.shape) coords = np.array(coords).reshape(-1); assert np.all(coords == goal) # find_smaller_than_fraction result = mpl.find_smaller_than_fraction(d, indices, search, fraction = 0.5, out = None, processes = None); coords = np.unravel_index(indices + search[result], d.shape) coords = np.array(coords).reshape(-1); assert np.all(coords == goal)