Public Member Functions | Public Attributes | Private Member Functions | List of all members
hough::MultiHough2P Class Reference

Hough transform based on 2-point combinations. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-03-07/HoughTransform/MultiHough2P.h"

Public Member Functions

 MultiHough2P (geo::View_t v, double rhosz, double thetasz, int rhosm, int thetasm, double rsqr, double pco, int loops, double rmdist)
 Construct the 2-point Hough transform. More...
 
 ~MultiHough2P ()
 
TH2F * MapToTH2F () const
 
void Map (const rb::HitList &h)
 
void MultiMap (rb::HitList h)
 

Public Attributes

geo::View_t fView
 Which detector view? More...
 
double fRhoSz
 Size of rho bins (cm) More...
 
double fThetaSz
 Size of theta bins (degrees) More...
 
int fXwinSz
 Smoothing size (bins) in x. More...
 
int fYwinSz
 Smoothing size (bins) in y. More...
 
double fSigma
 Assumed location error on points (cm) More...
 
double fRsqr
 Distance^2 cut on points (cm^2) More...
 
double fPco
 Number of sigma abouve average peak height to use as cutoff. More...
 
int fLoops
 Max number of Multi-Hough loops. More...
 
double fRmDst
 Distance cut for removing points from a Hough line. More...
 
LiteTH2fHspaceW
 Hough transform. More...
 
rb::HoughResult fH
 Results of Hough transformation;. More...
 

Private Member Functions

void BuildMap (double xlo, double xhi, double ylo, double yhi)
 
void RhoTheta (double x1, double y1, double x2, double y2, double *rho, double *theta, double *sigmarho, double *sigmatheta)
 
void SmoothMap ()
 
void RefinePeak (double &rho, double &theta, int xp, int yp)
 
void ReweightHits (double rho, double theta, rb::HitList &h)
 

Detailed Description

Hough transform based on 2-point combinations.

Definition at line 21 of file MultiHough2P.h.

Constructor & Destructor Documentation

hough::MultiHough2P::MultiHough2P ( geo::View_t  v,
double  rhosz,
double  thetasz,
int  rhosm,
int  thetasm,
double  rsqr,
double  pco,
int  loops,
double  rmdist 
)

Construct the 2-point Hough transform.

The transform proceeds by considering all paris of points that occur within a cut off distance (rsqr) of each other. For each pair an ellipse is filled in Hough space weighted by the quality of the line fit to the two spatial points. The algorithm is largely derived from

"Real-time line detection through an improved Hough transform voting scheme", L.A.F. Fernandes, M.M. Oliveira / Pattern Recognition 41 (2008) 299-314

Parameters
v- Which detector view? geo::kX or geo::kY
rhosz- Size of rho bins (cm)
thetasz- Size of theta bins (degrees)
rhosm- Number of bins to use in smoothing rho
thetasm- Number of bins to use in smoothing theta
rsqr- Cut on distance between two points
pco- Number of sigma above average peak height for threshold (peak cut-off)
loops- Max number of multi-Hough loops

Definition at line 109 of file MultiHough2P.cxx.

117  :
118  fView(v),
119  fRhoSz(rhosz),
120  fThetaSz(thetasz),
121  fXwinSz(rhosm),
122  fYwinSz(thetasm),
123  fSigma(3.0/sqrt(12.)),
124  fRsqr(rsqr),
125  fPco(pco),
126  fLoops(loops),
127  fRmDst(rmdst),
128  fHspaceW(0),
129  fH(v, 0, 0, 0, 0)
130  { }
double fRhoSz
Size of rho bins (cm)
Definition: MultiHough2P.h:79
LiteTH2 * fHspaceW
Hough transform.
Definition: MultiHough2P.h:90
double fRsqr
Distance^2 cut on points (cm^2)
Definition: MultiHough2P.h:84
T sqrt(T number)
Definition: d0nt_math.hpp:156
int fXwinSz
Smoothing size (bins) in x.
Definition: MultiHough2P.h:81
int fYwinSz
Smoothing size (bins) in y.
Definition: MultiHough2P.h:82
double fSigma
Assumed location error on points (cm)
Definition: MultiHough2P.h:83
double fThetaSz
Size of theta bins (degrees)
Definition: MultiHough2P.h:80
rb::HoughResult fH
Results of Hough transformation;.
Definition: MultiHough2P.h:91
double fPco
Number of sigma abouve average peak height to use as cutoff.
Definition: MultiHough2P.h:85
int fLoops
Max number of Multi-Hough loops.
Definition: MultiHough2P.h:86
double fRmDst
Distance cut for removing points from a Hough line.
Definition: MultiHough2P.h:87
geo::View_t fView
Which detector view?
Definition: MultiHough2P.h:78
hough::MultiHough2P::~MultiHough2P ( )

