peak_fit.py 12.7 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 37

import numpy as np
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
44
from ...io.FitH5 import BackgroundTypes
45
from ...util import gaussian, project
46
from .fitresults import FitResult
47

48 49 50 51

_logger = logging.getLogger(__name__)


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

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

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

    else:
        raise ValueError("Unsupported background mode")


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


87 88
class PeakFitter(object):
    """Class performing fit/com processing
Damien Naudet's avatar
Damien Naudet committed
89

90
    :param str qspace_f: path to the HDF5 file containing the qspace cubes
Thomas Vincent's avatar
Thomas Vincent committed
91
    :param FitTypes fit_type:
92 93
    :param Union[numpy.ndarray,None] indices:
        indices of the cubes (in the input HDF5 dataset) for which
94
        the dim0/dim1/dim2 peaks coordinates will be computed. E.g : if the array
Damien Naudet's avatar
Damien Naudet committed
95
        [1, 2, 3] is provided, only those cubes will be fitted.
Thomas Vincent's avatar
Thomas Vincent committed
96 97
    :param Union[int,None] n_proc:
        Number of process to use. If None, the config value is used.
98
    :param Union[List[List[int]],None] roi_indices: QSpace ROI to process
99
    :param BackgroundTypes background: The background subtraction to perform
Damien Naudet's avatar
Damien Naudet committed
100
    """
Damien Naudet's avatar
Damien Naudet committed
101

Damien Naudet's avatar
Damien Naudet committed
102 103 104 105
    READY, RUNNING, DONE, ERROR, CANCELED = __STATUSES = range(5)

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

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

        self.__qspace_f = qspace_f
        self.__fit_type = fit_type
118
        self.__background = background
Damien Naudet's avatar
Damien Naudet committed
119

120
        self.__n_proc = n_proc if n_proc else config.DEFAULT_PROCESS_NUMBER
Damien Naudet's avatar
Damien Naudet committed
121 122

        if roi_indices is not None:
Damien Naudet's avatar
Damien Naudet committed
123
            self.__roi_indices = np.array(roi_indices[:])
Damien Naudet's avatar
Damien Naudet committed
124 125 126 127 128
        else:
            self.__roi_indices = None

        if fit_type not in FitTypes.ALLOWED:
            self.__set_status(self.ERROR)
129 130 131 132 133
            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))
134

Damien Naudet's avatar
Damien Naudet committed
135 136 137 138 139 140 141
        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
142

143 144 145 146 147
        if indices is None:
            n_points = qdata_shape[0]
            self.__indices = np.arange(n_points)
        else:
            self.__indices = np.array(indices, copy=True)
Damien Naudet's avatar
Damien Naudet committed
148 149 150 151 152 153 154

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

157 158
    def peak_fit(self, blocking=True):
        """Run fit/com processing
Damien Naudet's avatar
Damien Naudet committed
159

160 161 162 163
        :param bool blocking:
            True for blocking until the end of the processing,
            False for yielding progress during the processing.
        """
Damien Naudet's avatar
Damien Naudet committed
164 165 166
        self.__results = None

        if blocking:
167 168
            for _ in self.__peak_fit():
                pass
Damien Naudet's avatar
Damien Naudet committed
169
        else:
170 171
            for progress, _ in enumerate(self.__peak_fit()):
                yield progress / len(self.__indices)
Damien Naudet's avatar
Damien Naudet committed
172 173 174 175

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

Thomas Vincent's avatar
Thomas Vincent committed
176
        pool = multiprocessing.Pool(self.__n_proc)
177

Thomas Vincent's avatar
Thomas Vincent committed
178 179 180 181 182 183 184 185 186
        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)
187
            yield result
Damien Naudet's avatar
Damien Naudet committed
188 189 190
        pool.close()
        pool.join()

191
        # Prepare FitResult object
Thomas Vincent's avatar
Thomas Vincent committed
192
        with QSpaceH5.QSpaceH5(self.__qspace_f) as qspace_h5:
193 194
            x_pos = qspace_h5.sample_x[self.__indices]
            y_pos = qspace_h5.sample_y[self.__indices]
195
            q_dim0, q_dim1, q_dim2 = qspace_h5.qspace_dimension_values
Damien Naudet's avatar
Damien Naudet committed
196

197 198 199 200
        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
201

202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220
        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
221

222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237
        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
238 239
        self.__set_status(self.DONE)

240

241 242 243 244
# Per process cache for fit process
_per_process_cache = None


245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266
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]]]
    """
267 268 269 270 271 272 273 274 275 276 277
    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)
278 279 280 281 282 283 284 285 286 287 288 289 290 291

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

293 294 295 296
    # Background subtraction
    if background_type != BackgroundTypes.NONE:
        for array in projections:
            array -= background_estimation(background_type, array)
Damien Naudet's avatar
Damien Naudet committed
297

298 299 300 301
    # 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
302

303
    return result
Damien Naudet's avatar
Damien Naudet committed
304

305

306
# Center of mass
Damien Naudet's avatar
Damien Naudet committed
307

308 309
def centroid(axis, signal):
    """Returns Center of mass and maximum information
310

311 312 313 314 315 316 317 318 319
    :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
320

321 322 323 324 325 326 327
    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
328

329

330
# Gaussian fit
331

332 333
def _gaussian_err(parameters, axis, signal):
    """Returns difference between signal and given gaussian
Damien Naudet's avatar
Damien Naudet committed
334

335 336 337 338 339 340 341
    :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
342

Damien Naudet's avatar
Damien Naudet committed
343

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

Damien Naudet's avatar
Damien Naudet committed
346

347 348
def gaussian_fit(axis, signal):
    """Returns gaussian fitting information
Damien Naudet's avatar
Damien Naudet committed
349

350 351
    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
352

353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376
    :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