Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  S-57 Reader
       4             :  * Purpose:  Implements polygon assembly from a bunch of arcs.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 1999, Frank Warmerdam
       9             :  * Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_port.h"
      15             : #include "ogr_api.h"
      16             : 
      17             : #include <cmath>
      18             : #include <cstddef>
      19             : #include <list>
      20             : #include <vector>
      21             : 
      22             : #include "ogr_core.h"
      23             : #include "ogr_geometry.h"
      24             : #include "cpl_conv.h"
      25             : #include "cpl_error.h"
      26             : 
      27             : /************************************************************************/
      28             : /*                            CheckPoints()                             */
      29             : /*                                                                      */
      30             : /*      Check if two points are closer than the current best            */
      31             : /*      distance.  Update the current best distance if they are.        */
      32             : /************************************************************************/
      33             : 
      34       28692 : static bool CheckPoints(OGRLineString *poLine1, int iPoint1,
      35             :                         OGRLineString *poLine2, int iPoint2,
      36             :                         double *pdfDistance)
      37             : 
      38             : {
      39       28692 :     if (pdfDistance == nullptr || *pdfDistance == 0)
      40             :     {
      41       35781 :         if (poLine1->getX(iPoint1) == poLine2->getX(iPoint2) &&
      42        7291 :             poLine1->getY(iPoint1) == poLine2->getY(iPoint2))
      43             :         {
      44        5332 :             if (pdfDistance)
      45        4954 :                 *pdfDistance = 0.0;
      46        5332 :             return true;
      47             :         }
      48       23158 :         return false;
      49             :     }
      50             : 
      51             :     const double dfDeltaX =
      52         202 :         std::abs(poLine1->getX(iPoint1) - poLine2->getX(iPoint2));
      53             : 
      54         202 :     if (dfDeltaX > *pdfDistance)
      55          30 :         return false;
      56             : 
      57             :     const double dfDeltaY =
      58         172 :         std::abs(poLine1->getY(iPoint1) - poLine2->getY(iPoint2));
      59             : 
      60         172 :     if (dfDeltaY > *pdfDistance)
      61           6 :         return false;
      62             : 
      63         166 :     const double dfDistance = sqrt(dfDeltaX * dfDeltaX + dfDeltaY * dfDeltaY);
      64             : 
      65         166 :     if (dfDistance < *pdfDistance)
      66             :     {
      67         166 :         *pdfDistance = dfDistance;
      68         166 :         return true;
      69             :     }
      70             : 
      71           0 :     return false;
      72             : }
      73             : 
      74             : /************************************************************************/
      75             : /*                           AddEdgeToRing()                            */
      76             : /************************************************************************/
      77             : 
      78        2717 : static void AddEdgeToRing(OGRLinearRing *poRing, OGRLineString *poLine,
      79             :                           bool bReverse, double dfTolerance)
      80             : 
      81             : {
      82             :     /* -------------------------------------------------------------------- */
      83             :     /*      Establish order and range of traverse.                          */
      84             :     /* -------------------------------------------------------------------- */
      85        2717 :     const int nVertToAdd = poLine->getNumPoints();
      86             : 
      87        2717 :     int iStart = bReverse ? nVertToAdd - 1 : 0;
      88        2717 :     const int iEnd = bReverse ? 0 : nVertToAdd - 1;
      89        2717 :     const int iStep = bReverse ? -1 : 1;
      90             : 
      91             :     /* -------------------------------------------------------------------- */
      92             :     /*      Should we skip a repeating vertex?                              */
      93             :     /* -------------------------------------------------------------------- */
      94        5083 :     if (poRing->getNumPoints() > 0 &&
      95        2366 :         CheckPoints(poRing, poRing->getNumPoints() - 1, poLine, iStart,
      96             :                     &dfTolerance))
      97             :     {
      98        2366 :         iStart += iStep;
      99             :     }
     100             : 
     101        2717 :     poRing->addSubLineString(poLine, iStart, iEnd);
     102        2717 : }
     103             : 
     104             : /************************************************************************/
     105             : /*                      OGRBuildPolygonFromEdges()                      */
     106             : /************************************************************************/
     107             : 
     108             : /**
     109             :  * Build a ring from a bunch of arcs.
     110             :  *
     111             :  * @param hLines handle to an OGRGeometryCollection (or OGRMultiLineString)
     112             :  * containing the line string geometries to be built into rings.
     113             :  * @param bBestEffort not yet implemented???.
     114             :  * @param bAutoClose indicates if the ring should be close when first and
     115             :  * last points of the ring are the same.
     116             :  * @param dfTolerance tolerance into which two arcs are considered
     117             :  * close enough to be joined.
     118             :  * @param peErr OGRERR_NONE on success, or OGRERR_FAILURE on failure.
     119             :  * @return a handle to the new geometry, a polygon.
     120             :  *
     121             :  */
     122             : 
     123         343 : OGRGeometryH OGRBuildPolygonFromEdges(OGRGeometryH hLines,
     124             :                                       CPL_UNUSED int bBestEffort,
     125             :                                       int bAutoClose, double dfTolerance,
     126             :                                       OGRErr *peErr)
     127             : 
     128             : {
     129         343 :     if (hLines == nullptr)
     130             :     {
     131           0 :         if (peErr != nullptr)
     132           0 :             *peErr = OGRERR_NONE;
     133           0 :         return nullptr;
     134             :     }
     135             : 
     136             :     /* -------------------------------------------------------------------- */
     137             :     /*      Check for the case of a geometrycollection that can be          */
     138             :     /*      promoted to MultiLineString.                                    */
     139             :     /* -------------------------------------------------------------------- */
     140         343 :     OGRGeometry *poGeom = OGRGeometry::FromHandle(hLines);
     141         343 :     if (wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection)
     142             :     {
     143        3055 :         for (auto &&poMember : poGeom->toGeometryCollection())
     144             :         {
     145        2715 :             if (wkbFlatten(poMember->getGeometryType()) != wkbLineString)
     146             :             {
     147           1 :                 if (peErr != nullptr)
     148           1 :                     *peErr = OGRERR_FAILURE;
     149           1 :                 CPLError(CE_Failure, CPLE_NotSupported,
     150             :                          "The geometry collection contains "
     151             :                          "non-line string geometries");
     152           1 :                 return nullptr;
     153             :             }
     154             :         }
     155             :     }
     156           2 :     else if (wkbFlatten(poGeom->getGeometryType()) != wkbMultiLineString)
     157             :     {
     158           1 :         if (peErr != nullptr)
     159           1 :             *peErr = OGRERR_FAILURE;
     160           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     161             :                  "The passed geometry is not an OGRGeometryCollection "
     162             :                  "(or OGRMultiLineString) "
     163             :                  "containing line string geometries");
     164           1 :         return nullptr;
     165             :     }
     166             : 
     167         341 :     bool bSuccess = true;
     168         341 :     OGRGeometryCollection *poLines = poGeom->toGeometryCollection();
     169         682 :     std::vector<OGRLinearRing *> apoRings;
     170             : 
     171             :     /* -------------------------------------------------------------------- */
     172             :     /*      Setup array of line markers indicating if they have been        */
     173             :     /*      added to a ring yet.                                            */
     174             :     /* -------------------------------------------------------------------- */
     175         341 :     const int nEdges = poLines->getNumGeometries();
     176         682 :     std::list<OGRLineString *> oListEdges;
     177        3060 :     for (int i = 0; i < nEdges; i++)
     178             :     {
     179        2719 :         OGRLineString *poLine = poLines->getGeometryRef(i)->toLineString();
     180        2719 :         if (poLine->getNumPoints() >= 2)
     181             :         {
     182        2717 :             oListEdges.push_back(poLine);
     183             :         }
     184             :     }
     185             : 
     186             :     /* ==================================================================== */
     187             :     /*      Loop generating rings.                                          */
     188             :     /* ==================================================================== */
     189         692 :     while (!oListEdges.empty())
     190             :     {
     191             :         /* --------------------------------------------------------------------
     192             :          */
     193             :         /*      Find the first unconsumed edge. */
     194             :         /* --------------------------------------------------------------------
     195             :          */
     196         351 :         OGRLineString *poLine = oListEdges.front();
     197         351 :         oListEdges.erase(oListEdges.begin());
     198             : 
     199             :         /* --------------------------------------------------------------------
     200             :          */
     201             :         /*      Start a new ring, copying in the current line directly */
     202             :         /* --------------------------------------------------------------------
     203             :          */
     204         351 :         OGRLinearRing *poRing = new OGRLinearRing();
     205             : 
     206         351 :         AddEdgeToRing(poRing, poLine, FALSE, 0);
     207             : 
     208             :         /* ====================================================================
     209             :          */
     210             :         /*      Loop adding edges to this ring until we make a whole pass */
     211             :         /*      within finding anything to add. */
     212             :         /* ====================================================================
     213             :          */
     214         351 :         bool bWorkDone = true;
     215         351 :         double dfBestDist = dfTolerance;
     216             : 
     217        2717 :         while (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1,
     218        2371 :                             nullptr) &&
     219        2717 :                !oListEdges.empty() && bWorkDone)
     220             :         {
     221        2366 :             bool bReverse = false;
     222             : 
     223        2366 :             bWorkDone = false;
     224        2366 :             dfBestDist = dfTolerance;
     225             : 
     226             :             // We consider linking the end to the beginning.  If this is
     227             :             // closer than any other option we will just close the loop.
     228             : 
     229             :             // CheckPoints(poRing, 0, poRing, poRing->getNumPoints()-1,
     230             :             //             &dfBestDist);
     231             : 
     232             :             // Find unused edge with end point closest to our loose end.
     233        2366 :             OGRLineString *poBestEdge = nullptr;
     234        2366 :             std::list<OGRLineString *>::iterator oBestIter;
     235       11604 :             for (auto oIter = oListEdges.begin(); oIter != oListEdges.end();
     236        9238 :                  ++oIter)
     237             :             {
     238       11592 :                 poLine = *oIter;
     239             : 
     240       11592 :                 if (CheckPoints(poLine, 0, poRing, poRing->getNumPoints() - 1,
     241             :                                 &dfBestDist))
     242             :                 {
     243        2031 :                     poBestEdge = poLine;
     244        2031 :                     oBestIter = oIter;
     245        2031 :                     bReverse = false;
     246             :                 }
     247       11592 :                 if (CheckPoints(poLine, poLine->getNumPoints() - 1, poRing,
     248       11592 :                                 poRing->getNumPoints() - 1, &dfBestDist))
     249             :                 {
     250         335 :                     poBestEdge = poLine;
     251         335 :                     oBestIter = oIter;
     252         335 :                     bReverse = true;
     253             :                 }
     254             : 
     255             :                 // If we found an exact match, jump now.
     256       11592 :                 if (dfBestDist == 0.0 && poBestEdge != nullptr)
     257        2354 :                     break;
     258             :             }
     259             : 
     260             :             // We found one within tolerance - add it.
     261        2366 :             if (poBestEdge)
     262             :             {
     263        2366 :                 AddEdgeToRing(poRing, poBestEdge, bReverse, dfTolerance);
     264             : 
     265        2366 :                 oListEdges.erase(oBestIter);
     266        2366 :                 bWorkDone = true;
     267             :             }
     268             :         }
     269             : 
     270             :         /* --------------------------------------------------------------------
     271             :          */
     272             :         /*      Did we fail to complete the ring? */
     273             :         /* --------------------------------------------------------------------
     274             :          */
     275         351 :         dfBestDist = dfTolerance;
     276             : 
     277         351 :         if (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1,
     278             :                          &dfBestDist))
     279             :         {
     280           0 :             CPLDebug("OGR",
     281             :                      "Failed to close ring %d.\n"
     282             :                      "End Points are: (%.8f,%.7f) and (%.7f,%.7f)",
     283           0 :                      static_cast<int>(apoRings.size()), poRing->getX(0),
     284           0 :                      poRing->getY(0), poRing->getX(poRing->getNumPoints() - 1),
     285           0 :                      poRing->getY(poRing->getNumPoints() - 1));
     286             : 
     287           0 :             bSuccess = false;
     288             :         }
     289             : 
     290             :         /* --------------------------------------------------------------------
     291             :          */
     292             :         /*      Do we need to auto-close this ring? */
     293             :         /* --------------------------------------------------------------------
     294             :          */
     295         351 :         dfBestDist = dfTolerance;
     296             : 
     297         351 :         if (bAutoClose)
     298             :         {
     299          37 :             if (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1,
     300             :                              &dfBestDist))
     301             :             {
     302           0 :                 poRing->addPoint(poRing->getX(0), poRing->getY(0),
     303             :                                  poRing->getZ(0));
     304             :             }
     305          37 :             else if (!CheckPoints(poRing, 0, poRing, poRing->getNumPoints() - 1,
     306             :                                   nullptr))
     307             :             {
     308             :                 // The endpoints are very close but do not exactly match.
     309             :                 // Alter the last one so it is equal to the first, to prevent
     310             :                 // invalid self-intersecting rings.
     311           5 :                 poRing->setPoint(poRing->getNumPoints() - 1, poRing->getX(0),
     312             :                                  poRing->getY(0), poRing->getZ(0));
     313             :             }
     314             :         }
     315             : 
     316         351 :         apoRings.push_back(poRing);
     317             :     }  // Next ring.
     318             : 
     319             :     /* -------------------------------------------------------------------- */
     320             :     /*      Identify exterior ring - it will be the largest.  #3610         */
     321             :     /* -------------------------------------------------------------------- */
     322         341 :     double maxarea = 0.0;
     323         341 :     int maxring = -1;
     324         341 :     OGREnvelope tenv;
     325             : 
     326         692 :     for (int rn = 0; rn < static_cast<int>(apoRings.size()); ++rn)
     327             :     {
     328         351 :         apoRings[rn]->getEnvelope(&tenv);
     329         351 :         const double tarea = (tenv.MaxX - tenv.MinX) * (tenv.MaxY - tenv.MinY);
     330         351 :         if (tarea > maxarea)
     331             :         {
     332         342 :             maxarea = tarea;
     333         342 :             maxring = rn;
     334             :         }
     335             :     }
     336             : 
     337         341 :     OGRPolygon *poPolygon = new OGRPolygon();
     338             : 
     339         341 :     if (maxring != -1)
     340             :     {
     341         341 :         poPolygon->addRingDirectly(apoRings[maxring]);
     342         692 :         for (int rn = 0; rn < static_cast<int>(apoRings.size()); ++rn)
     343             :         {
     344         351 :             if (rn == maxring)
     345         341 :                 continue;
     346          10 :             poPolygon->addRingDirectly(apoRings[rn]);
     347             :         }
     348             :     }
     349             :     else
     350             :     {
     351           0 :         for (auto &poRing : apoRings)
     352           0 :             delete poRing;
     353             :     }
     354             : 
     355         341 :     if (peErr != nullptr)
     356             :     {
     357         341 :         *peErr = bSuccess ? OGRERR_NONE : OGRERR_FAILURE;
     358             :     }
     359             : 
     360         341 :     return OGRGeometry::ToHandle(poPolygon);
     361             : }