LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
corner::CornerFinderAlg Class Reference

#include "CornerFinderAlg.h"

Public Member Functions

 CornerFinderAlg (fhicl::ParameterSet const &pset)
 
void GrabWires (std::vector< recob::Wire > const &wireVec, geo::Geometry const &)
 
void get_feature_points (std::vector< recob::EndPoint2D > &, geo::Geometry const &)
 
void get_feature_points_LineIntegralScore (std::vector< recob::EndPoint2D > &, geo::Geometry const &)
 
void get_feature_points_fast (std::vector< recob::EndPoint2D > &, geo::Geometry const &)
 
float line_integral (TH2F const &hist, int x1, float y1, int x2, float y2, float threshold) const
 
TH2F const & GetWireDataHist (unsigned int) const
 

Private Member Functions

void InitializeGeometry (geo::Geometry const &)
 
void create_image_histo (TH2F const &h_wire_data, TH2F &h_conversion) const
 
void create_derivative_histograms (TH2F const &h_conversion, TH2F &h_derivative_x, TH2F &h_derivative_y)
 
void create_cornerScore_histogram (TH2F const &h_derivative_x, TH2F const &h_derivative_y, TH2D &h_cornerScore)
 
std::vector< recob::EndPoint2Dperform_maximum_suppression (TH2D const &h_cornerScore, std::vector< geo::WireID > wireIDs, geo::View_t view, TH2D &h_maxSuppress, int startx=0, int starty=0) const
 
void calculate_line_integral_score (TH2F const &h_wire_data, std::vector< recob::EndPoint2D > const &corner_vector, std::vector< recob::EndPoint2D > &corner_lineIntegralScore_vector, TH2F &h_lineIntegralScore) const
 
void attach_feature_points (TH2F const &h_wire_data, std::vector< geo::WireID > const &wireIDs, geo::View_t view, std::vector< recob::EndPoint2D > &, int startx=0, int starty=0)
 
void attach_feature_points_LineIntegralScore (TH2F const &h_wire_data, std::vector< geo::WireID > const &wireIDs, geo::View_t view, std::vector< recob::EndPoint2D > &)
 
void create_smaller_histos (geo::Geometry const &)
 

Private Attributes

std::string fCalDataModuleLabel
 
std::string fConversion_algorithm
 
std::string fConversion_func
 
float fTrimming_threshold
 
int fTrimming_buffer
 
double fTrimming_totalThreshold
 
int fConversion_func_neighborhood
 
float fConversion_threshold
 
int fConversion_bins_per_input_x
 
int fConversion_bins_per_input_y
 
std::string fDerivative_method
 
int fDerivative_neighborhood
 
std::string fDerivative_BlurFunc
 
int fDerivative_BlurNeighborhood
 
int fCornerScore_neighborhood
 
std::string fCornerScore_algorithm
 
float fCornerScore_Noble_epsilon
 
float fCornerScore_Harris_kappa
 
int fMaxSuppress_neighborhood
 
int fMaxSuppress_threshold
 
float fIntegral_bin_threshold
 
float fIntegral_fraction_threshold
 
std::vector< TH2F > WireData_histos
 
std::vector< TH1D > WireData_histos_ProjectionX
 
std::vector< TH1D > WireData_histos_ProjectionY
 
std::vector< std::tuple< int, TH2F, int, int > > WireData_trimmed_histos
 
std::vector< std::vector< geo::WireID > > WireData_IDs
 
unsigned int event_number {}
 
unsigned int run_number {}
 

Detailed Description

Definition at line 33 of file CornerFinderAlg.h.

Constructor & Destructor Documentation

corner::CornerFinderAlg::CornerFinderAlg ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 40 of file CornerFinderAlg.cxx.

References fConversion_algorithm, fConversion_bins_per_input_x, fConversion_bins_per_input_y, fConversion_func, fConversion_func_neighborhood, fConversion_threshold, fCornerScore_algorithm, fCornerScore_Harris_kappa, fCornerScore_neighborhood, fCornerScore_Noble_epsilon, fDerivative_BlurFunc, fDerivative_BlurNeighborhood, fDerivative_method, fDerivative_neighborhood, fIntegral_bin_threshold, fIntegral_fraction_threshold, fMaxSuppress_neighborhood, fMaxSuppress_threshold, fTrimming_buffer, fTrimming_threshold, fTrimming_totalThreshold, and fhicl::ParameterSet::get().

41  : fCalDataModuleLabel{pset.get<std::string>("CalDataModuleLabel")}
42  , fConversion_algorithm{pset.get<std::string>("Conversion_algorithm")}
43  , fConversion_func{pset.get<std::string>("Conversion_function")}
44  , fTrimming_threshold{pset.get<float>("Trimming_threshold")}
45  , fTrimming_totalThreshold{pset.get<double>("Trimming_totalThreshold")}
46  , fConversion_func_neighborhood{pset.get<int>("Conversion_func_neighborhood")}
47  , fConversion_threshold{pset.get<float>("Conversion_threshold")}
48  , fConversion_bins_per_input_x{pset.get<int>("Conversion_bins_per_input_x")}
49  , fConversion_bins_per_input_y{pset.get<int>("Conversion_bins_per_input_y")}
50  , fDerivative_method{pset.get<std::string>("Derivative_method")}
51  , fDerivative_neighborhood{pset.get<int>("Derivative_neighborhood")}
52  , fDerivative_BlurFunc{pset.get<std::string>("Derivative_BlurFunc")}
53  , fDerivative_BlurNeighborhood{pset.get<int>("Derivative_BlurNeighborhood")}
54  , fCornerScore_neighborhood{pset.get<int>("CornerScore_neighborhood")}
55  , fCornerScore_algorithm{pset.get<std::string>("CornerScore_algorithm")}
56  , fCornerScore_Noble_epsilon{pset.get<float>("CornerScore_Noble_epsilon")}
57  , fCornerScore_Harris_kappa{pset.get<float>("CornerScore_Harris_kappa")}
58  , fMaxSuppress_neighborhood{pset.get<int>("MaxSuppress_neighborhood")}
59  , fMaxSuppress_threshold{pset.get<int>("MaxSuppress_threshold")}
60  , fIntegral_bin_threshold{pset.get<float>("Integral_bin_threshold")}
61  , fIntegral_fraction_threshold{pset.get<float>("Integral_fraction_threshold")}
62 {
68 }
std::string fDerivative_BlurFunc
std::string fCornerScore_algorithm
std::string fCalDataModuleLabel
std::string fConversion_algorithm
std::string fDerivative_method

Member Function Documentation

void corner::CornerFinderAlg::attach_feature_points ( TH2F const &  h_wire_data,
std::vector< geo::WireID > const &  wireIDs,
geo::View_t  view,
std::vector< recob::EndPoint2D > &  corner_vector,
int  startx = 0,
int  starty = 0 
)
private

Definition at line 441 of file CornerFinderAlg.cxx.

References create_cornerScore_histogram(), create_derivative_histograms(), create_image_histo(), event_number, fConversion_bins_per_input_x, fConversion_bins_per_input_y, perform_maximum_suppression(), run_number, x_max, and x_min.

Referenced by get_feature_points(), and get_feature_points_fast().

