SpatialEntityAlignment.hpp

Go to the documentation of this file.

/* 
 * SpatialEntityAlignment.hpp 
 *
 * This file is part of the Chemical Data Processing Toolkit
 *
 * Copyright (C) 2003 Thomas Seidel <thomas.seidel@univie.ac.at>
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; see the file COPYING. If not, write to
 * the Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */
 
#ifndef CDPL_CHEM_SPATIALENTITYALIGNMENT_HPP
#define CDPL_CHEM_SPATIALENTITYALIGNMENT_HPP
 
#include <cstddef>
#include <algorithm>
#include <vector>
#include <set>
#include <unordered_set>
#include <functional>
 
#include <boost/functional/hash.hpp>
 
#include "CDPL/Chem/TopologicalEntityAlignment.hpp"
#include "CDPL/Math/KabschAlgorithm.hpp"
#include "CDPL/Math/Matrix.hpp"
#include "CDPL/Math/Vector.hpp"
#include "CDPL/Util/ObjectStack.hpp"
 
 
namespace CDPL
{
 
    namespace Chem
    {
 
        template <typename T>
        class SpatialEntityAlignment
        {
 
            typedef TopologicalEntityAlignment<T> TopologicalAlignment;
 
          public:
            typedef T EntityType;
 
            typedef typename TopologicalAlignment::ConstEntityIterator ConstEntityIterator;
 
            typedef std::function<bool(const Util::STPairArray&)> TopologicalAlignmentConstraintFunction;
 
            typedef std::function<const Math::Vector3D&(const EntityType&)> Entity3DCoordinatesFunction;
 
            typedef std::function<double(const EntityType&)> EntityWeightFunction;
 
            typedef typename TopologicalAlignment::EntityMatchFunction EntityMatchFunction;
 
            typedef typename TopologicalAlignment::EntityPairMatchFunction EntityPairMatchFunction;
 
            SpatialEntityAlignment();
 
            virtual ~SpatialEntityAlignment() {}
 
            void setMinTopologicalMappingSize(std::size_t min_size);
 
            std::size_t getMinTopologicalMappingSize();
 
            void setEntity3DCoordinatesFunction(const Entity3DCoordinatesFunction& func);
 
            const Entity3DCoordinatesFunction& getEntity3DCoordinatesFunction() const;
 
            void setEntityWeightFunction(const EntityWeightFunction& func);
 
            const EntityWeightFunction& getEntityWeightFunction() const;
 
            void setTopAlignmentConstraintFunction(const TopologicalAlignmentConstraintFunction& func);
 
            const TopologicalAlignmentConstraintFunction& getTopAlignmentConstraintFunction() const;
 
            void setEntityMatchFunction(const EntityMatchFunction& func);
 
            const EntityMatchFunction& getEntityMatchFunction() const;
 
            void setEntityPairMatchFunction(const EntityPairMatchFunction& func);
 
            const EntityPairMatchFunction& getEntityPairMatchFunction() const;
 
            void performExhaustiveSearch(bool exhaustive);
 
            bool exhaustiveSearchPerformed() const;
 
            void addEntity(const EntityType& entity, bool first_set);
 
            void clearEntities(bool first_set);
 
            std::size_t getNumEntities(bool first_set) const;
 
            ConstEntityIterator getEntitiesBegin(bool first_set) const;
 
            ConstEntityIterator getEntitiesEnd(bool first_set) const;
 
            const EntityType& getEntity(std::size_t idx, bool first_set) const;
 
            bool nextAlignment();
 
            void reset();
 
            const Math::Matrix4D& getTransform() const;
 
            const Util::STPairArray& getTopologicalMapping() const;
 
          private:
            void init();
 
            Util::STPairArray* allocTopMapping();
 
            struct TopMappingCmpFunc
            {
 
                bool operator()(const Util::STPairArray* m1, const Util::STPairArray* m2) const
                {
                    return (m1->getSize() > m2->getSize());
                }
            };
 
            struct TopMappingHashFunc
            {
 
                std::size_t operator()(const Util::STPairArray* m) const
                {
                    return boost::hash_value(m->getData());
                }
            };
 
            struct TopMappingEqCmpFunc
            {
 
                bool operator()(const Util::STPairArray* m1, const Util::STPairArray* m2) const
                {
                    return (m1->getData() == m2->getData());
                }
            };
 
            typedef Util::ObjectStack<Util::STPairArray>                                                  TopMappingCache;
            typedef std::vector<Math::Vector3D>                                                           Vector3DArray;
            typedef std::vector<double>                                                                   DoubleArray;
            typedef std::multiset<Util::STPairArray*, TopMappingCmpFunc>                                  TopMappingSet;
            typedef typename TopMappingSet::iterator                                                      TopMappingSetIterator;
            typedef std::unordered_set<const Util::STPairArray*, TopMappingHashFunc, TopMappingEqCmpFunc> TopMappingHashSet;
 
