All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
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 27 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.

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 }
p
pdgs p
Definition: selectors.fcl:22

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 122 of file EndPointAlg.cxx.

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

Definition at line 60 of file EndPointAlg.cxx.

61 {
62  double const Norm = 1. / std::sqrt(2 * TMath::Pi() * cet::square(sigma));
63  return Norm * exp(-cet::sum_of_squares(x, y) / (2 * cet::square(sigma)));
64 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
double cluster::EndPointAlg::GaussianDerivativeX ( int  x,
int  y 
) const
private

Definition at line 68 of file EndPointAlg.cxx.

69 {
70  double const Norm = 1. / (std::sqrt(2 * TMath::Pi()) * cet::cube(fGsigma));
71  return Norm * (-x) * exp(-cet::sum_of_squares(x, y) / (2 * cet::square(fGsigma)));
72 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
cube
T cube(T side)
Definition: PointIsolationAlgStress_test.cc:52
double cluster::EndPointAlg::GaussianDerivativeY ( int  x,
int  y 
) const
private

Definition at line 76 of file EndPointAlg.cxx.

77 {
78  double const Norm = 1. / (std::sqrt(2 * TMath::Pi()) * cet::cube(fGsigma));
79  return Norm * (-y) * exp(-cet::sum_of_squares(x, y) / (2 * cet::square(fGsigma)));
80 }
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
make_calib_df.exp
def exp
Definition: make_calib_df.py:43
cube
T cube(T side)
Definition: PointIsolationAlgStress_test.cc:52
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 85 of file EndPointAlg.cxx.

86 {
87  std::ofstream bmpFile(fileName, std::ios::binary);
88  bmpFile.write("B", 1);
89  bmpFile.write("M", 1);
90  int bitsOffset = 54 + 256 * 4;
91  int size = bitsOffset + dx * dy; //header plus 256 entry LUT plus pixels
92  bmpFile.write((const char*)&size, 4);
93  int reserved = 0;
94  bmpFile.write((const char*)&reserved, 4);
95  bmpFile.write((const char*)&bitsOffset, 4);
96  int bmiSize = 40;
97  bmpFile.write((const char*)&bmiSize, 4);
98  bmpFile.write((const char*)&dx, 4);
99  bmpFile.write((const char*)&dy, 4);
100  short planes = 1;
101  bmpFile.write((const char*)&planes, 2);
102  short bitCount = 8;
103  bmpFile.write((const char*)&bitCount, 2);
104  int i, temp = 0;
105  for (i = 0; i < 6; i++)
106  bmpFile.write((const char*)&temp, 4); // zero out optional color info
107  // write a linear LUT
108  char lutEntry[4]; // blue,green,red
109  lutEntry[3] = 0; // reserved part
110  for (i = 0; i < 256; i++) {
111  lutEntry[0] = i;
112  lutEntry[1] = i + 1;
113  lutEntry[2] = i + 2;
114  bmpFile.write(lutEntry, sizeof lutEntry);
115  }
116  // write the actual pixels
117  bmpFile.write((const char*)pix, dx * dy);
118 }
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
planes
BEGIN_PROLOG Z planes
Definition: trackfinderalgorithms.fcl:26
fileName
services TFileService fileName
Definition: prodcorsika_genie_protononly_icarus_numi_nue_volDetEnclosure_tpc.fcl:18

Member Data Documentation

double cluster::EndPointAlg::fGsigma
private

Definition at line 51 of file EndPointAlg.h.

int cluster::EndPointAlg::fMaxCorners
private

Definition at line 50 of file EndPointAlg.h.

int cluster::EndPointAlg::fSaveVertexMap
private

Definition at line 54 of file EndPointAlg.h.

double cluster::EndPointAlg::fThreshold
private

Definition at line 53 of file EndPointAlg.h.

int cluster::EndPointAlg::fTimeBins
private

Definition at line 49 of file EndPointAlg.h.

int cluster::EndPointAlg::fWindow
private

Definition at line 52 of file EndPointAlg.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.5