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