HoughBaseAlg.h

Go to the documentation of this file.

// HoughBaseAlg.h
//
// HoughBaseAlg class
//
// Ben Carls (bcarls@fnal.gov)
// Some optimization by Gianluca Petrillo (petrillo@fnal.gov)
//
// Hough transform algorithm
// ============================================================================
//
// This implementation is a pattern recognition algorithm trying to detect
// straight lines in a image.
// In our application, the image is a hitmap on a wire plane, represented by
// wire number as one coordinate and drift time as the other.
// For each hit coordinate, all straight lines through that coordinate are
// recorded on counters, one for each line. If a counter rises to two, it
// mesans that there are two hits who can cast that line, i.e. there are two
// hits on that line.
// A line can be represented by two independent parameters, therefore we need a
// two-dimension parameter space to represent all of them (two degrees of
// freedom).
// Since there are infinite straight lines, each counter represents not just
// one line but a small region of parameter space.
// Passing by one point constraints one of the two parameters, therefore the
// parameter set of a generic line passing by a given point has one degree of
// freedom.
// We follow the custom of choosing as parameters of an arbitrary line passing
// through a point (x, y) the angle (from "x" axis) of the line, and the
// distance of the line from the chosen point. Fixing the point, we can assign
// freedom to the angle (that has a limited range by definition) and then
// the distance will be defined by some functional form d(angle).
// The shape (d vs. a) of the set of parameters of all lines passing by a point
// is a sinusoidal curve. We can define the angle to be between 0 and pi
// (half a period) and the distance potentially covers the whole real axis
// (but actually it will be never larger than the distance of the selected
// point from the origin).
//
// The implementation of the algorithm is based on a "accumulator", the set
// of counters of how many hits are passed by a given line (represented by
// its tow parameters). The accumulator is a two-dimensional container of
// counters, with the first dimension given by the angle and the second by the
// distance.
// In this scenario, all angles are sampled, therefore each angle will have
// at least one count per hit (somewhere at some distance d).
// Each angle will see one entry per hit.
// We choose to have a large number of sampled angles (the current standard
// configuration says 10800), therefore most of the counters will be empty.
// The sinusoidal shape can be steep enough that got example the angle a(1)
// has d(a1)=30 and the next angle (a2) has d(a2)=50. In this case we cover
// also the distances between 31 and 50, (assigning them, as an approximation,
// to a2). In this way we don't leave gaps and make sure that each two
// sinusoidal curves cross in at least one point, or, equivalently, that we
// can always find at least one straight line throw any two points.
// This translates in having very often the counts at a given angle clustered
// around some distances.
//
// We need to discretize the angles and the distances. Tests show that for a
// "natural" plane size of 9600 (TDC counts) x ~3000 (wires), we need to use
// O(10000) angles and to oversample the distance, that typically goes in the
// range [0-10000], by some factor (5 in the default parameters).
// The oversampling factor is just an artifact of the fact that we use integral
// distances but we want to have a resolution better than "1" -- we just
// multiply the distance by e.g. 5 and we get that the distace "1" actually
// means 1/5 = 0.2.
// Note that the algorithm is usually applied to subsets of that plane rather
// than on the full plane, making the typical distance range somehow smaller.
//
// The fastest data structure for the accumulator is a two-dimensional array.
// The proper size of this array would be some gigabyte, that makes this
// approach unfeasible. Since the dimension of angles has a very clear number
// of "bins" (covering always a fixed 0-pi range), but for each given angle
// the counters actually used are a few, a sparse structure (associative
// container, or "map") is used to describe all the counters at a given angle,
// with key the parameter r (discretized).
// Since all the angles always have data, the outer container is a vector
// whose index is the parameter a (also discretized).
// Therefore, for a given line (a, r), we can find how many hits pass though
// it by finding the associative container of the angle a from a vector,
// and in there looking for an entry with d as key: accum[a][d].
//
// Optimization
// ----------------------------------------------------------------------------
//
// Given the constraints and data structure described above, it turns out that
// a lot of time is spent creating new counters. On each hit, and each angle,
// one or more elements of the associative containers must be created.
// In total, millions of counter instances ar used at once.
// The standard C++ implementation of it, std::map, dynamically a new node
// each time a new counter is required, and in the end it frees them one by
// one. Both operations are very time-demanding.
// We used a custom memory allocator (BulkAllocator) that prepares memory
// for chunks of nodes and then returns a preallocated space at each request
// for a new node. The allocator is designed to be fast, giving up features:
// nodes are never really freed, that saves a lot of book-keeping.
// A lot of tricky aspects of it are documented in its own header file.
// Also freeing the memory is fast, since it implies the deletion of a few
// (very large) chunks rather than millions.
// The other aspect taking a lot of time is the lookup of a counter: where is
// counter (a,d)? the location of a is fast (constant time, from a vector).
// The location of d is the next best thing, a binary tree with log(N) access
// time. Unfortunately a balanced tree needs to be rebalanced often, and that
// takes a lot of time. Also, N may be "small" (O(1000)), but there are still
// million insertions and look ups.
// We use here a replacement of the plain map. Since it often happens that,
// to "fill the gaps", sequential counters are allocated and increased,
// we have blocks of couners allocated together (in the current version, 64
// of them). This can save memory in case of crowded spaces: the overhead
// for each node is 40 bytes, whether it is a single counter or a block of
// 64). Look up within the same block becomes constant-time.
// Also, to access sequences of counters, special code is used so that we
// take advantage of the result from the previous look up to perform the next
// one, including also the insertion of a new counter block after an existing
// one.
// Finally, the algorithm acts when it finds a relatively small number of
// aligned hits, afterward removing the hits already clustered. The hit count
// keeps small all the time: we use a signed char as basic data type for
// the counters, allowing a maximum of 127 aligned hits. This saves a lot of
// memory, at the cost of a small slow-down for large (e.g. 64-bit) bus
// architectures. No check is performed for overflow; that can also be
// implemented at a small cost.
//
//
#ifndef HOUGHBASEALG_H
#define HOUGHBASEALG_H

