peak_fit.py 13.6 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
#!/usr/bin/python
# coding: utf8
# /*##########################################################################
#
# Copyright (c) 2015-2016 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.
#
# ###########################################################################*/

Damien Naudet's avatar
WIP  
Damien Naudet committed
27 28
from __future__ import absolute_import

29 30 31 32
__authors__ = ["D. Naudet"]
__date__ = "01/06/2016"
__license__ = "MIT"

33
import logging
Thomas Vincent's avatar
Thomas Vincent committed
34
import functools
35
import ctypes
Thomas Vincent's avatar
Thomas Vincent committed
36
import multiprocessing
Damien Naudet's avatar
Damien Naudet committed
37
from threading import Thread
38
import multiprocessing.sharedctypes as mp_sharedctypes
Thomas Vincent's avatar
Thomas Vincent committed
39
import warnings
40 41

import numpy as np
42
from scipy.optimize import leastsq
43

Thomas Vincent's avatar
Thomas Vincent committed
44 45
from silx.math.fit import snip1d

Thomas Vincent's avatar
Thomas Vincent committed
46
from ... import config
Damien Naudet's avatar
Damien Naudet committed
47
from ...io import QSpaceH5
48
from ...io.FitH5 import BackgroundTypes
49 50
from ...util import gaussian
from .fitresults import FitResult
51

52 53 54 55

_logger = logging.getLogger(__name__)


56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
def background_estimation(mode, data):
    """Estimates a background of mode kind for data

    :param BackgroundTypes mode: The kind of background to compute
    :param numpy.ndarray data: The data array to process
    :return: The estimated background as an array with same shape as input
    :rtype: numpy.ndarray
    :raises ValueError: In case mode is not known
    """
    # Background subtraction
    if mode == BackgroundTypes.CONSTANT:
        # Shift data so that smallest value is 0
        return np.ones_like(data) * np.nanmin(data)

    elif mode == BackgroundTypes.LINEAR:
        # Simple linear background
        return np.linspace(data[0], data[-1], num=len(data), endpoint=True)

74 75 76 77
    elif mode == BackgroundTypes.SNIP:
        # Using snip background
        return snip1d(data, snip_width=len(data))

78 79 80 81 82 83 84
    elif mode == BackgroundTypes.NONE:
        return np.zeros_like(data)

    else:
        raise ValueError("Unsupported background mode")


85 86 87 88 89 90
class FitTypes(object):
    """Kind of fit available"""
    ALLOWED = range(2)
    GAUSSIAN, CENTROID = ALLOWED


Damien Naudet's avatar
Damien Naudet committed
91
class PeakFitter(Thread):
Damien Naudet's avatar
Damien Naudet committed
92
    """
Thomas Vincent's avatar
Thomas Vincent committed
93
    :param str qspace_f: path to the HDF5 file containing the qspace cubes
Damien Naudet's avatar
Damien Naudet committed
94

Thomas Vincent's avatar
Thomas Vincent committed
95
    :param FitTypes fit_type:
Damien Naudet's avatar
Damien Naudet committed
96 97

    :param indices: indices of the cubes (in the input HDF5 dataset) for which
98
        the dim0/dim1/dim2 peaks coordinates will be computed. E.g : if the array
Damien Naudet's avatar
Damien Naudet committed
99
        [1, 2, 3] is provided, only those cubes will be fitted.
Thomas Vincent's avatar
Thomas Vincent committed
100
    :type indices: *optional* `array_like`
Damien Naudet's avatar
Damien Naudet committed
101

Thomas Vincent's avatar
Thomas Vincent committed
102 103
    :param Union[int,None] n_proc:
        Number of process to use. If None, the config value is used.
104 105

    :param BackgroundTypes background: The background subtraction to perform
Damien Naudet's avatar
Damien Naudet committed
106
    """
Damien Naudet's avatar
Damien Naudet committed
107

Damien Naudet's avatar
Damien Naudet committed
108 109 110 111
    READY, RUNNING, DONE, ERROR, CANCELED = __STATUSES = range(5)

    def __init__(self,
                 qspace_f,
Damien Naudet's avatar
Damien Naudet committed
112
                 fit_type=FitTypes.GAUSSIAN,
Damien Naudet's avatar
Damien Naudet committed
113 114
                 indices=None,
                 n_proc=None,
115
                 roi_indices=None,
Thomas Vincent's avatar
Thomas Vincent committed
116
                 background=BackgroundTypes.NONE):
Damien Naudet's avatar
Damien Naudet committed
117 118 119 120
        super(PeakFitter, self).__init__()

        self.__results = None
        self.__thread = None
Damien Naudet's avatar
Damien Naudet committed
121
        self.__callback = None
Damien Naudet's avatar
Damien Naudet committed
122 123 124 125 126 127 128

        self.__status = self.READY

        self.__indices = None

        self.__qspace_f = qspace_f
        self.__fit_type = fit_type
129
        self.__background = background
Damien Naudet's avatar
Damien Naudet committed
130

131
        self.__n_proc = n_proc if n_proc else config.DEFAULT_PROCESS_NUMBER
Damien Naudet's avatar
Damien Naudet committed
132
        self.__shared_progress = mp_sharedctypes.RawArray(ctypes.c_int32,
133
                                                          self.__n_proc)
Damien Naudet's avatar
Damien Naudet committed
134 135

        if roi_indices is not None:
Damien Naudet's avatar
Damien Naudet committed
136
            self.__roi_indices = np.array(roi_indices[:])
Damien Naudet's avatar
Damien Naudet committed
137 138 139 140 141
        else:
            self.__roi_indices = None

        if fit_type not in FitTypes.ALLOWED:
            self.__set_status(self.ERROR)
142 143 144 145 146
            raise ValueError('Unknown fit type: {0}'.format(fit_type))

        if background not in BackgroundTypes.ALLOWED:
            self.__set_status(self.ERROR)
            raise ValueError('Unknown background type: {}'.format(background))
147

Damien Naudet's avatar
Damien Naudet committed
148 149 150 151
        try:
            with QSpaceH5.QSpaceH5(qspace_f) as qspace_h5:
                with qspace_h5.qspace_dset_ctx() as dset:
                    qdata_shape = dset.shape
Damien Naudet's avatar
Damien Naudet committed
152

Damien Naudet's avatar
Damien Naudet committed
153
                n_points = qdata_shape[0]
154

Damien Naudet's avatar
Damien Naudet committed
155
                if indices is None:
Damien Naudet's avatar
Damien Naudet committed
156
                    indices = list(range(n_points))
Damien Naudet's avatar
Damien Naudet committed
157 158 159 160 161
                else:
                    indices = indices[:]
        except IOError:
            self.__set_status(self.ERROR)
            raise
Damien Naudet's avatar
Damien Naudet committed
162

Damien Naudet's avatar
Damien Naudet committed
163
        self.__indices = np.array(indices)
Damien Naudet's avatar
Damien Naudet committed
164 165 166 167 168 169 170

    def __set_status(self, status):
        assert status in self.__STATUSES
        self.__status = status

    status = property(lambda self: self.__status)

Damien Naudet's avatar
Damien Naudet committed
171 172
    results = property(lambda self: self.__results)

173
    def peak_fit(self, blocking=True, callback=None):
Damien Naudet's avatar
Damien Naudet committed
174 175 176 177 178 179 180 181
        if self.__thread and self.__thread.is_alive():
            raise RuntimeError('A fit is already in progress.')

        self.__results = None

        if blocking:
            return self.__peak_fit()
        else:
Thomas Vincent's avatar
Thomas Vincent committed
182
            self.__thread = Thread(target=self.__peak_fit)
Damien Naudet's avatar
Damien Naudet committed
183
            self.__callback = callback
Thomas Vincent's avatar
Thomas Vincent committed
184
            self.__thread.start()
Damien Naudet's avatar
Damien Naudet committed
185 186 187 188 189 190 191 192 193

    def progress(self):
        return (100.0 *
                np.frombuffer(self.__shared_progress, dtype='int32').max() /
                (len(self.__indices) - 1))

    def __peak_fit(self):
        self.__set_status(self.RUNNING)

194
        # TODO
195
        progress = np.frombuffer(self.__shared_progress, dtype='int32')
Damien Naudet's avatar
Damien Naudet committed
196
        progress[:] = 0
197

Thomas Vincent's avatar
Thomas Vincent committed
198
        pool = multiprocessing.Pool(self.__n_proc)
199 200 201 202 203 204 205
        fit_results = pool.map(
            functools.partial(_fit_process,
                              qspace_f=self.__qspace_f,
                              fit_type=self.__fit_type,
                              background_type=self.__background,
                              roiIndices=self.__roi_indices),
            self.__indices)
Damien Naudet's avatar
Damien Naudet committed
206 207 208
        pool.close()
        pool.join()

209
        # Prepare FitResult object
Thomas Vincent's avatar
Thomas Vincent committed
210
        with QSpaceH5.QSpaceH5(self.__qspace_f) as qspace_h5:
211 212
            x_pos = qspace_h5.sample_x[self.__indices]
            y_pos = qspace_h5.sample_y[self.__indices]
213
            q_dim0, q_dim1, q_dim2 = qspace_h5.qspace_dimension_values
Damien Naudet's avatar
Damien Naudet committed
214

215 216 217 218
        if self.__roi_indices is not None:
            q_dim0 = q_dim0[self.__roi_indices[0][0]:self.__roi_indices[0][1]]
            q_dim1 = q_dim1[self.__roi_indices[1][0]:self.__roi_indices[1][1]]
            q_dim2 = q_dim2[self.__roi_indices[2][0]:self.__roi_indices[2][1]]
Damien Naudet's avatar
Damien Naudet committed
219

220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238
        if self.__fit_type == FitTypes.GAUSSIAN:
            fit_name = 'Gaussian'
            result_name = 'gauss_0'
            result_dtype = [('Area', np.float64),
                            ('Center', np.float64),
                            ('Sigma', np.float64),
                            ('Status', np.bool_)]

        elif self.__fit_type == FitTypes.CENTROID:
            fit_name = 'Centroid'
            result_name = 'centroid'
            result_dtype = [('COM', np.float64),
                            ('I_sum', np.float64),
                            ('I_max', np.float64),
                            ('Pos_max', np.float64),
                            ('Status', np.bool_)]

        else:
            raise RuntimeError('Unknown Fit Type')
Damien Naudet's avatar
Damien Naudet committed
239

240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
        results = FitResult(entry=fit_name,
                            sample_x=x_pos,
                            sample_y=y_pos,
                            q_x=q_dim0,
                            q_y=q_dim1,
                            q_z=q_dim2,
                            background_mode=self.__background)
        fit_results = np.array(fit_results, dtype=result_dtype)
        # From points x axes to axes x points
        fit_results = np.transpose(fit_results)
        for axis_index, array in enumerate(fit_results):
            for name, _ in result_dtype[:-1]:
                results._add_axis_result(result_name, axis_index, name, array[name])
            results._set_axis_status(axis_index, array['Status'])

        self.__results = results
Damien Naudet's avatar
Damien Naudet committed
256 257 258

        self.__set_status(self.DONE)

Damien Naudet's avatar
Damien Naudet committed
259 260 261
        if self.__callback:
            self.__callback()

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
        return results


def _fit_process(index,
                 qspace_f,
                 fit_type=FitTypes.GAUSSIAN,
                 background_type=BackgroundTypes.NONE,
                 roiIndices=None):
    """Run fit processing.

    It loads a QSpace, extracts a ROI from it, project to axes,
    and then for each axis, it subtracts a background and performs a fit/COM.

    This function is run through a multiprocessing.Pool

    :param int index: The index of the QSpace to process
    :param str qspace_f: Filename of the hdf5 file containing QSpace
    :param FitTypes fit_type: The kind of fit to perform
    :param BackgroundTypes background_type:
        The kind of background subtraction to perform
    :param Union[List[List[int]],None] roiIndices:
        Optional QSpace ROI start:end in the 3 dimensions
    :return: Fit results as a list of results for dim0, dim1 and dim2
    :rtype: List[List[Union[float,bool]]]
    """
    # Read data from file
    with QSpaceH5.QSpaceH5(qspace_f) as qspace_h5:
        axes = qspace_h5.qspace_dimension_values
        hits = qspace_h5.histo
        qspace = qspace_h5.qspace_slice(index)

    # apply Qspace ROI
    if roiIndices is not None:
        dim0Slice = slice(roiIndices[0][0], roiIndices[0][1], 1)
        dim1Slice = slice(roiIndices[1][0], roiIndices[1][1], 1)
        dim2Slice = slice(roiIndices[2][0], roiIndices[2][1], 1)

        axes = [axis[roi] for axis, roi in
                zip(axes, (dim0Slice, dim1Slice, dim2Slice))]
        hits = hits[dim0Slice, dim1Slice, dim2Slice]
        qspace = qspace[dim0Slice, dim1Slice, dim2Slice]

    # Normalize with hits and project to axes
    projections = project(qspace, hits)
306

307 308 309 310
    # Background subtraction
    if background_type != BackgroundTypes.NONE:
        for array in projections:
            array -= background_estimation(background_type, array)
Damien Naudet's avatar
Damien Naudet committed
311

312 313 314 315
    # Fit/COM
    fit = {FitTypes.CENTROID: centroid,
           FitTypes.GAUSSIAN: gaussian_fit}[fit_type]
    result = [fit(axis, values) for axis, values in zip(axes, projections)]
Damien Naudet's avatar
Damien Naudet committed
316

317
    return result
Damien Naudet's avatar
Damien Naudet committed
318

319

320 321
def project(data, hits=None):
    """Sum data along each axis
Damien Naudet's avatar
Damien Naudet committed
322

323 324 325 326 327 328 329
    :param numpy.ndarray data: 3D histogram
    :param Union[numpy.ndarray,None] hits:
        Number of bin count of the histogram or None to ignore
    :return: Projections on each axis of the dataset
    :rtype: List[numpy.ndarray]
    """
    if hits is not None:
Thomas Vincent's avatar
Thomas Vincent committed
330 331 332
        with warnings.catch_warnings():
            warnings.simplefilter('ignore', RuntimeWarning)
            data /= hits
333
        data[hits <= 0] = 0
334

335 336 337 338 339
    dim2_sum = data.sum(axis=0).sum(axis=0)
    data_sum_dim2 = data.sum(axis=2)
    dim1_sum = data_sum_dim2.sum(axis=0)
    dim0_sum = data_sum_dim2.sum(axis=1)
    return dim0_sum, dim1_sum, dim2_sum
340

341

342
# Center of mass
Damien Naudet's avatar
Damien Naudet committed
343

344 345
def centroid(axis, signal):
    """Returns Center of mass and maximum information
346

347 348 349 350 351 352 353 354 355
    :param numpy.ndarray axis: 1D x data
    :param numpy.ndarray signal: 1D y data
    :return: Center of mass, sum of signal, max of signal, position of max and status
        ('COM', 'I_sum', 'I_max', 'Pos_max', 'status')
    :rtype: List[Union[float,bool]]
    """
    signal_sum = signal.sum()
    if signal_sum == 0:
        return float('nan'), float('nan'), float('nan'), float('nan'), False
Damien Naudet's avatar
Damien Naudet committed
356

357 358 359 360 361 362 363
    else:
        max_idx = signal.argmax()
        return (float(np.dot(axis, signal) / signal_sum),
                float(signal_sum),
                float(signal[max_idx]),
                float(axis[max_idx]),
                True)
Damien Naudet's avatar
Damien Naudet committed
364

365

366
# Gaussian fit
367

368 369
def _gaussian_err(parameters, axis, signal):
    """Returns difference between signal and given gaussian
Damien Naudet's avatar
Damien Naudet committed
370

371 372 373 374 375 376 377
    :param List[float] parameters: area, center, sigma
    :param numpy.ndarray axis: 1D x data
    :param numpy.ndarray signal: 1D y data
    :return:
    """
    area, center, sigma = parameters
    return gaussian(axis, area, center, sigma) - signal
378

Damien Naudet's avatar
Damien Naudet committed
379

380
_SQRT_2_PI = np.sqrt(2 * np.pi)
Damien Naudet's avatar
Damien Naudet committed
381

Damien Naudet's avatar
Damien Naudet committed
382

383 384
def gaussian_fit(axis, signal):
    """Returns gaussian fitting information
Damien Naudet's avatar
Damien Naudet committed
385

386 387
    parameters: (a, c, s)
    and f(x) = (a / (sqrt(2 * pi) * s)) * exp(- 0.5 * ((x - c) / s)^2)
Damien Naudet's avatar
Damien Naudet committed
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
    :param axis: 1D x data
    :param signal: 1D y data
    :return: Area, center, sigma, status
        ('Area', 'Center', 'Sigma', 'status)
    :rtype: List[Union[float,bool]]
    """
    # compute guess
    area = signal.sum() * (axis[-1] - axis[0]) / len(axis)
    center = axis[signal.argmax()]
    sigma = area / (signal.max() * _SQRT_2_PI)

    # Fit a gaussian
    result = leastsq(_gaussian_err,
                     x0=(area, center, sigma),
                     args=(axis, signal),
                     maxfev=100000,
                     full_output=True)

    if result[4] not in [1, 2, 3, 4]:
        return float('nan'), float('nan'), float('nan'), False

    else:
        area, center, sigma = result[0]
        return float(area), float(center), float(sigma), True