File: Peaks.cpp

package info (click to toggle)
bornagain 23.0-4
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 103,936 kB
  • sloc: cpp: 423,131; python: 40,997; javascript: 11,167; awk: 630; sh: 318; ruby: 173; xml: 130; makefile: 51; ansic: 24
file content (50 lines) | stat: -rw-r--r-- 1,966 bytes parent folder | download | duplicates (2)
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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Device/Analyze/Peaks.cpp
//! @brief     Implements function FindPeaks namespace Analyze.
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************

#include "Device/Analyze/Peaks.h"
#include "Base/Axis/Scale.h"
#include "Device/Data/DataUtil.h"
#include "Device/Data/Datafield.h"
#include <tspectrum.h> // third-party code, extracted from CERN ROOT (class TSpectrum2)

std::vector<std::pair<double, double>>
Analyze::FindPeaks(const Datafield& data, double sigma, const std::string& option, double threshold)
{
    const double2d_t arr = data.values2D();
    tspectrum::Spectrum2D spec;
    const auto peaks = spec.find_peaks(arr, sigma, option, threshold);

    // coordinates of peaks in histogram axes units
    std::vector<std::pair<double, double>> result;

    for (const auto& p : peaks) {
        const double row_value = p.first;
        const double col_value = p.second;

        const auto xaxis_index = static_cast<size_t>(col_value);
        const size_t yaxis_index = data.yAxis().size() - 1 - static_cast<size_t>(row_value);

        const Bin1D xbin = data.xAxis().bin(xaxis_index);
        const Bin1D ybin = data.yAxis().bin(yaxis_index);

        const double dx = col_value - static_cast<size_t>(col_value);
        const double dy = -1.0 * (row_value - static_cast<size_t>(row_value));

        const double x = xbin.center() + xbin.binSize() * dx;
        const double y = ybin.center() + ybin.binSize() * dy;

        result.emplace_back(x, y);
    }
    return result;
}