SmoothNeighbours.cpp 30.8 KB
Newer Older
1
#include "MantidAlgorithms/SmoothNeighbours.h"
2
#include "MantidAPI/InstrumentValidator.h"
3
#include "MantidDataObjects/EventList.h"
4
#include "MantidDataObjects/EventWorkspace.h"
5
#include "MantidDataObjects/OffsetsWorkspace.h"
6
#include "MantidDataObjects/Workspace2D.h"
7
#include "MantidGeometry/ICompAssembly.h"
8
#include "MantidGeometry/IComponent.h"
9
#include "MantidGeometry/Instrument/RectangularDetector.h"
10
#include "MantidGeometry/Instrument/DetectorGroup.h"
11
#include "MantidKernel/BoundedValidator.h"
12
#include "MantidKernel/EnabledWhenProperty.h"
13
14
#include "MantidKernel/ListValidator.h"

15
#include <boost/algorithm/string.hpp>
16
#include <boost/regex.hpp>
17

18
19
20
21
using namespace Mantid::Kernel;
using namespace Mantid::Geometry;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
22
using std::map;
23

24
typedef std::vector<Mantid::Kernel::Property *> VecProperties;
25
26
typedef const VecProperties ConstVecProperties;

27
28
namespace Mantid {
namespace Algorithms {
29
30
31
32

// Register the class into the algorithm factory
DECLARE_ALGORITHM(SmoothNeighbours)

33
34
// Used in custom GUI. Make sure you change them in SmoothNeighboursDialog.cpp
// as well.
35
36
37
const std::string SmoothNeighbours::NON_UNIFORM_GROUP = "NonUniform Detectors";
const std::string SmoothNeighbours::RECTANGULAR_GROUP = "Rectangular Detectors";
const std::string SmoothNeighbours::INPUT_WORKSPACE = "InputWorkspace";
38

39
SmoothNeighbours::SmoothNeighbours()
40
41
42
43
    : API::Algorithm(), AdjX(0), AdjY(0), Edge(0), Radius(0.), nNeighbours(0),
      WeightedSum(new NullWeighting), PreserveEvents(false),
      expandSumAllPixels(false), outWI(0), inWS(), m_neighbours(),
      m_prog(NULL) {}
44

45
46
47
/** Initialisation method.
 *
 */
48
49
50
51
52
53
54
55
56
57
58
59
void SmoothNeighbours::init() {
  declareProperty(new WorkspaceProperty<MatrixWorkspace>(
                      INPUT_WORKSPACE, "", Direction::Input,
                      boost::make_shared<InstrumentValidator>()),
                  "The workspace containing the spectra to be averaged.");
  declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace", "",
                                                         Direction::Output),
                  "The name of the workspace to be created as the output of "
                  "the algorithm.");

  // Unsigned double
  auto mustBePositiveDouble = boost::make_shared<BoundedValidator<double>>();
60
61
  mustBePositiveDouble->setLower(0.0);

62
63
64
  // Unsigned int.
  auto mustBePositive = boost::make_shared<BoundedValidator<int>>();
  mustBePositive->setLower(0);
65

66
  std::vector<std::string> propOptions;
67
68
69
70
71
72
73
74
75
76
77
78
  propOptions.push_back("Flat");
  propOptions.push_back("Linear");
  propOptions.push_back("Parabolic");
  propOptions.push_back("Gaussian");
  declareProperty("WeightedSum", "Flat",
                  boost::make_shared<StringListValidator>(propOptions),
                  "What sort of Weighting scheme to use?\n"
                  "  Flat: Effectively no-weighting, all weights are 1.\n"
                  "  Linear: Linear weighting 1 - r/R from origin.\n"
                  "  Parabolic : Weighting as cutoff - x + cutoff - y + 1."
                  "  Gaussian : Uses the absolute distance x^2 + y^2 ... "
                  "normalised by the cutoff^2");
79

80
81
82
83
84
  declareProperty(
      "Sigma", 0.5, mustBePositiveDouble,
      "Sigma value for gaussian weighting schemes. Defaults to 0.5. ");
  setPropertySettings(
      "Sigma", new EnabledWhenProperty("WeightedSum", IS_EQUAL_TO, "Gaussian"));
85

86
87
88
  declareProperty(
      "IgnoreMaskedDetectors", true,
      "If true, do not consider masked detectors in the NN search.");
89

90
91
92
93
  declareProperty("PreserveEvents", true, "If the InputWorkspace is an "
                                          "EventWorkspace, this will preserve "
                                          "the full event list (warning: this "
                                          "will use much more memory!).");
94

95
96
  // -- Rectangular properties
  // ----------------------------------------------------------------------
97

98
99
100
101
  declareProperty(
      "AdjX", 1, mustBePositive,
      "The number of X (horizontal) adjacent pixels to average together. "
      "Only for instruments with RectangularDetectors. ");
Owen Arnold's avatar
Owen Arnold committed
102

103
104
105
106
  declareProperty(
      "AdjY", 1, mustBePositive,
      "The number of Y (vertical) adjacent pixels to average together. "
      "Only for instruments with RectangularDetectors. ");
107

108
109
110
111
112
  declareProperty(
      "SumPixelsX", 1, mustBePositive,
      "The total number of X (horizontal) adjacent pixels to sum together. "
      "Only for instruments with RectangularDetectors.  AdjX will be ignored "
      "if SumPixelsX > 1.");
113

114
115
116
117
118
  declareProperty(
      "SumPixelsY", 1, mustBePositive,
      "The total number of Y (vertical) adjacent pixels to sum together. "
      "Only for instruments with RectangularDetectors. AdjY will be ignored if "
      "SumPixelsY > 1");
119
120

