Skip to content

Commit

Permalink
refs #9194. Try replace custom DisjointElement
Browse files Browse the repository at this point in the history
  • Loading branch information
@OwenArnold
OwenArnold committed Apr 1, 2014
1 parent 8c71a11 commit 7b2bd12
Show file tree
Hide file tree
Showing 6 changed files with 206 additions and 310 deletions.
11 changes: 8 additions & 3 deletions Code/Mantid/Framework/Crystal/inc/MantidCrystal/ConnectedComponentLabeling.h
View file
Expand Up @@ -18,8 +18,11 @@ namespace Crystal
typedef std::map<size_t, SignalErrorSQPair > LabelIdIntensityMap;
typedef std::map<Mantid::Kernel::V3D, size_t> PositionToLabelIdMap;
typedef std::vector<size_t> VecIndexes;
typedef std::vector<DisjointElement> VecElements;
typedef std::vector<size_t> VecElements;
typedef std::set<size_t> SetIds;
typedef typename VecElements::value_type ElementType;
typedef std::map<ElementType,std::size_t> RankTypeMap;
typedef std::map<ElementType, ElementType> ElementParentMap;
}

class BackgroundStrategy;
Expand Down Expand Up @@ -72,9 +75,11 @@ namespace Crystal
int getNThreads() const;
/// Calculate the disjoint element tree across the image.
void calculateDisjointTree(Mantid::API::IMDHistoWorkspace_sptr ws,
BackgroundStrategy * const strategy, std::vector<DisjointElement>& neighbourElements,
BackgroundStrategy * const strategy, ConnectedComponentMappingTypes::VecElements& neighbourElements,
ConnectedComponentMappingTypes::LabelIdIntensityMap& labelMap,
ConnectedComponentMappingTypes::PositionToLabelIdMap& positionLabelMap) const;
ConnectedComponentMappingTypes::PositionToLabelIdMap& positionLabelMap,
ConnectedComponentMappingTypes::ElementParentMap& parentMap
) const;

/// Start labeling index
size_t m_startId;
Expand Down
23 changes: 23 additions & 0 deletions Code/Mantid/Framework/Crystal/inc/MantidCrystal/DisjointElement.h
View file
Expand Up @@ -57,6 +57,24 @@ namespace Crystal
DisjointElement(const DisjointElement& other);
/// Assignment operator.
DisjointElement& operator=(const DisjointElement& other);

bool operator==(const DisjointElement& other)
{
return other.getId() == this->getId();
}

bool operator!=(const DisjointElement& other)
{
return !(*this==other);
}

bool operator <(const DisjointElement& rhs)
{
return this->getId() < rhs.getId();
}
// friend bool operator< ( DisjointElement &a, DisjointElement &b);


private:
bool hasParent() const;
int compress();
Expand All @@ -72,6 +90,11 @@ namespace Crystal

};


bool operator< (const DisjointElement &a,const DisjointElement &b);



void unionElements(DisjointElement* a, DisjointElement* b);

} // namespace Crystal
Expand Down
146 changes: 83 additions & 63 deletions Code/Mantid/Framework/Crystal/src/ConnectedComponentLabeling.cpp
View file
Expand Up @@ -6,7 +6,7 @@
#include "MantidAPI/IMDIterator.h"
#include "MantidCrystal/ConnectedComponentLabeling.h"
#include "MantidCrystal/BackgroundStrategy.h"
#include "MantidCrystal/DisjointElement.h"
#include <boost/pending/disjoint_sets.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/tuple/tuple.hpp>
#include <boost/scoped_ptr.hpp>
Expand All @@ -27,7 +27,7 @@ namespace Mantid
{
const size_t ndims = ws->getNumDims();
size_t maxNeighbours = 3;
for(size_t i = 1; i < ndims; ++i)
for (size_t i = 1; i < ndims; ++i)
{
maxNeighbours *= 3;
}
Expand All @@ -50,37 +50,50 @@ namespace Mantid
}
return outWS;
}

