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
|
# -*- coding: utf-8 -*-
#cython: embedsignature=True, language_level=3
## This is for optimisation
#cython: boundscheck=False, wraparound=False, cdivision=True, initializedcheck=False,
## This is for developping:
##cython: profile=True, warn.undeclared=True, warn.unused=True, warn.unused_result=False, warn.unused_arg=True
#
# Project: Fable Input/Output
# https://github.com/silx-kit/fabio
#
# Copyright (C) 2020-2020 European Synchrotron Radiation Facility, Grenoble, France
#
# 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.
"""Densification of sparse frame format
"""
__author__ = "Jérôme Kieffer"
__date__ = "18/12/2020"
__contact__ = "Jerome.kieffer@esrf.fr"
__license__ = "MIT"
import time
import numpy
from libc.stdint cimport int8_t, uint8_t, \
uint16_t, int16_t,\
int32_t, uint32_t,\
int64_t, uint64_t
from libc.math cimport isfinite, log, sqrt, cos, M_PI
from libc.stdlib cimport rand, RAND_MAX, srand
ctypedef fused any_t:
double
float
int8_t
uint8_t
uint16_t
int16_t
int32_t
uint32_t
int64_t
uint64_t
cdef:
double EPS64 = numpy.finfo(numpy.float64).eps
double two_pi = 2.0*M_PI
cdef double random_normal(double mu, double sigma) nogil:
"""
Calculate the gaussian distribution using the Box–Muller algorithm
Credits:
https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
:param mu: the center of the distribution
:param sigma: the width of the distribution
"""
cdef:
double u1, u2, mag
u1 = u2 = 0.0
while (u1 < EPS64):
u1 = <double> rand() / RAND_MAX
u2 = <double> rand() / RAND_MAX
mag = sigma * sqrt(-2.0 * log(u1));
return mag * cos(two_pi * u2) + mu;
def densify(float[:,::1] mask,
float[::1] radius,
uint32_t[::1] index,
any_t[::1] intensity,
any_t dummy,
dtype,
float[::1] background,
float[::1] background_std=None):
"""
Densify a sparse representation to generate a normal frame
:param mask: 2D array with NaNs for mask and pixel radius for the valid pixels
:param radius: 1D array with the radial distance
:param background: 1D array with the background values at given distance from the center
:param index: position of non-background pixels
:param intensity: intensities of non background pixels (at index position)
:param dummy: numerical value for masked-out pixels in dense image
:param dtype: dtype of intensity.
:param background_std: 1D array with the background std at given distance from the center --> activates the noisy mode.
:return: dense frame as 2D array
"""
cdef:
Py_ssize_t i, j, size, pos, size_over, width, height
double value, fres, fpos, idelta, start, std
bint integral, noisy
any_t[:, ::1] dense
size = radius.shape[0]
assert background.shape[0] == size
size_over = index.shape[0]
assert intensity.shape[0] == size_over
integral = numpy.issubdtype(dtype, numpy.integer)
height =mask.shape[0]
width = mask.shape[1]
dense = numpy.zeros((height, width), dtype=dtype)
noisy = background_std is not None
if noisy:
try:
value = time.time_ns()
except Exception:
value = int(time.time()/EPS64)
srand(<unsigned int> (value%RAND_MAX))
with nogil:
start = radius[0]
idelta = (size - 1)/(radius[size-1] - start)
#Linear interpolation
for i in range(height):
for j in range(width):
fpos = (mask[i,j] - start)*idelta
if (fpos<0) or (fpos>=size) or (not isfinite(fpos)):
dense[i,j] = dummy
else:
pos = <uint32_t> fpos
if pos+1 == size:
value = background[pos]
fres = 0.0
else:
fres = fpos - pos
value = (1.0 - fres)*background[pos] + fres*background[pos+1]
if noisy:
std = (1.0 - fres)*background_std[pos] + fres*background_std[pos+1]
value = max(0.0, random_normal(value, std))
if integral:
dense[i,j] = <any_t>(value + 0.5) #this is rounding
else:
dense[i,j] = <any_t>(value)
# Assignment of outliers
for i in range(size_over):
j = index[i]
dense[j//width, j%width] = intensity[i]
return numpy.asarray(dense)
|
