cluster::HoughBaseAlg Class Reference

#include <HoughBaseAlg.h>

Classes
struct	ChargeInfo_t
	Data structure collecting charge information to be filled in cluster. More...

Public Member Functions
	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, unsigned int *nClusters, std::vector< protoTrack > *protoTracks)
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)
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)

Private Member Functions
void	HLSSaveBMPFile (char const *, unsigned char *, int, int)

Private Attributes
int	fMaxLines
	Max number of lines that can be found. More...
int	fMinHits
int	fSaveAccumulator
	Save bitmap image of accumulator for debugging? More...
int	fNumAngleCells
	that fall into the "correct" bin will be small and consistent with noise. More...
float	fMaxDistance
	Max distance that a hit can be from a line to be considered part of that line. More...
float	fMaxSlope
	Max slope a line can have. More...
int	fRhoZeroOutRange
	Range in rho over which to zero out area around previously found lines in the accumulator. More...
int	fThetaZeroOutRange
	Range in theta over which to zero out area around previously found lines in the accumulator. More...
float	fRhoResolutionFactor
	Factor determining the resolution in rho. More...
int	fPerCluster
int	fMissedHits
float	fMissedHitsDistance
	Distance between hits in a hough line before a hit is considered missed. More...
float	fMissedHitsToLineSize
	Ratio of missed hits to line size for a line to be considered a fake. More...

Friends
class	HoughTransformClus

Detailed Description

Definition at line 472 of file HoughBaseAlg.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 261 of file HoughBaseAlg.cxx.

{
  fMaxLines = pset.get<int>("MaxLines");
  fMinHits = pset.get<int>("MinHits");
  fSaveAccumulator = pset.get<int>("SaveAccumulator");
  fNumAngleCells = pset.get<int>("NumAngleCells");
  fMaxDistance = pset.get<float>("MaxDistance");
  fMaxSlope = pset.get<float>("MaxSlope");
  fRhoZeroOutRange = pset.get<int>("RhoZeroOutRange");
  fThetaZeroOutRange = pset.get<int>("ThetaZeroOutRange");
  fRhoResolutionFactor = pset.get<float>("RhoResolutionFactor");
  fPerCluster = pset.get<int>("HitsPerCluster");
  fMissedHits = pset.get<int>("MissedHits");
  fMissedHitsDistance = pset.get<float>("MissedHitsDistance");
  fMissedHitsToLineSize = pset.get<float>("MissedHitsToLineSize");
}

virtual cluster::HoughBaseAlg::~HoughBaseAlg ( )

virtualdefault

Member Function Documentation

size_t cluster::HoughBaseAlg::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
	)

Definition at line 886 of file HoughBaseAlg.cxx.

{
  std::vector<int> skip;

  art::FindManyP<recob::Hit> fmh(clusIn, evt, label);

  geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
  HoughTransform c;

  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
  auto const detProp =
    art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
  util::GeometryUtilities const gser{*geom, clockData, detProp};

  // prepare the algorithm to compute the cluster characteristics;
  // we use the "standard" one here; configuration would happen here,
  // but we are using the default configuration for that algorithm
  ClusterParamsImportWrapper<StandardClusterParamsAlg> ClusterParamAlgo;

  std::vector<art::Ptr<recob::Hit>> hit;

  for (auto view : geom->Views()) {

    MF_LOG_DEBUG("HoughBaseAlg") << "Analyzing view " << view;

    art::PtrVector<recob::Cluster>::const_iterator clusterIter = clusIn.begin();
    int clusterID = 0; //the unique ID of the cluster

    size_t cinctr = 0;
    while (clusterIter != clusIn.end()) {

      MF_LOG_DEBUG("HoughBaseAlg")
        << "Analyzing cinctr=" << cinctr << " which is in view " << (*clusterIter)->View();

      hit.clear();
      if (fPerCluster) {
        if ((*clusterIter)->View() == view) hit = fmh.at(cinctr);
      }
      else {
        while (clusterIter != clusIn.end()) {
          if ((*clusterIter)->View() == view) {

            hit = fmh.at(cinctr);
          } // end if cluster is in correct view
          clusterIter++;
          ++cinctr;
        } //end loop over clusters
      }   //end if not fPerCluster
      if (hit.size() == 0) {
        if (fPerCluster) {
          clusterIter++;
          ++cinctr;
        }
        continue;
      }

      MF_LOG_DEBUG("HoughBaseAlg") << "We have all the hits..." << hit.size();

      std::vector<double> slopevec;
      std::vector<ChargeInfo_t> totalQvec;
      std::vector<art::PtrVector<recob::Hit>> planeClusHitsOut;
      this->FastTransform(clockData, detProp, hit, planeClusHitsOut, engine, slopevec, totalQvec);

      MF_LOG_DEBUG("HoughBaseAlg") << "Made it through FastTransform" << planeClusHitsOut.size();

      for (size_t xx = 0; xx < planeClusHitsOut.size(); ++xx) {
        auto const& hits = planeClusHitsOut.at(xx);
        recob::Hit const& FirstHit = *hits.front();
        recob::Hit const& LastHit = *hits.back();
        const unsigned int sw = FirstHit.WireID().Wire;
        const unsigned int ew = LastHit.WireID().Wire;

        // feed the algorithm with all the cluster hits
        ClusterParamAlgo.ImportHits(gser, hits);

        // create the recob::Cluster directly in the vector;
        // NOTE usually we would use cluster::ClusterCreator to save some typing
        // and some mistakes. In this case, we don't want to pull in the
        // dependency on ClusterFinder, where ClusterCreator currently lives
        ccol.emplace_back(float(sw),                // start_wire
                          0.,                       // sigma_start_wire
                          FirstHit.PeakTime(),      // start_tick
                          FirstHit.SigmaPeakTime(), // sigma_start_tick
                          ClusterParamAlgo.StartCharge(gser).value(),
                          ClusterParamAlgo.StartAngle().value(),
                          ClusterParamAlgo.StartOpeningAngle().value(),
                          float(ew),               // end_wire
                          0.,                      // sigma_end_wire,
                          LastHit.PeakTime(),      // end_tick
                          LastHit.SigmaPeakTime(), // sigma_end_tick
                          ClusterParamAlgo.EndCharge(gser).value(),
                          ClusterParamAlgo.EndAngle().value(),
                          ClusterParamAlgo.EndOpeningAngle().value(),
                          ClusterParamAlgo.Integral().value(),
                          ClusterParamAlgo.IntegralStdDev().value(),
                          ClusterParamAlgo.SummedADC().value(),
                          ClusterParamAlgo.SummedADCStdDev().value(),
                          ClusterParamAlgo.NHits(),
                          ClusterParamAlgo.MultipleHitDensity(),
                          ClusterParamAlgo.Width(gser),
                          clusterID,
                          FirstHit.View(),
                          FirstHit.WireID().planeID(),
                          recob::Cluster::Sentry);

        ++clusterID;
        clusHitsOut.push_back(planeClusHitsOut.at(xx));
      }

      hit.clear();
      if (clusterIter != clusIn.end()) {
        clusterIter++;
        ++cinctr;
      }
      // listofxmax.clear();
      // listofymax.clear();
    } // end loop over clusters

  } // end loop over views

  return ccol.size();
}