447 {
448 
449  const int x_bins = h_wire_data.GetNbinsX();
450  const float x_min = h_wire_data.GetXaxis()->GetBinLowEdge(1);
451  const float x_max = h_wire_data.GetXaxis()->GetBinUpEdge(x_bins);
452 
453  const int y_bins = h_wire_data.GetNbinsY();
454  const float y_min = h_wire_data.GetYaxis()->GetBinLowEdge(1);
455  const float y_max = h_wire_data.GetYaxis()->GetBinUpEdge(y_bins);
456 
457  const int converted_y_bins = y_bins / fConversion_bins_per_input_y;
458  const int converted_x_bins = x_bins / fConversion_bins_per_input_x;
459 
460  std::stringstream conversion_name;
461  conversion_name << "h_conversion_" << view << "_" << run_number << "_" << event_number;
462  std::stringstream dx_name;
463  dx_name << "h_derivative_x_" << view << "_" << run_number << "_" << event_number;
464  std::stringstream dy_name;
465  dy_name << "h_derivative_y_" << view << "_" << run_number << "_" << event_number;
466  std::stringstream cornerScore_name;
467  cornerScore_name << "h_cornerScore_" << view << "_" << run_number << "_" << event_number;
468  std::stringstream maxSuppress_name;
469  maxSuppress_name << "h_maxSuppress_" << view << "_" << run_number << "_" << event_number;
470 
471  TH2F conversion_histo(conversion_name.str().c_str(),
472  "Image Conversion Histogram",
473  converted_x_bins,
474  x_min,
475  x_max,
476  converted_y_bins,
477  y_min,
478  y_max);
479 
480  TH2F derivativeX_histo(dx_name.str().c_str(),
481  "Partial Derivatives (x)",
482  converted_x_bins,
483  x_min,
484  x_max,
485  converted_y_bins,
486  y_min,
487  y_max);
488 
489  TH2F derivativeY_histo(dy_name.str().c_str(),
490  "Partial Derivatives (y)",
491  converted_x_bins,
492  x_min,
493  x_max,
494  converted_y_bins,
495  y_min,
496  y_max);
497 
498  TH2D cornerScore_histo(cornerScore_name.str().c_str(),
499  "Corner Score",
500  converted_x_bins,
501  x_min,
502  x_max,
503  converted_y_bins,
504  y_min,
505  y_max);
506 
507  TH2D maxSuppress_histo(maxSuppress_name.str().c_str(),
508  "Corner Points (Maximum Suppressed)",
509  converted_x_bins,
510  x_min,
511  x_max,
512  converted_y_bins,
513  y_min,
514  y_max);
515 
516  create_image_histo(h_wire_data, conversion_histo);
517  create_derivative_histograms(conversion_histo, derivativeX_histo, derivativeY_histo);
518  create_cornerScore_histogram(derivativeX_histo, derivativeY_histo, cornerScore_histo);
519  corner_vector = perform_maximum_suppression(
520  cornerScore_histo, wireIDs, view, maxSuppress_histo, startx, starty);
521 }
void create_image_histo(TH2F const &h_wire_data, TH2F &h_conversion) const
void create_cornerScore_histogram(TH2F const &h_derivative_x, TH2F const &h_derivative_y, TH2D &h_cornerScore)
std::vector< recob::EndPoint2D > perform_maximum_suppression(TH2D const &h_cornerScore, std::vector< geo::WireID > wireIDs, geo::View_t view, TH2D &h_maxSuppress, int startx=0, int starty=0) const
void create_derivative_histograms(TH2F const &h_conversion, TH2F &h_derivative_x, TH2F &h_derivative_y)
double x_min
Definition: berger.C:15
double x_max
Definition: berger.C:16
void corner::CornerFinderAlg::attach_feature_points_LineIntegralScore ( TH2F const &  h_wire_data,
std::vector< geo::WireID > const &  wireIDs,
geo::View_t  view,
std::vector< recob::EndPoint2D > &  corner_vector 
)
private

Definition at line 525 of file CornerFinderAlg.cxx.

References calculate_line_integral_score(), create_cornerScore_histogram(), create_derivative_histograms(), create_image_histo(), event_number, fConversion_bins_per_input_x, fConversion_bins_per_input_y, perform_maximum_suppression(), run_number, x_max, and x_min.

Referenced by get_feature_points_LineIntegralScore().

530 {
531  const int x_bins = h_wire_data.GetNbinsX();
532  const float x_min = h_wire_data.GetXaxis()->GetBinLowEdge(1);
533  const float x_max = h_wire_data.GetXaxis()->GetBinUpEdge(x_bins);
534 
535  const int y_bins = h_wire_data.GetNbinsY();
536  const float y_min = h_wire_data.GetYaxis()->GetBinLowEdge(1);
537  const float y_max = h_wire_data.GetYaxis()->GetBinUpEdge(y_bins);
538 
539  const int converted_y_bins = y_bins / fConversion_bins_per_input_y;
540  const int converted_x_bins = x_bins / fConversion_bins_per_input_x;
541 
542  std::stringstream conversion_name;
543  conversion_name << "h_conversion_" << view << "_" << run_number << "_" << event_number;
544  std::stringstream dx_name;
545  dx_name << "h_derivative_x_" << view << "_" << run_number << "_" << event_number;
546  std::stringstream dy_name;
547  dy_name << "h_derivative_y_" << view << "_" << run_number << "_" << event_number;
548  std::stringstream cornerScore_name;
549  cornerScore_name << "h_cornerScore_" << view << "_" << run_number << "_" << event_number;
550  std::stringstream maxSuppress_name;
551  maxSuppress_name << "h_maxSuppress_" << view << "_" << run_number << "_" << event_number;
552 
553  TH2F h_conversion(conversion_name.str().c_str(),
554  "Image Conversion Histogram",
555  converted_x_bins,
556  x_min,
557  x_max,
558  converted_y_bins,
559  y_min,
560  y_max);
561  TH2F h_derivative_x(dx_name.str().c_str(),
562  "Partial Derivatives (x)",
563  converted_x_bins,
564  x_min,
565  x_max,
566  converted_y_bins,
567  y_min,
568  y_max);
569  TH2F h_derivative_y(dy_name.str().c_str(),
570  "Partial Derivatives (y)",
571  converted_x_bins,
572  x_min,
573  x_max,
574  converted_y_bins,
575  y_min,
576  y_max);
577  TH2D h_cornerScore(cornerScore_name.str().c_str(),
578  "Feature Point Corner Score",
579  converted_x_bins,
580  x_min,
581  x_max,
582  converted_y_bins,
583  y_min,
584  y_max);
585  TH2D h_maxSuppress(maxSuppress_name.str().c_str(),
586  "Corner Points (Maximum Suppressed)",
587  converted_x_bins,
588  x_min,
589  x_max,
590  converted_y_bins,
591  y_min,
592  y_max);
593 
594  create_image_histo(h_wire_data, h_conversion);
595  create_derivative_histograms(h_conversion, h_derivative_x, h_derivative_y);
596  create_cornerScore_histogram(h_derivative_x, h_derivative_y, h_cornerScore);
597 
598  auto corner_vector_tmp = perform_maximum_suppression(h_cornerScore, wireIDs, view, h_maxSuppress);
599 
600  std::stringstream LI_name;
601  LI_name << "h_lineIntegralScore_" << view << "_" << run_number << "_" << event_number;
602  TH2F h_lineIntegralScore(
603  LI_name.str().c_str(), "Line Integral Score", x_bins, x_min, x_max, y_bins, y_min, y_max);
604  calculate_line_integral_score(h_wire_data, corner_vector_tmp, corner_vector, h_lineIntegralScore);
605 }
void create_image_histo(TH2F const &h_wire_data, TH2F &h_conversion) const
void create_cornerScore_histogram(TH2F const &h_derivative_x, TH2F const &h_derivative_y, TH2D &h_cornerScore)
std::vector< recob::EndPoint2D > perform_maximum_suppression(TH2D const &h_cornerScore, std::vector< geo::WireID > wireIDs, geo::View_t view, TH2D &h_maxSuppress, int startx=0, int starty=0) const
void create_derivative_histograms(TH2F const &h_conversion, TH2F &h_derivative_x, TH2F &h_derivative_y)
void calculate_line_integral_score(TH2F const &h_wire_data, std::vector< recob::EndPoint2D > const &corner_vector, std::vector< recob::EndPoint2D > &corner_lineIntegralScore_vector, TH2F &h_lineIntegralScore) const
double x_min
Definition: berger.C:15
double x_max
Definition: berger.C:16
void corner::CornerFinderAlg::calculate_line_integral_score ( TH2F const &  h_wire_data,
std::vector< recob::EndPoint2D > const &  corner_vector,
std::vector< recob::EndPoint2D > &  corner_lineIntegralScore_vector,
TH2F &  h_lineIntegralScore 
) const
private

Definition at line 1145 of file CornerFinderAlg.cxx.

References fIntegral_bin_threshold, fIntegral_fraction_threshold, and line_integral().

Referenced by attach_feature_points_LineIntegralScore().

