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 }
ValidateOpDetReco.end
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
Definition: ValidateOpDetReco.py:548
tss::Hit2D::Point2D
TVector2 const & Point2D() const
Definition: TssHit2D.h:36
truth_ana.hx
hx
Definition: truth_ana.py:105
tss::Segmentation2D::mergeClusters
int mergeClusters(std::vector< tss::Cluster2D > &group, const std::vector< size_t > &idxs) const
Definition: Segmentation2D.cxx:330
result
static QCString result
Definition: fortranscanner.cpp:56614
pma::Dist2
double Dist2(const TVector2 &v1, const TVector2 &v2)
Definition: Utilities.cxx:37
tss::Cluster2D::start
const tss::Hit2D * start(void) const
Definition: SimpleClustering.h:78
tss::Segmentation2D::mergeDenseParts
void mergeDenseParts(std::vector< tss::Cluster2D > &group) const
Definition: Segmentation2D.cxx:234
tss::Cluster2D::push_back
void push_back(const tss::Hit2D *hit)
Definition: SimpleClustering.h:66
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
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:211
tss::Cluster2D::closest
const Hit2D * closest(const TVector2 &p2d, size_t &idx) const
Definition: SimpleClustering.cxx:41
truth_ana.cx
cx
Definition: truth_ana.py:338
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
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
tss::Cluster2D::hits
const std::vector< const tss::Hit2D * > & hits(void) const
Definition: SimpleClustering.h:48
pma::GetProjectionToSegment
TVector2 GetProjectionToSegment(const TVector2 &p, const TVector2 &p0, const TVector2 &p1)
Definition: Utilities.cxx:149
cet::n
std::void_t< T > n
Definition: metaprogramming.h:28
tss::Segmentation2D::Cl2InsideCl1
bool Cl2InsideCl1(tss::Cluster2D &cl1, tss::Cluster2D &cl2) const
Definition: Segmentation2D.cxx:462
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
tss::SimpleClustering::run
std::vector< tss::Cluster2D > run(const std::vector< tss::Hit2D > &inp) const
Definition: SimpleClustering.cxx:247
tss::Cluster2D::end
const tss::Hit2D * end(void) const
Definition: SimpleClustering.h:83
test.p
p
Definition: test.py:223
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
tss::Cluster2D::isEM
bool isEM(void) const
Definition: SimpleClustering.h:103
tca::seeds
std::vector< TrajPoint > seeds
Definition: DataStructs.cxx:13
tss::SimpleClustering::hitsTouching
bool hitsTouching(const tss::Hit2D &h1, const tss::Hit2D &h2) const
Definition: SimpleClustering.cxx:178
Utilities.h
Implementation of the Projection Matching Algorithm.
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
wirecell.gen.depos.center
def center(depos, point)
Definition: depos.py:117
tss::Cluster2D::min
const TVector2 min(void) const
Definition: SimpleClustering.cxx:80
truth_ana.hy
hy
Definition: truth_ana.py:106
tss::Segmentation2D::fSimpleClustering
tss::SimpleClustering fSimpleClustering
Definition: Segmentation2D.h:70
ns
QAsciiDict< Entry > ns
Definition: tclscanner.cpp:1147
tss::Cluster2D::size
size_t size(void) const
Definition: SimpleClustering.h:44
truth_ana.cy
cy
Definition: truth_ana.py:341
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.11