            TopologicalAlignment                   topAlignment;
            TopMappingSet                          topMappings;
            TopMappingSetIterator                  nextTopMappingIter;
            Util::STPairArray*                     currTopMapping;
            TopologicalAlignmentConstraintFunction topAlignConstrFunc;
            Entity3DCoordinatesFunction            coordsFunc;
            EntityWeightFunction                   weightFunc;
            Math::KabschAlgorithm<double>          kabschAlgorithm;
            Vector3DArray                          firstSetCoords;
            DoubleArray                            firstSetWeights;
            Vector3DArray                          secondSetCoords;
            DoubleArray                            secondSetWeights;
            Math::DMatrix                          refPoints;
            Math::DMatrix                          alignedPoints;
            Math::DVector                          almntWeights;
            Math::Matrix4D                         transform;
            std::size_t                            minTopMappingSize;
            bool                                   changes;
            bool                                   exhaustiveMode;
            TopMappingHashSet                      seenTopMappings;
            TopMappingCache                        topMappingCache;
        };
    } // namespace Chem
} // namespace CDPL
 
 
// Implementation
 
template <typename T>
CDPL::Chem::SpatialEntityAlignment<T>::SpatialEntityAlignment():
    minTopMappingSize(3), changes(true), exhaustiveMode(true), topMappingCache(5000)
{
    currTopMapping = allocTopMapping();
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::setMinTopologicalMappingSize(std::size_t min_size)
{
    minTopMappingSize = min_size;
    changes           = true;
}
 
template <typename T>
std::size_t CDPL::Chem::SpatialEntityAlignment<T>::getMinTopologicalMappingSize()
{
    return minTopMappingSize;
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::setEntity3DCoordinatesFunction(const Entity3DCoordinatesFunction& func)
{
    coordsFunc = func;
    changes    = true;
}
 
template <typename T>
const typename CDPL::Chem::SpatialEntityAlignment<T>::Entity3DCoordinatesFunction&
CDPL::Chem::SpatialEntityAlignment<T>::getEntity3DCoordinatesFunction() const
{
    return coordsFunc;
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::setEntityWeightFunction(const EntityWeightFunction& func)
{
    weightFunc = func;
    changes    = true;
}
 
template <typename T>
const typename CDPL::Chem::SpatialEntityAlignment<T>::EntityWeightFunction&
CDPL::Chem::SpatialEntityAlignment<T>::getEntityWeightFunction() const
{
    return weightFunc;
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::setTopAlignmentConstraintFunction(const TopologicalAlignmentConstraintFunction& func)
{
    topAlignConstrFunc = func;
    changes            = true;
}
 
template <typename T>
const typename CDPL::Chem::SpatialEntityAlignment<T>::TopologicalAlignmentConstraintFunction&
CDPL::Chem::SpatialEntityAlignment<T>::getTopAlignmentConstraintFunction() const
{
    return topAlignConstrFunc;
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::setEntityMatchFunction(const EntityMatchFunction& func)
{
    topAlignment.setEntityMatchFunction(func);
}
 
template <typename T>
const typename CDPL::Chem::SpatialEntityAlignment<T>::EntityMatchFunction&
CDPL::Chem::SpatialEntityAlignment<T>::getEntityMatchFunction() const
{
    return topAlignment.getEntityMatchFunction();
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::setEntityPairMatchFunction(const EntityPairMatchFunction& func)
{
    topAlignment.setEntityPairMatchFunction(func);
}
 
template <typename T>
const typename CDPL::Chem::SpatialEntityAlignment<T>::EntityPairMatchFunction&
CDPL::Chem::SpatialEntityAlignment<T>::getEntityPairMatchFunction() const
{
    return topAlignment.getEntityPairMatchFunction();
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::performExhaustiveSearch(bool exhaustive)
{
    exhaustiveMode = exhaustive;
    changes        = true;
}
 
template <typename T>
bool CDPL::Chem::SpatialEntityAlignment<T>::exhaustiveSearchPerformed() const
{
    return exhaustiveMode;
}
 
template <typename T>
std::size_t CDPL::Chem::SpatialEntityAlignment<T>::getNumEntities(bool first_set) const
{
    return topAlignment.getNumEntities(first_set);
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::addEntity(const EntityType& entity, bool first_set)
{
    topAlignment.addEntity(entity, first_set);
    changes = true;
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::clearEntities(bool first_set)
{
    topAlignment.clearEntities(first_set);
    changes = true;
}
 
template <typename T>
typename CDPL::Chem::SpatialEntityAlignment<T>::ConstEntityIterator
CDPL::Chem::SpatialEntityAlignment<T>::getEntitiesBegin(bool first_set) const
{
    return topAlignment.getEntitiesBegin(first_set);
}
 
template <typename T>
typename CDPL::Chem::SpatialEntityAlignment<T>::ConstEntityIterator
CDPL::Chem::SpatialEntityAlignment<T>::getEntitiesEnd(bool first_set) const
{
    return topAlignment.getEntitiesEnd(first_set);
}
 
template <typename T>
const typename CDPL::Chem::SpatialEntityAlignment<T>::EntityType&
CDPL::Chem::SpatialEntityAlignment<T>::getEntity(std::size_t idx, bool first_set) const
{
    return topAlignment.getEntity(idx, first_set);
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::reset()
{
    changes = true;
}
 
template <typename T>
bool CDPL::Chem::SpatialEntityAlignment<T>::nextAlignment()
{
    if (changes)
        init();
 
    if (firstSetCoords.empty() || secondSetCoords.empty())
        return false;
 
    bool        have_weights     = !firstSetWeights.empty();
    std::size_t min_sub_mpg_size = (minTopMappingSize < 3 ? minTopMappingSize : std::size_t(3));
 
    while (nextTopMappingIter != topMappings.end()) {
        currTopMapping = *nextTopMappingIter;
 
        std::size_t num_points = currTopMapping->getSize();
 
        refPoints.resize(3, num_points, false);
        alignedPoints.resize(3, num_points, false);
 
        if (have_weights)
            almntWeights.resize(num_points, 1.0);
 
        std::size_t i = 0;
 
        for (Util::STPairArray::ConstElementIterator it = currTopMapping->getElementsBegin(), end = currTopMapping->getElementsEnd();
             it != end; ++it, i++) {
 
            std::size_t first_idx = it->first;
            std::size_t sec_idx   = it->second;
 
            column(refPoints, i)     = firstSetCoords[first_idx];
            column(alignedPoints, i) = secondSetCoords[sec_idx];
 
            if (have_weights)
                almntWeights(i) = std::max(firstSetWeights[first_idx], secondSetWeights[sec_idx]);
        }
 
        if (exhaustiveMode && (num_points > min_sub_mpg_size)) {
            Util::STPairArray* sub_mpg = 0;
 
            for (std::size_t j = 0; j < num_points; j++) {
                if (!sub_mpg)
                    sub_mpg = allocTopMapping();
 
                sub_mpg->clear();
 
                for (std::size_t k = 0; k < num_points; k++)
                    if (k != j)
                        sub_mpg->addElement(currTopMapping->getElement(k));
 
                if (!seenTopMappings.insert(sub_mpg).second) {
                    continue;
                }
 
                topMappings.insert(sub_mpg);
                sub_mpg = 0;
            }
 
            if (sub_mpg)
                topMappingCache.put();
        }
 
        if (!have_weights) {
            if (!kabschAlgorithm.align(alignedPoints, refPoints)) {
                ++nextTopMappingIter;
                continue;
            }
 
        } else if (!kabschAlgorithm.align(alignedPoints, refPoints, almntWeights)) {
            ++nextTopMappingIter;
            continue;
        }
 
        transform.assign(kabschAlgorithm.getTransform());
        ++nextTopMappingIter;
 
        return true;
    }
 
    return false;
}
 
template <typename T>
const CDPL::Math::Matrix4D& CDPL::Chem::SpatialEntityAlignment<T>::getTransform() const
{
    return transform;
}
 
template <typename T>
const CDPL::Util::STPairArray&
CDPL::Chem::SpatialEntityAlignment<T>::getTopologicalMapping() const
{
    return *currTopMapping;
}
 
template <typename T>
void CDPL::Chem::SpatialEntityAlignment<T>::init()
{
    firstSetCoords.clear();
    secondSetCoords.clear();
 
    firstSetWeights.clear();
    secondSetWeights.clear();
 
    if (coordsFunc) {
        for (ConstEntityIterator it = getEntitiesBegin(true), end = getEntitiesEnd(true); it != end; ++it) {
            const EntityType& entity = *it;
 
            firstSetCoords.push_back(coordsFunc(entity));
 
            if (weightFunc)
                firstSetWeights.push_back(weightFunc(entity));
        }
 
        for (ConstEntityIterator it = getEntitiesBegin(false), end = getEntitiesEnd(false); it != end; ++it) {
            const EntityType& entity = *it;
 
            secondSetCoords.push_back(coordsFunc(entity));
 
            if (weightFunc)
                secondSetWeights.push_back(weightFunc(entity));
        }
    }
 
    topAlignment.reset();
    topMappingCache.putAll();
    seenTopMappings.clear();
    topMappings.clear();
 
    currTopMapping = allocTopMapping();
 
    while (topAlignment.nextAlignment(*currTopMapping)) {
        if (currTopMapping->getSize() < minTopMappingSize)
            continue;
 
        if (topAlignConstrFunc && !topAlignConstrFunc(*currTopMapping))
            continue;
 
        topMappings.insert(currTopMapping);
        currTopMapping = allocTopMapping();
    }
 
    currTopMapping->clear();
 
    nextTopMappingIter = topMappings.begin();
    changes            = false;
}
 
template <typename T>
CDPL::Util::STPairArray*
CDPL::Chem::SpatialEntityAlignment<T>::allocTopMapping()
{
    return topMappingCache.get();
}
 
#endif // CDPL_CHEM_SPATIALENTITYALIGNMENT_HPP