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...
IslandsTH2 *	fHspaceI
	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.

                               :
    fView(v),
    fRhoSz(rhosz),
    fThetaSz(thetasz),
    fXwinSz(rhosm),
    fYwinSz(thetasm),
    fSigma(3.0/sqrt(12.)),
    fRsqr(rsqr),
    fPco(pco),
    fHspaceW(0),
    fHspaceI(0),
    fH(v, 0, 0, 0, 0)
  { }

hough::Hough2P::~Hough2P

(

)

Definition at line 72 of file Hough2P.cxx.

References fHspaceI, and fHspaceW.

  {
    if (fHspaceW) { delete fHspaceW; fHspaceW = 0; }
    if (fHspaceI) { delete fHspaceI; fHspaceI = 0; }
  }

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

  {
    double r;   // Size of hough map in rho
    int    nr;  // Number of bins to use in rho
    int    nth; // Number of bins to use in theta
    
    // 0.55 gives a little buffer beyond the minimum required 0.5
    r  = 0.55*util::pythag((xhi-xlo),(yhi-ylo));
    nr = (int)rint(r/fRhoSz+0.5);
    
    nth = (int)rint(180.0/fThetaSz+0.5);
    
    const char* tiw=0;
    const char* tii=0;
    if      (fView==geo::kX) { tiw = "fHough2P_Hx"; tii="Hough2P_Ix"; }
    else if (fView==geo::kY) { tiw = "fHough2P_Hy"; tii="Hough2P_Iy"; }
    else abort();
    
    // Allocate the Hough space histogram
    if (fHspaceW!=0) { delete fHspaceW; fHspaceW=0; }
    fHspaceW = new TH2F(tiw,
                        ";rho [cm];theta [radians]",
                        nr, -r,  r,
                        nth,-0.5*M_PI/180.0,179.5*M_PI/180.0);

    // The islands histogram is closely related to the Hough space
    // histogram. Rebuild it now to ensure they stay in sync.
    if (fHspaceI!=0) { delete fHspaceI; fHspaceI=0; }
    fHspaceI = new IslandsTH2(tii, 
                              ";rho [cm];theta [radians]",
                              fHspaceW,
                              fPco,
                              false,
                              true);
  }

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

  {

    fHspaceI->LabelIslands(fHspaceW);

    unsigned int ipeak;
    int ix, iy;
    double peak  = 0.0;
    double r     = 0.0;
    double theta = 0.0;
    double p;
    for (ipeak=0; ipeak<fHspaceI->Nislands(); ++ipeak) {
      peak  = 0.0;
      r     = 0.0;
      theta = 0.0;
      for (ix=1; ix<=fHspaceI->GetNbinsX(); ++ix) {
        for (iy=1; iy<=fHspaceI->GetNbinsY(); ++iy) {
          if (fHspaceI->GetBinContent(ix,iy)==ipeak+1) {
            p = fHspaceW->GetBinContent(ix,iy);
            if (p>peak) {
              peak  = p;
              r     = fHspaceW->GetXaxis()->GetBinCenter(ix);
              theta = fHspaceW->GetYaxis()->GetBinCenter(iy);
            }
          } // if peak matches
        } // loop on y bins
      } // loop on x bins
      
      // NOTE:  The last three arguments below are threshold, sigma-rho, and sigma-theta for the Hough line
      //        which this Hough transform currently does NOT calculate (the island space calculates the
      //        threshold but it does NOT pass that value back to Hough2P.)
      if(peak > 0.0) fH.fPeak.push_back(rb::HoughPeak(r,theta,peak,0.0,0.0,0.0));
      
    } // loop on peaks
    
    //
    // Assert that there must always be at least one peak
    //
    if (fH.fPeak.size()==0) {
      peak  = 0.0;
      r     = 0.0;
      theta = 0.0;
      for (ix=1; ix<=fHspaceW->GetNbinsX(); ++ix) {
        for (iy=1; iy<=fHspaceW->GetNbinsY(); ++iy) {
          p = fHspaceW->GetBinContent(ix,iy);
          if (p>peak) {
            peak = p;
            r     = fHspaceW->GetXaxis()->GetBinCenter(ix);
            theta = fHspaceW->GetYaxis()->GetBinCenter(iy);
          }
        }
      }
      
      // NOTE:  The last three arguments below are threshold, sigma-rho, and sigma-theta for the Hough line
      //        which this Hough transform currently does NOT calculate (the island space calculates the
      //        threshold but it does NOT pass that value back to Hough2P.)
      if(peak > 0.0) fH.fPeak.push_back(rb::HoughPeak(r,theta,peak,0.0,0.0,0.0));
    }
    
    fH.SortPeaks();
  }

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

  {
    unsigned int i;
    unsigned int j;
    
    //
    // Find a box that completely contains all the hits. A note on
    // coordinates: Since the literature on Hough transforms uses an
    // "x-y" coodinate system I map NOvA z to Hough x and NOvA X/Y to
    // Hough y so that the formalizm matches what one would find in the
    // literature.
    //
    double tlo=9E9, thi=-9e9;
    double xlo=9e9, xhi=-9e9;
    double ylo=9e9, yhi=-9e9;
    for (i=0; i<h.size(); ++i) {
      if (h[i].fView == fView) {
        if (h[i].fZ<xlo)  xlo = h[i].fZ;
        if (h[i].fZ>xhi)  xhi = h[i].fZ;
        if (h[i].fXY<ylo) ylo = h[i].fXY;
        if (h[i].fXY>yhi) yhi = h[i].fXY;
        if (h[i].fT<tlo)  tlo = h[i].fT;
        if (h[i].fT>tlo)  thi = h[i].fT;
      }
    }
    fH.fTNSlo = tlo;
    fH.fTNShi = thi;
    fH.fZ0  = 0.5*(xlo+xhi);
    fH.fXY0 = 0.5*(ylo+yhi);
    
    //
    // Fill the hough maps
    //
    this->BuildMap(xlo,xhi,ylo,yhi);
    TAxis* xaxis = fHspaceW->GetXaxis();
    TAxis* yaxis = fHspaceW->GetYaxis();
    for (i=0; i<h.size(); ++i) {
      if (h[i].fView != fView) continue;
      for (j=i+1; j<h.size(); ++j) {
        if (h[j].fView != fView) continue;
        
        // Check distance between hit pair
        double rsqr = util::sqr(h[i].fZ-h[j].fZ) + util::sqr(h[i].fXY-h[j].fXY);
        if (rsqr>fRsqr) continue;
        
        // don't let points in the same column vote unless they are close together
        // if ( abs(h[i].fZ - h[j].fZ) < 6.0 && rsqr > fRsqr/4.0 ) continue;
        
        // don't let points in the same row vote unless they are far apart
        // if ( abs(h[i].fXY - h[j].fXY) < 6.0 && rsqr < fRsqr/4.0 ) continue;
        
        double rho,      theta;
        double sigmarho, sigmatheta;
        this->RhoTheta(h[i].fZ-fH.fZ0, h[i].fXY-fH.fXY0,
                       h[j].fZ-fH.fZ0, h[j].fXY-fH.fXY0,
                       &rho,           &theta,
                       &sigmarho,      &sigmatheta);
        //
        // Fill region weighted by errors
        //
        int ixlo = xaxis->FindBin(rho-3.0*sigmarho);
        int ixhi = xaxis->FindBin(rho+3.0*sigmarho);
        int iylo = yaxis->FindBin(theta-3.0*sigmatheta);
        int iyhi = yaxis->FindBin(theta+3.0*sigmatheta);
        
        if (ixlo<1)                 ixlo = 1;
        if (iylo<1)                 iylo = 1;
        if (ixhi>xaxis->GetNbins()) ixhi = xaxis->GetNbins();
        if (iyhi>yaxis->GetNbins()) iyhi = yaxis->GetNbins();
        
        int ix, iy;
        for (ix=ixlo; ix<=ixhi; ++ix) {
          double x = xaxis->GetBinCenter(ix);
          for (iy=iylo; iy<=iyhi; ++iy) {
            double y = yaxis->GetBinCenter(iy);
            double w = 2.0*h[i].fW*h[j].fW/(h[i].fW+h[j].fW);
            w *= 
              TMath::Gaus(x-rho,   0, sigmarho,   kTRUE) *
              TMath::Gaus(y-theta, 0, sigmatheta, kTRUE);
            fHspaceW->Fill(x,y,w);
          }
        }
      }
    }
    if (fXwinSz>0 || fYwinSz>0) this->SmoothMap();
    this->FindPeaks();

  }

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

  {
    *theta = atan2(x1-x2,y2-y1);
    *rho   = 0.5*(cos(*theta)*(x1+x2)+sin(*theta)*(y1+y2));
    //
    // Keep rho and theta within range [0,pi) (-R,+R)
    //
    while (*theta<0) {
      *theta +=  M_PI;
      *rho    = -(*rho);
    }
    while (*theta>=M_PI) {
      *theta -=  M_PI;
      *rho    = -(*rho);
    }
    
    // Compute the error estimates
    double d = util::pythag((x2-x1),(y2-y1));
    *sigmarho   = fSigma;
    *sigmatheta = M_SQRT2*fSigma/d;
  }

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

  {
    
    int i, j, m, n;
    int nx, ny;
    int ix, iy;
    double sum;
    int wx, wy;
    
    wx = 2*fXwinSz + 1;
    wy = 2*fYwinSz + 1;
    
    double w[wx][wy];
    
    // fill the smoothing weight map
    for(int a = 0; a < wx; a++) {
      for(int b = 0; b < wy; b++) {
        
        w[a][b] = TMath::Gaus(a - fXwinSz, 0, fXwinSz/3.0, kTRUE)*TMath::Gaus(b - fYwinSz, 0, fYwinSz/3.0, kTRUE);
        
      } // end for b
    } // end for a
    
    TH2F* htmp = new TH2F(*fHspaceW);
    fHspaceW->Reset();
    
    nx = htmp->GetNbinsX();
    ny = htmp->GetNbinsY();
    for (i=1; i<=nx; ++i) {
      for(j=1; j<=ny; ++j) {
        sum = 0.0;
        for (m=-fXwinSz; m<=fXwinSz; m++) {
          ix = i+m;
          // Don't wrap around in rho
          if (ix>nx) continue;
          if (ix<1)  continue;
          for (n=-fYwinSz; n<=fYwinSz; ++n) {
            iy = j+n;
            // Wrap in theta
            if (iy>ny) iy = iy%ny;
            if (iy<1)  iy = ny+iy;
            
            sum += w[m + fXwinSz][n + fYwinSz]*htmp->GetBinContent(ix,iy);
          }
        }
        fHspaceW->SetBinContent(i,j,sum);
      }
    }
    delete htmp;
    htmp = 0;
    
  }

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:

• /cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-23/HoughTransform/Hough2P.h
• /cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-23/HoughTransform/Hough2P.cxx