LArSoft  v10_04_05
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::WireReadoutGeom const &)
 
void get_feature_points (std::vector< recob::EndPoint2D > &, geo::WireReadoutGeom const &)
 
void get_feature_points_LineIntegralScore (std::vector< recob::EndPoint2D > &, geo::WireReadoutGeom const &)
 
void get_feature_points_fast (std::vector< recob::EndPoint2D > &, geo::GeometryCore const &, geo::WireReadoutGeom 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::WireReadoutGeom 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::WireReadoutGeom 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 28 of file CornerFinderAlg.h.

Constructor & Destructor Documentation

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

Definition at line 42 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().

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

421 {
422  const int x_bins = h_wire_data.GetNbinsX();
423  const float x_min = h_wire_data.GetXaxis()->GetBinLowEdge(1);
424  const float x_max = h_wire_data.GetXaxis()->GetBinUpEdge(x_bins);
425 
426  const int y_bins = h_wire_data.GetNbinsY();
427  const float y_min = h_wire_data.GetYaxis()->GetBinLowEdge(1);
428  const float y_max = h_wire_data.GetYaxis()->GetBinUpEdge(y_bins);
429 
430  const int converted_y_bins = y_bins / fConversion_bins_per_input_y;
431  const int converted_x_bins = x_bins / fConversion_bins_per_input_x;
432 
433  std::stringstream conversion_name;
434  conversion_name << "h_conversion_" << view << "_" << run_number << "_" << event_number;
435  std::stringstream dx_name;
436  dx_name << "h_derivative_x_" << view << "_" << run_number << "_" << event_number;
437  std::stringstream dy_name;
438  dy_name << "h_derivative_y_" << view << "_" << run_number << "_" << event_number;
439  std::stringstream cornerScore_name;
440  cornerScore_name << "h_cornerScore_" << view << "_" << run_number << "_" << event_number;
441  std::stringstream maxSuppress_name;
442  maxSuppress_name << "h_maxSuppress_" << view << "_" << run_number << "_" << event_number;
443 
444  TH2F conversion_histo(conversion_name.str().c_str(),
445  "Image Conversion Histogram",
446  converted_x_bins,
447  x_min,
448  x_max,
449  converted_y_bins,
450  y_min,
451  y_max);
452 
453  TH2F derivativeX_histo(dx_name.str().c_str(),
454  "Partial Derivatives (x)",
455  converted_x_bins,
456  x_min,
457  x_max,
458  converted_y_bins,
459  y_min,
460  y_max);
461 
462  TH2F derivativeY_histo(dy_name.str().c_str(),
463  "Partial Derivatives (y)",
464  converted_x_bins,
465  x_min,
466  x_max,
467  converted_y_bins,
468  y_min,
469  y_max);
470 
471  TH2D cornerScore_histo(cornerScore_name.str().c_str(),
472  "Corner Score",
473  converted_x_bins,
474  x_min,
475  x_max,
476  converted_y_bins,
477  y_min,
478  y_max);
479 
480  TH2D maxSuppress_histo(maxSuppress_name.str().c_str(),
481  "Corner Points (Maximum Suppressed)",
482  converted_x_bins,
483  x_min,
484  x_max,
485  converted_y_bins,
486  y_min,
487  y_max);
488 
489  create_image_histo(h_wire_data, conversion_histo);
490  create_derivative_histograms(conversion_histo, derivativeX_histo, derivativeY_histo);
491  create_cornerScore_histogram(derivativeX_histo, derivativeY_histo, cornerScore_histo);
492  corner_vector = perform_maximum_suppression(
493  cornerScore_histo, wireIDs, view, maxSuppress_histo, startx, starty);
494 }
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 498 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().

503 {
504  const int x_bins = h_wire_data.GetNbinsX();
505  const float x_min = h_wire_data.GetXaxis()->GetBinLowEdge(1);
506  const float x_max = h_wire_data.GetXaxis()->GetBinUpEdge(x_bins);
507 
508  const int y_bins = h_wire_data.GetNbinsY();
509  const float y_min = h_wire_data.GetYaxis()->GetBinLowEdge(1);
510  const float y_max = h_wire_data.GetYaxis()->GetBinUpEdge(y_bins);
511 
512  const int converted_y_bins = y_bins / fConversion_bins_per_input_y;
513  const int converted_x_bins = x_bins / fConversion_bins_per_input_x;
514 
515  std::stringstream conversion_name;
516  conversion_name << "h_conversion_" << view << "_" << run_number << "_" << event_number;
517  std::stringstream dx_name;
518  dx_name << "h_derivative_x_" << view << "_" << run_number << "_" << event_number;
519  std::stringstream dy_name;
520  dy_name << "h_derivative_y_" << view << "_" << run_number << "_" << event_number;
521  std::stringstream cornerScore_name;
522  cornerScore_name << "h_cornerScore_" << view << "_" << run_number << "_" << event_number;
523  std::stringstream maxSuppress_name;
524  maxSuppress_name << "h_maxSuppress_" << view << "_" << run_number << "_" << event_number;
525 
526  TH2F h_conversion(conversion_name.str().c_str(),
527  "Image Conversion Histogram",
528  converted_x_bins,
529  x_min,
530  x_max,
531  converted_y_bins,
532  y_min,
533  y_max);
534  TH2F h_derivative_x(dx_name.str().c_str(),
535  "Partial Derivatives (x)",
536  converted_x_bins,
537  x_min,
538  x_max,
539  converted_y_bins,
540  y_min,
541  y_max);
542  TH2F h_derivative_y(dy_name.str().c_str(),
543  "Partial Derivatives (y)",
544  converted_x_bins,
545  x_min,
546  x_max,
547  converted_y_bins,
548  y_min,
549  y_max);
550  TH2D h_cornerScore(cornerScore_name.str().c_str(),
551  "Feature Point Corner Score",
552  converted_x_bins,
553  x_min,
554  x_max,
555  converted_y_bins,
556  y_min,
557  y_max);
558  TH2D h_maxSuppress(maxSuppress_name.str().c_str(),
559  "Corner Points (Maximum Suppressed)",
560  converted_x_bins,
561  x_min,
562  x_max,
563  converted_y_bins,
564  y_min,
565  y_max);
566 
567  create_image_histo(h_wire_data, h_conversion);
568  create_derivative_histograms(h_conversion, h_derivative_x, h_derivative_y);
569  create_cornerScore_histogram(h_derivative_x, h_derivative_y, h_cornerScore);
570 
571  auto corner_vector_tmp = perform_maximum_suppression(h_cornerScore, wireIDs, view, h_maxSuppress);
572 
573  std::stringstream LI_name;
574  LI_name << "h_lineIntegralScore_" << view << "_" << run_number << "_" << event_number;
575  TH2F h_lineIntegralScore(
576  LI_name.str().c_str(), "Line Integral Score", x_bins, x_min, x_max, y_bins, y_min, y_max);
577  calculate_line_integral_score(h_wire_data, corner_vector_tmp, corner_vector, h_lineIntegralScore);
578 }
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 1114 of file CornerFinderAlg.cxx.

References fIntegral_bin_threshold, fIntegral_fraction_threshold, and line_integral().

Referenced by attach_feature_points_LineIntegralScore().

1119 {
1120  for (auto const i_corner : corner_vector) {
1121 
1122  float score = 0;
1123  for (auto const j_corner : corner_vector) {
1124 
1125  if (line_integral(h_wire_data,
1126  i_corner.WireID().Wire,
1127  i_corner.DriftTime(),
1128  j_corner.WireID().Wire,
1129  j_corner.DriftTime(),
1131  score += 1.;
1132  }
1133  }
1134 
1135  recob::EndPoint2D corner(i_corner.DriftTime(),
1136  i_corner.WireID(),
1137  score,
1138  i_corner.ID(),
1139  i_corner.View(),
1140  i_corner.Charge());
1141 
1142  corner_lineIntegralScore_vector.push_back(corner);
1143 
1144  h_lineIntegralScore.SetBinContent(h_wire_data.GetXaxis()->FindBin(i_corner.WireID().Wire),
1145  h_wire_data.GetYaxis()->FindBin(i_corner.DriftTime()),
1146  score);
1147  }
1148 }
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 930 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().

933 {
934  const int x_bins = h_derivative_x.GetNbinsX();
935  const int y_bins = h_derivative_y.GetNbinsY();
936 
937  //the structure tensor elements
938  double st_xx = 0., st_xy = 0., st_yy = 0.;
939 
940  for (int iy = 1 + fCornerScore_neighborhood; iy <= (y_bins - fCornerScore_neighborhood); iy++) {
941  for (int ix = 1 + fCornerScore_neighborhood; ix <= (x_bins - fCornerScore_neighborhood); ix++) {
942 
943  if (ix == 1 + fCornerScore_neighborhood) {
944  st_xx = 0.;
945  st_xy = 0.;
946  st_yy = 0.;
947 
948  for (int jx = ix - fCornerScore_neighborhood; jx <= ix + fCornerScore_neighborhood; jx++) {
949  for (int jy = iy - fCornerScore_neighborhood; jy <= iy + fCornerScore_neighborhood;
950  jy++) {
951 
952  st_xx += h_derivative_x.GetBinContent(jx, jy) * h_derivative_x.GetBinContent(jx, jy);
953  st_yy += h_derivative_y.GetBinContent(jx, jy) * h_derivative_y.GetBinContent(jx, jy);
954  st_xy += h_derivative_x.GetBinContent(jx, jy) * h_derivative_y.GetBinContent(jx, jy);
955  }
956  }
957  }
958 
959  // we do it this way to reduce computation time
960  else {
961  for (int jy = iy - fCornerScore_neighborhood; jy <= iy + fCornerScore_neighborhood; jy++) {
962 
963  st_xx -= h_derivative_x.GetBinContent(ix - fCornerScore_neighborhood - 1, jy) *
964  h_derivative_x.GetBinContent(ix - fCornerScore_neighborhood - 1, jy);
965  st_xx += h_derivative_x.GetBinContent(ix + fCornerScore_neighborhood, jy) *
966  h_derivative_x.GetBinContent(ix + fCornerScore_neighborhood, jy);
967 
968  st_yy -= h_derivative_y.GetBinContent(ix - fCornerScore_neighborhood - 1, jy) *
969  h_derivative_y.GetBinContent(ix - fCornerScore_neighborhood - 1, jy);
970  st_yy += h_derivative_y.GetBinContent(ix + fCornerScore_neighborhood, jy) *
971  h_derivative_y.GetBinContent(ix + fCornerScore_neighborhood, jy);
972 
973  st_xy -= h_derivative_x.GetBinContent(ix - fCornerScore_neighborhood - 1, jy) *
974  h_derivative_y.GetBinContent(ix - fCornerScore_neighborhood - 1, jy);
975  st_xy += h_derivative_x.GetBinContent(ix + fCornerScore_neighborhood, jy) *
976  h_derivative_y.GetBinContent(ix + fCornerScore_neighborhood, jy);
977  }
978  }
979 
980  if (fCornerScore_algorithm.compare("Noble") == 0) {
981  h_cornerScore.SetBinContent(
982  ix, iy, (st_xx * st_yy - st_xy * st_xy) / (st_xx + st_yy + fCornerScore_Noble_epsilon));
983  }
984  else if (fCornerScore_algorithm.compare("Harris") == 0) {
985  h_cornerScore.SetBinContent(
986  ix,
987  iy,
988  (st_xx * st_yy - st_xy * st_xy) -
989  ((st_xx + st_yy) * (st_xx + st_yy) * fCornerScore_Harris_kappa));
990  }
991  else {
992  mf::LogError("CornerFinderAlg") << "BAD CORNER ALGORITHM: " << fCornerScore_algorithm;
993  return;
994  }
995 
996  } // end for loop over x bins
997  } // end for loop over y bins
998 }
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 649 of file CornerFinderAlg.cxx.

References fDerivative_BlurNeighborhood, fDerivative_method, and fDerivative_neighborhood.

Referenced by attach_feature_points(), and attach_feature_points_LineIntegralScore().

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

583 {
584  double temp_integral = 0;
585 
586  const TF2 fConversion_TF2("fConversion_func", fConversion_func.c_str(), -20, 20, -20, 20);
587 
588  for (int ix = 1; ix <= h_conversion.GetNbinsX(); ix++) {
589  for (int iy = 1; iy <= h_conversion.GetNbinsY(); iy++) {
590 
591  temp_integral = h_wire_data.Integral(ix, ix, iy, iy);
592 
593  if (temp_integral > fConversion_threshold) {
594 
595  if (fConversion_algorithm.compare("binary") == 0)
596  h_conversion.SetBinContent(ix, iy, 10 * fConversion_threshold);
597  else if (fConversion_algorithm.compare("standard") == 0)
598  h_conversion.SetBinContent(ix, iy, temp_integral);
599 
600  else if (fConversion_algorithm.compare("function") == 0) {
601 
602  temp_integral = 0;
603  for (int jx = ix - fConversion_func_neighborhood;
605  jx++) {
606  for (int jy = iy - fConversion_func_neighborhood;
608  jy++) {
609  temp_integral +=
610  h_wire_data.GetBinContent(jx, jy) * fConversion_TF2.Eval(ix - jx, iy - jy);
611  }
612  }
613  h_conversion.SetBinContent(ix, iy, temp_integral);
614  }
615 
616  else if (fConversion_algorithm.compare("skeleton") == 0) {
617 
618  if ((temp_integral > h_wire_data.GetBinContent(ix - 1, iy) &&
619  temp_integral > h_wire_data.GetBinContent(ix + 1, iy)) ||
620  (temp_integral > h_wire_data.GetBinContent(ix, iy - 1) &&
621  temp_integral > h_wire_data.GetBinContent(ix, iy + 1)))
622  h_conversion.SetBinContent(ix, iy, temp_integral);
623  else
624  h_conversion.SetBinContent(ix, iy, fConversion_threshold);
625  }
626  else if (fConversion_algorithm.compare("sk_bin") == 0) {
627 
628  if ((temp_integral > h_wire_data.GetBinContent(ix - 1, iy) &&
629  temp_integral > h_wire_data.GetBinContent(ix + 1, iy)) ||
630  (temp_integral > h_wire_data.GetBinContent(ix, iy - 1) &&
631  temp_integral > h_wire_data.GetBinContent(ix, iy + 1)))
632  h_conversion.SetBinContent(ix, iy, 10 * fConversion_threshold);
633  else
634  h_conversion.SetBinContent(ix, iy, fConversion_threshold);
635  }
636  else
637  h_conversion.SetBinContent(ix, iy, temp_integral);
638  }
639 
640  else
641  h_conversion.SetBinContent(ix, iy, fConversion_threshold);
642  }
643  }
644 }
std::string fConversion_algorithm
void corner::CornerFinderAlg::create_smaller_histos ( geo::WireReadoutGeom const &  wireReadoutGeom)
private

Definition at line 235 of file CornerFinderAlg.cxx.

References fTrimming_buffer, fTrimming_threshold, fTrimming_totalThreshold, geo::Iterable< IterationPolicy, Transform >::Iterate(), MF_LOG_DEBUG, WireData_histos, WireData_histos_ProjectionX, WireData_histos_ProjectionY, and WireData_trimmed_histos.

Referenced by get_feature_points_fast().

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

Definition at line 167 of file CornerFinderAlg.cxx.

References attach_feature_points(), geo::Iterable< IterationPolicy, Transform >::Iterate(), geo::WireReadoutGeom::Plane(), geo::PlaneGeo::View(), WireData_histos, and WireData_IDs.

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

169 {
170  for (auto const& pid : wireReadoutGeom.Iterate<geo::PlaneID>()) {
172  WireData_IDs.at(pid.Plane),
173  wireReadoutGeom.Plane(pid).View(),
174  corner_vector);
175  }
176 }
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:364
std::vector< TH2F > WireData_histos
void corner::CornerFinderAlg::get_feature_points_fast ( std::vector< recob::EndPoint2D > &  corner_vector,
geo::GeometryCore const &  my_geometry,
geo::WireReadoutGeom const &  wireReadoutGeom 
)

Definition at line 180 of file CornerFinderAlg.cxx.

References attach_feature_points(), create_smaller_histos(), geo::Iterable< IterationPolicy, Transform >::Iterate(), MF_LOG_DEBUG, geo::WireReadoutGeom::Plane(), WireData_IDs, and WireData_trimmed_histos.

183 {
184  create_smaller_histos(wireReadoutGeom);
185 
186  for (auto const& cryostat : my_geometry.Iterate<geo::CryostatGeo>()) {
187  for (unsigned int tpc = 0; tpc < cryostat.NTPC(); ++tpc) {
188  for (size_t histos = 0; histos != WireData_trimmed_histos.size(); histos++) {
189 
190  unsigned int plane = std::get<0>(WireData_trimmed_histos.at(histos));
191  int startx = std::get<2>(WireData_trimmed_histos.at(histos));
192  int starty = std::get<3>(WireData_trimmed_histos.at(histos));
193 
194  MF_LOG_DEBUG("CornerFinderAlg") << "Doing histogram " << histos << ", of plane " << plane
195  << " with start points " << startx << " " << starty;
196 
197  attach_feature_points(std::get<1>(WireData_trimmed_histos.at(histos)),
198  WireData_IDs.at(plane),
199  wireReadoutGeom.Plane({cryostat.ID().Cryostat, tpc, plane}).View(),
200  corner_vector,
201  startx,
202  starty);
203 
204  MF_LOG_DEBUG("CornerFinderAlg") << "Total feature points now is " << corner_vector.size();
205  }
206  }
207  }
208 }
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)
Geometry information for a single cryostat.
Definition: CryostatGeo.h:42
void create_smaller_histos(geo::WireReadoutGeom const &)
#define MF_LOG_DEBUG(id)
void corner::CornerFinderAlg::get_feature_points_LineIntegralScore ( std::vector< recob::EndPoint2D > &  corner_vector,
geo::WireReadoutGeom const &  wireReadoutGeom 
)