1150 {
1151 
1152  for (auto const i_corner : corner_vector) {
1153 
1154  float score = 0;
1155 
1156  for (auto const j_corner : corner_vector) {
1157 
1158  if (line_integral(h_wire_data,
1159  i_corner.WireID().Wire,
1160  i_corner.DriftTime(),
1161  j_corner.WireID().Wire,
1162  j_corner.DriftTime(),
1164  score += 1.;
1165  }
1166  }
1167 
1168  recob::EndPoint2D corner(i_corner.DriftTime(),
1169  i_corner.WireID(),
1170  score,
1171  i_corner.ID(),
1172  i_corner.View(),
1173  i_corner.Charge());
1174 
1175  corner_lineIntegralScore_vector.push_back(corner);
1176 
1177  h_lineIntegralScore.SetBinContent(h_wire_data.GetXaxis()->FindBin(i_corner.WireID().Wire),
1178  h_wire_data.GetYaxis()->FindBin(i_corner.DriftTime()),
1179  score);
1180  }
1181 }
float line_integral(TH2F const &hist, int x1, float y1, int x2, float y2, float threshold) const
void corner::CornerFinderAlg::create_cornerScore_histogram ( TH2F const &  h_derivative_x,
TH2F const &  h_derivative_y,
TH2D &  h_cornerScore 
)
private

Definition at line 959 of file CornerFinderAlg.cxx.

References fCornerScore_algorithm, fCornerScore_Harris_kappa, fCornerScore_neighborhood, and fCornerScore_Noble_epsilon.

Referenced by attach_feature_points(), and attach_feature_points_LineIntegralScore().

962 {
963 
964  const int x_bins = h_derivative_x.GetNbinsX();
965  const int y_bins = h_derivative_y.GetNbinsY();
966 
967  //the structure tensor elements
968  double st_xx = 0., st_xy = 0., st_yy = 0.;
969 
970  for (int iy = 1 + fCornerScore_neighborhood; iy <= (y_bins - fCornerScore_neighborhood); iy++) {
971  for (int ix = 1 + fCornerScore_neighborhood; ix <= (x_bins - fCornerScore_neighborhood); ix++) {
972 
973  if (ix == 1 + fCornerScore_neighborhood) {
974  st_xx = 0.;
975  st_xy = 0.;
976  st_yy = 0.;
977 
978  for (int jx = ix - fCornerScore_neighborhood; jx <= ix + fCornerScore_neighborhood; jx++) {
979  for (int jy = iy - fCornerScore_neighborhood; jy <= iy + fCornerScore_neighborhood;
980  jy++) {
981 
982  st_xx += h_derivative_x.GetBinContent(jx, jy) * h_derivative_x.GetBinContent(jx, jy);
983  st_yy += h_derivative_y.GetBinContent(jx, jy) * h_derivative_y.GetBinContent(jx, jy);
984  st_xy += h_derivative_x.GetBinContent(jx, jy) * h_derivative_y.GetBinContent(jx, jy);
985  }
986  }
987  }
988 
989  // we do it this way to reduce computation time
990  else {
991  for (int jy = iy - fCornerScore_neighborhood; jy <= iy + fCornerScore_neighborhood; jy++) {
992 
993  st_xx -= h_derivative_x.GetBinContent(ix - fCornerScore_neighborhood - 1, jy) *
994  h_derivative_x.GetBinContent(ix - fCornerScore_neighborhood - 1, jy);
995  st_xx += h_derivative_x.GetBinContent(ix + fCornerScore_neighborhood, jy) *
996  h_derivative_x.GetBinContent(ix + fCornerScore_neighborhood, jy);
997 
998  st_yy -= h_derivative_y.GetBinContent(ix - fCornerScore_neighborhood - 1, jy) *
999  h_derivative_y.GetBinContent(ix - fCornerScore_neighborhood - 1, jy);
1000  st_yy += h_derivative_y.GetBinContent(ix + fCornerScore_neighborhood, jy) *
1001  h_derivative_y.GetBinContent(ix + fCornerScore_neighborhood, jy);
1002 
1003  st_xy -= h_derivative_x.GetBinContent(ix - fCornerScore_neighborhood - 1, jy) *
1004  h_derivative_y.GetBinContent(ix - fCornerScore_neighborhood - 1, jy);
1005  st_xy += h_derivative_x.GetBinContent(ix + fCornerScore_neighborhood, jy) *
1006  h_derivative_y.GetBinContent(ix + fCornerScore_neighborhood, jy);
1007  }
1008  }
1009 
1010  if (fCornerScore_algorithm.compare("Noble") == 0) {
1011  h_cornerScore.SetBinContent(
1012  ix, iy, (st_xx * st_yy - st_xy * st_xy) / (st_xx + st_yy + fCornerScore_Noble_epsilon));
1013  }
1014  else if (fCornerScore_algorithm.compare("Harris") == 0) {
1015  h_cornerScore.SetBinContent(
1016  ix,
1017  iy,
1018  (st_xx * st_yy - st_xy * st_xy) -
1019  ((st_xx + st_yy) * (st_xx + st_yy) * fCornerScore_Harris_kappa));
1020  }
1021  else {
1022  mf::LogError("CornerFinderAlg") << "BAD CORNER ALGORITHM: " << fCornerScore_algorithm;
1023  return;
1024  }
1025 
1026  } // end for loop over x bins
1027  } // end for loop over y bins
1028 }
std::string fCornerScore_algorithm
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
void corner::CornerFinderAlg::create_derivative_histograms ( TH2F const &  h_conversion,
TH2F &  h_derivative_x,
TH2F &  h_derivative_y 
)
private

Definition at line 677 of file CornerFinderAlg.cxx.

References fDerivative_BlurNeighborhood, fDerivative_method, and fDerivative_neighborhood.

Referenced by attach_feature_points(), and attach_feature_points_LineIntegralScore().

680 {
681 
682  const int x_bins = h_conversion.GetNbinsX();
683  const int y_bins = h_conversion.GetNbinsY();
684 
685  for (int iy = 1 + fDerivative_neighborhood; iy <= (y_bins - fDerivative_neighborhood); iy++) {
686  for (int ix = 1 + fDerivative_neighborhood; ix <= (x_bins - fDerivative_neighborhood); ix++) {
687 
688  if (fDerivative_method.compare("Sobel") == 0) {
689 
690  if (fDerivative_neighborhood == 1) {
691  h_derivative_x.SetBinContent(ix,
692  iy,
693  0.5 * (h_conversion.GetBinContent(ix + 1, iy) -
694  h_conversion.GetBinContent(ix - 1, iy)) +
695  0.25 * (h_conversion.GetBinContent(ix + 1, iy + 1) -
696  h_conversion.GetBinContent(ix - 1, iy + 1)) +
697  0.25 * (h_conversion.GetBinContent(ix + 1, iy - 1) -
698  h_conversion.GetBinContent(ix - 1, iy - 1)));
699  h_derivative_y.SetBinContent(ix,
700  iy,
701  0.5 * (h_conversion.GetBinContent(ix, iy + 1) -
702  h_conversion.GetBinContent(ix, iy - 1)) +
703  0.25 * (h_conversion.GetBinContent(ix - 1, iy + 1) -
704  h_conversion.GetBinContent(ix - 1, iy - 1)) +
705  0.25 * (h_conversion.GetBinContent(ix + 1, iy + 1) -
706  h_conversion.GetBinContent(ix + 1, iy - 1)));
707  }
708  else if (fDerivative_neighborhood == 2) {
709  h_derivative_x.SetBinContent(
710  ix,
711  iy,
712  12 * (h_conversion.GetBinContent(ix + 1, iy) - h_conversion.GetBinContent(ix - 1, iy)) +
713  8 * (h_conversion.GetBinContent(ix + 1, iy + 1) -
714  h_conversion.GetBinContent(ix - 1, iy + 1)) +
715  8 * (h_conversion.GetBinContent(ix + 1, iy - 1) -
716  h_conversion.GetBinContent(ix - 1, iy - 1)) +
717  2 * (h_conversion.GetBinContent(ix + 1, iy + 2) -
718  h_conversion.GetBinContent(ix - 1, iy + 2)) +
719  2 * (h_conversion.GetBinContent(ix + 1, iy - 2) -
720  h_conversion.GetBinContent(ix - 1, iy - 2)) +
721  6 *
722  (h_conversion.GetBinContent(ix + 2, iy) - h_conversion.GetBinContent(ix - 2, iy)) +
723  4 * (h_conversion.GetBinContent(ix + 2, iy + 1) -
724  h_conversion.GetBinContent(ix - 2, iy + 1)) +
725  4 * (h_conversion.GetBinContent(ix + 2, iy - 1) -
726  h_conversion.GetBinContent(ix - 2, iy - 1)) +
727  1 * (h_conversion.GetBinContent(ix + 2, iy + 2) -
728  h_conversion.GetBinContent(ix - 2, iy + 2)) +
729  1 * (h_conversion.GetBinContent(ix + 2, iy - 2) -
730  h_conversion.GetBinContent(ix - 2, iy - 2)));
731  h_derivative_y.SetBinContent(
732  ix,
733  iy,
734  12 * (h_conversion.GetBinContent(ix, iy + 1) - h_conversion.GetBinContent(ix, iy - 1)) +
735  8 * (h_conversion.GetBinContent(ix - 1, iy + 1) -
736  h_conversion.GetBinContent(ix - 1, iy - 1)) +
737  8 * (h_conversion.GetBinContent(ix + 1, iy + 1) -
738  h_conversion.GetBinContent(ix + 1, iy - 1)) +
739  2 * (h_conversion.GetBinContent(ix - 2, iy + 1) -
740  h_conversion.GetBinContent(ix - 2, iy - 1)) +
741  2 * (h_conversion.GetBinContent(ix + 2, iy + 1) -
742  h_conversion.GetBinContent(ix + 2, iy - 1)) +
743  6 *
744  (h_conversion.GetBinContent(ix, iy + 2) - h_conversion.GetBinContent(ix, iy - 2)) +
745  4 * (h_conversion.GetBinContent(ix - 1, iy + 2) -
746  h_conversion.GetBinContent(ix - 1, iy - 2)) +
747  4 * (h_conversion.GetBinContent(ix + 1, iy + 2) -
748  h_conversion.GetBinContent(ix + 1, iy - 2)) +
749  1 * (h_conversion.GetBinContent(ix - 2, iy + 2) -
750  h_conversion.GetBinContent(ix - 2, iy - 2)) +
751  1 * (h_conversion.GetBinContent(ix + 2, iy + 2) -
752  h_conversion.GetBinContent(ix + 2, iy - 2)));
753  }
754  else {
755  mf::LogError("CornerFinderAlg")
756  << "Sobel derivative not supported for neighborhoods > 2.";
757  return;
758  }
759 
760  } //end if Sobel
761 
762  else if (fDerivative_method.compare("local") == 0) {
763 
764  if (fDerivative_neighborhood == 1) {
765  h_derivative_x.SetBinContent(
766  ix,
767  iy,
768  (h_conversion.GetBinContent(ix + 1, iy) - h_conversion.GetBinContent(ix - 1, iy)));
769  h_derivative_y.SetBinContent(
770  ix,
771  iy,
772  (h_conversion.GetBinContent(ix, iy + 1) - h_conversion.GetBinContent(ix, iy - 1)));
773  }
774  else {
775  mf::LogError("CornerFinderAlg")
776  << "Local derivative not yet supported for neighborhoods > 1.";
777  return;
778  }
779  } //end if local
780 
781  else {
782  mf::LogError("CornerFinderAlg") << "Bad derivative algorithm! " << fDerivative_method;
783  return;
784  }
785  }
786  }
787 
788  //this is just a double Gaussian
789  float func_blur[11][11];
790  func_blur[0][0] = 0.000000;
791  func_blur[0][1] = 0.000000;
792  func_blur[0][2] = 0.000000;
793  func_blur[0][3] = 0.000001;
794  func_blur[0][4] = 0.000002;
795  func_blur[0][5] = 0.000004;
796  func_blur[0][6] = 0.000002;
797  func_blur[0][7] = 0.000001;
798  func_blur[0][8] = 0.000000;
799  func_blur[0][9] = 0.000000;
800  func_blur[0][10] = 0.000000;
801  func_blur[1][0] = 0.000000;
802  func_blur[1][1] = 0.000000;
803  func_blur[1][2] = 0.000004;
804  func_blur[1][3] = 0.000045;
805  func_blur[1][4] = 0.000203;
806  func_blur[1][5] = 0.000335;
807  func_blur[1][6] = 0.000203;
808  func_blur[1][7] = 0.000045;
809  func_blur[1][8] = 0.000004;
810  func_blur[1][9] = 0.000000;
811  func_blur[1][10] = 0.000000;
812  func_blur[2][0] = 0.000000;
813  func_blur[2][1] = 0.000004;
814  func_blur[2][2] = 0.000123;
815  func_blur[2][3] = 0.001503;
816  func_blur[2][4] = 0.006738;
817  func_blur[2][5] = 0.011109;
818  func_blur[2][6] = 0.006738;
819  func_blur[2][7] = 0.001503;
820  func_blur[2][8] = 0.000123;
821  func_blur[2][9] = 0.000004;
822  func_blur[2][10] = 0.000000;
823  func_blur[3][0] = 0.000001;
824  func_blur[3][1] = 0.000045;
825  func_blur[3][2] = 0.001503;
826  func_blur[3][3] = 0.018316;
827  func_blur[3][4] = 0.082085;
828  func_blur[3][5] = 0.135335;
829  func_blur[3][6] = 0.082085;
830  func_blur[3][7] = 0.018316;
831  func_blur[3][8] = 0.001503;
832  func_blur[3][9] = 0.000045;
833  func_blur[3][10] = 0.000001;
834  func_blur[4][0] = 0.000002;
835  func_blur[4][1] = 0.000203;
836  func_blur[4][2] = 0.006738;
837  func_blur[4][3] = 0.082085;
838  func_blur[4][4] = 0.367879;
839  func_blur[4][5] = 0.606531;
840  func_blur[4][6] = 0.367879;
841  func_blur[4][7] = 0.082085;
842  func_blur[4][8] = 0.006738;
843  func_blur[4][9] = 0.000203;
844  func_blur[4][10] = 0.000002;
845  func_blur[5][0] = 0.000004;
846  func_blur[5][1] = 0.000335;
847  func_blur[5][2] = 0.011109;
848  func_blur[5][3] = 0.135335;
849  func_blur[5][4] = 0.606531;
850  func_blur[5][5] = 1.000000;
851  func_blur[5][6] = 0.606531;
852  func_blur[5][7] = 0.135335;
853  func_blur[5][8] = 0.011109;
854  func_blur[5][9] = 0.000335;
855  func_blur[5][10] = 0.000004;
856  func_blur[6][0] = 0.000002;
857  func_blur[6][1] = 0.000203;
858  func_blur[6][2] = 0.006738;
859  func_blur[6][3] = 0.082085;
860  func_blur[6][4] = 0.367879;
861  func_blur[6][5] = 0.606531;
862  func_blur[6][6] = 0.367879;
863  func_blur[6][7] = 0.082085;
864  func_blur[6][8] = 0.006738;
865  func_blur[6][9] = 0.000203;
866  func_blur[6][10] = 0.000002;
867  func_blur[7][0] = 0.000001;
868  func_blur[7][1] = 0.000045;
869  func_blur[7][2] = 0.001503;
870  func_blur[7][3] = 0.018316;
871  func_blur[7][4] = 0.082085;
872  func_blur[7][5] = 0.135335;
873  func_blur[7][6] = 0.082085;
874  func_blur[7][7] = 0.018316;
875  func_blur[7][8] = 0.001503;
876  func_blur[7][9] = 0.000045;
877  func_blur[7][10] = 0.000001;
878  func_blur[8][0] = 0.000000;
879  func_blur[8][1] = 0.000004;
880  func_blur[8][2] = 0.000123;
881  func_blur[8][3] = 0.001503;
882  func_blur[8][4] = 0.006738;
883  func_blur[8][5] = 0.011109;
884  func_blur[8][6] = 0.006738;
885  func_blur[8][7] = 0.001503;
886  func_blur[8][8] = 0.000123;
887  func_blur[8][9] = 0.000004;
888  func_blur[8][10] = 0.000000;
889  func_blur[9][0] = 0.000000;
890  func_blur[9][1] = 0.000000;
891  func_blur[9][2] = 0.000004;
892  func_blur[9][3] = 0.000045;
893  func_blur[9][4] = 0.000203;
894  func_blur[9][5] = 0.000335;
895  func_blur[9][6] = 0.000203;
896  func_blur[9][7] = 0.000045;
897  func_blur[9][8] = 0.000004;
898  func_blur[9][9] = 0.000000;
899  func_blur[9][10] = 0.000000;
900  func_blur[10][0] = 0.000000;
901  func_blur[10][1] = 0.000000;
902  func_blur[10][2] = 0.000000;
903  func_blur[10][3] = 0.000001;
904  func_blur[10][4] = 0.000002;
905  func_blur[10][5] = 0.000004;
906  func_blur[10][6] = 0.000002;
907  func_blur[10][7] = 0.000001;
908  func_blur[10][8] = 0.000000;
909  func_blur[10][9] = 0.000000;
910  func_blur[10][10] = 0.000000;
911 
912  double temp_integral_x = 0;
913  double temp_integral_y = 0;
914 
916 
917  if (fDerivative_BlurNeighborhood > 10) {
918  mf::LogWarning("CornerFinderAlg")
919  << "WARNING...BlurNeighborhoods>10 not currently allowed. Shrinking to 10.";
921  }
922 
923  TH2F* h_clone_derivative_x = (TH2F*)h_derivative_x.Clone("h_clone_derivative_x");
924  TH2F* h_clone_derivative_y = (TH2F*)h_derivative_y.Clone("h_clone_derivative_y");
925 
926  temp_integral_x = 0;
927  temp_integral_y = 0;
928 
929  for (int ix = 1; ix <= h_derivative_x.GetNbinsX(); ix++) {
930  for (int iy = 1; iy <= h_derivative_y.GetNbinsY(); iy++) {
931 
932  temp_integral_x = 0;
933  temp_integral_y = 0;
934 
935  for (int jx = ix - fDerivative_BlurNeighborhood; jx <= ix + fDerivative_BlurNeighborhood;
936  jx++) {
937  for (int jy = iy - fDerivative_BlurNeighborhood; jy <= iy + fDerivative_BlurNeighborhood;
938  jy++) {
939  temp_integral_x +=
940  h_clone_derivative_x->GetBinContent(jx, jy) * func_blur[(ix - jx) + 5][(iy - jy) + 5];
941  temp_integral_y +=
942  h_clone_derivative_y->GetBinContent(jx, jy) * func_blur[(ix - jx) + 5][(iy - jy) + 5];
943  }
944  }
945  h_derivative_x.SetBinContent(ix, iy, temp_integral_x);
946  h_derivative_y.SetBinContent(ix, iy, temp_integral_y);
947  }
948  }
949 
950  delete h_clone_derivative_x;
951  delete h_clone_derivative_y;
952 
953  } //end if blur
954 }
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::string fDerivative_method
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void corner::CornerFinderAlg::create_image_histo ( TH2F const &  h_wire_data,
TH2F &  h_conversion 
) const
private