  declareProperty("ZeroEdgePixels", 0, mustBePositive,
121
122
                  "The number of pixels to zero at edges. "
                  "Only for instruments with RectangularDetectors. ");
123
124
125
126
127
128
129

  setPropertyGroup("AdjX", RECTANGULAR_GROUP);
  setPropertyGroup("AdjY", RECTANGULAR_GROUP);
  setPropertyGroup("SumPixelsX", RECTANGULAR_GROUP);
  setPropertyGroup("SumPixelsY", RECTANGULAR_GROUP);
  setPropertyGroup("ZeroEdgePixels", RECTANGULAR_GROUP);

130
131
  // -- Non-uniform properties
  // ----------------------------------------------------------------------
132
133
134
135

  std::vector<std::string> radiusPropOptions;
  radiusPropOptions.push_back("Meters");
  radiusPropOptions.push_back("NumberOfPixels");
136
137
138
139
140
141
  declareProperty(
      "RadiusUnits", "Meters",
      boost::make_shared<StringListValidator>(radiusPropOptions),
      "Units used to specify the radius?\n"
      "  Meters : Radius is in meters.\n"
      "  NumberOfPixels : Radius is in terms of the number of pixels.");
142

143
144
  declareProperty(
      "Radius", 0.0, mustBePositiveDouble,
145
      "The radius around a pixel to look for nearest neighbours to average. \n"
146
147
      "If 0, will use the AdjX and AdjY parameters for rectangular detectors "
      "instead.");
148

149
150
151
  declareProperty("NumberOfNeighbours", 8, mustBePositive,
                  "Number of nearest neighbouring pixels.\n"
                  "Alternative to providing the radius. The default is 8.");
152

153
154
155
  declareProperty("SumNumberOfNeighbours", 1,
                  "Sum nearest neighbouring pixels with same parent.\n"
                  "Number of pixels will be reduced. The default is false.");
156

157
158
159
160
  declareProperty("ExpandSumAllPixels", false,
                  "OuputWorkspace will have same number of pixels as "
                  "InputWorkspace using SumPixelsX and SumPixelsY.  Individual "
                  "pixels will have averages.");
161

162
163
164
165
  setPropertyGroup("RadiusUnits", NON_UNIFORM_GROUP);
  setPropertyGroup("Radius", NON_UNIFORM_GROUP);
  setPropertyGroup("NumberOfNeighbours", NON_UNIFORM_GROUP);
  setPropertyGroup("SumNumberOfNeighbours", NON_UNIFORM_GROUP);
166
167
}

168
169
170
//--------------------------------------------------------------------------------------------
/** Fill the neighbours list given the AdjX AdjY parameters and an
 * instrument with rectangular detectors.
171
 */
