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 ()
size_t	FastTransform (const std::vector< art::Ptr< recob::Cluster > > &clusIn, std::vector< recob::Cluster > &ccol, std::vector< art::PtrVector< recob::Hit > > &clusHitsOut, art::Event const &evt, std::string const &label)
size_t	Transform (std::vector< art::Ptr< recob::Hit > > const &hits, std::vector< unsigned int > *fpointId_to_clusterId, unsigned int clusterId, unsigned int *nClusters, std::vector< protoTrack > *protoTracks)
size_t	FastTransform (std::vector< art::Ptr< recob::Hit >> &clusIn, std::vector< art::PtrVector< recob::Hit >> &clusHitsOut)
size_t	FastTransform (std::vector< art::Ptr< recob::Hit >> &clusIn, std::vector< art::PtrVector< recob::Hit >> &clusHitsOut, std::vector< double > &slope, std::vector< ChargeInfo_t > &totalQ)
size_t	Transform (std::vector< art::Ptr< recob::Hit > > const &hits)
size_t	Transform (std::vector< art::Ptr< recob::Hit > > const &hits, double &slope, double &intercept)
virtual void	reconfigure (fhicl::ParameterSet const &pset)

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

Constructor & Destructor Documentation

cluster::HoughBaseAlg::HoughBaseAlg

(

fhicl::ParameterSet const &

pset

)

Definition at line 211 of file HoughBaseAlg.cxx.

References reconfigure().

{
  this->reconfigure(pset);
}

cluster::HoughBaseAlg::~HoughBaseAlg ( )

virtual

Definition at line 217 of file HoughBaseAlg.cxx.

{
}

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,
		art::Event const &	evt,
		std::string const &	label
	)

Definition at line 1051 of file HoughBaseAlg.cxx.

References lar::util::StatCollector< T, W >::add(), cluster::HoughTransform::AddPointReturnMax(), art::PtrVector< T >::back(), art::PtrVector< T >::clear(), geo::GeometryCore::Cryostat(), art::PtrVector< T >::end(), fMaxDistance, fMaxLines, fMaxSlope, fMinHits, fMissedHits, fMissedHitsDistance, fMissedHitsToLineSize, fNumAngleCells, fPerCluster, fRhoResolutionFactor, fRhoZeroOutRange, fSaveAccumulator, fThetaZeroOutRange, cluster::HoughTransform::GetCell(), cluster::HoughTransform::GetEquation(), hits(), HLSSaveBMPFile(), cluster::ClusterParamsImportWrapper< Algo >::ImportHits(), cluster::HoughTransform::Init(), geo::kInduction, LOG_DEBUG, max, geo::GeometryCore::Nplanes(), geo::PlaneGeo::Nwires(), recob::Hit::PeakTime(), geo::TPCGeo::Plane(), geo::WireID::planeID(), art::PtrVector< T >::push_back(), lar::util::StatCollector< T, W >::RMS(), recob::Cluster::Sentry, cluster::HoughTransform::SetCell(), recob::Hit::SigmaPeakTime(), geo::GeometryCore::SignalType(), lar::util::StatCollector< T, W >::Sum(), sw, geo::CryostatGeo::TPC(), lar::dump::vector(), recob::Hit::View(), geo::GeometryCore::Views(), geo::WireID::Wire, recob::Hit::WireID(), geo::GeometryCore::WirePitch(), x, xx, and y.

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>();
  //  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