template<typename Rank, typename Parent>
boost::disjoint_sets<Rank, Parent>* makeDisjointSet(Rank& r, Parent& p, VecElements& elements)
{
boost::disjoint_sets<Rank, Parent>* dsets = new boost::disjoint_sets<Rank, Parent>(r, p);
for (auto e = elements.begin(); e != elements.end(); e++)
{
dsets->make_set(*e);
}
return dsets;
}
}

//----------------------------------------------------------------------------------------------
/** Constructor
*/
ConnectedComponentLabeling::ConnectedComponentLabeling(const size_t& startId, const bool runMultiThreaded)
: m_startId(startId), m_runMultiThreaded(runMultiThreaded)
*/
ConnectedComponentLabeling::ConnectedComponentLabeling(const size_t& startId,
const bool runMultiThreaded) :
m_startId(startId), m_runMultiThreaded(runMultiThreaded)
{

}

/**
* Set a custom start id. This has no bearing on the output of the process other than
* the initial id used.
* @param id: Id to start with
*/
* Set a custom start id. This has no bearing on the output of the process other than
* the initial id used.
* @param id: Id to start with
*/
void ConnectedComponentLabeling::startLabelingId(const size_t& id)
{
m_startId = id;
}

/**
@return: The start label id.
*/
@return: The start label id.
*/
size_t ConnectedComponentLabeling::getStartLabelId() const
{
return m_startId;
}

//----------------------------------------------------------------------------------------------
/** Destructor
*/
*/
ConnectedComponentLabeling::~ConnectedComponentLabeling()
{
}
Expand All @@ -90,63 +103,67 @@ namespace Mantid
return m_runMultiThreaded ? API::FrameworkManager::Instance().getNumOMPThreads() : 1;
}

void ConnectedComponentLabeling::calculateDisjointTree(IMDHistoWorkspace_sptr ws,
BackgroundStrategy * const strategy, VecElements& neighbourElements,
LabelIdIntensityMap& labelMap,
PositionToLabelIdMap& positionLabelMap
) const
void ConnectedComponentLabeling::calculateDisjointTree(IMDHistoWorkspace_sptr ws,
BackgroundStrategy * const strategy, VecElements& neighbourElements,
LabelIdIntensityMap& labelMap, PositionToLabelIdMap& positionLabelMap, ElementParentMap& parentMap) const
{

// ---- Pre-process to eliminate those cells that are below background. Record the indexes of those that are above. ----
VecIndexes allNonBackgroundIndexes;
allNonBackgroundIndexes.reserve(ws->getNPoints());

if(m_runMultiThreaded)
if (m_runMultiThreaded)
{

std::vector<API::IMDIterator*> iterators = ws->createIterators(getNThreads());
const int nthreads = static_cast<int>(iterators.size());
std::vector<VecIndexes> manyNonBackgroundIndexes(nthreads);

PARALLEL_FOR_NO_WSP_CHECK()
for(int i = 0; i < nthreads; ++i)
for (int i = 0; i < nthreads; ++i)
{
boost::scoped_ptr<BackgroundStrategy> strategyCopy(strategy->clone());
API::IMDIterator *iterator = iterators[i];
VecIndexes& nonBackgroundIndexes = manyNonBackgroundIndexes[i];
do
{
boost::scoped_ptr<BackgroundStrategy> strategyCopy(strategy->clone());
API::IMDIterator *iterator = iterators[i];
VecIndexes& nonBackgroundIndexes = manyNonBackgroundIndexes[i];
do
if (!strategyCopy->isBackground(iterator))
{
if(!strategyCopy->isBackground(iterator))
{
nonBackgroundIndexes.push_back( iterator->getLinearIndex() );
}
nonBackgroundIndexes.push_back(iterator->getLinearIndex());
}
while(iterator->next());
}
// Consolidate work from individual threads.
for(size_t i = 0; i < nthreads; ++i)
{
VecIndexes& source = manyNonBackgroundIndexes[i];
allNonBackgroundIndexes.insert(allNonBackgroundIndexes.end(), source.begin(), source.end());
}
} while (iterator->next());
}
// Consolidate work from individual threads.
for (size_t i = 0; i < nthreads; ++i)
{
VecIndexes& source = manyNonBackgroundIndexes[i];
allNonBackgroundIndexes.insert(allNonBackgroundIndexes.end(), source.begin(), source.end());
}
}
else
{
API::IMDIterator *iterator = ws->createIterator(NULL);
do
{
if(!strategy->isBackground(iterator))
if (!strategy->isBackground(iterator))
{
allNonBackgroundIndexes.push_back( iterator->getLinearIndex() );
allNonBackgroundIndexes.push_back(iterator->getLinearIndex());
}
}
while(iterator->next());
} while (iterator->next());
}

