src/operation/distance/DistanceOp.cpp

   1 /**********************************************************************
   2  *
   3  * GEOS - Geometry Engine Open Source
   4  * http://geos.osgeo.org
   5  *
   6  * Copyright (C) 2011 Sandro Santilli <strk@keybit.net>
   7  * Copyright (C) 2006 Refractions Research Inc.
   8  * Copyright (C) 2001-2002 Vivid Solutions Inc.
   9  *
  10  * This is free software; you can redistribute and/or modify it under
  11  * the terms of the GNU Lesser General Public Licence as published
  12  * by the Free Software Foundation.
  13  * See the COPYING file for more information.
  14  *
  15  **********************************************************************
  16  *
  17  * Last port: operation/distance/DistanceOp.java r335 (JTS-1.12-)
  18  *
  19  **********************************************************************/
  20
  21 #include <geos/operation/distance/DistanceOp.h>
  22 #include <geos/operation/distance/GeometryLocation.h>
  23 #include <geos/operation/distance/ConnectedElementLocationFilter.h>
  24 #include <geos/algorithm/PointLocator.h>
  25 #include <geos/algorithm/CGAlgorithms.h>
  26 #include <geos/geom/Coordinate.h>
  27 #include <geos/geom/CoordinateSequence.h>
  28 #include <geos/geom/CoordinateArraySequence.h>
  29 #include <geos/geom/LineString.h>
  30 #include <geos/geom/Point.h>
  31 #include <geos/geom/Polygon.h>
  32 #include <geos/geom/Envelope.h>
  33 #include <geos/geom/LineSegment.h>
  34 #include <geos/geom/util/PolygonExtracter.h>
  35 #include <geos/geom/util/LinearComponentExtracter.h>
  36 #include <geos/geom/util/PointExtracter.h>
  37 #include <geos/util/IllegalArgumentException.h>
  38
  39 #include <vector>
  40 #include <iostream>
  41
  42 #ifndef GEOS_DEBUG
  43 #define GEOS_DEBUG 0
  44 #endif
  45
  46 using namespace std;
  47 using namespace geos::geom;
  48 //using namespace geos::algorithm;
  49
  50 namespace geos {
  51 namespace operation { // geos.operation
  52 namespace distance { // geos.operation.distance
  53
  54 using namespace geom;
  55 //using namespace geom::util;
  56
  57 /*public static (deprecated)*/
  58 double
  59 DistanceOp::distance(const Geometry *g0, const Geometry *g1)
  60 {
  61         DistanceOp distOp(g0,g1);
  62         return distOp.distance();
  63 }
  64
  65 /*public static*/
  66 double
  67 DistanceOp::distance(const Geometry& g0, const Geometry& g1)
  68 {
  69         DistanceOp distOp(g0,g1);
  70         return distOp.distance();
  71 }
  72
  73 /*public static deprecated*/
  74 CoordinateSequence*
  75 DistanceOp::closestPoints(const Geometry *g0, const Geometry *g1)
  76 {
  77         DistanceOp distOp(g0,g1);
  78         return distOp.nearestPoints();
  79 }
  80
  81 /*public static*/
  82 CoordinateSequence*
  83 DistanceOp::nearestPoints(const Geometry *g0, const Geometry *g1)
  84 {
  85         DistanceOp distOp(g0,g1);
  86         return distOp.nearestPoints();
  87 }
  88
  89 DistanceOp::DistanceOp(const Geometry *g0, const Geometry *g1):
  90         geom(2),
  91         terminateDistance(0.0),
  92         minDistanceLocation(0),
  93         minDistance(DoubleMax)
  94 {
  95         geom[0] = g0;
  96         geom[1] = g1;
  97 }
  98
  99 DistanceOp::DistanceOp(const Geometry& g0, const Geometry& g1):
 100         geom(2),
 101         terminateDistance(0.0),
 102         minDistanceLocation(0),
 103         minDistance(DoubleMax)
 104 {
 105         geom[0] = &g0;
 106         geom[1] = &g1;
 107 }
 108
 109 DistanceOp::DistanceOp(const Geometry& g0, const Geometry& g1, double tdist)
 110         :
 111         geom(2),
 112         terminateDistance(tdist),
 113         minDistanceLocation(0),
 114         minDistance(DoubleMax)
 115 {
 116         geom[0] = &g0;
 117         geom[1] = &g1;
 118 }
 119
 120 DistanceOp::~DistanceOp()
 121 {
 122         size_t i;
 123         for (i=0; i<newCoords.size(); i++) delete newCoords[i];
 124         if ( minDistanceLocation )
 125         {
 126                 for (i=0; i<minDistanceLocation->size(); i++)
 127                 {
 128                         delete (*minDistanceLocation)[i];
 129                 }
 130                 delete minDistanceLocation;
 131         }
 132 }
 133
 134 /**
 135 * Report the distance between the closest points on the input geometries.
 136 *
 137 * @return the distance between the geometries
 138 */
 139 double
 140 DistanceOp::distance()
 141 {
 142         using geos::util::IllegalArgumentException;
 143
 144         if ( geom[0] == 0 || geom[1] == 0 )
 145                 throw IllegalArgumentException("null geometries are not supported");
 146         if ( geom[0]->isEmpty() || geom[1]->isEmpty() ) return 0.0;
 147         computeMinDistance();
 148         return minDistance;
 149 }
 150
 151 /* public */
 152 CoordinateSequence*
 153 DistanceOp::closestPoints()
 154 {
 155         return nearestPoints();
 156 }
 157
 158
 159 /* public */
 160 CoordinateSequence*
 161 DistanceOp::nearestPoints()
 162 {
 163         // lazily creates minDistanceLocation
 164         computeMinDistance();
 165
 166         assert(0 != minDistanceLocation);
 167         std::vector<GeometryLocation*>& locs = *minDistanceLocation;
 168
 169         // Empty input geometries result in this behaviour
 170         if ( locs[0] == 0 || locs[1] == 0 )
 171         {
 172                 // either both or none are set..
 173                 assert(locs[0] == 0 && locs[1] == 0);
 174
 175                 return 0;
 176         }
 177
 178         GeometryLocation* loc0 = locs[0];
 179         GeometryLocation* loc1 = locs[1];
 180         const Coordinate& c0 = loc0->getCoordinate();
 181         const Coordinate& c1 = loc1->getCoordinate();
 182
 183         CoordinateSequence* nearestPts = new CoordinateArraySequence();
 184         nearestPts->add(c0);
 185         nearestPts->add(c1);
 186
 187         return nearestPts;
 188 }
 189
 190 /*private, unused!*/
 191 vector<GeometryLocation*>* DistanceOp::nearestLocations(){
 192         computeMinDistance();
 193         return minDistanceLocation;
 194 }
 195
 196 void
 197 DistanceOp::updateMinDistance(vector<GeometryLocation*>& locGeom, bool flip)
 198 {
 199         assert(minDistanceLocation);
 200
 201         // if not set then don't update
 202         if (locGeom[0]==NULL)
 203         {
 204 #if GEOS_DEBUG
 205 std::cerr << "updateMinDistance called with loc[0] == null and loc[1] == " << locGeom[1] << std::endl;
 206 #endif
 207                 assert(locGeom[1] == NULL);
 208                 return;
 209         }
 210
 211         delete (*minDistanceLocation)[0];
 212         delete (*minDistanceLocation)[1];
 213         if (flip) {
 214                 (*minDistanceLocation)[0]=locGeom[1];
 215                 (*minDistanceLocation)[1]=locGeom[0];
 216         } else {
 217                 (*minDistanceLocation)[0]=locGeom[0];
 218                 (*minDistanceLocation)[1]=locGeom[1];
 219         }
 220 }
 221
 222 /*private*/
 223 void
 224 DistanceOp::computeMinDistance()
 225 {
 226         // only compute once!
 227         if (minDistanceLocation) return;
 228
 229 #if GEOS_DEBUG
 230         std::cerr << "---Start: " << geom[0]->toString() << " - " << geom[1]->toString() << std::endl;
 231 #endif
 232
 233         minDistanceLocation = new vector<GeometryLocation*>(2);
 234
 235         computeContainmentDistance();
 236
 237         if (minDistance <= terminateDistance)
 238         {
 239                 return;
 240         }
 241
 242         computeFacetDistance();
 243
 244 #if GEOS_DEBUG
 245         std::cerr << "---End " << std::endl;
 246 #endif
 247 }
 248
 249 /*private*/
 250 void
 251 DistanceOp::computeContainmentDistance()
 252 {
 253         using geom::util::PolygonExtracter;
 254
 255         Polygon::ConstVect polys1;
 256         PolygonExtracter::getPolygons(*(geom[1]), polys1);
 257
 258
 259 #if GEOS_DEBUG
 260         std::cerr << "PolygonExtracter found " << polys1.size() << " polygons in geometry 2" << std::endl;
 261 #endif
 262
 263         // NOTE:
 264         // Expected to fill minDistanceLocation items
 265         // if minDistance <= terminateDistance
 266
 267         vector<GeometryLocation*> *locPtPoly = new vector<GeometryLocation*>(2);
 268         // test if either geometry has a vertex inside the other
 269         if ( ! polys1.empty() )
 270         {
 271                 vector<GeometryLocation*> *insideLocs0 =
 272                         ConnectedElementLocationFilter::getLocations(geom[0]);
 273                 computeInside(insideLocs0, polys1, locPtPoly);
 274                 if (minDistance <= terminateDistance) {
 275                         assert( (*locPtPoly)[0] );
 276                         assert( (*locPtPoly)[1] );
 277                         (*minDistanceLocation)[0] = (*locPtPoly)[0];
 278                         (*minDistanceLocation)[1] = (*locPtPoly)[1];
 279                         delete locPtPoly;
 280                         for (size_t i=0; i<insideLocs0->size(); i++)
 281                         {
 282                                 GeometryLocation *l = (*insideLocs0)[i];
 283                                 if ( l != (*minDistanceLocation)[0] &&
 284                                         l != (*minDistanceLocation)[1] )
 285                                 {
 286                                         delete l;
 287                                 }
 288                         }
 289                         delete insideLocs0;
 290                         return;
 291                 }
 292                 for (size_t i=0; i<insideLocs0->size(); i++)
 293                         delete (*insideLocs0)[i];
 294                 delete insideLocs0;
 295         }
 296
 297         Polygon::ConstVect polys0;
 298         PolygonExtracter::getPolygons(*(geom[0]), polys0);
 299
 300 #if GEOS_DEBUG
 301         std::cerr << "PolygonExtracter found " << polys0.size() << " polygons in geometry 1" << std::endl;
 302 #endif
 303
 304
 305         if ( ! polys0.empty() )
 306         {
 307                 vector<GeometryLocation*> *insideLocs1 = ConnectedElementLocationFilter::getLocations(geom[1]);
 308                 computeInside(insideLocs1, polys0, locPtPoly);
 309                 if (minDistance <= terminateDistance) {
 310 // flip locations, since we are testing geom 1 VS geom 0
 311                         assert( (*locPtPoly)[0] );
 312                         assert( (*locPtPoly)[1] );
 313                         (*minDistanceLocation)[0] = (*locPtPoly)[1];
 314                         (*minDistanceLocation)[1] = (*locPtPoly)[0];
 315                         delete locPtPoly;
 316                         for (size_t i=0; i<insideLocs1->size(); i++)
 317                         {
 318                                 GeometryLocation *l = (*insideLocs1)[i];
 319                                 if ( l != (*minDistanceLocation)[0] &&
 320                                         l != (*minDistanceLocation)[1] )
 321                                 {
 322                                         delete l;
 323                                 }
 324                         }
 325                         delete insideLocs1;
 326                         return;
 327                 }
 328                 for (size_t i=0; i<insideLocs1->size(); i++)
 329                         delete (*insideLocs1)[i];
 330                 delete insideLocs1;
 331         }
 332
 333         delete locPtPoly;
 334
 335         // If minDistance <= terminateDistance we must have
 336         // set minDistanceLocations to some non-null item
 337         assert( minDistance > terminateDistance ||
 338                 ( (*minDistanceLocation)[0] && (*minDistanceLocation)[1] ) );
 339 }
 340
 341
 342 /*private*/
 343 void
 344 DistanceOp::computeInside(vector<GeometryLocation*> *locs,
 345                 const Polygon::ConstVect& polys,
 346                 vector<GeometryLocation*> *locPtPoly)
 347 {
 348         for (size_t i=0, ni=locs->size(); i<ni; ++i)
 349         {
 350                 GeometryLocation *loc=(*locs)[i];
 351                 for (size_t j=0, nj=polys.size(); j<nj; ++j)
 352                 {
 353                         computeInside(loc, polys[j], locPtPoly);
 354                         if (minDistance<=terminateDistance) return;
 355                 }
 356         }
 357 }
 358
 359 /*private*/
 360 void
 361 DistanceOp::computeInside(GeometryLocation *ptLoc,
 362                 const Polygon *poly,
 363                 vector<GeometryLocation*> *locPtPoly)
 364 {
 365         const Coordinate &pt=ptLoc->getCoordinate();
 366
 367         // if pt is not in exterior, distance to geom is 0
 368         if (Location::EXTERIOR!=ptLocator.locate(pt, static_cast<const Geometry *>(poly)))
 369         {
 370                 minDistance = 0.0;
 371                 (*locPtPoly)[0] = ptLoc;
 372                 GeometryLocation *locPoly = new GeometryLocation(poly, pt);
 373                 (*locPtPoly)[1] = locPoly;
 374                 return;
 375         }
 376 }
 377
 378 /*private*/
 379 void
 380 DistanceOp::computeFacetDistance()
 381 {
 382         using geom::util::LinearComponentExtracter;
 383         using geom::util::PointExtracter;
 384
 385         vector<GeometryLocation*> locGeom(2);
 386
 387         /**
 388          * Geometries are not wholely inside, so compute distance from lines
 389          * and points
 390          * of one to lines and points of the other
 391          */
 392         LineString::ConstVect lines0;
 393         LineString::ConstVect lines1;
 394         LinearComponentExtracter::getLines(*(geom[0]), lines0);
 395         LinearComponentExtracter::getLines(*(geom[1]), lines1);
 396
 397 #if GEOS_DEBUG
 398         std::cerr << "LinearComponentExtracter found "
 399                   << lines0.size() << " lines in geometry 1 and "
 400                   << lines1.size() << " lines in geometry 2 "
 401                   << std::endl;
 402 #endif
 403
 404         Point::ConstVect pts0;
 405         Point::ConstVect pts1;
 406         PointExtracter::getPoints(*(geom[0]), pts0);
 407         PointExtracter::getPoints(*(geom[1]), pts1);
 408
 409 #if GEOS_DEBUG
 410         std::cerr << "PointExtracter found "
 411                   << pts0.size() << " points in geometry 1 and "
 412                   << pts1.size() << " points in geometry 2 "
 413                   << std::endl;
 414 #endif
 415
 416         // exit whenever minDistance goes LE than terminateDistance
 417         computeMinDistanceLines(lines0, lines1, locGeom);
 418         updateMinDistance(locGeom, false);
 419         if (minDistance <= terminateDistance) {
 420 #if GEOS_DEBUG
 421                 std::cerr << "Early termination after line-line distance" << std::endl;
 422 #endif
 423                 return;
 424         };
 425
 426         locGeom[0]=NULL;
 427         locGeom[1]=NULL;
 428         computeMinDistanceLinesPoints(lines0, pts1, locGeom);
 429         updateMinDistance(locGeom, false);
 430         if (minDistance <= terminateDistance) {
 431 #if GEOS_DEBUG
 432                 std::cerr << "Early termination after lines0-points1 distance" << std::endl;
 433 #endif
 434                 return;
 435         };
 436
 437         locGeom[0]=NULL;
 438         locGeom[1]=NULL;
 439         computeMinDistanceLinesPoints(lines1, pts0, locGeom);
 440         updateMinDistance(locGeom, true);
 441         if (minDistance <= terminateDistance){
 442 #if GEOS_DEBUG
 443                 std::cerr << "Early termination after lines1-points0 distance" << std::endl;
 444 #endif
 445                 return;
 446         };
 447
 448         locGeom[0]=NULL;
 449         locGeom[1]=NULL;
 450         computeMinDistancePoints(pts0, pts1, locGeom);
 451         updateMinDistance(locGeom, false);
 452
 453 #if GEOS_DEBUG
 454         std::cerr << "termination after pts-pts distance" << std::endl;
 455 #endif
 456 }
 457
 458 /*private*/
 459 void
 460 DistanceOp::computeMinDistanceLines(
 461                 const LineString::ConstVect& lines0,
 462                 const LineString::ConstVect& lines1,
 463                 vector<GeometryLocation*>& locGeom)
 464 {
 465         for (size_t i=0, ni=lines0.size(); i<ni; ++i)
 466         {
 467                 const LineString *line0=lines0[i];
 468                 for (size_t j=0, nj=lines1.size(); j<nj; ++j) {
 469                         const LineString *line1=lines1[j];
 470                         computeMinDistance(line0, line1, locGeom);
 471                         if (minDistance<=terminateDistance) return;
 472                 }
 473         }
 474 }
 475
 476 /*private*/
 477 void
 478 DistanceOp::computeMinDistancePoints(
 479                 const Point::ConstVect& points0,
 480                 const Point::ConstVect& points1,
 481                 vector<GeometryLocation*>& locGeom)
 482 {
 483         for (size_t i=0, ni=points0.size(); i<ni; ++i)
 484         {
 485                 const Point *pt0 = points0[i];
 486                 for (size_t j=0, nj=points1.size(); j<nj; ++j)
 487                 {
 488                         const Point *pt1 = points1[j];
 489                         double dist = pt0->getCoordinate()->distance(*(pt1->getCoordinate()));
 490
 491 #if GEOS_DEBUG
 492         std::cerr << "Distance "
 493                   << pt0->toString() << " - "
 494                   << pt1->toString() << ": "
 495                   << dist << ", minDistance: " << minDistance
 496                   << std::endl;
 497 #endif
 498
 499                         if (dist < minDistance)
 500                         {
 501                                 minDistance = dist;
 502                                 // this is wrong - need to determine closest points on both segments!!!
 503                                 delete locGeom[0];
 504                                 locGeom[0] = new GeometryLocation(pt0, 0, *(pt0->getCoordinate()));
 505                                 delete locGeom[1];
 506                                 locGeom[1] = new GeometryLocation(pt1, 0, *(pt1->getCoordinate()));
 507                         }
 508
 509                         if (minDistance<=terminateDistance) return;
 510                 }
 511         }
 512 }
 513
 514 /*private*/
 515 void
 516 DistanceOp::computeMinDistanceLinesPoints(
 517                 const LineString::ConstVect& lines,
 518                 const Point::ConstVect& points,
 519                 vector<GeometryLocation*>& locGeom)
 520 {
 521         for (size_t i=0;i<lines.size();i++)
 522         {
 523                 const LineString *line=lines[i];
 524                 for (size_t j=0;j<points.size();j++)
 525                 {
 526                         const Point *pt=points[j];
 527                         computeMinDistance(line,pt,locGeom);
 528                         if (minDistance<=terminateDistance) return;
 529                 }
 530         }
 531 }
 532
 533 /*private*/
 534 void
 535 DistanceOp::computeMinDistance(
 536                 const LineString *line0,
 537                 const LineString *line1,
 538                 vector<GeometryLocation*>& locGeom)
 539 {
 540         using geos::algorithm::CGAlgorithms;
 541
 542         const Envelope *env0=line0->getEnvelopeInternal();
 543         const Envelope *env1=line1->getEnvelopeInternal();
 544         if (env0->distance(env1)>minDistance) {
 545                 return;
 546         }
 547
 548         const CoordinateSequence *coord0=line0->getCoordinatesRO();
 549         const CoordinateSequence *coord1=line1->getCoordinatesRO();
 550         size_t npts0=coord0->getSize();
 551         size_t npts1=coord1->getSize();
 552
 553         // brute force approach!
 554         for(size_t i=0; i<npts0-1; ++i)
 555         {
 556                 for(size_t j=0; j<npts1-1; ++j)
 557                 {
 558                         double dist=CGAlgorithms::distanceLineLine(coord0->getAt(i),coord0->getAt(i+1),
 559                                 coord1->getAt(j),coord1->getAt(j+1));
 560                         if (dist < minDistance) {
 561                                 minDistance = dist;
 562                                 LineSegment seg0(coord0->getAt(i), coord0->getAt(i + 1));
 563                                 LineSegment seg1(coord1->getAt(j), coord1->getAt(j + 1));
 564                                 CoordinateSequence* closestPt = seg0.closestPoints(seg1);
 565                                 Coordinate *c1 = new Coordinate(closestPt->getAt(0));
 566                                 Coordinate *c2 = new Coordinate(closestPt->getAt(1));
 567                                 newCoords.push_back(c1);
 568                                 newCoords.push_back(c2);
 569                                 delete closestPt;
 570
 571                                 delete locGeom[0];
 572                                 locGeom[0] = new GeometryLocation(line0, i, *c1);
 573                                 delete locGeom[1];
 574                                 locGeom[1] = new GeometryLocation(line1, j, *c2);
 575                         }
 576                         if (minDistance<=terminateDistance) return;
 577                 }
 578         }
 579 }
 580
 581 /*private*/
 582 void
 583 DistanceOp::computeMinDistance(const LineString *line,
 584                 const Point *pt,
 585                 vector<GeometryLocation*>& locGeom)
 586 {
 587         using geos::algorithm::CGAlgorithms;
 588
 589         const Envelope *env0=line->getEnvelopeInternal();
 590         const Envelope *env1=pt->getEnvelopeInternal();
 591         if (env0->distance(env1)>minDistance) {
 592                 return;
 593         }
 594         const CoordinateSequence *coord0=line->getCoordinatesRO();
 595         Coordinate *coord=new Coordinate(*(pt->getCoordinate()));
 596         newCoords.push_back(coord);
 597
 598         // brute force approach!
 599         size_t npts0=coord0->getSize();
 600         for(size_t i=0; i<npts0-1; ++i)
 601         {
 602                 double dist=CGAlgorithms::distancePointLine(*coord,coord0->getAt(i),coord0->getAt(i+1));
 603                 if (dist < minDistance) {
 604                         minDistance = dist;
 605                         LineSegment seg(coord0->getAt(i), coord0->getAt(i + 1));
 606                         Coordinate segClosestPoint;
 607                         seg.closestPoint(*coord, segClosestPoint);
 608
 609                         delete locGeom[0];
 610                         locGeom[0] = new GeometryLocation(line, i, segClosestPoint);
 611                         delete locGeom[1];
 612                         locGeom[1] = new GeometryLocation(pt, 0, *coord);
 613                 }
 614                 if (minDistance<=terminateDistance) return;
 615         }
 616 }
 617
 618 /* public static */
 619 bool
 620 DistanceOp::isWithinDistance(const geom::Geometry& g0,
 621                              const geom::Geometry& g1,
 622                              double distance)
 623 {
 624         DistanceOp distOp(g0, g1, distance);
 625         return distOp.distance() <= distance;
 626 }
 627
 628 } // namespace geos.operation.distance
 629 } // namespace geos.operation
 630 } // namespace geos
 631
 632 /**********************************************************************
 633  * $Log$
 634  * Revision 1.24  2006/06/12 11:29:24  strk
 635  * unsigned int => size_t
 636  *
 637  * Revision 1.23  2006/03/24 09:52:41  strk
 638  * USE_INLINE => GEOS_INLINE
 639  *
 640  * Revision 1.22  2006/03/23 12:12:01  strk
 641  * Fixes to allow build with -DGEOS_INLINE
 642  *
 643  * Revision 1.21  2006/03/21 17:55:01  strk
 644  * opDistance.h header split
 645  *
 646  * Revision 1.20  2006/03/03 10:46:21  strk
 647  * Removed 'using namespace' from headers, added missing headers in .cpp files, removed useless includes in headers (bug#46)
 648  **********************************************************************/
 649