//  lariov::ChannelStatusProvider const* channelStatus
//    = lar::providerFrom<lariov::ChannelStatusService>();
  HoughTransform c;

  // Get the random number generator
  art::ServiceHandle<art::RandomNumberGenerator> rng;
  CLHEP::HepRandomEngine & engine = rng -> getEngine();
  CLHEP::RandFlat flat(engine);

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

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

      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;
      }
      
      LOG_DEBUG("HoughBaseAlg") << "We have all the hits..." << hit.size();
      
      /*
      //factor to make x and y scale the same units
      double xyScale  = .001*detprop->DriftVelocity(detprop->Efield(),detprop->Temperature());
      xyScale        *= detprop->SamplingRate()/geom->WirePitch(p,t,cs);
      
      int x, y;
      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;
      double indcolscaling = 0.;       //a parameter to account for the different 
      //characteristic hit width of induction and collection plane
      double centerofmassx = 0;
      double centerofmassy = 0;
      double denom = 0; 
      double intercept=0.;
      double slope = 0.;
      //this array keeps track of the hits that have already been associated with a line. 
      int xMax = 0;
      int yMax = 0;
      double rho;
      double theta; 
      int accDx(0), accDy(0);
      
      
      //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());
      
      // 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(unsigned int i = 0; i < listofxmax.size(); ++i){
      int yClearStart = listofymax[i] - fRhoZeroOutRange;
      if (yClearStart < 0) yClearStart = 0;
      
      int yClearEnd = listofymax[i] + fRhoZeroOutRange;
      if (yClearEnd >= accDy) yClearEnd = accDy - 1;
      
      int xClearStart = listofxmax[i] - fThetaZeroOutRange;
      if (xClearStart < 0) xClearStart = 0;
      
      int xClearEnd = listofxmax[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[randInd]->WireID().Wire;
      xMax = 0;
      yMax = 0;
      
      // Add the randomly selected point to the accumulator
      maxCell = c.AddPointReturnMax(wireMax, (int)(hit[randInd]->PeakTime()), &yMax, &xMax, fMinHits);
      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){
      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);
      //c.GetEquation(yMax, xMax, rho, theta);
      slope = -1./tan(theta);    
      intercept = (rho/sin(theta));
      //mf::LogVerbatim("HoughBaseAlg") << std::endl;
      //mf::LogVerbatim("HoughBaseAlg") << "slope: " << slope << " intercept: " << intercept << std::endl; 
      //mf::LogInfo("HoughBaseAlg") << "slope: " << slope << " intercept: " << intercept;
      double distance;
      if(sigt == geo::kInduction)
      indcolscaling = 5.;
      else
      indcolscaling = 0.;
      // What is this?   
      indcolscaling = 0;
      
      if(!std::isinf(slope) && !std::isnan(slope)){
      sequenceHolder.clear();
      hitTemp.clear();
      for(size_t i = 0; i < hit.size(); ++i){
      distance = (TMath::Abs(hit[i]->PeakTime()-slope*(double)(hit[i]->WireID().Wire)-intercept)/(std::sqrt(pow(xyScale*slope,2)+1)));
      
      if(distance < fMaxDistance+hit[i]->RMS()+indcolscaling  && skip[i]!=1){
      hitTemp.push_back(i);
      sequenceHolder.push_back(hit[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[i]+j)) == true) j++;
      if(sequenceHolder[i+1]-sequenceHolder[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
      uint32_t     channel = hit[0]->Channel();
      double wirePitch = geom->WirePitch(geom->View(channel));
      double wire_dist = wirePitch;
      double tickToDist = detprop->DriftVelocity(detprop->Efield(),detprop->Temperature());
      tickToDist *= 1.e-3 * detprop->SamplingRate(); // 1e-3 is conversion of 1/us to 1/ns
      //std::cout << "New line" << std::endl;
      int missedHits=0;
      for(size_t i = 0; i < lastHits.size()-1; ++i) {
      //std::cout << hit[hitTemp[lastHits[i]]]->Channel() << std::endl;
      double pCorner0[2];
      pCorner0[0] = (hit[hitTemp[lastHits[i]]]->Channel())*wire_dist;
      pCorner0[1] = hit[hitTemp[lastHits[i]]]->PeakTime()*tickToDist;
      double pCorner1[2];
      pCorner1[0] = (hit[hitTemp[lastHits[i+1]]]->Channel())*wire_dist;
      pCorner1[1] = hit[hitTemp[lastHits[i+1]]]->PeakTime()*tickToDist;
      //std::cout << std::sqrt( pow(pCorner0[0]-pCorner1[0],2) + pow(pCorner0[1]-pCorner1[1],2)) << std::endl;
      if(std::sqrt( pow(pCorner0[0]-pCorner1[0],2) + pow(pCorner0[1]-pCorner1[1],2)) > fMissedHitsDistance             )
      missedHits++;
      }
      //std::cout << "missedHits " << missedHits << std::endl;
      //std::cout << "lastHits.size() " << lastHits.size() << std::endl;
      //std::cout << "missedHits/lastHits.size() " << (double)missedHits/(double)lastHits.size() << std::endl;
      if((double)missedHits/(double)lastHits.size() > fMissedHitsToLineSize)
      continue;
      
      
      
      
      
      clusterHits.clear();    
      double totalQ = 0.;
      if(lastHits.size() < 5) continue;
      
      
      for(size_t i = 0; i < lastHits.size(); ++i) {
      clusterHits.push_back(hit[hitTemp[lastHits[i]]]);
      totalQ += clusterHits.back()->Integral();
      skip[hitTemp[lastHits[i]]]=1;
      } 
      //protection against very steep uncorrelated hits
      if(std::abs(slope)>fMaxSlope 
      && std::abs((*clusterHits.begin())->Channel()-
      clusterHits[clusterHits.size()-1]->Channel())>=0
      )
      continue;
      
      
      
      unsigned int sw = (*clusterHits.begin())->WireID().Wire;
      unsigned int ew = (*(clusterHits.end()-1))->WireID().Wire;
      
      recob::Cluster cluster(sw, 0.,
      (*clusterHits.begin())->PeakTime(), 0.,
      ew, 0., 
      (clusterHits[clusterHits.size()-1])->PeakTime(), 0.,
      slope, 0., 
      -999., 0., 
      totalQ,
      geom->View((*clusterHits.begin())->Channel()),
      clusterID);             
      
      ++clusterID;
      ccol.push_back(cluster);
      clusHitsOut.push_back(clusterHits);
      
      //Turn off hit sharing. T. Yang 9/14/12
      //              //allow double assignment of first and last hits
      //              for(size_t i = 0; i < lastHits.size(); ++i){ 
      //                if(skip[hitTemp[lastHits[i]]] ==1){
      //                  channel = hit[hitTemp[lastHits[i]]]->Channel();  
      //                  if( channel == sc || channel == ec) skip[i] = 0;
      //                }
      //              }
      
      }// end if !std::isnan
      
      nLinesFound++;
      
      if(nLinesFound>(unsigned int)fMaxLines)
      break;
      
      
      }// end loop over hits*/
      
      std::vector<double> slopevec;
      std::vector<ChargeInfo_t> totalQvec;
      std::vector< art::PtrVector<recob::Hit> >   planeClusHitsOut;
      this->FastTransform(hit,planeClusHitsOut,slopevec,totalQvec );
      
      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;