// -------- Perform labeling -----------


RankTypeMap rankMap;
boost::associative_property_map<RankTypeMap> rankPmap(rankMap);
boost::associative_property_map<ElementParentMap> parentPmap(parentMap);


auto* disjointSet = makeDisjointSet(rankPmap, parentPmap, neighbourElements);

const size_t maxNeighbours = calculateMaxNeighbours(ws.get());
IMDIterator* iterator = ws->createIterator(NULL);
size_t currentLabelCount = m_startId;
for(size_t ii = 0; ii < allNonBackgroundIndexes.size(); ++ii)
for (size_t ii = 0; ii < allNonBackgroundIndexes.size(); ++ii)
{
size_t& currentIndex = allNonBackgroundIndexes[ii];
iterator->jumpTo(currentIndex);
Expand All @@ -171,9 +188,10 @@ namespace Mantid

if (nonEmptyNeighbourIndexes.empty())
{
neighbourElements[currentIndex] = DisjointElement(static_cast<int>(currentLabelCount)); // New leaf
neighbourElements[currentIndex] = currentLabelCount; // New leaf

labelMap[currentLabelCount] = 0; // Pre-fill the currentlabelcount.
const VMD& center = iterator->getCenter();
const VMD& center = iterator->getCenter();
positionLabelMap[V3D(center[0], center[1], center[2])] = currentLabelCount; // Get the position at this label.
++currentLabelCount;
}
Expand All @@ -183,44 +201,45 @@ namespace Mantid
}
else
{
// Choose the lowest neighbour index as the parent.
// Choose the lowest neighbour index as the parent. TODO: std::min instead!
size_t parentIndex = nonEmptyNeighbourIndexes[0];
for (size_t i = 1; i < nonEmptyNeighbourIndexes.size(); ++i)
{
size_t neighIndex = nonEmptyNeighbourIndexes[i];
if (neighbourElements[neighIndex].getId() < neighbourElements[parentIndex].getId())
if (neighbourElements[neighIndex] < neighbourElements[parentIndex])
{
parentIndex = i;
}
}
// Get the chosen parent
DisjointElement& parentElement = neighbourElements[parentIndex];

// Make this element a copy of the parent
neighbourElements[currentIndex] = parentElement;
neighbourElements[currentIndex] = neighbourElements[parentIndex]; // TODO. Not sure that this step is actually required.
// Union remainder parents with the chosen parent
for (size_t i = 0; i < nonEmptyNeighbourIndexes.size(); ++i)
{
size_t neighIndex = nonEmptyNeighbourIndexes[i];
if (neighIndex != parentIndex)
{
neighbourElements[neighIndex].unionWith(&parentElement);
disjointSet->union_set(neighbourElements[neighIndex], neighbourElements[parentIndex]);
}
}

}
}
disjointSet->compress_sets(neighbourElements.begin(), neighbourElements.end());
delete disjointSet;

}