size_t cluster::HoughBaseAlg::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
	)

Definition at line 1016 of file HoughBaseAlg.cxx.

{
  std::vector<double> slopevec;
  std::vector<ChargeInfo_t> totalQvec;
  return FastTransform(clockData, detProp, clusIn, clusHitsOut, engine, slopevec, totalQvec);
}

size_t cluster::HoughBaseAlg::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
	)

Todo:

explain where the 0.001 comes from

Todo:

: the collection plane's characteristic hit width's are,

Todo:

: on average, about 5 time samples wider than the induction plane's.

Todo:

: this is hard-coded for now.

Todo:

should it really be wire_pitch[0] in the if statements below, or the pitch for the plane of the hit?

Definition at line 1029 of file HoughBaseAlg.cxx.

{
  std::vector<int> skip;

  geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
  lariov::ChannelStatusProvider const* channelStatus =
    lar::providerFrom<lariov::ChannelStatusService>();

  CLHEP::RandFlat flat(engine);

  std::vector<art::Ptr<recob::Hit>> hit;

  size_t cinctr = 0;
  hit.clear();
  for (std::vector<art::Ptr<recob::Hit>>::const_iterator ii = clusIn.begin(); ii != clusIn.end();
       ii++)
    hit.push_back((*ii)); // this is new

  if (hit.size() == 0) {
    if (fPerCluster) { ++cinctr; }
    return -1;
  }

  unsigned int cs = hit.at(0)->WireID().Cryostat;
  unsigned int t = hit.at(0)->WireID().TPC;
  unsigned int p = hit.at(0)->WireID().Plane;
  geo::WireID const& wireid = hit.at(0)->WireID();

  geo::SigType_t sigt = geom->SignalType(wireid);

  if (hit.size() == 0) {
    if (fPerCluster) { ++cinctr; }
    return -1; // continue;
  }

  std::vector<double> wire_pitch(geom->Nplanes(t, cs), 0.);
  for (size_t p = 0; p < wire_pitch.size(); ++p) {
    wire_pitch[p] = geom->WirePitch(p);
  }

  //factor to make x and y scale the same units
  std::vector<double> xyScale(geom->Nplanes(t, cs), 0.);

  /// \todo explain where the 0.001 comes from
  double driftVelFactor = 0.001 * detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());

  for (size_t p = 0; p < xyScale.size(); ++p)
    xyScale[p] = driftVelFactor * sampling_rate(clockData) / wire_pitch[p];

  int x = 0, y = 0;
  int dx = geom->Cryostat(cs).TPC(t).Plane(p).Nwires(); // number of wires
  const int dy = detProp.ReadOutWindowSize();           // number of time samples.
  skip.clear();
  skip.resize(hit.size());
  std::vector<int> listofxmax;
  std::vector<int> listofymax;
  std::vector<int> hitTemp;        //indecies ofcandidate hits
  std::vector<int> sequenceHolder; //channels of hits in list
  std::vector<int> currentHits;    //working vector of hits
  std::vector<int> lastHits;       //best list of hits
  art::PtrVector<recob::Hit> clusterHits;
  float indcolscaling = 0.; //a parameter to account for the different
  //characteristic hit width of induction and collection plane
  float centerofmassx = 0;
  float centerofmassy = 0;
  float denom = 0;
  float intercept = 0.;
  float slope = 0.;
  //this array keeps track of the hits that have already been associated with a line.
  int xMax = 0;
  int yMax = 0;
  float rho;
  float theta;
  int accDx(0), accDy(0);

  HoughTransform c;
  //Init specifies the size of the two-dimensional accumulator
  //(based on the arguments, number of wires and number of time samples).
  //adds all of the hits (that have not yet been associated with a line) to the accumulator
  c.Init(dx, dy, fRhoResolutionFactor, fNumAngleCells);

  // count is how many points are left to randomly insert
  unsigned int count = hit.size();
  std::vector<unsigned int> accumPoints;
  accumPoints.resize(hit.size());
  int nAccum = 0;
  unsigned int nLinesFound = 0;

  for (; count > 0; count--) {

    // The random hit we are examining
    unsigned int randInd = (unsigned int)(flat.fire() * hit.size());

    MF_LOG_DEBUG("HoughBaseAlg") << "randInd=" << randInd << " and size is " << hit.size();

    // Skip if it's already in a line
    if (skip.at(randInd) == 1) continue;

    // If we have already accumulated the point, skip it
    if (accumPoints.at(randInd)) continue;
    accumPoints.at(randInd) = 1;

    // zeroes out the neighborhood of all previous lines
    for (unsigned int i = 0; i < listofxmax.size(); ++i) {
      int yClearStart = listofymax.at(i) - fRhoZeroOutRange;
      if (yClearStart < 0) yClearStart = 0;

      int yClearEnd = listofymax.at(i) + fRhoZeroOutRange;
      if (yClearEnd >= accDy) yClearEnd = accDy - 1;

      int xClearStart = listofxmax.at(i) - fThetaZeroOutRange;
      if (xClearStart < 0) xClearStart = 0;

      int xClearEnd = listofxmax.at(i) + fThetaZeroOutRange;
      if (xClearEnd >= accDx) xClearEnd = accDx - 1;

      for (y = yClearStart; y <= yClearEnd; ++y) {
        for (x = xClearStart; x <= xClearEnd; ++x) {
          c.SetCell(y, x, 0);
        }
      }
    } // end loop over size of listxmax

    //find the weightiest cell in the accumulator.
    int maxCell = 0;
    unsigned int wireMax = hit.at(randInd)->WireID().Wire;

    // Add the randomly selected point to the accumulator
    std::array<int, 3> max = c.AddPointReturnMax(wireMax, (int)(hit.at(randInd)->PeakTime()));
    maxCell = max.at(0);
    xMax = max.at(1);
    yMax = max.at(2);
    nAccum++;

    // Continue if the biggest maximum for the randomly selected point is smaller than fMinHits
    if (maxCell < fMinHits) continue;

    // Find the center of mass of the 3x3 cell system (with maxCell at the center).
    denom = centerofmassx = centerofmassy = 0;

    if (xMax > 0 && yMax > 0 && xMax + 1 < accDx && yMax + 1 < accDy) {
      for (int i = -1; i < 2; ++i) {
        for (int j = -1; j < 2; ++j) {
          int cell = c.GetCell(yMax + i, xMax + j);
          denom += cell;
          centerofmassx += j * cell;
          centerofmassy += i * cell;
        }
      }
      centerofmassx /= denom;
      centerofmassy /= denom;
    }
    else
      centerofmassx = centerofmassy = 0;

    //fill the list of cells that have already been found
    listofxmax.push_back(xMax);
    listofymax.push_back(yMax);

    // Find the lines equation
    c.GetEquation(yMax + centerofmassy, xMax + centerofmassx, rho, theta);
    MF_LOG_DEBUG("HoughBaseAlg") << "Transform(I) found maximum at (d=" << rho << " a=" << theta
                                 << ")"
                                    " from absolute maximum "
                                 << c.GetCell(yMax, xMax) << " at (d=" << yMax << ", a=" << xMax
                                 << ")";
    slope = -1. / tan(theta);
    intercept = (rho / sin(theta));
    double distance;
    /// \todo: the collection plane's characteristic hit width's are,
    /// \todo: on average, about 5 time samples wider than the induction
    /// plane's. \todo: this is hard-coded for now.
    if (sigt == geo::kInduction)
      indcolscaling = 5.;
    else
      indcolscaling = 0.;

    // note this doesn't work for wrapped wire planes!
    if (!std::isinf(slope) && !std::isnan(slope)) {
      sequenceHolder.clear();
      hitTemp.clear();
      for (size_t i = 0; i < hit.size(); ++i) {
        distance = (TMath::Abs(hit.at(i)->PeakTime() - slope * (double)(hit.at(i)->WireID().Wire) -
                               intercept) /
                    (std::sqrt(pow(xyScale[hit.at(i)->WireID().Plane] * slope, 2) + 1)));

        if (distance < fMaxDistance + hit.at(i)->RMS() + indcolscaling && skip.at(i) != 1) {
          hitTemp.push_back(i);
          sequenceHolder.push_back(hit.at(i)->Channel());
        }

      } // end loop over hits

      if (hitTemp.size() < 2) continue;
      currentHits.clear();
      lastHits.clear();
      int j;
      currentHits.push_back(0);
      for (size_t i = 0; i + 1 < sequenceHolder.size(); ++i) {
        j = 1;
        while ((channelStatus->IsBad(sequenceHolder.at(i) + j)) == true)
          j++;
        if (sequenceHolder.at(i + 1) - sequenceHolder.at(i) <= j + fMissedHits)
          currentHits.push_back(i + 1);
        else if (currentHits.size() > lastHits.size()) {
          lastHits = currentHits;
          currentHits.clear();
        }
        else
          currentHits.clear();
      }

      if (currentHits.size() > lastHits.size()) lastHits = currentHits;

      // Check if lastHits has hits with big gaps in it
      // lastHits[i] is ordered in increasing channel and then increasing peak time,
      // as a consequence, if the line has a negative slope and there are multiple hits in the line for a channel,
      // we have to go back to the first hit (in terms of lastHits[i]) of that channel to find the distance
      // between hits
      //std::cout << "New line" << std::endl;
      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
      int missedHits = 0;
      int lastHitsChannel = 0; //lastHits.at(0);
      int nHitsPerChannel = 1;

      MF_LOG_DEBUG("HoughBaseAlg") << "filling the pCorner arrays around here..."
                                   << "\n but first, lastHits size is " << lastHits.size()
                                   << " and lastHitsChannel=" << lastHitsChannel;

      double pCorner0[2];
      double pCorner1[2];
      unsigned int lastChannel = hit.at(hitTemp.at(lastHits.at(0)))->Channel();

      for (size_t i = 0; i < lastHits.size() - 1; ++i) {
        bool newChannel = false;
        if (slope < 0) {
          if (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel() != lastChannel) {
            newChannel = true;
          }
          if (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel() == lastChannel) nHitsPerChannel++;
        }

        /// \todo should it really be wire_pitch[0] in the if statements below,
        /// or the pitch for the plane of the hit?

        if (slope > 0 || (!newChannel && nHitsPerChannel <= 1)) {

          pCorner0[0] = (hit.at(hitTemp.at(lastHits.at(i)))->Channel()) * wire_pitch[0];
          pCorner0[1] = hit.at(hitTemp.at(lastHits.at(i)))->PeakTime() * tickToDist;
          pCorner1[0] = (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel()) * wire_pitch[0];
          pCorner1[1] = hit.at(hitTemp.at(lastHits.at(i + 1)))->PeakTime() * tickToDist;
          if (std::sqrt(pow(pCorner0[0] - pCorner1[0], 2) + pow(pCorner0[1] - pCorner1[1], 2)) >
              fMissedHitsDistance)
            missedHits++;
        }

        else if (slope < 0 && newChannel && nHitsPerChannel > 1) {

          pCorner0[0] =
            (hit.at(hitTemp.at(lastHits.at(lastHitsChannel)))->Channel()) * wire_pitch[0];
          pCorner0[1] = hit.at(hitTemp.at(lastHits.at(lastHitsChannel)))->PeakTime() * tickToDist;
          pCorner1[0] = (hit.at(hitTemp.at(lastHits.at(i + 1)))->Channel()) * wire_pitch[0];
          pCorner1[1] = hit.at(hitTemp.at(lastHits.at(i + 1)))->PeakTime() * tickToDist;
          if (std::sqrt(pow(pCorner0[0] - pCorner1[0], 2) + pow(pCorner0[1] - pCorner1[1], 2)) >
              fMissedHitsDistance)
            missedHits++;
          lastChannel = hit.at(hitTemp.at(lastHits.at(i)))->Channel();
          lastHitsChannel = i + 1;
          nHitsPerChannel = 0;
        }
      }

      if ((double)missedHits / ((double)lastHits.size() - 1) > fMissedHitsToLineSize) continue;

      clusterHits.clear();
      if (lastHits.size() < 5) continue;

      // reduce rounding errors by using double (RMS is very sensitive to them)
      lar::util::StatCollector<double> integralQ, summedQ;

      for (size_t i = 0; i < lastHits.size(); ++i) {
        clusterHits.push_back(hit.at(hitTemp.at(lastHits.at(i))));
        integralQ.add(clusterHits.back()->Integral());
        summedQ.add(clusterHits.back()->SummedADC());
        skip.at(hitTemp.at(lastHits.at(i))) = 1;
      }
      //protection against very steep uncorrelated hits
      if (std::abs(slope) > fMaxSlope) continue;

      clusHitsOut.push_back(clusterHits);
      slopevec.push_back(slope);
      totalQvec.emplace_back(integralQ.Sum(),
                             integralQ.RMS(), // TODO biased value; should unbias?
                             summedQ.Sum(),
                             summedQ.RMS() // TODO biased value; should unbias?
      );

    } // end if !std::isnan

    nLinesFound++;

    if (nLinesFound > (unsigned int)fMaxLines) break;

  } // end loop over hits

  // saves a bitmap image of the accumulator (useful for debugging),
  // with scaling based on the maximum cell value
  if (fSaveAccumulator) {
    //finds the maximum cell in the accumulator for image scaling
    int cell, pix = 0, maxCell = 0;
    for (y = 0; y < accDy; ++y) {
      for (x = 0; x < accDx; ++x) {
        cell = c.GetCell(y, x);
        if (cell > maxCell) maxCell = cell;
      }
    }

    std::unique_ptr<unsigned char[]> outPix(new unsigned char[accDx * accDy]);
    unsigned int PicIndex = 0;
    for (y = 0; y < accDy; ++y) {
      for (x = 0; x < accDx; ++x) {
        //scales the pixel weights based on the maximum cell value
        if (maxCell > 0) pix = (int)((1500 * c.GetCell(y, x)) / maxCell);
        outPix[PicIndex++] = pix;
      }
    }

    HLSSaveBMPFile("houghaccum.bmp", outPix.get(), accDx, accDy);
  } // end if saving accumulator

  hit.clear();
  lastHits.clear();
  listofxmax.clear();
  listofymax.clear();
  return clusHitsOut.size();
}