Definition at line 213 of file CornerFinderAlg.cxx.

References util::abs(), attach_feature_points_LineIntegralScore(), geo::Iterable< IterationPolicy, Transform >::Iterate(), geo::WireReadoutGeom::Plane(), geo::PlaneGeo::View(), WireData_histos, and WireData_IDs.

216 {
217  for (auto const& pid : wireReadoutGeom.Iterate<geo::PlaneID>()) {
219  WireData_IDs.at(pid.Plane),
220  wireReadoutGeom.Plane(pid).View(),
221  corner_vector);
222  }
223 }
std::vector< std::vector< geo::WireID > > WireData_IDs
The data type to uniquely identify a Plane.
Definition: geo_types.h:364
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 1150 of file CornerFinderAlg.cxx.

References WireData_histos.

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

1151 {
1152  return WireData_histos.at(i_plane);
1153 }
std::vector< TH2F > WireData_histos
void corner::CornerFinderAlg::GrabWires ( std::vector< recob::Wire > const &  wireVec,
geo::WireReadoutGeom const &  wireReadoutGeom 
)

Definition at line 100 of file CornerFinderAlg.cxx.

References geo::WireReadoutGeom::ChannelToWire(), fCalDataModuleLabel, InitializeGeometry(), geo::Iterable< IterationPolicy, Transform >::Iterate(), geo::WireReadoutGeom::Nplanes(), geo::WireReadoutGeom::Nwires(), geo::PlaneID::Plane, geo::WireID::Wire, WireData_histos, WireData_histos_ProjectionX, WireData_histos_ProjectionY, and WireData_IDs.

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

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

