diff --git a/DataStructures/EdgeBasedNode.h b/DataStructures/EdgeBasedNode.h index 542eda451..654de96de 100644 --- a/DataStructures/EdgeBasedNode.h +++ b/DataStructures/EdgeBasedNode.h @@ -4,6 +4,7 @@ #include "Coordinate.h" struct EdgeBasedNode { + EdgeBasedNode() : id(INT_MAX), lat1(INT_MAX), @@ -16,6 +17,55 @@ struct EdgeBasedNode { 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::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 { return other.id < id; } @@ -48,4 +98,4 @@ struct EdgeBasedNode { bool ignoreInGrid:1; }; -#endif //EDGE_BASED_NODE_H \ No newline at end of file +#endif //EDGE_BASED_NODE_H diff --git a/DataStructures/StaticRTree.h b/DataStructures/StaticRTree.h index 5feab4c7d..62ca89530 100644 --- a/DataStructures/StaticRTree.h +++ b/DataStructures/StaticRTree.h @@ -124,10 +124,10 @@ public: FixedPointCoordinate lower_left (other.min_lat, other.min_lon); return ( - Contains(upper_left) - || Contains(upper_right) - || Contains(lower_right) - || Contains(lower_left) + Contains(upper_left ) || + Contains(upper_right) || + Contains(lower_right) || + Contains(lower_left ) ); } @@ -188,32 +188,32 @@ public: min_max_dist = std::min( min_max_dist, std::max( - ApproximateDistance(location, upper_left ), - ApproximateDistance(location, upper_right) + ApproximateDistance(location, upper_left ), + ApproximateDistance(location, upper_right) ) ); min_max_dist = std::min( min_max_dist, std::max( - ApproximateDistance(location, upper_right), - ApproximateDistance(location, lower_right) + ApproximateDistance(location, upper_right), + ApproximateDistance(location, lower_right) ) ); min_max_dist = std::min( min_max_dist, std::max( - ApproximateDistance(location, lower_right), - ApproximateDistance(location, lower_left ) + ApproximateDistance(location, lower_right), + ApproximateDistance(location, lower_left ) ) ); min_max_dist = std::min( min_max_dist, std::max( - ApproximateDistance(location, lower_left ), - ApproximateDistance(location, upper_left ) + ApproximateDistance(location, lower_left ), + ApproximateDistance(location, upper_left ) ) ); return min_max_dist; @@ -318,8 +318,8 @@ public: uint64_t current_hilbert_value = HilbertCode::GetHilbertNumberForCoordinate(current_centroid); input_wrapper_vector[element_counter].m_hilbert_value = current_hilbert_value; - } + //open leaf file boost::filesystem::ofstream leaf_node_file(leaf_node_filename, std::ios::binary); leaf_node_file.write((char*) &m_element_count, sizeof(uint64_t)); @@ -334,6 +334,7 @@ public: LeafNode current_leaf; 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) { 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; @@ -391,6 +392,7 @@ public: //reverse and renumber tree to have root at index 0 std::reverse(m_search_tree.begin(), m_search_tree.end()); + #pragma omp parallel for schedule(guided) for(uint32_t i = 0; i < m_search_tree.size(); ++i) { TreeNode & current_tree_node = m_search_tree[i]; @@ -435,7 +437,7 @@ public: uint32_t tree_size = 0; 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); tree_node_file.read((char*)&m_search_tree[0], sizeof(TreeNode)*tree_size); tree_node_file.close(); @@ -526,12 +528,10 @@ public: } double current_ratio = 0.; - double current_perpendicular_distance = ComputePerpendicularDistance( + double current_perpendicular_distance = current_edge.ComputePerpendicularDistance( input_coordinate, - FixedPointCoordinate(current_edge.lat1, current_edge.lon1), - FixedPointCoordinate(current_edge.lat2, current_edge.lon2), nearest, - ¤t_ratio + current_ratio ); if( @@ -556,7 +556,7 @@ public: } else if( DoubleEpsilonCompare(current_perpendicular_distance, min_dist) && 1 == abs(current_edge.id - result_phantom_node.edgeBasedNode ) - && CoordinatesAreEquivalent( + && EdgesAreEquivalent( current_start_coordinate, FixedPointCoordinate( current_edge.lat1, @@ -736,7 +736,6 @@ public: return found_a_nearest_edge; } - bool FindPhantomNodeForCoordinate( const FixedPointCoordinate & input_coordinate, PhantomNode & result_phantom_node, @@ -746,7 +745,7 @@ public: bool ignore_tiny_components = (zoom_level <= 14); DataT nearest_edge; - uint32_t io_count = 0; + // uint32_t io_count = 0; uint32_t explored_tree_nodes_count = 0; //SimpleLogger().Write() << "searching for coordinate " << input_coordinate; double min_dist = std::numeric_limits::max(); @@ -758,9 +757,7 @@ public: //initialize queue with root element std::priority_queue traversal_queue; double current_min_dist = m_search_tree[0].minimum_bounding_rectangle.GetMinDist(input_coordinate); - traversal_queue.push( - QueryCandidate(0, current_min_dist) - ); + traversal_queue.push( QueryCandidate(0, current_min_dist) ); BOOST_ASSERT_MSG( std::numeric_limits::epsilon() > (0. - traversal_queue.top().min_dist), @@ -778,7 +775,7 @@ public: if (current_tree_node.child_is_on_disk) { LeafNode 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) { DataT & current_edge = current_leaf_node.objects[i]; if(ignore_tiny_components && current_edge.belongsToTinyComponent) { @@ -788,21 +785,20 @@ public: continue; } - double current_ratio = 0.; - double current_perpendicular_distance = ComputePerpendicularDistance( - input_coordinate, - FixedPointCoordinate(current_edge.lat1, current_edge.lon1), - FixedPointCoordinate(current_edge.lat2, current_edge.lon2), - nearest, - ¤t_ratio + double current_ratio = 0.; + double current_perpendicular_distance = current_edge.ComputePerpendicularDistance( + input_coordinate, + nearest, + current_ratio ); + BOOST_ASSERT( 0 <= current_perpendicular_distance ); if( - current_perpendicular_distance < min_dist - && !DoubleEpsilonCompare( - current_perpendicular_distance, - min_dist - ) + ( current_perpendicular_distance < min_dist ) && + !DoubleEpsilonCompare( + current_perpendicular_distance, + min_dist + ) ) { //found a new minimum min_dist = current_perpendicular_distance; result_phantom_node.edgeBasedNode = current_edge.id; @@ -816,10 +812,10 @@ public: current_end_coordinate.lon = current_edge.lon2; nearest_edge = current_edge; found_a_nearest_edge = true; - } else if( - DoubleEpsilonCompare(current_perpendicular_distance, min_dist) && - 1 == abs(current_edge.id - result_phantom_node.edgeBasedNode ) - && CoordinatesAreEquivalent( + } else + if( DoubleEpsilonCompare(current_perpendicular_distance, min_dist) && + ( 1 == abs(current_edge.id - result_phantom_node.edgeBasedNode ) ) && + EdgesAreEquivalent( current_start_coordinate, FixedPointCoordinate( current_edge.lat1, @@ -853,10 +849,10 @@ public: if( current_min_max_dist < 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; } - if (current_min_dist > min_dist) { //upward pruning + if( current_min_dist > min_dist ) { //upward pruning continue; } traversal_queue.push(QueryCandidate(child_id, current_min_dist)); @@ -866,10 +862,10 @@ public: } const double ratio = (found_a_nearest_edge ? - std::min(1., ApproximateDistance(current_start_coordinate, - result_phantom_node.location)/ApproximateDistance(current_start_coordinate, current_end_coordinate) - ) : 0 - ); + std::min(1., ApproximateDistance(current_start_coordinate, + result_phantom_node.location)/ApproximateDistance(current_start_coordinate, current_end_coordinate) + ) : 0 + ); result_phantom_node.weight1 *= ratio; if(INT_MAX != result_phantom_node.weight2) { result_phantom_node.weight2 *= (1.-ratio); @@ -885,7 +881,6 @@ public: } return found_a_nearest_edge; - } private: @@ -910,55 +905,7 @@ private: thread_local_rtree_stream->read((char *)&result_node, sizeof(LeafNode)); } - inline double ComputePerpendicularDistance( - 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::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 { + inline bool EdgesAreEquivalent(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); }