Definition at line 609 of file CornerFinderAlg.cxx.

References fConversion_algorithm, fConversion_func, fConversion_func_neighborhood, and fConversion_threshold.

Referenced by attach_feature_points(), and attach_feature_points_LineIntegralScore().

610 {
611 
612  double temp_integral = 0;
613 
614  const TF2 fConversion_TF2("fConversion_func", fConversion_func.c_str(), -20, 20, -20, 20);
615 
616  for (int ix = 1; ix <= h_conversion.GetNbinsX(); ix++) {
617  for (int iy = 1; iy <= h_conversion.GetNbinsY(); iy++) {
618 
619  temp_integral = h_wire_data.Integral(ix, ix, iy, iy);
620 
621  if (temp_integral > fConversion_threshold) {
622 
623  if (fConversion_algorithm.compare("binary") == 0)
624  h_conversion.SetBinContent(ix, iy, 10 * fConversion_threshold);
625  else if (fConversion_algorithm.compare("standard") == 0)
626  h_conversion.SetBinContent(ix, iy, temp_integral);
627 
628  else if (fConversion_algorithm.compare("function") == 0) {
629 
630  temp_integral = 0;
631  for (int jx = ix - fConversion_func_neighborhood;
633  jx++) {
634  for (int jy = iy - fConversion_func_neighborhood;
636  jy++) {
637  temp_integral +=
638  h_wire_data.GetBinContent(jx, jy) * fConversion_TF2.Eval(ix - jx, iy - jy);
639  }
640  }
641  h_conversion.SetBinContent(ix, iy, temp_integral);
642  }
643 
644  else if (fConversion_algorithm.compare("skeleton") == 0) {
645 
646  if ((temp_integral > h_wire_data.GetBinContent(ix - 1, iy) &&
647  temp_integral > h_wire_data.GetBinContent(ix + 1, iy)) ||
648  (temp_integral > h_wire_data.GetBinContent(ix, iy - 1) &&
649  temp_integral > h_wire_data.GetBinContent(ix, iy + 1)))
650  h_conversion.SetBinContent(ix, iy, temp_integral);
651  else
652  h_conversion.SetBinContent(ix, iy, fConversion_threshold);
653  }
654  else if (fConversion_algorithm.compare("sk_bin") == 0) {
655 
656  if ((temp_integral > h_wire_data.GetBinContent(ix - 1, iy) &&
657  temp_integral > h_wire_data.GetBinContent(ix + 1, iy)) ||
658  (temp_integral > h_wire_data.GetBinContent(ix, iy - 1) &&
659  temp_integral > h_wire_data.GetBinContent(ix, iy + 1)))
660  h_conversion.SetBinContent(ix, iy, 10 * fConversion_threshold);
661  else
662  h_conversion.SetBinContent(ix, iy, fConversion_threshold);
663  }
664  else
665  h_conversion.SetBinContent(ix, iy, temp_integral);
666  }
667 
668  else
669  h_conversion.SetBinContent(ix, iy, fConversion_threshold);
670  }
671  }
672 }
std::string fConversion_algorithm
void corner::CornerFinderAlg::create_smaller_histos ( geo::Geometry const &  my_geometry)
private

