LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
cluster::EndPointAlg Class Reference

Algorithm to find 2D end points. More...

#include "EndPointAlg.h"

Public Member Functions

 EndPointAlg (fhicl::ParameterSet const &pset)
 
void reconfigure (fhicl::ParameterSet const &pset)
 
size_t EndPoint (const art::PtrVector< recob::Cluster > &clusIn, std::vector< recob::EndPoint2D > &vtxcol, std::vector< art::PtrVector< recob::Hit >> &vtxHitsOut, art::Event const &evt, std::string const &label) const
 

Private Member Functions

double Gaussian (int x, int y, double sigma) const
 
double GaussianDerivativeX (int x, int y) const
 
double GaussianDerivativeY (int x, int y) const
 
void VSSaveBMPFile (const char *fileName, unsigned char *pix, int dx, int dy) const
 

Private Attributes

int fTimeBins
 
int fMaxCorners
 
double fGsigma
 
int fWindow
 
double fThreshold
 
int fSaveVertexMap
 

Detailed Description

Algorithm to find 2D end points.

Definition at line 29 of file EndPointAlg.h.

Constructor & Destructor Documentation

cluster::EndPointAlg::EndPointAlg ( fhicl::ParameterSet const &  p)
explicit

The algorithm is based on: C. Harris and M. Stephens (1988). "A combined corner and edge detector". Proceedings of the 4th Alvey Vision Conference. pp. 147-151. B. Morgan (2010). "Interest Point Detection for Reconstruction in High Granularity Tracking Detectors". arXiv:1006.3012v1 [physics.ins-det]

Definition at line 48 of file EndPointAlg.cxx.

References fGsigma, fMaxCorners, fSaveVertexMap, fThreshold, fTimeBins, fWindow, and fhicl::ParameterSet::get().

49 {
50  fTimeBins = p.get<int>("TimeBins");
51  fMaxCorners = p.get<int>("MaxCorners");
52  fGsigma = p.get<double>("Gsigma");
53  fWindow = p.get<int>("Window");
54  fThreshold = p.get<double>("Threshold");
55  fSaveVertexMap = p.get<int>("SaveVertexMap");
56 }

Member Function Documentation

size_t cluster::EndPointAlg::EndPoint ( const art::PtrVector< recob::Cluster > &  clusIn,
std::vector< recob::EndPoint2D > &  vtxcol,
std::vector< art::PtrVector< recob::Hit >> &  vtxHitsOut,
art::Event const &  evt,
std::string const &  label 
) const

Definition at line 120 of file EndPointAlg.cxx.

References util::abs(), geo::PlaneID::asPlaneID(), art::PtrVector< T >::begin(), art::PtrVector< T >::clear(), art::PtrVector< T >::end(), fGsigma, fMaxCorners, fSaveVertexMap, fThreshold, fTimeBins, fWindow, Gaussian(), GaussianDerivativeX(), GaussianDerivativeY(), n, geo::PlaneID::Plane, geo::GeometryCore::Plane(), art::PtrVector< T >::push_back(), detinfo::sampling_rate(), art::PtrVector< T >::size(), geo::GeometryCore::Views(), VSSaveBMPFile(), w, geo::GeometryCore::WirePitch(), x, and y.

Referenced by cluster::EndPointModule::produce().

