cluster::DBScanAlg Class Reference

#include <DBScanAlg.h>

Public Member Functions
	DBScanAlg (fhicl::ParameterSet const &pset)
void	InitScan (const detinfo::DetectorClocksData &clockData, const detinfo::DetectorPropertiesData &detProp, const std::vector< art::Ptr< recob::Hit >> &allhits, std::set< uint32_t > badChannels, const std::vector< geo::WireID > &wireids=std::vector< geo::WireID >())
double	getSimilarity (const std::vector< double > v1, const std::vector< double > v2)
std::vector< unsigned int >	findNeighbors (unsigned int pid, double threshold, double threshold2)
void	computeSimilarity ()
void	run_cluster ()
double	getSimilarity2 (const std::vector< double > v1, const std::vector< double > v2)
void	computeSimilarity2 ()
double	getWidthFactor (const std::vector< double > v1, const std::vector< double > v2)
void	computeWidthFactor ()

Public Attributes
std::vector< std::vector< unsigned int > >	fclusters
	collection of something More...
std::vector< std::vector< double > >	fps
	the collection of points we are working on More...
std::vector< unsigned int >	fpointId_to_clusterId
	mapping point_id -> clusterId More...
std::vector< std::vector< double > >	fsim
std::vector< std::vector< double > >	fsim2
std::vector< std::vector< double > >	fsim3
double	fMaxWidth
RTree	fRTree
std::vector< dbsPoint >	fRect

Private Member Functions
void	run_dbscan_cluster ()
void	run_FN_cluster ()
void	run_FN_naive_cluster ()
bool	ExpandCluster (unsigned int point, unsigned int clusterID)
std::set< unsigned int >	RegionQuery (unsigned int point)
std::vector< unsigned int >	RegionQuery_vector (unsigned int point)

Private Attributes
double	fEps
double	fEps2
unsigned int	fMinPts
unsigned int	fClusterMethod
	Which clustering method to use. More...
unsigned int	fDistanceMetric
	Which distance metric to use. More...
std::vector< bool >	fnoise
std::vector< bool >	fvisited
std::vector< double >	fWirePitch
	the pitch of the wires in each plane More...
std::set< uint32_t >	fBadChannels
	set of bad channels in this detector More...
std::vector< uint32_t >	fBadWireSum

Detailed Description

Definition at line 52 of file DBScanAlg.h.

Constructor & Destructor Documentation

cluster::DBScanAlg::DBScanAlg ( fhicl::ParameterSet const &  pset )

explicit

Definition at line 261 of file DBScanAlg.cxx.

{
  fEps = p.get<double>("eps");
  fEps2 = p.get<double>("epstwo");
  fMinPts = p.get<int>("minPts");
  fClusterMethod = p.get<int>("Method");
  fDistanceMetric = p.get<int>("Metric");
}

Member Function Documentation

void cluster::DBScanAlg::computeSimilarity

(

)

Definition at line 487 of file DBScanAlg.cxx.

{
  int size = fps.size();
  fsim.resize(size, std::vector<double>(size));
  for (int i = 0; i < size; i++) {
    for (int j = i + 1; j < size; j++) {
      fsim[j][i] = fsim[i][j] = getSimilarity(fps[i], fps[j]);
    }
  }
}

void cluster::DBScanAlg::computeSimilarity2

(

)

Definition at line 500 of file DBScanAlg.cxx.

{
  int size = fps.size();
  fsim2.resize(size, std::vector<double>(size));
  for (int i = 0; i < size; i++) {
    for (int j = i + 1; j < size; j++) {
      fsim2[j][i] = fsim2[i][j] = getSimilarity2(fps[i], fps[j]);
    }
  }
}

void cluster::DBScanAlg::computeWidthFactor

(

)

Definition at line 513 of file DBScanAlg.cxx.

{
  int size = fps.size();
  fsim3.resize(size, std::vector<double>(size));

  for (int i = 0; i < size; i++) {
    for (int j = i + 1; j < size; j++) {
      fsim3[j][i] = fsim3[i][j] = getWidthFactor(fps[i], fps[j]);
    }
  }
}

bool cluster::DBScanAlg::ExpandCluster ( unsigned int  point, unsigned int  clusterID  )