Definition at line 73 of file CornerFinderAlg.cxx.

References geo::Iterable< IterationPolicy, Transform >::Iterate(), geo::WireReadoutGeom::Nplanes(), geo::WireReadoutGeom::Nwires(), WireData_histos, WireData_histos_ProjectionX, WireData_histos_ProjectionY, WireData_IDs, and WireData_trimmed_histos.

Referenced by GrabWires().

74 {
75  // Reset containers
76  WireData_histos.clear();
79  WireData_IDs.clear();
80 
82 
83  // set the sizes of the WireData_histos and WireData_IDs
84  constexpr geo::TPCID tpcid{0, 0};
85  unsigned int nPlanes = wireReadoutGeom.Nplanes(tpcid);
86  WireData_histos.resize(nPlanes);
87  WireData_histos_ProjectionX.resize(nPlanes);
88  WireData_histos_ProjectionY.resize(nPlanes);
89 
90  /* For now, we need something to associate each wire in the histogram with a wire_id.
91  This is not a beautiful way of handling this, but for now it should work. */
92  WireData_IDs.resize(nPlanes);
93  for (auto const& planeid : wireReadoutGeom.Iterate<geo::PlaneID>(tpcid))
94  WireData_IDs[planeid.Plane].resize(wireReadoutGeom.Nwires(planeid));
95 
96  WireData_trimmed_histos.resize(0);
97 }
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:364
std::vector< TH1D > WireData_histos_ProjectionY
std::vector< TH2F > WireData_histos
The data type to uniquely identify a TPC.
Definition: geo_types.h:306
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 1054 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().

