40const std::string NON_UNIFORM_GROUP =
"NonUniform Detectors";
41const std::string RECTANGULAR_GROUP =
"Rectangular Detectors";
42const std::string INPUT_WORKSPACE =
"InputWorkspace";
47 using Callable = std::function<void()>;
48 CallOnExit(Callable &&callable) noexcept :
m_callable(std::move(callable)) {}
65 std::make_shared<InstrumentValidator>()),
66 "The workspace containing the spectra to be averaged.");
68 "The name of the workspace to be created as the output of "
72 auto mustBePositiveDouble = std::make_shared<BoundedValidator<double>>();
73 mustBePositiveDouble->setLower(0.0);
76 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
77 mustBePositive->setLower(0);
79 std::vector<std::string> propOptions{
"Flat",
"Linear",
"Parabolic",
"Gaussian"};
80 declareProperty(
"WeightedSum",
"Flat", std::make_shared<StringListValidator>(propOptions),
81 "What sort of Weighting scheme to use?\n"
82 " Flat: Effectively no-weighting, all weights are 1.\n"
83 " Linear: Linear weighting 1 - r/R from origin.\n"
84 " Parabolic : Weighting as cutoff - x + cutoff - y + 1."
85 " Gaussian : Uses the absolute distance x^2 + y^2 ... "
86 "normalised by the cutoff^2");
88 declareProperty(
"Sigma", 0.5, mustBePositiveDouble,
"Sigma value for gaussian weighting schemes. Defaults to 0.5. ");
91 declareProperty(
"IgnoreMaskedDetectors",
true,
"If true, do not consider masked detectors in the NN search.");
94 "If the InputWorkspace is an "
95 "EventWorkspace, this will preserve "
96 "the full event list (warning: this "
97 "will use much more memory!).");
102 "The number of X (horizontal) adjacent pixels to average together. "
103 "Only for instruments with RectangularDetectors. ");
106 "The number of Y (vertical) adjacent pixels to average together. "
107 "Only for instruments with RectangularDetectors. ");
110 "The total number of X (horizontal) adjacent pixels to sum together. "
111 "Only for instruments with RectangularDetectors. AdjX will be ignored "
112 "if SumPixelsX > 1.");
115 "The total number of Y (vertical) adjacent pixels to sum together. "
116 "Only for instruments with RectangularDetectors. AdjY will be ignored if "
120 "The number of pixels to zero at edges. "
121 "Only for instruments with RectangularDetectors. ");
131 std::vector<std::string> radiusPropOptions{
"Meters",
"NumberOfPixels"};
132 declareProperty(
"RadiusUnits",
"Meters", std::make_shared<StringListValidator>(radiusPropOptions),
133 "Units used to specify the radius.\n"
134 " Meters : Radius is in meters.\n"
135 " NumberOfPixels : Radius is in terms of the number of pixels.");
138 "The radius cut-off around a pixel to look for nearest neighbours to "
140 "This radius cut-off is applied to a set of nearest neighbours whose "
142 "defined in the NumberOfNeighbours property. See below for more details. "
144 "If 0, will use the AdjX and AdjY parameters for rectangular detectors "
148 "Number of nearest neighbouring pixels.\n"
149 "The default is 8.");
152 "Sum nearest neighbouring pixels with same parent.\n"
153 "Number of pixels will be reduced. The default is false.");
156 "OuputWorkspace will have same number of pixels as "
157 "InputWorkspace using SumPixelsX and SumPixelsY. Individual "
158 "pixels will have averages.");
171 g_log.
debug(
"SmoothNeighbours processing assuming rectangular detectors.");
180 Progress prog(
this, 0.0, 1.0, inst->nelements());
183 std::vector<std::shared_ptr<RectangularDetector>> detList;
184 for (
int i = 0; i < inst->nelements(); i++) {
185 std::shared_ptr<RectangularDetector> det;
186 std::shared_ptr<ICompAssembly> assem;
187 std::shared_ptr<ICompAssembly> assem2;
189 det = std::dynamic_pointer_cast<RectangularDetector>((*inst)[i]);
191 detList.emplace_back(det);
196 assem = std::dynamic_pointer_cast<ICompAssembly>((*inst)[i]);
198 for (
int j = 0; j < assem->nelements(); j++) {
199 det = std::dynamic_pointer_cast<RectangularDetector>((*assem)[j]);
201 detList.emplace_back(det);
208 assem2 = std::dynamic_pointer_cast<ICompAssembly>((*assem)[j]);
210 for (
int k = 0; k < assem2->nelements(); k++) {
211 det = std::dynamic_pointer_cast<RectangularDetector>((*assem2)[k]);
213 detList.emplace_back(det);
223 if (detList.empty()) {
240 bool sum = sumX * sumY > 1;
250 std::vector<std::pair<int, int>> idToIndexMap;
251 idToIndexMap.reserve(detList.size());
252 for (
int i = 0; i < static_cast<int>(detList.size()); i++)
253 idToIndexMap.emplace_back(detList[i]->getAtXY(0, 0)->getID(), i);
257 stable_sort(idToIndexMap.begin(), idToIndexMap.end());
260 for (
const auto &idIndex : idToIndexMap) {
261 std::shared_ptr<RectangularDetector> det = detList[idIndex.second];
262 const std::string det_name = det->getName();
263 for (
int j = 0; j < det->xpixels(); j += sumX) {
264 for (
int k = 0; k < det->ypixels(); k += sumY) {
265 double totalWeight = 0;
267 std::vector<weightedNeighbour> neighbours;
269 for (
int ix = startX; ix <= endX; ix++)
270 for (
int iy = startY; iy <= endY; iy++) {
276 if (j + ix >= det->xpixels() -
m_edge || j + ix <
m_edge)
278 if (k + iy >= det->ypixels() -
m_edge || k + iy <
m_edge)
280 int pixelID = det->getAtXY(j + ix, k + iy)->getID();
283 auto mapEntry = pixel_to_wi.find(pixelID);
284 if (mapEntry != pixel_to_wi.end()) {
285 size_t wi = mapEntry->second;
286 neighbours.emplace_back(wi, smweight);
288 totalWeight += smweight;
294 for (
auto &neighbour : neighbours)
295 neighbour.second /= totalWeight;
312 g_log.
debug(
"SmoothNeighbours processing NOT assuming rectangular detectors.");
315 this->
progress(0.2,
"Building Neighbour Map");
323 bool ignoreMaskedDetectors =
getProperty(
"IgnoreMaskedDetectors");
332 std::shared_ptr<const Geometry::IComponent> parent, neighbParent, grandparent, neighbGParent;
333 std::vector<bool> used(
m_inWS->getNumberHistograms(),
false);
334 const auto &detectorInfo =
m_inWS->detectorInfo();
335 for (
size_t wi = 0; wi <
m_inWS->getNumberHistograms(); wi++) {
342 const auto &dets =
m_inWS->getSpectrum(wi).getDetectorIDs();
343 const auto index = detectorInfo.indexOf(*dets.begin());
344 if (detectorInfo.isMonitor(
index))
346 if (detectorInfo.isMasked(
index)) {
354 const auto &det = detectorInfo.detector(
index);
355 parent = det.getParent();
357 grandparent = parent->getParent();
374 neighbSpectra[inSpec] =
V3D(0.0, 0.0, 0.0);
377 double totalWeight = 0;
379 std::vector<weightedNeighbour> neighbours;
382 for (
auto &specDistance : neighbSpectra) {
386 double weight =
m_weightedSum->weightAt(specDistance.second);
390 auto mapIt = spec2index.find(spec);
391 if (mapIt != spec2index.end()) {
392 size_t neighWI = mapIt->second;
395 const std::set<detid_t> &dets =
m_inWS->getSpectrum(neighWI).getDetectorIDs();
396 const auto &det = detectorInfo.detector(*dets.begin());
397 neighbParent = det.getParent();
398 neighbGParent = neighbParent->getParent();
399 if (noNeigh >= sum || neighbParent->getName() != parent->getName() ||
400 neighbGParent->getName() != grandparent->getName() || used[neighWI])
403 used[neighWI] =
true;
405 neighbours.emplace_back(neighWI, weight);
406 totalWeight += weight;
413 for (
auto &neighbour : neighbours)
414 neighbour.second /= totalWeight;
432 if (strategyName ==
"Flat") {
434 }
else if (strategyName ==
"Linear") {
436 }
else if (strategyName ==
"Parabolic") {
437 m_weightedSum = std::make_unique<ParabolicWeighting>(cutOff);
438 }
else if (strategyName ==
"Gaussian") {
451 double translatedRadius = 0;
452 if (radiusUnits ==
"Meters") {
454 translatedRadius = enteredRadius;
455 }
else if (radiusUnits ==
"NumberOfPixels") {
457 const auto &firstDet =
m_inWS->spectrumInfo().detector(0);
460 firstDet.getBoundingBox(bbox);
463 translatedRadius = bbox.
width().
norm() * enteredRadius;
465 const std::string message =
"SmoothNeighbours::translateToMeters, Unknown Unit: " + radiusUnits;
466 throw std::invalid_argument(message);
468 return translatedRadius;
492 m_progress = std::make_unique<Progress>(
this, 0.0, 0.2,
m_inWS->getNumberHistograms());
495 CallOnExit resetInWSOnExit([
this]() {
m_inWS.reset(); });
496 CallOnExit resetNeighboursOnExit([
this]() {
507 auto wsEvent = std::dynamic_pointer_cast<EventWorkspace>(
m_inWS);
522 const size_t numberOfSpectra =
m_outWI;
524 const size_t YLength =
m_inWS->blocksize();
528 if (std::dynamic_pointer_cast<OffsetsWorkspace>(
m_inWS)) {
532 outWS = std::dynamic_pointer_cast<MatrixWorkspace>(
541 for (
int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
544 auto &outSpec = outWS->getSpectrum(outWIi);
547 auto &outY = outSpec.mutableY();
549 auto &outE = outSpec.mutableE();
551 auto &outX = outSpec.mutableX();
554 std::vector<weightedNeighbour> &neighbours =
m_neighbours[outWIi];
555 std::vector<weightedNeighbour>::iterator it;
556 for (it = neighbours.begin(); it != neighbours.end(); ++it) {
557 size_t inWI = it->first;
558 double weight = it->second;
559 double weightSquared = weight * weight;
561 const auto &inSpec =
m_inWS->getSpectrum(inWI);
562 const auto &inY = inSpec.y();
563 const auto &inE = inSpec.e();
564 const auto &inX = inSpec.x();
566 for (
size_t i = 0; i < YLength; i++) {
568 outY[i] += inY[i] * weight;
571 double errorSquared = inE[i];
572 errorSquared *= errorSquared;
573 errorSquared *= weightSquared;
574 outE[i] += errorSquared;
578 if (
m_inWS->isHistogramData()) {
579 outX[YLength] = inX[YLength];
584 for (
size_t i = 0; i < YLength; i++)
585 outE[i] = sqrt(outE[i]);
605 for (
int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
613 const auto &inSpec =
m_inWS->getSpectrum(neighbor.first);
614 outSpec.addDetectorIDs(inSpec.getDetectorIDs());
623 const size_t numberOfSpectra =
m_inWS->getNumberHistograms();
625 const size_t YLength =
m_inWS->blocksize();
629 if (std::dynamic_pointer_cast<OffsetsWorkspace>(
m_inWS)) {
631 outws2 = std::make_shared<OffsetsWorkspace>(
m_inWS->getInstrument());
633 outws2 = std::dynamic_pointer_cast<MatrixWorkspace>(
640 for (
int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
641 const auto &inSpec =
m_inWS->getSpectrum(outWIi);
642 auto &outSpec2 = outws2->getSpectrum(outWIi);
643 outSpec2.mutableX() = inSpec.x();
644 outSpec2.addDetectorIDs(inSpec.getDetectorIDs());
646 outSpec2.clearData();
650 const size_t numberOfSpectra2 = outWS->getNumberHistograms();
651 for (
int outWIi = 0; outWIi < int(numberOfSpectra2); outWIi++) {
655 outws2->setHistogram(neighbor.first, outWS->histogram(outWIi));
668 const size_t numberOfSpectra =
m_outWI;
669 const auto YLength =
static_cast<int>(
m_inWS->blocksize());
673 outWS = std::dynamic_pointer_cast<EventWorkspace>(
680 this->
setProperty(
"OutputWorkspace", std::dynamic_pointer_cast<MatrixWorkspace>(outWS));
684 for (
int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
688 EventList &outEL = outWS->getSpectrum(outWIi);
691 std::vector<weightedNeighbour> &neighbours =
m_neighbours[outWIi];
692 std::vector<weightedNeighbour>::iterator it;
693 for (it = neighbours.begin(); it != neighbours.end(); ++it) {
694 size_t inWI = it->first;
696 double weight = it->second;
711 outWS->setAllX(
m_inWS->binEdges(0));
#define DECLARE_ALGORITHM(classname)
std::map< DeltaEMode::Type, std::string > index
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
const VecProperties ConstVecProperties
std::vector< Mantid::Kernel::Property * > VecProperties
Base class from which all concrete algorithm classes should be derived.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
void clearDetectorIDs()
Clear the detector IDs set.
Base MatrixWorkspace Abstract Class.
virtual ISpectrum & getSpectrum(const size_t index)=0
Return the underlying ISpectrum ptr at the given workspace index.
virtual std::size_t getNumberHistograms() const =0
Returns the number of histograms in the workspace.
Helper class for reporting progress from algorithms.
WorkspaceNearestNeighbourInfo provides easy access to nearest-neighbour information for a workspace.
std::map< specnum_t, Kernel::V3D > getNeighboursExact(specnum_t spec) const
Queries the WorkspaceNearestNeighbours object for the selected spectrum number.
A property class for workspaces.
SpectraDistanceMap apply(SpectraDistanceMap &unfiltered) const
Apply the filtering based on radius.
void execWorkspace2D()
Execute the algorithm for a Workspace2D/don't preserve events input.
std::vector< std::vector< weightedNeighbour > > m_neighbours
Vector of list of neighbours (with weight) for each workspace index.
void setWeightingStrategy(const std::string &strategyName, double &cutOff)
Sets the weighting stragegy.
int m_edge
Edge pixels to ignore.
void init() override
Virtual method - must be overridden by concrete algorithm.
std::unique_ptr< WeightingStrategy > m_weightedSum
Weight the neighbours during summing.
void findNeighboursUbiquitous()
Use NearestNeighbours to find the neighbours for any instrument.
size_t m_outWI
number of output workspace pixels
bool m_expandSumAllPixels
expand by pixel IDs
double m_radius
Radius to search nearest neighbours.
SmoothNeighbours()
Default constructor.
void spreadPixels(const API::MatrixWorkspace_sptr &outWS)
Build the instrument/detector setup in workspace.
void setupNewInstrument(API::MatrixWorkspace &outWS) const
Build the instrument/detector setup in workspace.
double translateToMeters(const std::string &radiusUnits, const double &enteredRadius) const
Translate the entered radius into meters.
Mantid::API::MatrixWorkspace_sptr m_inWS
Input workspace.
void exec() override
Executes the algorithm.
std::unique_ptr< Mantid::API::Progress > m_progress
Progress reporter.
void findNeighboursRectangular()
Fill the neighbours list given the AdjX AdjY parameters and an instrument with rectangular detectors.
void execEvent(Mantid::DataObjects::EventWorkspace_sptr &ws)
Execute the algorithm for a EventWorkspace input.
int m_nNeighbours
Number of neighbours.
bool m_preserveEvents
PreserveEvents.
An OffsetsWorkspace is a specialized Workspace2D where the Y value at each pixel is the offset to be ...
A simple structure that defines an axis-aligned cuboid shaped bounding box for a geometrical object.
Kernel::V3D width() const
Returns the width of the box.
Exception for when an item is not found in a collection.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
void debug(const std::string &msg)
Logs at debug level.
void information(const std::string &msg)
Logs at information level.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
double norm() const noexcept
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::map< specnum_t, Mantid::Kernel::V3D > SpectraDistanceMap
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
std::unordered_map< specnum_t, size_t > spec2index_map
Map with key = spectrum number, value = workspace index.
std::unordered_map< detid_t, size_t > detid2index_map
Map with key = detector ID, value = workspace index.
int32_t specnum_t
Typedef for a spectrum Number.
@ Input
An input workspace.
@ Output
An output workspace.