All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CornerFinderAlg.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // CornerFinderAlg class
4 //
5 // wketchum@fnal.gov
6 //
7 // CornerFinder is meant to use image-processing techniques (mainly Harris-Stephens
8 // corner-finding) to find "corners" using the information from calibrated wires.
9 //
10 // Conversion_algorithm options:
11 // standard --- basically a copy of the calibrated wires
12 // skeleton --- a thinned copy of the calibrated wires
13 // binary --- ticks above threshold get assigned a value 10*threshold, everything else = threshold
14 // function --- apply a function (like a double-Gaussian) to a neighborhood around each tick
15 //
16 // Derivative options:
17 // Sobel --- apply a Sobel mask (neighborhood of 1 or 2 supported)
18 // local --- take slope from immediately neighboring bins (neighborhood of 1 supported)
19 //
20 // Derivative_BlurFunc options: none. You're stuck with a double gaussian.
21 //
22 // CornerScore_algorithm options:
23 // Noble --- determinant / (trace + Noble_epsilon)
24 // Harris --- determinant - (trace)^2 * Harris_kappa
25 ////////////////////////////////////////////////////////////////////////
26 
27 
28 #include "larreco/RecoAlg/CornerFinderAlg.h"
29 
30 #include "messagefacility/MessageLogger/MessageLogger.h"
31 
32 #include "larcore/Geometry/Geometry.h"
33 #include "larcorealg/Geometry/CryostatGeo.h"
34 #include "larcorealg/Geometry/TPCGeo.h"
35 #include "larcorealg/Geometry/PlaneGeo.h"
36 
37 // NOTE: In the .h file I assumed this would belong in the cluster class....if
38 // we decide otherwise we will need to search and replace for this
39 
40 
41 //-----------------------------------------------------------------------------
42 corner::CornerFinderAlg::CornerFinderAlg(fhicl::ParameterSet const& pset)
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 {
65  fTrimming_buffer = std::max({ fConversion_func_neighborhood,
66  fDerivative_neighborhood,
67  fDerivative_BlurNeighborhood,
68  fCornerScore_neighborhood,
69  fMaxSuppress_neighborhood });
70 }
71 
72 //-----------------------------------------------------------------------------
73 void corner::CornerFinderAlg::InitializeGeometry(geo::Geometry const& my_geometry){
74  // Reset containers
75  WireData_histos.clear();
76  WireData_histos_ProjectionX.clear();
77  WireData_histos_ProjectionY.clear();
78  WireData_IDs.clear();
79 
80  WireData_trimmed_histos.clear();
81 
82  // set the sizes of the WireData_histos and WireData_IDs
83  unsigned int nPlanes = my_geometry.Nplanes();
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(unsigned int i_plane=0; i_plane < nPlanes; ++i_plane)
92  WireData_IDs.at(i_plane).resize(my_geometry.Nwires(i_plane));
93 
94  WireData_trimmed_histos.resize(0);
95 
96 }
97 
98 //-----------------------------------------------------------------------------
99 void corner::CornerFinderAlg::GrabWires( std::vector<recob::Wire> const& wireVec, geo::Geometry const& my_geometry){
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  for (unsigned int i_plane=0; i_plane < my_geometry.Nplanes(); i_plane++){
108 
109  std::stringstream ss_tmp_name,ss_tmp_title;
110  ss_tmp_name << "h_WireData_" << i_plane;
111  ss_tmp_title << fCalDataModuleLabel << " wire data for plane " << i_plane << ";Wire Number;Time Tick";
112 
113  if( (unsigned int)(WireData_histos.at(i_plane).GetNbinsX()) == (my_geometry.Nwires(i_plane)) ) {
114  WireData_histos.at(i_plane).Reset();
115  WireData_histos.at(i_plane).SetName(ss_tmp_name.str().c_str());
116  WireData_histos.at(i_plane).SetTitle(ss_tmp_title.str().c_str());
117  }
118  else
119  WireData_histos.at(i_plane) = TH2F(ss_tmp_name.str().c_str(),
120  ss_tmp_title.str().c_str(),
121  my_geometry.Nwires(i_plane),
122  0,
123  my_geometry.Nwires(i_plane),
124  nTimeTicks,
125  0,
126  nTimeTicks);
127  }
128 
129 
130  /* Now do the loop over the wires. */
131  for (std::vector<recob::Wire>::const_iterator iwire = wireVec.begin(); iwire < wireVec.end(); iwire++) {
132 
133  std::vector<geo::WireID> possible_wireIDs = my_geometry.ChannelToWire(iwire->Channel());
134  geo::WireID this_wireID;
135  try { this_wireID = possible_wireIDs.at(0);}
136  catch(cet::exception& excep) { mf::LogError("CornerFinderAlg") << "Bail out! No Possible Wires!\n"; }
137 
138  unsigned int i_plane = this_wireID.Plane;
139  unsigned int i_wire = this_wireID.Wire;
140 
141  WireData_IDs.at(i_plane).at(i_wire) = this_wireID;
142 
143  for(unsigned int i_time = 0; i_time < nTimeTicks; i_time++){
144  WireData_histos.at(i_plane).SetBinContent(i_wire,i_time,(iwire->Signal()).at(i_time));
145  }//<---End time loop
146 
147  }//<-- End loop over wires
148 
149 
150  for (unsigned int i_plane=0; i_plane < my_geometry.Nplanes(); i_plane++){
151  WireData_histos_ProjectionX.at(i_plane) = *(WireData_histos.at(i_plane).ProjectionX());
152  WireData_histos_ProjectionY.at(i_plane) = *(WireData_histos.at(i_plane).ProjectionY());
153  }
154 
155 
156 }
157 
158 
159 //-----------------------------------------------------------------------------------
160 // This gives us a vecotr of EndPoint2D objects that correspond to possible corners
161 void corner::CornerFinderAlg::get_feature_points(std::vector<recob::EndPoint2D> & corner_vector,
162  geo::Geometry const& my_geometry){
163  for(auto const& pid : my_geometry.IteratePlaneIDs()){
164  attach_feature_points(WireData_histos.at(pid.Plane),
165  WireData_IDs.at(pid.Plane),
166  my_geometry.View(pid),
167  corner_vector);
168  }
169 }
170 
171 //-----------------------------------------------------------------------------------
172 // This gives us a vector of EndPoint2D objects that correspond to possible corners, but quickly!
173 void corner::CornerFinderAlg::get_feature_points_fast(std::vector<recob::EndPoint2D> & corner_vector,
174  geo::Geometry const& my_geometry){
175  create_smaller_histos(my_geometry);
176 
177  for(unsigned int cstat = 0; cstat < my_geometry.Ncryostats(); ++cstat){
178  for(unsigned int tpc = 0; tpc < my_geometry.Cryostat(cstat).NTPC(); ++tpc){
179  for(size_t histos=0; histos!= WireData_trimmed_histos.size(); histos++){
180 
181  int plane = std::get<0>(WireData_trimmed_histos.at(histos));
182  int startx = std::get<2>(WireData_trimmed_histos.at(histos));
183  int starty = std::get<3>(WireData_trimmed_histos.at(histos));
184 
185  MF_LOG_DEBUG("CornerFinderAlg")
186  << "Doing histogram " << histos
187  << ", of plane " << plane
188  << " with start points " << startx << " " << starty;
189 
190  attach_feature_points(std::get<1>(WireData_trimmed_histos.at(histos)),
191  WireData_IDs.at(plane),
192  my_geometry.Cryostat(cstat).TPC(tpc).Plane(plane).View(),
193  corner_vector,
194  startx,
195  starty);
196 
197  MF_LOG_DEBUG("CornerFinderAlg") << "Total feature points now is " << corner_vector.size();
198  }
199  }
200  }
201 
202 }
203 
204 //-----------------------------------------------------------------------------------
205 // This gives us a vecotr of EndPoint2D objects that correspond to possible corners
206 // Uses line integral score as corner strength
207 void corner::CornerFinderAlg::get_feature_points_LineIntegralScore(std::vector<recob::EndPoint2D> & corner_vector,
208  geo::Geometry const& my_geometry){
209  for(auto const& pid : my_geometry.IteratePlaneIDs()){
210  attach_feature_points_LineIntegralScore(WireData_histos.at(pid.Plane),
211  WireData_IDs.at(pid.Plane),
212  my_geometry.View(pid),
213  corner_vector);
214  }
215 }
216 
217 struct compare_to_value{
218 
219  compare_to_value(int b) {this->b = b;}
220  bool operator() (int i, int j) {
221  return std::abs(b-i)<std::abs(b-j);
222  }
223 
224  int b;
225 
226 };
227 
228 struct compare_to_range{
229 
230  compare_to_range(int a, int b) {this->a = a; this->b = b;}
231  bool operator() (int i, int j) {
232 
233  int mid = (b-a)/2 + a;
234  if(i>=a && i<=b && j>=a && j<=b)
235  return std::abs(mid-i)<std::abs(mid-j);
236 
237  else if(j>=a && j<=b && (i<a || i>b) )
238  return false;
239 
240  else if(i>=a && i<=b && (j<a || j>b) )
241  return true;
242 
243  else
244  return true;
245  }
246 
247  int a;
248  int b;
249 
250 };
251 
252 //-----------------------------------------------------------------------------
253 // This looks for areas of the wires that are non-noise, to speed up evaluation
254 void corner::CornerFinderAlg::create_smaller_histos(geo::Geometry const& my_geometry){
255 
256  for(auto const& pid : my_geometry.IteratePlaneIDs() ){
257 
258  MF_LOG_DEBUG("CornerFinderAlg")
259  << "Working plane " << pid.Plane << ".";
260 
261  int x_bins = WireData_histos_ProjectionX.at(pid.Plane).GetNbinsX();
262  int y_bins = WireData_histos_ProjectionY.at(pid.Plane).GetNbinsX();
263 
264  std::vector<int> cut_points_x {0};
265  std::vector<int> cut_points_y {0};
266 
267  for (int ix=1; ix<=x_bins; ix++){
268 
269  float this_value = WireData_histos_ProjectionX.at(pid.Plane).GetBinContent(ix);
270 
271  if(ix<fTrimming_buffer || ix>(x_bins-fTrimming_buffer)) continue;
272 
273  int jx=ix-fTrimming_buffer;
274  while(this_value<fTrimming_threshold){
275  if(jx==ix+fTrimming_buffer) break;
276  this_value = WireData_histos_ProjectionX.at(pid.Plane).GetBinContent(jx);
277  jx++;
278  }
279  if(this_value<fTrimming_threshold){
280  cut_points_x.push_back(ix);
281  }
282 
283  }
284 
285  for (int iy=1; iy<=y_bins; iy++){
286 
287  float this_value = WireData_histos_ProjectionY.at(pid.Plane).GetBinContent(iy);
288 
289  if(iy<fTrimming_buffer || iy>(y_bins-fTrimming_buffer)) continue;
290 
291  int jy=iy-fTrimming_buffer;
292  while(this_value<fTrimming_threshold){
293  if(jy==iy+fTrimming_buffer) break;
294  this_value = WireData_histos_ProjectionY.at(pid.Plane).GetBinContent(jy);
295  jy++;
296  }
297  if(this_value<fTrimming_threshold){
298  cut_points_y.push_back(iy);
299  }
300 
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))
326  continue;
327 
328  double integral_low = WireData_histos.at(pid.Plane).Integral(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(cut_points_x.at(0),x_high.at(il),y_low.at(il),y_high.at(il));
330  if(integral_low > fTrimming_totalThreshold && integral_high > fTrimming_totalThreshold){
331  x_low.push_back(cut_points_x.at(0));
332  x_high.push_back(x_high.at(il));
333  y_low.push_back(y_low.at(il));
334  y_high.push_back(y_high.at(il));
335 
336  x_high[il] = cut_points_x.at(0);
337  x_change = true;
338  }
339  else if(integral_low > fTrimming_totalThreshold && integral_high < fTrimming_totalThreshold){
340  x_high[il] = cut_points_x.at(0);
341  x_change = true;
342  }
343  else if(integral_low < fTrimming_totalThreshold && integral_high > fTrimming_totalThreshold){
344  x_low[il] = cut_points_x.at(0);
345  x_change = true;
346  }
347  }
348 
349  current_size = x_low.size();
350 
351  for(size_t il=0; il<current_size; il++){
352 
353  int comp_value = (y_high.at(il) - y_low.at(il)) / 2;
354  std::sort(cut_points_y.begin(),cut_points_y.end(),compare_to_value(comp_value));
355 
356  if(cut_points_y.at(0) <= y_low.at(il) || cut_points_y.at(0) >= y_high.at(il))
357  continue;
358 
359  double integral_low = WireData_histos.at(pid.Plane).Integral(x_low.at(il),x_high.at(il),y_low.at(il),cut_points_y.at(0));
360  double integral_high = WireData_histos.at(pid.Plane).Integral(x_low.at(il),x_high.at(il),cut_points_y.at(0),y_high.at(il));
361  if(integral_low > fTrimming_totalThreshold && integral_high > fTrimming_totalThreshold){
362  y_low.push_back(cut_points_y.at(0));
363  y_high.push_back(y_high.at(il));
364  x_low.push_back(x_low.at(il));
365  x_high.push_back(x_high.at(il));
366 
367  y_high[il] = cut_points_y.at(0);
368  y_change = true;
369  }
370  else if(integral_low > fTrimming_totalThreshold && integral_high < fTrimming_totalThreshold){
371  y_high[il] = cut_points_y.at(0);
372  y_change = true;
373  }
374  else if(integral_low < fTrimming_totalThreshold && integral_high > fTrimming_totalThreshold){
375  y_low[il] = cut_points_y.at(0);
376  y_change = true;
377  }
378  }
379 
380  }
381 
382  MF_LOG_DEBUG("CornerFinderAlg")
383  << "First point in x is " << cut_points_x.at(0);
384 
385  std::sort(cut_points_x.begin(),cut_points_x.end(),compare_to_value(x_bins/2));
386 
387  MF_LOG_DEBUG("CornerFinderAlg")
388  << "Now the first point in x is " << cut_points_x.at(0);
389 
390  MF_LOG_DEBUG("CornerFinderAlg")
391  << "First point in y is " << cut_points_y.at(0);
392 
393  std::sort(cut_points_y.begin(),cut_points_y.end(),compare_to_value(y_bins/2));
394 
395  MF_LOG_DEBUG("CornerFinderAlg")
396  << "Now the first point in y is " << cut_points_y.at(0);
397 
398  MF_LOG_DEBUG("CornerFinderAlg")
399  << "\nIntegral on the SW side is "
400  << WireData_histos.at(pid.Plane).Integral(1,cut_points_x.at(0),1,cut_points_y.at(0))
401  << "\nIntegral on the SE side is "
402  << WireData_histos.at(pid.Plane).Integral(cut_points_x.at(0),x_bins,1,cut_points_y.at(0))
403  << "\nIntegral on the NW side is "
404  << WireData_histos.at(pid.Plane).Integral(1,cut_points_x.at(0),cut_points_y.at(0),y_bins)
405  << "\nIntegral on the NE side is "
406  << WireData_histos.at(pid.Plane).Integral(cut_points_x.at(0),x_bins,cut_points_y.at(0),y_bins);
407 
408 
409  for(size_t il=0; il<x_low.size(); il++){
410 
411  std::stringstream h_name;
412  h_name << "h_" << pid.Cryostat << "_" << pid.TPC << "_" << pid.Plane << "_sub" << il;
413  TH2F h_tmp( (h_name.str()).c_str(),"",
414  x_high.at(il)-x_low.at(il)+1,x_low.at(il),x_high.at(il),
415  y_high.at(il)-y_low.at(il)+1,y_low.at(il),y_high.at(il));
416 
417  for(int ix=1; ix<=(x_high.at(il)-x_low.at(il)+1); ix++){
418  for(int iy=1; iy<=(y_high.at(il)-y_low.at(il)+1); iy++){
419  h_tmp.SetBinContent(ix,iy,WireData_histos.at(pid.Plane).GetBinContent(x_low.at(il)+(ix-1),y_low.at(il)+(iy-1)));
420  }
421  }
422 
423  WireData_trimmed_histos.push_back(std::make_tuple(pid.Plane,h_tmp,x_low.at(il)-1,y_low.at(il)-1));
424  }
425 
426  }// end loop over PlaneIDs
427 
428 
429 }
430 
431 //-----------------------------------------------------------------------------
432 // This puts on all the feature points in a given view, using a given data histogram
433 void corner::CornerFinderAlg::attach_feature_points( TH2F const& h_wire_data,
434  std::vector<geo::WireID> const& wireIDs,
435  geo::View_t view,
436  std::vector<recob::EndPoint2D> & corner_vector,
437  int startx,
438  int starty){
439 
440 
441  const int x_bins = h_wire_data.GetNbinsX();
442  const float x_min = h_wire_data.GetXaxis()->GetBinLowEdge(1);
443  const float x_max = h_wire_data.GetXaxis()->GetBinUpEdge(x_bins);
444 
445  const int y_bins = h_wire_data.GetNbinsY();
446  const float y_min = h_wire_data.GetYaxis()->GetBinLowEdge(1);
447  const float y_max = h_wire_data.GetYaxis()->GetBinUpEdge(y_bins);
448 
449  const int converted_y_bins = y_bins/fConversion_bins_per_input_y;
450  const int converted_x_bins = x_bins/fConversion_bins_per_input_x;
451 
452  std::stringstream conversion_name; conversion_name << "h_conversion_" << view << "_" << run_number << "_" << event_number;
453  std::stringstream dx_name; dx_name << "h_derivative_x_" << view << "_" << run_number << "_" << event_number;
454  std::stringstream dy_name; dy_name << "h_derivative_y_" << view << "_" << run_number << "_" << event_number;
455  std::stringstream cornerScore_name; cornerScore_name << "h_cornerScore_" << view << "_" << run_number << "_" << event_number;
456  std::stringstream maxSuppress_name; maxSuppress_name << "h_maxSuppress_" << view << "_" << run_number << "_" << event_number;
457 
458  TH2F conversion_histo(conversion_name.str().c_str(),"Image Conversion Histogram",
459  converted_x_bins,x_min,x_max,
460  converted_y_bins,y_min,y_max);
461 
462  TH2F derivativeX_histo(dx_name.str().c_str(),"Partial Derivatives (x)",
463  converted_x_bins,x_min,x_max,
464  converted_y_bins,y_min,y_max);
465 
466  TH2F derivativeY_histo(dy_name.str().c_str(),"Partial Derivatives (y)",
467  converted_x_bins,x_min,x_max,
468  converted_y_bins,y_min,y_max);
469 
470  TH2D cornerScore_histo(cornerScore_name.str().c_str(),"Corner Score",
471  converted_x_bins,x_min,x_max,
472  converted_y_bins,y_min,y_max);
473 
474  TH2D maxSuppress_histo(maxSuppress_name.str().c_str(),"Corner Points (Maximum Suppressed)",
475  converted_x_bins,x_min,x_max,
476  converted_y_bins,y_min,y_max);
477 
478  create_image_histo(h_wire_data,conversion_histo);
479  create_derivative_histograms(conversion_histo,derivativeX_histo,derivativeY_histo);
480  create_cornerScore_histogram(derivativeX_histo,derivativeY_histo,cornerScore_histo);
481  corner_vector = perform_maximum_suppression(cornerScore_histo,wireIDs,view,maxSuppress_histo,startx,starty);
482 }
483 
484 
485 //-----------------------------------------------------------------------------
486 // This puts on all the feature points in a given view, using a given data histogram
487 void corner::CornerFinderAlg::attach_feature_points_LineIntegralScore(TH2F const& h_wire_data,
488  std::vector<geo::WireID> const& wireIDs,
489  geo::View_t view,
490  std::vector<recob::EndPoint2D> & corner_vector)
491 {
492  const int x_bins = h_wire_data.GetNbinsX();
493  const float x_min = h_wire_data.GetXaxis()->GetBinLowEdge(1);
494  const float x_max = h_wire_data.GetXaxis()->GetBinUpEdge(x_bins);
495 
496  const int y_bins = h_wire_data.GetNbinsY();
497  const float y_min = h_wire_data.GetYaxis()->GetBinLowEdge(1);
498  const float y_max = h_wire_data.GetYaxis()->GetBinUpEdge(y_bins);
499 
500  const int converted_y_bins = y_bins/fConversion_bins_per_input_y;
501  const int converted_x_bins = x_bins/fConversion_bins_per_input_x;
502 
503  std::stringstream conversion_name; conversion_name << "h_conversion_" << view << "_" << run_number << "_" << event_number;
504  std::stringstream dx_name; dx_name << "h_derivative_x_" << view << "_" << run_number << "_" << event_number;
505  std::stringstream dy_name; dy_name << "h_derivative_y_" << view << "_" << run_number << "_" << event_number;
506  std::stringstream cornerScore_name; cornerScore_name << "h_cornerScore_" << view << "_" << run_number << "_" << event_number;
507  std::stringstream maxSuppress_name; maxSuppress_name << "h_maxSuppress_" << view << "_" << run_number << "_" << event_number;
508 
509  TH2F h_conversion (conversion_name.str().c_str(),
510  "Image Conversion Histogram",
511  converted_x_bins,x_min,x_max,
512  converted_y_bins,y_min,y_max);
513  TH2F h_derivative_x(dx_name.str().c_str(),
514  "Partial Derivatives (x)",
515  converted_x_bins,x_min,x_max,
516  converted_y_bins,y_min,y_max);
517  TH2F h_derivative_y(dy_name.str().c_str(),
518  "Partial Derivatives (y)",
519  converted_x_bins,x_min,x_max,
520  converted_y_bins,y_min,y_max);
521  TH2D h_cornerScore (cornerScore_name.str().c_str(),
522  "Feature Point Corner Score",
523  converted_x_bins,x_min,x_max,
524  converted_y_bins,y_min,y_max);
525  TH2D h_maxSuppress (maxSuppress_name.str().c_str(),
526  "Corner Points (Maximum Suppressed)",
527  converted_x_bins,x_min,x_max,
528  converted_y_bins,y_min,y_max);
529 
530  create_image_histo(h_wire_data,h_conversion);
531  create_derivative_histograms(h_conversion,h_derivative_x,h_derivative_y);
532  create_cornerScore_histogram(h_derivative_x,h_derivative_y,h_cornerScore);
533 
534  auto corner_vector_tmp = perform_maximum_suppression(h_cornerScore,wireIDs,view,h_maxSuppress);
535 
536  std::stringstream LI_name; LI_name << "h_lineIntegralScore_" << view << "_" << run_number << "_" << event_number;
537  TH2F h_lineIntegralScore(LI_name.str().c_str(),
538  "Line Integral Score",
539  x_bins,x_min,x_max,
540  y_bins,y_min,y_max);
541  calculate_line_integral_score(h_wire_data,corner_vector_tmp,corner_vector,h_lineIntegralScore);
542 
543 }
544 
545 
546 //-----------------------------------------------------------------------------
547 // Convert to pixel
548 void corner::CornerFinderAlg::create_image_histo(TH2F const& h_wire_data, TH2F & h_conversion) const {
549 
550  double temp_integral=0;
551 
552  const TF2 fConversion_TF2("fConversion_func",fConversion_func.c_str(),-20,20,-20,20);
553 
554  for(int ix=1; ix<=h_conversion.GetNbinsX(); ix++){
555  for(int iy=1; iy<=h_conversion.GetNbinsY(); iy++){
556 
557  temp_integral = h_wire_data.Integral(ix,ix,iy,iy);
558 
559  if( temp_integral > fConversion_threshold){
560 
561  if(fConversion_algorithm.compare("binary")==0)
562  h_conversion.SetBinContent(ix,iy,10*fConversion_threshold);
563  else if(fConversion_algorithm.compare("standard")==0)
564  h_conversion.SetBinContent(ix,iy,temp_integral);
565 
566  else if(fConversion_algorithm.compare("function")==0){
567 
568  temp_integral = 0;
569  for(int jx=ix-fConversion_func_neighborhood; jx<=ix+fConversion_func_neighborhood; jx++){
570  for(int jy=iy-fConversion_func_neighborhood; jy<=iy+fConversion_func_neighborhood; jy++){
571  temp_integral += h_wire_data.GetBinContent(jx,jy)*fConversion_TF2.Eval(ix-jx,iy-jy);
572  }
573  }
574  h_conversion.SetBinContent(ix,iy,temp_integral);
575  }
576 
577  else if(fConversion_algorithm.compare("skeleton")==0){
578 
579  if( (temp_integral > h_wire_data.GetBinContent(ix-1,iy) && temp_integral > h_wire_data.GetBinContent(ix+1,iy))
580  || (temp_integral > h_wire_data.GetBinContent(ix,iy-1) && temp_integral > h_wire_data.GetBinContent(ix,iy+1)))
581  h_conversion.SetBinContent(ix,iy,temp_integral);
582  else
583  h_conversion.SetBinContent(ix,iy,fConversion_threshold);
584  }
585  else if(fConversion_algorithm.compare("sk_bin")==0){
586 
587  if( (temp_integral > h_wire_data.GetBinContent(ix-1,iy) && temp_integral > h_wire_data.GetBinContent(ix+1,iy))
588  || (temp_integral > h_wire_data.GetBinContent(ix,iy-1) && temp_integral > h_wire_data.GetBinContent(ix,iy+1)))
589  h_conversion.SetBinContent(ix,iy,10*fConversion_threshold);
590  else
591  h_conversion.SetBinContent(ix,iy,fConversion_threshold);
592  }
593  else
594  h_conversion.SetBinContent(ix,iy,temp_integral);
595  }
596 
597  else
598  h_conversion.SetBinContent(ix,iy,fConversion_threshold);
599 
600  }
601  }
602 
603 }
604 
605 //-----------------------------------------------------------------------------
606 // Derivative
607 
608 void corner::CornerFinderAlg::create_derivative_histograms(TH2F const& h_conversion, TH2F & h_derivative_x, TH2F & h_derivative_y){
609 
610  const int x_bins = h_conversion.GetNbinsX();
611  const int y_bins = h_conversion.GetNbinsY();
612 
613  for(int iy=1+fDerivative_neighborhood; iy<=(y_bins-fDerivative_neighborhood); iy++){
614  for(int ix=1+fDerivative_neighborhood; ix<=(x_bins-fDerivative_neighborhood); ix++){
615 
616  if(fDerivative_method.compare("Sobel")==0){
617 
618  if(fDerivative_neighborhood==1){
619  h_derivative_x.SetBinContent(ix,iy,
620  0.5*(h_conversion.GetBinContent(ix+1,iy)-h_conversion.GetBinContent(ix-1,iy))
621  + 0.25*(h_conversion.GetBinContent(ix+1,iy+1)-h_conversion.GetBinContent(ix-1,iy+1))
622  + 0.25*(h_conversion.GetBinContent(ix+1,iy-1)-h_conversion.GetBinContent(ix-1,iy-1)));
623  h_derivative_y.SetBinContent(ix,iy,
624  0.5*(h_conversion.GetBinContent(ix,iy+1)-h_conversion.GetBinContent(ix,iy-1))
625  + 0.25*(h_conversion.GetBinContent(ix-1,iy+1)-h_conversion.GetBinContent(ix-1,iy-1))
626  + 0.25*(h_conversion.GetBinContent(ix+1,iy+1)-h_conversion.GetBinContent(ix+1,iy-1)));
627  }
628  else if(fDerivative_neighborhood==2){
629  h_derivative_x.SetBinContent(ix,iy,
630  12*(h_conversion.GetBinContent(ix+1,iy)-h_conversion.GetBinContent(ix-1,iy))
631  + 8*(h_conversion.GetBinContent(ix+1,iy+1)-h_conversion.GetBinContent(ix-1,iy+1))
632  + 8*(h_conversion.GetBinContent(ix+1,iy-1)-h_conversion.GetBinContent(ix-1,iy-1))
633  + 2*(h_conversion.GetBinContent(ix+1,iy+2)-h_conversion.GetBinContent(ix-1,iy+2))
634  + 2*(h_conversion.GetBinContent(ix+1,iy-2)-h_conversion.GetBinContent(ix-1,iy-2))
635  + 6*(h_conversion.GetBinContent(ix+2,iy)-h_conversion.GetBinContent(ix-2,iy))
636  + 4*(h_conversion.GetBinContent(ix+2,iy+1)-h_conversion.GetBinContent(ix-2,iy+1))
637  + 4*(h_conversion.GetBinContent(ix+2,iy-1)-h_conversion.GetBinContent(ix-2,iy-1))
638  + 1*(h_conversion.GetBinContent(ix+2,iy+2)-h_conversion.GetBinContent(ix-2,iy+2))
639  + 1*(h_conversion.GetBinContent(ix+2,iy-2)-h_conversion.GetBinContent(ix-2,iy-2)));
640  h_derivative_y.SetBinContent(ix,iy,
641  12*(h_conversion.GetBinContent(ix,iy+1)-h_conversion.GetBinContent(ix,iy-1))
642  + 8*(h_conversion.GetBinContent(ix-1,iy+1)-h_conversion.GetBinContent(ix-1,iy-1))
643  + 8*(h_conversion.GetBinContent(ix+1,iy+1)-h_conversion.GetBinContent(ix+1,iy-1))
644  + 2*(h_conversion.GetBinContent(ix-2,iy+1)-h_conversion.GetBinContent(ix-2,iy-1))
645  + 2*(h_conversion.GetBinContent(ix+2,iy+1)-h_conversion.GetBinContent(ix+2,iy-1))
646  + 6*(h_conversion.GetBinContent(ix,iy+2)-h_conversion.GetBinContent(ix,iy-2))
647  + 4*(h_conversion.GetBinContent(ix-1,iy+2)-h_conversion.GetBinContent(ix-1,iy-2))
648  + 4*(h_conversion.GetBinContent(ix+1,iy+2)-h_conversion.GetBinContent(ix+1,iy-2))
649  + 1*(h_conversion.GetBinContent(ix-2,iy+2)-h_conversion.GetBinContent(ix-2,iy-2))
650  + 1*(h_conversion.GetBinContent(ix+2,iy+2)-h_conversion.GetBinContent(ix+2,iy-2)));
651  }
652  else{
653  mf::LogError("CornerFinderAlg") << "Sobel derivative not supported for neighborhoods > 2.";
654  return;
655  }
656 
657  } //end if Sobel
658 
659  else if(fDerivative_method.compare("local")==0){
660 
661  if(fDerivative_neighborhood==1){
662  h_derivative_x.SetBinContent(ix,iy,
663  (h_conversion.GetBinContent(ix+1,iy)-h_conversion.GetBinContent(ix-1,iy)));
664  h_derivative_y.SetBinContent(ix,iy,
665  (h_conversion.GetBinContent(ix,iy+1)-h_conversion.GetBinContent(ix,iy-1)));
666  }
667  else{
668  mf::LogError("CornerFinderAlg") << "Local derivative not yet supported for neighborhoods > 1.";
669  return;
670  }
671  } //end if local
672 
673  else{
674  mf::LogError("CornerFinderAlg") << "Bad derivative algorithm! " << fDerivative_method;
675  return;
676  }
677 
678  }
679  }
680 
681 
682  //this is just a double Gaussian
683  float func_blur[11][11];
684  func_blur[0][0] = 0.000000;
685  func_blur[0][1] = 0.000000;
686  func_blur[0][2] = 0.000000;
687  func_blur[0][3] = 0.000001;
688  func_blur[0][4] = 0.000002;
689  func_blur[0][5] = 0.000004;
690  func_blur[0][6] = 0.000002;
691  func_blur[0][7] = 0.000001;
692  func_blur[0][8] = 0.000000;
693  func_blur[0][9] = 0.000000;
694  func_blur[0][10] = 0.000000;
695  func_blur[1][0] = 0.000000;
696  func_blur[1][1] = 0.000000;
697  func_blur[1][2] = 0.000004;
698  func_blur[1][3] = 0.000045;
699  func_blur[1][4] = 0.000203;
700  func_blur[1][5] = 0.000335;
701  func_blur[1][6] = 0.000203;
702  func_blur[1][7] = 0.000045;
703  func_blur[1][8] = 0.000004;
704  func_blur[1][9] = 0.000000;
705  func_blur[1][10] = 0.000000;
706  func_blur[2][0] = 0.000000;
707  func_blur[2][1] = 0.000004;
708  func_blur[2][2] = 0.000123;
709  func_blur[2][3] = 0.001503;
710  func_blur[2][4] = 0.006738;
711  func_blur[2][5] = 0.011109;
712  func_blur[2][6] = 0.006738;
713  func_blur[2][7] = 0.001503;
714  func_blur[2][8] = 0.000123;
715  func_blur[2][9] = 0.000004;
716  func_blur[2][10] = 0.000000;
717  func_blur[3][0] = 0.000001;
718  func_blur[3][1] = 0.000045;
719  func_blur[3][2] = 0.001503;
720  func_blur[3][3] = 0.018316;
721  func_blur[3][4] = 0.082085;
722  func_blur[3][5] = 0.135335;
723  func_blur[3][6] = 0.082085;
724  func_blur[3][7] = 0.018316;
725  func_blur[3][8] = 0.001503;
726  func_blur[3][9] = 0.000045;
727  func_blur[3][10] = 0.000001;
728  func_blur[4][0] = 0.000002;
729  func_blur[4][1] = 0.000203;
730  func_blur[4][2] = 0.006738;
731  func_blur[4][3] = 0.082085;
732  func_blur[4][4] = 0.367879;
733  func_blur[4][5] = 0.606531;
734  func_blur[4][6] = 0.367879;
735  func_blur[4][7] = 0.082085;
736  func_blur[4][8] = 0.006738;
737  func_blur[4][9] = 0.000203;
738  func_blur[4][10] = 0.000002;
739  func_blur[5][0] = 0.000004;
740  func_blur[5][1] = 0.000335;
741  func_blur[5][2] = 0.011109;
742  func_blur[5][3] = 0.135335;
743  func_blur[5][4] = 0.606531;
744  func_blur[5][5] = 1.000000;
745  func_blur[5][6] = 0.606531;
746  func_blur[5][7] = 0.135335;
747  func_blur[5][8] = 0.011109;
748  func_blur[5][9] = 0.000335;
749  func_blur[5][10] = 0.000004;
750  func_blur[6][0] = 0.000002;
751  func_blur[6][1] = 0.000203;
752  func_blur[6][2] = 0.006738;
753  func_blur[6][3] = 0.082085;
754  func_blur[6][4] = 0.367879;
755  func_blur[6][5] = 0.606531;
756  func_blur[6][6] = 0.367879;
757  func_blur[6][7] = 0.082085;
758  func_blur[6][8] = 0.006738;
759  func_blur[6][9] = 0.000203;
760  func_blur[6][10] = 0.000002;
761  func_blur[7][0] = 0.000001;
762  func_blur[7][1] = 0.000045;
763  func_blur[7][2] = 0.001503;
764  func_blur[7][3] = 0.018316;
765  func_blur[7][4] = 0.082085;
766  func_blur[7][5] = 0.135335;
767  func_blur[7][6] = 0.082085;
768  func_blur[7][7] = 0.018316;
769  func_blur[7][8] = 0.001503;
770  func_blur[7][9] = 0.000045;
771  func_blur[7][10] = 0.000001;
772  func_blur[8][0] = 0.000000;
773  func_blur[8][1] = 0.000004;
774  func_blur[8][2] = 0.000123;
775  func_blur[8][3] = 0.001503;
776  func_blur[8][4] = 0.006738;
777  func_blur[8][5] = 0.011109;
778  func_blur[8][6] = 0.006738;
779  func_blur[8][7] = 0.001503;
780  func_blur[8][8] = 0.000123;
781  func_blur[8][9] = 0.000004;
782  func_blur[8][10] = 0.000000;
783  func_blur[9][0] = 0.000000;
784  func_blur[9][1] = 0.000000;
785  func_blur[9][2] = 0.000004;
786  func_blur[9][3] = 0.000045;
787  func_blur[9][4] = 0.000203;
788  func_blur[9][5] = 0.000335;
789  func_blur[9][6] = 0.000203;
790  func_blur[9][7] = 0.000045;
791  func_blur[9][8] = 0.000004;
792  func_blur[9][9] = 0.000000;
793  func_blur[9][10] = 0.000000;
794  func_blur[10][0] = 0.000000;
795  func_blur[10][1] = 0.000000;
796  func_blur[10][2] = 0.000000;
797  func_blur[10][3] = 0.000001;
798  func_blur[10][4] = 0.000002;
799  func_blur[10][5] = 0.000004;
800  func_blur[10][6] = 0.000002;
801  func_blur[10][7] = 0.000001;
802  func_blur[10][8] = 0.000000;
803  func_blur[10][9] = 0.000000;
804  func_blur[10][10] = 0.000000;
805 
806  double temp_integral_x = 0;
807  double temp_integral_y = 0;
808 
809  if(fDerivative_BlurNeighborhood>0){
810 
811  if(fDerivative_BlurNeighborhood>10){
812  mf::LogWarning("CornerFinderAlg") << "WARNING...BlurNeighborhoods>10 not currently allowed. Shrinking to 10.";
813  fDerivative_BlurNeighborhood=10;
814  }
815 
816  TH2F *h_clone_derivative_x = (TH2F*)h_derivative_x.Clone("h_clone_derivative_x");
817  TH2F *h_clone_derivative_y = (TH2F*)h_derivative_y.Clone("h_clone_derivative_y");
818 
819  temp_integral_x = 0;
820  temp_integral_y = 0;
821 
822  for(int ix=1; ix<=h_derivative_x.GetNbinsX(); ix++){
823  for(int iy=1; iy<=h_derivative_y.GetNbinsY(); iy++){
824 
825  temp_integral_x = 0;
826  temp_integral_y = 0;
827 
828  for(int jx=ix-fDerivative_BlurNeighborhood; jx<=ix+fDerivative_BlurNeighborhood; jx++){
829  for(int jy=iy-fDerivative_BlurNeighborhood; jy<=iy+fDerivative_BlurNeighborhood; jy++){
830  temp_integral_x += h_clone_derivative_x->GetBinContent(jx,jy)*func_blur[(ix-jx)+5][(iy-jy)+5];
831  temp_integral_y += h_clone_derivative_y->GetBinContent(jx,jy)*func_blur[(ix-jx)+5][(iy-jy)+5];
832  }
833  }
834  h_derivative_x.SetBinContent(ix,iy,temp_integral_x);
835  h_derivative_y.SetBinContent(ix,iy,temp_integral_y);
836 
837  }
838  }
839 
840  delete h_clone_derivative_x;
841  delete h_clone_derivative_y;
842 
843 
844  } //end if blur
845 
846 
847 }
848 
849 
850 //-----------------------------------------------------------------------------
851 // Corner Score
852 
853 void corner::CornerFinderAlg::create_cornerScore_histogram(TH2F const& h_derivative_x, TH2F const& h_derivative_y, TH2D & h_cornerScore){
854 
855  const int x_bins = h_derivative_x.GetNbinsX();
856  const int y_bins = h_derivative_y.GetNbinsY();
857 
858  //the structure tensor elements
859  double st_xx = 0., st_xy = 0., st_yy = 0.;
860 
861  for(int iy=1+fCornerScore_neighborhood; iy<=(y_bins-fCornerScore_neighborhood); iy++){
862  for(int ix=1+fCornerScore_neighborhood; ix<=(x_bins-fCornerScore_neighborhood); ix++){
863 
864  if(ix==1+fCornerScore_neighborhood){
865  st_xx=0.; st_xy=0.; st_yy=0.;
866 
867  for(int jx=ix-fCornerScore_neighborhood; jx<=ix+fCornerScore_neighborhood; jx++){
868  for(int jy=iy-fCornerScore_neighborhood; jy<=iy+fCornerScore_neighborhood; jy++){
869 
870  st_xx += h_derivative_x.GetBinContent(jx,jy)*h_derivative_x.GetBinContent(jx,jy);
871  st_yy += h_derivative_y.GetBinContent(jx,jy)*h_derivative_y.GetBinContent(jx,jy);
872  st_xy += h_derivative_x.GetBinContent(jx,jy)*h_derivative_y.GetBinContent(jx,jy);
873 
874  }
875  }
876  }
877 
878  // we do it this way to reduce computation time
879  else{
880  for(int jy=iy-fCornerScore_neighborhood; jy<=iy+fCornerScore_neighborhood; jy++){
881 
882  st_xx -= h_derivative_x.GetBinContent(ix-fCornerScore_neighborhood-1,jy)*h_derivative_x.GetBinContent(ix-fCornerScore_neighborhood-1,jy);
883  st_xx += h_derivative_x.GetBinContent(ix+fCornerScore_neighborhood,jy)*h_derivative_x.GetBinContent(ix+fCornerScore_neighborhood,jy);
884 
885  st_yy -= h_derivative_y.GetBinContent(ix-fCornerScore_neighborhood-1,jy)*h_derivative_y.GetBinContent(ix-fCornerScore_neighborhood-1,jy);
886  st_yy += h_derivative_y.GetBinContent(ix+fCornerScore_neighborhood,jy)*h_derivative_y.GetBinContent(ix+fCornerScore_neighborhood,jy);
887 
888  st_xy -= h_derivative_x.GetBinContent(ix-fCornerScore_neighborhood-1,jy)*h_derivative_y.GetBinContent(ix-fCornerScore_neighborhood-1,jy);
889  st_xy += h_derivative_x.GetBinContent(ix+fCornerScore_neighborhood,jy)*h_derivative_y.GetBinContent(ix+fCornerScore_neighborhood,jy);
890  }
891  }
892 
893  if( fCornerScore_algorithm.compare("Noble")==0 ) {
894  h_cornerScore.SetBinContent(ix,iy,
895  (st_xx*st_yy-st_xy*st_xy) / (st_xx+st_yy + fCornerScore_Noble_epsilon));
896  }
897  else if( fCornerScore_algorithm.compare("Harris")==0 ) {
898  h_cornerScore.SetBinContent(ix,iy,
899  (st_xx*st_yy-st_xy*st_xy) - ((st_xx+st_yy)*(st_xx+st_yy)*fCornerScore_Harris_kappa));
900  }
901  else{
902  mf::LogError("CornerFinderAlg") << "BAD CORNER ALGORITHM: " << fCornerScore_algorithm;
903  return;
904  }
905 
906  } // end for loop over x bins
907  } // end for loop over y bins
908 
909 }
910 
911 
912 //-----------------------------------------------------------------------------
913 // Max Supress
914 std::vector<recob::EndPoint2D>
915 corner::CornerFinderAlg::perform_maximum_suppression(TH2D const& h_cornerScore,
916  std::vector<geo::WireID> wireIDs,
917  geo::View_t view,
918  TH2D & h_maxSuppress,
919  int startx,
920  int starty) const
921 {
922  std::vector<recob::EndPoint2D> corner_vector;
923  const int x_bins = h_cornerScore.GetNbinsX();
924  const int y_bins = h_cornerScore.GetNbinsY();
925 
926  for(int iy=1; iy<=y_bins; iy++){
927  for(int ix=1; ix<=x_bins; ix++){
928 
929  if(h_cornerScore.GetBinContent(ix,iy) < fMaxSuppress_threshold)
930  continue;
931 
932  double temp_max = -1000;
933  double temp_center_bin = false;
934 
935  for(int jx=ix-fMaxSuppress_neighborhood; jx<=ix+fMaxSuppress_neighborhood; jx++){
936  for(int jy=iy-fMaxSuppress_neighborhood; jy<=iy+fMaxSuppress_neighborhood; jy++){
937 
938  if(h_cornerScore.GetBinContent(jx,jy) > temp_max){
939  temp_max = h_cornerScore.GetBinContent(jx,jy);
940  if(jx==ix && jy==iy) temp_center_bin=true;
941  else{ temp_center_bin=false; }
942  }
943 
944  }
945  }
946 
947  if(temp_center_bin){
948 
949  float time_tick = 0.5 * (float)((2*(iy+starty)) * fConversion_bins_per_input_y);
950  int wire_number = ( (2*(ix+startx))*fConversion_bins_per_input_x ) / 2;
951  double totalQ = 0;
952  int id = 0;
953  recob::EndPoint2D corner(time_tick,
954  wireIDs[wire_number],
955  h_cornerScore.GetBinContent(ix,iy),
956  id,
957  view,
958  totalQ);
959  corner_vector.push_back(corner);
960 
961  h_maxSuppress.SetBinContent(ix,iy,h_cornerScore.GetBinContent(ix,iy));
962  }
963 
964  }
965  }
966  return corner_vector;
967 }
968 
969 
970 /* Silly little function for doing a line integral type thing. Needs improvement. */
971 float corner::CornerFinderAlg::line_integral(TH2F const& hist, int begin_x, float begin_y, int end_x, float end_y, float threshold) const{
972 
973  int x1 = hist.GetXaxis()->FindBin( begin_x );
974  int y1 = hist.GetYaxis()->FindBin( begin_y );
975  int x2 = hist.GetXaxis()->FindBin( end_x );
976  int y2 = hist.GetYaxis()->FindBin( end_y );
977 
978  if(x1==x2 && abs(y1-y2)<1e-5)
979  return 0;
980 
981  if(x2<x1){
982  std::swap(x1, x2);
983  std::swap(y1, y2);
984  }
985 
986  float fraction = 0;
987  int bin_counter = 0;
988 
989  if(x2!=x1){
990 
991  float slope = (y2-y1)/((float)(x2-x1));
992 
993  for(int ix=x1; ix<=x2; ix++){
994 
995  int y_min,y_max;
996 
997  if(slope>=0){
998  y_min = y1 + slope*(ix-x1);
999  y_max = y1 + slope*(ix+1-x1);
1000  }
1001  else {
1002  y_max = (y1+1) + slope*(ix-x1);
1003  y_min = (y1+1) + slope*(ix+1-x1);
1004  }
1005 
1006  for(int iy=y_min; iy<=y_max; iy++){
1007  bin_counter++;
1008 
1009  if( hist.GetBinContent(ix,iy) > threshold )
1010  fraction += 1.;
1011  }
1012 
1013  }
1014  }
1015  else{
1016 
1017  auto const [y_min, y_max] = std::minmax(y1, y2);
1018  for(int iy=y_min; iy<=y_max; iy++){
1019  bin_counter++;
1020  if( hist.GetBinContent(x1,iy) > threshold)
1021  fraction += 1.;
1022  }
1023 
1024  }
1025 
1026  return fraction/bin_counter;
1027 }
1028 
1029 
1030 //-----------------------------------------------------------------------------
1031 // Do the silly little line integral score thing
1032 void corner::CornerFinderAlg::calculate_line_integral_score( TH2F const& h_wire_data,
1033  std::vector<recob::EndPoint2D> const & corner_vector,
1034  std::vector<recob::EndPoint2D> & corner_lineIntegralScore_vector,
1035  TH2F & h_lineIntegralScore) const {
1036 
1037  for(auto const i_corner : corner_vector){
1038 
1039  float score=0;
1040 
1041  for(auto const j_corner : corner_vector){
1042 
1043 
1044  if( line_integral(h_wire_data,
1045  i_corner.WireID().Wire,i_corner.DriftTime(),
1046  j_corner.WireID().Wire,j_corner.DriftTime(),
1047  fIntegral_bin_threshold) > fIntegral_fraction_threshold)
1048  {
1049  score+=1.;
1050  }
1051 
1052  }
1053 
1054  recob::EndPoint2D corner(i_corner.DriftTime(),
1055  i_corner.WireID(),
1056  score,
1057  i_corner.ID(),
1058  i_corner.View(),
1059  i_corner.Charge());
1060 
1061  corner_lineIntegralScore_vector.push_back(corner);
1062 
1063 
1064  h_lineIntegralScore.SetBinContent(h_wire_data.GetXaxis()->FindBin(i_corner.WireID().Wire),
1065  h_wire_data.GetYaxis()->FindBin(i_corner.DriftTime()),
1066  score);
1067 
1068  }
1069 }
1070 
1071 
1072 
1073 TH2F const& corner::CornerFinderAlg::GetWireDataHist(unsigned int i_plane) const {
1074  return WireData_histos.at(i_plane);
1075 }
corner::CornerFinderAlg::fDerivative_BlurFunc
std::string fDerivative_BlurFunc
Definition: CornerFinderAlg.h:74
corner::CornerFinderAlg::WireData_trimmed_histos
std::vector< std::tuple< int, TH2F, int, int > > WireData_trimmed_histos
Definition: CornerFinderAlg.h:89
corner::CornerFinderAlg::create_image_histo
void create_image_histo(TH2F const &h_wire_data, TH2F &h_conversion) const
Definition: CornerFinderAlg.cxx:548
corner::CornerFinderAlg::create_cornerScore_histogram
void create_cornerScore_histogram(TH2F const &h_derivative_x, TH2F const &h_derivative_y, TH2D &h_cornerScore)
Definition: CornerFinderAlg.cxx:853
CryostatGeo.h
Encapsulate the construction of a single cyostat.
lar_content::minmax
std::pair< float, float > minmax(const float a, const float b)
minmax
Definition: KDTreeLinkerToolsT.cc:14
corner::CornerFinderAlg::fCornerScore_algorithm
std::string fCornerScore_algorithm
Definition: CornerFinderAlg.h:77
corner::CornerFinderAlg::perform_maximum_suppression
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
Definition: CornerFinderAlg.cxx:915
geo::View_t
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
corner::CornerFinderAlg::attach_feature_points
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)
Definition: CornerFinderAlg.cxx:433
corner::CornerFinderAlg::create_smaller_histos
void create_smaller_histos(geo::Geometry const &)
Definition: CornerFinderAlg.cxx:254
score
BEGIN_PROLOG or score(default)}sbnd_crttrackmatchingalg_crID
Definition: crttrackmatchingalg_sbnd.fcl:9
corner::CornerFinderAlg::fCalDataModuleLabel
std::string fCalDataModuleLabel
Definition: CornerFinderAlg.h:62
corner::CornerFinderAlg::InitializeGeometry
void InitializeGeometry(geo::Geometry const &)
Definition: CornerFinderAlg.cxx:73
geo::GeometryCore::ChannelToWire
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition: GeometryCore.cxx:1196
compare_to_range::compare_to_range
compare_to_range(int a, int b)
Definition: CornerFinderAlg.cxx:230
geo::WireID::Wire
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
corner::CornerFinderAlg::fConversion_algorithm
std::string fConversion_algorithm
Definition: CornerFinderAlg.h:63
corner::CornerFinderAlg::GetWireDataHist
TH2F const & GetWireDataHist(unsigned int) const
Definition: CornerFinderAlg.cxx:1073
geo::GeometryCore::Ncryostats
unsigned int Ncryostats() const
Returns the number of cryostats in the detector.
Definition: GeometryCore.h:2002
corner::CornerFinderAlg::get_feature_points
void get_feature_points(std::vector< recob::EndPoint2D > &, geo::Geometry const &)
Definition: CornerFinderAlg.cxx:161
corner::CornerFinderAlg::create_derivative_histograms
void create_derivative_histograms(TH2F const &h_conversion, TH2F &h_derivative_x, TH2F &h_derivative_y)
Definition: CornerFinderAlg.cxx:608
compare_to_range::a
int a
Definition: CornerFinderAlg.cxx:247
geo::GeometryCore::Nwires
unsigned int Nwires(unsigned int p, unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wires in the specified plane.
Definition: GeometryCore.h:3345
compare_to_range::b
int b
Definition: CornerFinderAlg.cxx:248
corner::CornerFinderAlg::get_feature_points_fast
void get_feature_points_fast(std::vector< recob::EndPoint2D > &, geo::Geometry const &)
Definition: CornerFinderAlg.cxx:173
corner::CornerFinderAlg::calculate_line_integral_score
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
Definition: CornerFinderAlg.cxx:1032
a
process_name gaushit a
Definition: decode_signalprocess_icarus.fcl:37
geo::GeometryCore::IteratePlaneIDs
IteratorBox< plane_id_iterator,&GeometryCore::begin_plane_id,&GeometryCore::end_plane_id > IteratePlaneIDs() const
Enables ranged-for loops on all plane IDs of the detector.
Definition: GeometryCore.h:3130
corner::CornerFinderAlg::GrabWires
void GrabWires(std::vector< recob::Wire > const &wireVec, geo::Geometry const &)
Definition: CornerFinderAlg.cxx:99
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
geo::PlaneGeo::View
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:184
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2854
CornerFinderAlg.h
algorithm to find feature 2D points
corner::CornerFinderAlg::CornerFinderAlg
CornerFinderAlg(fhicl::ParameterSet const &pset)
Definition: CornerFinderAlg.cxx:42
corner::CornerFinderAlg::get_feature_points_LineIntegralScore
void get_feature_points_LineIntegralScore(std::vector< recob::EndPoint2D > &, geo::Geometry const &)
Definition: CornerFinderAlg.cxx:207
corner::CornerFinderAlg::WireData_histos_ProjectionY
std::vector< TH1D > WireData_histos_ProjectionY
Definition: CornerFinderAlg.h:88
corner::CornerFinderAlg::line_integral
float line_integral(TH2F const &hist, int x1, float y1, int x2, float y2, float threshold) const
Definition: CornerFinderAlg.cxx:971
geo::CryostatGeo::NTPC
unsigned int NTPC() const
Number of TPCs in this cryostat.
Definition: CryostatGeo.h:181
corner::CornerFinderAlg::WireData_histos
std::vector< TH2F > WireData_histos
Definition: CornerFinderAlg.h:86
geo::GeometryCore::Cryostat
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.cxx:321
compare_to_value
Definition: CornerFinderAlg.cxx:217
corner::CornerFinderAlg::fDerivative_method
std::string fDerivative_method
Definition: CornerFinderAlg.h:72
geo::Geometry
The geometry of one entire detector, as served by art.
Definition: Geometry.h:181
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
geo::GeometryCore::View
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Definition: GeometryCore.cxx:622
corner::CornerFinderAlg::WireData_IDs
std::vector< std::vector< geo::WireID > > WireData_IDs
Definition: CornerFinderAlg.h:90
corner::CornerFinderAlg::WireData_histos_ProjectionX
std::vector< TH1D > WireData_histos_ProjectionX
Definition: CornerFinderAlg.h:87
PlaneGeo.h
Encapsulate the construction of a single detector plane.
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
e
do i e
Definition: updateVersionList.sh:68
corner::CornerFinderAlg::attach_feature_points_LineIntegralScore
void attach_feature_points_LineIntegralScore(TH2F const &h_wire_data, std::vector< geo::WireID > const &wireIDs, geo::View_t view, std::vector< recob::EndPoint2D > &)
Definition: CornerFinderAlg.cxx:487
geo::TPCGeo::Plane
PlaneGeo const & Plane(geo::View_t view) const
Return the plane in the tpc with View_t view.
Definition: TPCGeo.cxx:263
compare_to_range
Definition: CornerFinderAlg.cxx:228
compare_to_value::b
int b
Definition: CornerFinderAlg.cxx:224
Geometry.h
art framework interface to geometry description
TPCGeo.h
Encapsulate the construction of a single detector plane.
compare_to_value::compare_to_value
compare_to_value(int b)
Definition: CornerFinderAlg.cxx:219
Generated by   doxygen 1.8.5