Definition at line 168 of file MultiHough2P.cxx.

References fHspaceW.

169  {
170  if (fHspaceW) { delete fHspaceW; fHspaceW = 0;}
171  }
LiteTH2 * fHspaceW
Hough transform.
Definition: MultiHough2P.h:90

Member Function Documentation

void hough::MultiHough2P::BuildMap ( double  xlo,
double  xhi,
double  ylo,
double  yhi 
)
private

Definition at line 175 of file MultiHough2P.cxx.

References fHspaceW, fRhoSz, fThetaSz, makeTrainCVSamples::int, M_PI, util::pythag(), and r().

Referenced by Map().

177  {
178  double r; // Size of hough map in rho
179  int nr; // Number of bins to use in rho
180  int nth; // Number of bins to use in theta
181 
182  // 0.55 gives a little buffer beyond the minimum required 0.5
183  r = 0.55*util::pythag((xhi-xlo),(yhi-ylo));
184  nr = (int)rint(r/fRhoSz+0.5);
185 
186  nth = (int)rint(180.0/fThetaSz+0.5);
187 
188  // Allocate the Hough space histogram
189  if (fHspaceW!=0) { delete fHspaceW; fHspaceW=0; }
190  fHspaceW = new LiteTH2(nr, -r, r,
191  nth,-0.5*M_PI/180.0,179.5*M_PI/180.0);
192  }
double fRhoSz
Size of rho bins (cm)
Definition: MultiHough2P.h:79
LiteTH2 * fHspaceW
Hough transform.
Definition: MultiHough2P.h:90
#define M_PI
Definition: SbMath.h:34
double fThetaSz
Size of theta bins (degrees)
Definition: MultiHough2P.h:80
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
TRandom3 r(0)
void hough::MultiHough2P::Map ( const rb::HitList h)

Make the Hough map for the hit list h.

Parameters
h- The list of hits (can be mixed X and Y view hits)

Definition at line 225 of file MultiHough2P.cxx.

References abs(), hough::LiteTH2::AddBinContent(), BuildMap(), fH, fHspaceW, hough::LiteTH2::FindBinX(), hough::LiteTH2::FindBinY(), fRsqr, rb::HoughResult::fTNShi, rb::HoughResult::fTNSlo, fView, fXwinSz, rb::HoughResult::fXY0, fYwinSz, rb::HoughResult::fZ0, hough::Gaus(), hough::LiteTH2::GetBinCenterX(), hough::LiteTH2::GetBinCenterY(), hough::LiteTH2::GetNbinsX(), hough::LiteTH2::GetNbinsY(), MECModelEnuComparisons::i, calib::j, RhoTheta(), SmoothMap(), util::sqr(), chisquared::theta, w, submit_syst::x, xhi, make_syst_table_plots::xlo, submit_syst::y, yhi, and ylo.

Referenced by MultiMap().

