FindPeaksMD.cpp 29.1 KB
Newer Older
Alex Buts's avatar
Alex Buts committed
1
#include "MantidMDAlgorithms/FindPeaksMD.h"
2
#include "MantidAPI/Run.h"
3
4
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidDataObjects/MDHistoWorkspace.h"
5
#include "MantidDataObjects/PeaksWorkspace.h"
Lynch, Vickie's avatar
Lynch, Vickie committed
6
#include "MantidGeometry/Crystal/EdgePixel.h"
7
#include "MantidGeometry/Objects/InstrumentRayTracer.h"
8
#include "MantidKernel/BoundedValidator.h"
9
10
#include "MantidKernel/EnabledWhenProperty.h"
#include "MantidKernel/ListValidator.h"
11
#include "MantidKernel/VMD.h"
12

13
#include <map>
14
#include <vector>
15
16
17
18

using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::DataObjects;
19
using namespace Mantid::DataObjects;
20

21
22
23
24
25
26
27
28
29
30
31
32
33
34
namespace Mantid {
namespace MDAlgorithms {
namespace {
// ---------- Template deduction of the event type
// --------------------------------
// See boost::type_traits documentation

/// Type trait to indicate that a general type is not a full MDEvent
template <typename MDE, size_t nd> struct IsFullEvent : boost::false_type {};

/// Specialization of type trait to indicate that a MDEvent is a full event
template <size_t nd> struct IsFullEvent<MDEvent<nd>, nd> : boost::true_type {};

/**
35
 * Specialization if isFullEvent for DataObjects
36
37
38
39
40
41
42
43
 * to return true
 */
template <typename MDE, size_t nd>
bool isFullMDEvent(const boost::true_type &) {
  return true;
}

/**
44
 * Specialization if isFullEvent for DataObjects
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
 * to return false
 */
template <typename MDE, size_t nd>
bool isFullMDEvent(const boost::false_type &) {
  return false;
}

/**
 * Returns true if the templated type is a full MDEvent
 */
template <typename MDE, size_t nd> bool isFullMDEvent() {
  return isFullMDEvent<MDE, nd>(IsFullEvent<MDE, nd>());
}

/**
 * Add the detectors from the given box as contributing detectors to the peak
 * @param peak :: The peak that relates to the box
 * @param box :: A reference to the box containing the peak
 */
template <typename MDE, size_t nd>
void addDetectors(DataObjects::Peak &peak, MDBoxBase<MDE, nd> &box,
                  const boost::true_type &) {
  if (box.getNumChildren() > 0) {
    std::cerr << "Box has children\n";
    addDetectors(peak, box, boost::true_type());
70
  }
71
72
73
74
  MDBox<MDE, nd> *mdBox = dynamic_cast<MDBox<MDE, nd> *>(&box);
  if (!mdBox) {
    throw std::invalid_argument("FindPeaksMD::addDetectors - Unexpected Box "
                                "type, cannot retrieve events");
75
  }
76
77
78
79
  const auto &events = mdBox->getConstEvents();
  auto itend = events.end();
  for (auto it = events.begin(); it != itend; ++it) {
    peak.addContributingDetID(it->getDetectorID());
80
  }
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
}

/// Add detectors based on lean events. Always throws as they do not know their
/// IDs
template <typename MDE, size_t nd>
void addDetectors(DataObjects::Peak &, MDBoxBase<MDE, nd> &,
                  const boost::false_type &) {
  throw std::runtime_error("FindPeaksMD - Workspace contains lean events, "
                           "cannot include detector information");
}

/**
 * Add the detectors from the given box as contributing detectors to the peak
 * @param peak :: The peak that relates to the box
 * @param box :: A reference to the box containing the peak
 */
template <typename MDE, size_t nd>
void addDetectors(DataObjects::Peak &peak, MDBoxBase<MDE, nd> &box) {
  // Compile time deduction of the correct function call
  addDetectors(peak, box, IsFullEvent<MDE, nd>());
}
102
} // namespace
103
104
105
106

// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(FindPeaksMD)

107
const std::string FindPeaksMD::volumeNormalization = "VolumeNormalization";
108
109
const std::string FindPeaksMD::numberOfEventsNormalization =
    "NumberOfEventsNormalization";
110

111
112
113
114
//----------------------------------------------------------------------------------------------
/** Constructor
 */
FindPeaksMD::FindPeaksMD()
115
    : peakWS(), peakRadiusSquared(), DensityThresholdFactor(0.0), m_maxPeaks(0),
116
      m_addDetectors(true), m_densityScaleFactor(1e-6), prog(nullptr), inst(),
117
      m_runNumber(-1), m_peakNumber(1), dimType(), m_goniometer() {}
118
119
120
121
122

//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
 */
void FindPeaksMD::init() {
123
  declareProperty(make_unique<WorkspaceProperty<IMDWorkspace> >(
124
                      "InputWorkspace", "", Direction::Input),
125
126
127
128
                  "An input MDEventWorkspace or MDHistoWorkspace with at least "
                  "3 dimensions.");

  declareProperty(
129
130
      make_unique<PropertyWithValue<double> >("PeakDistanceThreshold", 0.1,
                                              Direction::Input),
131
132
133
134
      "Threshold distance for rejecting peaks that are found to be too close "
      "from each other.\n"
      "This should be some multiple of the radius of a peak. Default: 0.1.");

135
136
  declareProperty(make_unique<PropertyWithValue<int64_t> >("MaxPeaks", 500,
                                                           Direction::Input),
137
                  "Maximum number of peaks to find. Default: 500.");
138

139
140
  std::vector<std::string> strategy = { volumeNormalization,
                                        numberOfEventsNormalization };
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
  declareProperty(
      "PeakFindingStrategy", volumeNormalization,
      boost::make_shared<StringListValidator>(strategy),
      "Strategy for finding peaks in an MD workspace."
      "1. VolumeNormalization: This is the default strategy. It will sort "
      "all boxes in the workspace by deacresing order of signal density "
      "(total weighted event sum divided by box volume).\n"
      "2.NumberOfEventsNormalization: This option is only valid for "
      "MDEventWorkspaces. It will use the total weighted event sum divided"
      "by the number of events. This can improve peak finding for "
      "histogram-based"
      "raw data which has been converted to an EventWorkspace. The threshold "
      "for"
      "peak finding can be controlled by the SingalThresholdFactor property "
      "which should"
      "be larger than 1. Note that this approach does not work for event-based "
      "raw data.\n");
158

159
  declareProperty(make_unique<PropertyWithValue<double> >(
160
                      "DensityThresholdFactor", 10.0, Direction::Input),
161
162
163
164
165
166
                  "The overall signal density of the workspace will be "
                  "multiplied by this factor \n"
                  "to get a threshold signal density below which boxes are NOT "
                  "considered to be peaks. See the help.\n"
                  "Default: 10.0");

167
168
  setPropertySettings("DensityThresholdFactor",
                      make_unique<EnabledWhenProperty>(
169
170
171
                          "PeakFindingStrategy",
                          Mantid::Kernel::ePropertyCriterion::IS_EQUAL_TO,
                          volumeNormalization));
Anton Piccardo-Selg's avatar
Anton Piccardo-Selg committed
172

173
  declareProperty(make_unique<PropertyWithValue<double> >(
174
                      "SignalThresholdFactor", 1.5, Direction::Input),
175
176
177
178
179
180
181
182
                  "The overal signal value (not density!) normalized by the "
                  "number of events is compared to the specified signal "
                  "threshold. Boxes which are below this threshold are NOT "
                  "considered to be peaks."
                  "This property is enabled when the PeakFindingStrategy has "
                  "been set to NumberOfEventsNormalization.\n"
                  "The value of boxes which contain peaks will be above 1. See "
                  "the below for more information.\n"
183
184
                  "Default: 1.50");

185
186
  setPropertySettings("SignalThresholdFactor",
                      make_unique<EnabledWhenProperty>(
187
188
189
                          "PeakFindingStrategy",
                          Mantid::Kernel::ePropertyCriterion::IS_EQUAL_TO,
                          numberOfEventsNormalization));
190

191
192
193
194
195
  declareProperty("CalculateGoniometerForCW", false,
                  "This will calculate the goniometer rotation (around y-axis "
                  "only) for a constant wavelength. This only works for Q "
                  "sample workspaces.");

196
  auto nonNegativeDbl = boost::make_shared<BoundedValidator<double> >();
197
198
199
200
201
202
203
204
205
  nonNegativeDbl->setLower(0);
  declareProperty("Wavelength", DBL_MAX, nonNegativeDbl,
                  "Wavelength to use when calculating goniometer angle");

  setPropertySettings("Wavelength",
                      make_unique<EnabledWhenProperty>(
                          "CalculateGoniometerForCW",
                          Mantid::Kernel::ePropertyCriterion::IS_NOT_DEFAULT));

206
  declareProperty(make_unique<WorkspaceProperty<PeaksWorkspace> >(
207
                      "OutputWorkspace", "", Direction::Output),
208
209
210
211
212
213
                  "An output PeaksWorkspace with the peaks' found positions.");

  declareProperty("AppendPeaks", false,
                  "If checked, then append the peaks in the output workspace "
                  "if it exists. \n"
                  "If unchecked, the output workspace is replaced (Default).");
Lynch, Vickie's avatar
Lynch, Vickie committed
214

215
  auto nonNegativeInt = boost::make_shared<BoundedValidator<int> >();
Lynch, Vickie's avatar
Lynch, Vickie committed
216
  nonNegativeInt->setLower(0);
217
  declareProperty("EdgePixels", 0, nonNegativeInt,
Lynch, Vickie's avatar
Lynch, Vickie committed
218
                  "Remove peaks that are at pixels this close to edge. ");
219
220
221
222
223
224
225
226
227
}

//----------------------------------------------------------------------------------------------
/** Extract needed data from the workspace's experiment info */
void FindPeaksMD::readExperimentInfo(const ExperimentInfo_sptr &ei,
                                     const IMDWorkspace_sptr &ws) {
  // Instrument associated with workspace
  inst = ei->getInstrument();
  // Find the run number
228
  m_runNumber = ei->getRunNumber();
229
230
231
232
233
234
235
236
237
238
239
240

  // Check that the workspace dimensions are in Q-sample-frame or Q-lab-frame.
  std::string dim0 = ws->getDimension(0)->getName();
  if (dim0 == "H") {
    dimType = HKL;
    throw std::runtime_error(
        "Cannot find peaks in a workspace that is already in HKL space.");
  } else if (dim0 == "Q_lab_x") {
    dimType = QLAB;
  } else if (dim0 == "Q_sample_x")
    dimType = QSAMPLE;
  else
241
242
    throw std::runtime_error(
        "Unexpected dimensions: need either Q_lab_x or Q_sample_x.");
243
244

  // Find the goniometer rotation matrix
245
  m_goniometer =
246
247
      Mantid::Kernel::Matrix<double>(3, 3, true); // Default IDENTITY matrix
  try {
248
    m_goniometer = ei->mutableRun().getGoniometerMatrix();
249
250
  }
  catch (std::exception &e) {
251
    g_log.warning() << "Error finding goniometer matrix. It will not be set in "
252
253
                       "the peaks found.\n";
    g_log.warning() << e.what() << '\n';
254
  }
255
256
257
258
259
260
261
}

//----------------------------------------------------------------------------------------------
/** Create and add a Peak to the output workspace
 *
 * @param Q :: Q_lab or Q_sample, depending on workspace
 * @param binCount :: bin count to give to the peak.
Owen Arnold's avatar
Owen Arnold committed
262
 * @param tracer :: Ray tracer to use for detector finding
263
 */
264
265
void FindPeaksMD::addPeak(const V3D &Q, const double binCount,
                          const Geometry::InstrumentRayTracer &tracer) {
266
  try {
267
    auto p = this->createPeak(Q, binCount, tracer);
Lynch, Vickie's avatar
Lynch, Vickie committed
268
    if (m_edge > 0) {
Lynch, Vickie's avatar
Lynch, Vickie committed
269
      if (edgePixel(inst, p->getBankName(), p->getCol(), p->getRow(), m_edge))
Lynch, Vickie's avatar
Lynch, Vickie committed
270
271
        return;
    }
272
273
    if (p->getDetectorID() != -1)
      peakWS->addPeak(*p);
274
275
  }
  catch (std::exception &e) {
276
    g_log.notice() << "Error creating peak at " << Q << " because of '"
277
                   << e.what() << "'. Peak will be skipped.\n";
278
  }
279
280
281
282
283
284
}

/**
 * Creates a Peak object from Q & bin count
 * */
boost::shared_ptr<DataObjects::Peak>
285
286
FindPeaksMD::createPeak(const Mantid::Kernel::V3D &Q, const double binCount,
                        const Geometry::InstrumentRayTracer &tracer) {
287
288
289
  boost::shared_ptr<DataObjects::Peak> p;
  if (dimType == QLAB) {
    // Build using the Q-lab-frame constructor
290
    p = boost::make_shared<Peak>(inst, Q);
291
    // Save gonio matrix for later
292
    p->setGoniometerMatrix(m_goniometer);
293
294
  } else if (dimType == QSAMPLE) {
    // Build using the Q-sample-frame constructor
295
296
297
298
299
300
301
302
303
304
305
306
307
308
    bool calcGoniometer = getProperty("CalculateGoniometerForCW");

    if (calcGoniometer) {
      // Calculate Q lab from Q sample and wavelength
      double wavelength = getProperty("Wavelength");
      double wv = 2.0 * M_PI / wavelength;
      double norm_q2 = Q.norm2();
      double theta = acos(1 - norm_q2 / (2 * wv * wv));
      double phi = asin(-Q[1] / wv * sin(theta));
      V3D Q_lab(-wv * sin(theta) * cos(phi), -wv * sin(theta) * sin(phi),
                wv * (1 - cos(theta)));

      // Solve to find rotation matrix, assuming only rotation around y-axis
      // A * X = B
309
      Matrix<double> A({ Q[0], Q[2], Q[2], -Q[0] }, 2, 2);
310
      A.Invert();
311
      std::vector<double> B{ Q_lab[0], Q_lab[2] };
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
      std::vector<double> X = A * B;
      double rot = atan2(X[1], X[0]);
      g_log.information() << "Found goniometer rotation to be "
                          << rot * 180 / M_PI
                          << " degrees for peak at Q sample = " << Q << "\n";

      Matrix<double> goniometer(3, 3, true);
      goniometer[0][0] = cos(rot);
      goniometer[0][2] = sin(rot);
      goniometer[2][0] = -sin(rot);
      goniometer[2][2] = cos(rot);
      p = boost::make_shared<Peak>(inst, Q, goniometer);

    } else {
      p = boost::make_shared<Peak>(inst, Q, m_goniometer);
    }
328
  } else {
329
    throw std::invalid_argument(
330
        "Cannot Integrate peaks unless the dimension is QLAB or QSAMPLE");
331
332
  }

333
  try { // Look for a detector
334
    p->findDetector(tracer);
335
336
  }
  catch (...) { /* Ignore errors in ray-tracer */
337
338
  }

339
340
  p->setBinCount(binCount);
  // Save the run number found before.
341
  p->setRunNumber(m_runNumber);
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
  return p;
}

//----------------------------------------------------------------------------------------------
/** Integrate the peaks of the workspace using parameters saved in the algorithm
 * class
 * @param ws ::  MDEventWorkspace to integrate
 */
template <typename MDE, size_t nd>
void FindPeaksMD::findPeaks(typename MDEventWorkspace<MDE, nd>::sptr ws) {
  if (nd < 3)
    throw std::invalid_argument("Workspace must have at least 3 dimensions.");

  if (isFullMDEvent<MDE, nd>()) {
    m_addDetectors = true;
  } else {
    m_addDetectors = false;
    g_log.warning("Workspace contains only lean events. Resultant "
                  "PeaksWorkspaces will not contain full detector "
                  "information.");
  }
363

364
365
366
367
368
369
  progress(0.01, "Refreshing Centroids");

  // TODO: This might be slow, progress report?
  // Make sure all centroids are fresh
  // ws->getBox()->refreshCentroid();

370
371
372
  uint16_t nexp = ws->getNumExperimentInfo();

  if (nexp == 0)
373
374
    throw std::runtime_error(
        "No instrument was found in the MDEventWorkspace. Cannot find peaks.");
375
376
377

  for (uint16_t iexp = 0; iexp < ws->getNumExperimentInfo(); iexp++) {
    ExperimentInfo_sptr ei = ws->getExperimentInfo(iexp);
378
    this->readExperimentInfo(ei, boost::dynamic_pointer_cast<IMDWorkspace>(ws));
379
380

    Geometry::InstrumentRayTracer tracer(inst);
381
382
383
384
385
    // Copy the instrument, sample, run to the peaks workspace.
    peakWS->copyExperimentInfoFrom(ei.get());

    // Calculate a threshold below which a box is too diffuse to be considered a
    // peak.
386
    signal_t threshold =
Anton Piccardo-Selg's avatar
Anton Piccardo-Selg committed
387
388
389
        m_useNumberOfEventsNormalization
            ? ws->getBox()->getSignalByNEvents() * m_signalThresholdFactor
            : ws->getBox()->getSignalNormalized() * DensityThresholdFactor;
390
391
392
    threshold *= m_densityScaleFactor;

    if (!std::isfinite(threshold)) {
393
394
      g_log.warning()
          << "Infinite or NaN overall density found. Your input data "
395
             "may be invalid. Using a 0 threshold instead.\n";
396
      threshold = 0;
397
    }
398
    g_log.information() << "Threshold signal density: " << threshold << '\n';
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421

    typedef API::IMDNode *boxPtr;
    // We will fill this vector with pointers to all the boxes (up to a given
    // depth)
    typename std::vector<API::IMDNode *> boxes;

    // Get all the MDboxes
    progress(0.10, "Getting Boxes");
    ws->getBox()->getBoxes(boxes, 1000, true);

    // This pair is the <density, ptr to the box>
    typedef std::pair<double, API::IMDNode *> dens_box;

    // Map that will sort the boxes by increasing density. The key = density;
    // value = box *.
    typename std::multimap<double, API::IMDNode *> sortedBoxes;

    // --------------- Sort and Filter by Density -----------------------------
    progress(0.20, "Sorting Boxes by Density");
    auto it1 = boxes.begin();
    auto it1_end = boxes.end();
    for (; it1 != it1_end; it1++) {
      auto box = *it1;
422
423
424
      double value = m_useNumberOfEventsNormalization
                         ? box->getSignalByNEvents()
                         : box->getSignalNormalized();
425
426
427
428
      value *= m_densityScaleFactor;
      // Skip any boxes with too small a signal value.
      if (value > threshold)
        sortedBoxes.insert(dens_box(value, box));
429
    }
430

431
432
433
434
    // --------------- Find Peak Boxes -----------------------------
    // List of chosen possible peak boxes.
    std::vector<API::IMDNode *> peakBoxes;

435
    prog = make_unique<Progress>(this, 0.30, 0.95, m_maxPeaks);
436
437
438
439
440
441
442
443

    // used for selecting method for calculating BinCount
    bool isMDEvent(ws->id().find("MDEventWorkspace") != std::string::npos);

    int64_t numBoxesFound = 0;
    // Now we go (backwards) through the map
    // e.g. from highest density down to lowest density.
    typename std::multimap<double, boxPtr>::reverse_iterator it2;
444
    auto it2_end = sortedBoxes.rend();
445
446
447
    for (it2 = sortedBoxes.rbegin(); it2 != it2_end; it2++) {
      signal_t density = it2->first;
      boxPtr box = it2->second;
448
#ifndef MDBOX_TRACK_CENTROID
449
450
      coord_t boxCenter[nd];
      box->calculateCentroid(boxCenter);
451
#else
452
      const coord_t *boxCenter = box->getCentroid();
453
#endif
454
455
456

      // Compare to all boxes already picked.
      bool badBox = false;
457
      for (auto &peakBoxe : peakBoxes) {
458

459
#ifndef MDBOX_TRACK_CENTROID
460
461
        coord_t otherCenter[nd];
        (*it3)->calculateCentroid(otherCenter);
462
#else
463
        const coord_t *otherCenter = peakBoxe->getCentroid();
464
#endif
465
466
467
468
469
470
471
472
473
474
475
476
477

        // Distance between this box and a box we already put in.
        coord_t distSquared = 0.0;
        for (size_t d = 0; d < nd; d++) {
          coord_t dist = otherCenter[d] - boxCenter[d];
          distSquared += (dist * dist);
        }

        // Reject this box if it is too close to another previously found box.
        if (distSquared < peakRadiusSquared) {
          badBox = true;
          break;
        }
478
      }
479

480
481
      // The box was not rejected for another reason.
      if (!badBox) {
482
        if (numBoxesFound++ >= m_maxPeaks) {
483
          g_log.notice() << "Number of peaks found exceeded the limit of "
Hahn, Steven's avatar
Hahn, Steven committed
484
                         << m_maxPeaks << ". Stopping peak finding.\n";
485
486
487
488
          break;
        }

        peakBoxes.push_back(box);
489
        g_log.debug() << "Found box at ";
490
        for (size_t d = 0; d < nd; d++)
491
          g_log.debug() << (d > 0 ? "," : "") << boxCenter[d];
492
        g_log.debug() << "; Density = " << density << '\n';
493
494
        // Report progres for each box found.
        prog->report("Finding Peaks");
495
496
497
      }
    }

498
    prog->resetNumSteps(numBoxesFound, 0.95, 1.0);
499

500
    // --- Convert the "boxes" to peaks ----
501
    for (auto box : peakBoxes) {
502
      //  If no events from this experimental contribute to the box then skip
503
504
505
      if (nexp > 1) {
        MDBox<MDE, nd> *mdbox = dynamic_cast<MDBox<MDE, nd> *>(box);
        typename std::vector<MDE> &events = mdbox->getEvents();
506
507
        if (std::none_of(events.cbegin(), events.cend(),
                         [&iexp, &nexp](MDE event) {
Owen Arnold's avatar
Owen Arnold committed
508
509
              return event.getRunIndex() == iexp || event.getRunIndex() >= nexp;
            }))
510
511
          continue;
      }
Owen Arnold's avatar
Owen Arnold committed
512
// The center of the box = Q in the lab frame
513

514
#ifndef MDBOX_TRACK_CENTROID
515
516
      coord_t boxCenter[nd];
      box->calculateCentroid(boxCenter);
517
#else
518
      const coord_t *boxCenter = box->getCentroid();
519
#endif
520
521
522
523
524
525

      // Q of the centroid of the box
      V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);

      // The "bin count" used will be the box density or the number of events in
      // the box
526
      double binCount = box->getSignalNormalized() * m_densityScaleFactor;
527
528
529
530
      if (isMDEvent)
        binCount = static_cast<double>(box->getNPoints());

      try {
531
        auto p = this->createPeak(Q, binCount, tracer);
532
533
534
535
536
537
538
        if (m_addDetectors) {
          auto mdBox = dynamic_cast<MDBoxBase<MDE, nd> *>(box);
          if (!mdBox) {
            throw std::runtime_error("Failed to cast box to MDBoxBase");
          }
          addDetectors(*p, *mdBox);
        }
Zhou, Wenduo's avatar
Zhou, Wenduo committed
539
        if (p->getDetectorID() != -1) {
Lynch, Vickie's avatar
Lynch, Vickie committed
540
          if (m_edge > 0) {
Lynch, Vickie's avatar
Lynch, Vickie committed
541
542
            if (!edgePixel(inst, p->getBankName(), p->getCol(), p->getRow(),
                           m_edge))
Lynch, Vickie's avatar
Lynch, Vickie committed
543
544
545
546
547
              peakWS->addPeak(*p);
            ;
          } else {
            peakWS->addPeak(*p);
          }
Zhou, Wenduo's avatar
Zhou, Wenduo committed
548
549
550
          g_log.information() << "Add new peak with Q-center = " << Q[0] << ", "
                              << Q[1] << ", " << Q[2] << "\n";
        }
551
552
      }
      catch (std::exception &e) {
553
        g_log.notice() << "Error creating peak at " << Q << " because of '"
554
                       << e.what() << "'. Peak will be skipped.\n";
555
      }
556

557
558
      // Report progress for each box found.
      prog->report("Adding Peaks");
559

560
    } // for each box found
561
  }
562
  g_log.notice() << "Number of peaks found: " << peakWS->getNumberPeaks()
563
                 << '\n';
564
565
566
567
568
569
570
}

//----------------------------------------------------------------------------------------------
/** Find peaks in the given MDHistoWorkspace
 *
 * @param ws :: MDHistoWorkspace
 */
571
572
void
FindPeaksMD::findPeaksHisto(Mantid::DataObjects::MDHistoWorkspace_sptr ws) {
573
574
575
576
577
578
579
580
581
582
  size_t nd = ws->getNumDims();
  if (nd < 3)
    throw std::invalid_argument("Workspace must have at least 3 dimensions.");

  g_log.warning("Workspace is an MDHistoWorkspace. Resultant PeaksWorkspaces "
                "will not contain full detector information.");

  if (ws->getNumExperimentInfo() == 0)
    throw std::runtime_error(
        "No instrument was found in the workspace. Cannot find peaks.");
583
584
585

  for (uint16_t iexp = 0; iexp < ws->getNumExperimentInfo(); iexp++) {
    ExperimentInfo_sptr ei = ws->getExperimentInfo(iexp);
586
    this->readExperimentInfo(ei, boost::dynamic_pointer_cast<IMDWorkspace>(ws));
587
    Geometry::InstrumentRayTracer tracer(inst);
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609

    // Copy the instrument, sample, run to the peaks workspace.
    peakWS->copyExperimentInfoFrom(ei.get());

    // This pair is the <density, box index>
    typedef std::pair<double, size_t> dens_box;

    // Map that will sort the boxes by increasing density. The key = density;
    // value = box index.
    std::multimap<double, size_t> sortedBoxes;

    size_t numBoxes = ws->getNPoints();

    // --------- Count the overall signal density -----------------------------
    progress(0.10, "Counting Total Signal");
    double totalSignal = 0;
    for (size_t i = 0; i < numBoxes; i++)
      totalSignal += ws->getSignalAt(i);
    // Calculate the threshold density
    double thresholdDensity =
        (totalSignal * ws->getInverseVolume() / double(numBoxes)) *
        DensityThresholdFactor * m_densityScaleFactor;
610
    if (!std::isfinite(thresholdDensity)) {
611
612
      g_log.warning()
          << "Infinite or NaN overall density found. Your input data "
613
             "may be invalid. Using a 0 threshold instead.\n";
614
615
      thresholdDensity = 0;
    }
616
617
    g_log.information() << "Threshold signal density: " << thresholdDensity
                        << '\n';
618
619
620
621
622
623
624
625

    // -------------- Sort and Filter by Density -----------------------------
    progress(0.20, "Sorting Boxes by Density");
    for (size_t i = 0; i < numBoxes; i++) {
      double density = ws->getSignalNormalizedAt(i) * m_densityScaleFactor;
      // Skip any boxes with too small a signal density.
      if (density > thresholdDensity)
        sortedBoxes.insert(dens_box(density, i));
626
627
    }

628
629
630
631
    // --------------- Find Peak Boxes -----------------------------
    // List of chosen possible peak boxes.
    std::vector<size_t> peakBoxes;

632
    prog = make_unique<Progress>(this, 0.30, 0.95, m_maxPeaks);
633
634
635
636
637

    int64_t numBoxesFound = 0;
    // Now we go (backwards) through the map
    // e.g. from highest density down to lowest density.
    std::multimap<double, size_t>::reverse_iterator it2;
638
    auto it2_end = sortedBoxes.rend();
639
640
641
642
643
644
645
646
    for (it2 = sortedBoxes.rbegin(); it2 != it2_end; ++it2) {
      signal_t density = it2->first;
      size_t index = it2->second;
      // Get the center of the box
      VMD boxCenter = ws->getCenter(index);

      // Compare to all boxes already picked.
      bool badBox = false;
647
      for (auto &peakBoxe : peakBoxes) {
648
        VMD otherCenter = ws->getCenter(peakBoxe);
649
650
651
652
653
654
655
656
657
658
659
660
661

        // Distance between this box and a box we already put in.
        coord_t distSquared = 0.0;
        for (size_t d = 0; d < nd; d++) {
          coord_t dist = otherCenter[d] - boxCenter[d];
          distSquared += (dist * dist);
        }

        // Reject this box if it is too close to another previously found box.
        if (distSquared < peakRadiusSquared) {
          badBox = true;
          break;
        }
662
      }
663

664
665
      // The box was not rejected for another reason.
      if (!badBox) {
666
        if (numBoxesFound++ >= m_maxPeaks) {
667
          g_log.notice() << "Number of peaks found exceeded the limit of "
Hahn, Steven's avatar
Hahn, Steven committed
668
                         << m_maxPeaks << ". Stopping peak finding.\n";
669
670
671
672
673
          break;
        }

        peakBoxes.push_back(index);
        g_log.debug() << "Found box at index " << index;
674
        g_log.debug() << "; Density = " << density << '\n';
675
676
677
        // Report progres for each box found.
        prog->report("Finding Peaks");
      }
678
    }
679
    // --- Convert the "boxes" to peaks ----
680
    for (auto index : peakBoxes) {
681
682
      // The center of the box = Q in the lab frame
      VMD boxCenter = ws->getCenter(index);
683

684
685
      // Q of the centroid of the box
      V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
686

687
688
      // The "bin count" used will be the box density.
      double binCount = ws->getSignalNormalizedAt(index) * m_densityScaleFactor;
689

690
      // Create the peak
691
      addPeak(Q, binCount, tracer);
692

693
694
      // Report progres for each box found.
      prog->report("Adding Peaks");
695

696
697
    } // for each box found
  }
698
  g_log.notice() << "Number of peaks found: " << peakWS->getNumberPeaks()
699
                 << '\n';
700
}
701
702
703
704
705

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
 */
void FindPeaksMD::exec() {
Owen Arnold's avatar
Owen Arnold committed
706

707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
  bool AppendPeaks = getProperty("AppendPeaks");

  // Output peaks workspace, create if needed
  peakWS = getProperty("OutputWorkspace");
  if (!peakWS || !AppendPeaks)
    peakWS = PeaksWorkspace_sptr(new PeaksWorkspace());

  // The MDEventWorkspace as input
  IMDWorkspace_sptr inWS = getProperty("InputWorkspace");
  MDHistoWorkspace_sptr inMDHW =
      boost::dynamic_pointer_cast<MDHistoWorkspace>(inWS);
  IMDEventWorkspace_sptr inMDEW =
      boost::dynamic_pointer_cast<IMDEventWorkspace>(inWS);

  // Other parameters
  double PeakDistanceThreshold = getProperty("PeakDistanceThreshold");
  peakRadiusSquared =
      static_cast<coord_t>(PeakDistanceThreshold * PeakDistanceThreshold);

  DensityThresholdFactor = getProperty("DensityThresholdFactor");
727
  m_signalThresholdFactor = getProperty("SignalThresholdFactor");
728
729
730

  std::string strategy = getProperty("PeakFindingStrategy");
  m_useNumberOfEventsNormalization = strategy == numberOfEventsNormalization;
731

732
  m_maxPeaks = getProperty("MaxPeaks");
733
  m_edge = this->getProperty("EdgePixels");
734
735
736
737
738
739
740
741
742
743
744
745
746

  // Execute the proper algo based on the type of workspace
  if (inMDHW) {
    this->findPeaksHisto(inMDHW);
  } else if (inMDEW) {
    CALL_MDEVENT_FUNCTION3(this->findPeaks, inMDEW);
  } else {
    throw std::runtime_error("This algorithm can only find peaks on a "
                             "MDHistoWorkspace or a MDEventWorkspace; it does "
                             "not work on a regular MatrixWorkspace.");
  }

  // Do a sort by bank name and then descending bin count (intensity)
747
  std::vector<std::pair<std::string, bool> > criteria;
748
  criteria.push_back(std::pair<std::string, bool>("RunNumber", true));
749
750
751
  criteria.push_back(std::pair<std::string, bool>("BankName", true));
  criteria.push_back(std::pair<std::string, bool>("bincount", false));
  peakWS->sort(criteria);
752

753
  for (auto i = 0; i != peakWS->getNumberPeaks(); ++i) {
754
    Peak &p = peakWS->getPeak(i);
Lynch, Vickie's avatar
Lynch, Vickie committed
755
    p.setPeakNumber(static_cast<int>(i + 1));
756
  }
757
758
759
  // Save the output
  setProperty("OutputWorkspace", peakWS);
}
760

Anton Piccardo-Selg's avatar
Anton Piccardo-Selg committed
761
762
763
764
765
766
//----------------------------------------------------------------------------------------------
/** Validate the inputs.
 */
std::map<std::string, std::string> FindPeaksMD::validateInputs() {
  std::map<std::string, std::string> result;
  // Check for number of event normalzation
767
768
  std::string strategy = getProperty("PeakFindingStrategy");

769
770
  const bool useNumberOfEventsNormalization =
      strategy == numberOfEventsNormalization;
Anton Piccardo-Selg's avatar
Anton Piccardo-Selg committed
771
772
773
774
775
  IMDWorkspace_sptr inWS = getProperty("InputWorkspace");
  IMDEventWorkspace_sptr inMDEW =
      boost::dynamic_pointer_cast<IMDEventWorkspace>(inWS);

  if (useNumberOfEventsNormalization && !inMDEW) {
776
777
778
    result["PeakFindingStrategy"] = "The NumberOfEventsNormalization selection "
                                    "can only be used with an MDEventWorkspace "
                                    "as the input.";
Anton Piccardo-Selg's avatar
Anton Piccardo-Selg committed
779
  }
780
781
782
783
784
785

  double wavelength = getProperty("Wavelength");
  if (getProperty("CalculateGoniometerForCW") && wavelength == DBL_MAX)
    result["Wavelength"] =
        "Must set wavelength when using CalculateGoniometerForCW option";

Anton Piccardo-Selg's avatar
Anton Piccardo-Selg committed
786
787
788
  return result;
}

789
} // namespace MDAlgorithms
790
} // namespace Mantid