void cluster::HoughBaseAlg::HLSSaveBMPFile ( char const *  fileName, unsigned char *  pix, int  dx, int  dy  )

private

Definition at line 851 of file HoughBaseAlg.cxx.

{
  std::ofstream bmpFile(fileName, std::ios::binary);
  bmpFile.write("B", 1);
  bmpFile.write("M", 1);
  int bitsOffset = 54 + 256 * 4;
  int size = bitsOffset + dx * dy; //header plus 256 entry LUT plus pixels
  bmpFile.write((const char*)&size, 4);
  int reserved = 0;
  bmpFile.write((const char*)&reserved, 4);
  bmpFile.write((const char*)&bitsOffset, 4);
  int bmiSize = 40;
  bmpFile.write((const char*)&bmiSize, 4);
  bmpFile.write((const char*)&dx, 4);
  bmpFile.write((const char*)&dy, 4);
  short planes = 1;
  bmpFile.write((const char*)&planes, 2);
  short bitCount = 8;
  bmpFile.write((const char*)&bitCount, 2);
  int i, temp = 0;
  for (i = 0; i < 6; i++)
    bmpFile.write((const char*)&temp, 4); // zero out optional color info
  // write a linear LUT
  char lutEntry[4]; // blue,green,red
  lutEntry[3] = 0;  // reserved part
  for (i = 0; i < 256; i++) {
    lutEntry[0] = lutEntry[1] = lutEntry[2] = i;
    bmpFile.write(lutEntry, sizeof lutEntry);
  }
  // write the actual pixels
  bmpFile.write((const char*)pix, dx * dy);
}