1060 {
1061  int x1 = hist.GetXaxis()->FindBin(begin_x);
1062  int y1 = hist.GetYaxis()->FindBin(begin_y);
1063  int x2 = hist.GetXaxis()->FindBin(end_x);
1064  int y2 = hist.GetYaxis()->FindBin(end_y);
1065 
1066  if (x1 == x2 && abs(y1 - y2) < 1e-5) return 0;
1067 
1068  if (x2 < x1) {
1069  std::swap(x1, x2);
1070  std::swap(y1, y2);
1071  }
1072 
1073  float fraction = 0;
1074  int bin_counter = 0;
1075 
1076  if (x2 != x1) {
1077 
1078  float slope = (y2 - y1) / ((float)(x2 - x1));
1079 
1080  for (int ix = x1; ix <= x2; ix++) {
1081 
1082  int y_min, y_max;
1083 
1084  if (slope >= 0) {
1085  y_min = y1 + slope * (ix - x1);
1086  y_max = y1 + slope * (ix + 1 - x1);
1087  }
1088  else {
1089  y_max = (y1 + 1) + slope * (ix - x1);
1090  y_min = (y1 + 1) + slope * (ix + 1 - x1);
1091  }
1092 
1093  for (int iy = y_min; iy <= y_max; iy++) {
1094  bin_counter++;
1095 
1096  if (hist.GetBinContent(ix, iy) > threshold) fraction += 1.;
1097  }
1098  }
1099  }
1100  else {
1101 
1102  auto const [y_min, y_max] = std::minmax(y1, y2);
1103  for (int iy = y_min; iy <= y_max; iy++) {
1104  bin_counter++;
1105  if (hist.GetBinContent(x1, iy) > threshold) fraction += 1.;
1106  }
1107  }
1108 
1109  return fraction / bin_counter;
1110 }
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 1002 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().