private

Definition at line 631 of file DBScanAlg.cxx.

{
  /* GetSetOfPoints for point*/
  std::set<unsigned int> seeds = RegionQuery(point);

  // not enough support -> mark as noise
  if (seeds.size() < fMinPts) {
    fpointId_to_clusterId[point] = kNOISE_CLUSTER;
    return false;
  }
  else {
    // Add to the currecnt cluster
    fpointId_to_clusterId[point] = clusterID;
    for (std::set<unsigned int>::iterator itr = seeds.begin(); itr != seeds.end(); itr++) {
      fpointId_to_clusterId[*itr] = clusterID;
    }
    seeds.erase(point);
    while (!seeds.empty()) {
      unsigned int currentP = *(seeds.begin());
      std::set<unsigned int> result = RegionQuery(currentP);

      if (result.size() >= fMinPts) {
        for (std::set<unsigned int>::iterator itr = result.begin(); itr != result.end(); itr++) {
          unsigned int resultP = *itr;
          // not already assigned to a cluster
          if (fpointId_to_clusterId[resultP] == kNO_CLUSTER ||
              fpointId_to_clusterId[resultP] == kNOISE_CLUSTER) {
            if (fpointId_to_clusterId[resultP] == kNO_CLUSTER) { seeds.insert(resultP); }
            fpointId_to_clusterId[resultP] = clusterID;
          } // unclassified or noise
        }   // for
      }     // enough support
      seeds.erase(currentP);
    } // while
    return true;
  }
}

std::vector< unsigned int > cluster::DBScanAlg::findNeighbors	(	unsigned int	pid,
		double	threshold,
		double	threshold2
	)

Definition at line 470 of file DBScanAlg.cxx.

{
  std::vector<unsigned int> ne;

  for (int unsigned j = 0; j < fsim.size(); j++) {
    if ((pid != j) &&
        (((fsim[pid][j]) / (threshold * threshold)) +
         ((fsim2[pid][j]) / (threshold2 * threshold2 * (fsim3[pid][j])))) < 1) { //ellipse
      ne.push_back(j);
    }
  } // end loop over fsim

  return ne;
}

double cluster::DBScanAlg::getSimilarity	(	const std::vector< double >	v1,
		const std::vector< double >	v2
	)

Todo:

this code assumes that all planes have the same wire pitch

Definition at line 365 of file DBScanAlg.cxx.

{

  //for Euclidean distance comment everything out except this-->>>
  // return std::sqrt((v2[1]-v1[1])*(v2[1]-v1[1])+(v2[0]-v1[0])*(v2[0]-v1[0]));
  //------------------------------------------------------------------------
  // return std::abs( v2[0]-v1[0]); //for rectangle
  //----------------------------------------------------------------------
  //Manhattan distance:
  //return std::abs(v1[0]-v2[0])+std::abs(v1[1]-v2[1]);

  /// \todo this code assumes that all planes have the same wire pitch
  double wire_dist = fWirePitch[0];

  unsigned int wire1 =
    (unsigned int)(v1[0] / wire_dist + 0.5); //to make sure to get desired integer
  unsigned int wire2 = (unsigned int)(v2[0] / wire_dist + 0.5);
  int wirestobridge = 0;

  if (wire1 > wire2) {
    unsigned int wire = wire1;
    wire1 = wire2;
    wire2 = wire;
  }

  for (unsigned int i = wire1; i < wire2; i++) {
    if (fBadChannels.find(i) != fBadChannels.end()) wirestobridge++;
  }

  double cmtobridge = wirestobridge * wire_dist;
  //---------------------------------------------------------------------
  return ((std::abs(v2[0] - v1[0]) - cmtobridge) *
          (std::abs(v2[0] - v1[0]) - cmtobridge)); //for ellipse
}

double cluster::DBScanAlg::getSimilarity2	(	const std::vector< double >	v1,
		const std::vector< double >	v2
	)

Todo:

this code assumes all planes have the same wire pitch

Definition at line 402 of file DBScanAlg.cxx.