size_t cluster::HoughBaseAlg::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,
		unsigned int *	nClusters,
		std::vector< protoTrack > *	protoTracks
	)

Todo:

provide comment about where the 0.001 comes from

characteristic hit width of induction and collection plane

Todo:

: the collection plane's characteristic hit width's are,

Todo:

: on average, about 5 time samples wider than the induction plane's.

Todo:

: this is hard-coded for now.

Outline of PPHT, J. Matas et. al.

LOOP over hits, picking a random one Enter the point into the accumulator IF it is already in the accumulator or part of a line, skip it Store it in a vector of points that have been chosen

Find max value in accumulator; IF above threshold, create a line Subtract points in line from accumulator

END LOOP over hits, picking a random one

Init specifies the size of the two-dimensional accumulator (based on the arguments, number of wires and number of time samples).

Adds all of the hits to the accumulator

count is how many points are left to randomly insert

The random hit we are examining unsigned int randInd = rand() % hits.size();

If the point isn't in the current fuzzy cluster, skip it

Skip if it's already in a line

If we have already accumulated the point, skip it

zeroes out the neighborhood of all previous lines

end loop over size of listxmax

Find the weightiest cell in the accumulator.

Add the randomly selected point to the accumulator

Find the center of mass of the 3x3 cell system (with maxCell at the center).

