peak_fit.py 12.9 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
Damien Naudet's avatar
Damien Naudet committed
36
from threading import Thread
37 38

import numpy as np
39
from scipy.optimize import leastsq
40

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

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

49 50 51 52

_logger = logging.getLogger(__name__)


53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70
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)

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

75 76 77 78 79 80 81
    elif mode == BackgroundTypes.NONE:
        return np.zeros_like(data)

    else:
        raise ValueError("Unsupported background mode")


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


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

Thomas Vincent's avatar
Thomas Vincent committed
92
    :param FitTypes fit_type:
Damien Naudet's avatar
Damien Naudet committed
93 94

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

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

    :param BackgroundTypes background: The background subtraction to perform
Damien Naudet's avatar
Damien Naudet committed
103
    """
Damien Naudet's avatar
Damien Naudet committed
104

Damien Naudet's avatar
Damien Naudet committed
105 106 107 108
    READY, RUNNING, DONE, ERROR, CANCELED = __STATUSES = range(5)

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

Thomas Vincent's avatar
Thomas Vincent committed
116
        self.__progress = 0
Damien Naudet's avatar
Damien Naudet committed
117 118
        self.__results = None
        self.__thread = None
Damien Naudet's avatar
Damien Naudet committed
119
        self.__callback = None
Damien Naudet's avatar
Damien Naudet committed
120 121 122 123 124 125 126

        self.__status = self.READY

        self.__indices = None

        self.__qspace_f = qspace_f
        self.__fit_type = fit_type
127
        self.__background = background
Damien Naudet's avatar
Damien Naudet committed
128

129
        self.__n_proc = n_proc if n_proc else config.DEFAULT_PROCESS_NUMBER
Damien Naudet's avatar
Damien Naudet committed
130 131

        if roi_indices is not None:
Damien Naudet's avatar
Damien Naudet committed
132
            self.__roi_indices = np.array(roi_indices[:])
Damien Naudet's avatar
Damien Naudet committed
133 134 135 136 137
        else:
            self.__roi_indices = None

        if fit_type not in FitTypes.ALLOWED:
            self.__set_status(self.ERROR)
138 139 140 141 142
            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))
143

Damien Naudet's avatar
Damien Naudet committed
144 145 146 147
        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
148

Damien Naudet's avatar
Damien Naudet committed
149
                n_points = qdata_shape[0]
150

Damien Naudet's avatar
Damien Naudet committed
151
                if indices is None:
Damien Naudet's avatar
Damien Naudet committed
152
                    indices = list(range(n_points))
Damien Naudet's avatar
Damien Naudet committed
153 154 155 156 157
                else:
                    indices = indices[:]
        except IOError:
            self.__set_status(self.ERROR)
            raise
Damien Naudet's avatar
Damien Naudet committed
158

Damien Naudet's avatar
Damien Naudet committed
159
        self.__indices = np.array(indices)
Damien Naudet's avatar
Damien Naudet committed
160 161 162 163 164 165 166

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

169
    def peak_fit(self, blocking=True, callback=None):
Damien Naudet's avatar
Damien Naudet committed
170 171 172 173 174 175 176 177
        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
178
            self.__thread = Thread(target=self.__peak_fit)
Damien Naudet's avatar
Damien Naudet committed
179
            self.__callback = callback
Thomas Vincent's avatar
Thomas Vincent committed
180
            self.__thread.start()
Damien Naudet's avatar
Damien Naudet committed
181 182

    def progress(self):
Thomas Vincent's avatar
Thomas Vincent committed
183
        return 100 * self.__progress / len(self.__indices)
Damien Naudet's avatar
Damien Naudet committed
184 185

    def __peak_fit(self):
Thomas Vincent's avatar
Thomas Vincent committed
186
        self.__progress = 0
Damien Naudet's avatar
Damien Naudet committed
187 188
        self.__set_status(self.RUNNING)

Thomas Vincent's avatar
Thomas Vincent committed
189
        pool = multiprocessing.Pool(self.__n_proc)
190

Thomas Vincent's avatar
Thomas Vincent committed
191 192 193 194 195 196 197 198 199 200
        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)
            self.__progress += 1
Damien Naudet's avatar
Damien Naudet committed
201 202 203
        pool.close()
        pool.join()

204
        # Prepare FitResult object
Thomas Vincent's avatar
Thomas Vincent committed
205
        with QSpaceH5.QSpaceH5(self.__qspace_f) as qspace_h5:
206 207
            x_pos = qspace_h5.sample_x[self.__indices]
            y_pos = qspace_h5.sample_y[self.__indices]
208
            q_dim0, q_dim1, q_dim2 = qspace_h5.qspace_dimension_values
Damien Naudet's avatar
Damien Naudet committed
209

210 211 212 213
        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
214

215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233
        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
234

235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250
        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
251 252 253

        self.__set_status(self.DONE)

Damien Naudet's avatar
Damien Naudet committed
254 255 256
        if self.__callback:
            self.__callback()

257 258 259
        return results


260 261 262 263
# Per process cache for fit process
_per_process_cache = None


264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
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]]]
    """
286 287 288 289 290 291 292 293 294 295 296
    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)
297 298 299 300 301 302 303 304 305 306 307 308 309 310

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

312 313 314 315
    # Background subtraction
    if background_type != BackgroundTypes.NONE:
        for array in projections:
            array -= background_estimation(background_type, array)
Damien Naudet's avatar
Damien Naudet committed
316

317 318 319 320
    # 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
321

322
    return result
Damien Naudet's avatar
Damien Naudet committed
323

324

325
# Center of mass
Damien Naudet's avatar
Damien Naudet committed
326

327 328
def centroid(axis, signal):
    """Returns Center of mass and maximum information
329

330 331 332 333 334 335 336 337 338
    :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
339

340 341 342 343 344 345 346
    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
347

348

349
# Gaussian fit
350

351 352
def _gaussian_err(parameters, axis, signal):
    """Returns difference between signal and given gaussian
Damien Naudet's avatar
Damien Naudet committed
353

354 355 356 357 358 359 360
    :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
361

Damien Naudet's avatar
Damien Naudet committed
362

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

Damien Naudet's avatar
Damien Naudet committed
365

366 367
def gaussian_fit(axis, signal):
    """Returns gaussian fitting information
Damien Naudet's avatar
Damien Naudet committed
368

369 370
    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
371

372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395
    :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