{

  //-------------------------------------------
  //return std::abs( v2[1]-v1[1]);//for rectangle
  //------------------------------------------

  /// \todo this code assumes all planes have the same wire pitch
  double wire_dist = fWirePitch[0];

  unsigned int wire1 =
    (unsigned int)(v1[0] / wire_dist + 0.5); //to make sure to get desired integer
  unsigned int wire2 = (unsigned int)(v2[0] / wire_dist + 0.5);
  int wirestobridge = 0;

  if (wire1 > wire2) {
    unsigned int wire = wire1;
    wire1 = wire2;
    wire2 = wire;
  }

  for (unsigned int i = wire1; i < wire2; i++) {
    if (fBadChannels.find(i) != fBadChannels.end()) wirestobridge++;
  }

  double cmtobridge = wirestobridge * wire_dist;

  if (std::abs(v2[0] - v1[0]) > 1e-10) {
    cmtobridge *= std::abs((v2[1] - v1[1]) / (v2[0] - v1[0]));
  }
  else
    cmtobridge = 0;

  return ((std::abs(v2[1] - v1[1]) - cmtobridge) *
          (std::abs(v2[1] - v1[1]) - cmtobridge)); //for ellipse
}

double cluster::DBScanAlg::getWidthFactor	(	const std::vector< double >	v1,
		const std::vector< double >	v2
	)

Definition at line 441 of file DBScanAlg.cxx.

{

  //double k=0.13; //this number was determined by looking at flat muon hits' widths.
  //The average width of these hits in cm is 0.505, so 4*2*(w1^2)=2.04
  //where w1=w2=0.505, e^2.044= 7.69. In order not to change the distance
  //in time direction of the ellipse we want to make it equal to 1 for
  //these hits. Thus the k factor is k=1/7.69=0.13//for coeff=4

  //..................................................
  double k = 0.1; //for 4.5 coeff
  double WFactor = (exp(4.6 * ((v1[2] * v1[2]) + (v2[2] * v2[2])))) * k;
  //........................................................
  // Let's try something different:
  if (WFactor > 1) {
    if (WFactor < 6.25)
      return WFactor; //remember that we are increasing the distance in
                      //eps2 as std::sqrt of this number (i.e std::sqrt(6.25))
    else
      return 6.25;
  }
  else
    return 1.0;
}

void cluster::DBScanAlg::InitScan	(	const detinfo::DetectorClocksData &	clockData,
		const detinfo::DetectorPropertiesData &	detProp,
		const std::vector< art::Ptr< recob::Hit >> &	allhits,
		std::set< uint32_t >	badChannels,
		const std::vector< geo::WireID > &	wireids = std::vector<geo::WireID>()
	)

get 2 first wires and find their spacing (wire_dist)

Definition at line 272 of file DBScanAlg.cxx.

