src/foam/interpolations/GGIInterpolation/GGIInterpolationPolygonIntersection.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 Description
  25     2D Polygon intersection algorithms
  26
  27 Author
  28     Martin Beaudoin, Hydro-Quebec, (2008)
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #include "boundBox.H"
  33 #include "plane.H"
  34
  35 // Point in polygon algorithm
  36 #include "HormannAgathos.H"
  37
  38 // Polygon clipping algorithm
  39 #include "SutherlandHodgman.H"
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45
  46 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  47
  48
  49 // Computation of the intersection area between 2 2D polygons.
  50 //
  51 // Two following algorithms are considered here:
  52 //   Sutherland-Hodgman algorithm : a classic, but polygons need to be convex
  53 //   Greiner-Hormann algorithm    : very nice algorithm; accept both convex
  54 //                                  and concave polygons
  55 //
  56 //   We also rely here on a very efficient algorithm for the point in
  57 //   polygon problem for arbitrary polygons. That algorithm was proposed
  58 //   by Hormann and Agathos.
  59
  60 template<class MasterPatch, class SlavePatch>
  61 scalar
  62 GGIInterpolation<MasterPatch, SlavePatch>::polygonIntersection
  63 (
  64     const List<point2D>& poly1,
  65     const List<point2D>& poly2
  66 ) const
  67 {
  68     // Using pointers because clipping and subject may be swapped
  69     // by the algorithm.  HJ, 24/Oct/2008
  70
  71      // The first polygon will be the clipping polygon
  72     const List<point2D>* clippingPolygon = &poly1;
  73
  74     // Yhe second polygon will be the subject polygon; the one that
  75     // is going to be clipped
  76     const List<point2D>* subjectPolygon = &poly2;
  77
  78     // Empty list so we can detect weird cases later
  79     List<point2D> clippedPolygon;
  80     scalar intersectionArea = 0.0;
  81
  82     // First, let's get rid of the obvious:
  83     //  1: Neither polygons intersect one another
  84     //     --> intersection area == 0.0 : that should not happen
  85     //  2: If both polygons completely overlap one another
  86     //     ---> subjectPolygon is the intersection polygon
  87     //  3: If clippingPolygon totally enclosed subjectPolygon
  88     //     ---> subjectPolygon is the intersecting polygon
  89     //  4: If subjectPolygon totally enclosed clippingPolygon
  90     //     --> clippingPolygon is the intersecting polygon
  91     //  5: Otherwise, we have partial or full intersection
  92     //
  93     //  For this, we first detect if the vertices of the subject polygon
  94     //  are inside or outside of the clipping polygon
  95     //  This is a quick intersection test...
  96
  97     // Keep track of who is inside or outside
  98     List<bool> subjectVertexInside(subjectPolygon->size());
  99
 100     insideOutside statusInOut =
 101         isVertexInsidePolygon
 102         (
 103             *clippingPolygon,
 104             *subjectPolygon,
 105             subjectVertexInside
 106         );
 107
 108     // We check right away if statusInOut == ALL_OUTSIDE
 109     // Sometimes, it is just that the clipping polygon is inside the
 110     // subject polygon instead So let's explore this situation, just
 111     // in case
 112     if (statusInOut == ALL_OUTSIDE)
 113     {
 114         // We need to check if subject is not completely or partially
 115         // enclosing clipping instead
 116
 117         clippingPolygon = &poly2;
 118         subjectPolygon  = &poly1;
 119
 120         subjectVertexInside.setSize(subjectPolygon->size());
 121         statusInOut =
 122             isVertexInsidePolygon
 123             (
 124                 *clippingPolygon,
 125                 *subjectPolygon,
 126                 subjectVertexInside
 127             );
 128
 129     }
 130
 131     switch(statusInOut)
 132     {
 133         case ALL_INSIDE:
 134         {
 135             clippedPolygon = *subjectPolygon;
 136             break;
 137         }
 138         case ALL_OUTSIDE:
 139         case PARTIALLY_OVERLAPPING:
 140         default:
 141         {
 142             // Compute the intersection
 143
 144             // If by any chance, we have reached a situation where the
 145             // intersection is really zero, it is because the quick
 146             // reject tests have missed something. The intersection
 147             // area will be 0.0, and the calling function should at
 148             // least check for this, and report loudly...
 149             clippedPolygon =
 150                 clipPolygon2DSutherlandHodgman
 151                 (
 152                     *clippingPolygon,
 153                     *subjectPolygon
 154                 );// For convex polygons
 155
 156 #if 0
 157             // This is the next candidate to code if the Sutherland
 158             // Hodgman algorithm encounters concave polygons...
 159             clippedPolygon =
 160                 clipPolygon2DGreinerHormann
 161                 (
 162                     *clippingPolygon,
 163                     *subjectPolygon,
 164                     subjectVertexInside
 165                 );   // For arbitrary polygons
 166 #endif
 167             break;
 168         }
 169     }
 170
 171     // Compute the area of clippedPolygon if we indeed do have a
 172     // clipped polygon; otherwise, the computed area stays at 0.0;
 173     if (clippedPolygon.size() > 2)
 174     {
 175         // We are only interested in the absolute value
 176         intersectionArea = mag(area2D(clippedPolygon));
 177     }
 178
 179     if (debug)
 180     {
 181         // Check against tolerances
 182         scalar clippingArea  = area2D(*clippingPolygon);
 183         scalar subjectArea   = area2D(*subjectPolygon);
 184
 185         if
 186         (
 187             mag(intersectionArea/clippingArea) < areaErrorTol_()
 188          || mag(intersectionArea/subjectArea)  < areaErrorTol_()
 189         )
 190         {
 191             WarningIn
 192             (
 193                 "GGIInterpolation<MasterPatch, SlavePatch>::"
 194                 "polygonIntersection"
 195             )   << "Intersection might be wrong.  Clipping side "
 196                 << intersectionArea/clippingArea << " subject: "
 197                 << intersectionArea/subjectArea << endl;
 198         }
 199     }
 200
 201     return intersectionArea;
 202 }
 203
 204
 205 // Compute the point in polygon problem using the computation of the
 206 // winding number. This technique is based on a paper by Hormann and
 207 // Agathos.
 208 //
 209 // This is based on the winding number technique, but optimized in
 210 // order to only evaluate quarter revolutions instead of the whole
 211 // arcos/sqrt basic algorithm.  When the GGI weighting factors will
 212 // have to be recomputed often for moving meshes, this performance
 213 // will be useful.  We can also go back to the classical winding
 214 // number algorithm if need be Notice: The list subjectVertexInside
 215 // will return a boolean marking if a point from the subject polygon
 216 // is inside or outside the clipping polygon. A point is "inside"
 217 // (value == true) if it is inside the clipping polygon, or on a
 218 // vertex or edge.  This information will be very useful for the
 219 // Sutherland Hodgman algo.
 220 template<class MasterPatch, class SlavePatch>
 221 typename Foam::GGIInterpolation<MasterPatch, SlavePatch>::insideOutside
 222 GGIInterpolation<MasterPatch, SlavePatch>::isVertexInsidePolygon
 223 (
 224     const List<point2D>& clippingPolygon,
 225     const List<point2D>& subjectPolygon,
 226     List<bool>& subjectVertexInside
 227 ) const
 228 {
 229     insideOutside retValue = ALL_OUTSIDE;
 230
 231     // The class HormannAgathos implements the algorithm
 232     // We use distErrorTol_ to evaluate a distance factor called
 233     // epsilon.  That epsilon factor will be used to detect if a point
 234     // is on a vertex or an edge.
 235     const scalar distErrorTol = sqrt(areaErrorTol_());
 236
 237     HormannAgathos pip(clippingPolygon, distErrorTol);
 238
 239     // Counter
 240     label nbrsOutside = 0;
 241
 242     // Iterate over all subject vertices, determining if they are
 243     // inside/outside the clipping polygon
 244     forAll(subjectPolygon, sPI)
 245     {
 246         switch (pip.evaluate(subjectPolygon[sPI]))
 247         {
 248             case HormannAgathos::POINT_OUTSIDE:
 249                 nbrsOutside++;
 250                 subjectVertexInside[sPI] = false;
 251                 break;
 252             case HormannAgathos::POINT_ON_VERTEX:
 253             case HormannAgathos::POINT_ON_EDGE:
 254             case HormannAgathos::POINT_INSIDE:
 255             default:
 256                 // Vertex, Edge or Inside: We are inside!
 257                 subjectVertexInside[sPI] = true;
 258                 break;
 259         }
 260     }
 261
 262     // Let's do the inventory now...
 263     if (nbrsOutside == 0)
 264     {
 265         retValue = ALL_INSIDE;
 266     }
 267     else if (nbrsOutside < subjectPolygon.size())
 268     {
 269         retValue = PARTIALLY_OVERLAPPING;
 270     }
 271     // else, all the points are outside, which is not necessarily a
 272     // problem if the subject Polygon is enclosing partially or
 273     // completely the clipping polygon instead
 274
 275     return retValue;
 276 }
 277
 278
 279 // This is an implementation of the Sutherland Hodgman algorithm:
 280 // Reentrant Polygon Clipping, Sutherland, Hodgman, Communications of
 281 // the ACM, 1974
 282 //
 283 // Wikipedia has a very simple Pseudo-Code example of this rather
 284 // simple algorithm http://en.wikipedia.org/wiki/Sutherland-Hodgeman.
 285 //
 286 // The subject polygon will be clipped by the clipping polygon. The
 287 // list of boolean values subjectVertexInside provide the bonus
 288 // information if a subject vertex is either "inside" or "outside"
 289 // the clipping polygon. A vertex is "inside" if it is inside, on a
 290 // vertex, or on an edge. Otherwise, the vertex is "outside"
 291 template<class MasterPatch, class SlavePatch>
 292 List<point2D>
 293 GGIInterpolation<MasterPatch, SlavePatch>::clipPolygon2DSutherlandHodgman
 294 (
 295     const List<point2D>& clippingPolygon,
 296     const List<point2D>& subjectPolygon
 297 ) const
 298 {
 299     // Create Sutherland-Hodgman intersector, with a distance tolerance
 300     // factor for intersection computation
 301     // We use a tolerance that is consistant with  the quick reject tests
 302     return SutherlandHodgman
 303     (
 304         clippingPolygon,
 305         subjectPolygon,
 306         sqrt(areaErrorTol_()) // = distErrorTol
 307     ).evaluate();
 308 }
 309
 310
 311 // Compute the area of a polygon
 312 // See:  http://mathworld.wolfram.com/PolygonArea.html
 313 template<class MasterPatch, class SlavePatch>
 314 scalar GGIInterpolation<MasterPatch, SlavePatch>::area2D
 315 (
 316     const List<point2D>& polygon
 317 ) const
 318 {
 319     // For a non-self-intersecting (simple) polygon with n vertices,
 320     // the area is :
 321     // A = 0.5 * sum(x(i)y(i+1) - x(i+1)y(i)  , i=1 to n; where n+1 = 0;
 322     scalar area = 0;
 323
 324     // We start with last term from Wolfram equation
 325     label indexI   = polygon.size() - 1;
 326     label indexIp1 = indexI;
 327
 328     for (label i = 0; i < polygon.size(); i++)
 329     {
 330         indexI = indexIp1;
 331         indexIp1 = i;
 332
 333         area += polygon[indexI].x()*polygon[indexIp1].y()
 334               - polygon[indexIp1].x()*polygon[indexI].y();
 335     }
 336
 337     // NB: THe area can be positive or negative
 338     return area/2.0;
 339 }
 340
 341
 342 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 343
 344 } // End namespace Foam
 345
 346 // ************************************************************************* //