peak_fit.py 17 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__)


Thomas Vincent's avatar
Thomas Vincent committed
51 52
class FitTypes(object):
    """Kind of fit available"""
53

Thomas Vincent's avatar
Thomas Vincent committed
54 55
    GAUSSIAN = 0
    """1 gaussian fit"""
56

Thomas Vincent's avatar
Thomas Vincent committed
57 58
    CENTROID = 1
    """Center of mass and max"""
59

Thomas Vincent's avatar
Thomas Vincent committed
60
    ALLOWED = GAUSSIAN, CENTROID
61 62


Thomas Vincent's avatar
Thomas Vincent committed
63
class FitStatus(object):
Thomas Vincent's avatar
Thomas Vincent committed
64 65
    """Enum for the fit status

Thomas Vincent's avatar
Thomas Vincent committed
66 67 68 69 70 71 72 73
    Starting at 1 for compatibility reasons.
    """
    UNKNOWN, OK, FAILED = range(0, 3)


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

Thomas Vincent's avatar
Thomas Vincent committed
74
    It also allows to save as hdf5 with :meth:`to_fit_h5`.
Thomas Vincent's avatar
Thomas Vincent committed
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

    :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
Thomas Vincent's avatar
Thomas Vincent committed
122
    def available_result_names(self, dimension=None):
Thomas Vincent's avatar
Thomas Vincent committed
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
        """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)):
Thomas Vincent's avatar
Thomas Vincent committed
176
                for name in self.available_result_names:
Thomas Vincent's avatar
Thomas Vincent committed
177 178 179 180 181 182 183
                    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)


184 185
class PeakFitter(object):
    """Class performing fit/com processing
Damien Naudet's avatar
Damien Naudet committed
186

187
    :param str qspace_f: path to the HDF5 file containing the qspace cubes
Thomas Vincent's avatar
Thomas Vincent committed
188
    :param FitTypes fit_type:
189 190
    :param Union[numpy.ndarray,None] indices:
        indices of the cubes (in the input HDF5 dataset) for which
191
        the dim0/dim1/dim2 peaks coordinates will be computed. E.g : if the array
Damien Naudet's avatar
Damien Naudet committed
192
        [1, 2, 3] is provided, only those cubes will be fitted.
Thomas Vincent's avatar
Thomas Vincent committed
193 194
    :param Union[int,None] n_proc:
        Number of process to use. If None, the config value is used.
195
    :param Union[List[List[int]],None] roi_indices: QSpace ROI to process
196
    :param BackgroundTypes background: The background subtraction to perform
Damien Naudet's avatar
Damien Naudet committed
197
    """
Damien Naudet's avatar
Damien Naudet committed
198

Thomas Vincent's avatar
Thomas Vincent committed
199
    # Different status values
Damien Naudet's avatar
Damien Naudet committed
200 201 202 203
    READY, RUNNING, DONE, ERROR, CANCELED = __STATUSES = range(5)

    def __init__(self,
                 qspace_f,
Damien Naudet's avatar
Damien Naudet committed
204
                 fit_type=FitTypes.GAUSSIAN,
Damien Naudet's avatar
Damien Naudet committed
205 206
                 indices=None,
                 n_proc=None,
207
                 roi_indices=None,
Thomas Vincent's avatar
Thomas Vincent committed
208
                 background=BackgroundTypes.NONE):
Damien Naudet's avatar
Damien Naudet committed
209 210 211 212 213 214 215
        super(PeakFitter, self).__init__()

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

        self.__qspace_f = qspace_f
        self.__fit_type = fit_type
216
        self.__background = background
Damien Naudet's avatar
Damien Naudet committed
217

218
        self.__n_proc = n_proc if n_proc else config.DEFAULT_PROCESS_NUMBER
Damien Naudet's avatar
Damien Naudet committed
219 220

        if roi_indices is not None:
Thomas Vincent's avatar
Thomas Vincent committed
221
            self.__roi_indices = numpy.array(roi_indices[:])
Damien Naudet's avatar
Damien Naudet committed
222 223 224 225 226
        else:
            self.__roi_indices = None

        if fit_type not in FitTypes.ALLOWED:
            self.__set_status(self.ERROR)
227 228 229 230 231
            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))
232

Damien Naudet's avatar
Damien Naudet committed
233 234 235 236 237 238 239
        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
240

241 242
        if indices is None:
            n_points = qdata_shape[0]
Thomas Vincent's avatar
Thomas Vincent committed
243
            self.__indices = numpy.arange(n_points)
244
        else:
Thomas Vincent's avatar
Thomas Vincent committed
245
            self.__indices = numpy.array(indices, copy=True)
Damien Naudet's avatar
Damien Naudet committed
246 247

    def __set_status(self, status):
Thomas Vincent's avatar
Thomas Vincent committed
248 249 250 251
        """Set the status of the processing

        :param int status:
        """
Damien Naudet's avatar
Damien Naudet committed
252 253 254
        assert status in self.__STATUSES
        self.__status = status

Thomas Vincent's avatar
Thomas Vincent committed
255 256
    status = property(lambda self: self.__status,
                      doc="Status of the fit/COM process (int)")
Damien Naudet's avatar
Damien Naudet committed
257

Thomas Vincent's avatar
Thomas Vincent committed
258 259
    results = property(lambda self: self.__results,
                       doc="The result of the fit/COM (FitResult or None)")
Damien Naudet's avatar
Damien Naudet committed
260

261 262
    def peak_fit(self):
        """Blocking execution of fit/com processing.
