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

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

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-23/HoughTransform/Hough2P.h"

Public Member Functions

 Hough2P (geo::View_t v, double rhosz, double thetasz, int rhosm, int thetasm, double rsqr, double pco)
 Construct the 2-point Hough transform. More...
 
 ~Hough2P ()
 
void Map (const 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...
 
TH2F * fHspaceW
 Hough transform. More...
 
IslandsTH2fHspaceI
 Islands for x view. 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 FindPeaks ()
 

Detailed Description

Hough transform based on 2-point combinations.

Definition at line 18 of file Hough2P.h.

Constructor & Destructor Documentation

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

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)

Definition at line 23 of file Hough2P.cxx.

29  :
30  fView(v),
31  fRhoSz(rhosz),
32  fThetaSz(thetasz),
33  fXwinSz(rhosm),
34  fYwinSz(thetasm),
35  fSigma(3.0/sqrt(12.)),
36  fRsqr(rsqr),
37  fPco(pco),
38  fHspaceW(0),
39  fHspaceI(0),
40  fH(v, 0, 0, 0, 0)
41  { }
geo::View_t fView
Which detector view?
Definition: Hough2P.h:66
double fRsqr
Distance^2 cut on points (cm^2)
Definition: Hough2P.h:72
int fYwinSz
Smoothing size (bins) in y.
Definition: Hough2P.h:70
double fThetaSz
Size of theta bins (degrees)
Definition: Hough2P.h:68
IslandsTH2 * fHspaceI
Islands for x view.
Definition: Hough2P.h:77
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH2F * fHspaceW
Hough transform.
Definition: Hough2P.h:76
double fRhoSz
Size of rho bins (cm)
Definition: Hough2P.h:67
rb::HoughResult fH
Results of Hough transformation;.
Definition: Hough2P.h:78
int fXwinSz
Smoothing size (bins) in x.
Definition: Hough2P.h:69
double fSigma
Assumed location error on points (cm)
Definition: Hough2P.h:71
double fPco
Number of sigma abouve average peak height to use as cutoff.
Definition: Hough2P.h:73
hough::Hough2P::~Hough2P ( )

Definition at line 72 of file Hough2P.cxx.

References fHspaceI, and fHspaceW.

73  {
74  if (fHspaceW) { delete fHspaceW; fHspaceW = 0; }
75  if (fHspaceI) { delete fHspaceI; fHspaceI = 0; }
76  }
IslandsTH2 * fHspaceI
Islands for x view.
Definition: Hough2P.h:77
TH2F * fHspaceW
Hough transform.
Definition: Hough2P.h:76

Member Function Documentation

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

Definition at line 80 of file Hough2P.cxx.

References fHspaceI, fHspaceW, fPco, fRhoSz, fThetaSz, fView, makeTrainCVSamples::int, geo::kX, geo::kY, M_PI, util::pythag(), and r().

Referenced by Map().

82  {
83  double r; // Size of hough map in rho
84  int nr; // Number of bins to use in rho
85  int nth; // Number of bins to use in theta
86 
87  // 0.55 gives a little buffer beyond the minimum required 0.5
88  r = 0.55*util::pythag((xhi-xlo),(yhi-ylo));
89  nr = (int)rint(r/fRhoSz+0.5);
90 
91  nth = (int)rint(180.0/fThetaSz+0.5);
92 
93  const char* tiw=0;
94  const char* tii=0;
95  if (fView==geo::kX) { tiw = "fHough2P_Hx"; tii="Hough2P_Ix"; }
96  else if (fView==geo::kY) { tiw = "fHough2P_Hy"; tii="Hough2P_Iy"; }
97  else abort();
98 
99  // Allocate the Hough space histogram
100  if (fHspaceW!=0) { delete fHspaceW; fHspaceW=0; }
101  fHspaceW = new TH2F(tiw,
102  ";rho [cm];theta [radians]",
103  nr, -r, r,
104  nth,-0.5*M_PI/180.0,179.5*M_PI/180.0);
105 
106  // The islands histogram is closely related to the Hough space
107  // histogram. Rebuild it now to ensure they stay in sync.
108  if (fHspaceI!=0) { delete fHspaceI; fHspaceI=0; }
109  fHspaceI = new IslandsTH2(tii,
110  ";rho [cm];theta [radians]",
111  fHspaceW,
112  fPco,
113  false,
114  true);
115  }
geo::View_t fView
Which detector view?
Definition: Hough2P.h:66
double fThetaSz
Size of theta bins (degrees)
Definition: Hough2P.h:68
IslandsTH2 * fHspaceI
Islands for x view.
Definition: Hough2P.h:77
Vertical planes which measure X.
Definition: PlaneGeo.h:28
#define M_PI
Definition: SbMath.h:34
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
TH2F * fHspaceW
Hough transform.
Definition: Hough2P.h:76
double fRhoSz
Size of rho bins (cm)
Definition: Hough2P.h:67
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
double fPco
Number of sigma abouve average peak height to use as cutoff.
Definition: Hough2P.h:73
TRandom3 r(0)
void hough::Hough2P::FindPeaks ( )
private