Definition at line 260 of file CornerFinderAlg.cxx.

References fTrimming_buffer, fTrimming_threshold, fTrimming_totalThreshold, geo::GeometryCore::Iterate(), MF_LOG_DEBUG, WireData_histos, WireData_histos_ProjectionX, WireData_histos_ProjectionY, and WireData_trimmed_histos.

Referenced by get_feature_points_fast().

261 {
262 
263  for (auto const& pid : my_geometry.Iterate<geo::PlaneID>()) {
264 
265  MF_LOG_DEBUG("CornerFinderAlg") << "Working plane " << pid.Plane << ".";
266 
267  int x_bins = WireData_histos_ProjectionX.at(pid.Plane).GetNbinsX();
268  int y_bins = WireData_histos_ProjectionY.at(pid.Plane).GetNbinsX();
269 
270  std::vector<int> cut_points_x{0};
271  std::vector<int> cut_points_y{0};
272 
273  for (int ix = 1; ix <= x_bins; ix++) {
274 
275  float this_value = WireData_histos_ProjectionX.at(pid.Plane).GetBinContent(ix);
276 
277  if (ix < fTrimming_buffer || ix > (x_bins - fTrimming_buffer)) continue;
278 
279  int jx = ix - fTrimming_buffer;
280  while (this_value < fTrimming_threshold) {
281  if (jx == ix + fTrimming_buffer) break;
282  this_value = WireData_histos_ProjectionX.at(pid.Plane).GetBinContent(jx);
283  jx++;
284  }
285  if (this_value < fTrimming_threshold) { cut_points_x.push_back(ix); }
286  }
287 
288  for (int iy = 1; iy <= y_bins; iy++) {
289 
290  float this_value = WireData_histos_ProjectionY.at(pid.Plane).GetBinContent(iy);
291 
292  if (iy < fTrimming_buffer || iy > (y_bins - fTrimming_buffer)) continue;
293 
294  int jy = iy - fTrimming_buffer;
295  while (this_value < fTrimming_threshold) {
296  if (jy == iy + fTrimming_buffer) break;
297  this_value = WireData_histos_ProjectionY.at(pid.Plane).GetBinContent(jy);
298  jy++;
299  }
300  if (this_value < fTrimming_threshold) { cut_points_y.push_back(iy); }
301  }
302 
303  MF_LOG_DEBUG("CornerFinderAlg")
304  << "We have a total of " << cut_points_x.size() << " x cut points."
305  << "\nWe have a total of " << cut_points_y.size() << " y cut points.";
306 
307  std::vector<int> x_low{1};
308  std::vector<int> x_high{x_bins};
309  std::vector<int> y_low{1};
310  std::vector<int> y_high{y_bins};
311  bool x_change = true;
312  bool y_change = true;
313  while (x_change || y_change) {
314 
315  x_change = false;
316  y_change = false;
317 
318  size_t current_size = x_low.size();
319 
320  for (size_t il = 0; il < current_size; il++) {
321 
322  int comp_value = (x_high.at(il) + x_low.at(il)) / 2;
323  std::sort(cut_points_x.begin(), cut_points_x.end(), compare_to_value(comp_value));
324 
325  if (cut_points_x.at(0) <= x_low.at(il) || cut_points_x.at(0) >= x_high.at(il)) continue;
326 
327  double integral_low = WireData_histos.at(pid.Plane).Integral(
328  x_low.at(il), cut_points_x.at(0), y_low.at(il), y_high.at(il));
329  double integral_high = WireData_histos.at(pid.Plane).Integral(
330  cut_points_x.at(0), x_high.at(il), y_low.at(il), y_high.at(il));
331  if (integral_low > fTrimming_totalThreshold && integral_high > fTrimming_totalThreshold) {
332  x_low.push_back(cut_points_x.at(0));
333  x_high.push_back(x_high.at(il));
334  y_low.push_back(y_low.at(il));
335  y_high.push_back(y_high.at(il));
336 
337  x_high[il] = cut_points_x.at(0);
338  x_change = true;
339  }
340  else if (integral_low > fTrimming_totalThreshold &&
341  integral_high < fTrimming_totalThreshold) {
342  x_high[il] = cut_points_x.at(0);
343  x_change = true;
344  }
345  else if (integral_low < fTrimming_totalThreshold &&
346  integral_high > fTrimming_totalThreshold) {
347  x_low[il] = cut_points_x.at(0);
348  x_change = true;
349  }
350  }
351 
352  current_size = x_low.size();
353 
354  for (size_t il = 0; il < current_size; il++) {
355 
356  int comp_value = (y_high.at(il) - y_low.at(il)) / 2;
357  std::sort(cut_points_y.begin(), cut_points_y.end(), compare_to_value(comp_value));
358 
359  if (cut_points_y.at(0) <= y_low.at(il) || cut_points_y.at(0) >= y_high.at(il)) continue;
360 
361  double integral_low = WireData_histos.at(pid.Plane).Integral(
362  x_low.at(il), x_high.at(il), y_low.at(il), cut_points_y.at(0));
363  double integral_high = WireData_histos.at(pid.Plane).Integral(
364  x_low.at(il), x_high.at(il), cut_points_y.at(0), y_high.at(il));
365  if (integral_low > fTrimming_totalThreshold && integral_high > fTrimming_totalThreshold) {
366  y_low.push_back(cut_points_y.at(0));
367  y_high.push_back(y_high.at(il));
368  x_low.push_back(x_low.at(il));
369  x_high.push_back(x_high.at(il));
370 
371  y_high[il] = cut_points_y.at(0);
372  y_change = true;
373  }
374  else if (integral_low > fTrimming_totalThreshold &&
375  integral_high < fTrimming_totalThreshold) {
376  y_high[il] = cut_points_y.at(0);
377  y_change = true;
378  }
379  else if (integral_low < fTrimming_totalThreshold &&
380  integral_high > fTrimming_totalThreshold) {
381  y_low[il] = cut_points_y.at(0);
382  y_change = true;
383  }
384  }
385  }
386 
387  MF_LOG_DEBUG("CornerFinderAlg") << "First point in x is " << cut_points_x.at(0);
388 
389  std::sort(cut_points_x.begin(), cut_points_x.end(), compare_to_value(x_bins / 2));
390 
391  MF_LOG_DEBUG("CornerFinderAlg") << "Now the first point in x is " << cut_points_x.at(0);
392 
393  MF_LOG_DEBUG("CornerFinderAlg") << "First point in y is " << cut_points_y.at(0);
394 
395  std::sort(cut_points_y.begin(), cut_points_y.end(), compare_to_value(y_bins / 2));
396 
397  MF_LOG_DEBUG("CornerFinderAlg") << "Now the first point in y is " << cut_points_y.at(0);
398 
399  MF_LOG_DEBUG("CornerFinderAlg")
400  << "\nIntegral on the SW side is "
401  << WireData_histos.at(pid.Plane).Integral(1, cut_points_x.at(0), 1, cut_points_y.at(0))
402  << "\nIntegral on the SE side is "
403  << WireData_histos.at(pid.Plane).Integral(cut_points_x.at(0), x_bins, 1, cut_points_y.at(0))
404  << "\nIntegral on the NW side is "
405  << WireData_histos.at(pid.Plane).Integral(1, cut_points_x.at(0), cut_points_y.at(0), y_bins)
406  << "\nIntegral on the NE side is "
407  << WireData_histos.at(pid.Plane).Integral(
408  cut_points_x.at(0), x_bins, cut_points_y.at(0), y_bins);
409 
410  for (size_t il = 0; il < x_low.size(); il++) {
411 
412  std::stringstream h_name;
413  h_name << "h_" << pid.Cryostat << "_" << pid.TPC << "_" << pid.Plane << "_sub" << il;
414  TH2F h_tmp((h_name.str()).c_str(),
415  "",
416  x_high.at(il) - x_low.at(il) + 1,
417  x_low.at(il),
418  x_high.at(il),
419  y_high.at(il) - y_low.at(il) + 1,
420  y_low.at(il),
421  y_high.at(il));
422 
423  for (int ix = 1; ix <= (x_high.at(il) - x_low.at(il) + 1); ix++) {
424  for (int iy = 1; iy <= (y_high.at(il) - y_low.at(il) + 1); iy++) {
425  h_tmp.SetBinContent(ix,
426  iy,
427  WireData_histos.at(pid.Plane).GetBinContent(x_low.at(il) + (ix - 1),
428  y_low.at(il) + (iy - 1)));
429  }
430  }
431 
432  WireData_trimmed_histos.push_back(
433  std::make_tuple(pid.Plane, h_tmp, x_low.at(il) - 1, y_low.at(il) - 1));
434  }
435 
436  } // end loop over PlaneIDs
437 }
std::vector< std::tuple< int, TH2F, int, int > > WireData_trimmed_histos
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
std::vector< TH1D > WireData_histos_ProjectionY
std::vector< TH2F > WireData_histos
std::vector< TH1D > WireData_histos_ProjectionX
#define MF_LOG_DEBUG(id)
void corner::CornerFinderAlg::get_feature_points ( std::vector< recob::EndPoint2D > &  corner_vector,
geo::Geometry const &  my_geometry 
)

