Hough2P.cxx

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
/// \brief   Perform a "2 point" Hough transform on a collection of hits
/// \authors messier@indiana.edu, mibaird@indiana.edu
/// \date
////////////////////////////////////////////////////////////////////////
#include "HoughTransform/Hough2P.h"

#include <algorithm>
#include <cstdlib>
#include <cmath>
#include <iostream>

#include "TH2F.h"
#include "TMath.h"

#include "HoughTransform/IslandsTH2.h"
#include "RecoBase/HitList.h"
#include "Utilities/func/MathUtil.h"

namespace hough
{ 
  //..................................................................... 
  Hough2P::Hough2P(geo::View_t v,
                   double rhosz,
                   double thetasz,
                   int rhosm,
                   int thetasm,
                   double rsqr,
                   double pco) :
    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)
  { }
  
  //......................................................................
  
  void Hough2P::RhoTheta(double x1, double y1,
                         double x2, double y2,
                         double* rho, double* theta,
                         double* sigmarho, double* sigmatheta)
  {
    *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;
  }
  
  //......................................................................
  
  Hough2P::~Hough2P() 
  {
    if (fHspaceW) { delete fHspaceW; fHspaceW = 0; }
    if (fHspaceI) { delete fHspaceI; fHspaceI = 0; }
  }
  
  //......................................................................
  
  void Hough2P::BuildMap(double xlo, double xhi,
                         double ylo, double yhi)
  {
    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 Hough2P::Map(const rb::HitList& h)
  {
    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 Hough2P::FindPeaks()
  {

    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 Hough2P::SmoothMap() 
  {
    
    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;
    
  }

} // end namespace hough

////////////////////////////////////////////////////////////////////////