DBScanAlg_tool.cc

Go to the documentation of this file.

// Framework Includes
#include "art/Utilities/ToolMacros.h"
#include "cetlib/search_path.h"
#include "cetlib/cpu_timer.h"

#include "larreco/RecoAlg/Cluster3DAlgs/IClusterAlg.h"

// LArSoft includes
#include "larreco/RecoAlg/Cluster3DAlgs/ClusterParamsBuilder.h"
#include "larreco/RecoAlg/Cluster3DAlgs/kdTree.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "lardataobj/RecoBase/Hit.h"

// std includes
#include <string>
#include <functional>
#include <iostream>
#include <memory>

//------------------------------------------------------------------------------------------------------------------------------------------
// implementation follows

namespace lar_cluster3d {

class DBScanAlg : virtual public IClusterAlg
{
public:
    explicit DBScanAlg(fhicl::ParameterSet const &pset);
    
    ~DBScanAlg();
    
    void configure(const fhicl::ParameterSet&) override;
    
    void Cluster3DHits(reco::HitPairList&           hitPairList,
                       reco::ClusterParametersList& clusterParametersList) const override;
    
    void Cluster3DHits(reco::HitPairListPtr&        hitPairList,
                       reco::ClusterParametersList& clusterParametersList) const override;

    float getTimeToExecute(IClusterAlg::TimeValues index) const override {return m_timeVector.at(index);}
    
private:
    
    void expandCluster(const kdTree::KdTreeNode&,
                       kdTree::CandPairList&,
                       reco::ClusterParameters&,
                       size_t) const;
    
    bool                       m_enableMonitoring;      
    size_t                     m_minPairPts;
    mutable std::vector<float> m_timeVector;            
    