172
void SmoothNeighbours::findNeighboursRectangular() {
173
174
  g_log.debug("SmoothNeighbours processing assuming rectangular detectors.");

175
  m_prog->resetNumSteps(inWS->getNumberHistograms(), 0.2, 0.5);
176

177
  Instrument_const_sptr inst = inWS->getInstrument();
178

179
180
181
  // To get the workspace index from the detector ID
  const detid2index_map pixel_to_wi =
      inWS->getDetectorIDToWorkspaceIndexMap(true);
182

183
184
  // std::cout << " inst->nelements() " << inst->nelements() << "\n";
  Progress prog(this, 0.0, 1.0, inst->nelements());
185

186
187
188
  // Build a list of Rectangular Detectors
  std::vector<boost::shared_ptr<RectangularDetector>> detList;
  for (int i = 0; i < inst->nelements(); i++) {
189
190
    boost::shared_ptr<RectangularDetector> det;
    boost::shared_ptr<ICompAssembly> assem;
191
    boost::shared_ptr<ICompAssembly> assem2;
192

193
194
    det = boost::dynamic_pointer_cast<RectangularDetector>((*inst)[i]);
    if (det) {
195
      detList.push_back(det);
196
197
198
199
200
201
202
203
204
    } else {
      // Also, look in the first sub-level for RectangularDetectors (e.g. PG3).
      // We are not doing a full recursive search since that will be very long
      // for lots of pixels.
      assem = boost::dynamic_pointer_cast<ICompAssembly>((*inst)[i]);
      if (assem) {
        for (int j = 0; j < assem->nelements(); j++) {
          det = boost::dynamic_pointer_cast<RectangularDetector>((*assem)[j]);
          if (det) {
205
            detList.push_back(det);
206

207
208
209
210
211
212
213
214
215
216
217
          } else {
            // Also, look in the second sub-level for RectangularDetectors (e.g.
            // PG3).
            // We are not doing a full recursive search since that will be very
            // long for lots of pixels.
            assem2 = boost::dynamic_pointer_cast<ICompAssembly>((*assem)[j]);
            if (assem2) {
              for (int k = 0; k < assem2->nelements(); k++) {
                det = boost::dynamic_pointer_cast<RectangularDetector>(
                    (*assem2)[k]);
                if (det) {
218
219
220
221
222
                  detList.push_back(det);
                }
              }
            }
          }
223
224
225
226
227
        }
      }
    }
  }

228
  if (detList.empty()) {
229
    // Not rectangular so use Nearest Neighbours
230
    Radius = translateToMeters("NumberOfPixels", std::max(AdjX, AdjY));
231
232
    setWeightingStrategy("Flat", Radius);
    nNeighbours = AdjX * AdjY - 1;
233
    findNeighboursUbiqutious();
234
    return;
235
  }
236

237
238
  // Resize the vector we are setting
  m_neighbours.resize(inWS->getNumberHistograms());
239
240
241
242
243
244
  int StartX = -AdjX;
  int StartY = -AdjY;
  int EndX = AdjX;
  int EndY = AdjY;
  int SumX = getProperty("SumPixelsX");
  int SumY = getProperty("SumPixelsY");
245
246
  bool sum = (SumX * SumY > 1) ? true : false;
  if (sum) {
247
248
249
250
    StartX = 0;
    StartY = 0;
    EndX = SumX - 1;
    EndY = SumY - 1;
251
252
253
254
  }

  outWI = 0;
  // Build a map to sort by the detectorID
255
256
257
  std::vector<std::pair<int, int>> v1;
  for (int i = 0; i < static_cast<int>(detList.size()); i++)
    v1.push_back(std::pair<int, int>(detList[i]->getAtXY(0, 0)->getID(), i));
258
259

  // To sort in descending order
260
261
  if (sum)
    stable_sort(v1.begin(), v1.end());
262

263
  std::vector<std::pair<int, int>>::iterator Iter1;
264

265
266
  // Loop through the RectangularDetector's we listed before.
  for (Iter1 = v1.begin(); Iter1 != v1.end(); ++Iter1) {
267
    int i = (*Iter1).second;
268
269
    boost::shared_ptr<RectangularDetector> det = detList[i];
    std::string det_name = det->getName();
270
271
272
    if (det) {
      for (int j = 0; j < det->xpixels(); j += SumX) {
        for (int k = 0; k < det->ypixels(); k += SumY) {
273
274
          double totalWeight = 0;
          // Neighbours and weights
275
          std::vector<weightedNeighbour> neighbours;
276

277
278
279
          for (int ix = StartX; ix <= EndX; ix++)
            for (int iy = StartY; iy <= EndY; iy++) {
              // Weights for corners=1; higher for center and adjacent pixels
280
              double smweight = WeightedSum->weightAt(AdjX, ix, AdjY, iy);
281

282
283
284
285
286
287
288
              // Find the pixel ID at that XY position on the rectangular
              // detector
              if (j + ix >= det->xpixels() - Edge || j + ix < Edge)
                continue;
              if (k + iy >= det->ypixels() - Edge || k + iy < Edge)
                continue;
              int pixelID = det->getAtXY(j + ix, k + iy)->getID();
289

290
              // Find the corresponding workspace index, if any
291
              auto mapEntry = pixel_to_wi.find(pixelID);
292
              if (mapEntry != pixel_to_wi.end()) {
293
                size_t wi = mapEntry->second;
294
                neighbours.push_back(weightedNeighbour(wi, smweight));
295
296
                // Count the total weight
                totalWeight += smweight;
297
298
              }
            }
299
300

          // Adjust the weights of each neighbour to normalize to unity
Lynch, Vickie's avatar
Lynch, Vickie committed
301
          if (!sum || expandSumAllPixels)
302
303
            for (size_t q = 0; q < neighbours.size(); q++)
              neighbours[q].second /= totalWeight;
304
305
306

          // Save the list of neighbours for this output workspace index.
          m_neighbours[outWI] = neighbours;
307
          outWI++;
308
309

          m_prog->report("Finding Neighbours");
310
311
        }
      }
312
    }
313
    prog.report(det_name);
314
315
316
317
318
319
  }
}