#include <array>
#include <map>
#include <memory> // std::allocator<>
#include <string>
#include <utility> // std::pair<>
#include <vector>

namespace art {
  class Event;
}
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
namespace fhicl {
  class ParameterSet;
}

#include "lardata/Utilities/CountersMap.h"

namespace CLHEP {
  class HepRandomEngine;
}
namespace detinfo {
  class DetectorClocksData;
  class DetectorPropertiesData;
}
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Hit.h"

struct houghCorner {
  double strength = 0;
  double p0 = 0;
  double p1 = 0;
  houghCorner(double strengthTemp = 0, double p0Temp = 0, double p1Temp = 0)
  {
    strength = strengthTemp;
    p0 = p0Temp;
    p1 = p1Temp;
  }

  bool operator<(const houghCorner& houghCornerComp) const
  {
    return (strength < houghCornerComp.strength);
  }
};

// This stores information about merged lines
struct mergedLines {
  double totalQ = 0;
  double pMin0 = 0;
  double pMin1 = 0;
  double pMax0 = 0;
  double pMax1 = 0;
  int clusterNumber = -999999;
  double showerLikeness = 0;
  mergedLines(double totalQTemp = 0,
              double pMin0Temp = 0,
              double pMin1Temp = 0,
              double pMax0Temp = 0,
              double pMax1Temp = 0,
              double clusterNumberTemp = -999999,
              double showerLikenessTemp = 0)
  {
    totalQ = totalQTemp;
    pMin0 = pMin0Temp;
    pMin1 = pMin1Temp;
    pMax0 = pMax0Temp;
    pMax1 = pMax1Temp;
    clusterNumber = clusterNumberTemp;
    showerLikeness = showerLikenessTemp;
  }
};

struct protoTrack {
  int clusterNumber = 999999;
  int planeNumber = 999999;
  int oldClusterNumber = 999999;
  float clusterSlope = 999999;
  float clusterIntercept = 999999;
  float totalQ = -999999;
  float pMin0 = 999999;
  float pMin1 = 999999;
  float pMax0 = -999999;
  float pMax1 = -999999;
  float iMinWire = 999999;
  float iMaxWire = -999999;
  float minWire = 999999;
  float maxWire = -999999;
  float isolation = -999999;
  float showerLikeness = -999999;
  bool merged = false;
  bool showerMerged = false;
  bool mergedLeft = false;
  bool mergedRight = false;
  std::vector<art::Ptr<recob::Hit>> hits;
  protoTrack() {}

  void Init(unsigned int num = 999999,
            unsigned int pnum = 999999,
            float slope = 999999,
            float intercept = 999999,
            float totalQTemp = -999999,
            float Min0 = 999999,
            float Min1 = 999999,
            float Max0 = -999999,
            float Max1 = -999999,
            int iMinWireTemp = 999999,
            int iMaxWireTemp = -999999,
            int minWireTemp = 999999,
            int maxWireTemp = -999999,
            std::vector<art::Ptr<recob::Hit>> hitsTemp = std::vector<art::Ptr<recob::Hit>>())
  {
    clusterNumber = num;
    planeNumber = pnum;
    oldClusterNumber = num;
    clusterSlope = slope;
    clusterIntercept = intercept;
    totalQ = totalQTemp;
    pMin0 = Min0;
    pMin1 = Min1;
    pMax0 = Max0;
    pMax1 = Max1;
    iMinWire = iMinWireTemp;
    iMaxWire = iMaxWireTemp;
    minWire = minWireTemp;
    maxWire = maxWireTemp;
    merged = false;
    showerMerged = false;
    showerLikeness = 0;
    hits.swap(hitsTemp);
  }
};

namespace cluster {

  template <typename KEY,
            typename COUNTER,
            size_t SIZE,
            typename ALLOC = std::allocator<std::pair<KEY, std::array<COUNTER, SIZE>>>,
            unsigned int SUBCOUNTERS = 1>
  class HoughTransformCounters : public lar::CountersMap<KEY, COUNTER, SIZE, ALLOC, SUBCOUNTERS> {
  public:
    using CounterMap_t = HoughTransformCounters<KEY, COUNTER, SIZE, ALLOC, SUBCOUNTERS>;
    using Base_t = lar::CountersMap<KEY, COUNTER, SIZE, ALLOC, SUBCOUNTERS>;

    // import useful types
    using BaseMap_t = typename Base_t::BaseMap_t;
    using Allocator_t = typename Base_t::Allocator_t;
    using Key_t = typename Base_t::Key_t;
    using SubCounter_t = typename Base_t::SubCounter_t;
    using CounterBlock_t = typename Base_t::CounterBlock_t;
    using const_iterator = typename Base_t::const_iterator;

    using PairValue_t = std::pair<const_iterator, SubCounter_t>;

    HoughTransformCounters() : Base_t() {}

    HoughTransformCounters(Allocator_t alloc) : Base_t(alloc) {}

    SubCounter_t set(Key_t key, SubCounter_t value) { return Base_t::set(key, value); }

    SubCounter_t increment(Key_t key) { return Base_t::increment(key); }

    SubCounter_t decrement(Key_t key) { return Base_t::decrement(key); }

    SubCounter_t set(Key_t key_begin, Key_t key_end, SubCounter_t value)
    {
      return unchecked_set_range(key_begin, key_end, value);
    }

    void increment(Key_t key_begin, Key_t key_end);

    PairValue_t increment_and_get_max(Key_t key_begin, Key_t key_end)
    {
      return unchecked_add_range_max(key_begin, key_end, +1);
    }

    PairValue_t increment_and_get_max(Key_t key_begin, Key_t key_end, SubCounter_t current_max)
    {
      return unchecked_add_range_max(key_begin, key_end, +1, current_max);
    }

    void decrement(Key_t key_begin, Key_t key_end);

    PairValue_t get_max(SubCounter_t current_max) const;

    PairValue_t get_max() const;

  protected:
    using CounterKey_t = typename Base_t::CounterKey_t;

  private:
    SubCounter_t unchecked_set_range(Key_t key_begin,
                                     Key_t key_end,
                                     SubCounter_t value,
                                     typename BaseMap_t::iterator start);
    SubCounter_t unchecked_set_range(Key_t key_begin, Key_t key_end, SubCounter_t value);
    PairValue_t unchecked_add_range_max(
      Key_t key_begin,
      Key_t key_end,
      SubCounter_t delta,
      typename BaseMap_t::iterator start,
      SubCounter_t min_max = std::numeric_limits<SubCounter_t>::min());
    PairValue_t unchecked_add_range_max(
      Key_t key_begin,
      Key_t key_end,
      SubCounter_t delta,
      SubCounter_t min_max = std::numeric_limits<SubCounter_t>::min());

  }; // class HoughTransformCounters

  class HoughBaseAlg {
  public:
    struct ChargeInfo_t {
      float integral = 0.0F;
      float integral_stddev = 0.0F;
      float summedADC = 0.0F;
      float summedADC_stddev = 0.0F;

      ChargeInfo_t(float in, float in_stdev, float sum, float sum_stdev)
        : integral(in), integral_stddev(in_stdev), summedADC(sum), summedADC_stddev(sum_stdev)
      {}
    }; // ChargeInfo_t

    explicit HoughBaseAlg(fhicl::ParameterSet const& pset);
    virtual ~HoughBaseAlg() = default;

    size_t FastTransform(const std::vector<art::Ptr<recob::Cluster>>& clusIn,
                         std::vector<recob::Cluster>& ccol,
                         std::vector<art::PtrVector<recob::Hit>>& clusHitsOut,
                         CLHEP::HepRandomEngine& engine,
                         art::Event const& evt,
                         std::string const& label);

    size_t Transform(const detinfo::DetectorClocksData& clockData,
                     detinfo::DetectorPropertiesData const& detProp,
                     std::vector<art::Ptr<recob::Hit>> const& hits,
                     CLHEP::HepRandomEngine& engine,
                     std::vector<unsigned int>* fpointId_to_clusterId,
                     unsigned int clusterId, // The id of the cluster we are examining
                     unsigned int* nClusters,
                     std::vector<protoTrack>* protoTracks);

    // interface to look for lines only on a set of hits,without slope and
    // totalQ arrays
    size_t FastTransform(detinfo::DetectorClocksData const& clockData,
                         detinfo::DetectorPropertiesData const& detProp,
                         std::vector<art::Ptr<recob::Hit>> const& clusIn,
                         std::vector<art::PtrVector<recob::Hit>>& clusHitsOut,
                         CLHEP::HepRandomEngine& engine);

    // interface to look for lines only on a set of hits
    size_t FastTransform(detinfo::DetectorClocksData const& clockData,
                         detinfo::DetectorPropertiesData const& detProp,
                         std::vector<art::Ptr<recob::Hit>> const& clusIn,
                         std::vector<art::PtrVector<recob::Hit>>& clusHitsOut,
                         CLHEP::HepRandomEngine& engine,
                         std::vector<double>& slope,
                         std::vector<ChargeInfo_t>& totalQ);

    size_t Transform(std::vector<art::Ptr<recob::Hit>> const& hits);

    size_t Transform(detinfo::DetectorPropertiesData const& detProp,
                     std::vector<art::Ptr<recob::Hit>> const& hits,
                     double& slope,
                     double& intercept);

    // FIXME: Is this necessary? Document or remove!
    //        2022-04-18 CHG
    friend class HoughTransformClus;

  private:
    void HLSSaveBMPFile(char const*, unsigned char*, int, int);

    int fMaxLines;        
    int fMinHits;         
    int fSaveAccumulator; 
    int fNumAngleCells;   
    float
      fMaxDistance; 
    float fMaxSlope; 
    int
      fRhoZeroOutRange; 
    int
      fThetaZeroOutRange; 
    float fRhoResolutionFactor; 
    int
      fPerCluster; 
    int
      fMissedHits; 
    float
      fMissedHitsDistance; 
    float
      fMissedHitsToLineSize; 
  };

} // namespace

#endif // HOUGHBASEALG_H