    ClusterParamsBuilder       m_clusterBuilder;        // Common cluster builder tool
    kdTree                     m_kdTree;                // For the kdTree
};

DBScanAlg::DBScanAlg(fhicl::ParameterSet const &pset) :
    m_clusterBuilder(pset.get<fhicl::ParameterSet>("ClusterParamsBuilder")),
    m_kdTree(pset.get<fhicl::ParameterSet>("kdTree"))
{
    this->configure(pset);
}

//------------------------------------------------------------------------------------------------------------------------------------------

DBScanAlg::~DBScanAlg()
{
}

//------------------------------------------------------------------------------------------------------------------------------------------

void DBScanAlg::configure(fhicl::ParameterSet const &pset)
{
    m_enableMonitoring         = pset.get<bool>  ("EnableMonitoring",  true  );
    m_minPairPts               = pset.get<size_t>("MinPairPts",        2     );
}
    
void DBScanAlg::Cluster3DHits(reco::HitPairList&           hitPairList,
                              reco::ClusterParametersList& clusterParametersList) const
{
    cet::cpu_timer theClockDBScan;
    
    m_timeVector.resize(NUMTIMEVALUES, 0.);
    
    // DBScan is driven of its "epsilon neighborhood". Computing adjacency within DBScan can be time
    // consuming so the idea is the prebuild the adjaceny map and then run DBScan.
    // We'll employ a kdTree to implement this scheme
    kdTree::KdTreeNodeList kdTreeNodeContainer;
    kdTree::KdTreeNode     topNode = m_kdTree.BuildKdTree(hitPairList, kdTreeNodeContainer);

    if (m_enableMonitoring) m_timeVector.at(BUILDHITTOHITMAP) = m_kdTree.getTimeToExecute();
    
    if (m_enableMonitoring) theClockDBScan.start();
    
    // Ok, here we go!
    // The idea is to loop through all of the input 3D hits and do the clustering
    for(const auto& hit : hitPairList)
    {
        // Check if the hit has already been visited
        if (hit->getStatusBits() & reco::ClusterHit3D::CLUSTERVISITED) continue;
        
        // Mark as visited
        hit->setStatusBit(reco::ClusterHit3D::CLUSTERVISITED);
        
        // Find the neighborhood for this hit
        kdTree::CandPairList candPairList;
        float                bestDistance(std::numeric_limits<float>::max());
        
        m_kdTree.FindNearestNeighbors(hit.get(), topNode, candPairList, bestDistance);
        
        if (candPairList.size() < m_minPairPts)
        {
            hit->setStatusBit(reco::ClusterHit3D::CLUSTERNOISE);
        }
        else
        {
            // "Create" a new cluster and get a reference to it
            clusterParametersList.push_back(reco::ClusterParameters());
            
            reco::ClusterParameters& curCluster = clusterParametersList.back();
            
            hit->setStatusBit(reco::ClusterHit3D::CLUSTERATTACHED);
            curCluster.addHit3D(hit.get());
            
            // expand the cluster
            expandCluster(topNode, candPairList, curCluster, m_minPairPts);
        }
    }

    if (m_enableMonitoring)
    {
        theClockDBScan.stop();
        
        m_timeVector[RUNDBSCAN] = theClockDBScan.accumulated_real_time();
    }
    
    // Initial clustering is done, now trim the list and get output parameters
    cet::cpu_timer theClockBuildClusters;
    
    // Start clocks if requested
    if (m_enableMonitoring) theClockBuildClusters.start();
    
    m_clusterBuilder.BuildClusterInfo(clusterParametersList);
    
    if (m_enableMonitoring)
    {
        theClockBuildClusters.stop();
        
        m_timeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
    }
    
    mf::LogDebug("Cluster3D") << ">>>>> DBScan done, found " << clusterParametersList.size() << " clusters" << std::endl;
    
    return;
}
    
void DBScanAlg::Cluster3DHits(reco::HitPairListPtr&        hitPairList,
                              reco::ClusterParametersList& clusterParametersList) const
{
    cet::cpu_timer theClockDBScan;
    
    m_timeVector.resize(NUMTIMEVALUES, 0.);
    
    // DBScan is driven of its "epsilon neighborhood". Computing adjacency within DBScan can be time
    // consuming so the idea is the prebuild the adjaceny map and then run DBScan.
    // We'll employ a kdTree to implement this scheme
    kdTree::KdTreeNodeList kdTreeNodeContainer;
    kdTree::KdTreeNode     topNode = m_kdTree.BuildKdTree(hitPairList, kdTreeNodeContainer);
    
    if (m_enableMonitoring) m_timeVector.at(BUILDHITTOHITMAP) = m_kdTree.getTimeToExecute();
    
    if (m_enableMonitoring) theClockDBScan.start();
    
    // Ok, here we go!
    // The idea is to loop through all of the input 3D hits and do the clustering
    for(const auto& hit : hitPairList)
    {
        // Check if the hit has already been visited
        if (hit->getStatusBits() & reco::ClusterHit3D::CLUSTERVISITED) continue;
        
        // Mark as visited
        hit->setStatusBit(reco::ClusterHit3D::CLUSTERVISITED);
        
        // Find the neighborhood for this hit
        kdTree::CandPairList candPairList;
        float                bestDistance(std::numeric_limits<float>::max());
        
        m_kdTree.FindNearestNeighbors(hit, topNode, candPairList, bestDistance);
        
        if (candPairList.size() < m_minPairPts)
        {
            hit->setStatusBit(reco::ClusterHit3D::CLUSTERNOISE);
        }
        else
        {
            // "Create" a new cluster and get a reference to it
            clusterParametersList.push_back(reco::ClusterParameters());
            
            reco::ClusterParameters& curCluster = clusterParametersList.back();
            
            hit->setStatusBit(reco::ClusterHit3D::CLUSTERATTACHED);
            curCluster.addHit3D(hit);
            
            // expand the cluster
            expandCluster(topNode, candPairList, curCluster, m_minPairPts);
        }
    }
    
    if (m_enableMonitoring)
    {
        theClockDBScan.stop();
        
        m_timeVector[RUNDBSCAN] = theClockDBScan.accumulated_real_time();
    }
    
    // Initial clustering is done, now trim the list and get output parameters
    cet::cpu_timer theClockBuildClusters;
    
    // Start clocks if requested
    if (m_enableMonitoring) theClockBuildClusters.start();
    
    m_clusterBuilder.BuildClusterInfo(clusterParametersList);
    
    if (m_enableMonitoring)
    {
        theClockBuildClusters.stop();
        
        m_timeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
    }
    
    mf::LogDebug("Cluster3D") << ">>>>> DBScan done, found " << clusterParametersList.size() << " clusters" << std::endl;
    
    return;
}

void DBScanAlg::expandCluster(const kdTree::KdTreeNode& topNode,
                              kdTree::CandPairList&     candPairList,
                              reco::ClusterParameters&  cluster,
                              size_t                    minPts) const
{
    // This is the main inside loop for the DBScan based clustering algorithm
    
    // Loop over added hits until list has been exhausted
    while(!candPairList.empty())
    {
        // Dereference the point so we can see in the debugger...
        const reco::ClusterHit3D* neighborHit = candPairList.front().second;
        
        // Process if we've not been here before
        if (!(neighborHit->getStatusBits() & reco::ClusterHit3D::CLUSTERVISITED))
        {
            // set as visited
            neighborHit->setStatusBit(reco::ClusterHit3D::CLUSTERVISITED);
            
            // get the neighborhood around this point
            kdTree::CandPairList neighborCandPairList;
            float                bestDistance(std::numeric_limits<float>::max());
            
            m_kdTree.FindNearestNeighbors(neighborHit, topNode, neighborCandPairList, bestDistance);

            // If the epsilon neighborhood of this point is large enough then add its points to our list
            if (neighborCandPairList.size() >= minPts)
            {
                std::copy(neighborCandPairList.begin(),neighborCandPairList.end(),std::back_inserter(candPairList));
            }
        }
        
        // If the point is not yet in a cluster then we now add
        if (!(neighborHit->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED))
        {
            neighborHit->setStatusBit(reco::ClusterHit3D::CLUSTERATTACHED);
            cluster.addHit3D(neighborHit);
        }
        
        candPairList.pop_front();
    }
    
    return;
}
    
//------------------------------------------------------------------------------------------------------------------------------------------
    
DEFINE_ART_CLASS_TOOL(DBScanAlg)
} // namespace lar_cluster3d