peak_fit.py 16.3 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.
#
# ###########################################################################*/

Thomas Vincent's avatar
Thomas Vincent committed
27
from __future__ import absolute_import, division
Damien Naudet's avatar
WIP  
Damien Naudet committed
28

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

33
import logging
Thomas Vincent's avatar
Thomas Vincent committed
34 35
import functools
import multiprocessing
36

Thomas Vincent's avatar
Thomas Vincent committed
37
import numpy
38
from scipy.optimize import leastsq
39

Thomas Vincent's avatar
Thomas Vincent committed
40 41
from silx.math.fit import snip1d

Thomas Vincent's avatar
Thomas Vincent committed
42
from ... import config
Damien Naudet's avatar
Damien Naudet committed
43
from ...io import QSpaceH5
Thomas Vincent's avatar
Thomas Vincent committed
44
from ...io.FitH5 import BackgroundTypes, FitH5Writer
45
from ...util import gaussian, project
46

47 48 49 50

_logger = logging.getLogger(__name__)


51 52 53 54 55 56 57 58 59 60 61 62
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
Thomas Vincent's avatar
Thomas Vincent committed
63
        return numpy.ones_like(data) * numpy.nanmin(data)
64 65 66

    elif mode == BackgroundTypes.LINEAR:
        # Simple linear background
Thomas Vincent's avatar
Thomas Vincent committed
67
        return numpy.linspace(data[0], data[-1], num=len(data), endpoint=True)
68

69 70 71 72
    elif mode == BackgroundTypes.SNIP:
        # Using snip background
        return snip1d(data, snip_width=len(data))

73
    elif mode == BackgroundTypes.NONE:
Thomas Vincent's avatar
Thomas Vincent committed
74
        return numpy.zeros_like(data)
75 76 77 78 79

    else:
        raise ValueError("Unsupported background mode")


80 81 82 83 84 85
class FitTypes(object):
    """Kind of fit available"""
    ALLOWED = range(2)
    GAUSSIAN, CENTROID = ALLOWED


Thomas Vincent's avatar
Thomas Vincent committed
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
class FitStatus(object):
    """
    Enum for the fit status
    Starting at 1 for compatibility reasons.
    """
    UNKNOWN, OK, FAILED = range(0, 3)


class FitResult(object):
    """Object storing fit/com results

    It also allows to save as hdf5.

    :param numpy.ndarray sample_x: N X sample position of the results
    :param numpy.ndarray sample_y: N Y sample position of the results
    :param List[numpy.ndarray] q_dim_values:
        Values along each axis of the QSpace
    :param List[str] q_dim_names:
        Name of axes for each dimension of the QSpace
    :param FitTypes fit_mode: Kind of fit
    :param BackgroundTypes background_mode: Kind of background subtraction
    :param numpy.ndarray fit_results:
        The fit/com results as a N (points) x 3 (axes) array of struct
        containing the results.

        Warning: This array is used as is and not copied.
    """

    def __init__(self,
                 sample_x, sample_y,
                 q_dim_values,
                 q_dim_names,
                 fit_mode, background_mode,
                 fit_results):

        super(FitResult, self).__init__()

        self.sample_x = sample_x
        """X position on the sample of each fit result (numpy.ndarray)"""

        self.sample_y = sample_y
        """Y position on the sample of each fit result (numpy.ndarray)"""

        self.qspace_dimension_values = q_dim_values
        """QSpace axis values (List[numpy.ndarray])"""

        self.qspace_dimension_names = q_dim_names
        """QSpace axis names (List[str])"""

        self.fit_mode = fit_mode
        """Fit type (FitTypes)"""

        self.background_mode = background_mode
        """Background type (BackgroundTypes)"""

        # transpose from N (points) x 3 (axes) to 3 (axes) x N (points)
        self._fit_results = numpy.transpose(fit_results)

    @property
    def available_results(self, dimension=None):
        """Returns the available result names

        :param Union[int,None] dimension:
        :rtype: List[str]
        """
        if dimension is None:
            dimension = 0
        return self._fit_results[dimension].dtype.names

    def get_results(self, dimension, parameter, copy=True):
        """Returns a given parameter of the result

        :param int dimension: QSpace dimension from which to return result
        :param str parameter: Name of the result to return
        :param bool copy: True to return a copy, False to return internal data
        :return: A 1D array
        :rtype: numpy.ndarray
        """
        return numpy.array(self._fit_results[dimension][parameter], copy=copy)

    def to_fit_h5(self, fit_h5, mode=None):
        """Write fit results to an HDF5 file

        :param str fit_h5: Filename where to save fit results
        :param Union[None,str] mode: HDF5 file opening mode
        """
        if self.fit_mode == FitTypes.GAUSSIAN:
            fit_name = 'Gaussian'
            result_name = 'gauss_0'
        elif self.fit_mode == FitTypes.CENTROID:
            fit_name = 'Centroid'
            result_name = 'centroid'
        else:
            raise RuntimeError('Unknown Fit Type')

        with FitH5Writer(fit_h5, mode=mode) as fitH5:
            fitH5.create_entry(fit_name)

            fitH5.set_scan_x(fit_name, self.sample_x)
            fitH5.set_scan_y(fit_name, self.sample_y)

            q_dim0, q_dim1, q_dim2 = self.qspace_dimension_values
            fitH5.set_qx(fit_name, q_dim0)
            fitH5.set_qy(fit_name, q_dim1)
            fitH5.set_qz(fit_name, q_dim2)

            fitH5.set_background_mode(fit_name, self.background_mode)
            fitH5.create_process(fit_name, result_name)

            for array, func, axis in zip(
                    self._fit_results,
                    (fitH5.set_qx_result, fitH5.set_qy_result, fitH5.set_qz_result),
                    (0, 1, 2)):
                for name in self.available_results:
                    results = self.get_results(axis, name, copy=False)
                    if name == 'Status':
                        fitH5.set_status(fit_name, axis, results)
                    else:
                        func(fit_name, result_name, name, results)