{
  if (wireids.size() && wireids.size() != allhits.size()) {
    throw cet::exception("DBScanAlg") << "allhits size = " << allhits.size()
                                      << " wireids size = " << wireids.size() << " do not match\n";
  }
  // clear all the data member vectors for the new set of hits
  fps.clear();
  fpointId_to_clusterId.clear();
  fnoise.clear();
  fvisited.clear();
  fsim.clear();
  fsim2.clear();
  fsim3.clear();
  fclusters.clear();
  fWirePitch.clear();

  fBadChannels = badChannels;
  fBadWireSum.clear();

  // Clear the RTree
  fRTree.Remove(RTree::AcceptAny(), RTree::RemoveLeaf());
  // and the bounds list
  fRect.clear();

  //------------------------------------------------------------------
  // Determine spacing between wires (different for each detector)
  ///get 2 first wires and find their spacing (wire_dist)

  art::ServiceHandle<geo::Geometry const> geom;

  for (size_t p = 0; p < geom->Nplanes(); ++p)
    fWirePitch.push_back(geom->WirePitch(p));

  // Collect the bad wire list into a useful form
  if (fClusterMethod) { // Using the R*-tree
    fBadWireSum.resize(geom->Nchannels());
    unsigned int count = 0;
    for (unsigned int i = 0; i < fBadWireSum.size(); ++i) {
      count += fBadChannels.count(i);
      fBadWireSum[i] = count;
    }
  }

  // Collect the hits in a useful form,
  // and take note of the maximum time width
  fMaxWidth = 0.0;
  for (unsigned int j = 0; j < allhits.size(); ++j) {
    int dims = 3; //our point is defined by 3 elements:wire#,center of the hit, and the hit width
    std::vector<double> p(dims);

    double tickToDist = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
    tickToDist *= 1.e-3 * sampling_rate(clockData); // 1e-3 is conversion of 1/us to 1/ns
    if (!wireids.size())
      p[0] = (allhits[j]->WireID().Wire) * fWirePitch[allhits[j]->WireID().Plane];
    else
      p[0] = (wireids[j].Wire) * fWirePitch[allhits[j]->WireID().Plane];
    p[1] = allhits[j]->PeakTime() * tickToDist;
    p[2] = 2. * allhits[j]->RMS() * tickToDist; //width of a hit in cm

    // check on the maximum width condition
    if (p[2] > fMaxWidth) fMaxWidth = p[2];

    fps.push_back(p);

    if (fClusterMethod) { // Using the R*-tree
      // Convert these same values into dbsPoints to feed into the R*-tree
      dbsPoint pp(p[0], p[1], 0.0, p[2] / 2.0); // note dividing by two
      fRTree.Insert(j, pp.bounds());
      // Keep a parallel list already made up. We could use fps instead, but...
      fRect.push_back(pp);
    }
  }

  fpointId_to_clusterId.resize(fps.size(), kNO_CLUSTER); // Not zero as before!
  fnoise.resize(fps.size(), false);
  fvisited.resize(fps.size(), false);

  if (fClusterMethod) { // Using the R*-tree
    Visitor visitor = fRTree.Query(RTree::AcceptAny(), Visitor());
    mf::LogInfo("DBscan") << "InitScan: hits RTree loaded with " << visitor.count << " items.";
  }
  mf::LogInfo("DBscan") << "InitScan: hits vector size is " << fps.size();

  return;
}

std::set< unsigned int > cluster::DBScanAlg::RegionQuery ( unsigned int  point )

private

Definition at line 593 of file DBScanAlg.cxx.

{
  dbsPoint region(fRect[point]);
  Visitor visitor = fRTree.Query(AcceptFindNeighbors(region.bounds(),
                                                     fEps,
                                                     fEps2,
                                                     fMaxWidth,
                                                     fWirePitch[0], //\todo
                                                     fBadWireSum),  // assumes
                                 Visitor());                        // equal
                                                                    // pitch
  return visitor.sResult;
}

std::vector< unsigned int > cluster::DBScanAlg::RegionQuery_vector ( unsigned int  point )

private

Definition at line 609 of file DBScanAlg.cxx.

{
  dbsPoint region(fRect[point]);
  Visitor visitor = fRTree.Query(AcceptFindNeighbors(region.bounds(),
                                                     fEps,
                                                     fEps2,
                                                     fMaxWidth,
                                                     fWirePitch[0], //\todo
                                                     fBadWireSum),  // assumes
                                 Visitor());                        // equal
                                                                    // pitch
  std::vector<unsigned int>& v = visitor.vResult;
  // find neighbors insures that the called point is not in the
  // returned and this is intended as a drop-in replacement, so insure
  // this condition
  v.erase(std::remove(v.begin(), v.end(), point), v.end());
  return v;
}

void cluster::DBScanAlg::run_cluster

(

)

Definition at line 530 of file DBScanAlg.cxx.

{
  switch (fClusterMethod) {
  case 2: return run_dbscan_cluster();
  case 1: return run_FN_cluster();
  default:
    computeSimilarity();  // watch out for this, they are *slow*
    computeSimilarity2(); // "
    computeWidthFactor(); // "
    return run_FN_naive_cluster();
  }
}

void cluster::DBScanAlg::run_dbscan_cluster ( )

private

Definition at line 549 of file DBScanAlg.cxx.

