35 size_t maxNeighbours = 3;
36 for (
size_t i = 1; i < ndims; ++i) {
51size_t calculateMaxClusters(
IMDHistoWorkspace const *
const ws,
size_t nIterators) {
52 size_t maxClusters = 1;
53 for (
size_t i = 0; i < ws->
getNumDims(); ++i) {
56 maxClusters /= nIterators;
57 if (maxClusters == 0) {
73 for (
int i = 0; i < static_cast<
int>(outWS->getNPoints()); ++i) {
74 outWS->setSignalAt(i, 0);
75 outWS->setErrorSquaredAt(i, 0);
87template <
typename T> T reportEvery(
const T &maxReports,
const T &maxIterations) {
88 T frequency = maxReports;
89 if (maxIterations >= maxReports) {
90 frequency = maxIterations / maxReports;
96template <
typename Container>
97bool does_contain_key(
const Container &container,
const typename Container::key_type &
value) {
98 return container.find(
value) != container.end();
101using EdgeIndexPair = boost::tuple<size_t, size_t>;
102using VecEdgeIndexPair = std::vector<EdgeIndexPair>;
122 size_t startLabelId, VecEdgeIndexPair &edgeIndexVec) {
123 size_t currentLabelCount = startLabelId;
133 nonEmptyNeighbourIndexes.reserve(maxNeighbours);
136 for (
auto neighIndex : neighbourIndexes) {
145 edgeIndexVec.emplace_back(currentIndex, neighIndex);
149 const DisjointElement &neighbourElement = neighbourElements[neighIndex];
151 if (!neighbourElement.
isEmpty()) {
152 nonEmptyNeighbourIndexes.emplace_back(neighIndex);
153 neighbourIds.insert(neighbourElement.
getId());
157 if (nonEmptyNeighbourIndexes.empty()) {
159 element.
setId(
static_cast<int>(currentLabelCount));
161 }
else if (neighbourIds.size() == 1)
163 neighbourElements[currentIndex] = neighbourElements[nonEmptyNeighbourIndexes.front()];
168 size_t candidateSourceParentIndex = nonEmptyNeighbourIndexes[0];
169 for (
size_t i = 1; i < nonEmptyNeighbourIndexes.size(); ++i) {
170 size_t neighIndex = nonEmptyNeighbourIndexes[i];
171 if (neighbourElements[neighIndex].getRoot() < neighbourElements[candidateSourceParentIndex].getRoot()) {
172 candidateSourceParentIndex = neighIndex;
176 DisjointElement &parentElement = neighbourElements[candidateSourceParentIndex];
178 for (
auto neighIndex : nonEmptyNeighbourIndexes) {
179 if (neighIndex != candidateSourceParentIndex) {
180 neighbourElements[neighIndex].
unionWith(&parentElement);
184 neighbourElements[currentIndex].unionWith(&parentElement);
187 }
while (iterator->
next());
189 return currentLabelCount;
194void memoryCheck(
size_t nPoints) {
195 size_t sizeOfElement = (3 *
sizeof(
signal_t)) +
sizeof(bool);
198 const size_t freeMemory = memoryStats.
availMem();
199 const size_t memoryCost = sizeOfElement * nPoints / 1000;
200 if (memoryCost > freeMemory) {
201 std::string basicMessage =
"CCL requires more free memory than you have available.";
202 std::stringstream sstream;
203 sstream << basicMessage <<
" Requires " << memoryCost <<
" KB of contiguous memory.";
205 throw std::runtime_error(basicMessage);
216 : m_startId(startId), m_nThreads(
std::move(nThreads)) {
218 throw std::invalid_argument(
"Cannot request that CCL runs with less than one thread!");
265 std::map<size_t, std::shared_ptr<ICluster>> clusterMap;
268 const size_t maxNeighbours = calculateMaxNeighbours(ws.get());
270 progress.
doReport(
"Identifying clusters");
271 auto frequency = reportEvery<size_t>(10000, ws->getNPoints());
278 if (nThreadsToUse > 1) {
279 auto iterators = ws->createIterators(nThreadsToUse);
280 const size_t maxClustersPossible = calculateMaxClusters(ws.get(), nThreadsToUse);
282 std::vector<VecEdgeIndexPair> parallelEdgeVec(nThreadsToUse);
284 std::vector<std::map<size_t, std::shared_ptr<Cluster>>> parallelClusterMapVec(nThreadsToUse);
289 for (
int i = 0; i < nThreadsToUse; ++i) {
291 boost::scoped_ptr<BackgroundStrategy> strategy(baseStrategy->
clone());
292 VecEdgeIndexPair &edgeVec = parallelEdgeVec[i];
294 const size_t startLabel =
m_startId + (i * maxClustersPossible);
297 const size_t endLabel = doConnectedComponentLabeling(iterator, strategy.get(), neighbourElements, progress,
298 maxNeighbours, startLabel, edgeVec);
301 std::map<size_t, std::shared_ptr<Cluster>> &localClusterMap = parallelClusterMapVec[i];
302 for (
size_t labelId = startLabel; labelId != endLabel; ++labelId) {
303 auto cluster = std::make_shared<Cluster>(labelId);
304 localClusterMap[labelId] = cluster;
315 const size_t &labelAtIndex = neighbourElements[currentIndex].getRoot();
316 localClusterMap[labelAtIndex]->addIndex(currentIndex);
318 }
while (iterator->
next());
325 for (
const auto ¶llelClusterMap : parallelClusterMapVec) {
326 for (
const auto &cluster : parallelClusterMap) {
327 clusterRegister.
add(cluster.first, cluster.second);
333 g_log.
debug(
"Percolate minimum label across boundaries");
335 for (
auto &indexPairVec : parallelEdgeVec) {
336 for (
auto &iit : indexPairVec) {
339 clusterRegister.
merge(a, b);
342 clusterMap = clusterRegister.
clusters(neighbourElements);
345 auto iterator = ws->createIterator(
nullptr);
346 VecEdgeIndexPair edgeIndexPair;
348 size_t endLabelId = doConnectedComponentLabeling(iterator.get(), baseStrategy, neighbourElements, progress,
349 maxNeighbours,
m_startId, edgeIndexPair);
352 for (
size_t labelId =
m_startId; labelId != endLabelId; ++labelId) {
353 auto cluster = std::make_shared<Cluster>(labelId);
354 clusterMap[labelId] = cluster;
363 const int labelAtIndex = neighbourElements[currentIndex].getRoot();
364 clusterMap[labelAtIndex]->addIndex(currentIndex);
366 }
while (iterator->
next());
398 memoryCheck(ws->getNPoints());
407 g_log.
debug(
"Finish cloning input workspace");
411 keys.reserve(clusters.size());
412 std::transform(clusters.cbegin(), clusters.cend(), std::back_inserter(keys),
413 [](
const auto &cluster) { return cluster.first; });
416 for (
int i = 0; i < static_cast<int>(keys.size()); ++i) {
417 clusters[keys[i]]->writeTo(outWS);
double value
The value of the point.
#define PARALLEL_FOR_NO_WSP_CHECK()
Abstract interface to MDHistoWorkspace, for use in exposing to Python.
This is an interface to an iterator of an IMDWorkspace.
virtual std::vector< size_t > findNeighbourIndexes() const =0
Find neighbouring indexes vertex touching.
virtual bool isWithinBounds(size_t index) const =0
Is index reachable by the iterator.
virtual bool next()=0
Advance to the next cell.
virtual size_t getLinearIndex() const =0
Get the linear index.
virtual void jumpTo(size_t index)=0
Jump to the index^th cell.
virtual uint64_t getNPoints() const =0
Get the number of points associated with the workspace.
virtual std::shared_ptr< const Mantid::Geometry::IMDDimension > getDimension(size_t index) const
Get a dimension.
virtual size_t getNumDims() const
Helper class for reporting progress from algorithms.
void doReport(const std::string &msg="") override
Actually do the reporting, without changing the loop counter.
BackgroundStrategy : Abstract class used for identifying elements of a IMDWorkspace that are not cons...
virtual void configureIterator(Mantid::API::IMDIterator *const iterator) const =0
virtual BackgroundStrategy * clone() const =0
virtual bool isBackground(Mantid::API::IMDIterator *const iterator) const =0
ClusterRegister : A fly-weight ICluster regeister.
void add(const size_t &label, const std::shared_ptr< ICluster > &cluster)
Add clusters.
MapCluster clusters(std::vector< DisjointElement > &elements) const
Get all combined clusters.
void merge(const DisjointElement &a, const DisjointElement &b) const
Merge clusters on the basis of known pairs of disjoint elements.
const boost::optional< int > m_nThreads
Run multithreaded.
void startLabelingId(const size_t &id)
Setter for the label id.
size_t getStartLabelId() const
Getter for the start label id.
ConnectedComponentLabeling(const size_t &startId=1, boost::optional< int > nThreads=boost::none)
Constructor.
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.
std::shared_ptr< Mantid::API::IMDHistoWorkspace > execute(Mantid::API::IMDHistoWorkspace_sptr ws, BackgroundStrategy *const strategy, Mantid::API::Progress &progress) const
Execute and return clusters.
ConnectedComponentMappingTypes::ClusterMap calculateDisjointTree(const Mantid::API::IMDHistoWorkspace_sptr &ws, BackgroundStrategy *const baseStrategy, Mantid::API::Progress &progress) const
Calculate the disjoint element tree across the image.
virtual ~ConnectedComponentLabeling()
Destructor.
size_t m_startId
Start labeling index.
int getNThreads() const
Get the number of threads to use.
DisjointElement : Cluster item used in a disjoint-set data structure.
bool isEmpty() const
Is empty.
void unionWith(DisjointElement *other)
Union with other.
void setId(int id)
Set the id.
Models a Container is used to hold a sample in the beam.
The Logger class is in charge of the publishing messages from the framework through various channels.
void debug(const std::string &msg)
Logs at debug level.
void notice(const std::string &msg)
Logs at notice level.
This class is responsible for memory statistics.
std::size_t availMem() const
Returns the available memory of the system in kiB.
void resetNumSteps(int64_t nsteps, double start, double end)
Change the number of steps between start/end.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< IMDHistoWorkspace > IMDHistoWorkspace_sptr
shared pointer to Mantid::API::IMDHistoWorkspace
Namespace containing useful typedefs.
std::map< size_t, std::shared_ptr< Mantid::Crystal::ICluster > > ClusterMap
std::unordered_set< size_t > SetIds
std::vector< size_t > VecIndexes
std::vector< DisjointElement > VecElements
boost::tuple< Mantid::API::IMDHistoWorkspace_sptr, ClusterMap > ClusterTuple
double signal_t
Typedef for the signal recorded in a MDBox, etc.