//--------------------------------------------------------------------------------------------
/** Use NearestNeighbours to find the neighbours for any instrument
 */
320
321
322
void SmoothNeighbours::findNeighboursUbiqutious() {
  g_log.debug(
      "SmoothNeighbours processing NOT assuming rectangular detectors.");
323
  /*
324
325
    This will cause the Workspace to rebuild the nearest neighbours map, so that
    we can pick-up any of the properties specified
326
327
328
329
    for this algorithm in the constructor for the NearestNeighboursObject.
  */
  inWS->rebuildNearestNeighbours();

330
331
332
333
  m_prog->resetNumSteps(inWS->getNumberHistograms(), 0.2, 0.5);
  this->progress(0.2, "Building Neighbour Map");

  Instrument_const_sptr inst = inWS->getInstrument();
334
  const spec2index_map spec2index = inWS->getSpectrumToWorkspaceIndexMap();
335
336
337

  // Resize the vector we are setting
  m_neighbours.resize(inWS->getNumberHistograms());
338

339
340
  bool ignoreMaskedDetectors = getProperty("IgnoreMaskedDetectors");

341
  // Cull by radius
342
343
  RadiusFilter radiusFilter(Radius);

Owen Arnold's avatar
Owen Arnold committed
344
  IDetector_const_sptr det;
345
  // Go through every input workspace pixel
346
  outWI = 0;
347
  int sum = getProperty("SumNumberOfNeighbours");
348
349
  boost::shared_ptr<const Geometry::IComponent> parent, neighbParent,
      grandparent, neighbGParent;
350
  auto used = new bool[inWS->getNumberHistograms()];
351
352
  if (sum > 1) {
    for (size_t wi = 0; wi < inWS->getNumberHistograms(); wi++)
353
354
      used[wi] = false;
  }
355
356
357
358
  for (size_t wi = 0; wi < inWS->getNumberHistograms(); wi++) {
    if (sum > 1)
      if (used[wi])
        continue;
Owen Arnold's avatar
Owen Arnold committed
359
    // We want to skip monitors
360
    try {
361
      // Get the list of detectors in this pixel
362
      const std::set<detid_t> &dets = inWS->getSpectrum(wi)->getDetectorIDs();
363
      det = inst->getDetector(*dets.begin());
364
365
366
367
368
369
370
371
      if (det->isMonitor())
        continue; // skip monitor
      if (det->isMasked()) {
        // Calibration masks many detectors, but there should be 0s after
        // smoothing
        if (sum == 1)
          outWI++;
        continue; // skip masked detectors
372
      }
373
      if (sum > 1) {
374
        parent = det->getParent();
375
376
        if (parent)
          grandparent = parent->getParent();
377
      }
378
379
    } catch (Kernel::Exception::NotFoundError &) {
      continue; // skip missing detector
Owen Arnold's avatar
Owen Arnold committed
380
381
    }

382
    specid_t inSpec = inWS->getSpectrum(wi)->getSpectrumNo();
383

384
385
386
    // Step one - Get the number of specified neighbours
    SpectraDistanceMap insideGrid =
        inWS->getNeighboursExact(inSpec, nNeighbours, ignoreMaskedDetectors);
387

388
    // Step two - Filter the results by the radius cut off.
389
    SpectraDistanceMap neighbSpectra = radiusFilter.apply(insideGrid);
390

391
    // Force the central pixel to always be there
392
393
    // There seems to be a bug in nearestNeighbours, returns distance != 0.0 for
    // the central pixel. So we force distance = 0
394
    neighbSpectra[inSpec] = V3D(0.0, 0.0, 0.0);
395
396
397

    // Neighbours and weights list
    double totalWeight = 0;
398
    int noNeigh = 0;
399
    std::vector<weightedNeighbour> neighbours;
400
401

    // Convert from spectrum numbers to workspace indices
402
    for (auto it = neighbSpectra.begin();
403
         it != neighbSpectra.end(); ++it) {
404
405
      specid_t spec = it->first;

406
407
      // Use the weighting strategy to calculate the weight.
      double weight = WeightedSum->weightAt(it->second);
408

409
      if (weight > 0) {
410
        // Find the corresponding workspace index
411
        auto mapIt = spec2index.find(spec);
412
        if (mapIt != spec2index.end()) {
413
          size_t neighWI = mapIt->second;
414
415
416
417
          if (sum > 1) {
            // Get the list of detectors in this pixel
            const std::set<detid_t> &dets =
                inWS->getSpectrum(neighWI)->getDetectorIDs();
418
419
            det = inst->getDetector(*dets.begin());
            neighbParent = det->getParent();
420
            neighbGParent = neighbParent->getParent();
421
422
423
424
425
            if (noNeigh >= sum ||
                neighbParent->getName().compare(parent->getName()) != 0 ||
                neighbGParent->getName().compare(grandparent->getName()) != 0 ||
                used[neighWI])
              continue;
426
427
428
            noNeigh++;
            used[neighWI] = true;
          }
429
          neighbours.push_back(weightedNeighbour(neighWI, weight));
430
431
432
          totalWeight += weight;
        }
      }
433
    }
434
435

    // Adjust the weights of each neighbour to normalize to unity
436
437
438
    if (sum == 1)
      for (size_t q = 0; q < neighbours.size(); q++)
        neighbours[q].second /= totalWeight;
439
440

    // Save the list of neighbours for this output workspace index.
441
    m_neighbours[outWI] = neighbours;
442
    outWI++;
443
444
445
446

    m_prog->report("Finding Neighbours");
  } // each workspace index

447
  delete[] used;
448
449
}