fill the list of cells that have already been found

end iterator over hits

Subtract points from the accumulator that have already been used

end if !std::isnan

end loop over hits

Definition at line 280 of file HoughBaseAlg.cxx.

{
  int nClustersTemp = *nClusters;

  geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
  lariov::ChannelStatusProvider const* channelStatus =
    lar::providerFrom<lariov::ChannelStatusService>();

  unsigned int wire = 0;
  unsigned int wireMax = 0;
  unsigned int cs = hits[0]->WireID().Cryostat;
  unsigned int t = hits[0]->WireID().TPC;
  geo::WireID const& wireid = hits[0]->WireID();

  geo::SigType_t sigt = geom->SignalType(wireid);
  std::vector<int> skip;

  std::vector<double> wire_pitch(geom->Nplanes(t, cs), 0.);
  for (size_t p = 0; p < wire_pitch.size(); ++p) {
    wire_pitch[p] = geom->WirePitch(p);
  }

  //factor to make x and y scale the same units
  std::vector<double> xyScale(geom->Nplanes(t, cs), 0.);

  /// \todo provide comment about where the 0.001 comes from
  double driftVelFactor = 0.001 * detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());

  for (size_t p = 0; p < xyScale.size(); ++p)
    xyScale[p] = driftVelFactor * sampling_rate(clockData) / wire_pitch[p];

  float tickToDist = driftVelFactor * sampling_rate(clockData);

  int x, y;
  // there must be a better way to find which plane a cluster comes from
  const int dx = geom->Cryostat(hits[0]->WireID().Cryostat)
                   .TPC(hits[0]->WireID().TPC)
                   .Plane(hits[0]->WireID().Plane)
                   .Nwires();                 // number of wires
  const int dy = detProp.NumberTimeSamples(); // number of time samples.
  skip.clear();
  skip.resize(hits.size());
  std::vector<int> listofxmax;
  std::vector<int> listofymax;
  std::vector<int> hitsTemp;       //indecies ofcandidate hits
  std::vector<int> sequenceHolder; //wire numbers of hits in list
  std::vector<int> channelHolder;  //channels of hits in list
  std::vector<int> currentHits;    //working vector of hits
  std::vector<int> lastHits;       //best list of hits
  std::vector<art::Ptr<recob::Hit>> clusterHits;
  float indcolscaling = 0.; //a parameter to account for the different
                            ////characteristic hit width of induction and collection plane
  /// \todo: the collection plane's characteristic hit width's are,
  /// \todo: on average, about 5 time samples wider than the induction plane's.
  /// \todo: this is hard-coded for now.
  if (sigt == geo::kInduction)
    indcolscaling = 0.;
  else
    indcolscaling = 1.;
  float centerofmassx = 0;
  float centerofmassy = 0;
  float denom = 0;
  float intercept = 0.;
  float slope = 0.;
  //this array keeps track of the hits that have already been associated with a line.
  int xMax = 0;
  int yMax = 0;
  int yClearStart;
  int yClearEnd;
  int xClearStart;
  int xClearEnd;
  int maxCell = 0;
  float rho;
  float theta;
  int accDx(0), accDy(0);
  float pCornerMin[2];
  float pCornerMax[2];

  unsigned int fMaxWire = 0;
  int iMaxWire = 0;
  unsigned int fMinWire = 99999999;
  int iMinWire = -1;

  /// Outline of PPHT, J. Matas et. al.
  /// ---------------------------------------
  ///
  ///LOOP over hits, picking a random one
  ///  Enter the point into the accumulator
  ///  IF it is already in the accumulator or part of a line, skip it
  ///  Store it in a vector of points that have been chosen
  ///
  ///  Find max value in accumulator; IF above threshold, create a line
  ///    Subtract points in line from accumulator
  ///
  ///
  ///END LOOP over hits, picking a random one
  ///

  mf::LogInfo("HoughBaseAlg") << "dealing with " << hits.size() << " hits";

  HoughTransform c;

  ///Init specifies the size of the two-dimensional accumulator
  ///(based on the arguments, number of wires and number of time samples).
  c.Init(dx, dy, fRhoResolutionFactor, fNumAngleCells);
  /// Adds all of the hits to the accumulator

  c.GetAccumSize(accDy, accDx);

  // Define the prototrack object
  protoTrack protoTrackToLoad;

  // The number of lines we've found
  unsigned int nLinesFound = 0;
  std::vector<unsigned int> accumPoints(hits.size(), 0);
  int nAccum = 0;

  /// count is how many points are left to randomly insert
  int count = 0;
  for (auto fpointId_to_clusterIdItr = fpointId_to_clusterId->begin();
       fpointId_to_clusterIdItr != fpointId_to_clusterId->end();
       fpointId_to_clusterIdItr++)
    if (*fpointId_to_clusterIdItr == clusterId) count++;

  unsigned int randInd;

  CLHEP::RandFlat flat(engine);
  TStopwatch w;

  for (; count > 0;) {

    /// The random hit we are examining
    ///unsigned int randInd = rand() % hits.size();
    randInd = (unsigned int)(flat.fire() * hits.size());

    /// If the point isn't in the current fuzzy cluster, skip it
    if (fpointId_to_clusterId->at(randInd) != clusterId) continue;

    --count;

    /// Skip if it's already in a line
    if (skip[randInd] == 1) { continue; }

    /// If we have already accumulated the point, skip it
    if (accumPoints[randInd]) continue;
    accumPoints[randInd] = 1;

    /// zeroes out the neighborhood of all previous lines
    for (auto listofxmaxItr = listofxmax.begin(); listofxmaxItr != listofxmax.end();
         ++listofxmaxItr) {
      yClearStart = listofymax[listofxmaxItr - listofxmax.begin()] - fRhoZeroOutRange;
      if (yClearStart < 0) yClearStart = 0;

      yClearEnd = listofymax[listofxmaxItr - listofxmax.begin()] + fRhoZeroOutRange;
      if (yClearEnd >= accDy) yClearEnd = accDy - 1;

      xClearStart = *listofxmaxItr - fThetaZeroOutRange;
      if (xClearStart < 0) xClearStart = 0;

      xClearEnd = *listofxmaxItr + fThetaZeroOutRange;
      if (xClearEnd >= accDx) xClearEnd = accDx - 1;

      for (y = yClearStart; y <= yClearEnd; ++y) {
        for (x = xClearStart; x <= xClearEnd; ++x) {
          c.SetCell(y, x, 0);
        }
      }
    } /// end loop over size of listxmax

    /// Find the weightiest cell in the accumulator.
    uint32_t channel = hits[randInd]->Channel();
    wireMax = hits[randInd]->WireID().Wire;

    /// Add the randomly selected point to the accumulator
    std::array<int, 3> max = c.AddPointReturnMax(wireMax, (int)(hits[randInd]->PeakTime()));
    maxCell = max[0];
    xMax = max[1];
    yMax = max[2];
    ++nAccum;

    // The threshold calculation using a Binomial distribution
    double binomProbSum = TMath::BinomialI(1 / (double)accDx, nAccum, maxCell);
    if (binomProbSum > 1e-21) continue;

    /// Find the center of mass of the 3x3 cell system (with maxCell at the
    /// center).
    denom = centerofmassx = centerofmassy = 0;

    if (xMax > 0 && yMax > 0 && xMax + 1 < accDx && yMax + 1 < accDy) {
      for (int i = -1; i < 2; ++i) {
        for (int j = -1; j < 2; ++j) {
          denom += c.GetCell(yMax + i, xMax + j);
          centerofmassx += j * c.GetCell(yMax + i, xMax + j);
          centerofmassy += i * c.GetCell(yMax + i, xMax + j);
        }
      }
      centerofmassx /= denom;
      centerofmassy /= denom;
    }
    else
      centerofmassx = centerofmassy = 0;

    ///fill the list of cells that have already been found
    listofxmax.push_back(xMax);
    listofymax.push_back(yMax);

    // Find the lines equation
    c.GetEquation(yMax + centerofmassy, xMax + centerofmassx, rho, theta);
    MF_LOG_DEBUG("HoughBaseAlg") << "Transform(II) found maximum at (d=" << rho << " a=" << theta
                                 << ")"
                                    " from absolute maximum "
                                 << c.GetCell(yMax, xMax) << " at (d=" << yMax << ", a=" << xMax
                                 << ")";
    slope = -1. / tan(theta);
    intercept = (rho / sin(theta));
    float distance;

    if (!std::isinf(slope) && !std::isnan(slope)) {
      channelHolder.clear();
      sequenceHolder.clear();
      fMaxWire = 0;
      iMaxWire = 0;
      fMinWire = 99999999;
      iMinWire = -1;
      hitsTemp.clear();
      for (auto hitsItr = hits.cbegin(); hitsItr != hits.cend(); ++hitsItr) {
        wire = (*hitsItr)->WireID().Wire;
        if (fpointId_to_clusterId->at(hitsItr - hits.begin()) != clusterId) continue;
        channel = (*hitsItr)->Channel();
        distance = (std::abs((*hitsItr)->PeakTime() - slope * (float)((*hitsItr)->WireID().Wire) -
                             intercept) /
                    (std::sqrt(sqr(xyScale[(*hitsItr)->WireID().Plane] * slope) + 1.)));
        if (distance < fMaxDistance + (*hitsItr)->RMS() + indcolscaling &&
            skip[hitsItr - hits.begin()] != 1) {
          hitsTemp.push_back(hitsItr - hits.begin());
          sequenceHolder.push_back(wire);
          channelHolder.push_back(channel);
        }
      } /// end iterator over hits

      if (hitsTemp.size() < 5) continue;
      currentHits.clear();
      lastHits.clear();
      int j;
      currentHits.push_back(0);
      for (auto sequenceHolderItr = sequenceHolder.begin();
           sequenceHolderItr + 1 != sequenceHolder.end();
           ++sequenceHolderItr) {
        j = 1;
        while (channelStatus->IsBad(sequenceHolderItr - sequenceHolder.begin() + j))
          j++;
        if (sequenceHolder[sequenceHolderItr - sequenceHolder.begin() + 1] -
              sequenceHolder[sequenceHolderItr - sequenceHolder.begin()] <=
            j + fMissedHits)
          currentHits.push_back(sequenceHolderItr - sequenceHolder.begin() + 1);
        else if (currentHits.size() > lastHits.size()) {
          lastHits = currentHits;
          currentHits.clear();
        }
        else
          currentHits.clear();
      }
      if (currentHits.size() > lastHits.size()) lastHits = currentHits;

      clusterHits.clear();

      if (lastHits.size() < (unsigned int)fMinHits) continue;

      if (std::abs(slope) > fMaxSlope) { continue; }

      // Add new line to list of lines
      float totalQ = 0;
      std::vector<art::Ptr<recob::Hit>> lineHits;
      nClustersTemp++;
      for (auto lastHitsItr = lastHits.begin(); lastHitsItr != lastHits.end(); ++lastHitsItr) {
        fpointId_to_clusterId->at(hitsTemp[(*lastHitsItr)]) = nClustersTemp - 1;
        totalQ += hits[hitsTemp[(*lastHitsItr)]]->Integral();
        wire = hits[hitsTemp[(*lastHitsItr)]]->WireID().Wire;

        if (!accumPoints[hitsTemp[(*lastHitsItr)]]) count--;

        skip[hitsTemp[(*lastHitsItr)]] = 1;

        lineHits.push_back(hits[hitsTemp[(*lastHitsItr)]]);

        /// Subtract points from the accumulator that have already been used
        if (accumPoints[hitsTemp[(*lastHitsItr)]])
          c.SubtractPoint(wire, (int)(hits[hitsTemp[(*lastHitsItr)]]->PeakTime()));

        if (wire < fMinWire) {
          fMinWire = wire;
          iMinWire = hitsTemp[(*lastHitsItr)];
        }
        if (wire > fMaxWire) {
          fMaxWire = wire;
          iMaxWire = hitsTemp[(*lastHitsItr)];
        }
      }
      int pnum = hits[iMinWire]->WireID().Plane;
      pCornerMin[0] = (hits[iMinWire]->WireID().Wire) * wire_pitch[pnum];
      pCornerMin[1] = hits[iMinWire]->PeakTime() * tickToDist;
      pCornerMax[0] = (hits[iMaxWire]->WireID().Wire) * wire_pitch[pnum];
      pCornerMax[1] = hits[iMaxWire]->PeakTime() * tickToDist;

      protoTrackToLoad.Init(nClustersTemp - 1,
                            pnum,
                            slope,
                            intercept,
                            totalQ,
                            pCornerMin[0],
                            pCornerMin[1],
                            pCornerMax[0],
                            pCornerMax[1],
                            iMinWire,
                            iMaxWire,
                            fMinWire,
                            fMaxWire,
                            lineHits);
      linesFound->push_back(protoTrackToLoad);

    } /// end if !std::isnan

    nLinesFound++;

    if (nLinesFound > (unsigned int)fMaxLines) break;

  } /// end loop over hits

  lastHits.clear();

  listofxmax.clear();
  listofymax.clear();

  // saves a bitmap image of the accumulator (useful for debugging),
  // with scaling based on the maximum cell value
  if (fSaveAccumulator) {
    unsigned char* outPix = new unsigned char[accDx * accDy];
    //finds the maximum cell in the accumulator for image scaling
    int cell, pix = 0, maxCell = 0;
    for (int y = 0; y < accDy; ++y) {
      for (int x = 0; x < accDx; ++x) {
        cell = c.GetCell(y, x);
        if (cell > maxCell) maxCell = cell;
      }
    }
    for (int y = 0; y < accDy; ++y) {
      for (int x = 0; x < accDx; ++x) {
        //scales the pixel weights based on the maximum cell value
        if (maxCell > 0) pix = (int)((1500 * c.GetCell(y, x)) / maxCell);
        outPix[y * accDx + x] = pix;
      }
    }

    HLSSaveBMPFile("houghaccum.bmp", outPix, accDx, accDy);
    delete[] outPix;
  } // end if saving accumulator

  *nClusters = nClustersTemp;

  return 1;
}

