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