450
/**
451
452
453
454
Attempts to reset the Weight based on the strategyName provided. Note that if
these conditional
statements fail to override the existing WeightedSum member, it should stay as a
NullWeighting, which
455
will throw during usage.
456
457
@param strategyName : The name of the weighting strategy to use
@param cutOff : The cutoff distance
458
*/
459
460
461
void SmoothNeighbours::setWeightingStrategy(const std::string strategyName,
                                            double &cutOff) {
  if (strategyName == "Flat") {
462
463
    boost::scoped_ptr<WeightingStrategy> flatStrategy(new FlatWeighting);
    WeightedSum.swap(flatStrategy);
464
    g_log.information("Smoothing with Flat Weighting");
465
466
467
  } else if (strategyName == "Linear") {
    boost::scoped_ptr<WeightingStrategy> linearStrategy(
        new LinearWeighting(cutOff));
468
    WeightedSum.swap(linearStrategy);
469
    g_log.information("Smoothing with Linear Weighting");
470
471
472
  } else if (strategyName == "Parabolic") {
    boost::scoped_ptr<WeightingStrategy> parabolicStrategy(
        new ParabolicWeighting(cutOff));
473
    WeightedSum.swap(parabolicStrategy);
474
    g_log.information("Smoothing with Parabolic Weighting");
475
  } else if (strategyName == "Gaussian") {
476
    double sigma = getProperty("Sigma");
477
478
    boost::scoped_ptr<WeightingStrategy> gaussian1DStrategy(
        new GaussianWeightingnD(cutOff, sigma));
479
    WeightedSum.swap(gaussian1DStrategy);
480
    g_log.information("Smoothing with Gaussian Weighting");
481
  }
482
483
}

484
485
486
/*
Fetch the instrument associated with the workspace
@return const shared pointer to the instrument,
487
*/
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
Instrument_const_sptr SmoothNeighbours::fetchInstrument() const {
  Instrument_const_sptr instrument;
  EventWorkspace_sptr wsEvent =
      boost::dynamic_pointer_cast<EventWorkspace>(inWS);
  MatrixWorkspace_sptr wsMatrix =
      boost::dynamic_pointer_cast<MatrixWorkspace>(inWS);
  if (wsEvent) {
    instrument =
        boost::dynamic_pointer_cast<EventWorkspace>(inWS)->getInstrument();
  } else if (wsMatrix) {
    instrument =
        boost::dynamic_pointer_cast<MatrixWorkspace>(inWS)->getInstrument();
  } else {
    throw std::invalid_argument("Neither a Matrix Workspace or an "
                                "EventWorkpace provided to SmoothNeighbours.");
  }
  return instrument;
505
506
507
508
509
}

/**
Translate the radius into meters.
@param radiusUnits : The name of the radius units
510
511
@param enteredRadius : The numerical value of the radius in whatever units have
been specified
512
*/
513
514
double SmoothNeighbours::translateToMeters(const std::string radiusUnits,
                                           const double &enteredRadius) {
515
  double translatedRadius = 0;
516
  if (radiusUnits == "Meters") {
517
    // Nothing more to do.
518
519
    translatedRadius = enteredRadius;
  } else if (radiusUnits == "NumberOfPixels") {
520
521
    // Fetch the instrument.
    Instrument_const_sptr instrument = fetchInstrument();
522

523
524
    // Get the first idetector from the workspace index 0.
    IDetector_const_sptr firstDet = inWS->getDetector(0);
525
526
527
    // Find the bounding box of that detector
    BoundingBox bbox;
    firstDet->getBoundingBox(bbox);
528
529
    // Multiply (meters/pixels) by number of pixels, note that enteredRadius
    // takes on meaning of the number of pixels.
530
    translatedRadius = bbox.width().norm() * enteredRadius;
531
532
533
  } else {
    std::string message =
        "SmoothNeighbours::translateToMeters, Unknown Unit: " + radiusUnits;
534
535
536
537
538
    throw std::invalid_argument(message);
  }
  return translatedRadius;
}