125 {
126 
128  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
129  auto const det_prop =
131 
132  //Point to a collection of vertices to output.
133  std::vector<art::Ptr<recob::Hit>> hit;
134 
135  art::FindManyP<recob::Hit> fmh(clusIn, evt, label);
136 
137  int flag = 0;
138  int windex = 0; // the wire index to make sure the end point finder does not
139  // fall off the edge of the hit map
140  int tindex = 0; // the time index to make sure the end point finder does not
141  // fall off the edge of the hit map
142  int n = 0; // index of window cell. There are 49 cells in the 7X7 Gaussian and
143  // Gaussian derivative windows
144  unsigned int numberwires = 0;
145  const unsigned int numbertimesamples = det_prop.ReadOutWindowSize();
146  const float TicksPerBin = numbertimesamples / fTimeBins;
147 
148  double MatrixAAsum = 0.;
149  double MatrixBBsum = 0.;
150  double MatrixCCsum = 0.;
151  std::vector<double> Cornerness2;
152  //gaussian window definitions. The cell weights are calculated here to help the algorithm's speed
153  double w[49] = {0.};
154  double wx[49] = {0.};
155  double wy[49] = {0.};
156  int ctr = 0;
157  for (int i = -3; i < 4; ++i) {
158  for (int j = 3; j > -4; --j) {
159  w[ctr] = Gaussian(i, j, fGsigma);
160  wx[ctr] = GaussianDerivativeX(i, j);
161  wy[ctr] = GaussianDerivativeY(i, j);
162  ++ctr;
163  }
164  }
165 
166  unsigned int wire = 0;
167  unsigned int wire2 = 0;
168  for (auto view : geom->Views()) {
171  hit.clear();
172  size_t cinctr = 0;
173  while (clusterIter != clusIn.end()) {
174  if ((*clusterIter)->View() == view) {
175  hit = fmh.at(cinctr);
176  } //end if cluster is in the correct view
177  clusterIter++;
178  ++cinctr;
179  }
180 
181  if (hit.size() == 0) continue;
182 
183  geo::WireID wid = hit[0]->WireID();
184  numberwires = geom->Plane(wid.asPlaneID()).Nwires();
185  mf::LogInfo("EndPointAlg") << " --- endpoints check " << numberwires << " " << numbertimesamples
186  << " " << fTimeBins;
187 
188  std::vector<std::vector<double>> MatrixAsum(numberwires);
189  std::vector<std::vector<double>> MatrixBsum(numberwires);
190  std::vector<std::vector<double>> hit_map(numberwires); //the map of hits
191 
192  //the index of the hit that corresponds to the potential corner
193  std::vector<std::vector<int>> hit_loc(numberwires);
194 
195  std::vector<std::vector<double>> Cornerness(numberwires); //the "weight" of a corner
196 
197  for (unsigned int wi = 0; wi < numberwires; ++wi) {
198  hit_map[wi].resize(fTimeBins, 0);
199  hit_loc[wi].resize(fTimeBins, -1);
200  Cornerness[wi].resize(fTimeBins, 0);
201  MatrixAsum[wi].resize(fTimeBins, 0);
202  MatrixBsum[wi].resize(fTimeBins, 0);
203  }
204  for (unsigned int i = 0; i < hit.size(); ++i) {
205  wire = hit[i]->WireID().Wire;
206  const float center = hit[i]->PeakTime(), sigma = hit[i]->RMS();
207  const int iFirstBin = int((center - 3 * sigma) / TicksPerBin),
208  iLastBin = int((center + 3 * sigma) / TicksPerBin);
209  for (int iBin = iFirstBin; iBin <= iLastBin; ++iBin) {
210  const float bin_center = iBin * TicksPerBin;
211  hit_map[wire][iBin] += Gaussian(bin_center, center, sigma);
212  }
213  }
214 
215  // Gaussian derivative convolution
216  for (unsigned int wire = 1; wire < numberwires - 1; ++wire) {
217 
218  for (int timebin = 1; timebin < fTimeBins - 1; ++timebin) {
219  MatrixAsum[wire][timebin] = 0.;
220  MatrixBsum[wire][timebin] = 0.;
221  n = 0;
222  for (int i = -3; i <= 3; ++i) {
223  windex = wire + i;
224  if (windex < 0) windex = 0;
225  // this is ok, because the line before makes sure it's not negative
226  else if ((unsigned int)windex >= numberwires)
227  windex = numberwires - 1;
228 
229  for (int j = -3; j <= 3; ++j) {
230  tindex = timebin + j;
231  if (tindex < 0)
232  tindex = 0;
233  else if (tindex >= fTimeBins)
234  tindex = fTimeBins - 1;
235 
236  MatrixAsum[wire][timebin] += wx[n] * hit_map[windex][tindex];
237  MatrixBsum[wire][timebin] += wy[n] * hit_map[windex][tindex];
238  ++n;
239  } // end loop over j
240  } // end loop over i
241  } // end loop over time bins
242  } // end loop over wires
243 
244  //calculate the cornerness of each pixel while making sure not to fall off the hit map.
245  for (unsigned int wire = 1; wire < numberwires - 1; ++wire) {
246 
247  for (int timebin = 1; timebin < fTimeBins - 1; ++timebin) {
248  MatrixAAsum = 0.;
249  MatrixBBsum = 0.;
250  MatrixCCsum = 0.;
251  //Gaussian smoothing convolution
252  n = 0;
253  for (int i = -3; i <= 3; ++i) {
254  windex = wire + i;
255  if (windex < 0) windex = 0;
256  // this is ok, because the line before makes sure it's not negative
257  else if ((unsigned int)windex >= numberwires)
258  windex = numberwires - 1;
259 
260  for (int j = -3; j <= 3; ++j) {
261  tindex = timebin + j;
262  if (tindex < 0)
263  tindex = 0;
264  else if (tindex >= fTimeBins)
265  tindex = fTimeBins - 1;
266 
267  MatrixAAsum += w[n] * pow(MatrixAsum[windex][tindex], 2);
268  MatrixBBsum += w[n] * pow(MatrixBsum[windex][tindex], 2);
269  MatrixCCsum += w[n] * MatrixAsum[windex][tindex] * MatrixBsum[windex][tindex];
270  ++n;
271  } // end loop over j
272  } // end loop over i
273 
274  if ((MatrixAAsum + MatrixBBsum) > 0)
275  Cornerness[wire][timebin] =
276  (MatrixAAsum * MatrixBBsum - pow(MatrixCCsum, 2)) / (MatrixAAsum + MatrixBBsum);
277  else
278  Cornerness[wire][timebin] = 0;
279 
280  if (Cornerness[wire][timebin] > 0) {
281  for (unsigned int i = 0; i < hit.size(); ++i) {
282  wire2 = hit[i]->WireID().Wire;
283  //make sure the end point candidate coincides with an actual hit.
284  if (wire == wire2 && std::abs(hit[i]->TimeDistanceAsRMS(
285  timebin * (numbertimesamples / fTimeBins))) < 1.) {
286  //this index keeps track of the hit number
287  hit_loc[wire][timebin] = i;
288  Cornerness2.push_back(Cornerness[wire][timebin]);
289  break;
290  }
291  } // end loop over hits
292  } // end if cornerness > 0
293  } // end loop over time bins
294  } // end wire loop
295 
296  std::sort(Cornerness2.rbegin(), Cornerness2.rend());
297 
298  for (int vertexnum = 0; vertexnum < fMaxCorners; ++vertexnum) {
299  flag = 0;
300  for (unsigned int wire = 0; wire < numberwires && flag == 0; ++wire) {
301  for (int timebin = 0; timebin < fTimeBins && flag == 0; ++timebin) {
302  if (Cornerness2.size() > (unsigned int)vertexnum)
303  if (Cornerness[wire][timebin] == Cornerness2[vertexnum] &&
304  Cornerness[wire][timebin] > 0. && hit_loc[wire][timebin] > -1) {
305  ++flag;
306 
307  //thresholding
308  if (Cornerness2.size())
309  if (Cornerness[wire][timebin] < (fThreshold * Cornerness2[0]))
310  vertexnum = fMaxCorners;
311  vHits.push_back(hit[hit_loc[wire][timebin]]);
312 
313  // get the total charge from the associated hits
314  double totalQ = 0.;
315  for (size_t vh = 0; vh < vHits.size(); ++vh)
316  totalQ += vHits[vh]->Integral();
317 
318  recob::EndPoint2D endpoint(hit[hit_loc[wire][timebin]]->PeakTime(),
319  hit[hit_loc[wire][timebin]]->WireID(),
320  Cornerness[wire][timebin],
321  vtxcol.size(),
322  view,
323  totalQ);
324  vtxcol.push_back(endpoint);
325  vtxHitsOut.push_back(vHits);
326  vHits.clear();
327 
328  // non-maximal suppression on a square window. The wire coordinate units are
329  // converted to time ticks so that the window is truly square.
330  // Note that there are 1/0.0743=13.46 time samples per 4.0 mm (wire pitch in ArgoNeuT),
331  // assuming a 1.5 mm/us drift velocity for a 500 V/cm E-field
332 
333  double drifttick = det_prop.DriftVelocity(det_prop.Efield(), det_prop.Temperature());
334  drifttick *= sampling_rate(clock_data) * 1.e-3;
335  double wirepitch = geom->WirePitch();
336  double corrfactor = drifttick / wirepitch;
337 
338  for (int wireout =
339  (int)wire -
340  (int)((fWindow * (numbertimesamples / fTimeBins) * corrfactor) + .5);
341  wireout <=
342  (int)wire + (int)((fWindow * (numbertimesamples / fTimeBins) * corrfactor) + .5);
343  ++wireout) {
344  for (int timebinout = timebin - fWindow; timebinout <= timebin + fWindow;
345  timebinout++) {
346  if (std::sqrt(pow(wire - wireout, 2) + pow(timebin - timebinout, 2)) <
347  fWindow) //circular window
348  Cornerness[wireout][timebinout] = 0;
349  }
350  }
351  }
352  }
353  }
354  }
355  Cornerness2.clear();
356  hit.clear();
357  if (clusterIter != clusIn.end()) clusterIter++;
358 
359  if ((int)wid.Plane == fSaveVertexMap) {
360  unsigned char* outPix = new unsigned char[fTimeBins * numberwires];
361  //finds the maximum cell in the map for image scaling
362  int cell = 0;
363  int pix = 0;
364  int maxCell = 0;
365  //int xmaxx, ymaxx;
366  for (int y = 0; y < fTimeBins; ++y) {
367  for (unsigned int x = 0; x < numberwires; ++x) {
368  cell = (int)(hit_map[x][y] * 1000);
369  if (cell > maxCell) { maxCell = cell; }
370  }
371  }
372  for (int y = 0; y < fTimeBins; ++y) {
373  for (unsigned int x = 0; x < numberwires; ++x) {
374  //scales the pixel weights based on the maximum cell value
375  if (maxCell > 0) pix = (int)((1500000 * hit_map[x][y]) / maxCell);
376  outPix[y * numberwires + x] = pix;
377  }
378  }
379 
380  // add 3x3 pixel squares to with the harris vertex finders to the .bmp file
381 
382  for (unsigned int ii = 0; ii < vtxcol.size(); ++ii) {
383  if (vtxcol[ii].View() == (unsigned int)view) {
384  pix = (int)(255);
385  outPix[(int)(vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) * numberwires +
386  vtxcol[ii].WireID().Wire] = pix;
387  outPix[(int)(vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) * numberwires +
388  vtxcol[ii].WireID().Wire - 1] = pix;
389  outPix[(int)(vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) * numberwires +
390  vtxcol[ii].WireID().Wire + 1] = pix;
391  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) - 1) *
392  numberwires +
393  vtxcol[ii].WireID().Wire] = pix;
394  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) + 1) *
395  numberwires +
396  vtxcol[ii].WireID().Wire] = pix;
397  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) - 1) *
398  numberwires +
399  vtxcol[ii].WireID().Wire - 1] = pix;
400  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) - 1) *
401  numberwires +
402  vtxcol[ii].WireID().Wire - 2] = pix;
403  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) + 1) *
404  numberwires +
405  vtxcol[ii].WireID().Wire + 1] = pix;
406  outPix[(int)((vtxcol[ii].DriftTime() * (fTimeBins / numbertimesamples)) + 1) *
407  numberwires +
408  vtxcol[ii].WireID().Wire + 2] = pix;
409  }
410  }
411 
412  VSSaveBMPFile(Form("harrisvertexmap_%d_%d.bmp", (*clusIn.begin())->ID(), wid.Plane),
413  outPix,
414  numberwires,
415  fTimeBins);
416  delete[] outPix;
417  }
418  } // end loop over views
419 
420  return vtxcol.size();
421 }
Float_t x
Definition: compare.C:6
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
iterator begin()
Definition: PtrVector.h:217
Float_t y
Definition: compare.C:6
constexpr auto abs(T v)
Returns the absolute value of the argument.
PlaneGeo const & Plane(PlaneID const &planeid) const
Returns the specified wire.
void VSSaveBMPFile(const char *fileName, unsigned char *pix, int dx, int dy) const
Definition: EndPointAlg.cxx:81
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:435
iterator end()
Definition: PtrVector.h:231
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Definition: geo_types.h:520
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
size_type size() const
Definition: PtrVector.h:302
double GaussianDerivativeY(int x, int y) const
Definition: EndPointAlg.cxx:73
Detector simulation of raw signals on wires.
double Gaussian(int x, int y, double sigma) const
Definition: EndPointAlg.cxx:59
double GaussianDerivativeX(int x, int y) const
Definition: EndPointAlg.cxx:66
std::set< View_t > const & Views() const
Returns a list of possible views in the detector.
Char_t n[5]
TCEvent evt
Definition: DataStructs.cxx:8
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
void clear()
Definition: PtrVector.h:533
Float_t w
Definition: plot.C:20
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
double cluster::EndPointAlg::Gaussian ( int  x,
int  y,
double  sigma 
) const
private