226  {
227  unsigned int i;
228  unsigned int j;
229 
230  //
231  // Find a box that completely contains all the hits. A note on
232  // coordinates: Since the literature on Hough transforms uses an
233  // "x-y" coodinate system I map NOvA z to Hough x and NOvA X/Y to
234  // Hough y so that the formalizm matches what one would find in the
235  // literature.
236  //
237  double tlo=9E9, thi=-9e9;
238  double xlo=9e9, xhi=-9e9;
239  double ylo=9e9, yhi=-9e9;
240  for (i=0; i<h.size(); ++i) {
241  if (h[i].fView == fView) {
242  if (h[i].fZ<xlo) xlo = h[i].fZ;
243  if (h[i].fZ>xhi) xhi = h[i].fZ;
244  if (h[i].fXY<ylo) ylo = h[i].fXY;
245  if (h[i].fXY>yhi) yhi = h[i].fXY;
246  if (h[i].fT<tlo) tlo = h[i].fT;
247  if (h[i].fT>tlo) thi = h[i].fT;
248  }
249  }
250  fH.fTNSlo = tlo;
251  fH.fTNShi = thi;
252  fH.fZ0 = 0.5*(xlo+xhi);
253  fH.fXY0 = 0.5*(ylo+yhi);
254 
255 
256 
257  //
258  // Fill the hough map
259  //
260  this->BuildMap(xlo,xhi,ylo,yhi);
261 
262  for (i=0; i<h.size(); ++i) {
263  if (h[i].fView != fView) continue;
264  for (j=i+1; j<h.size(); ++j) {
265  if (h[j].fView != fView) continue;
266  const bool hasweight = h[i].fW != 0.0 && h[j].fW != 0.0;
267  if(!hasweight) continue;
268 
269  // Check distance between hit pair
270  double rsqr = util::sqr(h[i].fZ-h[j].fZ) + util::sqr(h[i].fXY-h[j].fXY);
271 
272  if (rsqr>fRsqr) continue;
273 
274  // don't let points in the same column vote unless they are close together
275  // if ( abs(h[i].fZ - h[j].fZ) < 5.9 && rsqr > fRsqr/4.0 ) continue;
276 
277  // don't let points in the same row vote unless they are far apart
278  if ( abs(h[i].fXY - h[j].fXY) < 3.9 && rsqr < fRsqr/4.0 ) continue;
279 
280  double rho, theta;
281  double sigmarho, sigmatheta;
282  this->RhoTheta(h[i].fZ-fH.fZ0, h[i].fXY-fH.fXY0,
283  h[j].fZ-fH.fZ0, h[j].fXY-fH.fXY0,
284  &rho, &theta,
285  &sigmarho, &sigmatheta);
286  //
287  // Fill region weighted by errors
288  //
289  int ixlo = fHspaceW->FindBinX(rho-3.0*sigmarho);
290  int ixhi = fHspaceW->FindBinX(rho+3.0*sigmarho);
291  int iylo = fHspaceW->FindBinY(theta-3.0*sigmatheta);
292  int iyhi = fHspaceW->FindBinY(theta+3.0*sigmatheta);
293 
294  if (ixlo<1) ixlo = 1;
295  if (iylo<1) iylo = 1;
296  if (ixhi>fHspaceW->GetNbinsX()) ixhi = fHspaceW->GetNbinsX();
297  if (iyhi>fHspaceW->GetNbinsY()) iyhi = fHspaceW->GetNbinsY();
298 
299  int ix, iy;
300 
301  // Precalcuate gaussian weights since they will be the same every loop
302  double yprob[iyhi]; // won't use first elements: use memory to save time
303  for (iy=iylo; iy<=iyhi; ++iy) {
304  const double y = fHspaceW->GetBinCenterY(iy);
305  yprob[iy] = Gaus(y-theta, sigmatheta);
306  }
307 
308  for (ix=ixlo; ix<=ixhi; ++ix) {
309  double x = fHspaceW->GetBinCenterX(ix);
310  const double xprob = Gaus(x-rho, sigmarho);
311  for (iy=iylo; iy<=iyhi; ++iy) {
312  double w = 0.0;
313 
314  // if(h[i].fW > 0.1 && h[j].fW > 0.1) w = 2.0*h[i].fW*h[j].fW/(h[i].fW+h[j].fW);
315  w = 2.0*h[i].fW*h[j].fW/(h[i].fW+h[j].fW) *
316  xprob * yprob[iy];
317  fHspaceW->AddBinContent(ix,iy,w);
318 
319  }
320  }
321  }
322  }
323  if (fXwinSz > 0 || fYwinSz > 0) this->SmoothMap();
324 
325  }
double fXY0
X/Y offset applied to hits.
Definition: HoughResult.h:59
int FindBinY(float y) const
double fZ0
Z offset applied to hits.
Definition: HoughResult.h:60
LiteTH2 * fHspaceW
Hough transform.
Definition: MultiHough2P.h:90
double fRsqr
Distance^2 cut on points (cm^2)
Definition: MultiHough2P.h:84
int fXwinSz
Smoothing size (bins) in x.
Definition: MultiHough2P.h:81
int fYwinSz
Smoothing size (bins) in y.
Definition: MultiHough2P.h:82
void abs(TH1 *hist)
float GetBinCenterX(int ix) const
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
int FindBinX(float x) const
rb::HoughResult fH
Results of Hough transformation;.
Definition: MultiHough2P.h:91
const double j
Definition: BetheBloch.cxx:29
static Double_t Gaus(const Double_t x, const Double_t sigma)
float GetBinCenterY(int iy) const
double fTNShi
Upper edge of time slice transformed.
Definition: HoughResult.h:58
int GetNbinsX() const
int GetNbinsY() const
double fTNSlo
Low edge of time slice transformed.
Definition: HoughResult.h:57
void BuildMap(double xlo, double xhi, double ylo, double yhi)
void RhoTheta(double x1, double y1, double x2, double y2, double *rho, double *theta, double *sigmarho, double *sigmatheta)
void AddBinContent(int ix, int iy, float val)
geo::View_t fView
Which detector view?
Definition: MultiHough2P.h:78
Float_t w
Definition: plot.C:20
TH2F * hough::MultiHough2P::MapToTH2F ( ) const