//      ChargeInfo_t const& charge_info = totalQvec.at(xx); // delegating to algos
        
        // feed the algorithm with all the cluster hits
        ClusterParamAlgo.ImportHits(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().value(),       // start_charge
          ClusterParamAlgo.StartAngle().value(),        // start_angle
          ClusterParamAlgo.StartOpeningAngle().value(), // start_opening
          float(ew),                                    // end_wire
          0.,                                           // sigma_end_wire,
          LastHit.PeakTime(),                           // end_tick
          LastHit.SigmaPeakTime(),                      // sigma_end_tick
          ClusterParamAlgo.EndCharge().value(),         // end_charge
          ClusterParamAlgo.EndAngle().value(),          // end_angle
          ClusterParamAlgo.EndOpeningAngle().value(),   // end_opening
          ClusterParamAlgo.Integral().value(),          // integral
          ClusterParamAlgo.IntegralStdDev().value(),    // integral_stddev
          ClusterParamAlgo.SummedADC().value(),         // summedADC
          ClusterParamAlgo.SummedADCStdDev().value(),   // summedADC_stddev
          ClusterParamAlgo.NHits(),                     // n_hits
          ClusterParamAlgo.MultipleHitDensity(),           // multiple_hit_density
          ClusterParamAlgo.Width(),                     // width
          clusterID,                                    // ID
          FirstHit.View(),                              // view
          FirstHit.WireID().planeID(),                  // plane
          recob::Cluster::Sentry                        // sentry
          );
        
        ++clusterID;
        clusHitsOut.push_back(planeClusHitsOut.at(xx));
      }
      
      
      hit.clear();
      //  lastHits.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	(	std::vector< art::Ptr< recob::Hit >> &	clusIn,
		std::vector< art::PtrVector< recob::Hit >> &	clusHitsOut
	)