Definition at line 59 of file EndPointAlg.cxx.

Referenced by EndPoint().

60 {
61  double const Norm = 1. / std::sqrt(2 * TMath::Pi() * cet::square(sigma));
62  return Norm * exp(-cet::sum_of_squares(x, y) / (2 * cet::square(sigma)));
63 }
Float_t x
Definition: compare.C:6
Float_t y
Definition: compare.C:6
double cluster::EndPointAlg::GaussianDerivativeX ( int  x,
int  y 
) const
private

Definition at line 66 of file EndPointAlg.cxx.

References fGsigma, and x.

Referenced by EndPoint().

67 {
68  double const Norm = 1. / (std::sqrt(2 * TMath::Pi()) * cet::cube(fGsigma));
69  return Norm * (-x) * exp(-cet::sum_of_squares(x, y) / (2 * cet::square(fGsigma)));
70 }
Float_t x
Definition: compare.C:6
Float_t y
Definition: compare.C:6
double cluster::EndPointAlg::GaussianDerivativeY ( int  x,
int  y 
) const
private

Definition at line 73 of file EndPointAlg.cxx.

References fGsigma, and y.

Referenced by EndPoint().

74 {
75  double const Norm = 1. / (std::sqrt(2 * TMath::Pi()) * cet::cube(fGsigma));
76  return Norm * (-y) * exp(-cet::sum_of_squares(x, y) / (2 * cet::square(fGsigma)));
77 }
Float_t x
Definition: compare.C:6
Float_t y
Definition: compare.C:6
void cluster::EndPointAlg::reconfigure ( fhicl::ParameterSet const &  pset)
void cluster::EndPointAlg::VSSaveBMPFile ( const char *  fileName,
unsigned char *  pix,
int  dx,
int  dy 
) const
private

