GEOS  3.10.1
TemplateSTRtreeDistance.h
1 /**********************************************************************
2  *
3  * GEOS - Geometry Engine Open Source
4  * http://geos.osgeo.org
5  *
6  * Copyright (C) 2020-2021 Daniel Baston
7  *
8  * This is free software; you can redistribute and/or modify it under
9  * the terms of the GNU Lesser General Public Licence as published
10  * by the Free Software Foundation.
11  * See the COPYING file for more information.
12  *
13  **********************************************************************/
14 
15 #pragma once
16 
17 #include <geos/index/strtree/TemplateSTRNode.h>
18 #include <geos/index/strtree/TemplateSTRNodePair.h>
19 
20 #include <queue>
21 #include <vector>
22 
23 namespace geos {
24 namespace index {
25 namespace strtree {
26 
27 template<typename ItemType, typename BoundsType, typename ItemDistance>
28 class TemplateSTRtreeDistance {
29  using Node = TemplateSTRNode<ItemType, BoundsType>;
30  using NodePair = TemplateSTRNodePair<ItemType, BoundsType, ItemDistance>;
31  using ItemPair = std::pair<ItemType, ItemType>;
32 
33  struct PairQueueCompare {
34  bool operator()(const NodePair& a, const NodePair& b) {
35  return a.getDistance() > b.getDistance();
36  }
37  };
38 
39  using PairQueue = std::priority_queue<NodePair, std::vector<NodePair>, PairQueueCompare>;
40 
41 public:
42  explicit TemplateSTRtreeDistance(ItemDistance& id) : m_id(id) {}
43 
44  ItemPair nearestNeighbour(const Node& root1, const Node& root2) {
45  NodePair initPair(root1, root2, m_id);
46  return nearestNeighbour(initPair);
47  }
48 
49  ItemPair nearestNeighbour(NodePair& initPair) {
50  return nearestNeighbour(initPair, DoubleInfinity);
51  }
52 
53 private:
54 
55  ItemPair nearestNeighbour(NodePair& initPair, double maxDistance) {
56  double distanceLowerBound = maxDistance;
57  std::unique_ptr<NodePair> minPair;
58 
59  PairQueue priQ;
60  priQ.push(initPair);
61 
62  while (!priQ.empty() && distanceLowerBound > 0) {
63  NodePair pair = priQ.top();
64  priQ.pop();
65  double currentDistance = pair.getDistance();
66 
67  /*
68  * If the distance for the first node in the queue
69  * is >= the current minimum distance, all other nodes
70  * in the queue must also have a greater distance.
71  * So the current minDistance must be the true minimum,
72  * and we are done.
73  */
74  if (minPair && currentDistance >= distanceLowerBound) {
75  break;
76  }
77 
78  /*
79  * If the pair members are leaves
80  * then their distance is the exact lower bound.
81  * Update the distanceLowerBound to reflect this
82  * (which must be smaller, due to the test
83  * immediately prior to this).
84  */
85  if (pair.isLeaves()) {
86  distanceLowerBound = currentDistance;
87  if (minPair) {
88  *minPair = pair;
89  } else {
90  minPair = detail::make_unique<NodePair>(pair);
91  }
92  } else {
93  /*
94  * Otherwise, expand one side of the pair,
95  * (the choice of which side to expand is heuristically determined)
96  * and insert the new expanded pairs into the queue
97  */
98  expandToQueue(pair, priQ, distanceLowerBound);
99  }
100  }
101 
102  if (!minPair) {
103  throw util::GEOSException("Error computing nearest neighbor");
104  }
105 
106  return minPair->getItems();
107  }
108 
109  void expandToQueue(const NodePair& pair, PairQueue& priQ, double minDistance) {
110  const Node& node1 = pair.getFirst();
111  const Node& node2 = pair.getSecond();
112 
113  bool isComp1 = node1.isComposite();
114  bool isComp2 = node2.isComposite();
115 
121  if (isComp1 && isComp2) {
122  if (node1.getSize() > node2.getSize()) {
123  expand(node1, node2, false, priQ, minDistance);
124  return;
125  } else {
126  expand(node2, node1, true, priQ, minDistance);
127  return;
128  }
129  } else if (isComp1) {
130  expand(node1, node2, false, priQ, minDistance);
131  return;
132  } else if (isComp2) {
133  expand(node2, node1, true, priQ, minDistance);
134  return;
135  }
136 
137  throw util::IllegalArgumentException("neither boundable is composite");
138 
139  }
140 
141  void expand(const Node &nodeComposite, const Node &nodeOther, bool isFlipped, PairQueue& priQ,
142  double minDistance) {
143  for (const auto *child = nodeComposite.beginChildren();
144  child < nodeComposite.endChildren(); ++child) {
145  NodePair sp = isFlipped ? NodePair(nodeOther, *child, m_id) : NodePair(*child, nodeOther, m_id);
146 
147  // only add to queue if this pair might contain the closest points
148  if (minDistance == DoubleInfinity || sp.getDistance() < minDistance) {
149  priQ.push(sp);
150  }
151  }
152  }
153 
154  ItemDistance& m_id;
155 };
156 }
157 }
158 }
Basic namespace for all GEOS functionalities.
Definition: Angle.h:26