size_t cluster::HoughBaseAlg::FastTransform	(	std::vector< art::Ptr< recob::Hit >> &	clusIn,
		std::vector< art::PtrVector< recob::Hit >> &	clusHitsOut,
		std::vector< double > &	slope,
		std::vector< ChargeInfo_t > &	totalQ
	)

void cluster::HoughBaseAlg::HLSSaveBMPFile ( char const *  , unsigned char *  , int  , int    )

protected

Definition at line 1015 of file HoughBaseAlg.cxx.

Referenced by FastTransform(), and Transform().

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

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

virtual

Definition at line 222 of file HoughBaseAlg.cxx.

References fMaxDistance, fMaxLines, fMaxSlope, fMinHits, fMissedHits, fMissedHitsDistance, fMissedHitsToLineSize, fNumAngleCells, fPerCluster, fRhoResolutionFactor, fRhoZeroOutRange, fSaveAccumulator, fThetaZeroOutRange, and fhicl::ParameterSet::get().

Referenced by HoughBaseAlg(), cluster::fuzzyClusterAlg::reconfigure(), and cluster::HoughLineFinder::reconfigure().

{ 
  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"           );
  return;
}

size_t cluster::HoughBaseAlg::Transform	(	std::vector< art::Ptr< recob::Hit > > const &	hits,
		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

Get the random number generator

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

The threshold calculation, see http://www.via.cornell.edu/ece547/projects/Hough/webpage/statistics.html accDx is the number of rho bins,m_rowLength

The threshold calculation using a Poisson distribution instead

Continue if the biggest maximum for the randomly selected point is smaller than fMinHits

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

std::cout << "nClusters: " << *nClusters << std::endl;

Subtract points from the accumulator that have already been used

std::cout << std::endl; std::cout << "pCornerMin[0]: " << pCornerMin[0] << " pCornerMin[1]: " << pCornerMin[1] << std::endl; std::cout << "pCornerMax[0]: " << pCornerMax[0] << " pCornerMax[1]: " << pCornerMax[1] << std::endl;

end if !std::isnan

end loop over hits

Definition at line 248 of file HoughBaseAlg.cxx.

References cluster::HoughTransform::AddPointReturnMax(), geo::GeometryCore::Cryostat(), e, fMaxDistance, fMaxLines, fMaxSlope, fMinHits, fMissedHits, fNumAngleCells, fRhoResolutionFactor, fRhoZeroOutRange, fSaveAccumulator, fThetaZeroOutRange, cluster::HoughTransform::GetAccumSize(), cluster::HoughTransform::GetCell(), cluster::HoughTransform::GetEquation(), hits(), HLSSaveBMPFile(), protoTrack::Init(), cluster::HoughTransform::Init(), geo::kInduction, LOG_DEBUG, max, geo::GeometryCore::Nplanes(), geo::PlaneGeo::Nwires(), geo::TPCGeo::Plane(), cluster::HoughTransform::SetCell(), geo::GeometryCore::SignalType(), sqr(), cluster::HoughTransform::SubtractPoint(), geo::CryostatGeo::TPC(), w, geo::GeometryCore::WirePitch(), x, and y.

Referenced by cluster::fuzzyClusterAlg::run_fuzzy_cluster().

{

  int nClustersTemp = *nClusters;
  
  geo::GeometryCore const* geom = lar::providerFrom<geo::Geometry>();
  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
  lariov::ChannelStatusProvider const* channelStatus
    = lar::providerFrom<lariov::ChannelStatusService>();

  //  uint32_t     channel = hits[0]->Channel();
  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();


  //mf::LogInfo("HoughBaseAlg") << "nClusters is: " << *nClusters;


  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.);
  
  double driftVelFactor = 0.001*detprop->DriftVelocity(detprop->Efield(),detprop->Temperature());
  
  for(size_t p = 0; p < xyScale.size(); ++p)
    xyScale[p] = driftVelFactor * detprop->SamplingRate()/wire_pitch[p];

  float tickToDist = driftVelFactor * detprop->SamplingRate();
  //tickToDist *= 1.e-3 * detprop->SamplingRate(); // 1e-3 is conversion of 1/us to 1/ns


  //mf::LogInfo("HoughBaseAlg") << "xyScale: " << xyScale;
  //mf::LogInfo("HoughBaseAlg") << "tickToDist: " << tickToDist;
  
  int x, y;
  //unsigned int channel, plane, wire, tpc, cstat;
  //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->ReadOutWindowSize(); // number of time samples. 
  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 
  if(sigt == geo::kInduction)
    indcolscaling = 0.;
  else
    indcolscaling = 1.;
  //indcolscaling = 0;
  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];
  //float pCorner0[2];
  //float pCorner1[2];
  //bool newChannel = false;
  //unsigned int lastChannel;

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





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

  c.Init(dx,dy,fRhoResolutionFactor,fNumAngleCells);
  //mf::LogInfo("HoughBaseAlg") << "Beginning PPHT";

  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;

  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;

  art::ServiceHandle<art::RandomNumberGenerator> rng;
  CLHEP::HepRandomEngine & engine = rng -> getEngine();
  CLHEP::RandFlat flat(engine);
  TStopwatch w;
  //float timeTotal = 0;

  for( ; count > 0; ){

    randInd = (unsigned int)(flat.fire()*hits.size());

    //std::cout << count << " " << randInd << std::endl;
    
    if(fpointId_to_clusterId->at(randInd) != clusterId)
      continue;

    --count;

    if(skip[randInd]==1){
      continue;
    }

    if(accumPoints[randInd])
      continue;
    accumPoints[randInd]=1;
   
    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);
        }
      }
    }
  
     uint32_t channel = hits[randInd]->Channel();
    wireMax = hits[randInd]->WireID().Wire;

    //w.Start();
    std::array<int, 3> max = c.AddPointReturnMax(wireMax, (int)(hits[randInd]->PeakTime()));
    maxCell = max[0];
    xMax    = max[1];
    yMax    = max[2];
    //w.Stop();
    //std::cout << "Real Time: " << w.RealTime() << std::endl;
    //timeTotal += w.CpuTime();
    ++nAccum; 

    //mf::LogVerbatim("HoughBaseAlg") << "cout: " << count << " maxCell: " << maxCell << std::endl;
    //mf::LogVerbatim("HoughBaseAlg") << "xMax: " << xMax << " yMax: " << yMax << std::endl;

    //TF1 *threshGaus = new TF1("threshGaus","(1/([0]*std::sqrt(2*TMath::Pi())))*exp(-0.5*pow(((x-[1])/[0]),2))");
    //double sigma = std::sqrt(((double)nAccum/(double)accDx)*(1-1/(double)accDx));
    //double mean = (double)nAccum/(double)accDx;
    //threshGaus->SetParameter(0,sigma);
    //threshGaus->SetParameter(1,mean);
    //mf::LogVerbatim("HoughBaseAlg") << "threshGaus mean: " << mean << " sigma: " << sigma << " accDx: " << accDx << std::endl;
    //mf::LogVerbatim("HoughBaseAlg") << "nAccum: " << nAccum << std::endl;
    //mf::LogVerbatim("HoughBaseAlg") << "threshGaus integral range: " << mean-2*sigma << " to " << maxCell << std::endl;
    //mf::LogVerbatim("HoughBaseAlg") << "threshGaus integral: " << threshGaus->Integral(mean-2*sigma,maxCell) << std::endl;
    //mf::LogVerbatim("HoughBaseAlg") << "threshGaus integral: " << threshGaus->Integral(0,maxCell) << std::endl;


    //double poisProbSum = 0;
    //for(int j = 0; j <= maxCell; j++){
    //double poisProb = TMath::Poisson(j,mean);
    //poisProbSum+=poisProb;
    //mf::LogVerbatim("HoughBaseAlg") << "Poisson: " << poisProb << std::endl;
    //}
    //mf::LogVerbatim("HoughBaseAlg") << "Poisson prob sum: " << poisProbSum << std::endl;
    //mf::LogVerbatim("HoughBaseAlg") << "Probability it is higher: " << 1-poisProbSum << std::endl;

    // Continue if the probability of finding a point, (1-poisProbSum) is the probability of finding a 
    // value of maxCell higher than what it currently is
    //if( (1-poisProbSum) > 1e-13)
    //continue;


    // The threshold calculation using a Binomial distribution instead
    double binomProbSum = TMath::BinomialI(1/(double)accDx,nAccum,maxCell);
    //std::cout << "nAccum: " << nAccum << std::endl;
    //std::cout << "maxCell: " << maxCell << std::endl;
    //std::cout << "BinomialI: " << binomProbSum << std::endl;
    //std::cout << std::endl;
    //mf::LogVerbatim("HoughBaseAlg") << "Probability it is higher: " << 1-binomProbSum << std::endl;
    //Continue if the probability of finding a point, (1-poisProbSum) is the probability of finding a
    //value of maxCell higher than what it currently is
    if( binomProbSum > 1e-21)
      continue;



    //if (maxCell < fMinHits) 
      //continue;


    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;
  
    listofxmax.push_back(xMax);
    listofymax.push_back(yMax);

    // Find the lines equation
    c.GetEquation(yMax+centerofmassy, xMax+centerofmassx, rho, theta);
    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);
        }
      }
   
      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() < 5) continue;
      if(lastHits.size() < (unsigned int)fMinHits) continue;







      //int missedHits=0;
      //int lastHitsChannel = 0;
      //fMaxWire = 0;
      //iMaxWire = 0;
      //fMinWire = 99999999;
      //iMinWire = -1;
      //newChannel = false;
      //lastChannel = hits[hitsTemp[lastHits[0]]]->Channel();
      //for(auto lastHitsItr = lastHits.begin(); lastHitsItr != lastHits.end()-1; ++lastHitsItr) {

        //newChannel = false;
        //if(slope < 0){
          //if(hits[hitsTemp[*lastHitsItr+1]]->Channel() != lastChannel){
            //newChannel = true;
          //}           
        //}
        //if(slope > 0 || !newChannel){

          //pCorner0[0] = (hits[hitsTemp[*lastHitsItr]]->Channel())*wire_dist;
          //pCorner0[1] = hits[hitsTemp[*lastHitsItr]]->PeakTime()*tickToDist;
          //pCorner1[0] = (hits[hitsTemp[*lastHitsItr+1]]->Channel())*wire_dist;
          //pCorner1[1] = hits[hitsTemp[*lastHitsItr+1]]->PeakTime()*tickToDist;
          //if((pCorner0[0]-pCorner1[0])*(pCorner0[0]-pCorner1[0]) + (pCorner0[1]-pCorner1[1])*(pCorner0[1]-pCorner1[1]) > fMissedHitsDistance*fMissedHitsDistance)
            //missedHits++;
        //} else if (slope < 0 && newChannel){
          //pCorner0[0] = (hits[hitsTemp[lastHits[lastHitsChannel]]]->Channel())*wire_dist;
          //pCorner0[1] = hits[hitsTemp[lastHits[lastHitsChannel]]]->PeakTime()*tickToDist;
          //pCorner1[0] = (hits[hitsTemp[*lastHitsItr+1]]->Channel())*wire_dist;
          //pCorner1[1] = hits[hitsTemp[*lastHitsItr+1]]->PeakTime()*tickToDist;
          //if((pCorner0[0]-pCorner1[0])*(pCorner0[0]-pCorner1[0]) + (pCorner0[1]-pCorner1[1])*(pCorner0[1]-pCorner1[1]) > fMissedHitsDistance*fMissedHitsDistance)
            //missedHits++;
          //lastChannel=hits[hitsTemp[*lastHitsItr+1]]->Channel();
          //lastHitsChannel=lastHitsItr-lastHits.begin()+1;
        //}
      //}
      //if((float)missedHits/((float)lastHits.size()-1) > fMissedHitsToLineSize)
        //continue;







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

      //std::cout << std::endl;
      //std::cout  << "Found line!" << std::endl
                                       //<< "Slope: " << slope << std::endl
                                       //<< "Intercept: " << intercept << std::endl
                                       //<< "Number of hits: " << lastHits.size() << std::endl
                                       //<< "Wire: " << fMinWire << " Peak time: " 
                                       //<< hits[iMinWire]->PeakTime() << std::endl
                                       //<< "Wire: " << fMaxWire << " Peak time: " 
                                       //<< hits[iMaxWire]->PeakTime() << std::endl;


      // 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;
        //clusterHits.push_back(hits[hitsTemp[(*lastHitsItr)]]);
        //totalQ += clusterHits.back()->Integral();
        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)]]);


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


    nLinesFound++;

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

  }

  //std::cout << "Average cpu time: " << timeTotal/nAccum << std::endl;
  //std::cout << "Total cpu time: " << timeTotal << std::endl;





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

