All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Segmentation2D.cxx
Go to the documentation of this file.
1 /**
2  * @file Segmentation2D.cxx
3  *
4  * @author D.Stefan and R.Sulej
5  *
6  * @brief Split into linear clusters.
7  */
8 
9 #include "Segmentation2D.h"
10 
11 #include "fhiclcpp/ParameterSet.h"
12 #include "messagefacility/MessageLogger/MessageLogger.h"
13 #include "larreco/RecoAlg/PMAlg/Utilities.h"
14 
15 void tss::Segmentation2D::reconfigure(const fhicl::ParameterSet& p)
16 {
17  fRadiusMin = p.get< double >("RadiusMin");
18  fRadiusMax = p.get< double >("RadiusMax");
19  fMaxLineDist = p.get< double >("MaxLineDist");
20 
21  fDenseVtxRadius = p.get< double >("DenseVtxRadius");
22  fDenseMinN = p.get< unsigned int >("DenseMinNVtx");
23 
24  fDenseHitRadius = p.get< double >("DenseHitRadius");
25  fDenseMinH = p.get< unsigned int >("DenseMinNHits");
26 }
27 
28 std::vector< tss::Cluster2D > tss::Segmentation2D::run(tss::Cluster2D & inp) const
29 {
30  std::vector< tss::Cluster2D > result;
31  while (inp.size() > 1)
32  {
33  size_t idx;
34  const tss::Hit2D* hFirst = inp.outermost(idx);
35  if (!hFirst) break;
36 
37  std::vector< TVector2 > centers;
38  centers.emplace_back(hFirst->Point2D());
39 
40  while (centers.size())
41  {
42  run(inp, result, centers);
43  }
44  }
45 
46  tagDenseEnds(result);
47  mergeDenseParts(result);
48 
49  return result;
50 }
51 // ------------------------------------------------------
52 
53 void tss::Segmentation2D::run(
54  tss::Cluster2D & inp,
55  std::vector< tss::Cluster2D > & result,
56  std::vector< TVector2 > & centers) const
57 {
58  if (!centers.size()) return;
59 
60  TVector2 center(centers.front());
61  centers.erase(centers.begin());
62 
63  size_t idx;
64  const tss::Hit2D* hFirst = inp.closest(center, idx);
65 
66  const double dmax2 = 0.5 * 0.5; // does not look like startpoint selected before
67  if (!hFirst || (pma::Dist2(hFirst->Point2D(), center) > dmax2)) return;
68  center = hFirst->Point2D();
69 
70  tss::Cluster2D ring = selectRing(inp, center);
71  if (ring.size())
72  {
73  std::vector< tss::Cluster2D > seeds = fSimpleClustering.run(ring);
74  while (seeds.size())
75  {
76  size_t seedIdx = 0, hitIdx, h;
77  double d2, min_d2 = seeds.front().dist2(center, hitIdx);
78  for (size_t i = 1; i < seeds.size(); i++)
79  {
80  d2 = seeds[i].dist2(center, h);
81  if (d2 < min_d2) { min_d2 = d2; seedIdx = i; hitIdx = h; }
82  }
83 
84  tss::Cluster2D segment = buildSegment(inp, center, seeds[seedIdx][hitIdx].Point2D());
85  if (segment.size())
86  {
87  result.emplace_back(segment);
88  if (segment.size() > 1)
89  {
90  const tss::Hit2D* hEnd = segment.end();
91  if (hEnd) centers.emplace_back(hEnd->Point2D());
92  }
93  }
94 
95  seeds.erase(seeds.begin() + seedIdx);
96  }
97  }
98  else
99  {
100  result.emplace_back(tss::Cluster2D());
101  result.back().push_back(hFirst);
102  inp.release(hFirst);
103  }
104 }
105 // ------------------------------------------------------
106 
107 tss::Cluster2D tss::Segmentation2D::buildSegment(tss::Cluster2D & inp, TVector2 center, TVector2 end) const
108 {
109  const double max_d2 = fMaxLineDist * fMaxLineDist;
110  TVector2 segDir = end - center;
111 
112  double dc, min_dc = 1.0e9;
113  size_t firstIdx = 0;
114 
115  tss::Cluster2D candidates;
116  for (auto h : inp.hits())
117  {
118  TVector2 proj = pma::GetProjectionToSegment(h->Point2D(), center, end);
119  if (pma::Dist2(h->Point2D(), proj) < max_d2)
120  {
121  TVector2 hDir = h->Point2D() - center;
122  dc = hDir.Mod();
123  if ((hDir * segDir >= 0.0) || (dc < 0.1))
124  {
125  candidates.push_back(h);
126  if (dc < min_dc)
127  {
128  min_dc = dc; firstIdx = candidates.size() - 1;
129  }
130  }
131  }
132  }
133  if (candidates.size() > 1)
134  {
135  const tss::Hit2D* hFirst = candidates.hits()[firstIdx];
136  candidates.hits()[firstIdx] = candidates.hits()[0];
137  candidates.hits()[0] = hFirst;
138  candidates.sort();
139  }
140 
141  tss::Cluster2D segment;
142  if (candidates.size())
143  {
144  segment.push_back(candidates.start());
145  if (!inp.release(segment.start()))
146  {
147  mf::LogError("Segmentation2D") << "Hit not found in the input cluster.";
148  }
149 
150  size_t i = 1;
151  while ((i < candidates.size()) &&
152  fSimpleClustering.hitsTouching(segment, candidates[i]))
153  {
154  segment.push_back(candidates.hits()[i++]);
155  if (!inp.release(segment.end()))
156  {
157  mf::LogError("Segmentation2D") << "Hit not found in the input cluster.";
158  }
159  }
160  }
161 
162  return segment;
163 }
164 // ------------------------------------------------------
165 
166 tss::Cluster2D tss::Segmentation2D::selectRing(const tss::Cluster2D & inp, TVector2 center) const
167 {
168  double d2_min = fRadiusMin * fRadiusMin;
169  double d2_max = fRadiusMax * fRadiusMax;
170 
171  tss::Cluster2D ring;
172  for (size_t h = 0; h < inp.size(); h++)
173  {
174  double d2 = pma::Dist2(center, inp[h].Point2D());
175  if ((d2 >= d2_min) && (d2 <= d2_max))
176  ring.push_back(inp.hits()[h]);
177  }
178  return ring;
179 }
180 // ------------------------------------------------------
181 
182 void tss::Segmentation2D::tagDenseEnds(std::vector< tss::Cluster2D > & group) const
183 {
184  const double rad2 = fDenseVtxRadius * fDenseVtxRadius;
185 
186  for (size_t i = 0; i < group.size(); i++)
187  {
188  bool denseStart = false, denseEnd = false;
189  TVector2 start0(group[i].start()->Point2D()), end0(group[i].end()->Point2D());
190 
191  for (size_t j = 0; j < group.size(); j++)
192  {
193  if (i == j) continue;
194 
195  TVector2 start1(group[j].start()->Point2D()), end1(group[j].end()->Point2D());
196 
197  if (!group[j].isDenseStart())
198  {
199 
200  if (pma::Dist2(start1, start0) < rad2)
201  {
202  group[j].tagDenseStart(true);
203  denseStart = true;
204  }
205  if (pma::Dist2(start1, end0) < rad2)
206  {
207  group[j].tagDenseStart(true);
208  denseEnd = true;
209  }
210  }
211 
212  if (!group[j].isDenseEnd())
213  {
214  if (pma::Dist2(end1, start0) < rad2)
215  {
216  group[j].tagDenseEnd(true);
217  denseStart = true;
218  }
219  if (pma::Dist2(end1, end0) < rad2)
220  {
221  group[j].tagDenseEnd(true);
222  denseEnd = true;
223  }
224  }
225  }
226 
227  if (denseStart) group[i].tagDenseStart(true);
228  if (denseEnd) group[i].tagDenseEnd(true);
229 
230  }
231 }
232 // ------------------------------------------------------
233 
234 void tss::Segmentation2D::mergeDenseParts(std::vector< tss::Cluster2D > & group) const
235 {
236  const double rad2 = fDenseVtxRadius * fDenseVtxRadius;
237 
238  bool merged = true;
239  while (merged)
240  {
241  merged = false;
242 
243  size_t maxS = fDenseMinN, maxE = fDenseMinN;
244  std::vector< size_t > toMergeS, toMergeE;
245  int idxMaxS = -1, idxMaxE = -1;
246 
247  for (size_t i = 0; i < group.size(); i++)
248  {
249 
250  if (group[i].isEM()) continue;
251 
252  if (group[i].isDenseStart())
253  {
254 
255  size_t ns = 0;
256  std::vector< size_t > toMerge;
257  TVector2 start0(group[i].start()->Point2D());
258  for (size_t j = 0; j < group.size(); j++)
259  {
260  if (group[j].isEM()) continue;
261 
262  if (i == j)
263  {
264  if ((group[j].size() > 1) &&
265  (group[j].length2() < rad2)) ns++;
266  }
267  else
268  {
269  bool tagged = false;
270  if (pma::Dist2(start0, group[j].start()->Point2D()) < rad2)
271  {
272  ns++; toMerge.push_back(j); tagged = true;
273  }
274  if ((group[j].size() > 1) && (pma::Dist2(start0, group[j].end()->Point2D()) < rad2))
275  {
276  ns++; if (!tagged) toMerge.push_back(j);
277  }
278  }
279  }
280  if (ns > maxS) { maxS = ns; idxMaxS = i; toMergeS = toMerge; }
281  }
282 
283  if ((group[i].size() > 1) && group[i].isDenseEnd())
284  {
285  size_t ne = 0;
286  std::vector< size_t > toMerge;
287  TVector2 end0(group[i].end()->Point2D());
288  for (size_t j = 0; j < group.size(); j++)
289  {
290  if (group[j].isEM()) continue;
291 
292  if (i == j)
293  {
294  if ((group[j].size() > 1) &&
295  (group[j].length2() < rad2)) ne++;
296  }
297  else
298  {
299  bool tagged = false;
300  if (pma::Dist2(end0, group[j].start()->Point2D()) < rad2)
301  {
302  ne++; toMerge.push_back(j); tagged = true;
303  }
304  if ((group[j].size() > 1) && (pma::Dist2(end0, group[j].end()->Point2D()) < rad2))
305  {
306  ne++; if (!tagged) toMerge.push_back(j);
307  }
308  }
309  }
310  if (ne > maxE) { maxE = ne; idxMaxE = i; toMergeE = toMerge; }
311  }
312  }
313 
314  int idx = idxMaxS;
315  std::vector< size_t > toMergeIdxs = toMergeS;
316  if (idxMaxE > idx) { idx = idxMaxE; toMergeIdxs = toMergeE; }
317  if (idx > -1)
318  {
319  toMergeIdxs.push_back(idx);
320  idx = mergeClusters(group, toMergeIdxs);
321 
322  if (idx > -1) group[idx].tagEM(true);
323 
324  merged = true;
325  }
326  }
327 }
328 // ------------------------------------------------------
329 
330 int tss::Segmentation2D::mergeClusters(std::vector< tss::Cluster2D > & group, const std::vector< size_t > & idxs) const
331 {
332  if (idxs.size() < 2) return 0;
333 
334  for (const auto i : idxs)
335  {
336  if (i < group.size()) group[i].setTag(true);
337  else
338  {
339  mf::LogError("Segmentation2D") << "Merging index out of group range.";
340  return -1;
341  }
342  }
343 
344  size_t k = idxs.front(), i = idxs.front() + 1;
345  while (i < group.size())
346  {
347  if (group[i].isTagged())
348  {
349  group[k].merge(group[i]);
350  group.erase(group.begin() + i);
351  }
352  else ++i;
353  }
354  group[k].setTag(false);
355 
356  return k;
357 }
358 // ------------------------------------------------------
359 
360 void tss::Segmentation2D::splitHits(
361  const std::vector< tss::Cluster2D > & inp,
362  std::vector< const tss::Hit2D* > & trackHits,
363  std::vector< const tss::Hit2D* > & emHits) const
364 {
365 
366  trackHits.clear();
367  emHits.clear();
368 
369  for (const auto & cx : inp)
370  {
371  if (!cx.size()) continue;
372 
373  if (cx.isEM())
374  {
375  for (auto h : cx.hits()) emHits.push_back(h);
376  }
377  else
378  {
379  for (auto h : cx.hits()) trackHits.push_back(h);
380  }
381  }
382 
383 }
384 // ------------------------------------------------------
385 
386 void tss::Segmentation2D::splitHitsNaive(
387  const tss::Cluster2D & inp,
388  std::vector< const tss::Hit2D* > & trackHits,
389  std::vector< const tss::Hit2D* > & emHits) const
390 {
391  const double rad2 = fDenseHitRadius * fDenseHitRadius;
392 
393  trackHits.clear();
394  emHits.clear();
395 
396  for (const auto hx : inp.hits())
397  {
398  size_t n = 0;
399  for (const auto hy : inp.hits())
400  {
401  if (hx->Hit2DPtr() == hy->Hit2DPtr()) continue;
402 
403  if (pma::Dist2(hx->Point2D(), hy->Point2D()) < rad2) n++;
404  }
405 
406  if (n > fDenseMinH)
407  {
408  emHits.push_back(hx);
409  }
410  else
411  {
412  trackHits.push_back(hx);
413  }
414  }
415 }
416 // ------------------------------------------------------
417 
418 void tss::Segmentation2D::splitHitsNaive(
419  const std::vector< tss::Cluster2D > & inp,
420  std::vector< const tss::Hit2D* > & trackHits,
421  std::vector< const tss::Hit2D* > & emHits) const
422 {
423  const double rad2 = fDenseHitRadius * fDenseHitRadius;
424 
425  trackHits.clear();
426  emHits.clear();
427 
428  for (const auto & cx : inp)
429  {
430  if (!cx.size()) continue;
431 
432  for (const auto hx : cx.hits())
433  {
434  size_t n = 0;
435  for (const auto & cy : inp)
436  {
437  if (!cy.size()) continue;
438 
439  for (const auto hy : cy.hits())
440  {
441  if (hx->Hit2DPtr() == hy->Hit2DPtr()) continue;
442 
443  if (pma::Dist2(hx->Point2D(), hy->Point2D()) < rad2) n++;
444  }
445  }
446 
447  if (n > fDenseMinH)
448  {
449  emHits.push_back(hx);
450  }
451  else
452  {
453  trackHits.push_back(hx);
454  }
455  }
456  }
457 }
458 // ------------------------------------------------------
459 
460 // ------------------------------------------------------
461 
462 bool tss::Segmentation2D::Cl2InsideCl1(tss::Cluster2D& cl1, tss::Cluster2D& cl2) const
463 {
464  bool clover = false; bool clunder = false;
465  bool clleft = false; bool clright = false;
466 
467  if (cl2.isEM()) return false;
468 
469  float shift = 5; // shift!
470 
471  TVector2 point((cl2.max()).X(), (cl2.min()).Y());
472  float width = (cl2.max()).X() - (cl2.min()).X();
473  float height = (cl2.max()).Y() - (cl2.min()).Y();
474 
475  for (unsigned int h = 0; h < cl1.size(); h++)
476  {
477  float wire = cl1[h].Point2D().X(); float drift = cl1[h].Point2D().Y();
478 
479  if ( (wire <= (point.X() + shift)) && (wire >= (point.X() - width - shift)) )
480  {
481  if (drift < point.Y())
482  {
483  clover = true;
484  }
485  else if (drift > (point.Y() + height))
486  {
487  clunder = true;
488  }
489 
490  if (wire > point.X())
491  {
492  clleft = true;
493  }
494  else if (wire < (point.X() - width))
495  {
496  clright = true;
497  }
498  }
499  }
500 
501  if (clover && clunder && clleft && clright) return true;
502  else return false;
503 }
tss::Hit2D::Point2D
TVector2 const & Point2D() const
Definition: TssHit2D.h:36
tss::Segmentation2D::mergeClusters
int mergeClusters(std::vector< tss::Cluster2D > &group, const std::vector< size_t > &idxs) const
Definition: Segmentation2D.cxx:330
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:40
tss::Segmentation2D::mergeDenseParts
void mergeDenseParts(std::vector< tss::Cluster2D > &group) const
Definition: Segmentation2D.cxx:234
p
pdgs p
Definition: selectors.fcl:22
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
tss::Segmentation2D::splitHitsNaive
void splitHitsNaive(const tss::Cluster2D &inp, std::vector< const tss::Hit2D * > &trackHits, std::vector< const tss::Hit2D * > &emHits) const
Definition: Segmentation2D.cxx:386
Segmentation2D.h
Split into linear clusters.
tss::Cluster2D::outermost
const Hit2D * outermost(size_t &idx) const
Definition: SimpleClustering.cxx:57
X
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
Definition: find_global_symbol.sh:107
tss::Cluster2D::closest
const Hit2D * closest(const TVector2 &p2d, size_t &idx) const
Definition: SimpleClustering.cxx:41
shift
shift
Definition: fcl_checks.sh:26
tss::Segmentation2D::selectRing
tss::Cluster2D selectRing(const tss::Cluster2D &inp, TVector2 center) const
Definition: Segmentation2D.cxx:166
tss::Segmentation2D::reconfigure
void reconfigure(const fhicl::ParameterSet &p)
Definition: Segmentation2D.cxx:15
h
while getopts h
Definition: getStatistics.sh:13
tss::Cluster2D::hits
const std::vector< const tss::Hit2D * > & hits(void) const
Definition: larreco/larreco/ClusterFinder/TrackShowerSplitter/Segmentation2D/SimpleClustering.h:48
group
then echo group
Definition: sbndpoms_genfclwithrunnumber_maker.sh:220
pma::GetProjectionToSegment
TVector2 GetProjectionToSegment(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: larreco/larreco/RecoAlg/PMAlg/Utilities.cxx:152
tss::Segmentation2D::Cl2InsideCl1
bool Cl2InsideCl1(tss::Cluster2D &cl1, tss::Cluster2D &cl2) const
Definition: Segmentation2D.cxx:462
icarus::ns::util::end
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
tss::Segmentation2D::run
std::vector< tss::Cluster2D > run(tss::Cluster2D &inp) const
Definition: Segmentation2D.cxx:28
tss::Segmentation2D::buildSegment
tss::Cluster2D buildSegment(tss::Cluster2D &inp, TVector2 center, TVector2 end) const
Definition: Segmentation2D.cxx:107
tca::seeds
std::vector< TrajPoint > seeds
Definition: DataStructs.cxx:14
tss::Segmentation2D::tagDenseEnds
void tagDenseEnds(std::vector< tss::Cluster2D > &group) const
Definition: Segmentation2D.cxx:182
tss::Cluster2D::release
bool release(const tss::Hit2D *hit)
Definition: SimpleClustering.cxx:31
tss::Cluster2D::max
const TVector2 max(void) const
Definition: SimpleClustering.cxx:98
tss::Cluster2D::min
const TVector2 min(void) const
Definition: SimpleClustering.cxx:80
k
pdgs k
Definition: selectors.fcl:22
tss::Segmentation2D::splitHits
void splitHits(const std::vector< tss::Cluster2D > &inp, std::vector< const tss::Hit2D * > &trackHits, std::vector< const tss::Hit2D * > &emHits) const
Definition: Segmentation2D.cxx:360
Generated by   doxygen 1.8.5