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 455 of file HoughBaseAlg.h.

Constructor & Destructor Documentation

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

explicit

Definition at line 252 of file HoughBaseAlg.cxx.

References util::abs(), cluster::HoughTransform::AddPointReturnMax(), geo::WireID::asPlaneID(), geo::PlaneID::asTPCID(), detinfo::DetectorPropertiesData::DriftVelocity(), e, detinfo::DetectorPropertiesData::Efield(), Get, fhicl::ParameterSet::get(), cluster::HoughTransform::GetAccumSize(), cluster::HoughTransform::GetCell(), cluster::HoughTransform::GetEquation(), hits(), cluster::HoughTransform::Init(), protoTrack::Init(), geo::kInduction, MF_LOG_DEBUG, detinfo::DetectorPropertiesData::NumberTimeSamples(), detinfo::sampling_rate(), cluster::HoughTransform::SetCell(), sqr(), cluster::HoughTransform::SubtractPoint(), detinfo::DetectorPropertiesData::Temperature(), Transform(), lar::dump::vector(), w, x, and y.

{
  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 868 of file HoughBaseAlg.cxx.

References art::PtrVector< T >::end(), Get, hits(), cluster::ClusterParamsImportWrapper< Algo >::ImportHits(), MF_LOG_DEBUG, recob::Hit::PeakTime(), geo::WireID::planeID(), recob::Cluster::Sentry, recob::Hit::SigmaPeakTime(), sw, recob::Hit::View(), geo::WireID::Wire, recob::Hit::WireID(), and xx.

Referenced by cluster::HoughLineFinder::produce().

{
  std::vector<int> skip;

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

  geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
  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, wireReadoutGeom, 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 : wireReadoutGeom.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 998 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 1010 of file HoughBaseAlg.cxx.

References util::abs(), lar::util::StatCollector< T, W >::add(), cluster::HoughTransform::AddPointReturnMax(), geo::WireID::asPlaneID(), geo::PlaneID::asTPCID(), art::PtrVector< T >::back(), art::PtrVector< T >::clear(), detinfo::DetectorPropertiesData::DriftVelocity(), detinfo::DetectorPropertiesData::Efield(), Get, cluster::HoughTransform::GetCell(), cluster::HoughTransform::GetEquation(), cluster::HoughTransform::Init(), geo::kInduction, MF_LOG_DEBUG, art::PtrVector< T >::push_back(), detinfo::DetectorPropertiesData::ReadOutWindowSize(), lar::util::StatCollector< T, W >::RMS(), detinfo::sampling_rate(), cluster::HoughTransform::SetCell(), lar::util::StatCollector< T, W >::Sum(), detinfo::DetectorPropertiesData::Temperature(), lar::dump::vector(), x, and y.

{
  std::vector<int> skip;

  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
  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;
  }

  geo::WireID const& wireid = hit.at(0)->WireID();
  auto const& tpcid = wireid.asPlaneID().asTPCID();

  geo::SigType_t sigt = wireReadoutGeom.SignalType(wireid);

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

  auto const num_planes = wireReadoutGeom.Nplanes(tpcid);
  std::vector<double> wire_pitch(num_planes, 0.);
  for (auto const& plane : wireReadoutGeom.Iterate<geo::PlaneGeo>(tpcid)) {
    auto const p = plane.ID().Plane;
    wire_pitch[p] = plane.WirePitch();
  }

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

  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 = wireReadoutGeom.Plane(wireid.asPlaneID()).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;
    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++;
        }


        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()->ROISummedADC());
        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 *  , unsigned char *  , int  , int    )

private

Definition at line 834 of file HoughBaseAlg.cxx.

References util::size().

{
  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
	)

Referenced by HoughBaseAlg().

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 1355 of file HoughBaseAlg.cxx.

References cluster::HoughTransform::AddPointReturnMax(), Get, cluster::HoughTransform::GetAccumSize(), cluster::HoughTransform::GetCell(), cluster::HoughTransform::GetEquation(), cluster::HoughTransform::GetMax(), hits(), cluster::HoughTransform::Init(), MF_LOG_DEBUG, and detinfo::DetectorPropertiesData::ReadOutWindowSize().

{
  HoughTransform c;

  auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();

  int dx = wireReadoutGeom.Nwires(geo::PlaneID{0, 0, 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);


  return hits.size();
}

Friends And Related Function Documentation

friend class HoughTransformClus

friend

Definition at line 514 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 528 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fMaxLines

private

Max number of lines that can be found.

Definition at line 519 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fMaxSlope

private

Max slope a line can have.

Definition at line 529 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 520 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 539 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 542 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 544 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 523 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 536 of file HoughBaseAlg.h.

float cluster::HoughBaseAlg::fRhoResolutionFactor

private

Factor determining the resolution in rho.

Definition at line 534 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 531 of file HoughBaseAlg.h.

int cluster::HoughBaseAlg::fSaveAccumulator

private

Save bitmap image of accumulator for debugging?

Definition at line 522 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 533 of file HoughBaseAlg.h.

---

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

• larreco/v10_01_07/source/larreco/RecoAlg/HoughBaseAlg.h
• larreco/v10_01_07/source/larreco/RecoAlg/HoughBaseAlg.cxx