Hough2P.cxx
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 /// \brief Perform a "2 point" Hough transform on a collection of hits
3 /// \authors messier@indiana.edu, mibaird@indiana.edu
4 /// \date
5 ////////////////////////////////////////////////////////////////////////
7 
8 #include <algorithm>
9 #include <cstdlib>
10 #include <cmath>
11 #include <iostream>
12 
13 #include "TH2F.h"
14 #include "TMath.h"
15 
17 #include "RecoBase/HitList.h"
18 #include "Utilities/func/MathUtil.h"
19 
20 namespace hough
21 {
22  //.....................................................................
24  double rhosz,
25  double thetasz,
26  int rhosm,
27  int thetasm,
28  double rsqr,
29  double pco) :
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  { }
42 
43  //......................................................................
44 
45  void Hough2P::RhoTheta(double x1, double y1,
46  double x2, double y2,
47  double* rho, double* theta,
48  double* sigmarho, double* sigmatheta)
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  }
69 
70  //......................................................................
71 
73  {
74  if (fHspaceW) { delete fHspaceW; fHspaceW = 0; }
75  if (fHspaceI) { delete fHspaceI; fHspaceI = 0; }
76  }
77 
78  //......................................................................
79 
80  void Hough2P::BuildMap(double xlo, double xhi,
81  double ylo, double yhi)
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  }
116 
117  //......................................................................
118 
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  }
207 
208  //......................................................................
209 
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  }
271 
272  //......................................................................
273 
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  }
326 
327 } // end namespace hough
328 
329 ////////////////////////////////////////////////////////////////////////
geo::View_t fView
Which detector view?
Definition: Hough2P.h:66
TAxis * xaxis
Definition: PlotUnfolding.C:64
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
Perform a "2 point" Hough transform on a collection of hits.
Definition: Hough2P.cxx:20
Float_t y1[n_points_granero]
Definition: compare.C:5
double fThetaSz
Size of theta bins (degrees)
Definition: Hough2P.h:68
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
IslandsTH2 * fHspaceI
Islands for x view.
Definition: Hough2P.h:77
Float_t x1[n_points_granero]
Definition: compare.C:5
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
Vertical planes which measure X.
Definition: PlaneGeo.h:28
TAxis * yaxis
Definition: PlotUnfolding.C:65
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
#define M_PI
Definition: SbMath.h:34
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
TH2F * fHspaceW
Hough transform.
Definition: Hough2P.h:76
void SortPeaks()
Sort Hough peaks from best (highest) to worst.
Definition: HoughResult.cxx:83
Just the essential information required for track fitting.
Definition: HitList.h:40
Data for a single peak in Hough space.
Definition: HoughResult.h:18
void FindPeaks()
Definition: Hough2P.cxx:210
const double a
double fRhoSz
Size of rho bins (cm)
Definition: Hough2P.h:67
Float_t d
Definition: plot.C:236
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
void LabelIslands(const TH2 *h)
Definition: IslandsTH2.cxx:37
T sin(T number)
Definition: d0nt_math.hpp:132
double fTNSlo
Low edge of time slice transformed.
Definition: HoughResult.h:57
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
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
unsigned int Nislands() const
Definition: IslandsTH2.h:20
double fPco
Number of sigma abouve average peak height to use as cutoff.
Definition: Hough2P.h:73
const hit & b
Definition: hits.cxx:21
void Map(const rb::HitList &h)
Definition: Hough2P.cxx:119
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
void SmoothMap()
Definition: Hough2P.cxx:274
Double_t sum
Definition: plot.C:31
void BuildMap(double xlo, double xhi, double ylo, double yhi)
Definition: Hough2P.cxx:80
Float_t w
Definition: plot.C:20
Hough2P(geo::View_t v, double rhosz, double thetasz, int rhosm, int thetasm, double rsqr, double pco)
Construct the 2-point Hough transform.
Definition: Hough2P.cxx:23
void RhoTheta(double x1, double y1, double x2, double y2, double *rho, double *theta, double *sigmarho, double *sigmatheta)
Definition: Hough2P.cxx:45
T atan2(T number)
Definition: d0nt_math.hpp:72