[geos-devel] Geometry/Rectangle Intersection: line touching rectangle

[Mika Heiskanen]

On 09/12/2014 03:23 AM, Martin Davis wrote:
>
> I agree that having a built-in way to extract only the outer polygons
> from a polygonized dataset would be useful.  Someone else asked for this
> just yesterday. I haven't come across this use case before, so hadn't
> realized it was needed.  Have you made any attempt to share the code you
> wrote to solve this problem?  You can do so on the JTS list, or via the
> usual ways of sharing code (ie GitHub).  If not, please don't complain
> that it isn't present in JTS.

Actually I initially reported it as a bug:

   http://lists.osgeo.org/pipermail/geos-devel/2014-March/006821.html

and was explained it was not. I was then informed of PostGIS BuildArea
function. The function does a final overlay operation, which I thought 
we would not need. I then offered our a bit different solution, but
there was no interest then:

   http://lists.osgeo.org/pipermail/geos-devel/2014-March/006831.html

I'm including the code below so you can check if you see any
use for the code. The main code has called normalize for all the input
polygons, and then calls buildGeometry to build the final result.


    * \brief Remove superfluous polygons from the list
    *
    * We need to remove those polygons whose outer rings
    * are holes in another polygon. For the original algorithm
    * see PostGIS liblwgeom/lwgeom_geos.c source code. However,
    * we can improve on it since we know for a fact that all
    * edges in the input appeared only once. Hence we can just
    * take a sample edge from all rings, and see if they appear
    * both in an outer ring and in a hole. If that happens
    * we remove the polygon where the edge is in the outer
    * ring. However, there may be multiple nesting, we must
    * only remove the polygon if its number of parents is
    * odd.
    *
    * For this to be efficient each polygon should be in normalized
    * form, which is actually good also for the stability of the
    * regression tests. In the normalized form the inner and
    * outer rings will have the same first coordinate, but
    * the order of the remaining coordinates is reversed for
    * holes. Hence we can check whether the first and second coordinates
    * match the first and second to last coordinates of the hole.
    *
    * PostGIS compares the full ring and also does an
    * extra overlay on the final result to dissolve duplicate
    * edges.
    */



   void 
GeosBuilder::remove_superfluous_polygons(std::vector<geos::geom::Polygon 
*>  & polys)
   {
     namespace gg = geos::geom;

     // No need for the extra check if there are no multiple polygons
     // or if there are no holes

     if(polys.size() < 2 || !has_holes(polys))
       return;

     // Child polygon indices
     std::vector<int> children(polys.size(), -1);


     for(std::size_t p=0; p<polys.size(); p++)
       {
         for(std::size_t h=0; h<polys[p]->getNumInteriorRing(); h++)
           {
             const gg::LineString * hole = polys[p]->getInteriorRingN(h);
             if(hole != NULL)
               {
                 const gg::CoordinateSequence * cs = 
hole->getCoordinatesRO();
                 const gg::Coordinate & p1 = cs->getAt(0);
                 const gg::Coordinate & p2 = cs->getAt(cs->size()-2);

                 // Now find if there is a matching outer ring

                 for(std::size_t pp=0; pp<polys.size(); pp++)
                   {
                     if(p != pp && polys[pp] != NULL)
                       {
                         const gg::LineString * shell = 
polys[pp]->getExteriorRing();
                         if(shell != NULL)
                           {
                             const gg::CoordinateSequence * cs2 = 
shell->getCoordinatesRO();
                             const gg::Coordinate & sp1 = cs2->getAt(0);
                             const gg::Coordinate & sp2 = cs2->getAt(1);

                             if(p1 == sp1 && p2 == sp2)
                               {
                                 children[pp] = static_cast<int>(p);
                                 break;
                               }
                           }
                       }
                       }
               }
           }
       }

     // Now delete polygons which have an even number of children

     std::size_t outpos = 0;
     for(std::size_t p=0; p<polys.size(); p++)
       {
         int count = 0;
         int pos = static_cast<int>(p);
         while(children[pos] >= 0)
           {
             ++count;
             pos = children[pos];
           }
         if(count % 2 == 1)
           {
             delete polys[p];
           }
         else
           {
             polys[outpos++] = polys[p];
           }
       }
     polys.resize(outpos);

   }

   bool GeosBuilder::has_holes(const std::vector<geos::geom::Polygon *> 
  & polys)
   {
     if(polys.empty())
       return false;

     for(std::size_t i=0; i<polys.size(); i++)
       {
         if( polys[i]->getNumInteriorRing() != 0)
           return true;
       }

     return false;
   }