Definition at line 196 of file MultiHough2P.cxx.

References fHspaceW, fView, geo::kX, geo::kY, runNovaSAM::ret, and hough::LiteTH2::ToTH2F().

Referenced by hough::MultiHoughT::WriteHoughHistos().

197  {
198  if(!fHspaceW) return 0;
199  TH2F* ret = fHspaceW->ToTH2F();
200  ret->SetTitle(";rho [cm];theta [radians]");
201 
202  const char* tiw = 0;
203  if(fView == geo::kX) tiw = "fHough2P_Hx";
204  else if (fView == geo::kY) tiw = "fHough2P_Hy";
205  else abort();
206  ret->SetTitle(tiw);
207 
208  return ret;
209  }
LiteTH2 * fHspaceW
Hough transform.
Definition: MultiHough2P.h:90
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
TH2F * ToTH2F() const
geo::View_t fView
Which detector view?
Definition: MultiHough2P.h:78
void hough::MultiHough2P::MultiMap ( rb::HitList  h)

Definition at line 329 of file MultiHough2P.cxx.

References a, b, stan::math::fabs(), fH, fHspaceW, fLoops, fPco, rb::HoughResult::fPeak, hough::LiteTH2::GetBinContent(), hough::LiteTH2::GetNbinsX(), hough::LiteTH2::GetNbinsY(), MECModelEnuComparisons::i, Map(), genie::units::nb, r(), RefinePeak(), ReweightHits(), rb::HoughResult::SortPeaks(), std::sqrt(), sr, and chisquared::theta.

Referenced by hough::MultiHoughT::produce().

