Skip to content

Commit

Permalink
Merge pull request #37159 from mantidproject/37107_findsxpeaks_peak_s…
Browse files Browse the repository at this point in the history 
…ize_validation

37107 findsxpeaks peak size validation
  • Loading branch information
@robertapplin
robertapplin committed Apr 26, 2024
2 parents 3815977 + e18a017 commit 4abe275
Show file tree
Hide file tree
Showing 7 changed files with 274 additions and 17 deletions.
14 changes: 14 additions & 0 deletions Framework/Crystal/inc/MantidCrystal/FindSXPeaksHelper.h
View file
Expand Up @@ -70,6 +70,8 @@ struct MANTID_CRYSTAL_DLL SXPeak {

detid_t getDetectorId() const;

const std::vector<int> &getPeakSpectras() const;

private:
/// TOF for the peak centre
double m_tof;
Expand Down Expand Up @@ -157,8 +159,11 @@ class MANTID_CRYSTAL_DLL PeakFindingStrategy {
virtual ~PeakFindingStrategy() = default;
PeakList findSXPeaks(const HistogramData::HistogramX &x, const HistogramData::HistogramY &y,
const HistogramData::HistogramE &e, const int workspaceIndex) const;
void setMinNBinsPerPeak(int minBinsPerPeak);

protected:
void filterPeaksForMinBins(std::vector<std::unique_ptr<PeakContainer>> &inputPeakList) const;

BoundsIterator getBounds(const HistogramData::HistogramX &x) const;
double calculatePhi(size_t workspaceIndex) const;
double getXValue(const HistogramData::HistogramX &x, const size_t peakLocation) const;
Expand All @@ -174,6 +179,7 @@ class MANTID_CRYSTAL_DLL PeakFindingStrategy {
const double m_maxValue = EMPTY_DBL();
const API::SpectrumInfo &m_spectrumInfo;
const XAxisUnit m_units;
int m_minNBinsPerPeak = EMPTY_INT();
};

class MANTID_CRYSTAL_DLL StrongestPeaksStrategy : public PeakFindingStrategy {
Expand Down Expand Up @@ -263,14 +269,22 @@ class MANTID_CRYSTAL_DLL ReducePeakListStrategy {
virtual std::vector<SXPeak> reduce(const std::vector<SXPeak> &peaks,
Mantid::Kernel::ProgressBase &progress) const = 0;

void setMinNSpectraPerPeak(int minSpectrasForPeak);
void setMaxNSpectraPerPeak(int maxSpectrasForPeak);

protected:
const CompareStrategy *m_compareStrategy;
int m_minNSpectraPerPeak = EMPTY_INT();
int m_maxNSpectraPerPeak = EMPTY_INT();
};

class MANTID_CRYSTAL_DLL SimpleReduceStrategy : public ReducePeakListStrategy {
public:
SimpleReduceStrategy(const CompareStrategy *compareStrategy);
std::vector<SXPeak> reduce(const std::vector<SXPeak> &peaks, Mantid::Kernel::ProgressBase &progress) const override;

private:
void reducePeaksFromNumberOfSpectras(std::vector<SXPeak> &inputPeaks) const;
};

class MANTID_CRYSTAL_DLL FindMaxReduceStrategy : public ReducePeakListStrategy {
Expand Down
86 changes: 80 additions & 6 deletions Framework/Crystal/src/FindSXPeaks.cpp
View file
Expand Up @@ -83,7 +83,8 @@ void FindSXPeaks::init() {
"spectrum. This is slower than StrongestPeakOnly. Note that the "
"recommended ResolutionStrategy in this mode is AbsoluteResolution.\n"
"3. AllPeaksNSigma: This stratergy will look for peaks by bins that are"
" more than nsigma different in intensity.");
" more than nsigma different in intensity. Note that the "
"recommended ResolutionStrategy in this mode is AbsoluteResolution.\n");

// Declare
declareProperty("SignalBackground", 10.0, mustBePositiveDouble,
Expand Down Expand Up @@ -178,10 +179,6 @@ void FindSXPeaks::init() {
"ResolutionStrategy", Mantid::Kernel::ePropertyCriterion::IS_EQUAL_TO,
absoluteResolutionPeaksStrategy));

declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
"The name of the PeaksWorkspace in which to store the list "
"of peaks found");

// Group
const std::string resolutionGroup = "Resolution Settings";
setPropertyGroup("ResolutionStrategy", resolutionGroup);
Expand All @@ -190,6 +187,26 @@ void FindSXPeaks::init() {
setPropertyGroup("PhiResolution", resolutionGroup);
setPropertyGroup("TwoThetaResolution", resolutionGroup);

// Declare
declareProperty("MinNBinsPerPeak", EMPTY_INT(), mustBePositive,
"Minimum number of bins contributing to a peak in an individual spectrum");

declareProperty("MinNSpectraPerPeak", EMPTY_INT(), mustBePositive,
"Minimum number of spectra contributing to a peak after they are grouped");

declareProperty("MaxNSpectraPerPeak", EMPTY_INT(), mustBePositive,
"Maximum number of spectra contributing to a peak after they are grouped");

// Group
const std::string peakValidationGroup = "Peak Validation Settings";
setPropertyGroup("MinNBinsPerPeak", peakValidationGroup);
setPropertyGroup("MinNSpectraPerPeak", peakValidationGroup);
setPropertyGroup("MaxNSpectraPerPeak", peakValidationGroup);

declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
"The name of the PeaksWorkspace in which to store the list "
"of peaks found");

// Create the output peaks workspace
m_peaks.reset(new PeaksWorkspace);
}
Expand All @@ -211,6 +228,47 @@ std::map<std::string, std::string> FindSXPeaks::validateInputs() {
validationOutput["XResolution"] = "XResolution must be set to a value greater than 0";
}

const int minNSpectraPerPeak = getProperty("MinNSpectraPerPeak");
const int maxNSpectraPerPeak = getProperty("MaxNSpectraPerPeak");
if (!isEmpty(minNSpectraPerPeak) && !isEmpty(maxNSpectraPerPeak)) {
if (maxNSpectraPerPeak < minNSpectraPerPeak) {
validationOutput["MaxNSpectraPerPeak"] = "MaxNSpectraPerPeak must be greater than MinNSpectraPerPeak";
validationOutput["MinNSpectraPerPeak"] = "MinNSpectraPerPeak must be lower than MaxNSpectraPerPeak";
}
}

MatrixWorkspace_const_sptr inputWorkspace = getProperty("InputWorkspace");
if (inputWorkspace) {
const int minWsIndex = getProperty("StartWorkspaceIndex");
const int maxWsIndex = getProperty("EndWorkspaceIndex");
size_t numberOfSpectraToConsider =
!isEmpty(minWsIndex) ? (!isEmpty(maxWsIndex) ? (maxWsIndex - minWsIndex + 1)
: (inputWorkspace->getNumberHistograms() - minWsIndex))
: (!isEmpty(maxWsIndex) ? (maxWsIndex + 1) : (inputWorkspace->getNumberHistograms()));

if (!isEmpty(minNSpectraPerPeak)) {
if (numberOfSpectraToConsider < minNSpectraPerPeak) {
validationOutput["MinNSpectraPerPeak"] =
"MinNSpectraPerPeak must be less than the number of spectrums considered in InputWorkspace";
}
}

if (!isEmpty(maxNSpectraPerPeak)) {
if (numberOfSpectraToConsider < maxNSpectraPerPeak) {
validationOutput["MaxNSpectraPerPeak"] =
"MaxNSpectraPerPeak must be less than the number of spectrums considered in InputWorkspace";
}
}

const int minNBinsPerPeak = getProperty("MinNBinsPerPeak");
if (!isEmpty(minNBinsPerPeak)) {
if (minNBinsPerPeak > inputWorkspace->getMaxNumberBins()) {
validationOutput["MinNBinsPerPeak"] =
"MinNBinsPerPeak must be less than the number of bins in the InputWorkspace";
}
}
}

return validationOutput;
}