Definition at line 165 of file CornerFinderAlg.cxx.

References attach_feature_points(), geo::GeometryCore::Iterate(), geo::GeometryCore::View(), WireData_histos, and WireData_IDs.

Referenced by trkf::FeatureTracker::produce().

167 {
168  for (auto const& pid : my_geometry.Iterate<geo::PlaneID>()) {
170  WireData_IDs.at(pid.Plane),
171  my_geometry.View(pid),
172  corner_vector);
173  }
174 }
std::vector< std::vector< geo::WireID > > WireData_IDs
void attach_feature_points(TH2F const &h_wire_data, std::vector< geo::WireID > const &wireIDs, geo::View_t view, std::vector< recob::EndPoint2D > &, int startx=0, int starty=0)
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
std::vector< TH2F > WireData_histos
void corner::CornerFinderAlg::get_feature_points_fast ( std::vector< recob::EndPoint2D > &  corner_vector,
geo::Geometry const &  my_geometry 
)

Definition at line 178 of file CornerFinderAlg.cxx.

References attach_feature_points(), create_smaller_histos(), geo::GeometryCore::Iterate(), MF_LOG_DEBUG, WireData_IDs, and WireData_trimmed_histos.

180 {
181  create_smaller_histos(my_geometry);
182 
183  for (auto const& cryostat : my_geometry.Iterate<geo::CryostatGeo>()) {
184  for (unsigned int tpc = 0; tpc < cryostat.NTPC(); ++tpc) {
185  for (size_t histos = 0; histos != WireData_trimmed_histos.size(); histos++) {
186 
187  int plane = std::get<0>(WireData_trimmed_histos.at(histos));
188  int startx = std::get<2>(WireData_trimmed_histos.at(histos));
189  int starty = std::get<3>(WireData_trimmed_histos.at(histos));
190 
191  MF_LOG_DEBUG("CornerFinderAlg") << "Doing histogram " << histos << ", of plane " << plane
192  << " with start points " << startx << " " << starty;
193 
194  attach_feature_points(std::get<1>(WireData_trimmed_histos.at(histos)),
195  WireData_IDs.at(plane),
196  cryostat.TPC(tpc).Plane(plane).View(),
197  corner_vector,
198  startx,
199  starty);
200 
201  MF_LOG_DEBUG("CornerFinderAlg") << "Total feature points now is " << corner_vector.size();
202  }
203  }
204  }
205 }
std::vector< std::vector< geo::WireID > > WireData_IDs
std::vector< std::tuple< int, TH2F, int, int > > WireData_trimmed_histos
void attach_feature_points(TH2F const &h_wire_data, std::vector< geo::WireID > const &wireIDs, geo::View_t view, std::vector< recob::EndPoint2D > &, int startx=0, int starty=0)
void create_smaller_histos(geo::Geometry const &)
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
#define MF_LOG_DEBUG(id)
void corner::CornerFinderAlg::get_feature_points_LineIntegralScore ( std::vector< recob::EndPoint2D > &  corner_vector,
geo::Geometry const &  my_geometry 
)