539
540
541
542
/*
Check whether the properties provided are all in their default state.
@param properties : Vector of mantid property pointers
@return True only if they are all default, otherwise False.
543
*/
544
bool areAllDefault(ConstVecProperties &properties) {
545
  bool areAllDefault = false;
546
  for (auto it = properties.cbegin(); it != properties.cend(); ++it) {
547
    if (!(*it)->isDefault()) {
548
      return areAllDefault;
549
550
    }
  }
551
552
  areAllDefault = true;
  return areAllDefault;
553
554
}

555
556
557
558
//--------------------------------------------------------------------------------------------
/** Executes the algorithm
 *
 */
559
void SmoothNeighbours::exec() {
560
  inWS = getProperty("InputWorkspace");
561
562

  PreserveEvents = getProperty("PreserveEvents");
563

564
565
  expandSumAllPixels = getProperty("ExpandSumAllPixels");

566
  // Use the unit type to translate the entered radius into meters.
567
  Radius = translateToMeters(getProperty("RadiusUnits"), getProperty("Radius"));
568

569
  setWeightingStrategy(getProperty("WeightedSum"), Radius);
570

Arturs Bekasovs's avatar
Arturs Bekasovs committed
571
572
573
574
575
576
  AdjX = getProperty("AdjX");
  AdjY = getProperty("AdjY");
  Edge = getProperty("ZeroEdgePixels");

  nNeighbours = getProperty("NumberOfNeighbours");

577
578
579
  // Progress reporting, first for the sorting
  m_prog = new Progress(this, 0.0, 0.2, inWS->getNumberHistograms());

Arturs Bekasovs's avatar
Arturs Bekasovs committed
580
  // Run the appropriate method depending on the type of the instrument
581
582
  if (inWS->getInstrument()->containsRectDetectors() ==
      Instrument::ContainsState::Full)
583
584
    findNeighboursRectangular();
  else
585
    findNeighboursUbiqutious();
586

587
588
  EventWorkspace_sptr wsEvent =
      boost::dynamic_pointer_cast<EventWorkspace>(inWS);
589
590
591
592
  if (wsEvent)
    wsEvent->sortAll(TOF_SORT, m_prog);

  if (!wsEvent || !PreserveEvents)
593
    this->execWorkspace2D();
594
595
596
  else if (wsEvent)
    this->execEvent(wsEvent);
  else
597
598
    throw std::runtime_error("This algorithm requires a Workspace2D or "
                             "EventWorkspace as its input.");
599
600
601
602
}

//--------------------------------------------------------------------------------------------
/** Execute the algorithm for a Workspace2D/don't preserve events input */
603
void SmoothNeighbours::execWorkspace2D() {
604
605
  m_prog->resetNumSteps(inWS->getNumberHistograms(), 0.5, 1.0);

606
  // Get some stuff from the input workspace
607
  const size_t numberOfSpectra = outWI;
608

609
610
611
  const size_t YLength = inWS->blocksize();

  MatrixWorkspace_sptr outWS;
612
613
  // Make a brand new Workspace2D
  if (boost::dynamic_pointer_cast<OffsetsWorkspace>(inWS)) {
614
    g_log.information() << "Creating new OffsetsWorkspace\n";
615
    outWS = MatrixWorkspace_sptr(new OffsetsWorkspace(inWS->getInstrument()));
616
617
618
619
  } else {
    outWS = boost::dynamic_pointer_cast<MatrixWorkspace>(
        API::WorkspaceFactory::Instance().create("Workspace2D", numberOfSpectra,
                                                 YLength + 1, YLength));
620
  }
621
  this->setProperty("OutputWorkspace", outWS);
622
623
624

  setupNewInstrument(outWS);

625
  // Copy geometry over.
626
  // API::WorkspaceFactory::Instance().initializeFromParent(inWS, outWS, false);
627
628

  // Go through all the output workspace
629
  PARALLEL_FOR2(inWS, outWS)
630
  for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
631
632
    PARALLEL_START_INTERUPT_REGION

633
    ISpectrum *outSpec = outWS->getSpectrum(outWIi);
634
635
    // Reset the Y and E vectors
    outSpec->clearData();
636
    MantidVec &outY = outSpec->dataY();
637
    // We will temporarily carry the squared error
638
    MantidVec &outE = outSpec->dataE();
Owen Arnold's avatar
Owen Arnold committed
639
    // tmp to carry the X Data.
640
    MantidVec &outX = outSpec->dataX();
641
642

    // Which are the neighbours?
643
644
645
    std::vector<weightedNeighbour> &neighbours = m_neighbours[outWIi];
    std::vector<weightedNeighbour>::iterator it;
    for (it = neighbours.begin(); it != neighbours.end(); ++it) {
646
647
      size_t inWI = it->first;
      double weight = it->second;
648
      double weightSquared = weight * weight;
649

650
      const ISpectrum *inSpec = inWS->getSpectrum(inWI);
651
      inSpec->lockData();
652
653
654
      const MantidVec &inY = inSpec->readY();
      const MantidVec &inE = inSpec->readE();
      const MantidVec &inX = inSpec->readX();
655

656
      for (size_t i = 0; i < YLength; i++) {
657
658
        // Add the weighted signal
        outY[i] += inY[i] * weight;
659
660
        // Square the error, scale by weight (which you have to square too),
        // then add in quadrature
661
662
        double errorSquared = inE[i];
        errorSquared *= errorSquared;
663
        errorSquared *= weightSquared;
664
        outE[i] += errorSquared;
665
        // Copy the X values as well
Owen Arnold's avatar
Owen Arnold committed
666
667
        outX[i] = inX[i];
      }
668
      if (inWS->isHistogramData()) {
Owen Arnold's avatar
Owen Arnold committed
669
        outX[YLength] = inX[YLength];
670
      }
671
672

      inSpec->unlockData();
673
674
    } //(each neighbour)

675
    // Now un-square the error, since we summed it in quadrature
676
    for (size_t i = 0; i < YLength; i++)
677
678
      outE[i] = sqrt(outE[i]);

679
680
    // Copy the single detector ID (of the center) and spectrum number from the
    // input workspace
681
    // outSpec->copyInfoFrom(*inWS->getSpectrum(outWIi));
682
683
684
685
686

    m_prog->report("Summing");
    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION
687
688
  if (expandSumAllPixels)
    spreadPixels(outWS);
689
690
}