330  {
331 
332  int flag = 0; // flag to jump out of multi hough loop
333  int mh = 1; // count of number of multihough loops
334 
335  // make the initial Hough Map
336  this->Map(h);
337 
338  // calculate the threshold (needs to be done once at the beginning so that it is
339  // the same for all subsequent maps)
340  double Pave = 0.0; // average peak height
341  double Ps = 0.0; // standard deviation of average peak height
342  double th = 0.0; // threshold to use for peak cutoff
343  int nb = 0; // number of bins counted for avergae peak height
344 
345  // calculate average peak height to determine threshold
346  // empty bins are currently counted
347  for (int a = 1; a <= fHspaceW->GetNbinsX(); ++a) {
348  for (int b = 1; b <= fHspaceW->GetNbinsY(); ++b) {
349  double hab = fHspaceW->GetBinContent(a, b);
350  if (hab != 0) {
351  Pave += hab;
352  Ps += hab*hab;
353  } // end if
354  nb++;
355  } // end for b
356  } // end for a
357 
358  Pave = Pave/nb;
359  Ps = sqrt(Ps/nb - Pave*Pave);
360  th = Pave + fPco*Ps;
361 
362  // If the threshold is zero, then the initial Hough map is empty so there
363  // is no reason to proceed.
364  if(th == 0.0) return;
365 
366  // find the tallest peak
367  int ix, iy;
368  int xpeak = 0, ypeak = 0;
369  double peak = 0.0;
370  double r = 0.0;
371  double theta = 0.0;
372  double p;
373  for(ix = 1; ix <= fHspaceW->GetNbinsX(); ++ix) {
374  for(iy = 1; iy <= fHspaceW->GetNbinsY(); ++iy) {
375  p = fHspaceW->GetBinContent(ix,iy);
376  if (p>peak) {
377  peak = p;
378  xpeak = ix;
379  ypeak = iy;
380  } // end if p>peak
381  } // loop on y bins
382  } // loop on x bins
383 
384  // Average bins around the peak to get a better rho and theta value
385  this->RefinePeak(r, theta, xpeak, ypeak);
386 
387  // AS OF RIGHT NOW, SIGMA-RHO AND SIGMA-THETA ARE NOT CALCULATED!!!
388  double sr = 0.0;
389  double st = 0.0;
390 
391  // put tallest peak into Hough Results
392  if(peak > 0.0) fH.fPeak.push_back(rb::HoughPeak(r,theta,peak,th,sr,st));
393 
394  // remove hits on the line by reweighting them
395  this->ReweightHits(r, theta, h);
396 
397 
398 
399 
400  //
401  // start the MultiHough loop
402  //
403 
404  // if either the number of desired loops is only one or the tallest peak was already below threshold,
405  // then skip the multi-hough loop because we are already done
406  if(peak < th) flag = 1;
407  if(mh >= fLoops) flag = 1;
408 
409  while(flag == 0) {
410 
411  mh++;
412  if(mh >= fLoops) flag = 1;
413 
414  // repeat everything done above EXCEPT: calculating the threshold
415 
416  // build the new map
417  this->Map(h);
418 
419 
420 
421  // find the tallest peak
422  xpeak = 0;
423  ypeak = 0;
424  peak = 0.0;
425  r = 0.0;
426  theta = 0.0;
427  p = 0.0;
428 
429  for (ix=1; ix<=fHspaceW->GetNbinsX(); ++ix) {
430  for (iy=1; iy<=fHspaceW->GetNbinsY(); ++iy) {
431  p = fHspaceW->GetBinContent(ix,iy);
432  if (p>peak) {
433  peak = p;
434  xpeak = ix;
435  ypeak = iy;
436  } // end if p>peak
437  } // loop on y bins
438  } // loop on x bins
439 
440  this->RefinePeak(r, theta, xpeak, ypeak);
441 
442  sr = 0.0;
443  st = 0.0;
444 
445 
446 
447  // Run parallel line test: (If two lines have almost the same rho and theta values,
448  // they are assumed to be part of the same line and that line is NOT put into the
449  // Hough results.)
450  //
451  // NOTE: This is a little sensitive to the numbers 15.0 and 0.02.
452  int parflag = 0;
453  for(unsigned int i = 0; i < fH.fPeak.size(); ++i) {
454  if(fabs(r - fH.fPeak[i].fRho) < 15.0 && fabs(theta - fH.fPeak[i].fTheta) < 0.02) parflag = 1;
455  }
456 
457  // put only the tallest peak into Houghresults
458  if(parflag != 1 && peak >= th) fH.fPeak.push_back(rb::HoughPeak(r,theta,peak,th,sr,st));
459  if(peak < th) flag = 1;
460 
461  this->ReweightHits(r, theta, h);
462 
463  } // end MultiHough loop
464 
465  fH.SortPeaks();
466 
467  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
LiteTH2 * fHspaceW
Hough transform.
Definition: MultiHough2P.h:90
std::vector< HoughPeak > fPeak
List of peaks found in Hough space.
Definition: HoughResult.h:61
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SortPeaks()
Sort Hough peaks from best (highest) to worst.
Definition: HoughResult.cxx:83
void ReweightHits(double rho, double theta, rb::HitList &h)
rb::HoughResult fH
Results of Hough transformation;.
Definition: MultiHough2P.h:91
Data for a single peak in Hough space.
Definition: HoughResult.h:18
double fPco
Number of sigma abouve average peak height to use as cutoff.
Definition: MultiHough2P.h:85
const double a
int fLoops
Max number of Multi-Hough loops.
Definition: MultiHough2P.h:86
void Map(const rb::HitList &h)
caf::StandardRecord * sr
float GetBinContent(int ix, int iy) const
int GetNbinsX() const
int GetNbinsY() const
const hit & b
Definition: hits.cxx:21
TRandom3 r(0)
void RefinePeak(double &rho, double &theta, int xp, int yp)
static const double nb
Definition: Units.h:89
void hough::MultiHough2P::RefinePeak ( double &  rho,
double &  theta,
int  xp,
int  yp 
)
private

Definition at line 532 of file MultiHough2P.cxx.

References dist, stan::math::fabs(), fHspaceW, hough::LiteTH2::GetBinCenterX(), hough::LiteTH2::GetBinCenterY(), hough::LiteTH2::GetBinContent(), hough::LiteTH2::GetNbinsX(), hough::LiteTH2::GetNbinsY(), and util::pythag().

Referenced by MultiMap().

533  {
534  // make these fcl parameters???
535  int nx = 3, ny = 3; // +/- number of bins around peak to use for averaging
536  double peaktotal = 0.0;
537  double BinValue = 0.0;
538  int nXbins = fHspaceW->GetNbinsX();
539  int nYbins = fHspaceW->GetNbinsY();
540 
541 
542  // determine if bins to be averaged over will be off the map and constrain nx and ny accordingly
543  if((xp+nx) >= nXbins) nx = (nXbins-1) - xp;
544  if((xp-nx) < 0) nx = xp;
545 
546  if((yp+ny) >= nYbins) ny = (nYbins-1) - yp;
547  if((yp-ny) < 0) ny = yp;
548 
549 
550 
551  double xd = 0.0, yd = 0.0;
552  double dist = 0.0;
553 
554  for(int xi = -nx; xi <= nx; ++xi) {
555  xd = fabs(xi);
556  for(int yi = -ny; yi <= ny; ++yi) {
557 
558  yd = fabs(yi);
559  dist = util::pythag(xd,yd);
560  if(dist == 0.0) dist = 0.707107;
561 
562  BinValue = (fHspaceW->GetBinContent(xp+xi,yp+yi))/dist;
563  peaktotal += BinValue;
564  rho += BinValue*fHspaceW->GetBinCenterX(xp+xi);
565  theta += BinValue*fHspaceW->GetBinCenterY(yp+yi);
566 
567  }
568  }
569 
570  rho = rho/peaktotal;
571  theta = theta/peaktotal;
572 
573  }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
LiteTH2 * fHspaceW
Hough transform.
Definition: MultiHough2P.h:90
float GetBinCenterX(int ix) const
double dist
Definition: runWimpSim.h:113
float GetBinContent(int ix, int iy) const
float GetBinCenterY(int iy) const
int GetNbinsX() const
int GetNbinsY() const
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void hough::MultiHough2P::ReweightHits ( double  rho,
double  theta,
rb::HitList h 
)
private

Definition at line 579 of file MultiHough2P.cxx.

References a, b, plot_validation_datamc::c, std::cos(), d, geo::DsqrToLine(), fH, fRmDst, fView, rb::HoughResult::fXY0, rb::HoughResult::fZ0, MECModelEnuComparisons::i, calib::j, m, std::sin(), std::sqrt(), and make_true_q0q3_plots::zmax.

Referenced by MultiMap().

580  {
581  //
582  // This function does several things:
583  // 1. It sets the weights of all hits within 6 cm of the line to zero.
584  // 2. Then it resets the weights of the end points of the line to one.
585  // 3. Then it applies more scrubbing to remove stray hits that have no neighbors.
586  //
587 
588  // These variable were used in an old scheme in which I used a Fermi distribution function
589  // to determine the new hit weights based on their distance to the line. I left it in for
590  // potential future use.
591  // double dc = 4.0; // distance cutoff
592  // double ds = 1.0; // distance spread
593 
594 
595  // loop over all points, determine which are on that line and set their weights to zero
596  double m = 0.0;
597  double b = 0.0;
598  double s = sin(theta);
599  double c = cos(theta);
600 
601  // calculate slope/intercept in detector coordinates
602  if (s == 0.0) s = 1.0E-9;
603  m = -c/s;
604  b = r/s - m*fH.fZ0 + fH.fXY0;
605 
606  double d = 0.0; // perpendicular distance of point to line
607 
608  // numbers needed for determining which points are end points
609  double zmin = 1.0e9, zmax = -1.0e9;
610  unsigned int amax = 999999;
611  unsigned int amin = 999999;
612  // amax2 & amin2 are used for the case that there are 2 end points (when there are two points
613  // in the same plane within fRmDist (cm) of the line.)
614  // NOTE: There could be three points that fit this, but I only allow for keeping two.
615  unsigned int amax2 = 999999;
616  unsigned int amin2 = 999999;
617 
618  for(unsigned int a = 0; a < h.size(); ++a) {
619 
620  // if a point is on the line, set its weight to zero
621  if(h[a].fView == fView) {
622 
623  d = geo::DsqrToLine(h[a].fZ, h[a].fXY, 0.0, b, 100.0, 100.0*m + b);
624  d = sqrt(d);
625  // h[a].fW *= 1.0 - 1.0/(exp((d - dc)/ds) + 1.0); // The Fermi weight function
626  if(d < fRmDst) {
627  h[a].fW = 0.0;
628 
629  // determine end points
630  if(h[a].fZ == zmax) { amax2 = a; }
631  if(h[a].fZ == zmin) { amin2 = a; }
632  if(h[a].fZ > zmax) { zmax = h[a].fZ; amax = a; }
633  if(h[a].fZ < zmin) { zmin = h[a].fZ; amin = a; }
634  }
635 
636  } // end if fView
637  } // end for a
638 
639  // Reset the weights of the line end points to one.
640  if(amax != 999999) h[amax].fW = 1.0;
641  if(amin != 999999) h[amin].fW = 1.0;
642  if(amax2 != 999999) h[amax2].fW = 1.0;
643  if(amin2 != 999999) h[amin2].fW = 1.0;
644 
645 
646 
647  // Rescrub the hit list to remove points with no neighbors within sqrt(ScrubDistSq) (cm).
648  double dminsq = 0.0, dsq = 0.0;
649  double ScrubDistSq = 900.0;
650 
651  // We will keep track of which points are within sqrt(ScrubDistSq) (cm) of
652  // another point (i.e. - which points have a buddy) to make this more efficient.
653  std::vector<int> BuddyList;
654  for(unsigned int a = 0; a < h.size(); ++a) BuddyList.push_back(0);
655 
656  // For each point, loop over all others to determine if it has a "buddy" within sqrt(ScrubDistSq).
657  // If it has no "buddies," scrub it! (Set the weight to zero.)
658  for(unsigned int i = 0; i < h.size(); ++i) {
659  if(BuddyList[i] > 0) continue;
660  if(h[i].fView != fView) continue;
661  if(h[i].fW == 0.0) continue;
662  dminsq = 1.0e9;
663  for(unsigned int j = 0; j < h.size(); ++j) {
664  if(i == j) continue;
665  if(h[j].fView != fView) continue;
666  if(h[j].fW == 0.0) continue;
667  dsq = (h[i].fZ - h[j].fZ)*(h[i].fZ - h[j].fZ) + (h[i].fXY - h[j].fXY)*(h[i].fXY - h[j].fXY);
668  if(dsq < ScrubDistSq) {
669  BuddyList[i] += 1;
670  BuddyList[j] += 1;
671  }
672  if(dsq < dminsq) dminsq = dsq;
673  }
674  if(dminsq >= ScrubDistSq) h[i].fW = 0.0;
675  }
676 
677  }
double fXY0
X/Y offset applied to hits.
Definition: HoughResult.h:59
double fZ0
Z offset applied to hits.
Definition: HoughResult.h:60
T sqrt(T number)
Definition: d0nt_math.hpp:156
const XML_Char * s
Definition: expat.h:262
rb::HoughResult fH
Results of Hough transformation;.
Definition: MultiHough2P.h:91
const double a
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
Definition: Geo.cxx:338
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
double fRmDst
Distance cut for removing points from a Hough line.
Definition: MultiHough2P.h:87
T sin(T number)
Definition: d0nt_math.hpp:132
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
geo::View_t fView
Which detector view?
Definition: MultiHough2P.h:78
void hough::MultiHough2P::RhoTheta ( double  x1,
double  y1,
double  x2,
double  y2,
double *  rho,
double *  theta,
double *  sigmarho,
double *  sigmatheta 
)
private

Definition at line 134 of file MultiHough2P.cxx.

References std::atan2(), std::cos(), d, fSigma, M_PI, util::pythag(), std::sin(), submit_syst::x2, and submit_syst::y2.

Referenced by Map().

138  {
139  *theta = atan2(x1-x2,y2-y1);
140  *rho = 0.5*(cos(*theta)*(x1+x2)+sin(*theta)*(y1+y2));
141  //
142  // Keep rho and theta within range [0,pi) (-R,+R)
143  //
144  while (*theta<0) {
145  *theta += M_PI;
146  *rho = -(*rho);
147  }
148  while (*theta>=M_PI) {
149  *theta -= M_PI;
150  *rho = -(*rho);
151  }
152 
153  // Compute the error estimates
154  double d = util::pythag((x2-x1),(y2-y1));
155  *sigmarho = fSigma;
156  *sigmatheta = M_SQRT2*fSigma/d;
157 
158  // these were failed attempts at making sigmatheta level off at a fixed value
159  //*sigmatheta = M_SQRT2*fSigma/d + 0.1;
160  /*
161  if(d < 20.0) *sigmatheta = M_SQRT2*fSigma/d;
162  else *sigmatheta = M_SQRT2*fSigma/20.0;
163  */
164  }
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
double fSigma
Assumed location error on points (cm)
Definition: MultiHough2P.h:83
#define M_PI
Definition: SbMath.h:34
Float_t d
Definition: plot.C:236
T sin(T number)
Definition: d0nt_math.hpp:132
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
T cos(T number)
Definition: d0nt_math.hpp:78
T atan2(T number)
Definition: d0nt_math.hpp:72
void hough::MultiHough2P::SmoothMap ( )
private

Definition at line 471 of file MultiHough2P.cxx.

References a, b, fHspaceW, fXwinSz, fYwinSz, hough::Gaus(), hough::LiteTH2::GetBinContent(), hough::LiteTH2::GetNbinsX(), hough::LiteTH2::GetNbinsY(), MECModelEnuComparisons::i, calib::j, m, getGoodRuns4SAM::n, hough::LiteTH2::SetBinContent(), sum, and w.

Referenced by Map().

472  {
473 
474  // NOTE: Since most of the Hough space map is empty, we only apply smoothing to the non-zero
475  // bins in order to save lots of calculation time.
476 
477  int i, j, m, n;
478  int nx, ny;
479  int ix, iy;
480  double sum;
481  int wx, wy;
482 
483  wx = 2*fXwinSz + 1;
484  wy = 2*fYwinSz + 1;
485 
486  double w[wx][wy];
487 
488  // fill the smoothing weight map
489  for(int a = 0; a < wx; a++) {
490  for(int b = 0; b < wy; b++) {
491 
492  w[a][b] = Gaus(a - fXwinSz, fXwinSz/3.0)*Gaus(b - fYwinSz, fYwinSz/3.0);
493 
494  } // end for b
495  } // end for a
496 
497  LiteTH2* htmp = new LiteTH2(*fHspaceW);
498 
499  nx = htmp->GetNbinsX();
500  ny = htmp->GetNbinsY();
501  for (i=1; i<=nx; ++i) {
502  for(j=1; j<=ny; ++j) {
503  sum = 0.0;
504  // skip smoothing if bin is empty
505  if(htmp->GetBinContent(i,j) == 0.0) continue;
506  for (m=-fXwinSz; m<=fXwinSz; m++) {
507  ix = i+m;
508  // Don't wrap around in rho
509  if (ix>nx) continue;
510  if (ix<1) continue;
511  for (n=-fYwinSz; n<=fYwinSz; ++n) {
512  iy = j+n;
513  // Wrap in theta
514  if (iy>ny) iy = iy%ny;
515  if (iy<1) iy = ny+iy;
516 
517  sum += w[m + fXwinSz][n + fYwinSz]*htmp->GetBinContent(ix,iy);
518  }
519  }
520  fHspaceW->SetBinContent(i,j,sum);
521  }
522  }
523  delete htmp;
524  htmp = 0;
525 
526  }
LiteTH2 * fHspaceW
Hough transform.
Definition: MultiHough2P.h:90
int fXwinSz
Smoothing size (bins) in x.
Definition: MultiHough2P.h:81
int fYwinSz
Smoothing size (bins) in y.
Definition: MultiHough2P.h:82
const double a
const double j
Definition: BetheBloch.cxx:29
void SetBinContent(int ix, int iy, float val)
static Double_t Gaus(const Double_t x, const Double_t sigma)
const hit & b
Definition: hits.cxx:21
Double_t sum
Definition: plot.C:31
Float_t w
Definition: plot.C:20

Member Data Documentation

rb::HoughResult hough::MultiHough2P::fH

Results of Hough transformation;.

Definition at line 91 of file MultiHough2P.h.

Referenced by Map(), MultiMap(), hough::MultiHoughT::produce(), and ReweightHits().

LiteTH2* hough::MultiHough2P::fHspaceW

Hough transform.

Definition at line 90 of file MultiHough2P.h.

Referenced by BuildMap(), Map(), MapToTH2F(), MultiMap(), RefinePeak(), SmoothMap(), and ~MultiHough2P().

int hough::MultiHough2P::fLoops

Max number of Multi-Hough loops.

Definition at line 86 of file MultiHough2P.h.

Referenced by MultiMap().

double hough::MultiHough2P::fPco

Number of sigma abouve average peak height to use as cutoff.

Definition at line 85 of file MultiHough2P.h.

Referenced by MultiMap().

double hough::MultiHough2P::fRhoSz

Size of rho bins (cm)

Definition at line 79 of file MultiHough2P.h.

Referenced by BuildMap().

double hough::MultiHough2P::fRmDst

Distance cut for removing points from a Hough line.

Definition at line 87 of file MultiHough2P.h.

Referenced by ReweightHits().

double hough::MultiHough2P::fRsqr

Distance^2 cut on points (cm^2)

Definition at line 84 of file MultiHough2P.h.

Referenced by Map().

double hough::MultiHough2P::fSigma

Assumed location error on points (cm)

Definition at line 83 of file MultiHough2P.h.

Referenced by RhoTheta().

double hough::MultiHough2P::fThetaSz

Size of theta bins (degrees)

Definition at line 80 of file MultiHough2P.h.

Referenced by BuildMap().

geo::View_t hough::MultiHough2P::fView

Which detector view?

Definition at line 78 of file MultiHough2P.h.

Referenced by Map(), MapToTH2F(), and ReweightHits().

int hough::MultiHough2P::fXwinSz

Smoothing size (bins) in x.

Definition at line 81 of file MultiHough2P.h.

Referenced by Map(), and SmoothMap().

int hough::MultiHough2P::fYwinSz

Smoothing size (bins) in y.

Definition at line 82 of file MultiHough2P.h.

Referenced by Map(), and SmoothMap().


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