207 208
class PeakFitter(object):
    """Class performing fit/com processing
Damien Naudet's avatar
Damien Naudet committed
209

210
    :param str qspace_f: path to the HDF5 file containing the qspace cubes
Thomas Vincent's avatar
Thomas Vincent committed
211
    :param FitTypes fit_type:
212 213
    :param Union[numpy.ndarray,None] indices:
        indices of the cubes (in the input HDF5 dataset) for which
214
        the dim0/dim1/dim2 peaks coordinates will be computed. E.g : if the array
Damien Naudet's avatar
Damien Naudet committed
215
        [1, 2, 3] is provided, only those cubes will be fitted.
Thomas Vincent's avatar
Thomas Vincent committed
216 217
    :param Union[int,None] n_proc:
        Number of process to use. If None, the config value is used.
218
    :param Union[List[List[int]],None] roi_indices: QSpace ROI to process
219
    :param BackgroundTypes background: The background subtraction to perform
Damien Naudet's avatar
Damien Naudet committed
220
    """
Damien Naudet's avatar
Damien Naudet committed
221

Damien Naudet's avatar
Damien Naudet committed
222 223 224 225
    READY, RUNNING, DONE, ERROR, CANCELED = __STATUSES = range(5)

    def __init__(self,
                 qspace_f,
Damien Naudet's avatar
Damien Naudet committed
226
                 fit_type=FitTypes.GAUSSIAN,
Damien Naudet's avatar
Damien Naudet committed
227 228
                 indices=None,
                 n_proc=None,
229
                 roi_indices=None,
Thomas Vincent's avatar
Thomas Vincent committed
230
                 background=BackgroundTypes.NONE):
Damien Naudet's avatar
Damien Naudet committed
231 232 233 234 235 236 237
        super(PeakFitter, self).__init__()

        self.__results = None
        self.__status = self.READY

        self.__qspace_f = qspace_f
        self.__fit_type = fit_type
238
        self.__background = background
Damien Naudet's avatar
Damien Naudet committed
239

240
        self.__n_proc = n_proc if n_proc else config.DEFAULT_PROCESS_NUMBER
Damien Naudet's avatar
Damien Naudet committed
241 242

        if roi_indices is not None:
Thomas Vincent's avatar
Thomas Vincent committed
243
            self.__roi_indices = numpy.array(roi_indices[:])
Damien Naudet's avatar
Damien Naudet committed
244 245 246 247 248
        else:
            self.__roi_indices = None

        if fit_type not in FitTypes.ALLOWED:
            self.__set_status(self.ERROR)
249 250 251 252 253
            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))
254

Damien Naudet's avatar
Damien Naudet committed
255 256 257 258 259 260 261
        try:
            with QSpaceH5.QSpaceH5(qspace_f) as qspace_h5:
                with qspace_h5.qspace_dset_ctx() as dset:
                    qdata_shape = dset.shape
        except IOError:
            self.__set_status(self.ERROR)
            raise
Damien Naudet's avatar
Damien Naudet committed
262

263 264
        if indices is None:
            n_points = qdata_shape[0]
Thomas Vincent's avatar
Thomas Vincent committed
265
            self.__indices = numpy.arange(n_points)
266
        else:
Thomas Vincent's avatar
Thomas Vincent committed
267
            self.__indices = numpy.array(indices, copy=True)
Damien Naudet's avatar
Damien Naudet committed
268 269 270 271 272 273 274

    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
275 276
    results = property(lambda self: self.__results)

277 278
    def peak_fit(self):
        """Blocking execution of fit/com processing.
Damien Naudet's avatar
Damien Naudet committed
279

280
        It returns the fit/com result
281
        """
282 283 284
        for _ in self.__peak_fit():
            pass
        return self.results
Damien Naudet's avatar
Damien Naudet committed
285

286 287 288 289 290 291 292
    def peak_fit_iterator(self):
        """Run fit/com processing as a generator.

        It yields a progress indicator in [0, 1].
        """
        for progress, _ in enumerate(self.__peak_fit()):
            yield progress / len(self.__indices)
Damien Naudet's avatar
Damien Naudet committed
293 294

    def __peak_fit(self):