References cluster::HoughTransform::AddPointReturnMax(), fNumAngleCells, fRhoResolutionFactor, cluster::HoughTransform::GetAccumSize(), cluster::HoughTransform::GetCell(), cluster::HoughTransform::GetEquation(), cluster::HoughTransform::GetMax(), hits(), cluster::HoughTransform::Init(), LOG_DEBUG, geo::GeometryCore::Nwires(), and detinfo::DetectorProperties::ReadOutWindowSize().

{
  HoughTransform c;

  art::ServiceHandle<geo::Geometry> geom;
  const detinfo::DetectorProperties* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
  
  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);
  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 635 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 621 of file HoughBaseAlg.h.

Referenced by FastTransform(), reconfigure(), and Transform().

int cluster::HoughBaseAlg::fMaxLines

private

Max number of lines that can be found.

Definition at line 613 of file HoughBaseAlg.h.

Referenced by FastTransform(), reconfigure(), and Transform().

float cluster::HoughBaseAlg::fMaxSlope

private

Max slope a line can have.

Definition at line 622 of file HoughBaseAlg.h.

Referenced by FastTransform(), reconfigure(), and Transform().

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

Referenced by FastTransform(), reconfigure(), and Transform().

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

Referenced by FastTransform(), reconfigure(), and Transform().