Definition at line 210 of file Hough2P.cxx.

References fH, fHspaceI, fHspaceW, rb::HoughResult::fPeak, hough::IslandsTH2::LabelIslands(), hough::IslandsTH2::Nislands(), r(), rb::HoughResult::SortPeaks(), and chisquared::theta.

Referenced by Map().

211  {
212 
214 
215  unsigned int ipeak;
216  int ix, iy;
217  double peak = 0.0;
218  double r = 0.0;
219  double theta = 0.0;
220  double p;
221  for (ipeak=0; ipeak<fHspaceI->Nislands(); ++ipeak) {
222  peak = 0.0;
223  r = 0.0;
224  theta = 0.0;
225  for (ix=1; ix<=fHspaceI->GetNbinsX(); ++ix) {
226  for (iy=1; iy<=fHspaceI->GetNbinsY(); ++iy) {
227  if (fHspaceI->GetBinContent(ix,iy)==ipeak+1) {
228  p = fHspaceW->GetBinContent(ix,iy);
229  if (p>peak) {
230  peak = p;
231  r = fHspaceW->GetXaxis()->GetBinCenter(ix);
232  theta = fHspaceW->GetYaxis()->GetBinCenter(iy);
233  }
234  } // if peak matches
235  } // loop on y bins
236  } // loop on x bins
237 
238  // NOTE: The last three arguments below are threshold, sigma-rho, and sigma-theta for the Hough line
239  // which this Hough transform currently does NOT calculate (the island space calculates the
240  // threshold but it does NOT pass that value back to Hough2P.)
241  if(peak > 0.0) fH.fPeak.push_back(rb::HoughPeak(r,theta,peak,0.0,0.0,0.0));
242 
243  } // loop on peaks
244 
245  //
246  // Assert that there must always be at least one peak
247  //
248  if (fH.fPeak.size()==0) {
249  peak = 0.0;
250  r = 0.0;
251  theta = 0.0;
252  for (ix=1; ix<=fHspaceW->GetNbinsX(); ++ix) {
253  for (iy=1; iy<=fHspaceW->GetNbinsY(); ++iy) {
254  p = fHspaceW->GetBinContent(ix,iy);
255  if (p>peak) {
256  peak = p;
257  r = fHspaceW->GetXaxis()->GetBinCenter(ix);
258  theta = fHspaceW->GetYaxis()->GetBinCenter(iy);
259  }
260  }
261  }
262 
263  // NOTE: The last three arguments below are threshold, sigma-rho, and sigma-theta for the Hough line
264  // which this Hough transform currently does NOT calculate (the island space calculates the
265  // threshold but it does NOT pass that value back to Hough2P.)
266  if(peak > 0.0) fH.fPeak.push_back(rb::HoughPeak(r,theta,peak,0.0,0.0,0.0));
267  }
268 
269  fH.SortPeaks();
270  }
IslandsTH2 * fHspaceI
Islands for x view.
Definition: Hough2P.h:77
std::vector< HoughPeak > fPeak
List of peaks found in Hough space.
Definition: HoughResult.h:61
const char * p
Definition: xmltok.h:285
TH2F * fHspaceW
Hough transform.
Definition: Hough2P.h:76
void SortPeaks()
Sort Hough peaks from best (highest) to worst.
Definition: HoughResult.cxx:83
Data for a single peak in Hough space.
Definition: HoughResult.h:18
void LabelIslands(const TH2 *h)
Definition: IslandsTH2.cxx:37
rb::HoughResult fH
Results of Hough transformation;.
Definition: Hough2P.h:78
unsigned int Nislands() const
Definition: IslandsTH2.h:20
TRandom3 r(0)
void hough::Hough2P::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 119 of file Hough2P.cxx.

