SkeletonAlg.cxx
Go to the documentation of this file.
1 /**
2  * @file SkeletonAlg.cxx
3  *
4  * @brief This algorithm accepts a list of 3D space points and attempts to find the "medial skeleton"
5  * describing the "best" path through them. Here, "medial" is defined by the best time match
6  * between associated 2D hits along the "best" directions.
7  *
8  */
9 
10 // Framework Includes
11 #include "fhiclcpp/ParameterSet.h"
12 
14 
15 // LArSoft includes
18 #include "lardata/RecoObjects/Cluster3D.h"
19 
20 // std includes
21 #include <string>
22 #include <functional>
23 #include <iostream>
24 #include <memory>
25 
26 //------------------------------------------------------------------------------------------------------------------------------------------
27 
28 namespace lar_cluster3d
29 {
30 
32 {
33  reconfigure(pset);
34 }
35 
36 //------------------------------------------------------------------------------------------------------------------------------------------
37 
39 {
40 }
41 
42 //------------------------------------------------------------------------------------------------------------------------------------------
43 
45 {
46  m_minimumDeltaTicks = pset.get<double>("MinimumDeltaTicks", 0.05);
47  m_maximumDeltaTicks = pset.get<double>("MaximumDeltaTicks", 10.0 );
48 }
49 
50 double SkeletonAlg::FindFirstAndLastWires(std::vector<const reco::ClusterHit3D*>& hitVec,
51  int viewToCheck,
52  int referenceWire,
53  double referenceTicks,
54  int& firstWire,
55  int& lastWire) const
56 {
57  // In the simple case the first and last wires are simply the front and back of the input vector
58  firstWire = hitVec.front()->getHits()[viewToCheck]->getHit().WireID().Wire;
59  lastWire = hitVec.back()->getHits()[viewToCheck]->getHit().WireID().Wire;
60 
61  double maxDeltaTicks = referenceTicks - hitVec.front()->getHits()[viewToCheck]->getTimeTicks();
62  double minDeltaTicks = referenceTicks - hitVec.back()->getHits()[viewToCheck]->getTimeTicks();
63 
64  if (minDeltaTicks > maxDeltaTicks) std::swap(maxDeltaTicks, minDeltaTicks);
65 
66  double bestDeltaTicks = 1000000.;
67 
68  // Can't have a gap if only one element
69  if (hitVec.size() > 1)
70  {
71  // The issue is that there may be a gap in wires which we need to find.
72  // Reset the last wire
73  lastWire = firstWire;
74 
75  // Keep track of all the deltas
76  double nextBestDeltaTicks = bestDeltaTicks;
77 
78  for(const auto& hitPair : hitVec)
79  {
80  int curWire = hitPair->getHits()[viewToCheck]->getHit().WireID().Wire;
81  double deltaTicks = referenceTicks - hitPair->getHits()[viewToCheck]->getTimeTicks();
82 
83  maxDeltaTicks = std::max(maxDeltaTicks, deltaTicks);
84  minDeltaTicks = std::min(minDeltaTicks, deltaTicks);
85 
86  if (bestDeltaTicks > fabs(deltaTicks))
87  {
88  // By definition, the "next best" is now the current best
89  nextBestDeltaTicks = bestDeltaTicks;
90  bestDeltaTicks = fabs(deltaTicks);
91  }
92  else if (nextBestDeltaTicks > fabs(deltaTicks))
93  {
94  nextBestDeltaTicks = fabs(deltaTicks);
95  }
96 
97  // if gap detected, take action depending on where the input wire is
98  // But remember we need to be willing to accept a 1 wire gap... for efficiency
99  if (fabs(curWire - lastWire) > 2000)
100  {
101  if (referenceWire <= lastWire) break;
102 
103  // Stepped over gap, reset everything
104  firstWire = curWire;
105  maxDeltaTicks = deltaTicks;
106  minDeltaTicks = deltaTicks;
107  bestDeltaTicks = fabs(deltaTicks);
108  nextBestDeltaTicks = bestDeltaTicks;
109  }
110 
111  lastWire = curWire;
112  }
113 
114  bestDeltaTicks = nextBestDeltaTicks;
115  }
116 
117  bestDeltaTicks = std::max(std::min(bestDeltaTicks,m_maximumDeltaTicks), m_minimumDeltaTicks);
118 
119  if (minDeltaTicks * maxDeltaTicks > 0. && bestDeltaTicks > 30.) bestDeltaTicks = 0.;
120 
121  return bestDeltaTicks;
122 }
123 
124 class OrderHitsAlongWire
125 {
126 public:
127  OrderHitsAlongWire(int view = 0) : m_view(view) {}
128 
130  {
131  int viewToCheck = (m_view + 1) % 3;
132 
133  return left->getHits()[viewToCheck]->getHit().WireID().Wire < right->getHits()[viewToCheck]->getHit().WireID().Wire;
134  }
135 private:
136  int m_view;
137 };
138 
139 struct OrderBestViews
140 {
141  bool operator()(const std::pair<size_t,size_t>& left, const std::pair<size_t,size_t>& right)
142  {
143  return left.second < right.second;
144  }
145 
146 };
147 
149 {
150  // Our mission is to try to find the medial skeletion of the input list of hits
151  // We define that as the set of hit pairs where the pairs share the same hit in a given direction
152  // and the selected medial hit is equal distance from the edges.
153  // The first step in trying to do this is to build a map which relates a give 2D hit to an ordered list of hitpairs using it
154  typedef std::map<const reco::ClusterHit2D*, std::vector<const reco::ClusterHit3D*> > Hit2DtoHit3DMapU;
155 
156  Hit2DtoHit3DMapU hit2DToHit3DMap[3];
157 
158  for(const auto& hitPair : hitPairList)
159  {
160  // Don't consider points "rejected" earlier
161  if (hitPair->bitsAreSet(reco::ClusterHit3D::REJECTEDHIT)) continue;
162 
163  for(const auto& hit2D : hitPair->getHits())
164  {
165  size_t view = hit2D->getHit().View();
166 
167  hit2DToHit3DMap[view][hit2D].push_back(hitPair);
168  }
169  }
170 
171  // Do an explicit loop through to sort?
172  for(size_t idx = 0; idx < 3; idx++)
173  {
174  for(auto& mapPair : hit2DToHit3DMap[idx])
175  {
176  size_t numHitPairs = mapPair.second.size();
177 
178  if (numHitPairs > 1) std::sort(mapPair.second.begin(), mapPair.second.end(), OrderHitsAlongWire(idx));
179  }
180  }
181 
182  // Keep a count of the number of skeleton points to be returned
183  int nSkeletonPoints(0);
184 
185  // The idea is go through all the hits again and determine if they could be "skeleton" elements
186  for(const auto& hitPair : hitPairList)
187  {
188  // Don't consider points "rejected" earlier
189  if (hitPair->bitsAreSet(reco::ClusterHit3D::REJECTEDHIT)) continue;
190 
191  // If a hit pair we skip for now
192  if (hitPair->getHits().size() < 3) continue;
193 
194  // Hopefully I am not confusing myself here.
195  // The goal is to know, for a given 3D hit, how many other 3D hits share the 2D hits that it is made of
196  // We want this to be a bit more precise in that we don't count other 3D hits if there is a gap between
197  // the current hit and the one we are comparing.
198  // How do we do this?
199  // We consider each 2D hit comprising the 3D hit in turn... for each 2D hit we have a list of the 3D
200  // hits which share this hit. Note that for the 3D hits to be unique, there must be an equal number of 2D
201  // hits from the other two views. So, to count "wires" to the boundary, we can simply pick one of the other views
202 
203  // so the idea is to through each of the views to determine
204  // a) the difference in distance to the last hit on this wire (in each direction)
205  // b) the number of 3D hits using the 2D hit in the given 3D hit
206  size_t numHitPairs[3] = {0, 0, 0}; // This keeps track of the number of 3D hits a given 2D hit is associated with
207  int deltaWires[3] = {0, 0, 0};
208  double viewDeltaT[3] = {0., 0., 0.};
209  double bestDeltaTicks[3] = {0., 0., 0.};
210  int wireNumByView[3] = {int(hitPair->getHits()[0]->getHit().WireID().Wire),
211  int(hitPair->getHits()[1]->getHit().WireID().Wire),
212  int(hitPair->getHits()[2]->getHit().WireID().Wire)};
213 
214  size_t bestViewIdx(0);
215 
216  // Initiate the loop over views
217  for(size_t viewIdx = 0; viewIdx < 3; viewIdx++)
218  {
219  const reco::ClusterHit2D* hit2D = hitPair->getHits()[viewIdx];
220  double hit2DTimeTicks = hit2D->getTimeTicks();
221 
222  std::vector<const reco::ClusterHit3D*>& hitVec(hit2DToHit3DMap[viewIdx][hit2D]);
223 
224  numHitPairs[viewIdx] = hitVec.size();
225 
226  if (numHitPairs[viewIdx] > 1)
227  {
228  int viewToCheck = (viewIdx + 1) % 3;
229  int firstWire;
230  int lastWire;
231 
232  bestDeltaTicks[viewIdx] = FindFirstAndLastWires(hitVec,
233  viewToCheck,
234  wireNumByView[viewToCheck],
235  hit2DTimeTicks,
236  firstWire,
237  lastWire);
238 
239  int deltaFirst = wireNumByView[viewToCheck] - firstWire;
240  int deltaLast = wireNumByView[viewToCheck] - lastWire;
241 
242  // If either distance to the first or last wire is zero then we are an edge point
243  if (deltaFirst == 0 || deltaLast == 0) hitPair->setStatusBit(reco::ClusterHit3D::EDGEHIT);
244 
245  deltaWires[viewIdx] = deltaFirst + deltaLast;
246  numHitPairs[viewIdx] = lastWire - firstWire + 1;
247  viewDeltaT[viewIdx] = fabs(hit2DTimeTicks - hitPair->getHits()[viewToCheck]->getTimeTicks());
248  }
249  // Otherwise, by definition, it is both a skeleton point and an edge point
250  else hitPair->setStatusBit(reco::ClusterHit3D::SKELETONHIT | reco::ClusterHit3D::EDGEHIT);
251 
252  if (numHitPairs[viewIdx] < numHitPairs[bestViewIdx])
253  {
254  bestViewIdx = viewIdx;
255  }
256  }
257 
258  // Make sure a real index was found (which should always happen)
259  if (bestViewIdx < 3 && numHitPairs[bestViewIdx] > 0)
260  {
261  // Make a sliding value for the width of the core region
262  int maxBestWires = 0.05 * numHitPairs[bestViewIdx] + 1;
263 
264  maxBestWires = 5000;
265 
266  // Check condition for a "skeleton" point
267  if (fabs(deltaWires[bestViewIdx]) < maxBestWires + 1)
268  {
269 
270  // If exactly in the middle then we are done
271  // Set a bit in the ClusterHit3D word to signify
272 // if (fabs(deltaWires[bestViewIdx]) < maxBestWires || numHitPairs[bestViewIdx] < 3)
273 // hitPair->setStatusBit(reco::ClusterHit3D::SKELETONHIT);
274 
275  // Otherwise we can try to look at the other view
276 // else
277  {
278  // Find the next best view
279  int nextBestIdx = (bestViewIdx + 1) % 3;
280 
281  for(size_t idx = 0; idx < 3; idx++)
282  {
283  if (idx == bestViewIdx) continue;
284 
285  if (numHitPairs[idx] < numHitPairs[nextBestIdx]) nextBestIdx = idx;
286  }
287 
288  if (nextBestIdx > 2)
289  {
290  std::cout << "***** invalid next best view: " << nextBestIdx << " *******" << std::endl;
291  continue;
292  }
293 
294  if (viewDeltaT[bestViewIdx] < 1.01*bestDeltaTicks[bestViewIdx] && viewDeltaT[nextBestIdx] < 6.01*bestDeltaTicks[nextBestIdx])
295  hitPair->setStatusBit(reco::ClusterHit3D::SKELETONHIT);
296  }
297  }
298  }
299 
300  // We want to keep count of "pure" skeleton points only
301  if (hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONHIT) && !hitPair->bitsAreSet(reco::ClusterHit3D::EDGEHIT)) nSkeletonPoints++;
302  }
303 
304  return nSkeletonPoints;
305 }
306 
307 void SkeletonAlg::GetSkeletonHits(const reco::HitPairListPtr& inputHitList, reco::HitPairListPtr& skeletonHitList) const
308 {
309  for(const auto& hit3D : inputHitList) if (hit3D->bitsAreSet(reco::ClusterHit3D::SKELETONHIT)) skeletonHitList.emplace_back(hit3D);
310 
311  return;
312 }
313 
315 {
316  // NOTE: This method assumes the list being given to it is comprised of skeleton hits
317  // YMMV if you send in a complete hit collection!
318 
319  // Define a mapping between 2D hits and the 3D hits they make
320  typedef std::map<const reco::ClusterHit2D*, std::vector<const reco::ClusterHit3D*> > Hit2DtoHit3DMap;
321 
322  // We want to keep track of this mapping by view
323  Hit2DtoHit3DMap hit2DToHit3DMap[3];
324 
325  // Keep count of the number of skeleton hits selected
326  unsigned int nSkeletonHits(0);
327 
328  // Execute the first loop through the hits to build the map
329  for(const auto& hitPair : skeletonHitList)
330  {
331  // Don't consider points "rejected" earlier
332  if (hitPair->bitsAreSet(reco::ClusterHit3D::REJECTEDHIT)) continue;
333 
334  // Count only those skeleton hits which have not been averaged
335  if (!hitPair->bitsAreSet(reco::ClusterHit3D::SKELETONPOSAVE)) nSkeletonHits++;
336 
337  for(const auto& hit2D : hitPair->getHits())
338  {
339  size_t view = hit2D->getHit().View();
340 
341  hit2DToHit3DMap[view][hit2D].push_back(hitPair);
342  }
343  }
344 
345  // Exit early if no skeleton hits to process
346  if (!nSkeletonHits) return;
347 
348  // The list of 3D hits associated to any given 2D hit is most useful to us if it is ordered
349  // So this will loop through entries in the map and do the ordering
350  for(size_t idx = 0; idx < 3; idx++)
351  {
352  for(auto& mapPair : hit2DToHit3DMap[idx])
353  {
354  size_t numHitPairs = mapPair.second.size();
355 
356  if (numHitPairs > 1) std::sort(mapPair.second.begin(), mapPair.second.end(), OrderHitsAlongWire(idx));
357  }
358  }
359 
360  // Ok, so the basic strategy is not entirely different from that used to build the skeleton hits in the first
361  // place. The idea is to loop through all skeleton hits and then determine the average position for the hit
362  // based on the wire direction with the fewest hits to the edge.
363  // Currently this is a pretty simple minded approach and it would appear that some effort could be spent
364  // here to improve this.
365  // NOTE: the averaging being performed is ONLY in the Y-Z plane since this is where the hit ambiguity
366  // issue arises.
367  reco::HitPairListPtr::iterator hitPairItr = skeletonHitList.begin();
368 
369  for(int bestViewVecIdx = 0; bestViewVecIdx < 2; bestViewVecIdx++)
370  {
371  std::list<reco::ClusterHit3D> tempHitPairList;
372  reco::HitPairListPtr tempHitPairListPtr;
373 
374  std::map<const reco::ClusterHit3D*, const reco::ClusterHit3D*> hit3DToHit3DMap;
375 
376  while(hitPairItr != skeletonHitList.end())
377  {
378  const reco::ClusterHit3D* hit3D = *hitPairItr++;
379 
380  std::vector<std::pair<size_t,size_t> > bestViewVec;
381 
382  for(const auto& hit2D : hit3D->getHits())
383  {
384  bestViewVec.push_back(std::pair<size_t,size_t>(hit2D->getHit().View(),hit2DToHit3DMap[hit2D->getHit().View()][hit2D].size()));
385  }
386 
387  std::sort(bestViewVec.begin(), bestViewVec.end(), OrderBestViews());
388 
389  size_t bestViewIdx = bestViewVec[bestViewVecIdx].first;
390  size_t bestViewCnt = bestViewVec[bestViewVecIdx].second;
391 
392  if (bestViewCnt > 5) continue;
393 
394  std::vector<const reco::ClusterHit3D*>& hit3DVec = hit2DToHit3DMap[bestViewIdx][hit3D->getHits()[bestViewIdx]];
395 
396  double avePosition[3] = {hit3D->getPosition()[0],0.,0.};
397 
398  for(const auto& tempHit3D : hit3DVec)
399  {
400  avePosition[1] += tempHit3D->getPosition()[1];
401  avePosition[2] += tempHit3D->getPosition()[2];
402  }
403 
404  avePosition[1] *= 1./double(hit3DVec.size());
405  avePosition[2] *= 1./double(hit3DVec.size());
406 
407  tempHitPairList.emplace_back(reco::ClusterHit3D(hit3D->getID(),
408  hit3D->getStatusBits(),
409  avePosition,
410  hit3D->getTotalCharge(),
411  hit3D->getAvePeakTime(),
412  hit3D->getDeltaPeakTime(),
413  hit3D->getSigmaPeakTime(),
414  hit3D->getDocaToAxis(),
415  hit3D->getArclenToPoca(),
416  hit3D->getOverlapFraction(),
417  hit3D->getHits()));
418 
419  tempHitPairListPtr.push_back(&tempHitPairList.back());
420 
421  hit3DToHit3DMap[tempHitPairListPtr.back()] = hit3D;
422  }
423 
424  for(const auto& pair : hit3DToHit3DMap)
425  {
426  pair.second->setPosition(pair.first->getPosition());
427  pair.second->setStatusBit(reco::ClusterHit3D::SKELETONPOSAVE);
428  }
429  }
430 
431  return;
432 }
433 
434 
435 } // namespace lar_cluster3d
intermediate_table::iterator iterator
float getTimeTicks() const
Definition: Cluster3D.h:76
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:158
Hit is an "edge" hit.
Definition: Cluster3D.h:101
float getTotalCharge() const
Definition: Cluster3D.h:162
Hit has been rejected for any reason.
Definition: Cluster3D.h:99
float getSigmaPeakTime() const
Definition: Cluster3D.h:165
float getOverlapFraction() const
Definition: Cluster3D.h:167
unsigned int getStatusBits() const
Definition: Cluster3D.h:157
void GetSkeletonHits(const reco::HitPairListPtr &inputHitList, reco::HitPairListPtr &skeletonHitList) const
Return the skeleton hits from the input list.
float getDocaToAxis() const
Definition: Cluster3D.h:169
SkeletonAlg(fhicl::ParameterSet const &pset)
Constructor.
Definition: SkeletonAlg.cxx:31
art framework interface to geometry description
void reconfigure(fhicl::ParameterSet const &pset)
a handler for the case where the algorithm control parameters are to be reset
Definition: SkeletonAlg.cxx:44
float getAvePeakTime() const
Definition: Cluster3D.h:163
bool operator()(const std::pair< size_t, size_t > &left, const std::pair< size_t, size_t > &right)
void swap(Handle< T > &a, Handle< T > &b)
int FindMedialSkeleton(reco::HitPairListPtr &hitPairList) const
This is intended to find the medial skeleton given a list of input hit pairs.
T get(std::string const &key) const
Definition: ParameterSet.h:271
void AverageSkeletonPositions(reco::HitPairListPtr &skeletonHitList) const
Modifies the position of input skeleton hits by averaging along the "best" wire direction.
size_t getID() const
Definition: Cluster3D.h:156
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right)
static int max(int a, int b)
double FindFirstAndLastWires(std::vector< const reco::ClusterHit3D * > &hitVec, int viewToCheck, int referenceWire, double referenceTicks, int &firstWire, int &lastWire) const
A function to find the bounding wires in a given view.
Definition: SkeletonAlg.cxx:50
Declaration of signal hit object.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
virtual ~SkeletonAlg()
Destructor.
Definition: SkeletonAlg.cxx:38
float getArclenToPoca() const
Definition: Cluster3D.h:170
float getDeltaPeakTime() const
Definition: Cluster3D.h:164
Skeleton hit position averaged.
Definition: Cluster3D.h:106
Hit is a "skeleton" hit.
Definition: Cluster3D.h:100
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:171
void setPosition(const Eigen::Vector3f &pos) const
Definition: Cluster3D.h:186
QTextStream & endl(QTextStream &s)