Expand Down Expand Up @@ -248,7 +306,7 @@ void FindSXPeaks::exec() {
m_MaxWsIndex = numberOfSpectra - 1;
if (m_MaxWsIndex > numberOfSpectra - 1 || m_MaxWsIndex < m_MinWsIndex) {
g_log.warning("EndSpectrum out of range! Set to max detector number");
m_MaxWsIndex = numberOfSpectra;
m_MaxWsIndex = numberOfSpectra - 1;
}
if (m_MinRange > m_MaxRange) {
g_log.warning("Range_upper is less than Range_lower. Will integrate up to "
Expand All @@ -269,6 +327,11 @@ void FindSXPeaks::exec() {
auto peakFindingStrategy =
getPeakFindingStrategy(backgroundStrategy.get(), spectrumInfo, m_MinRange, m_MaxRange, xUnit);

const int minNBinsPerPeak = getProperty("MinNBinsPerPeak");
if (!isEmpty(minNBinsPerPeak)) {
peakFindingStrategy->setMinNBinsPerPeak(minNBinsPerPeak);
}

peakvector entries;
entries.reserve(m_MaxWsIndex - m_MinWsIndex);
// Count the peaks so that we can resize the peak vector at the end.
Expand Down Expand Up @@ -318,6 +381,17 @@ void FindSXPeaks::reducePeakList(const peakvector &pcv, Progress &progress) {
auto &goniometerMatrix = localworkspace->run().getGoniometer().getR();
auto compareStrategy = getCompareStrategy();
auto reductionStrategy = getReducePeakListStrategy(compareStrategy.get());

const int minNSpectraPerPeak = getProperty("MinNSpectraPerPeak");
if (!isEmpty(minNSpectraPerPeak)) {
reductionStrategy->setMinNSpectraPerPeak(minNSpectraPerPeak);
}

const int maxNSpectraPerPeak = getProperty("MaxNSpectraPerPeak");
if (!isEmpty(maxNSpectraPerPeak)) {
reductionStrategy->setMaxNSpectraPerPeak(maxNSpectraPerPeak);
}

auto finalv = reductionStrategy->reduce(pcv, progress);

for (auto &finalPeak : finalv) {
Expand Down
69 changes: 67 additions & 2 deletions Framework/Crystal/src/FindSXPeaksHelper.cpp
View file
Expand Up @@ -201,6 +201,11 @@ Getter for the detector Id.
*/
detid_t SXPeak::getDetectorId() const { return m_detId; }

/**
Getter for spectrum indexes of a peak
*/
const std::vector<int> &SXPeak::getPeakSpectras() const { return m_spectra; }

PeakContainer::PeakContainer(const HistogramData::HistogramY &y)
: m_y(y), m_startIndex(0), m_stopIndex(m_y.size() - 1), m_maxIndex(0) {}

Expand Down Expand Up @@ -233,8 +238,10 @@ size_t PeakContainer::getNumberOfPointsInPeak() const {
if (m_startIndex >= m_y.size()) {
return 0;
}

return m_stopIndex - m_startIndex;
if (m_stopIndex >= m_startIndex) {
return m_stopIndex - m_startIndex + 1;
}
return 0;
}

yIt PeakContainer::getMaxIterator() const { return m_y.begin() + m_maxIndex; }
Expand Down Expand Up @@ -296,6 +303,22 @@ PeakList PeakFindingStrategy::findSXPeaks(const HistogramData::HistogramX &x, co
return dofindSXPeaks(x, y, e, lowit, highit, workspaceIndex);
}

void PeakFindingStrategy::setMinNBinsPerPeak(int minNBinsPerPeak) { m_minNBinsPerPeak = minNBinsPerPeak; }

void PeakFindingStrategy::filterPeaksForMinBins(std::vector<std::unique_ptr<PeakContainer>> &inputPeakList) const {
if (m_minNBinsPerPeak == EMPTY_INT() || inputPeakList.empty()) {
return;
}

for (auto inputIt = inputPeakList.begin(); inputIt != inputPeakList.end();) {
if ((*inputIt)->getNumberOfPointsInPeak() < m_minNBinsPerPeak) {
inputIt = inputPeakList.erase(inputIt);
} else {
++inputIt;
}
}
}

BoundsIterator PeakFindingStrategy::getBounds(const HistogramData::HistogramX &x) const {
// Find the range [min,max]
auto lowit = (m_minValue == EMPTY_DBL()) ? x.begin() : std::lower_bound(x.begin(), x.end(), m_minValue);
Expand Down Expand Up @@ -435,6 +458,9 @@ PeakList AllPeaksStrategy::dofindSXPeaks(const HistogramData::HistogramX &x, con
// Get all peaks from the container
auto foundPeaks = getAllPeaks(x, y, low, high, m_backgroundStrategy);

// Filter the found peaks having the mininum number of bins
filterPeaksForMinBins(foundPeaks);

// Convert the found peaks to SXPeaks
auto peaks = convertToSXPeaks(x, y, foundPeaks, workspaceIndex);

Expand Down Expand Up @@ -514,6 +540,10 @@ PeakList NSigmaPeaksStrategy::dofindSXPeaks(const HistogramData::HistogramX &x,
const HistogramData::HistogramE &e, Bound low, Bound high,
const int workspaceIndex) const {
auto nsigmaPeaks = getAllNSigmaPeaks(x, y, e, low, high);

// Filter the found peaks having the mininum number of bins
filterPeaksForMinBins(nsigmaPeaks);

auto sxPeaks = convertToSXPeaks(x, y, nsigmaPeaks, workspaceIndex);
return sxPeaks;
}
Expand Down Expand Up @@ -596,6 +626,14 @@ std::vector<std::unique_ptr<PeakContainer>> NSigmaPeaksStrategy::getAllNSigmaPea
ReducePeakListStrategy::ReducePeakListStrategy(const CompareStrategy *compareStrategy)
: m_compareStrategy(compareStrategy) {}

void ReducePeakListStrategy::setMinNSpectraPerPeak(int minNSpectraPerPeak) {
m_minNSpectraPerPeak = minNSpectraPerPeak;
}

void ReducePeakListStrategy::setMaxNSpectraPerPeak(int maxSpectrasForPeak) {
m_maxNSpectraPerPeak = maxSpectrasForPeak;
}

SimpleReduceStrategy::SimpleReduceStrategy(const CompareStrategy *compareStrategy)
: ReducePeakListStrategy(compareStrategy) {}

Expand All @@ -621,9 +659,26 @@ std::vector<SXPeak> SimpleReduceStrategy::reduce(const std::vector<SXPeak> &peak
}
}

reducePeaksFromNumberOfSpectras(finalPeaks);

return finalPeaks;
}

void SimpleReduceStrategy::reducePeaksFromNumberOfSpectras(std::vector<SXPeak> &inputPeaks) const {
if (m_minNSpectraPerPeak == EMPTY_INT() && m_maxNSpectraPerPeak == EMPTY_INT()) {
return;
}

for (auto peakIt = inputPeaks.begin(); peakIt != inputPeaks.end();) {
if (((m_minNSpectraPerPeak != EMPTY_INT()) && ((*peakIt).getPeakSpectras().size() < m_minNSpectraPerPeak)) ||
((m_maxNSpectraPerPeak != EMPTY_INT()) && ((*peakIt).getPeakSpectras().size() > m_maxNSpectraPerPeak))) {
peakIt = inputPeaks.erase(peakIt);
} else {
++peakIt;
}
}
}

FindMaxReduceStrategy::FindMaxReduceStrategy(const CompareStrategy *compareStrategy)
: ReducePeakListStrategy(compareStrategy) {}

Expand Down Expand Up @@ -733,8 +788,17 @@ std::vector<SXPeak> FindMaxReduceStrategy::getFinalPeaks(const std::vector<std::
// Currently we select the peak with the largest signal (this strategy could
// be changed to something like a weighted mean or similar)
for (const auto &group : peakGroups) {
// When MinNSpectraPerPeak or maxNSpectraPerPeak parameters are provided,
// a group will be ignored if it does not satisfy the minimum or the maximum number of spectrums
// required to identify as a peak.
if ((m_minNSpectraPerPeak != EMPTY_INT() && group.size() < m_minNSpectraPerPeak) ||
(m_maxNSpectraPerPeak != EMPTY_INT() && group.size() > m_maxNSpectraPerPeak)) {
continue;
}

SXPeak *maxPeak = nullptr;
double maxIntensity = std::numeric_limits<double>::min();

for (auto *element : group) {
if (element->getIntensity() > maxIntensity) {
maxIntensity = element->getIntensity();
Expand All @@ -748,6 +812,7 @@ std::vector<SXPeak> FindMaxReduceStrategy::getFinalPeaks(const std::vector<std::
peaks.emplace_back(*maxPeak);
}
}

return peaks;
}

Expand Down

0 comments on commit 4abe275

Please sign in to comment.