References BuildMap(), fH, fHspaceW, FindPeaks(), fRsqr, rb::HoughResult::fTNShi, rb::HoughResult::fTNSlo, fView, fXwinSz, rb::HoughResult::fXY0, fYwinSz, rb::HoughResult::fZ0, Gaus(), MECModelEnuComparisons::i, calib::j, RhoTheta(), SmoothMap(), util::sqr(), chisquared::theta, w, submit_syst::x, xaxis, xhi, make_syst_table_plots::xlo, submit_syst::y, yaxis, yhi, and ylo.

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

120  {
121  unsigned int i;
122  unsigned int j;
123 
124  //
125  // Find a box that completely contains all the hits. A note on
126  // coordinates: Since the literature on Hough transforms uses an
127  // "x-y" coodinate system I map NOvA z to Hough x and NOvA X/Y to
128  // Hough y so that the formalizm matches what one would find in the
129  // literature.
130  //
131  double tlo=9E9, thi=-9e9;
132  double xlo=9e9, xhi=-9e9;
133  double ylo=9e9, yhi=-9e9;
134  for (i=0; i<h.size(); ++i) {
135  if (h[i].fView == fView) {
136  if (h[i].fZ<xlo) xlo = h[i].fZ;
137  if (h[i].fZ>xhi) xhi = h[i].fZ;
138  if (h[i].fXY<ylo) ylo = h[i].fXY;
139  if (h[i].fXY>yhi) yhi = h[i].fXY;
140  if (h[i].fT<tlo) tlo = h[i].fT;
141  if (h[i].fT>tlo) thi = h[i].fT;
142  }
143  }
144  fH.fTNSlo = tlo;
145  fH.fTNShi = thi;
146  fH.fZ0 = 0.5*(xlo+xhi);
147  fH.fXY0 = 0.5*(ylo+yhi);
148 
149  //
150  // Fill the hough maps
151  //
152  this->BuildMap(xlo,xhi,ylo,yhi);
153  TAxis* xaxis = fHspaceW->GetXaxis();
154  TAxis* yaxis = fHspaceW->GetYaxis();
155  for (i=0; i<h.size(); ++i) {
156  if (h[i].fView != fView) continue;
157  for (j=i+1; j<h.size(); ++j) {
158  if (h[j].fView != fView) continue;
159 
160  // Check distance between hit pair
161  double rsqr = util::sqr(h[i].fZ-h[j].fZ) + util::sqr(h[i].fXY-h[j].fXY);
162  if (rsqr>fRsqr) continue;
163 
164  // don't let points in the same column vote unless they are close together
165  // if ( abs(h[i].fZ - h[j].fZ) < 6.0 && rsqr > fRsqr/4.0 ) continue;
166 
167  // don't let points in the same row vote unless they are far apart
168  // if ( abs(h[i].fXY - h[j].fXY) < 6.0 && rsqr < fRsqr/4.0 ) continue;
169 
170  double rho, theta;
171  double sigmarho, sigmatheta;
172  this->RhoTheta(h[i].fZ-fH.fZ0, h[i].fXY-fH.fXY0,
173  h[j].fZ-fH.fZ0, h[j].fXY-fH.fXY0,
174  &rho, &theta,
175  &sigmarho, &sigmatheta);
176  //
177  // Fill region weighted by errors
178  //
179  int ixlo = xaxis->FindBin(rho-3.0*sigmarho);
180  int ixhi = xaxis->FindBin(rho+3.0*sigmarho);
181  int iylo = yaxis->FindBin(theta-3.0*sigmatheta);
182  int iyhi = yaxis->FindBin(theta+3.0*sigmatheta);
183 
184  if (ixlo<1) ixlo = 1;
185  if (iylo<1) iylo = 1;
186  if (ixhi>xaxis->GetNbins()) ixhi = xaxis->GetNbins();
187  if (iyhi>yaxis->GetNbins()) iyhi = yaxis->GetNbins();
188 
189  int ix, iy;
190  for (ix=ixlo; ix<=ixhi; ++ix) {
191  double x = xaxis->GetBinCenter(ix);
192  for (iy=iylo; iy<=iyhi; ++iy) {
193  double y = yaxis->GetBinCenter(iy);
194  double w = 2.0*h[i].fW*h[j].fW/(h[i].fW+h[j].fW);
195  w *=
196  TMath::Gaus(x-rho, 0, sigmarho, kTRUE) *
197  TMath::Gaus(y-theta, 0, sigmatheta, kTRUE);
198  fHspaceW->Fill(x,y,w);
199  }
200  }
201  }
202  }
203  if (fXwinSz>0 || fYwinSz>0) this->SmoothMap();
204  this->FindPeaks();
205 
206  }
geo::View_t fView
Which detector view?
Definition: Hough2P.h:66
double fRsqr
Distance^2 cut on points (cm^2)
Definition: Hough2P.h:72
int fYwinSz
Smoothing size (bins) in y.
Definition: Hough2P.h:70
double fXY0
X/Y offset applied to hits.
Definition: HoughResult.h:59
double fZ0
Z offset applied to hits.
Definition: HoughResult.h:60
std::string yaxis
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
std::string xaxis
TH2F * fHspaceW
Hough transform.
Definition: Hough2P.h:76
void FindPeaks()
Definition: Hough2P.cxx:210
double Gaus(TH1D *h, double &err, bool isTruth)
Definition: absCal.cxx:489
const double j
Definition: BetheBloch.cxx:29
double fTNShi
Upper edge of time slice transformed.
Definition: HoughResult.h:58
double fTNSlo
Low edge of time slice transformed.
Definition: HoughResult.h:57
rb::HoughResult fH
Results of Hough transformation;.
Definition: Hough2P.h:78
int fXwinSz
Smoothing size (bins) in x.
Definition: Hough2P.h:69
void SmoothMap()
Definition: Hough2P.cxx:274
void BuildMap(double xlo, double xhi, double ylo, double yhi)
Definition: Hough2P.cxx:80
Float_t w
Definition: plot.C:20
void RhoTheta(double x1, double y1, double x2, double y2, double *rho, double *theta, double *sigmarho, double *sigmatheta)
Definition: Hough2P.cxx:45
void hough::Hough2P::RhoTheta ( double  x1,
double  y1,
double  x2,
double  y2,
double *  rho,
double *  theta,
double *  sigmarho,
double *  sigmatheta 
)
private

