Mantid
Loading...
Searching...
No Matches
SmoothNeighbours.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
23
24using namespace Mantid::Kernel;
25using namespace Mantid::Geometry;
26using namespace Mantid::API;
27using namespace Mantid::DataObjects;
28
29using VecProperties = std::vector<Mantid::Kernel::Property *>;
31
32namespace Mantid::Algorithms {
33
34// Register the class into the algorithm factory
35DECLARE_ALGORITHM(SmoothNeighbours)
36
37namespace {
38// Used in custom GUI. Make sure you change them in SmoothNeighboursDialog.cpp
39// as well.
40const std::string NON_UNIFORM_GROUP = "NonUniform Detectors";
41const std::string RECTANGULAR_GROUP = "Rectangular Detectors";
42const std::string INPUT_WORKSPACE = "InputWorkspace";
43
46struct CallOnExit {
47 using Callable = std::function<void()>;
48 CallOnExit(Callable &&callable) noexcept : m_callable(std::move(callable)) {}
49 ~CallOnExit() {
50 try {
51 m_callable();
52 } catch (...) {
53 };
54 }
55
56private:
57 Callable m_callable;
58};
59} // namespace
60
61SmoothNeighbours::SmoothNeighbours() : API::Algorithm(), m_weightedSum(std::make_unique<NullWeighting>()) {}
62
64 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>(INPUT_WORKSPACE, "", Direction::Input,
65 std::make_shared<InstrumentValidator>()),
66 "The workspace containing the spectra to be averaged.");
67 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
68 "The name of the workspace to be created as the output of "
69 "the algorithm.");
70
71 // Unsigned double
72 auto mustBePositiveDouble = std::make_shared<BoundedValidator<double>>();
73 mustBePositiveDouble->setLower(0.0);
74
75 // Unsigned int.
76 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
77 mustBePositive->setLower(0);
78
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");
87
88 declareProperty("Sigma", 0.5, mustBePositiveDouble, "Sigma value for gaussian weighting schemes. Defaults to 0.5. ");
89 setPropertySettings("Sigma", std::make_unique<EnabledWhenProperty>("WeightedSum", IS_EQUAL_TO, "Gaussian"));
90
91 declareProperty("IgnoreMaskedDetectors", true, "If true, do not consider masked detectors in the NN search.");
92
93 declareProperty("PreserveEvents", true,
94 "If the InputWorkspace is an "
95 "EventWorkspace, this will preserve "
96 "the full event list (warning: this "
97 "will use much more memory!).");
98
99 // -- Rectangular properties --
100
101 declareProperty("AdjX", 1, mustBePositive,
102 "The number of X (horizontal) adjacent pixels to average together. "
103 "Only for instruments with RectangularDetectors. ");
104
105 declareProperty("AdjY", 1, mustBePositive,
106 "The number of Y (vertical) adjacent pixels to average together. "
107 "Only for instruments with RectangularDetectors. ");
108
109 declareProperty("SumPixelsX", 1, mustBePositive,
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.");
113
114 declareProperty("SumPixelsY", 1, mustBePositive,
115 "The total number of Y (vertical) adjacent pixels to sum together. "
116 "Only for instruments with RectangularDetectors. AdjY will be ignored if "
117 "SumPixelsY > 1");
118
119 declareProperty("ZeroEdgePixels", 0, mustBePositive,
120 "The number of pixels to zero at edges. "
121 "Only for instruments with RectangularDetectors. ");
122
123 setPropertyGroup("AdjX", RECTANGULAR_GROUP);
124 setPropertyGroup("AdjY", RECTANGULAR_GROUP);
125 setPropertyGroup("SumPixelsX", RECTANGULAR_GROUP);
126 setPropertyGroup("SumPixelsY", RECTANGULAR_GROUP);
127 setPropertyGroup("ZeroEdgePixels", RECTANGULAR_GROUP);
128
129 // -- Non-uniform properties --
130
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.");
136
137 declareProperty("Radius", 0.0, mustBePositiveDouble,
138 "The radius cut-off around a pixel to look for nearest neighbours to "
139 "average. \n"
140 "This radius cut-off is applied to a set of nearest neighbours whose "
141 "number is "
142 "defined in the NumberOfNeighbours property. See below for more details. "
143 "\n"
144 "If 0, will use the AdjX and AdjY parameters for rectangular detectors "
145 "instead.");
146
147 declareProperty("NumberOfNeighbours", 8, mustBePositive,
148 "Number of nearest neighbouring pixels.\n"
149 "The default is 8.");
150
151 declareProperty("SumNumberOfNeighbours", 1,
152 "Sum nearest neighbouring pixels with same parent.\n"
153 "Number of pixels will be reduced. The default is false.");
154
155 declareProperty("ExpandSumAllPixels", false,
156 "OuputWorkspace will have same number of pixels as "
157 "InputWorkspace using SumPixelsX and SumPixelsY. Individual "
158 "pixels will have averages.");
159
160 setPropertyGroup("RadiusUnits", NON_UNIFORM_GROUP);
161 setPropertyGroup("Radius", NON_UNIFORM_GROUP);
162 setPropertyGroup("NumberOfNeighbours", NON_UNIFORM_GROUP);
163 setPropertyGroup("SumNumberOfNeighbours", NON_UNIFORM_GROUP);
164}
165
166//--------------------------------------------------------------------------------------------
171 g_log.debug("SmoothNeighbours processing assuming rectangular detectors.");
172
173 m_progress->resetNumSteps(m_inWS->getNumberHistograms(), 0.2, 0.5);
174
175 Instrument_const_sptr inst = m_inWS->getInstrument();
176
177 // To get the workspace index from the detector ID
178 const detid2index_map pixel_to_wi = m_inWS->getDetectorIDToWorkspaceIndexMap(true, true);
179
180 Progress prog(this, 0.0, 1.0, inst->nelements());
181
182 // Build a list of Rectangular Detectors
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;
188
189 det = std::dynamic_pointer_cast<RectangularDetector>((*inst)[i]);
190 if (det) {
191 detList.emplace_back(det);
192 } else {
193 // Also, look in the first sub-level for RectangularDetectors (e.g. PG3).
194 // We are not doing a full recursive search since that will be very long
195 // for lots of pixels.
196 assem = std::dynamic_pointer_cast<ICompAssembly>((*inst)[i]);
197 if (assem) {
198 for (int j = 0; j < assem->nelements(); j++) {
199 det = std::dynamic_pointer_cast<RectangularDetector>((*assem)[j]);
200 if (det) {
201 detList.emplace_back(det);
202
203 } else {
204 // Also, look in the second sub-level for RectangularDetectors (e.g.
205 // PG3).
206 // We are not doing a full recursive search since that will be very
207 // long for lots of pixels.
208 assem2 = std::dynamic_pointer_cast<ICompAssembly>((*assem)[j]);
209 if (assem2) {
210 for (int k = 0; k < assem2->nelements(); k++) {
211 det = std::dynamic_pointer_cast<RectangularDetector>((*assem2)[k]);
212 if (det) {
213 detList.emplace_back(det);
214 }
215 }
216 }
217 }
218 }
219 }
220 }
221 }
222
223 if (detList.empty()) {
224 // Not rectangular so use Nearest Neighbours
225 m_radius = translateToMeters("NumberOfPixels", std::max(m_adjX, m_adjY));
229 return;
230 }
231
232 // Resize the vector we are setting
233 m_neighbours.resize(m_inWS->getNumberHistograms());
234 int startX = -m_adjX;
235 int startY = -m_adjY;
236 int endX = m_adjX;
237 int endY = m_adjY;
238 const int sumX = getProperty("SumPixelsX");
239 const int sumY = getProperty("SumPixelsY");
240 bool sum = sumX * sumY > 1;
241 if (sum) {
242 startX = 0;
243 startY = 0;
244 endX = sumX - 1;
245 endY = sumY - 1;
246 }
247
248 m_outWI = 0;
249 // Build a map to sort by the detectorID
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);
254
255 // To sort in descending order
256 if (sum)
257 stable_sort(idToIndexMap.begin(), idToIndexMap.end());
258
259 // Loop through the RectangularDetector's we listed before.
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;
266 // Neighbours and weights
267 std::vector<weightedNeighbour> neighbours;
268
269 for (int ix = startX; ix <= endX; ix++)
270 for (int iy = startY; iy <= endY; iy++) {
271 // Weights for corners=1; higher for center and adjacent pixels
272 double smweight = m_weightedSum->weightAt(m_adjX, ix, m_adjY, iy);
273
274 // Find the pixel ID at that XY position on the rectangular
275 // detector
276 if (j + ix >= det->xpixels() - m_edge || j + ix < m_edge)
277 continue;
278 if (k + iy >= det->ypixels() - m_edge || k + iy < m_edge)
279 continue;
280 int pixelID = det->getAtXY(j + ix, k + iy)->getID();
281
282 // Find the corresponding workspace index, if any
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);
287 // Count the total weight
288 totalWeight += smweight;
289 }
290 }
291
292 // Adjust the weights of each neighbour to normalize to unity
293 if (!sum || m_expandSumAllPixels)
294 for (auto &neighbour : neighbours)
295 neighbour.second /= totalWeight;
296
297 // Save the list of neighbours for this output workspace index.
298 m_neighbours[m_outWI] = neighbours;
299 m_outWI++;
300
301 m_progress->report("Finding Neighbours");
302 }
303 }
304 prog.report(det_name);
305 }
306}
307
308//--------------------------------------------------------------------------------------------
312 g_log.debug("SmoothNeighbours processing NOT assuming rectangular detectors.");
313
314 m_progress->resetNumSteps(m_inWS->getNumberHistograms(), 0.2, 0.5);
315 this->progress(0.2, "Building Neighbour Map");
316
317 Instrument_const_sptr inst = m_inWS->getInstrument();
318 const spec2index_map spec2index = m_inWS->getSpectrumToWorkspaceIndexMap();
319
320 // Resize the vector we are setting
321 m_neighbours.resize(m_inWS->getNumberHistograms());
322
323 bool ignoreMaskedDetectors = getProperty("IgnoreMaskedDetectors");
324 WorkspaceNearestNeighbourInfo neighbourInfo(*m_inWS, ignoreMaskedDetectors, m_nNeighbours);
325
326 // Cull by radius
327 RadiusFilter radiusFilter(m_radius);
328
329 // Go through every input workspace pixel
330 m_outWI = 0;
331 int sum = getProperty("SumNumberOfNeighbours");
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++) {
336 if (sum > 1)
337 if (used[wi])
338 continue;
339 // We want to skip monitors
340 try {
341 // Get the list of detectors in this pixel
342 const auto &dets = m_inWS->getSpectrum(wi).getDetectorIDs();
343 const auto index = detectorInfo.indexOf(*dets.begin());
344 if (detectorInfo.isMonitor(index))
345 continue; // skip monitor
346 if (detectorInfo.isMasked(index)) {
347 // Calibration masks many detectors, but there should be 0s after
348 // smoothing
349 if (sum == 1)
350 m_outWI++;
351 continue; // skip masked detectors
352 }
353 if (sum > 1) {
354 const auto &det = detectorInfo.detector(index);
355 parent = det.getParent();
356 if (parent)
357 grandparent = parent->getParent();
358 }
360 continue; // skip missing detector
361 }
362
363 specnum_t inSpec = m_inWS->getSpectrum(wi).getSpectrumNo();
364
365 // Step one - Get the number of specified neighbours
366 SpectraDistanceMap insideGrid = neighbourInfo.getNeighboursExact(inSpec);
367
368 // Step two - Filter the results by the radius cut off.
369 SpectraDistanceMap neighbSpectra = radiusFilter.apply(insideGrid);
370
371 // Force the central pixel to always be there
372 // There seems to be a bug in nearestNeighbours, returns distance != 0.0 for
373 // the central pixel. So we force distance = 0
374 neighbSpectra[inSpec] = V3D(0.0, 0.0, 0.0);
375
376 // Neighbours and weights list
377 double totalWeight = 0;
378 int noNeigh = 0;
379 std::vector<weightedNeighbour> neighbours;
380
381 // Convert from spectrum numbers to workspace indices
382 for (auto &specDistance : neighbSpectra) {
383 specnum_t spec = specDistance.first;
384
385 // Use the weighting strategy to calculate the weight.
386 double weight = m_weightedSum->weightAt(specDistance.second);
387
388 if (weight > 0) {
389 // Find the corresponding workspace index
390 auto mapIt = spec2index.find(spec);
391 if (mapIt != spec2index.end()) {
392 size_t neighWI = mapIt->second;
393 if (sum > 1) {
394 // Get the list of detectors in this pixel
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])
401 continue;
402 noNeigh++;
403 used[neighWI] = true;
404 }
405 neighbours.emplace_back(neighWI, weight);
406 totalWeight += weight;
407 }
408 }
409 }
410
411 // Adjust the weights of each neighbour to normalize to unity
412 if (sum == 1)
413 for (auto &neighbour : neighbours)
414 neighbour.second /= totalWeight;
415
416 // Save the list of neighbours for this output workspace index.
417 m_neighbours[m_outWI] = neighbours;
418 m_outWI++;
419
420 m_progress->report("Finding Neighbours");
421 } // each workspace index
422}
423
431void SmoothNeighbours::setWeightingStrategy(const std::string &strategyName, double &cutOff) {
432 if (strategyName == "Flat") {
433 m_weightedSum = std::make_unique<FlatWeighting>();
434 } else if (strategyName == "Linear") {
435 m_weightedSum = std::make_unique<LinearWeighting>(cutOff);
436 } else if (strategyName == "Parabolic") {
437 m_weightedSum = std::make_unique<ParabolicWeighting>(cutOff);
438 } else if (strategyName == "Gaussian") {
439 m_weightedSum = std::make_unique<GaussianWeightingnD>(cutOff, getProperty("Sigma"));
440 }
441 g_log.information() << "Smoothing with " << strategyName << " Weighting\n";
442}
443
450double SmoothNeighbours::translateToMeters(const std::string &radiusUnits, const double &enteredRadius) const {
451 double translatedRadius = 0;
452 if (radiusUnits == "Meters") {
453 // Nothing more to do.
454 translatedRadius = enteredRadius;
455 } else if (radiusUnits == "NumberOfPixels") {
456 // Get the first idetector from the workspace index 0.
457 const auto &firstDet = m_inWS->spectrumInfo().detector(0);
458 // Find the bounding box of that detector
459 BoundingBox bbox;
460 firstDet.getBoundingBox(bbox);
461 // Multiply (meters/pixels) by number of pixels, note that enteredRadius
462 // takes on meaning of the number of pixels.
463 translatedRadius = bbox.width().norm() * enteredRadius;
464 } else {
465 const std::string message = "SmoothNeighbours::translateToMeters, Unknown Unit: " + radiusUnits;
466 throw std::invalid_argument(message);
467 }
468 return translatedRadius;
469}
470
475 m_inWS = getProperty("InputWorkspace");
476
477 m_preserveEvents = getProperty("PreserveEvents");
478 m_expandSumAllPixels = getProperty("ExpandSumAllPixels");
479
480 // Use the unit type to translate the entered radius into meters.
481 m_radius = translateToMeters(getProperty("RadiusUnits"), getProperty("Radius"));
482
484
485 m_adjX = getProperty("AdjX");
486 m_adjY = getProperty("AdjY");
487 m_edge = getProperty("ZeroEdgePixels");
488
489 m_nNeighbours = getProperty("NumberOfNeighbours");
490
491 // Progress reporting, first for the sorting
492 m_progress = std::make_unique<Progress>(this, 0.0, 0.2, m_inWS->getNumberHistograms());
493
494 // Clean up when we are done
495 CallOnExit resetInWSOnExit([this]() { m_inWS.reset(); });
496 CallOnExit resetNeighboursOnExit([this]() {
497 m_neighbours.clear();
498 m_neighbours.shrink_to_fit();
499 });
500
501 // Run the appropriate method depending on the type of the instrument
502 if (m_inWS->getInstrument()->containsRectDetectors() == Instrument::ContainsState::Full)
504 else
506
507 auto wsEvent = std::dynamic_pointer_cast<EventWorkspace>(m_inWS);
508 if (wsEvent)
509 wsEvent->sortAll(TOF_SORT, m_progress.get());
510 if (wsEvent && m_preserveEvents)
511 this->execEvent(wsEvent);
512 else
513 this->execWorkspace2D();
514}
515
516//--------------------------------------------------------------------------------------------
519 m_progress->resetNumSteps(m_inWS->getNumberHistograms(), 0.5, 1.0);
520
521 // Get some stuff from the input workspace
522 const size_t numberOfSpectra = m_outWI;
523
524 const size_t YLength = m_inWS->blocksize();
525
527 // Make a brand new Workspace2D
528 if (std::dynamic_pointer_cast<OffsetsWorkspace>(m_inWS)) {
529 g_log.information() << "Creating new OffsetsWorkspace\n";
530 outWS = MatrixWorkspace_sptr(new OffsetsWorkspace(m_inWS->getInstrument()));
531 } else {
532 outWS = std::dynamic_pointer_cast<MatrixWorkspace>(
533 API::WorkspaceFactory::Instance().create("Workspace2D", numberOfSpectra, YLength + 1, YLength));
534 }
535 this->setProperty("OutputWorkspace", outWS);
536
537 setupNewInstrument(*outWS);
538
539 // Go through all the output workspace
541 for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
543
544 auto &outSpec = outWS->getSpectrum(outWIi);
545 // Reset the Y and E vectors
546 outSpec.clearData();
547 auto &outY = outSpec.mutableY();
548 // We will temporarily carry the squared error
549 auto &outE = outSpec.mutableE();
550 // tmp to carry the X Data.
551 auto &outX = outSpec.mutableX();
552
553 // Which are the neighbours?
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;
560
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();
565
566 for (size_t i = 0; i < YLength; i++) {
567 // Add the weighted signal
568 outY[i] += inY[i] * weight;
569 // Square the error, scale by weight (which you have to square too),
570 // then add in quadrature
571 double errorSquared = inE[i];
572 errorSquared *= errorSquared;
573 errorSquared *= weightSquared;
574 outE[i] += errorSquared;
575 // Copy the X values as well
576 outX[i] = inX[i];
577 }
578 if (m_inWS->isHistogramData()) {
579 outX[YLength] = inX[YLength];
580 }
581 } //(each neighbour)
582
583 // Now un-square the error, since we summed it in quadrature
584 for (size_t i = 0; i < YLength; i++)
585 outE[i] = sqrt(outE[i]);
586
587 m_progress->report("Summing");
589 }
592 spreadPixels(outWS);
593}
594
595//--------------------------------------------------------------------------------------------
599 // Copy geometry over.
600 API::WorkspaceFactory::Instance().initializeFromParent(*m_inWS, outWS, false);
601
602 // Go through all the output workspace
603 size_t numberOfSpectra = outWS.getNumberHistograms();
604
605 for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
606 auto &outSpec = outWS.getSpectrum(outWIi);
607
608 // Reset detectors
609 outSpec.clearDetectorIDs();
610
611 // Which are the neighbours?
612 for (const auto &neighbor : m_neighbours[outWIi]) {
613 const auto &inSpec = m_inWS->getSpectrum(neighbor.first);
614 outSpec.addDetectorIDs(inSpec.getDetectorIDs());
615 }
616 }
617}
618//--------------------------------------------------------------------------------------------
622 // Get some stuff from the input workspace
623 const size_t numberOfSpectra = m_inWS->getNumberHistograms();
624
625 const size_t YLength = m_inWS->blocksize();
626
628 // Make a brand new Workspace2D
629 if (std::dynamic_pointer_cast<OffsetsWorkspace>(m_inWS)) {
630 g_log.information() << "Creating new OffsetsWorkspace\n";
631 outws2 = std::make_shared<OffsetsWorkspace>(m_inWS->getInstrument());
632 } else {
633 outws2 = std::dynamic_pointer_cast<MatrixWorkspace>(
634 API::WorkspaceFactory::Instance().create("Workspace2D", numberOfSpectra, YLength + 1, YLength));
635 }
636
637 // Copy geometry over.
638 API::WorkspaceFactory::Instance().initializeFromParent(*m_inWS, *outws2, false);
639 // Go through all the input workspace
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());
645 // Zero the Y and E vectors
646 outSpec2.clearData();
647 }
648
649 // Go through all the output workspace
650 const size_t numberOfSpectra2 = outWS->getNumberHistograms();
651 for (int outWIi = 0; outWIi < int(numberOfSpectra2); outWIi++) {
652
653 // Which are the neighbours?
654 for (const auto &neighbor : m_neighbours[outWIi]) {
655 outws2->setHistogram(neighbor.first, outWS->histogram(outWIi));
656 }
657 }
658 this->setProperty("OutputWorkspace", outws2);
659}
660//--------------------------------------------------------------------------------------------
665 m_progress->resetNumSteps(m_inWS->getNumberHistograms(), 0.5, 1.0);
666
667 // Get some stuff from the input workspace
668 const size_t numberOfSpectra = m_outWI;
669 const auto YLength = static_cast<int>(m_inWS->blocksize());
670
672 // Make a brand new EventWorkspace
673 outWS = std::dynamic_pointer_cast<EventWorkspace>(
674 API::WorkspaceFactory::Instance().create("EventWorkspace", numberOfSpectra, YLength + 1, YLength));
675 // Copy geometry over.
676 API::WorkspaceFactory::Instance().initializeFromParent(*ws, *outWS, false);
677 // Ensure thread-safety
678 outWS->sortAll(TOF_SORT, nullptr);
679
680 this->setProperty("OutputWorkspace", std::dynamic_pointer_cast<MatrixWorkspace>(outWS));
681
682 // Go through all the output workspace
684 for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
686
687 // Create the output event list (empty)
688 EventList &outEL = outWS->getSpectrum(outWIi);
689
690 // Which are the neighbours?
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;
695 // if(sum)outEL.copyInfoFrom(*ws->getSpectrum(inWI));
696 double weight = it->second;
697 // Copy the event list
698 EventList tmpEL = ws->getSpectrum(inWI);
699 // Scale it
700 tmpEL *= weight;
701 // Add it
702 outEL += tmpEL;
703 }
704
705 m_progress->report("Summing");
707 }
709
710 // Give the 0-th X bins to all the output spectra.
711 outWS->setAllX(m_inWS->binEdges(0));
713 spreadPixels(outWS);
714}
715
716} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#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...
Definition: MultiThreaded.h:94
#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
Callable m_callable
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
void clearDetectorIDs()
Clear the detector IDs set.
Definition: ISpectrum.cpp:117
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.
Definition: Progress.h:25
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.
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.
A class for holding :
Definition: EventList.h:56
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.
Definition: BoundingBox.h:34
Kernel::V3D width() const
Returns the width of the box.
Definition: BoundingBox.h:98
Exception for when an item is not found in a collection.
Definition: Exception.h:145
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.
Definition: Logger.cpp:114
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
double norm() const noexcept
Definition: V3D.h:263
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.
Definition: MultiThreaded.h:22
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.
Definition: IDTypes.h:16
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54