boost::shared_ptr<Mantid::API::IMDHistoWorkspace> ConnectedComponentLabeling::execute(
IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy) const
IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy) const
{
VecElements neighbourElements(ws->getNPoints());

// Perform the bulk of the connected component analysis, but don't collapse the elements yet.
LabelIdIntensityMap labelMap; // This will not get used.
PositionToLabelIdMap positionLabelMap; // This will not get used.
calculateDisjointTree(ws, strategy, neighbourElements, labelMap, positionLabelMap);
ElementParentMap parentMap;
calculateDisjointTree(ws, strategy, neighbourElements, labelMap, positionLabelMap, parentMap);

// Create the output workspace from the input workspace
IMDHistoWorkspace_sptr outWS = cloneInputWorkspace(ws);
Expand All @@ -230,9 +249,9 @@ namespace Mantid
for (int i = 0; i < static_cast<int>(neighbourElements.size()); ++i)
{
//std::cout << "Element\t" << i << " Id: \t" << neighbourElements[i].getId() << " This location:\t"<< &neighbourElements[i] << " Root location:\t" << neighbourElements[i].getParent() << " Root Id:\t" << neighbourElements[i].getRoot() << std::endl;
if(!neighbourElements[i].isEmpty())
if (neighbourElements[i] > m_startId)
{
outWS->setSignalAt(i, neighbourElements[i].getRoot());
outWS->setSignalAt(i, parentMap[neighbourElements[i]]);
}
else
{
Expand All @@ -245,34 +264,35 @@ namespace Mantid
}

boost::shared_ptr<Mantid::API::IMDHistoWorkspace> ConnectedComponentLabeling::executeAndIntegrate(
IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy, LabelIdIntensityMap& labelMap,
PositionToLabelIdMap& positionLabelMap) const
IMDHistoWorkspace_sptr ws, BackgroundStrategy * const strategy, LabelIdIntensityMap& labelMap,
PositionToLabelIdMap& positionLabelMap) const
{
VecElements neighbourElements(ws->getNPoints());

// Perform the bulk of the connected component analysis, but don't collapse the elements yet.
calculateDisjointTree(ws, strategy, neighbourElements, labelMap, positionLabelMap);
ElementParentMap parentMap;
calculateDisjointTree(ws, strategy, neighbourElements, labelMap, positionLabelMap, parentMap);

// Create the output workspace from the input workspace
IMDHistoWorkspace_sptr outWS = cloneInputWorkspace(ws);

// Set each pixel to the root of each disjointed element.
for (size_t i = 0; i < neighbourElements.size(); ++i)
{
if(!neighbourElements[i].isEmpty())
if (!neighbourElements[i] > m_startId)
{
const double& signal = ws->getSignalAt(i); // Intensity value at index
double errorSQ = ws->getErrorAt(i);
errorSQ *=errorSQ; // Error squared at index
const size_t& labelId = neighbourElements[i].getRoot();
errorSQ *= errorSQ; // Error squared at index
const size_t& labelId = parentMap[neighbourElements[i]];
// Set the output cluster workspace signal value
outWS->setSignalAt(i, labelId);

SignalErrorSQPair current = labelMap[labelId];
// Sum labels. This is integration!
labelMap[labelId] = SignalErrorSQPair(current.get<0>() + signal, current.get<1>() + errorSQ);
outWS->setSignalAt(i, neighbourElements[i].getRoot());

outWS->setSignalAt(i, labelId);
}
else
{
Expand Down
5 changes: 5 additions & 0 deletions Code/Mantid/Framework/Crystal/src/DisjointElement.cpp
View file
Expand Up @@ -6,6 +6,11 @@ namespace Mantid
namespace Crystal
{

bool operator< (const DisjointElement &a,const DisjointElement &b)
{
return a.getId() < b.getId();
}

/**
* Default constructor. Creates an 'empty' disjoint element.
*/
Expand Down

0 comments on commit 7b2bd12

Please sign in to comment.