691
692
693
//--------------------------------------------------------------------------------------------
/** Build the instrument/detector setup in workspace
  */
694
695
void SmoothNeighbours::setupNewInstrument(MatrixWorkspace_sptr outws) {
  // Copy geometry over.
696
  API::WorkspaceFactory::Instance().initializeFromParent(inWS, outws, false);
697

698
699
  // Go through all the output workspace
  size_t numberOfSpectra = outws->getNumberHistograms();
700

701
702
  for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
    ISpectrum *outSpec = outws->getSpectrum(outWIi);
703
704
705
706
707
708
709
    /*
    g_log.notice() << "[DBx555] Original spectrum number for wsindex " << outWIi
                   << " = " << outSpec->getSpectrumNo() << std::endl;
    outSpec->setSpectrumNo(outWIi+1);
    */

    // Reset detectors
710
711
    outSpec->clearDetectorIDs();
    ;
712
713

    // Which are the neighbours?
714
715
716
    std::vector<weightedNeighbour> &neighbours = m_neighbours[outWIi];
    std::vector<weightedNeighbour>::iterator it;
    for (it = neighbours.begin(); it != neighbours.end(); ++it) {
717
718
      size_t inWI = it->first;

719
      const ISpectrum *inSpec = inWS->getSpectrum(inWI);
720
721
722
723
724
725
726
727
728

      std::set<detid_t> thesedetids = inSpec->getDetectorIDs();
      outSpec->addDetectorIDs(thesedetids);

    } //(each neighbour)
  }

  return;
}
729
730
731
//--------------------------------------------------------------------------------------------
/** Spread the average over all the pixels
  */
732
733
void SmoothNeighbours::spreadPixels(MatrixWorkspace_sptr outws) {
  // Get some stuff from the input workspace
734
735
736
737
738
  const size_t numberOfSpectra = inWS->getNumberHistograms();

  const size_t YLength = inWS->blocksize();

  MatrixWorkspace_sptr outws2;
739
740
  // Make a brand new Workspace2D
  if (boost::dynamic_pointer_cast<OffsetsWorkspace>(inWS)) {
741
742
    g_log.information() << "Creating new OffsetsWorkspace\n";
    outws2 = MatrixWorkspace_sptr(new OffsetsWorkspace(inWS->getInstrument()));
743
744
745
746
  } else {
    outws2 = boost::dynamic_pointer_cast<MatrixWorkspace>(
        API::WorkspaceFactory::Instance().create("Workspace2D", numberOfSpectra,
                                                 YLength + 1, YLength));
747
748
  }

749
  // Copy geometry over.
750
  API::WorkspaceFactory::Instance().initializeFromParent(inWS, outws2, false);
751
  // Go through all the input workspace
752
753
754
  for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
    ISpectrum *inSpec = inWS->getSpectrum(outWIi);
    MantidVec &inX = inSpec->dataX();
755
756

    std::set<detid_t> thesedetids = inSpec->getDetectorIDs();
757
758
    ISpectrum *outSpec2 = outws2->getSpectrum(outWIi);
    MantidVec &outX = outSpec2->dataX();
759
760
761
762
    outX = inX;
    outSpec2->addDetectorIDs(thesedetids);
    // Zero the Y and E vectors
    outSpec2->clearData();
763
764
    outSpec2->dataY().assign(YLength, 0.0);
    outSpec2->dataE().assign(YLength, 0.0);
765
766
  }

