DBScan.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
*/
#ifndef TRIGTOOLS_TRIG_VSI_DBSCAN
#define TRIGTOOLS_TRIG_VSI_DBSCAN
 
#include "TrigVrtSecInclusive/HashTools.h"
 
#include <vector>
#include <functional>
#include <unordered_set>
#include <unordered_map>
 
namespace TrigVSI {
 
template<typename T>
class Cluster {
   public:
      Cluster(){ m_points.clear(); };
      Cluster( const std::vector<T>&  points ) : m_points( points ){};
      Cluster( std::vector<T>&&       points ) : m_points( std::move(points) ){};
 
      inline T getPoint(size_t ipt) const { return m_points[ipt]; };
      inline size_t nPoints() const { return m_points.size(); };
 
      inline std::vector<T>&  Points() { return m_points; };
      inline const std::vector<T>& getPoints() const { return m_points; };
 
   protected :
      std::vector<T> m_points;
};
 
 
template<typename pointType>
class DBScan {
   public :
      using RegionFunc = std::function<std::vector<pointType>(const pointType&, double)>; 
      DBScan( const std::unordered_set<pointType>&, RegionFunc);
      DBScan();
      ~DBScan(){};
 
      DBScan& operator= (DBScan&&) noexcept = default;  
      DBScan(DBScan&&) noexcept = default;              
 
      inline std::vector<Cluster<pointType>> getClusters() const;
      inline size_t nClusters() const { return m_clusters.size() + m_noisesCluster.size(); };
      inline const Cluster<pointType>& getCluster(size_t) const;
 
      void clusterize(double, size_t);
 
   private :
      RegionFunc  m_regionQuery;
      std::unordered_map<pointType, bool>   m_pointsVisited;
      std::vector<Cluster<pointType>>       m_clusters;
      std::unordered_set<pointType>         m_noises;
      std::vector<Cluster<pointType>>       m_noisesCluster;
};
 
 
//============================================================================
 
template<typename pointType>
DBScan<pointType>::DBScan( const std::unordered_set<pointType>& set, RegionFunc regionQuery):
m_regionQuery(regionQuery)
{
   m_pointsVisited.clear();
   for (const auto& point : set) {
      m_pointsVisited.emplace(point, false);
   }
}
 
 
template<typename pointType>
DBScan<pointType>::DBScan()
{
 
}
 
 
template<typename pointType>
std::vector<Cluster<pointType>> DBScan<pointType>::getClusters() const
{
   auto tmp = m_clusters;
   const auto& tmp_noise = m_noisesCluster;
   tmp.insert( tmp.end(), tmp_noise.begin(), tmp_noise.end() );
   return tmp;
}
 
 
template<typename pointType>
const Cluster<pointType>& DBScan<pointType>::getCluster(size_t icls) const
{
   const Cluster<pointType>& cls_tmp = ( icls < m_clusters.size() )? m_clusters[icls] : m_noisesCluster[icls-m_clusters.size()];
   return cls_tmp;
}
 
 
template<typename pointType>
void DBScan<pointType>::clusterize( double eps, size_t minN )
{
   m_clusters.clear();
   m_noisesCluster.clear();
   for (auto& pair : m_pointsVisited) {
      const pointType& point = pair.first;
 
      if ( pair.second ) continue;  // Continue when the point is visited
      pair.second = true;           // Set the point as visited
      std::vector<pointType> neighbor_points = m_regionQuery(point, eps);
      if ( neighbor_points.size() < minN ) {
         m_noises.emplace(point);
         continue;
      }
 
      m_clusters.emplace_back( std::vector<pointType>{point} );
      Cluster<pointType>& cls = *(--m_clusters.end());
 
      for (size_t i = 0; i < neighbor_points.size(); i++) {
         const auto& q = neighbor_points[i];
 
         // Check if q is labeled as noise
         // If so, remove q from noise lists and add q to cluster
         auto q_itr_noise = m_noises.find(q);
         if ( q_itr_noise != m_noises.end() ) {
            m_noises.erase(q_itr_noise);
            cls.Points().emplace_back(q);
         }
 
         // Check if q had been visited before
         // If visited, q is already labeled as a noise or added to other cluster so skip it
         // If not visited, label q as visited and add it to cluster
         auto q_itr = m_pointsVisited.find(q);
         if ( q_itr == m_pointsVisited.end() ) continue;
         if ( q_itr->second ) continue;
         q_itr->second = true;
         cls.Points().emplace_back(q);
 
         std::vector<pointType> q_neighbor_points = m_regionQuery(q, eps);
         if ( q_neighbor_points.size() >= minN ) neighbor_points.insert( neighbor_points.end(), q_neighbor_points.begin(), q_neighbor_points.end() );
      }
 
   }
 
   // Convert noises to one point clusters
   for (const auto& q : m_noises) {
      m_noisesCluster.emplace_back( std::vector<pointType>{q} );
   }
}
 
 
} // end of namespace TrigVSI
 
#endif