xmsgeom  1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
GmPtSearch.cpp
Go to the documentation of this file.
1 //------------------------------------------------------------------------------
7 //------------------------------------------------------------------------------
8 
9 //----- Included files ---------------------------------------------------------
10 
11 // 1. Precompiled header
12 
13 // 2. My header
15 
16 // 3. Standard Library Headers
17 
18 // 4. External Library Headers
19 #pragma warning(push)
20 #pragma warning(disable : 4512) // boost code: no assignment operator
21 #include <boost/geometry.hpp>
22 #include <boost/geometry/geometries/register/point.hpp>
23 #include <boost/geometry/index/rtree.hpp>
24 #include <boost/iterator/counting_iterator.hpp>
25 #pragma warning(pop)
26 
27 // 5. Shared Headers
28 #include <xmscore/points/pt.h>
29 #include <xmscore/stl/utility.h>
30 #include <xmscore/misc/XmError.h>
31 #include <xmscore/misc/DynBitset.h>
32 #include <xmscore/misc/XmConst.h>
33 
34 // 6. Non-shared Headers
35 
36 //----- Forward declarations ---------------------------------------------------
37 
38 //----- External globals -------------------------------------------------------
39 
40 //----- Namespace declaration --------------------------------------------------
41 namespace
42 {
43 class bPt
44 {
45 public:
46  bPt()
47  : x(0)
48  , y(0)
49  , z(0)
50  {
51  }
52  bPt(double a_x, double a_y, double a_z)
53  : x(a_x)
54  , y(a_y)
55  , z(a_z)
56  {
57  }
58  bPt(const xms::Pt3d& a_)
59  : x(a_.x)
60  , y(a_.y)
61  , z(a_.z)
62  {
63  }
64  bPt& operator=(const xms::Pt3d a_)
65  {
66  x = a_.x;
67  y = a_.y;
68  z = a_.z;
69  return (*this);
70  }
71 
72  double x, y, z;
73 };
74 } // unnamed namespace
78 BOOST_GEOMETRY_REGISTER_POINT_3D(bPt, double, cs::cartesian, x, y, z);
79 using namespace xms;
80 namespace xms
81 {
82 namespace bg = boost::geometry;
83 namespace bgi = bg::index;
84 
86 typedef bg::model::box<bPt> box;
88 typedef size_t value;
90 typedef bgi::quadratic<8> qRtree;
91 
92 //----- Constants / Enumerations -----------------------------------------------
93 
94 //----- Classes / Structs ------------------------------------------------------
98 {
99 public:
102  fSatisfies(size_t a_nVals)
103  : m_bits(a_nVals, 0)
104  {
105  }
106 
108 
112  template <typename value>
113  bool operator()(value const& a_) const
114  {
115  if (m_bits.test(a_))
116  {
117  return false;
118  }
119  return true;
120  } // operator()
121 };
122 
125 class idx_pt
126 {
127 public:
128  typedef const bPt result_type;
129 
133  explicit idx_pt(const Pt3d* a_, bool a_2d)
134  : m_(a_)
135  , m_2d(a_2d)
136  {
137  m_v.reset();
138  }
139 
143  explicit idx_pt(BSHP<std::vector<Pt3d>> a_, bool a_2d)
144  : m_(0)
145  , m_2d(a_2d)
146  , m_v(a_)
147  {
148  }
149 
153  result_type operator()(size_t i) const
154  {
155  if (m_)
156  return bPt(m_[i].x, m_[i].y, m_2d ? 0 : m_[i].z);
157  else
158  return bPt((*m_v)[i].x, (*m_v)[i].y, m_2d ? 0 : (*m_v)[i].z);
159  }
160 
161  const Pt3d* m_;
162  bool m_2d;
163  BSHP<std::vector<Pt3d>> m_v;
164 };
165 
170 {
171 public:
172  GmPtSearchImpl(bool a_2dSearch);
173  virtual ~GmPtSearchImpl();
174 
175  virtual void PtsToSearch(BSHP<std::vector<Pt3d>> a_pts) override;
176 
177  virtual void VectorThatGrowsToSearch(BSHP<std::vector<Pt3d>> a_) override;
178  virtual bool AddPtToVectorIfUnique(const Pt3d& a_, double a_tol, int& a_ptIdx) override;
179 
180  virtual void NearestPtsToPt(const Pt3d& a_pt,
181  int a_numPtsToFind,
182  bool a_quad_oct_Search,
183  std::vector<int>& a_nearest) const override;
184 
185  virtual void NearestPtsToPtInRtree(int a_ptIdx,
186  const Pt3d& a_pt,
187  int a_numPtsToFind,
188  bool a_quad_oct_Search,
189  std::vector<int>& a_nearest) const override;
190 
191  virtual bool PtInRTree(const Pt3d& a_pt, const double a_tol) override;
192 
193  virtual void PtsWithinDistanceToPtInRtree(int a_ptIdx,
194  const Pt3d& a_pt,
195  double a_dist,
196  std::vector<int>& a_nearest) const override;
197 
198  virtual void NearestPtsToPt(const Pt3d& a_pt,
199  int a_numPtsToFind,
200  bool a_quad_oct_Search,
201  std::vector<int>& a_nearest,
202  fSatisfies* a_fsat) const;
203 
204  virtual void SetActivity(DynBitset& a_activity) override;
205  virtual DynBitset GetActivity() override;
206 
209  virtual const BSHP<VecPt3d> GetPointsPtr() const override { return m_bshpPt3d; }
212  virtual bool Is2D() const override { return m_2dSearch; }
213  virtual std::string ToString() const override;
214 
215  void UpdateMinMax(const Pt3d* a_pts, size_t a_npts);
216  void CreateOctants(const Pt3d& a_pt, std::vector<box>& a_boxes) const;
217 
218  bool m_2dSearch;
222  bgi::rtree<value, qRtree, idx_pt>* m_rTree;
223  BSHP<VecPt3d> m_bshpPt3d;
224 }; // class GmPtSearchImpl
225 
226 //----- Internal functions -----------------------------------------------------
227 
228 //----- Class / Function definitions -------------------------------------------
229 
230 //------------------------------------------------------------------------------
235 //------------------------------------------------------------------------------
236 BSHP<GmPtSearch> GmPtSearch::New(bool a_2dSearch)
237 {
238  BSHP<GmPtSearch> ptr(new GmPtSearchImpl(a_2dSearch));
239  return ptr;
240 } // PtSearch::New
241 //------------------------------------------------------------------------------
243 //------------------------------------------------------------------------------
244 GmPtSearch::GmPtSearch()
245 {
246 } // PtSearch::PtSearch
247 //------------------------------------------------------------------------------
249 //------------------------------------------------------------------------------
250 GmPtSearch::~GmPtSearch()
251 {
252 } // PtSearch::~PtSearch
253 
257 //------------------------------------------------------------------------------
260 //------------------------------------------------------------------------------
262 : m_2dSearch(a_2dSearch)
263 , m_min(XM_FLT_HIGHEST)
264 , m_max(XM_FLT_LOWEST)
265 , m_activity()
266 , m_rTree(nullptr)
267 {
268 } // GmPtSearchImpl::GmPtSearchImpl
269 //------------------------------------------------------------------------------
271 //------------------------------------------------------------------------------
272 GmPtSearchImpl::~GmPtSearchImpl()
273 {
274  if (m_rTree)
275  delete (m_rTree);
276 } // GmPtSearchImpl::~GmPtSearchImpl
277 //------------------------------------------------------------------------------
280 //------------------------------------------------------------------------------
281 void GmPtSearchImpl::PtsToSearch(BSHP<std::vector<Pt3d>> a_pts)
282 {
283  if (m_rTree)
284  delete (m_rTree);
285 
286  m_bshpPt3d = a_pts;
287  qRtree rtree;
288  idx_pt getter(a_pts, m_2dSearch);
289  m_rTree = new bgi::rtree<size_t, qRtree, idx_pt>(
290  boost::counting_iterator<size_t>(0), boost::counting_iterator<size_t>(m_bshpPt3d->size()),
291  rtree, getter);
292 
293  UpdateMinMax(&(*m_bshpPt3d)[0], m_bshpPt3d->size());
294 } // GmPtSearchImpl::PtsToSearch
295 //------------------------------------------------------------------------------
299 //------------------------------------------------------------------------------
300 void GmPtSearchImpl::UpdateMinMax(const Pt3d* a_pts, size_t a_nPts)
301 {
302  for (size_t i = 0; i < a_nPts; ++i)
303  {
304  Pt3d pt = a_pts[i];
305  if (pt.x < m_min.x)
306  m_min.x = pt.x;
307  if (pt.y < m_min.y)
308  m_min.y = pt.y;
309  if (pt.z < m_min.z)
310  m_min.z = pt.z;
311  if (pt.x > m_max.x)
312  m_max.x = pt.x;
313  if (pt.y > m_max.y)
314  m_max.y = pt.y;
315  if (pt.z > m_max.z)
316  m_max.z = pt.z;
317  if (m_2dSearch)
318  pt.z = 0;
319  }
320  m_min -= 1;
321  m_max += 1;
322 } // GmPtSearchImpl::UpdateMinMax
323 //------------------------------------------------------------------------------
326 //------------------------------------------------------------------------------
327 void GmPtSearchImpl::VectorThatGrowsToSearch(BSHP<std::vector<Pt3d>> a_)
328 {
329  if (m_rTree)
330  delete (m_rTree);
331  m_bshpPt3d = a_;
332  qRtree rtree;
333  idx_pt getter(a_, m_2dSearch);
334  m_rTree = new bgi::rtree<size_t, qRtree, idx_pt>(boost::counting_iterator<size_t>(0),
335  boost::counting_iterator<size_t>(a_->size()),
336  rtree, getter);
337 
338  if (!m_bshpPt3d->empty())
339  UpdateMinMax(&(*m_bshpPt3d)[0], m_bshpPt3d->size());
340 } // GmPtSearchImpl::VectorThatGrowsToSearch
341 //------------------------------------------------------------------------------
349 //------------------------------------------------------------------------------
350 bool GmPtSearchImpl::AddPtToVectorIfUnique(const Pt3d& a_pt, double a_tol, int& a_ptIdx)
351 {
352  std::vector<int> n;
353  PtsWithinDistanceToPtInRtree(-1, a_pt, a_tol, n);
354  if (!n.empty())
355  {
356  a_ptIdx = *std::min_element(n.begin(), n.end());
357  return false;
358  }
359  value val(m_bshpPt3d->size());
360  a_ptIdx = (int)m_bshpPt3d->size();
361  m_bshpPt3d->push_back(a_pt);
362  m_rTree->insert(val);
363  // update min/max
364  if (a_pt.x < m_min.x)
365  m_min.x = a_pt.x;
366  if (a_pt.y < m_min.y)
367  m_min.y = a_pt.y;
368  if (a_pt.z < m_min.z)
369  m_min.z = a_pt.z;
370  if (a_pt.x > m_max.x)
371  m_max.x = a_pt.x;
372  if (a_pt.y > m_max.y)
373  m_max.y = a_pt.y;
374  if (a_pt.z > m_max.z)
375  m_max.z = a_pt.z;
376  // if (!m_bshpPt3d->empty()) UpdateMinMax(&(*m_bshpPt3d)[0], m_bshpPt3d->size());
377  return true;
378 } // GmPtSearchImpl::AddPtToVectorIfUnique
379 //------------------------------------------------------------------------------
387 //------------------------------------------------------------------------------
389  int a_numPtsToFind,
390  bool a_quad_oct_Search,
391  std::vector<int>& a_nearest) const
392 {
393  fSatisfies *fs(0), fs1(1);
394  if (!m_activity.empty())
395  {
396  fs1.m_bits = m_activity;
397  fs = &fs1;
398  }
399  NearestPtsToPt(a_pt, a_numPtsToFind, a_quad_oct_Search, a_nearest, fs);
400 } // NearestPtsToPt
401 //------------------------------------------------------------------------------
415 //------------------------------------------------------------------------------
417  int a_numPtsToFind,
418  bool a_quad_oct_Search,
419  std::vector<int>& a_nearest,
420  fSatisfies* a_fsat) const
421 {
422  std::vector<value> nearest;
423  a_nearest.resize(0);
424  if (!m_bshpPt3d || m_bshpPt3d->empty())
425  {
426  XM_LOG(xmlog::error, "Unable to find nearest points; no points exist.");
427  }
428 
429  bPt tmpPt(a_pt);
430  if (!a_quad_oct_Search)
431  {
432  if (a_fsat)
433  m_rTree->query(bgi::satisfies(*a_fsat) && bgi::nearest(tmpPt, a_numPtsToFind),
434  std::back_inserter(nearest));
435  else
436  m_rTree->query(bgi::nearest(tmpPt, a_numPtsToFind), std::back_inserter(nearest));
437  a_nearest.reserve(nearest.size());
438  for (size_t i = 0; i < nearest.size(); ++i)
439  {
440  a_nearest.push_back((int)nearest[i]);
441  }
442  }
443  else
444  {
445  fSatisfies *fPtr(a_fsat), tmpFsat(m_rTree->size());
446  if (!fPtr)
447  fPtr = &tmpFsat;
448  std::vector<box> boxes;
449  CreateOctants(a_pt, boxes);
450  for (size_t i = 0; i < boxes.size(); ++i)
451  {
452  std::vector<value> foundPts;
453  m_rTree->query(
454  bgi::covered_by(boxes[i]) && bgi::satisfies(*fPtr) && bgi::nearest(tmpPt, a_numPtsToFind),
455  std::back_inserter(foundPts));
456  a_nearest.reserve(a_nearest.size() + foundPts.size());
457  for (size_t j = 0; j < foundPts.size(); ++j)
458  {
459  a_nearest.push_back((int)foundPts[j]);
460  fPtr->m_bits.set(foundPts[j]);
461  }
462  }
463  }
464  std::sort(a_nearest.begin(), a_nearest.end());
465 } // GmPtSearchImpl::NearestPtsToPt
466 //------------------------------------------------------------------------------
478 //------------------------------------------------------------------------------
480  const Pt3d& a_pt,
481  int a_numPtsToFind,
482  bool a_quad_oct_Search,
483  std::vector<int>& a_nearest) const
484 {
485  fSatisfies fsat(m_rTree->size());
486  if (!m_activity.empty())
487  fsat.m_bits = m_activity;
488  fsat.m_bits.set(a_ptIdx);
489  NearestPtsToPt(a_pt, a_numPtsToFind, a_quad_oct_Search, a_nearest, &fsat);
490 } // GmPtSearchImpl::NearestPtsToPtInRtree
491 //------------------------------------------------------------------------------
497 //------------------------------------------------------------------------------
498 bool GmPtSearchImpl::PtInRTree(const Pt3d& a_pt, const double a_tol)
499 {
500  std::vector<int> n;
501  PtsWithinDistanceToPtInRtree(-1, a_pt, a_tol, n);
502  return !n.empty();
503 } // GmPtSearchImpl::PtInRTree
504 //------------------------------------------------------------------------------
516 //------------------------------------------------------------------------------
518  const Pt3d& a_pt,
519  double a_distance,
520  std::vector<int>& a_nearest) const
521 {
522  fSatisfies fsat(m_rTree->size());
523  if (!m_activity.empty())
524  fsat.m_bits = m_activity;
525  if (a_ptIdx > -1)
526  fsat.m_bits.set(a_ptIdx);
527 
528  Pt3d bMin(a_pt - a_distance), bMax(a_pt + a_distance);
529  if (m_2dSearch)
530  {
531  bMin.z = -1;
532  bMax.z = 1;
533  }
534  bPt boxMin, boxMax;
535  box aBox(bMin, bMax);
536 
537  std::vector<value> nearest;
538  a_nearest.resize(0);
539 
540  m_rTree->query(bgi::satisfies(fsat) && bgi::covered_by(aBox), std::back_inserter(nearest));
541  a_nearest.reserve(nearest.size());
542  for (size_t i = 0; i < nearest.size(); ++i)
543  {
544  a_nearest.push_back((int)nearest[i]);
545  }
546 } // GmPtSearchImpl::PtsWithinDistanceToPtInRtree
547 //------------------------------------------------------------------------------
552 //------------------------------------------------------------------------------
554 {
555  if (a_activity.size() == m_rTree->size())
556  {
557  m_activity = a_activity;
558  m_activity.flip();
559  }
560 } // SetActivity
561 //------------------------------------------------------------------------------
564 //------------------------------------------------------------------------------
566 {
567  DynBitset act(m_activity);
568  act.flip();
569  return act;
570 } // GmPtSearchImpl::GetActivity
571 //------------------------------------------------------------------------------
574 //------------------------------------------------------------------------------
575 std::string GmPtSearchImpl::ToString() const
576 {
577  std::stringstream ss;
578  ss << m_2dSearch << "=2dSearch " << m_min << "=min " << m_max << "=max "
579  << "\n";
580 
581  VecToStream(ss, m_activity, "activity");
582 
583  if (m_bshpPt3d)
584  VecToStream(ss, *m_bshpPt3d, "bshpPt3d");
585 
586  return ss.str();
587 } // GmPtSearchImpl::ToString
588 //------------------------------------------------------------------------------
592 //------------------------------------------------------------------------------
593 void GmPtSearchImpl::CreateOctants(const Pt3d& a_pt, std::vector<box>& a_boxes) const
594 {
595  a_boxes.resize(0);
596  Pt3d bmin(m_min), bmax(m_max), pt(a_pt);
597  box bound(bmin, bmax);
598  if (!m_2dSearch)
599  bound.min_corner().z = pt.z;
600  // here are the quadrants that we create
601  // ---------
602  // | 1 | 2 |
603  // ---------
604  // | 3 | 4 |
605  // ---------
606  box aBox(bound);
607  // 1
608  aBox.max_corner().x = pt.x;
609  aBox.min_corner().y = pt.y;
610  a_boxes.push_back(aBox);
611  // 2
612  aBox = bound;
613  aBox.min_corner() = pt;
614  if (m_2dSearch)
615  aBox.min_corner().z = -1;
616  a_boxes.push_back(aBox);
617  // 3
618  aBox = bound;
619  aBox.max_corner() = pt;
620  aBox.max_corner().z = m_max.z;
621  if (m_2dSearch)
622  aBox.max_corner().z = 1;
623  a_boxes.push_back(aBox);
624  // 4
625  aBox = bound;
626  aBox.min_corner().x = pt.x;
627  aBox.max_corner().y = pt.y;
628  a_boxes.push_back(aBox);
629  if (m_2dSearch)
630  return;
631  // do the lower octants if this is 3D
632  for (int i = 0; i < 4; ++i)
633  {
634  a_boxes.push_back(a_boxes[i]);
635  a_boxes.back().max_corner().z = pt.z;
636  a_boxes.back().min_corner().z = m_min.z;
637  }
638 } // GmPtSearchImpl::CreateOctants
639 
640 } // namespace xms
641 
642 #ifdef CXX_TEST
643 
646 #include <boost/assign.hpp>
648 
649 // namespace xms {
650 using namespace xms;
651 
652 //------------------------------------------------------------------------------
657 //------------------------------------------------------------------------------
658 void iGetPoints2(std::vector<Pt3d>& a_pts, std::vector<float>& a_scalar, std::vector<Pt3d>& a_iPts)
659 {
660  a_pts = {{-75.0, -16.0, 0.0}, {-60.0, 32.0, 0.0}, {-34.0, 50.0, 0.0}, {-34.0, 27.0, 0.0},
661  {-8.0, 40.0, 0.0}, {16.0, 38.0, 0.0}, {-25.0, 14.0, 43.64}, {10.0, 18.0, 44.16},
662  {27.0, 26.0, 0.0}, {63.0, 35.0, 0.0}, {-32.0, 0.0, 59.04}, {-7.0, 7.0, 90.2},
663  {26.0, 6.0, 67.2}, {75.0, 7.0, 0.0}, {-37.0, -15.0, 9.24}, {-7.0, -13.0, 71.0},
664  {2.0, -3.0, 98.4}, {31.0, -15.0, 25.56}, {60.0, -13.0, 0.0}, {-50.0, -30.0, 0.0},
665  {-30.0, -28.0, 0.0}, {43.0, -22.0, 0.0}, {-32.0, -50.0, 0.0}, {27.0, -37.0, 0.0},
666  {60.0, -33.0, 0.0}};
667 
668  a_scalar.resize(a_pts.size());
669  for (size_t i = 0; i < a_pts.size(); ++i)
670  a_scalar[i] = (float)a_pts[i].z;
671  a_iPts = {{-90, 60, 0}, {90, 60, 0}, {-90, -60, 0}, {90, -60, 0}, {60, -33, 0}};
672 } // iGetPoints
673 //------------------------------------------------------------------------------
675 #define _TS_ASSERT_POINTS_FOUND(a_file, a_line, a_required, a_optional, a_found) \
676  iAssertPointsFound(a_file, a_line, a_required, a_optional, a_found)
677 #define TS_ASSERT_POINTS_FOUND(a_required, a_optional, a_found) \
679  _TS_ASSERT_POINTS_FOUND(__FILE__, __LINE__, a_required, a_optional, a_found)
680 //------------------------------------------------------------------------------
687 //------------------------------------------------------------------------------
688 void iAssertPointsFound(const char* a_file,
689  int a_line,
690  const std::vector<int>& a_required,
691  const std::vector<int>& a_optional,
692  const std::vector<int>& a_found)
693 {
694  // look for required points
695  std::vector<bool> handled(a_found.size(), false);
696  for (size_t i = 0; i < a_required.size(); ++i)
697  {
698  int idx = a_required[i];
699  auto it = std::find(a_found.begin(), a_found.end(), idx);
700  bool found = it != a_found.end();
701  if (found)
702  {
703  handled[it - a_found.begin()] = found;
704  }
705  else
706  {
707  std::stringstream ss;
708  ss << "Result failed to include " << idx;
709  _TS_FAIL(a_file, a_line, ss.str().c_str());
710  }
711  }
712 
713  // ensure all other points are optional
714  for (size_t i = 0; i < handled.size(); ++i)
715  {
716  if (!handled[i])
717  {
718  int idx = a_found[i];
719  bool found = std::find(a_optional.begin(), a_optional.end(), idx) != a_optional.end();
720  if (!found)
721  {
722  std::stringstream ss;
723  ss << "Failed to find " << idx << " in optional points";
724  _TS_FAIL(a_file, a_line, ss.str().c_str());
725  }
726  }
727  }
728 } // iAssertPointsFound
733 //------------------------------------------------------------------------------
735 //------------------------------------------------------------------------------
737 {
738  BSHP<GmPtSearch> pts = GmPtSearch::New(true);
739  TS_ASSERT(pts);
740 }
741 //------------------------------------------------------------------------------
745 //------------------------------------------------------------------------------
746 static void iSetupPts(std::vector<Pt3d>& pts, bool a_2d)
747 {
748  pts.resize(0);
749  pts.reserve(300);
750  for (int i = 1; i < 11; ++i)
751  {
752  for (int j = 1; j < 11; ++j)
753  {
754  pts.push_back(Pt3d(i, j, 0));
755  if (i % 2)
756  pts.push_back(Pt3d(-i, j, 0));
757  pts.push_back(Pt3d(i, -j, 0));
758  if (i % 2)
759  pts.push_back(Pt3d(-i, -j, 0));
760  }
761  }
762  if (a_2d)
763  return;
764  std::vector<Pt3d> pts1(pts);
765  for (size_t i = 0; i < pts1.size(); ++i)
766  pts1[i].z = .1 * (double)(i + 1);
767  pts.insert(pts.end(), pts1.begin(), pts1.end());
768  for (size_t i = 0; i < pts1.size(); ++i)
769  pts1[i].z = -.1 * (double)(i + 1);
770  pts.insert(pts.end(), pts1.begin(), pts1.end());
771 } // iSetupPts
772 //------------------------------------------------------------------------------
774 //------------------------------------------------------------------------------
776 {
777  BSHP<std::vector<Pt3d>> pts(new std::vector<Pt3d>());
778  iSetupPts(*pts, true);
779 
780  BSHP<GmPtSearch> iPts = GmPtSearch::New(true);
781  iPts->PtsToSearch(pts);
782  TS_ASSERT_EQUALS(true, iPts->Is2D());
783 
784  Pt3d aPt(.5, 0, 0);
785  std::vector<int> foundPts, requiredPts, optionalPts;
786  iPts->NearestPtsToPt(aPt, 8, false, foundPts);
787  requiredPts = {0, 1, 2, 3, 4, 6, 40, 41};
788  std::sort(foundPts.begin(), foundPts.end());
789  TS_ASSERT_EQUALS(8, foundPts.size());
790  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
791 
792  iPts->NearestPtsToPt(aPt, 4, true, foundPts);
793  requiredPts = {0, 1, 2, 3, 4, 5, 6, 7, 9, 11, 40, 41, 42, 43, 61, 63};
794  std::sort(foundPts.begin(), foundPts.end());
795  TS_ASSERT_EQUALS(16, foundPts.size());
796  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
797 
798  aPt.x = -3.6;
799  aPt.y = 2.5;
800  iPts->NearestPtsToPt(aPt, 8, false, foundPts);
801  requiredPts = {61, 65, 69, 73, 121, 125, 129, 133};
802  std::sort(foundPts.begin(), foundPts.end());
803  TS_ASSERT_EQUALS(8, foundPts.size());
804  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
805 
806  iPts->NearestPtsToPt(aPt, 4, true, foundPts);
807  requiredPts = {1, 5, 9, 61, 65, 69, 73, 77, 121, 125, 129, 133, 137, 181, 185, 189};
808  std::sort(foundPts.begin(), foundPts.end());
809  TS_ASSERT_EQUALS(16, foundPts.size());
810  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
811 
812  aPt.x = -1;
813  aPt.y = 0;
814  iPts->NearestPtsToPt(aPt, 4, true, foundPts);
815  requiredPts = {0, 1, 2, 3, 4, 5, 6, 7, 9, 11, 40, 41, 61, 63, 65, 67};
816  std::sort(foundPts.begin(), foundPts.end());
817  TS_ASSERT_EQUALS(16, foundPts.size());
818  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
819 
820  iPts->NearestPtsToPt(aPt, 4, true, foundPts);
821  requiredPts = {0, 1, 2, 3, 4, 5, 6, 7, 9, 11, 40, 41, 61, 63, 65, 67};
822  std::sort(foundPts.begin(), foundPts.end());
823  TS_ASSERT_EQUALS(16, foundPts.size());
824  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
825 } // PtSearchUnitTests::testTest2d
826 //------------------------------------------------------------------------------
828 //------------------------------------------------------------------------------
830 {
831  BSHP<std::vector<Pt3d>> pts(new std::vector<Pt3d>());
832  std::vector<Pt3d> ipts;
833  std::vector<float> scalar, s, vBase;
834  iGetPoints2(*pts, scalar, ipts);
835 
836  BSHP<GmPtSearch> iPts = GmPtSearch::New(true);
837  iPts->PtsToSearch(pts);
838 
839  std::vector<int> foundPts, requiredPts, optionalPts;
840  iPts->NearestPtsToPt(ipts[0], 16, false, foundPts);
841  requiredPts = {0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 14, 15, 16, 19, 20};
842  std::sort(foundPts.begin(), foundPts.end());
843  TS_ASSERT_EQUALS(16, foundPts.size());
844  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
845 } // PtSearchUnitTests::testTest2dTutData
846 //------------------------------------------------------------------------------
848 //------------------------------------------------------------------------------
850 {
851  BSHP<std::vector<Pt3d>> pts(new std::vector<Pt3d>());
852  iSetupPts(*pts, false);
853 
854  BSHP<GmPtSearch> iPts = GmPtSearch::New(false);
855  iPts->PtsToSearch(pts);
856 
857  Pt3d aPt(0, 0, 0);
858  std::vector<int> foundPts, requiredPts, optionalPts;
859  iPts->NearestPtsToPt(aPt, 16, false, foundPts);
860  requiredPts = {0, 1, 2, 3, 300, 301, 302, 303, 600, 601, 602, 603};
861  optionalPts = {4, 5, 6, 7, 40, 41};
862  TS_ASSERT_EQUALS(16, foundPts.size());
863  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
864 }
865 //------------------------------------------------------------------------------
868 //------------------------------------------------------------------------------
869 static void iSetupPtsOctant(std::vector<Pt3d>& pts)
870 {
871  // points along the axis
872  pts.push_back(Pt3d(1, 0, 0));
873  pts.push_back(Pt3d(2, 0, 0));
874  pts.push_back(Pt3d(-1, 0, 0));
875  pts.push_back(Pt3d(-2, 0, 0));
876  pts.push_back(Pt3d(0, 1, 0));
877  pts.push_back(Pt3d(0, 2, 0));
878  pts.push_back(Pt3d(0, -1, 0));
879  pts.push_back(Pt3d(0, -2, 0));
880  pts.push_back(Pt3d(0, 0, 1));
881  pts.push_back(Pt3d(0, 0, 2));
882  pts.push_back(Pt3d(0, 0, -1));
883  pts.push_back(Pt3d(0, 0, -2));
884 
885  // points in the octants
886  // positive z octants
887  pts.push_back(Pt3d(1, 1, 2));
888  pts.push_back(Pt3d(2, 1, 2));
889  pts.push_back(Pt3d(1, 2, 2));
890  pts.push_back(Pt3d(1, 1, 3));
891 
892  pts.push_back(Pt3d(-1, 1, 2));
893  pts.push_back(Pt3d(-2, 1, 2));
894  pts.push_back(Pt3d(-1, 2, 2));
895  pts.push_back(Pt3d(-1, 1, 3));
896 
897  pts.push_back(Pt3d(1, -1, 2));
898  pts.push_back(Pt3d(2, -1, 2));
899  pts.push_back(Pt3d(1, -2, 2));
900  pts.push_back(Pt3d(1, -1, 3));
901 
902  pts.push_back(Pt3d(-1, -1, 2));
903  pts.push_back(Pt3d(-2, -1, 2));
904  pts.push_back(Pt3d(-1, -2, 2));
905  pts.push_back(Pt3d(-1, -1, 3));
906 
907  // negative z octants
908  pts.push_back(Pt3d(1, 1, -2));
909  pts.push_back(Pt3d(2, 1, -2));
910  pts.push_back(Pt3d(1, 2, -2));
911  pts.push_back(Pt3d(1, 1, -3));
912 
913  pts.push_back(Pt3d(-1, 1, -2));
914  pts.push_back(Pt3d(-2, 1, -2));
915  pts.push_back(Pt3d(-1, 2, -2));
916  pts.push_back(Pt3d(-1, 1, -3));
917 
918  pts.push_back(Pt3d(1, -1, -2));
919  pts.push_back(Pt3d(2, -1, -2));
920  pts.push_back(Pt3d(1, -2, -2));
921  pts.push_back(Pt3d(1, -1, -3));
922 
923  pts.push_back(Pt3d(-1, -1, -2));
924  pts.push_back(Pt3d(-2, -1, -2));
925  pts.push_back(Pt3d(-1, -2, -2));
926  pts.push_back(Pt3d(-1, -1, -3));
927 
928 } // iSetupPtsOctant
929 //------------------------------------------------------------------------------
931 //------------------------------------------------------------------------------
933 {
934  BSHP<std::vector<Pt3d>> pts(new std::vector<Pt3d>());
935  iSetupPtsOctant(*pts);
936 
937  BSHP<GmPtSearch> iPts = GmPtSearch::New(false);
938  iPts->PtsToSearch(pts);
939 
940  Pt3d aPt;
941  std::vector<int> foundPts;
942  iPts->NearestPtsToPt(aPt, 4, true, foundPts);
943  std::vector<int> requiredPts = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 20, 21, 22,
944  23, 24, 28, 29, 30, 31, 32, 36, 37, 38, 39, 40, 41, 42, 43};
945  std::vector<int> optionalPts = {33, 34};
946  std::sort(foundPts.begin(), foundPts.end());
947  TS_ASSERT_EQUALS(32, foundPts.size());
948  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
949 
950  iPts->NearestPtsToPt(aPt, 4, true, foundPts);
951  std::sort(foundPts.begin(), foundPts.end());
952  TS_ASSERT_EQUALS(32, foundPts.size());
953  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
954 } // PtSearchUnitTests::testTest3d
955 //------------------------------------------------------------------------------
957 //------------------------------------------------------------------------------
959 {
960  BSHP<std::vector<Pt3d>> pts(new std::vector<Pt3d>());
961  iSetupPts(*pts, true);
962 
963  BSHP<GmPtSearch> iPts = GmPtSearch::New(true);
964  iPts->PtsToSearch(pts);
965 
966  DynBitset activity_wrongSize;
967  // this should do nothing because it is sized wrong
968  iPts->SetActivity(activity_wrongSize);
969 
970  Pt3d aPt(.5, 0, 0);
971  std::vector<int> foundPts, requiredPts, optionalPts;
972  iPts->NearestPtsToPt(aPt, 8, false, foundPts);
973  requiredPts = {0, 1, 2, 3, 4, 6, 40, 41};
974  std::sort(foundPts.begin(), foundPts.end());
975  TS_ASSERT_EQUALS(8, foundPts.size());
976  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
977 
978  DynBitset act(300);
979  act.flip();
980  act.set(2, 0);
981  act.set(40, 0);
982  iPts->SetActivity(act);
983  iPts->NearestPtsToPt(aPt, 8, false, foundPts);
984  requiredPts = {0, 1, 3, 4, 6, 41, 42};
985  optionalPts = {7, 43};
986  std::sort(foundPts.begin(), foundPts.end());
987  TS_ASSERT_EQUALS(8, foundPts.size());
988  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
989 } // PtSearchUnitTests::testActivity2d
990 //------------------------------------------------------------------------------
992 //------------------------------------------------------------------------------
994 {
995  BSHP<std::vector<Pt3d>> pts(new std::vector<Pt3d>());
996  iSetupPtsOctant(*pts);
997 
998  BSHP<GmPtSearch> iPts = GmPtSearch::New(false);
999  iPts->PtsToSearch(pts);
1000 
1001  DynBitset activity_wrongSize;
1002  // this should do nothing because it is sized wrong
1003  iPts->SetActivity(activity_wrongSize);
1004 
1005  Pt3d aPt;
1006  std::vector<int> foundPts;
1007  iPts->NearestPtsToPt(aPt, 4, true, foundPts);
1008  std::vector<int> requiredPts = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 20, 21, 22,
1009  23, 24, 28, 29, 30, 31, 32, 36, 37, 38, 39, 40, 41, 42, 43};
1010  std::vector<int> optionalPts = {33, 34};
1011  std::sort(foundPts.begin(), foundPts.end());
1012  TS_ASSERT_EQUALS(32, foundPts.size());
1013  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
1014 
1015  DynBitset act(44);
1016  act.flip();
1017  act.set(2, 0);
1018  act.set(40, 0);
1019  iPts->SetActivity(act);
1020  iPts->NearestPtsToPt(aPt, 4, true, foundPts);
1021  requiredPts = {0, 1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 20, 21, 22,
1022  23, 24, 28, 29, 30, 31, 32, 36, 37, 38, 39, 41, 42, 43};
1023  optionalPts = {13, 14, 25, 33, 34};
1024  std::sort(foundPts.begin(), foundPts.end());
1025  TS_ASSERT_EQUALS(31, foundPts.size());
1026  TS_ASSERT_POINTS_FOUND(requiredPts, optionalPts, foundPts);
1027 } // PtSearchUnitTests::testActivity3d
1028 //------------------------------------------------------------------------------
1030 //------------------------------------------------------------------------------
1032 {
1033  BSHP<std::vector<Pt3d>> pts(new std::vector<Pt3d>());
1034  *pts = {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {1, 1, 0}, {1, -1, 0}};
1035 
1036  GmPtSearchImpl p(true);
1037  p.PtsToSearch(pts);
1038 
1039  std::vector<int> n;
1040  p.PtsWithinDistanceToPtInRtree(1, (*pts)[1], .99, n);
1041  std::vector<int> base;
1042  TS_ASSERT_EQUALS_VEC(base, n);
1043  base = {0, 2, 3, 4};
1044  p.PtsWithinDistanceToPtInRtree(1, (*pts)[1], 1.0, n);
1045  TS_ASSERT_EQUALS_VEC(base, n);
1046  base = {1, 3, 4};
1047  p.PtsWithinDistanceToPtInRtree(2, (*pts)[2], 1.0, n);
1048  TS_ASSERT_EQUALS_VEC(base, n);
1049 }
1050 //------------------------------------------------------------------------------
1052 //------------------------------------------------------------------------------
1054 {
1055  BSHP<std::vector<Pt3d>> pts(new std::vector<Pt3d>());
1056  *pts = {{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {1, 1, 0}, {1, -1, 0}};
1057 
1058  GmPtSearchImpl p(true);
1059  p.VectorThatGrowsToSearch(pts);
1060 
1061  // add a point
1062  int ptIdx;
1063  Pt3d pt(0, 1, 0);
1064  TS_ASSERT_EQUALS(false, p.AddPtToVectorIfUnique(pt, 1, ptIdx));
1065  TS_ASSERT_EQUALS(0, ptIdx);
1066  TS_ASSERT_EQUALS(true, p.AddPtToVectorIfUnique(pt, 1e-9, ptIdx));
1067  TS_ASSERT_EQUALS(5, ptIdx);
1068 } // PtSearchUnitTests::testVectorThatGrows
1069 
1070 //} // namespace xms
1071 #endif // CXX_TEST
virtual bool PtInRTree(const Pt3d &a_pt, const double a_tol) override
Checks if the point is in the Rtree within tolerance.
Definition: GmPtSearch.cpp:498
#define XM_LOG(A, B)
Spatial index for searching points.
Definition: GmPtSearch.h:27
Pt3< double > Pt3d
DynBitset m_activity
point activity
Definition: GmPtSearch.cpp:221
void iGetPoints2(std::vector< Pt3d > &a_pts, std::vector< float > &a_scalar, std::vector< Pt3d > &a_iPts)
Get GMS tutorial data.
Definition: GmPtSearch.cpp:658
bool m_2dSearch
flag specifying 2d searching only
Definition: GmPtSearch.cpp:218
BOOST_GEOMETRY_REGISTER_POINT_3D(bPt, double, cs::cartesian, x, y, z)
result_type operator()(size_t i) const
operator to return the point location based on an index
Definition: GmPtSearch.cpp:153
bool operator()(value const &a_) const
Determines if a point can be returned. If the bit is not set yet.
Definition: GmPtSearch.cpp:113
static void iSetupPts(std::vector< Pt3d > &pts, bool a_2d)
helper function for testing
Definition: GmPtSearch.cpp:746
virtual bool Is2D() const override
Returns true if class only searches in 2D.
Definition: GmPtSearch.cpp:212
void testTest3d()
testing in 3d
Definition: GmPtSearch.cpp:849
const bPt result_type
created to follow boost examples
Definition: GmPtSearch.cpp:128
virtual void PtsToSearch(BSHP< std::vector< Pt3d >> a_pts) override
Adds the point locations to the class.
Definition: GmPtSearch.cpp:281
GmPtSearchImpl(bool a_2dSearch)
Constructor.
Definition: GmPtSearch.cpp:261
virtual bool AddPtToVectorIfUnique(const Pt3d &a_, double a_tol, int &a_ptIdx) override
Adds the point to the R-tree if the point location is unique based on the passed in tolerance...
Definition: GmPtSearch.cpp:350
fSatisfies(size_t a_nVals)
Constructor.
Definition: GmPtSearch.cpp:102
void VecToStream(std::stringstream &a_ss, const T &a_v, std::string a_label)
Pt3d m_min
minimum extents of points
Definition: GmPtSearch.cpp:219
static void iSetupPtsOctant(std::vector< Pt3d > &pts)
helper function for testing
Definition: GmPtSearch.cpp:869
void testActivity2d()
testing searching with activity
Definition: GmPtSearch.cpp:958
#define TS_ASSERT_POINTS_FOUND(a_required, a_optional, a_found)
macro for testing
Definition: GmPtSearch.cpp:678
BSHP< std::vector< Pt3d > > m_v
vector of point locations
Definition: GmPtSearch.cpp:163
DynBitset m_bits
bitset that is the size of the points
Definition: GmPtSearch.cpp:107
boost::dynamic_bitset< size_t > DynBitset
void testCreateClass()
tests creating the class
Definition: GmPtSearch.cpp:736
const Pt3d * m_
array of point locations
Definition: GmPtSearch.cpp:161
virtual void NearestPtsToPtInRtree(int a_ptIdx, const Pt3d &a_pt, int a_numPtsToFind, bool a_quad_oct_Search, std::vector< int > &a_nearest) const override
Finds the nearest points to the input a_pt with an option to search quadrants/octants. Similar to NearestPtsToPt but this method will always pass the fSatisfies class to the RTree. This method is used to find the nearest points to a point that is included in the RTree.
Definition: GmPtSearch.cpp:479
virtual DynBitset GetActivity() override
Returns the point activity.
Definition: GmPtSearch.cpp:565
void CreateOctants(const Pt3d &a_pt, std::vector< box > &a_boxes) const
Creates octants (or quadrants for 2d) to be used in the rtree query.
Definition: GmPtSearch.cpp:593
void testTest2dTutData()
testing tutorial data from GMS in 2d
Definition: GmPtSearch.cpp:829
virtual void PtsWithinDistanceToPtInRtree(int a_ptIdx, const Pt3d &a_pt, double a_dist, std::vector< int > &a_nearest) const override
Finds the nearest points to the input a_pt with an option to search quadrants/octants. Similar to NearestPtsToPt but this method will always pass the fSatisfies class to the RTree. This method is used to find the nearest points to a point that is included in the RTree.
Definition: GmPtSearch.cpp:517
virtual void NearestPtsToPt(const Pt3d &a_pt, int a_numPtsToFind, bool a_quad_oct_Search, std::vector< int > &a_nearest) const override
Finds the nearest points to the input a_pt with an option to search quadrants/octants.
Definition: GmPtSearch.cpp:388
static const float XM_FLT_HIGHEST
class for indexing the points
Definition: GmPtSearch.cpp:125
#define BSHP
virtual void VectorThatGrowsToSearch(BSHP< std::vector< Pt3d >> a_) override
Adds the point locations to the class.
Definition: GmPtSearch.cpp:327
void UpdateMinMax(const Pt3d *a_pts, size_t a_npts)
Updates the m_min, m_max variables.
Definition: GmPtSearch.cpp:300
a class used with the boost::geometry::index::satisfies function
Definition: GmPtSearch.cpp:97
void testTest3dOct()
testing octant searching in 3d
Definition: GmPtSearch.cpp:932
bool m_2d
flag specifying 2d searching only
Definition: GmPtSearch.cpp:162
idx_pt(const Pt3d *a_, bool a_2d)
constructor
Definition: GmPtSearch.cpp:133
static BSHP< GmPtSearch > New(bool a_2dSearch)
Creates an PtSearch class.
Definition: GmPtSearch.cpp:236
idx_pt(BSHP< std::vector< Pt3d >> a_, bool a_2d)
constructor
Definition: GmPtSearch.cpp:143
Pt3d m_max
maximum extents of points
Definition: GmPtSearch.cpp:220
void testVectorThatGrows()
testing VectorThatGrows functionality
bgi::rtree< value, qRtree, idx_pt > * m_rTree
spatial index for searching the points
Definition: GmPtSearch.cpp:222
virtual void SetActivity(DynBitset &a_activity) override
Sets activity on the points in the Rtree so that points can be ignored when interpolating.
Definition: GmPtSearch.cpp:553
virtual const BSHP< VecPt3d > GetPointsPtr() const override
Returns shared point to the point locations in the RTree.
Definition: GmPtSearch.cpp:209
Implementation of GmPtSearch. Generic class for searching location data. Uses a boost R-tree to query...
Definition: GmPtSearch.cpp:169
BSHP< VecPt3d > m_bshpPt3d
Vector of point locations.
Definition: GmPtSearch.cpp:223
void testActivity3d()
testing searching with activity
Definition: GmPtSearch.cpp:993
#define TS_ASSERT_EQUALS_VEC(a, b)
void testPtsWithinDist()
testing Point within a distance to a specified point
void testTest2d()
testing in 2d
Definition: GmPtSearch.cpp:775
void iAssertPointsFound(const char *a_file, int a_line, const std::vector< int > &a_required, const std::vector< int > &a_optional, const std::vector< int > &a_found)
helper function for testing
Definition: GmPtSearch.cpp:688
virtual std::string ToString() const override
Write the internals to a string.
Definition: GmPtSearch.cpp:575
static const float XM_FLT_LOWEST