1009 {
1010  std::vector<recob::EndPoint2D> corner_vector;
1011  const int x_bins = h_cornerScore.GetNbinsX();
1012  const int y_bins = h_cornerScore.GetNbinsY();
1013 
1014  for (int iy = 1; iy <= y_bins; iy++) {
1015  for (int ix = 1; ix <= x_bins; ix++) {
1016 
1017  if (h_cornerScore.GetBinContent(ix, iy) < fMaxSuppress_threshold) continue;
1018 
1019  double temp_max = -1000;
1020  double temp_center_bin = false;
1021 
1022  for (int jx = ix - fMaxSuppress_neighborhood; jx <= ix + fMaxSuppress_neighborhood; jx++) {
1023  for (int jy = iy - fMaxSuppress_neighborhood; jy <= iy + fMaxSuppress_neighborhood; jy++) {
1024 
1025  if (h_cornerScore.GetBinContent(jx, jy) > temp_max) {
1026  temp_max = h_cornerScore.GetBinContent(jx, jy);
1027  if (jx == ix && jy == iy)
1028  temp_center_bin = true;
1029  else {
1030  temp_center_bin = false;
1031  }
1032  }
1033  }
1034  }
1035 
1036  if (temp_center_bin) {
1037 
1038  float time_tick = 0.5 * (float)((2 * (iy + starty)) * fConversion_bins_per_input_y);
1039  int wire_number = ((2 * (ix + startx)) * fConversion_bins_per_input_x) / 2;
1040  double totalQ = 0;
1041  int id = 0;
1043  time_tick, wireIDs[wire_number], h_cornerScore.GetBinContent(ix, iy), id, view, totalQ);
1044  corner_vector.push_back(corner);
1045 
1046  h_maxSuppress.SetBinContent(ix, iy, h_cornerScore.GetBinContent(ix, iy));
1047  }
1048  }
1049  }
1050  return corner_vector;
1051 }