295
        self.__results = None
Damien Naudet's avatar
Damien Naudet committed
296 297
        self.__set_status(self.RUNNING)

Thomas Vincent's avatar
Thomas Vincent committed
298
        pool = multiprocessing.Pool(self.__n_proc)
299

Thomas Vincent's avatar
Thomas Vincent committed
300 301 302 303 304 305 306 307 308
        fit_results = []
        for result in pool.imap(
                functools.partial(_fit_process,
                                  qspace_f=self.__qspace_f,
                                  fit_type=self.__fit_type,
                                  background_type=self.__background,
                                  roiIndices=self.__roi_indices),
                self.__indices):
            fit_results.append(result)
309
            yield result
Damien Naudet's avatar
Damien Naudet committed
310 311 312
        pool.close()
        pool.join()

313
        # Prepare FitResult object
Thomas Vincent's avatar
Thomas Vincent committed
314
        with QSpaceH5.QSpaceH5(self.__qspace_f) as qspace_h5:
315 316
            x_pos = qspace_h5.sample_x[self.__indices]
            y_pos = qspace_h5.sample_y[self.__indices]
317
            q_dim0, q_dim1, q_dim2 = qspace_h5.qspace_dimension_values
Thomas Vincent's avatar
Thomas Vincent committed
318
            q_dim_names = qspace_h5.qspace_dimension_names
Damien Naudet's avatar
Damien Naudet committed
319

320 321 322 323
        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
324

325
        if self.__fit_type == FitTypes.GAUSSIAN:
Thomas Vincent's avatar
Thomas Vincent committed
326 327 328 329
            result_dtype = [('Area', numpy.float64),
                            ('Center', numpy.float64),
                            ('Sigma', numpy.float64),
                            ('Status', numpy.bool_)]
330 331

        elif self.__fit_type == FitTypes.CENTROID:
Thomas Vincent's avatar
Thomas Vincent committed
332 333 334 335 336
            result_dtype = [('COM', numpy.float64),
                            ('I_sum', numpy.float64),
                            ('I_max', numpy.float64),
                            ('Pos_max', numpy.float64),
                            ('Status', numpy.bool_)]
337 338 339

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

Thomas Vincent's avatar
Thomas Vincent committed
341 342 343 344 345 346 347 348 349
        self.__results = FitResult(
            sample_x=x_pos,
            sample_y=y_pos,
            q_dim_values=(q_dim0, q_dim1, q_dim2),
            q_dim_names=q_dim_names,
            fit_mode=self.__fit_type,
            background_mode=self.__background,
            fit_results=numpy.array(fit_results, dtype=result_dtype))

Damien Naudet's avatar
Damien Naudet committed
350 351
        self.__set_status(self.DONE)

352

353 354 355 356
# Per process cache for fit process
_per_process_cache = None


357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378
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]]]
    """
379 380 381 382 383 384 385 386 387 388 389
    global _per_process_cache
    if _per_process_cache is None:  # Initialize per process cache
        qspace_h5 = QSpaceH5.QSpaceH5(qspace_f)
        _per_process_cache = (
            qspace_h5, qspace_h5.qspace_dimension_values, qspace_h5.histo)

    # Retrieve data/file from cache
    qspace_h5, axes, hits = _per_process_cache

    # Load qspace
    qspace = qspace_h5.qspace_slice(index)
390 391 392 393 394 395 396 397 398 399 400 401 402 403

    # 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)
404

405 406 407 408
    # Background subtraction
    if background_type != BackgroundTypes.NONE:
        for array in projections:
            array -= background_estimation(background_type, array)
Damien Naudet's avatar
Damien Naudet committed
409

410 411 412 413
    # 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
414

415
    return result
Damien Naudet's avatar
Damien Naudet committed
416

417

418
# Center of mass
Damien Naudet's avatar
Damien Naudet committed
419

420 421
def centroid(axis, signal):
    """Returns Center of mass and maximum information
422

423 424 425 426 427 428 429 430 431
    :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
432

433 434
    else:
        max_idx = signal.argmax()
Thomas Vincent's avatar
Thomas Vincent committed
435
        return (float(numpy.dot(axis, signal) / signal_sum),
436 437 438 439
                float(signal_sum),
                float(signal[max_idx]),
                float(axis[max_idx]),
                True)
Damien Naudet's avatar
Damien Naudet committed
440

441

442
# Gaussian fit
443

444 445
def _gaussian_err(parameters, axis, signal):
    """Returns difference between signal and given gaussian
Damien Naudet's avatar
Damien Naudet committed
446

447 448 449 450 451 452 453
    :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
454

Damien Naudet's avatar
Damien Naudet committed
455

Thomas Vincent's avatar
Thomas Vincent committed
456
_SQRT_2_PI = numpy.sqrt(2 * numpy.pi)
Damien Naudet's avatar
Damien Naudet committed
457

Damien Naudet's avatar
Damien Naudet committed
458

459 460
def gaussian_fit(axis, signal):
    """Returns gaussian fitting information
Damien Naudet's avatar
Damien Naudet committed
461

462 463
    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
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
    :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