Definition at line 45 of file Hough2P.cxx.

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

Referenced by Map().

49  {
50  *theta = atan2(x1-x2,y2-y1);
51  *rho = 0.5*(cos(*theta)*(x1+x2)+sin(*theta)*(y1+y2));
52  //
53  // Keep rho and theta within range [0,pi) (-R,+R)
54  //
55  while (*theta<0) {
56  *theta += M_PI;
57  *rho = -(*rho);
58  }
59  while (*theta>=M_PI) {
60  *theta -= M_PI;
61  *rho = -(*rho);
62  }
63 
64  // Compute the error estimates
65  double d = util::pythag((x2-x1),(y2-y1));
66  *sigmarho = fSigma;
67  *sigmatheta = M_SQRT2*fSigma/d;
68  }
Float_t y1[n_points_granero]
Definition: compare.C:5
Float_t x1[n_points_granero]
Definition: compare.C:5
#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
double fSigma
Assumed location error on points (cm)
Definition: Hough2P.h:71
T cos(T number)
Definition: d0nt_math.hpp:78
T atan2(T number)
Definition: d0nt_math.hpp:72
void hough::Hough2P::SmoothMap ( )
private

Definition at line 274 of file Hough2P.cxx.

References a, b, fHspaceW, fXwinSz, fYwinSz, Gaus(), MECModelEnuComparisons::i, calib::j, m, getGoodRuns4SAM::n, sum, and w.