float cluster::HoughBaseAlg::fMissedHitsDistance

private

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

Definition at line 630 of file HoughBaseAlg.h.

Referenced by FastTransform(), and reconfigure().

float cluster::HoughBaseAlg::fMissedHitsToLineSize

private

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

Definition at line 631 of file HoughBaseAlg.h.

Referenced by FastTransform(), and reconfigure().

int cluster::HoughBaseAlg::fNumAngleCells

private

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 that fall into the "correct" bin will be small and consistent with noise.

Definition at line 617 of file HoughBaseAlg.h.

Referenced by FastTransform(), reconfigure(), and Transform().

int cluster::HoughBaseAlg::fPerCluster

private

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

Definition at line 626 of file HoughBaseAlg.h.

Referenced by FastTransform(), and reconfigure().

float cluster::HoughBaseAlg::fRhoResolutionFactor

private

Factor determining the resolution in rho.

Definition at line 625 of file HoughBaseAlg.h.

Referenced by FastTransform(), reconfigure(), and Transform().

int cluster::HoughBaseAlg::fRhoZeroOutRange

private

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

Definition at line 623 of file HoughBaseAlg.h.

Referenced by FastTransform(), reconfigure(), and Transform().

int cluster::HoughBaseAlg::fSaveAccumulator

private

Save bitmap image of accumulator for debugging?

Definition at line 616 of file HoughBaseAlg.h.

Referenced by FastTransform(), reconfigure(), and Transform().

int cluster::HoughBaseAlg::fThetaZeroOutRange

private

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

Definition at line 624 of file HoughBaseAlg.h.

Referenced by FastTransform(), reconfigure(), and Transform().

---

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

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