i need store 1 billion spatial polygons (specified using minimum bounding rectangles) in quad-tree. doing so, wrote following code. however, turns out 1 billion points code running slowly. there way may improve code may run bit faster. if yes, can please same
//--------------------------------------------------------------------------- struct mbr { double xright, xleft, ybottom, ytop; mbr *zero,*first,*second,*third; unsigned level=0; vector<unsigned> result; //contains resulting intersecting spatial ids }; bool intersects(mbr& spatialid,mbr& mbr) { if (mbr.ybottom > spatialid.ytop || mbr.ytop < spatialid.ybottom) return false; if (mbr.xleft > spatialid.xright || mbr.xright < spatialid.xleft) return false; return true; } //--------------------------------------------------------------------------- bool contains(mbr& spatialid,mbr& mbr) { if (mbr.ybottom > spatialid.ybottom || mbr.ytop < spatialid.ytop) return false; if (mbr.xleft > spatialid.xleft || mbr.xright < spatialid.xright) return false; return true; } //--------------------------------------------------------------------------- bool touches(mbr& spatialid,mbr& mbr) { if ( (mbr.ybottom >= spatialid.ybottom + std::numeric_limits<double>::epsilon() && mbr.ybottom <= spatialid.ybottom - std::numeric_limits<double>::epsilon()) || (mbr.ytop >= spatialid.ytop + std::numeric_limits<double>::epsilon() && mbr.ytop <= spatialid.ytop - std::numeric_limits<double>::epsilon())) return true; if ( (mbr.xleft >= spatialid.xleft + std::numeric_limits<double>::epsilon() && mbr.xleft <= spatialid.xleft - std::numeric_limits<double>::epsilon()) || (mbr.xright >= spatialid.xright + std::numeric_limits<double>::epsilon() && mbr.xright <= spatialid.xright - std::numeric_limits<double>::epsilon())) return true; return false; } //--------------------------------------------------------------------------- mbr mbr1,mbr2,mbr3,mbr4; vector<unsigned> spatialids; //contain 1 billion spatial identifiers intersected mbr1, mbr2, mbr3, mbr4 //mbr1, mbr2, mbr3, mbr4 again specified using minimum bounding rectangles stack<mbr**> stackquadtree; mbr *root=new mbr(); (*root).ybottom=-90; (*root).ytop=90; (*root).xleft=-180; (*root).xright=180; (*root).level=0; stackquadtree.push(&root); while(!stackquadtree.empty()) { mbr** node=&(*stackquadtree.front()); if((*node)->level==50) break; (*node)->zero=new mbr(); (*node)->first=new mbr(); (*node)->second=new mbr(); (*node)->third=new mbr(); (*node)->zero->ybottom=(*node)->ybottom; (*node)->zero->ytop=((*node)->ybottom+(*node)->ytop)/2; (*node)->zero->xleft=(*node)->xleft; (*node)->zero->xright=((*node)->xleft+(*node)->xright)/2; (*node)->zero->level=(*node)->level+1; (*node)->first->ybottom=((*node)->ybottom+(*node)->ytop)/2; (*node)->first->ytop=(*node)->ytop; (*node)->first->xleft=(*node)->xleft; (*node)->first->xright=((*node)->xleft+(*node)->xright)/2; (*node)->first->level=(*node)->level+1; (*node)->second->ybottom=(*node)->ybottom; (*node)->second->ytop=((*node)->ybottom+(*node)->ytop)/2; (*node)->second->xleft=((*node)->xleft+(*node)->xright)/2; (*node)->second->xright=(*node)->xright; (*node)->second->level=(*node)->level+1; (*node)->third->ybottom=((*node)->ybottom+(*node)->ytop)/2; (*node)->third->ytop=(*node)->ytop; (*node)->third->xleft=((*node)->xleft+(*node)->xright)/2; (*node)->third->xright=(*node)->xright; (*node)->third->level=(*node)->level+1; mbr* node=*stackquadtree.top(); stackquadtree.pop(); for(vector<mbr>::iterator itspatialid=spatialids.begin(),lspatialid=spatialids.end();itspatialid!=lspatialid;++itspatialid) { if(intersects((*itspatialid),&(*node)->zero)||contains((*itspatialid),&(*node)->zero)||touches((*itspatialid),&(*node)->zero)) { (*node)->zero->result.push_back((*itspatialid)); stackquadtree.push(*(*node)->zero); } if(intersects((*itspatialid),&(*node)->first)||contains((*itspatialid),&(*node)->first)||touches((*itspatialid),&(*node)->first)) { (*node)->first->result.push_back((*itspatialid)); stackquadtree.push(*(*node)->first); } if(intersects((*itspatialid),&(*node)->second)||contains((*itspatialid),&(*node)->second)||touches((*itspatialid),&(*node)->second)) { (*node)->second->result.push_back((*itspatialid)); stackquadtree.push(*(*node)->second); } if(intersects((*itspatialid),&(*node)->third)||contains((*itspatialid),&(*node)->third)||touches((*itspatialid),&(*node)->third)) { (*node)->third->result.push_back((*itspatialid)); stackquadtree.push(*(*node)->third); } } }
log2(1 billion) approximately 30
in order reduce operations count, consider data structure more rapidly gets right neighborhood.
for example, if know objects located within grid 10km x 10km, consider breaking down 10x10 grid (1kmx1km) then, object falls within particular grid can jump 1% of objects (assuming evenly distributed).
obviously, can have recursive structure, don't need more 1 or 2 10x10 splits. why not try top level, each subgrid quadtree?
i notice have odd construction:
mbr.ytop >= spatialid.ytop + std::numeric_limits<double>::epsilon() i think can simplify to:
spatialid.ytop < mbr.ytop that kind of definition of floating point resolution, isn't it? if i'm right, should speed code lot that's constant factor.
class grid { private: stackquadtree qt[100]; // 10x10 array of quadtrees public: grid(double xmin, double xmax, double ymin, double ymax) { // store 10x10 grid in grid object gridsizex = 1/10 of size of grid in x gridsizey = 1/10 of size of grid in y } stackquadtree& findgrid(mbr& mbr) { int = (mbr.xleft - xmin) / gridsizex; int j = (mbr.ytop - ymin) / gridsizey; return qt[i*10+j]; } void add(mbr& mbr) { // add right quadtree , add object in } }; 
Comments
Post a Comment