Referenced by Map().

275  {
276 
277  int i, j, m, n;
278  int nx, ny;
279  int ix, iy;
280  double sum;
281  int wx, wy;
282 
283  wx = 2*fXwinSz + 1;
284  wy = 2*fYwinSz + 1;
285 
286  double w[wx][wy];
287 
288  // fill the smoothing weight map
289  for(int a = 0; a < wx; a++) {
290  for(int b = 0; b < wy; b++) {
291 
292  w[a][b] = TMath::Gaus(a - fXwinSz, 0, fXwinSz/3.0, kTRUE)*TMath::Gaus(b - fYwinSz, 0, fYwinSz/3.0, kTRUE);
293 
294  } // end for b
295  } // end for a
296 
297  TH2F* htmp = new TH2F(*fHspaceW);
298  fHspaceW->Reset();
299 
300  nx = htmp->GetNbinsX();
301  ny = htmp->GetNbinsY();
302  for (i=1; i<=nx; ++i) {
303  for(j=1; j<=ny; ++j) {
304  sum = 0.0;
305  for (m=-fXwinSz; m<=fXwinSz; m++) {
306  ix = i+m;
307  // Don't wrap around in rho
308  if (ix>nx) continue;
309  if (ix<1) continue;
310  for (n=-fYwinSz; n<=fYwinSz; ++n) {
311  iy = j+n;
312  // Wrap in theta
313  if (iy>ny) iy = iy%ny;
314  if (iy<1) iy = ny+iy;
315 
316  sum += w[m + fXwinSz][n + fYwinSz]*htmp->GetBinContent(ix,iy);
317  }
318  }
319  fHspaceW->SetBinContent(i,j,sum);
320  }
321  }
322  delete htmp;
323  htmp = 0;
324 
325  }
int fYwinSz
Smoothing size (bins) in y.
Definition: Hough2P.h:70
TH2F * fHspaceW
Hough transform.
Definition: Hough2P.h:76
const double a
double Gaus(TH1D *h, double &err, bool isTruth)
Definition: absCal.cxx:489
const double j
Definition: BetheBloch.cxx:29
int fXwinSz
Smoothing size (bins) in x.
Definition: Hough2P.h:69
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::Hough2P::fH

Results of Hough transformation;.

Definition at line 78 of file Hough2P.h.

Referenced by FindPeaks(), Map(), and hough::HoughT::produce().

IslandsTH2* hough::Hough2P::fHspaceI

Islands for x view.

Definition at line 77 of file Hough2P.h.

Referenced by BuildMap(), FindPeaks(), hough::HoughT::WriteHoughHistos(), and ~Hough2P().

TH2F* hough::Hough2P::fHspaceW

Hough transform.

Definition at line 76 of file Hough2P.h.

Referenced by BuildMap(), FindPeaks(), Map(), SmoothMap(), hough::HoughT::WriteHoughHistos(), and ~Hough2P().

double hough::Hough2P::fPco

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

Definition at line 73 of file Hough2P.h.

Referenced by BuildMap().

double hough::Hough2P::fRhoSz

Size of rho bins (cm)

Definition at line 67 of file Hough2P.h.

Referenced by BuildMap().

double hough::Hough2P::fRsqr

Distance^2 cut on points (cm^2)

Definition at line 72 of file Hough2P.h.

Referenced by Map().

double hough::Hough2P::fSigma

Assumed location error on points (cm)

Definition at line 71 of file Hough2P.h.

Referenced by RhoTheta().

double hough::Hough2P::fThetaSz

Size of theta bins (degrees)

Definition at line 68 of file Hough2P.h.

Referenced by BuildMap().

geo::View_t hough::Hough2P::fView

Which detector view?

Definition at line 66 of file Hough2P.h.

Referenced by BuildMap(), and Map().

int hough::Hough2P::fXwinSz

Smoothing size (bins) in x.

Definition at line 69 of file Hough2P.h.

Referenced by Map(), and SmoothMap().

int hough::Hough2P::fYwinSz

Smoothing size (bins) in y.

Definition at line 70 of file Hough2P.h.

Referenced by Map(), and SmoothMap().


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