Damien Naudet's avatar
Damien Naudet committed
263

264
        It returns the fit/com result
Thomas Vincent's avatar
Thomas Vincent committed
265 266

        :rtype: FitResult
267
        """
268 269 270
        for _ in self.__peak_fit():
            pass
        return self.results
Damien Naudet's avatar
Damien Naudet committed
271

272 273 274 275 276 277 278
    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
279 280

    def __peak_fit(self):
Thomas Vincent's avatar
Thomas Vincent committed
281
        """Run the fit/COM processing and set result and status"""
282
        self.__results = None
Damien Naudet's avatar
Damien Naudet committed
283 284
        self.__set_status(self.RUNNING)

Thomas Vincent's avatar
Thomas Vincent committed
285
        pool = multiprocessing.Pool(self.__n_proc)
286

Thomas Vincent's avatar
Thomas Vincent committed
287 288 289 290 291 292
        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,
Thomas Vincent's avatar
Thomas Vincent committed
293
                                  roi_indices=self.__roi_indices),
Thomas Vincent's avatar
Thomas Vincent committed
294 295
                self.__indices):
            fit_results.append(result)
296
            yield result
Damien Naudet's avatar
Damien Naudet committed
297 298 299
        pool.close()
        pool.join()

300
        # Prepare FitResult object
Thomas Vincent's avatar
Thomas Vincent committed
301
        with QSpaceH5.QSpaceH5(self.__qspace_f) as qspace_h5:
302 303
            x_pos = qspace_h5.sample_x[self.__indices]
            y_pos = qspace_h5.sample_y[self.__indices]
304
            q_dim0, q_dim1, q_dim2 = qspace_h5.qspace_dimension_values
Thomas Vincent's avatar
Thomas Vincent committed
305
            q_dim_names = qspace_h5.qspace_dimension_names
Damien Naudet's avatar
Damien Naudet committed
306

307 308 309 310
        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
311

312
        if self.__fit_type == FitTypes.GAUSSIAN:
Thomas Vincent's avatar
Thomas Vincent committed
313 314 315 316
            result_dtype = [('Area', numpy.float64),
                            ('Center', numpy.float64),
                            ('Sigma', numpy.float64),
                            ('Status', numpy.bool_)]
317 318

        elif self.__fit_type == FitTypes.CENTROID:
Thomas Vincent's avatar
Thomas Vincent committed
319 320 321 322 323
            result_dtype = [('COM', numpy.float64),
                            ('I_sum', numpy.float64),
                            ('I_max', numpy.float64),
                            ('Pos_max', numpy.float64),
                            ('Status', numpy.bool_)]
324 325 326

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

Thomas Vincent's avatar
Thomas Vincent committed
328 329 330 331 332 333 334 335 336
        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
337 338
        self.__set_status(self.DONE)

339

Thomas Vincent's avatar
Thomas Vincent committed
340 341
# Fit/COM function run through multiprocessing

342 343 344 345
# Per process cache for fit process
_per_process_cache = None


346 347 348 349
def _fit_process(index,
                 qspace_f,
                 fit_type=FitTypes.GAUSSIAN,
                 background_type=BackgroundTypes.NONE,
Thomas Vincent's avatar
Thomas Vincent committed
350
                 roi_indices=None):
351 352 353 354 355 356 357 358 359 360 361 362
    """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
