src/foam/algorithms/polygon/clipping/SutherlandHodgman.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 Author
  25     Martin Beaudoin, Hydro-Quebec, (2008)
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "SutherlandHodgman.H"
  30
  31 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  32
  33 // Compute the intersection between two line segments
  34 // See:
  35 //     http://mathworld.wolfram.com/Line-LineIntersection.html
  36 //     http://local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
  37 //
  38 //  We are only interested in the actual interpolation where ua and ub are
  39 //  both in the range [0, 1]. If outside of this range, we simply return false.
  40 bool Foam::SutherlandHodgman::lineSegmentIntersection
  41 (
  42     const point2D& p1,  // Segment 1 [p1, p2]
  43     const point2D& p2,
  44     const point2D& p3,  // Segment 2 [p3, p4]
  45     const point2D& p4,
  46     point2D& intersectionPoint
  47 ) const
  48 {
  49     bool segmentDoIntersect = false;
  50
  51     intersectionPoint[0] = 0.0;
  52     intersectionPoint[1] = 0.0;
  53
  54     vector2D delta21(p2 - p1);
  55     vector2D delta43(p4 - p3);
  56     vector2D delta13(p1 - p3);
  57
  58     scalar denom = delta43[1]*delta21[0] - delta43[0]*delta21[1];
  59
  60     // Comparison with 0 and 1...  Floating point precision
  61     // issue. Cannot compare with 0 and 1 exactly.  VSMALL is too
  62     // strict... we need to take into account that we tolerate an
  63     // epsilon error all along so in the worst case, if a point P is
  64     // separated by only epsilon from an edge, the smallest
  65     // determinant should be of the order of 2*(2*epsilon*2epsilon) =
  66     // 8*epsilon^2
  67     //
  68     // Epsilon is adjusted to the smallest edge length between [p1,p2]
  69     // and [p3,p4]
  70
  71     // Accounting for double precision tolerance.  HJ, 16/Jan/2009
  72     scalar epsilon =
  73         Foam::max
  74         (
  75             Foam::min(mag(delta21), mag(delta43))*intersectSegDistTol_,
  76             SMALL
  77         );
  78
  79     scalar tolFactor = 8.0*sqr(epsilon);
  80
  81     if (mag(denom) > tolFactor)
  82     {
  83         scalar ua = (delta43[0]*delta13[1] - delta43[1]*delta13[0])/denom;
  84         scalar ub = (delta21[0]*delta13[1] - delta21[1]*delta13[0])/denom;
  85
  86         //- Adjust tolFactor because of division:
  87         //  (y +- epsilon)/(x +- epsilon) = y/x +- 2*epsilon (rough estimation)
  88         tolFactor *= 2.0;
  89
  90         // Check for end points outside of range 0..1
  91         //if (ua >= 0.0 && ua <= 1.0 || ub >= 0.0 && ub <= 1.0)
  92         // One of the intersection must be inside of one of the segment
  93         // We don't want both of them at the same time, because each edge
  94         // should be considered as infinite line here. But still, the
  95         // intersection must be within the interval [0,1] for one of
  96         // the edge
  97         if (ua >= -tolFactor && ua <= (1.0 + tolFactor))
  98         {
  99             // This is the intersection we want
 100             intersectionPoint = p1 + ua*delta21;
 101             segmentDoIntersect = true;
 102
 103             if (mag(intersectionPoint - (p3 + ub*delta43)) > epsilon)
 104             {
 105                 WarningIn
 106                 (
 107                     "Foam::SutherlandHodgman::lineSegmentIntersection()"
 108                 )   << "ua does not match with ub: delta: "
 109                     << mag(intersectionPoint - (p3 + ub*delta43))
 110                     << " : epsilon: " << epsilon
 111                     << endl;
 112             }
 113         }
 114         else if (ub >= -tolFactor && ub <= (1.0 + tolFactor))
 115         {
 116             // This is the intersection we want
 117             intersectionPoint = p3 + ub*delta43;
 118             segmentDoIntersect = true;
 119
 120             if (mag(intersectionPoint - (p1 + ua*delta21)) > epsilon)
 121             {
 122                 WarningIn
 123                 (
 124                     "Foam::SutherlandHodgman::lineSegmentIntersection()"
 125                 )   << "ub does not match with ua: delta: "
 126                     << mag(intersectionPoint - (p1 + ua*delta21))
 127                     << " : epsilon: " << epsilon
 128                     << endl;
 129             }
 130         }
 131     }
 132
 133     // Parallel lines. We don't go any further
 134     return segmentDoIntersect;
 135 }
 136
 137
 138 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 139
 140 Foam::SutherlandHodgman::SutherlandHodgman
 141 (
 142     const List<point2D>& clippingPolygon,
 143     const List<point2D>& subjectPolygon,
 144     const scalar& intersectSegDistTol
 145 )
 146 :
 147     subjectPolygon_(subjectPolygon),
 148     clippingPolygon_(clippingPolygon),
 149     currentClipEdgeP1_(clippingPolygon.size() - 1),
 150     currentClipEdgeP2_(0),
 151     intersectSegDistTol_(intersectSegDistTol)
 152 {}
 153
 154
 155
 156 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 157
 158 Foam::List<Foam::point2D> Foam::SutherlandHodgman::evaluate()
 159 {
 160     DynamicList<point2D, 8> clippedPolygon(0);
 161     DynamicList<bool, 8> clippedVertexInside(0);
 162     point2D interSectPt;
 163
 164     label S = subjectPolygon_.size() - 1;
 165
 166     forAll (subjectPolygon_, sPI)
 167     {
 168         label P  = sPI;
 169
 170         // if P is inside or visible from cutting edge
 171         if (isPointVisibleFromCuttingSegment(subjectPolygon_[P]))
 172         {
 173             // if S is inside or visible from cutting edge
 174             if (isPointVisibleFromCuttingSegment(subjectPolygon_[S]))
 175             {
 176                 clippedPolygon.append(subjectPolygon_[P]);  // Output P
 177             }
 178             else
 179             {
 180                 // Output intersection between S and P
 181                 if
 182                 (
 183                     clipSegment
 184                     (
 185                         subjectPolygon_[S],
 186                         subjectPolygon_[P],
 187                         interSectPt
 188                     )
 189                 )
 190                 {
 191                     // We do intersect this edge
 192
 193                     // Output intersection point
 194                     clippedPolygon.append(interSectPt);
 195
 196                     // Build the inside/outside list
 197                     clippedVertexInside.append(true);
 198                 }
 199                 else
 200                 {
 201 //                     WarningIn("List<point2D> SutherlandHodgman::evaluate()")
 202 //                         << "missed intersection... recheck algorithm, 1"
 203 //                         << endl;
 204                 }
 205
 206                 // Output P
 207                 clippedPolygon.append(subjectPolygon_[P]);
 208             }
 209         }
 210         // if S is inside or visible from cutting edge
 211         else if (isPointVisibleFromCuttingSegment(subjectPolygon_[S]))
 212         {
 213             if
 214             (
 215                 clipSegment
 216                 (
 217                     subjectPolygon_[P],
 218                     subjectPolygon_[S],
 219                     interSectPt
 220                 )
 221             )
 222             {
 223                 // We do intersect this edge
 224
 225                 // Output intersection point
 226                 clippedPolygon.append(interSectPt);
 227
 228                 // Build the inside/outside list
 229                 clippedVertexInside.append(true);
 230             }
 231             else
 232             {
 233 //                 WarningIn("List<point2D> SutherlandHodgman::evaluate()")
 234 //                     << "missed intersection... recheck algorithm, 2"
 235 //                     << endl;
 236             }
 237         }
 238
 239         S = P;
 240     }
 241
 242     // Next clipping edge
 243     currentClipEdgeP1_ = currentClipEdgeP2_;
 244     currentClipEdgeP2_++;
 245
 246     // Update, clean up and return
 247     subjectPolygon_.transfer(clippedPolygon.shrink());
 248
 249     // Are we done?
 250     if (currentClipEdgeP2_ < clippingPolygon_.size())
 251     {
 252         // Clip against next edge, re-entrant.
 253         evaluate();
 254     }
 255
 256     // We are done
 257     return subjectPolygon_; // Which is now totally clipped
 258 }
 259
 260
 261 // ************************************************************************* //