{
  unsigned int cid = 0;
  // foreach pid
  for (size_t pid = 0; pid < fps.size(); pid++) {
    // not already visited
    if (fpointId_to_clusterId[pid] == kNO_CLUSTER) {
      if (ExpandCluster(pid, cid)) { cid++; }
    } // if (!visited
  }   // for
  //  END DBSCAN

  // Construct clusters, count noise, etc..
  int noise = 0;
  fclusters.resize(cid);
  for (size_t y = 0; y < fpointId_to_clusterId.size(); ++y) {
    if (fpointId_to_clusterId[y] == kNO_CLUSTER) {
      // This shouldn't happen...all points should be clasified by now!
      mf::LogWarning("DBscan") << "Unclassified point!";
    }
    else if (fpointId_to_clusterId[y] == kNOISE_CLUSTER) {
      ++noise;
    }
    else {
      unsigned int c = fpointId_to_clusterId[y];
      if (c >= cid) {
        mf::LogWarning("DBscan") << "Point in cluster " << c << " when only " << cid
                                 << " clusters wer found [0-" << cid - 1 << "]";
      }
      fclusters[c].push_back(y);
    }
  }
  mf::LogInfo("DBscan") << "DWM (R*-tree): Found " << cid << " clusters...";
  for (unsigned int c = 0; c < cid; ++c) {
    mf::LogVerbatim("DBscan") << "\t"
                              << "Cluster " << c << ":\t" << fclusters[c].size();
  }
  mf::LogVerbatim("DBscan") << "\t"
                            << "...and " << noise << " noise points.";
}

void cluster::DBScanAlg::run_FN_cluster ( )

private

Definition at line 676 of file DBScanAlg.cxx.

{

  unsigned int cid = 0;
  // foreach pid
  for (size_t pid = 0; pid < fps.size(); pid++) {
    // not already visited
    if (!fvisited[pid]) {

      fvisited[pid] = true;
      // get the neighbors
      std::vector<unsigned int> ne = RegionQuery_vector(pid);

      // not enough support -> mark as noise
      if (ne.size() < fMinPts) { fnoise[pid] = true; }
      else {
        // Add p to current cluster

        std::vector<unsigned int> c; // a new cluster

        c.push_back(pid); // assign pid to cluster
        fpointId_to_clusterId[pid] = cid;
        // go to neighbors
        for (size_t i = 0; i < ne.size(); ++i) {
          unsigned int nPid = ne[i];

          // not already visited
          if (!fvisited[nPid]) {
            fvisited[nPid] = true;
            // go to neighbors
            std::vector<unsigned int> ne1 = RegionQuery_vector(nPid);
            // enough support
            if (ne1.size() >= fMinPts) {

              // join

              for (size_t i = 0; i < ne1.size(); ++i) {
                // join neighbord
                ne.push_back(ne1[i]);
              }
            }
          }

          // not already assigned to a cluster
          if (fpointId_to_clusterId[nPid] == kNO_CLUSTER) {
            c.push_back(nPid);
            fpointId_to_clusterId[nPid] = cid;
          }
        }

        fclusters.push_back(c);

        cid++;
      }
    } // if (!visited
  }   // for

  int noise = 0;
  // no_hits=fnoise.size();

  for (size_t y = 0; y < fpointId_to_clusterId.size(); ++y) {
    if (fpointId_to_clusterId[y] == kNO_CLUSTER) noise++;
  }
  mf::LogInfo("DBscan") << "FindNeighbors (R*-tree): Found " << cid << " clusters...";
  for (unsigned int c = 0; c < cid; ++c) {
    mf::LogVerbatim("DBscan") << "\t"
                              << "Cluster " << c << ":\t" << fclusters[c].size();
  }
  mf::LogVerbatim("DBscan") << "\t"
                            << "...and " << noise << " noise points.";
}

void cluster::DBScanAlg::run_FN_naive_cluster ( )

private

Definition at line 754 of file DBScanAlg.cxx.