Definition at line 210 of file CornerFinderAlg.cxx.

References attach_feature_points_LineIntegralScore(), geo::GeometryCore::Iterate(), geo::GeometryCore::View(), WireData_histos, and WireData_IDs.

213 {
214  for (auto const& pid : my_geometry.Iterate<geo::PlaneID>()) {
216  WireData_IDs.at(pid.Plane),
217  my_geometry.View(pid),
218  corner_vector);
219  }
220 }
std::vector< std::vector< geo::WireID > > WireData_IDs
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
std::vector< TH2F > WireData_histos
void attach_feature_points_LineIntegralScore(TH2F const &h_wire_data, std::vector< geo::WireID > const &wireIDs, geo::View_t view, std::vector< recob::EndPoint2D > &)
TH2F const & corner::CornerFinderAlg::GetWireDataHist ( unsigned int  i_plane) const

Definition at line 1183 of file CornerFinderAlg.cxx.

References WireData_histos.

Referenced by trkf::FeatureTracker::GetValidLines().

1184 {
1185  return WireData_histos.at(i_plane);
1186 }
std::vector< TH2F > WireData_histos
void corner::CornerFinderAlg::GrabWires ( std::vector< recob::Wire > const &  wireVec,
geo::Geometry const &  my_geometry 
)

Definition at line 98 of file CornerFinderAlg.cxx.

References geo::GeometryCore::ChannelToWire(), fCalDataModuleLabel, InitializeGeometry(), geo::GeometryCore::Iterate(), geo::GeometryCore::Nplanes(), geo::GeometryCore::Nwires(), geo::PlaneID::Plane, geo::WireID::Wire, WireData_histos, WireData_histos_ProjectionX, WireData_histos_ProjectionY, and WireData_IDs.

Referenced by trkf::FeatureTracker::produce().

100 {
101  InitializeGeometry(my_geometry);
102 
103  const unsigned int nTimeTicks = wireVec.at(0).NSignal();
104 
105  // Initialize the histograms.
106  // All of this should eventually be changed to not need to use histograms...
107  constexpr geo::TPCID tpcid{0, 0};
108  for (auto const& planeid : my_geometry.Iterate<geo::PlaneID>(tpcid)) {
109  auto const i_plane = planeid.Plane;
110 
111  std::stringstream ss_tmp_name, ss_tmp_title;
112  ss_tmp_name << "h_WireData_" << i_plane;
113  ss_tmp_title << fCalDataModuleLabel << " wire data for plane " << i_plane
114  << ";Wire Number;Time Tick";
115 
116  auto const num_wires = my_geometry.Nwires(planeid);
117  if (static_cast<unsigned int>(WireData_histos[i_plane].GetNbinsX()) == num_wires) {
118  WireData_histos[i_plane].Reset();
119  WireData_histos[i_plane].SetName(ss_tmp_name.str().c_str());
120  WireData_histos[i_plane].SetTitle(ss_tmp_title.str().c_str());
121  }
122  else
123  WireData_histos[i_plane] = TH2F(ss_tmp_name.str().c_str(),
124  ss_tmp_title.str().c_str(),
125  num_wires,
126  0,
127  num_wires,
128  nTimeTicks,
129  0,
130  nTimeTicks);
131  }
132 
133  /* Now do the loop over the wires. */
134  for (std::vector<recob::Wire>::const_iterator iwire = wireVec.begin(); iwire < wireVec.end();
135  iwire++) {
136 
137  std::vector<geo::WireID> possible_wireIDs = my_geometry.ChannelToWire(iwire->Channel());
138  geo::WireID this_wireID;
139  try {
140  this_wireID = possible_wireIDs.at(0);
141  }
142  catch (cet::exception& excep) {
143  mf::LogError("CornerFinderAlg") << "Bail out! No Possible Wires!\n";
144  }
145 
146  unsigned int i_plane = this_wireID.Plane;
147  unsigned int i_wire = this_wireID.Wire;
148 
149  WireData_IDs.at(i_plane).at(i_wire) = this_wireID;
150 
151  for (unsigned int i_time = 0; i_time < nTimeTicks; i_time++) {
152  WireData_histos.at(i_plane).SetBinContent(i_wire, i_time, (iwire->Signal()).at(i_time));
153  } //<---End time loop
154 
155  } //<-- End loop over wires
156 
157  for (unsigned int i_plane = 0; i_plane < my_geometry.Nplanes(); i_plane++) {
158  WireData_histos_ProjectionX.at(i_plane) = *(WireData_histos.at(i_plane).ProjectionX());
159  WireData_histos_ProjectionY.at(i_plane) = *(WireData_histos.at(i_plane).ProjectionY());
160  }
161 }
std::vector< std::vector< geo::WireID > > WireData_IDs
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
std::string fCalDataModuleLabel
void InitializeGeometry(geo::Geometry const &)
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:563
intermediate_table::const_iterator const_iterator
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
std::vector< TH1D > WireData_histos_ProjectionY
std::vector< TH2F > WireData_histos
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
std::vector< TH1D > WireData_histos_ProjectionX
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
void corner::CornerFinderAlg::InitializeGeometry ( geo::Geometry const &  my_geometry)
private

Definition at line 71 of file CornerFinderAlg.cxx.

References geo::GeometryCore::Iterate(), geo::GeometryCore::Nplanes(), geo::GeometryCore::Nwires(), WireData_histos, WireData_histos_ProjectionX, WireData_histos_ProjectionY, WireData_IDs, and WireData_trimmed_histos.

Referenced by GrabWires().

72 {
73  // Reset containers
74  WireData_histos.clear();
77  WireData_IDs.clear();
78 
80 
81  // set the sizes of the WireData_histos and WireData_IDs
82  constexpr geo::TPCID tpcid{0, 0};
83  unsigned int nPlanes = my_geometry.Nplanes(tpcid);
84  WireData_histos.resize(nPlanes);
85  WireData_histos_ProjectionX.resize(nPlanes);
86  WireData_histos_ProjectionY.resize(nPlanes);
87 
88  /* For now, we need something to associate each wire in the histogram with a wire_id.
89  This is not a beautiful way of handling this, but for now it should work. */
90  WireData_IDs.resize(nPlanes);
91  for (auto const& planeid : my_geometry.Iterate<geo::PlaneID>(tpcid))
92  WireData_IDs[planeid.Plane].resize(my_geometry.Nwires(planeid));
93 
94  WireData_trimmed_histos.resize(0);
95 }
std::vector< std::vector< geo::WireID > > WireData_IDs
std::vector< std::tuple< int, TH2F, int, int > > WireData_trimmed_histos
The data type to uniquely identify a Plane.
Definition: geo_types.h:463
std::vector< TH1D > WireData_histos_ProjectionY
std::vector< TH2F > WireData_histos
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
std::vector< TH1D > WireData_histos_ProjectionX
float corner::CornerFinderAlg::line_integral ( TH2F const &  hist,
int  x1,
float  y1,
int  x2,
float  y2,
float  threshold 
) const

Definition at line 1084 of file CornerFinderAlg.cxx.

References util::abs(), e, lar_content::minmax(), x1, x2, y1, and y2.

Referenced by calculate_line_integral_score(), and trkf::FeatureTracker::GetValidLines().

