|
| 1 | +#include "spatialindexpoints.h" |
| 2 | + |
| 3 | +#include <vector> |
| 4 | +#include <boost/geometry.hpp> |
| 5 | +#include <boost/geometry/index/rtree.hpp> |
| 6 | + |
| 7 | +#include "domain/pointset.h" |
| 8 | +#include "domain/application.h" |
| 9 | + |
| 10 | +namespace bg = boost::geometry; |
| 11 | +namespace bgi = boost::geometry::index; |
| 12 | + |
| 13 | +typedef bg::model::point<double, 3, bg::cs::cartesian> Point3D; |
| 14 | +typedef bg::model::box<Point3D> Box; |
| 15 | +typedef std::pair<Box, size_t> Value; |
| 16 | + |
| 17 | +// create the R* variant of the rtree |
| 18 | +bgi::rtree< Value, bgi::rstar<16> > rtree; |
| 19 | + |
| 20 | +//the query pointset |
| 21 | +PointSet* pointset = nullptr; |
| 22 | + |
| 23 | +//the GEO-EAS columns of the X, Y, Z variables of the point set file |
| 24 | +int iX; |
| 25 | +int iY; |
| 26 | +int iZ; |
| 27 | + |
| 28 | +void setPointSet( PointSet* ps ){ |
| 29 | + pointset = ps; |
| 30 | + //loads the PointSet data. |
| 31 | + ps->loadData(); |
| 32 | + //get the GEO-EAS indexes -1 for the X, Y and Z coordinates |
| 33 | + iX = ps->getXindex() - 1; |
| 34 | + iY = ps->getYindex() - 1; |
| 35 | + iZ = ps->getZindex() - 1; |
| 36 | +} |
| 37 | + |
| 38 | +SpatialIndexPoints::SpatialIndexPoints() |
| 39 | +{ |
| 40 | + |
| 41 | +} |
| 42 | + |
| 43 | +void SpatialIndexPoints::fill(PointSet *ps, double tolerance) |
| 44 | +{ |
| 45 | + //first clear the index. |
| 46 | + clear(); |
| 47 | + |
| 48 | + setPointSet( ps ); |
| 49 | + |
| 50 | + //for each data line... |
| 51 | + uint totlines = ps->getDataLineCount(); |
| 52 | + for( uint iLine = 0; iLine < totlines; ++iLine){ |
| 53 | + //...make a Point3D for the index |
| 54 | + double x = ps->data( iLine, iX ); |
| 55 | + double y = ps->data( iLine, iY ); |
| 56 | + double z = 0.0; //put 2D data in the z==0.0 plane |
| 57 | + if( iZ >= 0 ) |
| 58 | + z = ps->data( iLine, iZ ); |
| 59 | + //make a bounding box around the point. |
| 60 | + Box box( Point3D(x-tolerance, y-tolerance, z-tolerance), |
| 61 | + Point3D(x+tolerance, y+tolerance, z+tolerance)); |
| 62 | + //insert the box representing the point into the spatial index. |
| 63 | + rtree.insert( std::make_pair(box, iLine) ); |
| 64 | + } |
| 65 | +} |
| 66 | + |
| 67 | +QList<uint> SpatialIndexPoints::getNearest(uint index, uint n) |
| 68 | +{ |
| 69 | + QList<uint> result; |
| 70 | + |
| 71 | + //get the location of the point. |
| 72 | + double x = pointset->data( index, iX ); |
| 73 | + double y = pointset->data( index, iY ); |
| 74 | + double z = 0.0; //put 2D data in the z==0.0 plane |
| 75 | + if( iZ >= 0 ) |
| 76 | + z = pointset->data( index, iZ ); |
| 77 | + |
| 78 | + // find n nearest values to a point |
| 79 | + std::vector<Value> result_n; |
| 80 | + rtree.query(bgi::nearest(Point3D(x, y, z), n), std::back_inserter(result_n)); |
| 81 | + |
| 82 | + // collect the point indexes |
| 83 | + std::vector<Value>::iterator it = result_n.begin(); |
| 84 | + for(; it != result_n.end(); ++it){ |
| 85 | + //do not return itself |
| 86 | + if( index != (*it).second ) |
| 87 | + result.push_back( (*it).second ); |
| 88 | + } |
| 89 | + |
| 90 | + //return the point indexes |
| 91 | + return result; |
| 92 | +} |
| 93 | + |
| 94 | +QList<uint> SpatialIndexPoints::getNearestWithin(uint index, uint n, double distance ) |
| 95 | +{ |
| 96 | + QList<uint> result; |
| 97 | + |
| 98 | + //get the location of the query point. |
| 99 | + double qx = pointset->data( index, iX ); |
| 100 | + double qy = pointset->data( index, iY ); |
| 101 | + double qz = 0.0; //put 2D data in the z==0.0 plane |
| 102 | + if( iZ >= 0 ) |
| 103 | + qz = pointset->data( index, iZ ); |
| 104 | + Point3D qPoint(qx, qy, qz); |
| 105 | + |
| 106 | + //get the n-nearest points |
| 107 | + QList<uint> nearestSamples = SpatialIndexPoints::getNearest( index, n ); |
| 108 | + |
| 109 | + //test the distance to each of the n-nearest points |
| 110 | + QList<uint>::iterator it = nearestSamples.begin(); |
| 111 | + for(; it != nearestSamples.end(); ++it){ |
| 112 | + //get the location of a near point. |
| 113 | + uint nIndex = *it; |
| 114 | + double nx = pointset->data( nIndex, iX ); |
| 115 | + double ny = pointset->data( nIndex, iY ); |
| 116 | + double nz = 0.0; //put 2D data in the z==0.0 plane |
| 117 | + if( iZ >= 0 ) |
| 118 | + nz = pointset->data( nIndex, iZ ); |
| 119 | + //compute the distance between the query point and a nearest point |
| 120 | + double dist = boost::geometry::distance( qPoint, Point3D(nx, ny, nz) ); |
| 121 | + if( dist < distance ){ |
| 122 | + result.push_back( nIndex ); |
| 123 | + } |
| 124 | + } |
| 125 | + return result; |
| 126 | +} |
| 127 | + |
| 128 | +void SpatialIndexPoints::clear() |
| 129 | +{ |
| 130 | + rtree.clear(); |
| 131 | + pointset = nullptr; |
| 132 | +} |
0 commit comments