DBScanAlg_tool.cc

Go to the documentation of this file.

// Framework Includes
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Utilities/ToolMacros.h"
#include "art/Utilities/make_tool.h"
#include "cetlib/cpu_timer.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// LArSoft includes
#include "larcore/Geometry/WireReadout.h"
#include "larcorealg/Geometry/PlaneGeo.h"
#include "larreco/RecoAlg/Cluster3DAlgs/Cluster3D.h"
#include "larreco/RecoAlg/Cluster3DAlgs/IClusterAlg.h"
#include "larreco/RecoAlg/Cluster3DAlgs/IClusterParamsBuilder.h"
#include "larreco/RecoAlg/Cluster3DAlgs/kdTree.h"

// std includes
#include <memory>

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

namespace lar_cluster3d {

  class DBScanAlg : public IClusterAlg {
  public:
    explicit DBScanAlg(fhicl::ParameterSet const& pset);

    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[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; 

    std::unique_ptr<lar_cluster3d::IClusterParametersBuilder>
      m_clusterBuilder; 
    kdTree m_kdTree;    // For the kdTree
  };

  DBScanAlg::DBScanAlg(fhicl::ParameterSet const& pset)
  {
    m_enableMonitoring = pset.get<bool>("EnableMonitoring", true);
    m_minPairPts = pset.get<size_t>("MinPairPts", 2);

    m_clusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
      pset.get<fhicl::ParameterSet>("ClusterParamsBuilder"));

    // Recover the parameter set for the kdTree
    fhicl::ParameterSet kdTreeParams(pset.get<fhicl::ParameterSet>("kdTree"));

    // Now work out the maximum wire pitch
    auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();

    // Returns the wire pitch per plane assuming they will be the same for all TPCs
    constexpr geo::TPCID tpcid{0, 0};
    std::vector<double> const wirePitchVec{
      wireReadoutGeom.Plane(geo::PlaneID{tpcid, 0}).WirePitch(),
      wireReadoutGeom.Plane(geo::PlaneID{tpcid, 1}).WirePitch(),
      wireReadoutGeom.Plane(geo::PlaneID{tpcid, 2}).WirePitch()};

    float maxBestDist = 1.99 * *std::max_element(wirePitchVec.begin(), wirePitchVec.end());

    kdTreeParams.put_or_replace<float>("RefLeafBestDist", maxBestDist);

    m_kdTree = kdTree(kdTreeParams);
  }

  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[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;
  }

  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[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;
  }

  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();
    }
  }

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

  DEFINE_ART_CLASS_TOOL(DBScanAlg)
} // namespace lar_cluster3d