1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
|
# coding: utf-8
# /*##########################################################################
#
# Copyright (c) 2018-2021 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
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ############################################################################*/
"""This module provides :func:`cmap` which applies a colormap to a dataset.
"""
__authors__ = ["T. Vincent"]
__license__ = "MIT"
__date__ = "16/05/2018"
import os
cimport cython
from cython.parallel import prange
cimport numpy as cnumpy
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']
_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)
cdef int USE_OPENMP_THRESHOLD = 1000
"""OpenMP is not used for arrays with less elements than this threshold"""
# Supported data types
ctypedef fused data_types:
cnumpy.uint8_t
cnumpy.int8_t
cnumpy.uint16_t
cnumpy.int16_t
cnumpy.uint32_t
cnumpy.int32_t
cnumpy.uint64_t
cnumpy.int64_t
float
double
long double
# Data types using a LUT to apply the colormap
ctypedef fused lut_types:
cnumpy.uint8_t
cnumpy.int8_t
cnumpy.uint16_t
cnumpy.int16_t
# Data types using default colormap implementation
ctypedef fused default_types:
cnumpy.uint32_t
cnumpy.int32_t
cnumpy.uint64_t
cnumpy.int64_t
float
double
long double
# Supported colors/output types
ctypedef fused image_types:
cnumpy.uint8_t
float
# 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)
cdef class SqrtNormalization(Normalization):
"""Square root normalization"""
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 class PowerNormalization(Normalization):
"""Gamma correction:
Linear normalization to [0, 1] followed by power normalization.
:param gamma: Gamma correction factor
"""
cdef:
readonly double gamma
def __cinit__(self, double gamma):
self.gamma = gamma
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
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
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef image_types[:, ::1] compute_cmap(
default_types[:] data,
image_types[:, ::1] colors,
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 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 normalization: Normalization to apply
:return: Data converted to colors
"""
cdef image_types[:, ::1] output
cdef double scale, value, normalized_vmin, normalized_vmax
cdef int length, nb_channels, nb_colors
cdef int channel, index, lut_index, num_threads
nb_colors = <int> colors.shape[0]
nb_channels = <int> colors.shape[1]
length = <int> data.size
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)
if length < USE_OPENMP_THRESHOLD:
num_threads = 1
else:
num_threads = min(
DEFAULT_NUM_THREADS,
int(os.environ.get("OMP_NUM_THREADS", DEFAULT_NUM_THREADS)))
with nogil:
for index in prange(length, num_threads=num_threads):
value = normalization.apply_double(
<double> data[index], vmin, vmax)
# Handle NaN
if isnan(value):
for channel in range(nb_channels):
output[index, channel] = nan_color[channel]
continue
if value <= normalized_vmin:
lut_index = 0
elif value >= normalized_vmax:
lut_index = nb_colors - 1
else:
lut_index = <int>((value - normalized_vmin) * scale)
# Index can overflow of 1
if lut_index >= nb_colors:
lut_index = nb_colors - 1
for channel in range(nb_channels):
output[index, channel] = colors[lut_index, channel]
return output
@cython.wraparound(False)
@cython.boundscheck(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef image_types[:, ::1] compute_cmap_with_lut(
lut_types[:] data,
image_types[:, ::1] colors,
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 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 normalization: Normalization to apply
:return: The generated image
"""
cdef image_types[:, ::1] output
cdef double[:] values
cdef image_types[:, ::1] lut
cdef int type_min, type_max
cdef int nb_channels, length
cdef int channel, index, lut_index, num_threads
length = <int> data.size
nb_channels = <int> colors.shape[1]
if lut_types is cnumpy.int8_t:
type_min = -128
type_max = 127
elif lut_types is cnumpy.uint8_t:
type_min = 0
type_max = 255
elif lut_types is cnumpy.int16_t:
type_min = -32768
type_max = 32767
else: # uint16_t
type_min = 0
type_max = 65535
colors_dtype = numpy.array(colors).dtype
values = numpy.arange(type_min, type_max + 1, dtype=numpy.float64)
lut = compute_cmap(
values, colors, normalization, vmin, vmax, nan_color)
output = numpy.empty((length, nb_channels), dtype=colors_dtype)
if length < USE_OPENMP_THRESHOLD:
num_threads = 1
else:
num_threads = min(
DEFAULT_NUM_THREADS,
int(os.environ.get("OMP_NUM_THREADS", DEFAULT_NUM_THREADS)))
with nogil:
# Apply LUT
for index in prange(length, num_threads=num_threads):
lut_index = data[index] - type_min
for channel in range(nb_channels):
output[index, channel] = lut[lut_index, channel]
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,
Normalization normalization,
double vmin,
double vmax,
image_types[::1] nan_color):
"""Implementation of colormap.
Use :func:`cmap`.
:param data: Input data
:param colors: Colors look-up-table
: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 image_types[:, ::1] output
# 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, normalization, vmin, vmax, nan_color)
elif data_types in default_types: # Use default implementation
output = compute_cmap(
data, colors, normalization, vmin, vmax, nan_color)
else:
raise ValueError('Unsupported data type')
return numpy.array(output, copy=False)
def cmap(data not None,
colors not None,
double vmin,
double vmax,
normalization='linear',
nan_color=None):
"""Convert data to colors with provided colors look-up table.
:param numpy.ndarray data: The input data
:param numpy.ndarray colors: Color look-up table as a 2D array.
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 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
:return: Array of colors. The shape of the
returned array is that of data array + the last dimension of colors.
The dtype of the returned array is that of the colors array.
:rtype: numpy.ndarray
:raises ValueError: If data of colors dtype is not supported
"""
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)
if data.dtype.kind not in ('b', 'i', 'u', 'f'):
raise ValueError("Unsupported data dtype: %s" % data.dtype)
native_endian_dtype = data.dtype.newbyteorder('N')
if native_endian_dtype.kind == 'f' and native_endian_dtype.itemsize == 2:
native_endian_dtype = "=f4" # Use native float32 instead of float16
data = numpy.array(data, copy=False, dtype=native_endian_dtype)
# Make colors a contiguous array of native endian type
colors = numpy.array(colors, copy=False)
if colors.dtype.kind == 'f':
colors_dtype = numpy.dtype('float32')
elif colors.dtype.kind in ('b', 'i', 'u'):
colors_dtype = numpy.dtype('uint8')
else:
raise ValueError("Unsupported colors dtype: %s" % colors.dtype)
if (colors_dtype.kind != colors.dtype.kind or
colors_dtype.itemsize != colors.dtype.itemsize):
# Do not warn if only endianness has changed
_logger.warning("Casting colors from %s to %s", colors.dtype, colors_dtype)
nb_channels = colors.shape[colors.ndim - 1]
colors = numpy.ascontiguousarray(colors, dtype=colors_dtype)
# 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)
else:
nan_color = numpy.ascontiguousarray(
nan_color, dtype=colors.dtype).reshape(-1)
assert nan_color.shape == (nb_channels,)
image = _cmap(
data.reshape(-1),
colors.reshape(-1, nb_channels),
norm,
vmin,
vmax,
nan_color)
image.shape = data.shape + (nb_channels,)
return image
|