Thomas Vincent's avatar
Thomas Vincent committed
363
    :param Union[List[List[int]],None] roi_indices:
364 365 366 367
        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]]]
    """
368 369 370 371 372 373 374 375 376 377 378
    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)
379 380

    # apply Qspace ROI
Thomas Vincent's avatar
Thomas Vincent committed
381 382 383 384
    if roi_indices is not None:
        dim0_slice = slice(roi_indices[0][0], roi_indices[0][1], 1)
        dim1_slice = slice(roi_indices[1][0], roi_indices[1][1], 1)
        dim2_slice = slice(roi_indices[2][0], roi_indices[2][1], 1)
385 386

        axes = [axis[roi] for axis, roi in
Thomas Vincent's avatar
Thomas Vincent committed
387 388 389
                zip(axes, (dim0_slice, dim1_slice, dim2_slice))]
        hits = hits[dim0_slice, dim1_slice, dim2_slice]
        qspace = qspace[dim0_slice, dim1_slice, dim2_slice]
390 391 392

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

394 395 396 397
    # Background subtraction
    if background_type != BackgroundTypes.NONE:
        for array in projections:
            array -= background_estimation(background_type, array)
Damien Naudet's avatar
Damien Naudet committed
398

Thomas Vincent's avatar
Thomas Vincent committed
399 400 401 402 403 404 405
    # Select Fit/COM function
    if fit_type == FitTypes.CENTROID:
        fit_func = centroid
    elif fit_type == FitTypes.GAUSSIAN:
        fit_func = gaussian_fit
    else:
        raise RuntimeError("Unsupported fit type: %s" % fit_type)
Damien Naudet's avatar
Damien Naudet committed
406

Thomas Vincent's avatar
Thomas Vincent committed
407 408
    # Fit/COM for projections on each axes
    return [fit_func(axis, signal) for axis, signal in zip(axes, projections)]
Damien Naudet's avatar
Damien Naudet committed
409

410

411
# Center of mass
Damien Naudet's avatar
Damien Naudet committed
412

413 414
def centroid(axis, signal):
    """Returns Center of mass and maximum information
415

416 417 418 419 420 421 422 423 424
    :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
425

426 427
    else:
        max_idx = signal.argmax()
Thomas Vincent's avatar
Thomas Vincent committed
428
        return (float(numpy.dot(axis, signal) / signal_sum),
429 430 431 432
                float(signal_sum),
                float(signal[max_idx]),
                float(axis[max_idx]),
                True)
Damien Naudet's avatar
Damien Naudet committed
433

434

435
# Gaussian fit
436

437 438
def _gaussian_err(parameters, axis, signal):
    """Returns difference between signal and given gaussian
Damien Naudet's avatar
Damien Naudet committed
439

440 441 442 443 444 445 446
    :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
447

Damien Naudet's avatar
Damien Naudet committed
448

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

Damien Naudet's avatar
Damien Naudet committed
451

452 453
def gaussian_fit(axis, signal):
    """Returns gaussian fitting information
Damien Naudet's avatar
Damien Naudet committed
454

455 456
    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
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
    :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
Thomas Vincent's avatar
Thomas Vincent committed
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


# background

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 numpy.ones_like(data) * numpy.nanmin(data)

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

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

    elif mode == BackgroundTypes.NONE:
        return numpy.zeros_like(data)

    else:
        raise ValueError("Unsupported background mode")