Definition at line 81 of file EndPointAlg.cxx.

References util::size().

Referenced by EndPoint().

85 {
86  std::ofstream bmpFile(fileName, std::ios::binary);
87  bmpFile.write("B", 1);
88  bmpFile.write("M", 1);
89  int bitsOffset = 54 + 256 * 4;
90  int size = bitsOffset + dx * dy; //header plus 256 entry LUT plus pixels
91  bmpFile.write((const char*)&size, 4);
92  int reserved = 0;
93  bmpFile.write((const char*)&reserved, 4);
94  bmpFile.write((const char*)&bitsOffset, 4);
95  int bmiSize = 40;
96  bmpFile.write((const char*)&bmiSize, 4);
97  bmpFile.write((const char*)&dx, 4);
98  bmpFile.write((const char*)&dy, 4);
99  short planes = 1;
100  bmpFile.write((const char*)&planes, 2);
101  short bitCount = 8;
102  bmpFile.write((const char*)&bitCount, 2);
103  int i, temp = 0;
104  for (i = 0; i < 6; i++)
105  bmpFile.write((const char*)&temp, 4); // zero out optional color info
106  // write a linear LUT
107  char lutEntry[4]; // blue,green,red
108  lutEntry[3] = 0; // reserved part
109  for (i = 0; i < 256; i++) {
110  lutEntry[0] = i;
111  lutEntry[1] = i + 1;
112  lutEntry[2] = i + 2;
113  bmpFile.write(lutEntry, sizeof lutEntry);
114  }
115  // write the actual pixels
116  bmpFile.write((const char*)pix, dx * dy);
117 }
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101

Member Data Documentation

double cluster::EndPointAlg::fGsigma
private

Definition at line 50 of file EndPointAlg.h.

Referenced by EndPoint(), EndPointAlg(), GaussianDerivativeX(), and GaussianDerivativeY().

int cluster::EndPointAlg::fMaxCorners
private

Definition at line 49 of file EndPointAlg.h.

Referenced by EndPoint(), and EndPointAlg().

int cluster::EndPointAlg::fSaveVertexMap
private

Definition at line 53 of file EndPointAlg.h.

Referenced by EndPoint(), and EndPointAlg().

double cluster::EndPointAlg::fThreshold
private

Definition at line 52 of file EndPointAlg.h.

Referenced by EndPoint(), and EndPointAlg().

int cluster::EndPointAlg::fTimeBins
private

Definition at line 48 of file EndPointAlg.h.

Referenced by EndPoint(), and EndPointAlg().

int cluster::EndPointAlg::fWindow
private

Definition at line 51 of file EndPointAlg.h.

Referenced by EndPoint(), and EndPointAlg().


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