size_t cluster::HoughBaseAlg::Transform

(

std::vector< art::Ptr< recob::Hit >> const &

hits

)

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

Todo:

could eventually refine this method to throw out hits that are

Todo:

far from the hough line and refine the fit

Definition at line 1375 of file HoughBaseAlg.cxx.

{
  HoughTransform c;

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

  int dx = geom->Nwires(0);                   // number of wires
  const int dy = detProp.ReadOutWindowSize(); // number of time samples.

  c.Init(dx, dy, fRhoResolutionFactor, fNumAngleCells);

  for (unsigned int i = 0; i < hits.size(); ++i) {
    c.AddPointReturnMax(hits[i]->WireID().Wire, (int)(hits[i]->PeakTime()));
  } // end loop over hits

  //gets the actual two-dimensional size of the accumulator
  int accDx = 0;
  int accDy = 0;
  c.GetAccumSize(accDy, accDx);

  //find the weightiest cell in the accumulator.
  int xMax = 0;
  int yMax = 0;
  c.GetMax(yMax, xMax);

  //find the center of mass of the 3x3 cell system (with maxCell at the center).
  float centerofmassx = 0.;
  float centerofmassy = 0.;
  float denom = 0.;

  if (xMax > 0 && yMax > 0 && xMax + 1 < accDx && yMax + 1 < accDy) {
    for (int i = -1; i < 2; ++i) {
      for (int j = -1; j < 2; ++j) {
        denom += c.GetCell(yMax + i, xMax + j);
        centerofmassx += j * c.GetCell(yMax + i, xMax + j);
        centerofmassy += i * c.GetCell(yMax + i, xMax + j);
      }
    }
    centerofmassx /= denom;
    centerofmassy /= denom;
  }
  else
    centerofmassx = centerofmassy = 0;

  float rho = 0.;
  float theta = 0.;
  c.GetEquation(yMax + centerofmassy, xMax + centerofmassx, rho, theta);
  MF_LOG_DEBUG("HoughBaseAlg") << "Transform(III) found maximum at (d=" << rho << " a=" << theta
                               << ")"
                                  " from absolute maximum "
                               << c.GetCell(yMax, xMax) << " at (d=" << yMax << ", a=" << xMax
                               << ")";
  slope = -1. / tan(theta);
  intercept = rho / sin(theta);

  ///\todo could eventually refine this method to throw out hits that are
  ///\todo far from the hough line and refine the fit

  return hits.size();
}

Friends And Related Function Documentation

friend class HoughTransformClus

friend

Definition at line 531 of file HoughBaseAlg.h.

Member Data Documentation

float cluster::HoughBaseAlg::fMaxDistance

private

Max distance that a hit can be from a line to be considered part of that line.

Definition at line 545 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fMaxLines

private

Max number of lines that can be found.

Definition at line 536 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fMaxSlope

private

Max slope a line can have.

Definition at line 546 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fMinHits

private

Min number of hits in the accumulator to consider (number of hits required to be considered a line).

Definition at line 537 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fMissedHits

private

Number of wires that are allowed to be missed before a line is broken up into segments

Definition at line 556 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fMissedHitsDistance

private

Distance between hits in a hough line before a hit is considered missed.

Definition at line 559 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fMissedHitsToLineSize

private

Ratio of missed hits to line size for a line to be considered a fake.

Definition at line 561 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fNumAngleCells

private

that fall into the "correct" bin will be small and consistent with noise.

Number of angle cells in the accumulator (a measure of the angular resolution of the line finder). If this number is too large than the number of votes

Definition at line 540 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fPerCluster

private

Tells the original Hough algorithm to look at clusters individually, or all hits at once

Definition at line 553 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fRhoResolutionFactor

private

Factor determining the resolution in rho.

Definition at line 551 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fRhoZeroOutRange

private

Range in rho over which to zero out area around previously found lines in the accumulator.

Definition at line 548 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fSaveAccumulator

private

Save bitmap image of accumulator for debugging?

Definition at line 539 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fThetaZeroOutRange

private

Range in theta over which to zero out area around previously found lines in the accumulator.

Definition at line 550 of file HoughBaseAlg.h.

---

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

• HoughBaseAlg.h
• HoughBaseAlg.cxx