1090 {
1091 
1092  int x1 = hist.GetXaxis()->FindBin(begin_x);
1093  int y1 = hist.GetYaxis()->FindBin(begin_y);
1094  int x2 = hist.GetXaxis()->FindBin(end_x);
1095  int y2 = hist.GetYaxis()->FindBin(end_y);
1096 
1097  if (x1 == x2 && abs(y1 - y2) < 1e-5) return 0;
1098 
1099  if (x2 < x1) {
1100  std::swap(x1, x2);
1101  std::swap(y1, y2);
1102  }
1103 
1104  float fraction = 0;
1105  int bin_counter = 0;
1106 
1107  if (x2 != x1) {
1108 
1109  float slope = (y2 - y1) / ((float)(x2 - x1));
1110 
1111  for (int ix = x1; ix <= x2; ix++) {
1112 
1113  int y_min, y_max;
1114 
1115  if (slope >= 0) {
1116  y_min = y1 + slope * (ix - x1);
1117  y_max = y1 + slope * (ix + 1 - x1);
1118  }
1119  else {
1120  y_max = (y1 + 1) + slope * (ix - x1);
1121  y_min = (y1 + 1) + slope * (ix + 1 - x1);
1122  }
1123 
1124  for (int iy = y_min; iy <= y_max; iy++) {
1125  bin_counter++;
1126 
1127  if (hist.GetBinContent(ix, iy) > threshold) fraction += 1.;
1128  }
1129  }
1130  }
1131  else {
1132 
1133  auto const [y_min, y_max] = std::minmax(y1, y2);
1134  for (int iy = y_min; iy <= y_max; iy++) {
1135  bin_counter++;
1136  if (hist.GetBinContent(x1, iy) > threshold) fraction += 1.;
1137  }
1138  }
1139 
1140  return fraction / bin_counter;
1141 }
Float_t y1[n_points_granero]
Definition: compare.C:5
std::pair< float, float > minmax(const float a, const float b)
minmax
Float_t x1[n_points_granero]
Definition: compare.C:5
constexpr auto abs(T v)
Returns the absolute value of the argument.
Float_t y2[n_points_geant4]
Definition: compare.C:26
TH2F * hist
Definition: plot.C:134
void swap(lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &a, lar::deep_const_fwd_iterator_nested< CITER, INNERCONTEXTRACT > &b)
Float_t x2[n_points_geant4]
Definition: compare.C:26
Float_t e
Definition: plot.C:35
std::vector< recob::EndPoint2D > corner::CornerFinderAlg::perform_maximum_suppression ( TH2D const &  h_cornerScore,
std::vector< geo::WireID wireIDs,
geo::View_t  view,
TH2D &  h_maxSuppress,
int  startx = 0,
int  starty = 0 
) const
private

Definition at line 1032 of file CornerFinderAlg.cxx.

References fConversion_bins_per_input_x, fConversion_bins_per_input_y, fMaxSuppress_neighborhood, and fMaxSuppress_threshold.

Referenced by attach_feature_points(), and attach_feature_points_LineIntegralScore().

1039 {
1040  std::vector<recob::EndPoint2D> corner_vector;
1041  const int x_bins = h_cornerScore.GetNbinsX();
1042  const int y_bins = h_cornerScore.GetNbinsY();
1043 
1044  for (int iy = 1; iy <= y_bins; iy++) {
1045  for (int ix = 1; ix <= x_bins; ix++) {
1046 
1047  if (h_cornerScore.GetBinContent(ix, iy) < fMaxSuppress_threshold) continue;
1048 
1049  double temp_max = -1000;
1050  double temp_center_bin = false;
1051 
1052  for (int jx = ix - fMaxSuppress_neighborhood; jx <= ix + fMaxSuppress_neighborhood; jx++) {
1053  for (int jy = iy - fMaxSuppress_neighborhood; jy <= iy + fMaxSuppress_neighborhood; jy++) {
1054 
1055  if (h_cornerScore.GetBinContent(jx, jy) > temp_max) {
1056  temp_max = h_cornerScore.GetBinContent(jx, jy);
1057  if (jx == ix && jy == iy)
1058  temp_center_bin = true;
1059  else {
1060  temp_center_bin = false;
1061  }
1062  }
1063  }
1064  }
1065 
1066  if (temp_center_bin) {
1067 
1068  float time_tick = 0.5 * (float)((2 * (iy + starty)) * fConversion_bins_per_input_y);
1069  int wire_number = ((2 * (ix + startx)) * fConversion_bins_per_input_x) / 2;
1070  double totalQ = 0;
1071  int id = 0;
1073  time_tick, wireIDs[wire_number], h_cornerScore.GetBinContent(ix, iy), id, view, totalQ);
1074  corner_vector.push_back(corner);
1075 
1076  h_maxSuppress.SetBinContent(ix, iy, h_cornerScore.GetBinContent(ix, iy));
1077  }
1078  }
1079  }
1080  return corner_vector;
1081 }

Member Data Documentation

unsigned int corner::CornerFinderAlg::event_number {}
private
std::string corner::CornerFinderAlg::fCalDataModuleLabel
private

Definition at line 62 of file CornerFinderAlg.h.

Referenced by GrabWires().

std::string corner::CornerFinderAlg::fConversion_algorithm
private

Definition at line 63 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_image_histo().

int corner::CornerFinderAlg::fConversion_bins_per_input_x
private
int corner::CornerFinderAlg::fConversion_bins_per_input_y
private
std::string corner::CornerFinderAlg::fConversion_func
private

Definition at line 64 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_image_histo().

int corner::CornerFinderAlg::fConversion_func_neighborhood
private

Definition at line 68 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_image_histo().

float corner::CornerFinderAlg::fConversion_threshold
private

Definition at line 69 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_image_histo().

std::string corner::CornerFinderAlg::fCornerScore_algorithm
private

Definition at line 77 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_cornerScore_histogram().

float corner::CornerFinderAlg::fCornerScore_Harris_kappa
private

Definition at line 79 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_cornerScore_histogram().

int corner::CornerFinderAlg::fCornerScore_neighborhood
private

Definition at line 76 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_cornerScore_histogram().

float corner::CornerFinderAlg::fCornerScore_Noble_epsilon
private

Definition at line 78 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_cornerScore_histogram().

std::string corner::CornerFinderAlg::fDerivative_BlurFunc
private

Definition at line 74 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg().

int corner::CornerFinderAlg::fDerivative_BlurNeighborhood
private

Definition at line 75 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_derivative_histograms().

std::string corner::CornerFinderAlg::fDerivative_method
private

Definition at line 72 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_derivative_histograms().

int corner::CornerFinderAlg::fDerivative_neighborhood
private

Definition at line 73 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_derivative_histograms().

float corner::CornerFinderAlg::fIntegral_bin_threshold
private

Definition at line 82 of file CornerFinderAlg.h.

Referenced by calculate_line_integral_score(), and CornerFinderAlg().

float corner::CornerFinderAlg::fIntegral_fraction_threshold
private

Definition at line 83 of file CornerFinderAlg.h.

Referenced by calculate_line_integral_score(), and CornerFinderAlg().

int corner::CornerFinderAlg::fMaxSuppress_neighborhood
private

Definition at line 80 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and perform_maximum_suppression().

int corner::CornerFinderAlg::fMaxSuppress_threshold
private

Definition at line 81 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and perform_maximum_suppression().

int corner::CornerFinderAlg::fTrimming_buffer
private

Definition at line 66 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_smaller_histos().

float corner::CornerFinderAlg::fTrimming_threshold
private

Definition at line 65 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_smaller_histos().

double corner::CornerFinderAlg::fTrimming_totalThreshold
private

Definition at line 67 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_smaller_histos().

unsigned int corner::CornerFinderAlg::run_number {}
private
std::vector<TH2F> corner::CornerFinderAlg::WireData_histos
private
std::vector<TH1D> corner::CornerFinderAlg::WireData_histos_ProjectionX
private

Definition at line 87 of file CornerFinderAlg.h.

Referenced by create_smaller_histos(), GrabWires(), and InitializeGeometry().

std::vector<TH1D> corner::CornerFinderAlg::WireData_histos_ProjectionY
private

Definition at line 88 of file CornerFinderAlg.h.

Referenced by create_smaller_histos(), GrabWires(), and InitializeGeometry().

std::vector<std::vector<geo::WireID> > corner::CornerFinderAlg::WireData_IDs
private
std::vector<std::tuple<int, TH2F, int, int> > corner::CornerFinderAlg::WireData_trimmed_histos
private

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