767
768
  // Go through all the output workspace
  const size_t numberOfSpectra2 = outws->getNumberHistograms();
769
770
  for (int outWIi = 0; outWIi < int(numberOfSpectra2); outWIi++) {
    const ISpectrum *outSpec = outws->getSpectrum(outWIi);
771
772

    // Which are the neighbours?
773
774
775
    std::vector<weightedNeighbour> &neighbours = m_neighbours[outWIi];
    std::vector<weightedNeighbour>::iterator it;
    for (it = neighbours.begin(); it != neighbours.end(); ++it) {
776
777
      size_t inWI = it->first;

778
      ISpectrum *outSpec2 = outws2->getSpectrum(inWI);
779
780
      // Reset the Y and E vectors
      outSpec2->clearData();
781
782
783
784
785
786
      MantidVec &out2Y = outSpec2->dataY();
      MantidVec &out2E = outSpec2->dataE();
      MantidVec &out2X = outSpec2->dataX();
      const MantidVec &outY = outSpec->dataY();
      const MantidVec &outE = outSpec->dataE();
      const MantidVec &outX = outSpec->dataX();
787
788
789
790
791
792
793
794
      out2Y = outY;
      out2E = outE;
      out2X = outX;
    } //(each neighbour)
  }
  this->setProperty("OutputWorkspace", outws2);
  return;
}
795
796
797
798
//--------------------------------------------------------------------------------------------
/** Execute the algorithm for a EventWorkspace input
 * @param ws :: EventWorkspace
 */
799
void SmoothNeighbours::execEvent(Mantid::DataObjects::EventWorkspace_sptr ws) {
800
801
  m_prog->resetNumSteps(inWS->getNumberHistograms(), 0.5, 1.0);

802
  // Get some stuff from the input workspace
803
  const size_t numberOfSpectra = outWI;
804
805
806
  const int YLength = static_cast<int>(inWS->blocksize());

  EventWorkspace_sptr outWS;
807
808
809
810
811
  // Make a brand new EventWorkspace
  outWS = boost::dynamic_pointer_cast<EventWorkspace>(
      API::WorkspaceFactory::Instance().create(
          "EventWorkspace", numberOfSpectra, YLength + 1, YLength));
  // Copy geometry over.
812
813
814
815
  API::WorkspaceFactory::Instance().initializeFromParent(ws, outWS, false);
  // Ensure thread-safety
  outWS->sortAll(TOF_SORT, NULL);

816
817
  this->setProperty("OutputWorkspace",
                    boost::dynamic_pointer_cast<MatrixWorkspace>(outWS));
818
819
820

  // Go through all the output workspace
  PARALLEL_FOR2(ws, outWS)
821
  for (int outWIi = 0; outWIi < int(numberOfSpectra); outWIi++) {
822
823
824
    PARALLEL_START_INTERUPT_REGION

    // Create the output event list (empty)
825
    EventList &outEL = outWS->getOrAddEventList(outWIi);
826
827

    // Which are the neighbours?
828
829
830
    std::vector<weightedNeighbour> &neighbours = m_neighbours[outWIi];
    std::vector<weightedNeighbour>::iterator it;
    for (it = neighbours.begin(); it != neighbours.end(); ++it) {
831
      size_t inWI = it->first;
832
      // if(sum)outEL.copyInfoFrom(*ws->getSpectrum(inWI));
833
834
835
836
837
838
839
840
841
      double weight = it->second;
      // Copy the event list
      EventList tmpEL = ws->getEventList(inWI);
      // Scale it
      tmpEL *= weight;
      // Add it
      outEL += tmpEL;
    }

842
843
844
    // Copy the single detector ID (of the center) and spectrum number from the
    // input workspace
    // if (!sum) outEL.copyInfoFrom(*ws->getSpectrum(outWIi));
845
846

    m_prog->report("Summing");
847
    PARALLEL_END_INTERUPT_REGION
848
  }
849
  PARALLEL_CHECK_INTERUPT_REGION
850

851
  // Give the 0-th X bins to all the output spectra.
852
  Kernel::cow_ptr<MantidVec> outX = inWS->refX(0);
853
854
855
  outWS->setAllX(outX);
  if (expandSumAllPixels)
    spreadPixels(outWS);
856
}
857
858
859

} // namespace Algorithms
} // namespace Mantid