Open-Harbor/osrm-backend
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 3 Wiki Activity

move distance computation from r-tree to element type class

Browse Source
This commit is contained in:
Dennis Luxen 2013-11-22 18:05:47 +01:00
parent 5ef7ea794a
commit 50d6b10be4
2 changed files with 95 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
DataStructures/EdgeBasedNode.h
View File

@ -4,6 +4,7 @@
#include "Coordinate.h" #include "Coordinate.h"
struct EdgeBasedNode { struct EdgeBasedNode {
EdgeBasedNode() : EdgeBasedNode() :
id(INT_MAX), id(INT_MAX),
lat1(INT_MAX), lat1(INT_MAX),
@ -16,6 +17,55 @@ struct EdgeBasedNode {
ignoreInGrid(false) ignoreInGrid(false)
{ } { }
double ComputePerpendicularDistance(
const FixedPointCoordinate& inputPoint,
FixedPointCoordinate & nearest_location,
double & r
) const {
const double x = inputPoint.lat/COORDINATE_PRECISION;
const double y = inputPoint.lon/COORDINATE_PRECISION;
const double a = lat1/COORDINATE_PRECISION;
const double b = lon1/COORDINATE_PRECISION;
const double c = lat2/COORDINATE_PRECISION;
const double d = lon2/COORDINATE_PRECISION;
double p,q,mX,nY;
if(std::fabs(a-c) > std::numeric_limits<double>::epsilon() ){
const double m = (d-b)/(c-a); // slope
// Projection of (x,y) on line joining (a,b) and (c,d)
p = ((x + (m*y)) + (m*m*a - m*b))/(1. + m*m);
q = b + m*(p - a);
} else {
p = c;
q = y;
}
nY = (d*p - c*q)/(a*d - b*c);
mX = (p - nY*a)/c;// These values are actually n/m+n and m/m+n , we need
// not calculate the explicit values of m an n as we
// are just interested in the ratio
if(std::isnan(mX)) {
r = ((lat2 == inputPoint.lat) && (lon2 == inputPoint.lon)) ? 1. : 0.;
} else {
r = mX;
}
if(r<=0.){
nearest_location.lat = lat1;
nearest_location.lon = lon1;
return ((b - y)*(b - y) + (a - x)*(a - x));
// return std::sqrt(((b - y)*(b - y) + (a - x)*(a - x)));
} else if(r >= 1.){
nearest_location.lat = lat2;
nearest_location.lon = lon2;
return ((d - y)*(d - y) + (c - x)*(c - x));
// return std::sqrt(((d - y)*(d - y) + (c - x)*(c - x)));
}
// point lies in between
nearest_location.lat = p*COORDINATE_PRECISION;
nearest_location.lon = q*COORDINATE_PRECISION;
// return std::sqrt((p-x)*(p-x) + (q-y)*(q-y));
return (p-x)*(p-x) + (q-y)*(q-y);
}
bool operator<(const EdgeBasedNode & other) const { bool operator<(const EdgeBasedNode & other) const {
return other.id < id; return other.id < id;
} }
@ -48,4 +98,4 @@ struct EdgeBasedNode {
bool ignoreInGrid:1; bool ignoreInGrid:1;
}; };
#endif //EDGE_BASED_NODE_H #endif //EDGE_BASED_NODE_H

141
DataStructures/StaticRTree.h
View File

@ -124,10 +124,10 @@ public:
FixedPointCoordinate lower_left (other.min_lat, other.min_lon); FixedPointCoordinate lower_left (other.min_lat, other.min_lon);
return ( return (
Contains(upper_left) Contains(upper_left ) ||
|| Contains(upper_right) Contains(upper_right) ||
|| Contains(lower_right) Contains(lower_right) ||
|| Contains(lower_left) Contains(lower_left )
); );
} }
@ -188,32 +188,32 @@ public:
min_max_dist = std::min( min_max_dist = std::min(
min_max_dist, min_max_dist,
std::max( std::max(
ApproximateDistance(location, upper_left ), ApproximateDistance(location, upper_left ),
ApproximateDistance(location, upper_right) ApproximateDistance(location, upper_right)
) )
); );
min_max_dist = std::min( min_max_dist = std::min(
min_max_dist, min_max_dist,
std::max( std::max(
ApproximateDistance(location, upper_right), ApproximateDistance(location, upper_right),
ApproximateDistance(location, lower_right) ApproximateDistance(location, lower_right)
) )
); );
min_max_dist = std::min( min_max_dist = std::min(
min_max_dist, min_max_dist,
std::max( std::max(
ApproximateDistance(location, lower_right), ApproximateDistance(location, lower_right),
ApproximateDistance(location, lower_left ) ApproximateDistance(location, lower_left )
) )
); );
min_max_dist = std::min( min_max_dist = std::min(
min_max_dist, min_max_dist,
std::max( std::max(
ApproximateDistance(location, lower_left ), ApproximateDistance(location, lower_left ),
ApproximateDistance(location, upper_left ) ApproximateDistance(location, upper_left )
) )
); );
return min_max_dist; return min_max_dist;
@ -318,8 +318,8 @@ public:
uint64_t current_hilbert_value = HilbertCode::GetHilbertNumberForCoordinate(current_centroid); uint64_t current_hilbert_value = HilbertCode::GetHilbertNumberForCoordinate(current_centroid);
input_wrapper_vector[element_counter].m_hilbert_value = current_hilbert_value; input_wrapper_vector[element_counter].m_hilbert_value = current_hilbert_value;
} }
//open leaf file //open leaf file
boost::filesystem::ofstream leaf_node_file(leaf_node_filename, std::ios::binary); boost::filesystem::ofstream leaf_node_file(leaf_node_filename, std::ios::binary);
leaf_node_file.write((char*) &m_element_count, sizeof(uint64_t)); leaf_node_file.write((char*) &m_element_count, sizeof(uint64_t));
@ -334,6 +334,7 @@ public:
LeafNode current_leaf; LeafNode current_leaf;
TreeNode current_node; TreeNode current_node;
//SimpleLogger().Write() << "reading " << tree_size << " tree nodes in " << (sizeof(TreeNode)*tree_size) << " bytes";
for(uint32_t current_element_index = 0; RTREE_LEAF_NODE_SIZE > current_element_index; ++current_element_index) { for(uint32_t current_element_index = 0; RTREE_LEAF_NODE_SIZE > current_element_index; ++current_element_index) {
if(m_element_count > (processed_objects_count + current_element_index)) { if(m_element_count > (processed_objects_count + current_element_index)) {
uint32_t index_of_next_object = input_wrapper_vector[processed_objects_count + current_element_index].m_array_index; uint32_t index_of_next_object = input_wrapper_vector[processed_objects_count + current_element_index].m_array_index;
@ -391,6 +392,7 @@ public:
//reverse and renumber tree to have root at index 0 //reverse and renumber tree to have root at index 0
std::reverse(m_search_tree.begin(), m_search_tree.end()); std::reverse(m_search_tree.begin(), m_search_tree.end());
#pragma omp parallel for schedule(guided) #pragma omp parallel for schedule(guided)
for(uint32_t i = 0; i < m_search_tree.size(); ++i) { for(uint32_t i = 0; i < m_search_tree.size(); ++i) {
TreeNode & current_tree_node = m_search_tree[i]; TreeNode & current_tree_node = m_search_tree[i];
@ -435,7 +437,7 @@ public:
uint32_t tree_size = 0; uint32_t tree_size = 0;
tree_node_file.read((char*)&tree_size, sizeof(uint32_t)); tree_node_file.read((char*)&tree_size, sizeof(uint32_t));
//SimpleLogger().Write() << "reading " << tree_size << " tree nodes in " << (sizeof(TreeNode)*tree_size) << " bytes";
m_search_tree.resize(tree_size); m_search_tree.resize(tree_size);
tree_node_file.read((char*)&m_search_tree[0], sizeof(TreeNode)*tree_size); tree_node_file.read((char*)&m_search_tree[0], sizeof(TreeNode)*tree_size);
tree_node_file.close(); tree_node_file.close();
@ -526,12 +528,10 @@ public:
} }
double current_ratio = 0.; double current_ratio = 0.;
double current_perpendicular_distance = ComputePerpendicularDistance( double current_perpendicular_distance = current_edge.ComputePerpendicularDistance(
input_coordinate, input_coordinate,
FixedPointCoordinate(current_edge.lat1, current_edge.lon1),
FixedPointCoordinate(current_edge.lat2, current_edge.lon2),
nearest, nearest,
&current_ratio current_ratio
); );
if( if(
@ -556,7 +556,7 @@ public:
} else if( } else if(
DoubleEpsilonCompare(current_perpendicular_distance, min_dist) && DoubleEpsilonCompare(current_perpendicular_distance, min_dist) &&
1 == abs(current_edge.id - result_phantom_node.edgeBasedNode ) 1 == abs(current_edge.id - result_phantom_node.edgeBasedNode )
&& CoordinatesAreEquivalent( && EdgesAreEquivalent(
current_start_coordinate, current_start_coordinate,
FixedPointCoordinate( FixedPointCoordinate(
current_edge.lat1, current_edge.lat1,
@ -736,7 +736,6 @@ public:
return found_a_nearest_edge; return found_a_nearest_edge;
} }
bool FindPhantomNodeForCoordinate( bool FindPhantomNodeForCoordinate(
const FixedPointCoordinate & input_coordinate, const FixedPointCoordinate & input_coordinate,
PhantomNode & result_phantom_node, PhantomNode & result_phantom_node,
@ -746,7 +745,7 @@ public:
bool ignore_tiny_components = (zoom_level <= 14); bool ignore_tiny_components = (zoom_level <= 14);
DataT nearest_edge; DataT nearest_edge;
uint32_t io_count = 0; // uint32_t io_count = 0;
uint32_t explored_tree_nodes_count = 0; uint32_t explored_tree_nodes_count = 0;
//SimpleLogger().Write() << "searching for coordinate " << input_coordinate; //SimpleLogger().Write() << "searching for coordinate " << input_coordinate;
double min_dist = std::numeric_limits<double>::max(); double min_dist = std::numeric_limits<double>::max();
@ -758,9 +757,7 @@ public:
//initialize queue with root element //initialize queue with root element
std::priority_queue<QueryCandidate> traversal_queue; std::priority_queue<QueryCandidate> traversal_queue;
double current_min_dist = m_search_tree[0].minimum_bounding_rectangle.GetMinDist(input_coordinate); double current_min_dist = m_search_tree[0].minimum_bounding_rectangle.GetMinDist(input_coordinate);
traversal_queue.push( traversal_queue.push( QueryCandidate(0, current_min_dist) );
QueryCandidate(0, current_min_dist)
);
BOOST_ASSERT_MSG( BOOST_ASSERT_MSG(
std::numeric_limits<double>::epsilon() > (0. - traversal_queue.top().min_dist), std::numeric_limits<double>::epsilon() > (0. - traversal_queue.top().min_dist),
@ -778,7 +775,7 @@ public:
if (current_tree_node.child_is_on_disk) { if (current_tree_node.child_is_on_disk) {
LeafNode current_leaf_node; LeafNode current_leaf_node;
LoadLeafFromDisk(current_tree_node.children[0], current_leaf_node); LoadLeafFromDisk(current_tree_node.children[0], current_leaf_node);
++io_count; // ++io_count;
for(uint32_t i = 0; i < current_leaf_node.object_count; ++i) { for(uint32_t i = 0; i < current_leaf_node.object_count; ++i) {
DataT & current_edge = current_leaf_node.objects[i]; DataT & current_edge = current_leaf_node.objects[i];
if(ignore_tiny_components && current_edge.belongsToTinyComponent) { if(ignore_tiny_components && current_edge.belongsToTinyComponent) {
@ -788,21 +785,20 @@ public:
continue; continue;
} }
double current_ratio = 0.; double current_ratio = 0.;
double current_perpendicular_distance = ComputePerpendicularDistance( double current_perpendicular_distance = current_edge.ComputePerpendicularDistance(
input_coordinate, input_coordinate,
FixedPointCoordinate(current_edge.lat1, current_edge.lon1), nearest,
FixedPointCoordinate(current_edge.lat2, current_edge.lon2), current_ratio
nearest,
&current_ratio
); );
BOOST_ASSERT( 0 <= current_perpendicular_distance );
if( if(
current_perpendicular_distance < min_dist ( current_perpendicular_distance < min_dist ) &&
&& !DoubleEpsilonCompare( !DoubleEpsilonCompare(
current_perpendicular_distance, current_perpendicular_distance,
min_dist min_dist
) )
) { //found a new minimum ) { //found a new minimum
min_dist = current_perpendicular_distance; min_dist = current_perpendicular_distance;
result_phantom_node.edgeBasedNode = current_edge.id; result_phantom_node.edgeBasedNode = current_edge.id;
@ -816,10 +812,10 @@ public:
current_end_coordinate.lon = current_edge.lon2; current_end_coordinate.lon = current_edge.lon2;
nearest_edge = current_edge; nearest_edge = current_edge;
found_a_nearest_edge = true; found_a_nearest_edge = true;
} else if( } else
DoubleEpsilonCompare(current_perpendicular_distance, min_dist) && if( DoubleEpsilonCompare(current_perpendicular_distance, min_dist) &&
1 == abs(current_edge.id - result_phantom_node.edgeBasedNode ) ( 1 == abs(current_edge.id - result_phantom_node.edgeBasedNode ) ) &&
&& CoordinatesAreEquivalent( EdgesAreEquivalent(
current_start_coordinate, current_start_coordinate,
FixedPointCoordinate( FixedPointCoordinate(
current_edge.lat1, current_edge.lat1,
@ -853,10 +849,10 @@ public:
if( current_min_max_dist < min_max_dist ) { if( current_min_max_dist < min_max_dist ) {
min_max_dist = current_min_max_dist; min_max_dist = current_min_max_dist;
} }
if (current_min_dist > min_max_dist) { if( current_min_dist > min_max_dist ) {
continue; continue;
} }
if (current_min_dist > min_dist) { //upward pruning if( current_min_dist > min_dist ) { //upward pruning
continue; continue;
} }
traversal_queue.push(QueryCandidate(child_id, current_min_dist)); traversal_queue.push(QueryCandidate(child_id, current_min_dist));
@ -866,10 +862,10 @@ public:
} }
const double ratio = (found_a_nearest_edge ? const double ratio = (found_a_nearest_edge ?
std::min(1., ApproximateDistance(current_start_coordinate, std::min(1., ApproximateDistance(current_start_coordinate,
result_phantom_node.location)/ApproximateDistance(current_start_coordinate, current_end_coordinate) result_phantom_node.location)/ApproximateDistance(current_start_coordinate, current_end_coordinate)
) : 0 ) : 0
); );
result_phantom_node.weight1 *= ratio; result_phantom_node.weight1 *= ratio;
if(INT_MAX != result_phantom_node.weight2) { if(INT_MAX != result_phantom_node.weight2) {
result_phantom_node.weight2 *= (1.-ratio); result_phantom_node.weight2 *= (1.-ratio);
@ -885,7 +881,6 @@ public:
} }
return found_a_nearest_edge; return found_a_nearest_edge;
} }
private: private:
@ -910,55 +905,7 @@ private:
thread_local_rtree_stream->read((char *)&result_node, sizeof(LeafNode)); thread_local_rtree_stream->read((char *)&result_node, sizeof(LeafNode));
} }
inline double ComputePerpendicularDistance( inline bool EdgesAreEquivalent(const FixedPointCoordinate & a, const FixedPointCoordinate & b, const FixedPointCoordinate & c, const FixedPointCoordinate & d) const {
const FixedPointCoordinate& inputPoint,
const FixedPointCoordinate& source,
const FixedPointCoordinate& target,
FixedPointCoordinate& nearest, double *r) const {
const double x = inputPoint.lat/COORDINATE_PRECISION;
const double y = inputPoint.lon/COORDINATE_PRECISION;
const double a = source.lat/COORDINATE_PRECISION;
const double b = source.lon/COORDINATE_PRECISION;
const double c = target.lat/COORDINATE_PRECISION;
const double d = target.lon/COORDINATE_PRECISION;
double p,q,mX,nY;
if(std::fabs(a-c) > std::numeric_limits<double>::epsilon() ){
const double m = (d-b)/(c-a); // slope
// Projection of (x,y) on line joining (a,b) and (c,d)
p = ((x + (m*y)) + (m*m*a - m*b))/(1. + m*m);
q = b + m*(p - a);
} else {
p = c;
q = y;
}
nY = (d*p - c*q)/(a*d - b*c);
mX = (p - nY*a)/c;// These values are actually n/m+n and m/m+n , we need
// not calculate the explicit values of m an n as we
// are just interested in the ratio
if(std::isnan(mX)) {
*r = (target == inputPoint) ? 1. : 0.;
} else {
*r = mX;
}
if(*r<=0.){
nearest.lat = source.lat;
nearest.lon = source.lon;
return ((b - y)*(b - y) + (a - x)*(a - x));
// return std::sqrt(((b - y)*(b - y) + (a - x)*(a - x)));
} else if(*r >= 1.){
nearest.lat = target.lat;
nearest.lon = target.lon;
return ((d - y)*(d - y) + (c - x)*(c - x));
// return std::sqrt(((d - y)*(d - y) + (c - x)*(c - x)));
}
// point lies in between
nearest.lat = p*COORDINATE_PRECISION;
nearest.lon = q*COORDINATE_PRECISION;
// return std::sqrt((p-x)*(p-x) + (q-y)*(q-y));
return (p-x)*(p-x) + (q-y)*(q-y);
}
inline bool CoordinatesAreEquivalent(const FixedPointCoordinate & a, const FixedPointCoordinate & b, const FixedPointCoordinate & c, const FixedPointCoordinate & d) const {
return (a == b && c == d) || (a == c && b == d) || (a == d && b == c); return (a == b && c == d) || (a == c && b == d) || (a == d && b == c);
} }