Mantid
Loading...
Searching...
No Matches
IntegratePeaksHybrid.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 +
7/*WIKI*
8
9 Integrates arbitrary shaped single crystal peaks defined on an
10 [[MDHistoWorkspace]] using connected component analysis to determine
11 regions of interest around each peak of the [[PeaksWorkspace]]. The output is
12 an integrated [[PeaksWorkspace]] as well as an image
13 containing the labels assigned to each cluster for diagnostic and visualisation
14 purposes.
15
16 This algorithm is very similar to, [[IntegratePeaksUsingClusters]] but breaks
17 the integration into a series of local image domains rather than integrating an
18 image in one-shot.
19 The advantages of this approach are that you can locally define a background
20 rather than using a global setting, and are therefore better able to capture
21 the peak shape. A further advantage
22 is that the memory requirement is reduced, given that [[MDHistoWorkspaces]] are
23 generated in the region of the peak, and therefore high resolution can be
24 achieved in the region of the peaks without
25 an overall high n-dimensional image cost.
26
27 Unlike [[IntegratePeaksUsingClusters]] you do not need to specify at Threshold
28 for background detection. You do however need to specify a
29 BackgroundOuterRadius in a similar fashion to
30 [[IntegratePeaksMD]]. This is used to determine the region in which to make an
31 [[MDHistoWorkspace]] around each peak. A liberal estimate is a good idea.
32
33 At present, you will need to provide NumberOfBins for binning (via [[BinMD]]
34 axis-aligned). By default, the algorithm will create a 20 by 20 by 20 grid
35 around each peak. The same
36 number of bins is applied in each dimension.
37
38 == Warnings and Logging ==
39 See [[IntegratePeaksUsingClusters]] for notes on logs and warnings.
40
41 *WIKI*/
42
44
49
59
60#include <boost/format.hpp>
61#include <cmath>
62
63using namespace Mantid::API;
64using namespace Mantid::Kernel;
65using namespace Mantid::DataObjects;
66
67namespace {
76std::string extractFormattedPropertyFromDimension(Mantid::Geometry::IMDDimension_const_sptr &dimension,
77 const double &min, const double &max, const int &nBins) {
78 std::string id = dimension->getDimensionId();
79 return boost::str(boost::format("%s, %f, %f, %d") % id % min % max % nBins);
80}
81} // namespace
82
83namespace Mantid::Crystal {
84using namespace ConnectedComponentMappingTypes;
85
86// Register the algorithm into the AlgorithmFactory
87DECLARE_ALGORITHM(IntegratePeaksHybrid)
88
89//----------------------------------------------------------------------------------------------
91const std::string IntegratePeaksHybrid::name() const { return "IntegratePeaksHybrid"; }
92
94int IntegratePeaksHybrid::version() const { return 1; }
95
97const std::string IntegratePeaksHybrid::category() const { return "MDAlgorithms\\Peaks;Crystal\\Integration"; }
98
99//----------------------------------------------------------------------------------------------
103 declareProperty(std::make_unique<WorkspaceProperty<IMDEventWorkspace>>("InputWorkspace", "", Direction::Input),
104 "Input md workspace.");
105 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeaksWorkspace", "", Direction::Input),
106 "A PeaksWorkspace containing the peaks to integrate.");
107
108 auto positiveIntValidator = std::make_shared<BoundedValidator<int>>();
109 positiveIntValidator->setExclusive(true);
110 positiveIntValidator->setLower(0);
111
112 declareProperty(std::make_unique<PropertyWithValue<int>>("NumberOfBins", 20, positiveIntValidator, Direction::Input),
113 "Number of bins to use while creating each local image. "
114 "Defaults to 20. Increase to reduce pixelation");
115
116 auto compositeValidator = std::make_shared<CompositeValidator>();
117 auto positiveDoubleValidator = std::make_shared<BoundedValidator<double>>();
118 positiveDoubleValidator->setExclusive(true);
119 positiveDoubleValidator->setLower(0);
120 compositeValidator->add(positiveDoubleValidator);
121 compositeValidator->add(std::make_shared<MandatoryValidator<double>>());
122
124 std::make_unique<PropertyWithValue<double>>("BackgroundOuterRadius", 0.0, compositeValidator, Direction::Input),
125 "Background outer radius estimate. Choose liberal value.");
126
127 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
128 "An output integrated peaks workspace.");
129
130 declareProperty(std::make_unique<WorkspaceProperty<WorkspaceGroup>>("OutputWorkspaces", "", Direction::Output),
131 "MDHistoWorkspaces containing the labeled clusters used by "
132 "the algorithm.");
133}
134
135const std::string IntegratePeaksHybrid::summary() const {
136 return "Integrate single crystal peaks using connected component analysis. "
137 "Binning invididual to each peak.";
138}
139
140//----------------------------------------------------------------------------------------------
144 IMDEventWorkspace_sptr mdWS = getProperty("InputWorkspace");
145 PeaksWorkspace_sptr inPeakWS = getProperty("PeaksWorkspace");
146 PeaksWorkspace_sptr peakWS = getProperty("OutputWorkspace");
147 const int numBins = getProperty("NumberOfBins");
148 const double peakOuterRadius = getProperty("BackgroundOuterRadius");
149 const double halfPeakOuterRadius = peakOuterRadius / 2;
150 if (peakWS != inPeakWS) {
151 peakWS = inPeakWS->clone();
152 }
153
154 {
155 const SpecialCoordinateSystem mdCoordinates = mdWS->getSpecialCoordinateSystem();
156 if (mdCoordinates == None) {
157 throw std::invalid_argument("The coordinate system of the input "
158 "MDWorkspace cannot be established. Create "
159 "your workspace with the an MDFrame which is "
160 "not a General Frame or Unknown Frame.");
161 }
162 }
163
164 PeakClusterProjection projection(mdWS);
165 auto outImageResults = std::make_shared<WorkspaceGroup>();
166
167 Progress progress(this, 0.0, 1.0, peakWS->getNumberPeaks());
168
170 for (int i = 0; i < peakWS->getNumberPeaks(); ++i) {
171
173 Geometry::IPeak &peak = peakWS->getPeak(i);
174 const V3D center = projection.peakCenter(peak);
175
176 auto binMDAlg = this->createChildAlgorithm("BinMD");
177 binMDAlg->setProperty("InputWorkspace", mdWS);
178 binMDAlg->setPropertyValue("OutputWorkspace", "output_ws");
179 binMDAlg->setProperty("AxisAligned", true);
180
181 for (int j = 0; j < static_cast<int>(mdWS->getNumDims()); ++j) {
182 std::stringstream propertyName;
183 propertyName << "AlignedDim" << j;
184
185 auto dimension = mdWS->getDimension(j);
186
187 double min = center[j] - halfPeakOuterRadius;
188 double max = center[j] + halfPeakOuterRadius;
189
190 binMDAlg->setPropertyValue(propertyName.str(),
191 extractFormattedPropertyFromDimension(dimension, min, max, numBins));
192 }
193 binMDAlg->execute();
194 Workspace_sptr temp = binMDAlg->getProperty("OutputWorkspace");
195 IMDHistoWorkspace_sptr localImage = std::dynamic_pointer_cast<IMDHistoWorkspace>(temp);
197 auto iterator = localImage->createIterator();
198 iterator->setNormalization(normalization);
199 signal_t cumulative = 0;
200 do {
201 cumulative += iterator->getSignal();
202 } while (iterator->next());
203 const double threshold = cumulative / signal_t(localImage->getNPoints());
204
205 HardThresholdBackground backgroundStrategy(threshold, normalization);
206 // CCL. Multi-processor version.
207 const size_t startId = 1;
208 const size_t nThreads = 1;
209 ConnectedComponentLabeling analysis(startId, nThreads); // CCL executed single threaded.
210
211 Progress dummyProgress;
212 // Perform CCL.
213 ClusterTuple clusters = analysis.executeAndFetchClusters(localImage, &backgroundStrategy, dummyProgress);
214 // Extract the clusters
215 ConnectedComponentMappingTypes::ClusterMap &clusterMap = clusters.get<1>();
216 // Extract the labeled image
217 IMDHistoWorkspace_sptr outHistoWS = clusters.get<0>();
218 outImageResults->addWorkspace(outHistoWS);
219
220 PeakClusterProjection localProjection(outHistoWS);
221 const Mantid::signal_t signalValue =
222 localProjection.signalAtPeakCenter(peak); // No normalization when extracting label ids!
223
224 if (std::isnan(signalValue)) {
225 g_log.warning() << "Warning: image for integration is off edge of detector for peak " << i << '\n';
226 } else if (signalValue < static_cast<Mantid::signal_t>(analysis.getStartLabelId())) {
227 g_log.information() << "Peak: " << i
228 << " Has no corresponding cluster/blob detected on "
229 "the image. This could be down to your Threshold "
230 "settings.\n";
231 } else {
232 const auto labelIdAtPeak = static_cast<size_t>(signalValue);
233 ICluster *const cluster = clusterMap[labelIdAtPeak].get();
234 ICluster::ClusterIntegratedValues integratedValues = cluster->integrate(localImage);
235 peak.setIntensity(integratedValues.get<0>());
236 peak.setSigmaIntensity(std::sqrt(integratedValues.get<1>()));
237 }
238 progress.report();
240 }
242
243 setProperty("OutputWorkspace", peakWS);
244 setProperty("OutputWorkspaces", outImageResults);
245}
246
247} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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.
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
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
Definition: Algorithm.cpp:842
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
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
ConnectedComponentLabelling : Implements connected component labeling on MDHistoWorkspaces.
size_t getStartLabelId() const
Getter for the start label id.
ConnectedComponentMappingTypes::ClusterTuple executeAndFetchClusters(Mantid::API::IMDHistoWorkspace_sptr ws, BackgroundStrategy *const strategy, Mantid::API::Progress &progress) const
Execute and return clusters, as well as maps to integratable clusters.
HardThresholdBackground : Implementation of BackgroundStrategy using a fixed background signal value ...
ICluster : Abstract cluster.
Definition: ICluster.h:25
boost::tuple< double, double > ClusterIntegratedValues
Definition: ICluster.h:27
virtual ClusterIntegratedValues integrate(std::shared_ptr< const Mantid::API::IMDHistoWorkspace > ws) const =0
integrate the cluster
IntegratePeaksHybrid : Integrate single crystal peaks algorithm.
const std::string summary() const override
function returns a summary message that will be displayed in the default GUI, and in the help.
void exec() override
Execute the algorithm.
void init() override
Initialize the algorithm's properties.
int version() const override
Algorithm's version for identification.
const std::string category() const override
Algorithm's category for identification.
PeakClusterProjection : Maps peaks onto IMDHistoWorkspaces and returns the signal value at the peak c...
Mantid::signal_t signalAtPeakCenter(const Mantid::Geometry::IPeak &peak, Mantid::API::MDNormalization normalization=Mantid::API::NoNormalization) const
Get the signal value at the peak center.
Mantid::Kernel::V3D peakCenter(const Mantid::Geometry::IPeak &peak) const
Get the peak center.
Structure describing a single-crystal peak.
Definition: IPeak.h:26
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Validator to check that a property is not left empty.
The concrete, templated class for properties.
Class for 3D vectors.
Definition: V3D.h:34
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< IMDHistoWorkspace > IMDHistoWorkspace_sptr
shared pointer to Mantid::API::IMDHistoWorkspace
MDNormalization
Enum describing different ways to normalize the signal in a MDWorkspace.
Definition: IMDIterator.h:25
@ NoNormalization
Don't normalize = return raw counts.
Definition: IMDIterator.h:27
std::map< size_t, std::shared_ptr< Mantid::Crystal::ICluster > > ClusterMap
boost::tuple< Mantid::API::IMDHistoWorkspace_sptr, ClusterMap > ClusterTuple
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< const IMDDimension > IMDDimension_const_sptr
Shared Pointer to const IMDDimension.
Definition: IMDDimension.h:101
SpecialCoordinateSystem
Special coordinate systems for Q3D.
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
double signal_t
Typedef for the signal recorded in a MDBox, etc.
Definition: MDTypes.h:36
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54