Member Data Documentation

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

Definition at line 58 of file CornerFinderAlg.h.

Referenced by GrabWires().

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

Definition at line 59 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 60 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_image_histo().

int corner::CornerFinderAlg::fConversion_func_neighborhood
private

Definition at line 64 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_image_histo().

float corner::CornerFinderAlg::fConversion_threshold
private

Definition at line 65 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_image_histo().

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

Definition at line 73 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_cornerScore_histogram().

float corner::CornerFinderAlg::fCornerScore_Harris_kappa
private

Definition at line 75 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_cornerScore_histogram().

int corner::CornerFinderAlg::fCornerScore_neighborhood
private

Definition at line 72 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_cornerScore_histogram().

float corner::CornerFinderAlg::fCornerScore_Noble_epsilon
private

Definition at line 74 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_cornerScore_histogram().

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

Definition at line 70 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg().

int corner::CornerFinderAlg::fDerivative_BlurNeighborhood
private

Definition at line 71 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_derivative_histograms().

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

Definition at line 68 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_derivative_histograms().

int corner::CornerFinderAlg::fDerivative_neighborhood
private

Definition at line 69 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_derivative_histograms().

float corner::CornerFinderAlg::fIntegral_bin_threshold
private

Definition at line 78 of file CornerFinderAlg.h.

Referenced by calculate_line_integral_score(), and CornerFinderAlg().

float corner::CornerFinderAlg::fIntegral_fraction_threshold
private

Definition at line 79 of file CornerFinderAlg.h.

Referenced by calculate_line_integral_score(), and CornerFinderAlg().

int corner::CornerFinderAlg::fMaxSuppress_neighborhood
private

Definition at line 76 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and perform_maximum_suppression().

int corner::CornerFinderAlg::fMaxSuppress_threshold
private

Definition at line 77 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and perform_maximum_suppression().

int corner::CornerFinderAlg::fTrimming_buffer
private

Definition at line 62 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_smaller_histos().

float corner::CornerFinderAlg::fTrimming_threshold
private

Definition at line 61 of file CornerFinderAlg.h.

Referenced by CornerFinderAlg(), and create_smaller_histos().

double corner::CornerFinderAlg::fTrimming_totalThreshold
private

Definition at line 63 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 83 of file CornerFinderAlg.h.

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

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

Definition at line 84 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: