diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SeparateBackgroundFromSignal.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SeparateBackgroundFromSignal.h index 9c58cf46cf2f06c6215f49c06959ba4c597cc9a7..bf04c9efb6823ac1fb5c30983f9f8a933dd8eeea 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SeparateBackgroundFromSignal.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/SeparateBackgroundFromSignal.h @@ -54,6 +54,19 @@ namespace Algorithms /// Implement abstract Algorithm methods void exec(); double moment(MantidVec& X, size_t n, double mean, int k); + struct cont_peak + { + size_t start; + size_t stop; + double maxY; + }; + struct by_len + { + bool operator()(cont_peak const &a, cont_peak const &b) + { + return a.maxY > b.maxY; + } + }; }; diff --git a/Code/Mantid/Framework/Algorithms/src/SeparateBackgroundFromSignal.cpp b/Code/Mantid/Framework/Algorithms/src/SeparateBackgroundFromSignal.cpp index c39e697b60f2459086de8ba1a5a5b60512bda9ca..ed70b108fbfb1f56c13156fd6f4a7b9e2435866e 100644 --- a/Code/Mantid/Framework/Algorithms/src/SeparateBackgroundFromSignal.cpp +++ b/Code/Mantid/Framework/Algorithms/src/SeparateBackgroundFromSignal.cpp @@ -147,54 +147,58 @@ namespace Algorithms } while (std::abs(Ymean-Yvariance) > k * Ysigma); - // remove single outliers - - if (mask[1] == mask[2] && mask[2] == mask[3]) - mask[0] = mask[1]; - if (mask[0] == mask[2] && mask[2] == mask[3]) - mask[1] = mask[2]; - for (size_t l = 2; l < n-2; ++l) - { - if (mask[l-1] == mask[l+1] && (mask[l-1] == mask[l-2] || mask[l+1] == mask[l+2])) - { - mask[l] = mask[l+1]; - } - } - if (mask[n-2] == mask[n-3] && mask[n-3] == mask[n-4]) - mask[n-1] = mask[n-2]; - if (mask[n-1] == mask[n-3] && mask[n-3] == mask[n-4]) - mask[n-2] = mask[n-1]; - - // mask regions not connected to largest region - vector<size_t> cont_start; - vector<size_t> cont_stop; - for (size_t l = 1; l < n-1; ++l) - { - if (mask[l] != mask[l-1] && mask[l] == 1) + if(n > 5) + { + // remove single outliers + if (mask[1] == mask[2] && mask[2] == mask[3]) + mask[0] = mask[1]; + if (mask[0] == mask[2] && mask[2] == mask[3]) + mask[1] = mask[2]; + for (size_t l = 2; l < n-3; ++l) { - cont_start.push_back(l); + if (mask[l-1] == mask[l+1] && (mask[l-1] == mask[l-2] || mask[l+1] == mask[l+2])) + { + mask[l] = mask[l+1]; + } } - if (mask[l] != mask[l-1] && mask[l] == 0) + if (mask[n-2] == mask[n-3] && mask[n-3] == mask[n-4]) + mask[n-1] = mask[n-2]; + if (mask[n-1] == mask[n-3] && mask[n-3] == mask[n-4]) + mask[n-2] = mask[n-1]; + + // mask regions not connected to largest region + // for loop can start > 1 for multiple peaks + vector<cont_peak> peaks; + for (size_t l = 1; l < n; ++l) { - cont_stop.push_back(l-1); + if (mask[l] != mask[l-1] && mask[l] == 1) + { + peaks.push_back(cont_peak()); + peaks[peaks.size()-1].start = l; + } + if (peaks.size() > 0) + { + size_t ipeak = peaks.size()-1; + if (mask[l] != mask[l-1] && mask[l] == 0) + { + peaks[ipeak].stop = l-1; + } + if (inpY[l] > peaks[ipeak].maxY) peaks[ipeak].maxY = inpY[l]; + } } - } - if(cont_start.size() > cont_stop.size()) cont_stop.push_back(n-1); - vector<size_t> cont_len; - for (size_t l = 0; l < cont_start.size(); ++l) - { - cont_len.push_back(cont_stop[l] - cont_start[l] + 1); - } - std::vector<size_t>::iterator ic = std::max_element(cont_len.begin(), cont_len.end()); - const size_t c = ic - cont_len.begin(); - for (size_t l = 0; l < cont_len.size(); ++l) - { - if (l != c) for (size_t j = cont_start[l]; j <= cont_stop[l]; ++j) + if(peaks.size()> 0) { - mask[j] = 0; + if(peaks[peaks.size()-1].stop == 0) peaks[peaks.size()-1].stop = n-1; + std::sort(peaks.begin(), peaks.end(), by_len()); + for (size_t l = 1; l < peaks.size(); ++l) + { + for (size_t j = peaks[l].start; j <= peaks[l].stop; ++j) + { + mask[j] = 0; + } + } } - } - + } // save output of mask * Y vecX = inpX; std::transform(mask.begin(), mask.end(), inpY.begin(), vecY.begin(), std::multiplies<double>());