{

  unsigned int cid = 0;
  // foreach pid
  for (size_t pid = 0; pid < fps.size(); ++pid) {
    // not already visited
    if (!fvisited[pid]) {

      fvisited[pid] = true;
      // get the neighbors
      std::vector<unsigned int> ne = findNeighbors(pid, fEps, fEps2);

      // not enough support -> mark as noise
      if (ne.size() < fMinPts) { fnoise[pid] = true; }
      else {
        // Add p to current cluster

        std::vector<unsigned int> c; // a new cluster

        c.push_back(pid); // assign pid to cluster
        fpointId_to_clusterId[pid] = cid;
        // go to neighbors
        for (size_t i = 0; i < ne.size(); ++i) {
          unsigned int nPid = ne[i];

          // not already visited
          if (!fvisited[nPid]) {
            fvisited[nPid] = true;
            // go to neighbors
            std::vector<unsigned int> ne1 = findNeighbors(nPid, fEps, fEps2);
            // enough support
            if (ne1.size() >= fMinPts) {

              // join

              for (unsigned int i = 0; i < ne1.size(); i++) {
                // join neighbord
                ne.push_back(ne1[i]);
              }
            }
          }

          // not already assigned to a cluster
          if (fpointId_to_clusterId[nPid] == kNO_CLUSTER) {
            c.push_back(nPid);
            fpointId_to_clusterId[nPid] = cid;
          }
        }

        fclusters.push_back(c);

        cid++;
      }
    } // if (!visited
  }   // for

  int noise = 0;

  for (size_t y = 0; y < fpointId_to_clusterId.size(); ++y) {
    if (fpointId_to_clusterId[y] == kNO_CLUSTER) ++noise;
  }
  mf::LogInfo("DBscan") << "FindNeighbors (naive): Found " << cid << " clusters...";
  for (unsigned int c = 0; c < cid; ++c) {
    mf::LogVerbatim("DBscan") << "\t"
                              << "Cluster " << c << ":\t" << fclusters[c].size() << " points";
  }
  mf::LogVerbatim("DBscan") << "\t"
                            << "...and " << noise << " noise points.";
}

Member Data Documentation

std::set<uint32_t> cluster::DBScanAlg::fBadChannels

private

set of bad channels in this detector

Definition at line 98 of file DBScanAlg.h.

std::vector<uint32_t> cluster::DBScanAlg::fBadWireSum

private

running total of bad channels. Used for fast intervening dead wire counting ala fBadChannelSum[m]-fBadChannelSum[n].

Definition at line 99 of file DBScanAlg.h.

unsigned int cluster::DBScanAlg::fClusterMethod

private

Which clustering method to use.

Definition at line 91 of file DBScanAlg.h.

std::vector<std::vector<unsigned int> > cluster::DBScanAlg::fclusters

collection of something

Definition at line 71 of file DBScanAlg.h.

unsigned int cluster::DBScanAlg::fDistanceMetric

private

Which distance metric to use.

Definition at line 92 of file DBScanAlg.h.

double cluster::DBScanAlg::fEps

private

Definition at line 86 of file DBScanAlg.h.

double cluster::DBScanAlg::fEps2

private

Definition at line 87 of file DBScanAlg.h.

double cluster::DBScanAlg::fMaxWidth

Definition at line 77 of file DBScanAlg.h.

unsigned int cluster::DBScanAlg::fMinPts

private

Definition at line 89 of file DBScanAlg.h.

std::vector<bool> cluster::DBScanAlg::fnoise

private

Definition at line 95 of file DBScanAlg.h.

std::vector<unsigned int> cluster::DBScanAlg::fpointId_to_clusterId

mapping point_id -> clusterId

Definition at line 73 of file DBScanAlg.h.

std::vector<std::vector<double> > cluster::DBScanAlg::fps

the collection of points we are working on

Definition at line 72 of file DBScanAlg.h.

std::vector<dbsPoint> cluster::DBScanAlg::fRect

Definition at line 80 of file DBScanAlg.h.

RTree cluster::DBScanAlg::fRTree

Definition at line 79 of file DBScanAlg.h.

std::vector<std::vector<double> > cluster::DBScanAlg::fsim

Definition at line 74 of file DBScanAlg.h.

std::vector<std::vector<double> > cluster::DBScanAlg::fsim2

Definition at line 75 of file DBScanAlg.h.

std::vector<std::vector<double> > cluster::DBScanAlg::fsim3

Definition at line 76 of file DBScanAlg.h.

std::vector<bool> cluster::DBScanAlg::fvisited

private

Definition at line 96 of file DBScanAlg.h.

std::vector<double> cluster::DBScanAlg::fWirePitch

private

the pitch of the wires in each plane

Definition at line 97 of file DBScanAlg.h.

---

The documentation for this class was generated from the following files:

• larreco/larreco/RecoAlg/DBScanAlg.h
• larreco/larreco/RecoAlg/DBScanAlg.cxx