mapping.py 6.29 KB
Newer Older
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
# coding: utf-8
# /*##########################################################################
#
# Copyright (c) 2016-2017 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.
#
# ###########################################################################*/


__authors__ = ["J. Garriga"]
__license__ = "MIT"
29
__date__ = "04/06/2021"
30
31

import numpy
32
33
34
35
try:
    import tqdm
except ImportError:
    tqdm = None
36
37
38
39
from functools import partial
import multiprocessing
from multiprocessing import Pool
from scipy.optimize import curve_fit
40
from scipy.stats import skew, kurtosis
41
from silx.math.medianfilter import medfilt2d
42
43


44
def gaussian(x, a, b, c, d):
45
46
47
48
49
50
51
52
53
54
55
56
    """
    Function to calculate the Gaussian with constants a, b, and c

    :param float x: Value to evaluate
    :param float a: height of the curve's peak
    :param float b: position of the center of the peak
    :param float c: standard deviation
    :param float d: lowest value of the curve (value of the limits)

    :returns: result of the function on x
    :rtype: float
    """
57
58
    return d + a * numpy.exp(-numpy.power(x - b, 2) / (2 * numpy.power(c, 2)))

59

60
61
62
def generator(data, moments=None):
    """
    Generator that returns the rocking curve for every pixel
63

64
65
66
67
    :param ndarray data: data to analyse
    :param moments: array of same shape as data with the moments values per pixel and image, optional
    :type moments: Union[None, ndarray]
    """
68
69
70
71
72
    for i in range(data.shape[1]):
        for j in range(data.shape[2]):
            if moments is not None:
                yield data[:, i, j], moments[:, i, j]
            yield data[:, i, j], None
73

74

75
def fit_rocking_curve(y_values, values=None, num_points=None, int_thresh=None):
76
77
78
    """
    Fit rocking curve.

79
    :param tuple y_values: the first element is the dependent data and the second element are
80
81
82
83
84
85
86
87
88
89
90
        the moments to use as starting values for the fit
    :param values: The independent variable where the data is measured, optional
    :type values: Union[None, list]
    :param int num_points: Number of points to evaluate the data on, optional
    :param float int_thresh: Intensity threshold. If not None, only the rocking curves with
        higher ptp (range of values) are fitted, others are assumed to be noise and not important
        data. This parameter is used to accelerate the fit. Optional.

    :returns: If curve was fitted, the fitted curve, else item[0]
    :rtype: list
    """
91
    y, moments = y_values
92
    y = numpy.asanyarray(y)
93
94
    if int_thresh is not None and y.ptp() < int_thresh:
        return y
95
96
97
98
99
100
101
    x = values if values is not None else numpy.arange(len(y))
    if moments is not None:
        p0 = [y.ptp(), moments[0], moments[1], min(y)]
    else:
        mean = sum(x * y) / sum(y)
        sigma = numpy.sqrt(sum(y * (x - mean)**2) / sum(y))
        p0 = [y.ptp(), mean, sigma, min(y)]
102
    if numpy.isclose(p0[2], 0):
103
        return y
104
    if num_points is None:
105
106
107
        num_points = len(y)
    try:
        pars, cov = curve_fit(f=gaussian, xdata=x, ydata=y, p0=p0)
108
109
110
        y_gauss = gaussian(numpy.linspace(x[0], x[-1], num_points), *pars)
        y_gauss[numpy.isnan(y_gauss)] = 0
        return y_gauss
111
112
    except RuntimeError:
        return y
113
114


115
116
117
118
119
120
121
122
123
124
def fit_data(data, moments=None, values=None, int_thresh=15, _tqdm=False):
    """
    Fit data in axis 0 of data

    :param bool _tqdm: If True, execut fitting under tqdm library.
    :returns: fitted data
    """
    g = generator(data, moments)
    cpus = multiprocessing.cpu_count()
    with Pool(cpus - 1) as p:
125
        if _tqdm and tqdm is not None:
126
            new_data = list(tqdm.tqdm(p.imap(partial(fit_rocking_curve, values=values, int_thresh=15), g), total=data.shape[1] * data.shape[2]))
127
        else:
128
            new_data = p.map(partial(fit_rocking_curve, values=values, int_thresh=15), g)
129
    return numpy.array(new_data).T.reshape(data.shape)
130
131


132
def compute_moments(values, data):
133
134
135
136
137
    """
    Compute first, second, third and fourth moment of data on values.

    :returns: The four moments
    """
138
    stack = [values[i] * data[i] for i in range(len(data))]
139
140
    zsum = numpy.sum(data, axis=0)
    com = numpy.sum(stack, axis=0) / zsum
141
    com = medfilt2d(com.astype(numpy.float64))
142
    com[numpy.isnan(com)] = numpy.min(com[~numpy.isnan(com)])
143
    repeat_values = numpy.repeat(values, com.shape[0] * com.shape[1]).reshape(len(values), com.shape[0], com.shape[1])
144
    std = numpy.sqrt(numpy.sum((data[i] * ((repeat_values[i] - com)**2) for i in range(len(data))), axis=0)) / zsum
145
    std = medfilt2d(std.astype(numpy.float64))
146
    std[numpy.isnan(std)] = numpy.min(std[~numpy.isnan(std)])
147
148
149
    skews = skew(stack, axis=0)
    kurt = kurtosis(stack, axis=0)

150
    return com, std, skews, kurt
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172


def compute_rsm(H, W, d, ffz, mainx):

    pix_arr = numpy.meshgrid(numpy.arange(H), numpy.arange(W))
    pix_arr[0] = (pix_arr[0] - W / 2) * d
    pix_arr[1] = (H / 2 - pix_arr[1]) * d
    pix_arr[0] = numpy.arctan2((ffz - W / 2) * pix_arr[0], mainx)
    pix_arr[1] = numpy.arctan2(pix_arr[1], numpy.sqrt(ffz * ffz + mainx * mainx))
    return pix_arr


def compute_magnification(H, W, d, obx, obpitch, mainx):

    pix_arr = numpy.meshgrid(numpy.arange(H), numpy.arange(W))
    d1 = obx / numpy.cos(obpitch)
    d2 = mainx / numpy.cos(obpitch) - d1
    M = d2 / d1
    d /= M
    pix_arr[0] = (pix_arr[0] - W / 2) * d
    pix_arr[1] = (H / 2 - pix_arr[1]) * d
    return pix_arr