diff options
Diffstat (limited to 'silx/math/colormap.pyx')
| -rw-r--r-- | silx/math/colormap.pyx | 336 |
1 files changed, 244 insertions, 92 deletions
diff --git a/silx/math/colormap.pyx b/silx/math/colormap.pyx index c5d3e09..dde249d 100644 --- a/silx/math/colormap.pyx +++ b/silx/math/colormap.pyx @@ -1,7 +1,7 @@ # coding: utf-8 # /*########################################################################## # -# Copyright (c) 2018 European Synchrotron Radiation Facility +# Copyright (c) 2018-2020 European Synchrotron Radiation Facility # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal @@ -30,13 +30,16 @@ __license__ = "MIT" __date__ = "16/05/2018" +import os cimport cython from cython.parallel import prange cimport numpy as cnumpy -from libc.math cimport frexp, sqrt +from libc.math cimport frexp, sinh, sqrt from .math_compatibility cimport asinh, isnan, isfinite, lrint, INFINITY, NAN import logging +import numbers + import numpy __all__ = ['cmap'] @@ -44,6 +47,16 @@ __all__ = ['cmap'] _logger = logging.getLogger(__name__) +cdef int DEFAULT_NUM_THREADS +if hasattr(os, 'sched_getaffinity'): + DEFAULT_NUM_THREADS = min(4, len(os.sched_getaffinity(0))) +elif os.cpu_count() is not None: + DEFAULT_NUM_THREADS = min(4, os.cpu_count()) +else: # Fallback + DEFAULT_NUM_THREADS = 1 +# Number of threads to use for the computation (initialized to up to 4) + + # Supported data types ctypedef fused data_types: cnumpy.uint8_t @@ -84,60 +97,193 @@ ctypedef fused image_types: float -# Type of scale function -ctypedef double (*scale_function)(double) nogil +# Normalization +ctypedef double (*NormalizationFunction)(double) nogil + + +cdef class Normalization: + """Base class for colormap normalization""" + + def apply(self, data, double vmin, double vmax): + """Apply normalization. + + :param Union[float,numpy.ndarray] data: + :param float vmin: Lower bound of the range + :param float vmax: Upper bound of the range + :rtype: Union[float,numpy.ndarray] + """ + cdef int length + cdef double[:] result + + if isinstance(data, numbers.Real): + return self.apply_double(<double> data, vmin, vmax) + else: + data = numpy.array(data, copy=False) + length = <int> data.size + result = numpy.empty(length, dtype=numpy.float64) + data1d = numpy.ravel(data) + for index in range(length): + result[index] = self.apply_double( + <double> data1d[index], vmin, vmax) + return numpy.array(result).reshape(data.shape) + + def revert(self, data, double vmin, double vmax): + """Revert normalization. + + :param Union[float,numpy.ndarray] data: + :param float vmin: Lower bound of the range + :param float vmax: Upper bound of the range + :rtype: Union[float,numpy.ndarray] + """ + cdef int length + cdef double[:] result + + if isinstance(data, numbers.Real): + return self.revert_double(<double> data, vmin, vmax) + else: + data = numpy.array(data, copy=False) + length = <int> data.size + result = numpy.empty(length, dtype=numpy.float64) + data1d = numpy.ravel(data) + for index in range(length): + result[index] = self.revert_double( + <double> data1d[index], vmin, vmax) + return numpy.array(result).reshape(data.shape) + + cdef double apply_double(self, double value, double vmin, double vmax) nogil: + """Apply normalization to a floating point value + + Override in subclass + + :param float value: + :param float vmin: Lower bound of the range + :param float vmax: Upper bound of the range + """ + return value + + cdef double revert_double(self, double value, double vmin, double vmax) nogil: + """Apply inverse of normalization to a floating point value + + Override in subclass + + :param float value: + :param float vmin: Lower bound of the range + :param float vmax: Upper bound of the range + """ + return value + + +cdef class LinearNormalization(Normalization): + """Linear normalization""" + + cdef double apply_double(self, double value, double vmin, double vmax) nogil: + return value + + cdef double revert_double(self, double value, double vmin, double vmax) nogil: + return value + + +cdef class LogarithmicNormalization(Normalization): + """Logarithmic normalization using a fast log approximation""" + cdef: + readonly int lutsize + readonly double[::1] lut # LUT used for fast log approximation + + def __cinit__(self, int lutsize=4096): + # Initialize log approximation LUT + self.lutsize = lutsize + self.lut = numpy.log2( + numpy.linspace(0.5, 1., lutsize + 1, + endpoint=True).astype(numpy.float64)) + # index_lut can overflow of 1 + self.lut[lutsize] = self.lut[lutsize - 1] + + def __dealloc__(self): + self.lut = None + + @cython.wraparound(False) + @cython.boundscheck(False) + @cython.nonecheck(False) + @cython.cdivision(True) + cdef double apply_double(self, double value, double vmin, double vmax) nogil: + """Return log10(value) fast approximation based on LUT""" + cdef double result = NAN # if value < 0.0 or value == NAN + cdef int exponent, index_lut + cdef double mantissa # in [0.5, 1) unless value == 0 NaN or +/-inf + + if value <= 0.0 or not isfinite(value): + if value == 0.0: + result = - INFINITY + elif value > 0.0: # i.e., value = +INFINITY + result = value # i.e. +INFINITY + else: + mantissa = frexp(value, &exponent) + index_lut = lrint(self.lutsize * 2 * (mantissa - 0.5)) + # 1/log2(10) = 0.30102999566398114 + result = 0.30102999566398114 * (<double> exponent + + self.lut[index_lut]) + return result + + cdef double revert_double(self, double value, double vmin, double vmax) nogil: + return 10**value + + +cdef class ArcsinhNormalization(Normalization): + """Inverse hyperbolic sine normalization""" + + cdef double apply_double(self, double value, double vmin, double vmax) nogil: + return asinh(value) + + cdef double revert_double(self, double value, double vmin, double vmax) nogil: + return sinh(value) -# Normalization +cdef class SqrtNormalization(Normalization): + """Square root normalization""" -cdef double linear_scale(double value) nogil: - """No-Op scaling function""" - return value + cdef double apply_double(self, double value, double vmin, double vmax) nogil: + return sqrt(value) + cdef double revert_double(self, double value, double vmin, double vmax) nogil: + return value**2 -cdef double asinh_scale(double value) nogil: - """asinh scaling function - Wraps asinh as it is defined as a macro for Windows support. +cdef class PowerNormalization(Normalization): + """Gamma correction: + + Linear normalization to [0, 1] followed by power normalization. + + :param gamma: Gamma correction factor """ - return asinh(value) + cdef: + readonly double gamma -DEF LOG_LUT_SIZE = 4096 -cdef double[::1] _log_lut -"""LUT used for fast log approximation""" + def __cinit__(self, double gamma): + self.gamma = gamma -# Initialize log approximation LUT -_log_lut = numpy.log2( - numpy.linspace(0.5, 1., LOG_LUT_SIZE + 1, - endpoint=True).astype(numpy.float64)) -# index_lut can overflow of 1 -_log_lut[LOG_LUT_SIZE] = _log_lut[LOG_LUT_SIZE - 1] + def __init__(self, gamma): + # Needed for multiple inheritance to work + pass + cdef double apply_double(self, double value, double vmin, double vmax) nogil: + if vmin == vmax: + return 0. + elif value <= vmin: + return 0. + elif value >= vmax: + return 1. + else: + return ((value - vmin) / (vmax - vmin))**self.gamma -@cython.wraparound(False) -@cython.boundscheck(False) -@cython.nonecheck(False) -@cython.cdivision(True) -cdef double fast_log10(double value) nogil: - """Return log10(value) fast approximation based on LUT""" - cdef double result = NAN # if value < 0.0 or value == NAN - cdef int exponent, index_lut - cdef double mantissa # in [0.5, 1) unless value == 0 NaN or +/-inf - - if value <= 0.0 or not isfinite(value): - if value == 0.0: - result = - INFINITY - elif value > 0.0: # i.e., value = +INFINITY - result = value # i.e. +INFINITY - else: - mantissa = frexp(value, &exponent) - index_lut = lrint(LOG_LUT_SIZE * 2 * (mantissa - 0.5)) - # 1/log2(10) = 0.30102999566398114 - result = 0.30102999566398114 * (<double> exponent + - _log_lut[index_lut]) - return result + cdef double revert_double(self, double value, double vmin, double vmax) nogil: + if value <= 0.: + return vmin + elif value >= 1.: + return vmax + else: + return vmin + (vmax - vmin) * value**(1.0/self.gamma) # Colormap @@ -149,22 +295,22 @@ cdef double fast_log10(double value) nogil: cdef image_types[:, ::1] compute_cmap( default_types[:] data, image_types[:, ::1] colors, - double normalized_vmin, - double normalized_vmax, - image_types[::1] nan_color, - scale_function scale_func): + Normalization normalization, + double vmin, + double vmax, + image_types[::1] nan_color): """Apply colormap to data. :param data: Input data :param colors: Colors look-up-table - :param normalized_vmin: Normalized lower bound of the colormap range - :param normalized_vmax: Normalized upper bound of the colormap range + :param vmin: Lower bound of the colormap range + :param vmax: Upper bound of the colormap range :param nan_color: Color to use for NaN value - :param scale_func: The function to use to scale data + :param normalization: Normalization to apply :return: Data converted to colors """ cdef image_types[:, ::1] output - cdef double scale, value + cdef double scale, value, normalized_vmin, normalized_vmax cdef int length, nb_channels, nb_colors cdef int channel, index, lut_index @@ -175,14 +321,21 @@ cdef image_types[:, ::1] compute_cmap( output = numpy.empty((length, nb_channels), dtype=numpy.array(colors, copy=False).dtype) + normalized_vmin = normalization.apply_double(vmin, vmin, vmax) + normalized_vmax = normalization.apply_double(vmax, vmin, vmax) + + if not isfinite(normalized_vmin) or not isfinite(normalized_vmax): + raise ValueError('Colormap range is not valid') + if normalized_vmin == normalized_vmax: scale = 0. else: scale = nb_colors / (normalized_vmax - normalized_vmin) with nogil: - for index in prange(length): - value = scale_func(<double> data[index]) + for index in prange(length, num_threads=DEFAULT_NUM_THREADS): + value = normalization.apply_double( + <double> data[index], vmin, vmax) # Handle NaN if isnan(value): @@ -212,20 +365,20 @@ cdef image_types[:, ::1] compute_cmap( cdef image_types[:, ::1] compute_cmap_with_lut( lut_types[:] data, image_types[:, ::1] colors, - double normalized_vmin, - double normalized_vmax, - image_types[::1] nan_color, - scale_function scale_func): + Normalization normalization, + double vmin, + double vmax, + image_types[::1] nan_color): """Convert data to colors using look-up table to speed the process. Only supports data of types: uint8, uint16, int8, int16. :param data: Input data :param colors: Colors look-up-table - :param normalized_vmin: Normalized lower bound of the colormap range - :param normalized_vmax: Normalized upper bound of the colormap range + :param vmin: Lower bound of the colormap range + :param vmax: Upper bound of the colormap range :param nan_color: Color to use for NaN values - :param scale_func: The function to use for scaling data + :param normalization: Normalization to apply :return: The generated image """ cdef image_types[:, ::1] output @@ -255,14 +408,13 @@ cdef image_types[:, ::1] compute_cmap_with_lut( values = numpy.arange(type_min, type_max + 1, dtype=numpy.float64) lut = compute_cmap( - values, colors, normalized_vmin, normalized_vmax, - nan_color, scale_func) + values, colors, normalization, vmin, vmax, nan_color) output = numpy.empty((length, nb_channels), dtype=colors_dtype) with nogil: # Apply LUT - for index in prange(length): + for index in prange(length, num_threads=DEFAULT_NUM_THREADS): lut_index = data[index] - type_min for channel in range(nb_channels): output[index, channel] = lut[lut_index, channel] @@ -270,13 +422,22 @@ cdef image_types[:, ::1] compute_cmap_with_lut( return output +# Normalizations without parameters +_BASIC_NORMALIZATIONS = { + 'linear': LinearNormalization(), + 'log': LogarithmicNormalization(), + 'arcsinh': ArcsinhNormalization(), + 'sqrt': SqrtNormalization(), + } + + @cython.wraparound(False) @cython.boundscheck(False) @cython.nonecheck(False) @cython.cdivision(True) def _cmap(data_types[:] data, image_types[:, ::1] colors, - str normalization, + Normalization normalization, double vmin, double vmax, image_types[::1] nan_color): @@ -286,42 +447,20 @@ def _cmap(data_types[:] data, :param data: Input data :param colors: Colors look-up-table - :param normalization: Kind of scaling to apply on data + :param normalization: Normalization object to apply :param vmin: Lower bound of the colormap range :param vmax: Upper bound of the colormap range :param nan_color: Color to use for NaN value. :return: The generated image """ - cdef double normalized_vmin, normalized_vmax - cdef scale_function scale_func - - if normalization == 'linear': - scale_func = linear_scale - elif normalization == 'log': - scale_func = fast_log10 - elif normalization == 'arcsinh': - scale_func = asinh_scale - elif normalization == 'sqrt': - scale_func = sqrt - else: - raise ValueError('Unsupported normalization %s' % normalization) - - normalized_vmin = scale_func(vmin) - normalized_vmax = scale_func(vmax) - - if not isfinite(normalized_vmin) or not isfinite(normalized_vmax): - raise ValueError('Colormap range is not valid') - # Proxy for calling the right implementation depending on data type if data_types in lut_types: # Use LUT implementation output = compute_cmap_with_lut( - data, colors, normalized_vmin, normalized_vmax, - nan_color, scale_func) + data, colors, normalization, vmin, vmax, nan_color) elif data_types in default_types: # Use default implementation output = compute_cmap( - data, colors, normalized_vmin, normalized_vmax, - nan_color, scale_func) + data, colors, normalization, vmin, vmax, nan_color) else: raise ValueError('Unsupported data type') @@ -342,12 +481,14 @@ def cmap(data, It MUST be of type uint8 or float32 :param vmin: Data value to map to the beginning of colormap. :param vmax: Data value to map to the end of the colormap. - :param str normalization: The normalization to apply: + :param Union[str,Normalization] normalization: + Either a :class:`Normalization` instance or a str in: - 'linear' (default) - 'log' - 'arcsinh' - 'sqrt' + - 'gamma' :param nan_color: Color to use for NaN value. Default: A color with all channels set to 0 @@ -357,6 +498,7 @@ def cmap(data, :rtype: numpy.ndarray """ cdef int nb_channels + cdef Normalization norm # Make data a numpy array of native endian type (no need for contiguity) data = numpy.array(data, copy=False) @@ -371,6 +513,14 @@ def cmap(data, colors = numpy.ascontiguousarray(colors, dtype=colors.dtype.newbyteorder('N')) + # Make normalization a Normalization object + if isinstance(normalization, str): + norm = _BASIC_NORMALIZATIONS.get(normalization, None) + if norm is None: + raise ValueError('Unsupported normalization %s' % normalization) + else: + norm = normalization + # Check nan_color if nan_color is None: nan_color = numpy.zeros((nb_channels,), dtype=colors.dtype) @@ -382,8 +532,10 @@ def cmap(data, image = _cmap( data.reshape(-1), colors.reshape(-1, nb_channels), - str(normalization), - vmin, vmax, nan_color) + norm, + vmin, + vmax, + nan_color) image.shape = data.shape + (nb_channels,) return image |