MeshKit  1.0
AF2FreeZone.cpp
Go to the documentation of this file.
00001 #include "meshkit/AF2FreeZone.hpp"
00002 
00003 #include <algorithm>
00004 #include <math.h>
00005 AF2FreeZone::AF2FreeZone(
00006     std::list<AF2Point2D> const & bndryPoints) :
00007     vertices(bndryPoints.begin(), bndryPoints.end())
00008 {
00009   typedef std::list<AF2Point2D>::const_iterator ItrType;
00010 
00011   // determine the bounding box
00012   minX = 0;
00013   maxX = 0;
00014   minY = 0;
00015   maxY = 0;
00016   bool isFirst = true;
00017   for (ItrType itr = vertices.begin(); itr != vertices.end(); ++itr)
00018   {
00019     if (isFirst)
00020     {
00021       minX = itr->getX();
00022       maxX = itr->getX();
00023       minY = itr->getY();
00024       maxY = itr->getY();
00025       isFirst = false;
00026     }
00027     else
00028     {
00029       minX = std::min(minX, itr->getX());
00030       maxX = std::max(maxX, itr->getX());
00031       minY = std::min(minY, itr->getY());
00032       maxY = std::max(maxY, itr->getY());
00033     }
00034   }
00035 
00036   // determine the scale
00037   scale = std::max(maxX - minX, maxY - minY);
00038   scale = std::max(scale, fabs(maxX));
00039   scale = std::max(scale, fabs(minX));
00040   scale = std::max(scale, fabs(maxY));
00041   scale = std::max(scale, fabs(minY));
00042 }
00043 
00044 bool AF2FreeZone::nearContains(AF2Point2D const & testPnt,
00045     bool const & containsBndry) const
00046 {
00047   // Remark: This method could execute faster if the class precomputed
00048   // and stored the coefficients needed for the comparisons.  That is
00049   // how netgen version 5.2 did it.
00050   // The test point x-coefficient would be -rayYDiff.
00051   // The test point y-coefficient would be rayXDiff.
00052   // The constant would be ( -rayXDiff * itr->getY() ) + . . .
00053   //     ( rayYDiff * itr->getX() )
00054 
00055   typedef std::list<AF2Point2D>::const_iterator ItrType;
00056 
00057   // check whether the test point lies outside of the free zone bounding box
00058   if (testPnt.getX() < minX || testPnt.getX() > maxX ||
00059       testPnt.getY() < minY || testPnt.getY() > maxY)
00060   {
00061     return false;
00062   }
00063 
00064   // check whether the test point is clockwise (i.e., outside) relative
00065   // to one of the bounding line segments of the free zone
00066   for (ItrType itr = vertices.begin(); itr != vertices.end(); ++itr)
00067   {
00068     // if the parameter is defined such that the free zone should not
00069     // be considered to contain points approximately equal to its boundary
00070     // points, explicitly test whether the test point is approximately
00071     // equal to the boundary point, and, if so, count it as clockwise
00072     if (!containsBndry && nearEqual(testPnt, *itr))
00073     {
00074       return false;
00075     }
00076  
00077     // get an iterator pointing to the point after the current point,
00078     // wrapping around if necessary
00079     ItrType next = itr;
00080     ++next;
00081     if (next == vertices.end())
00082     {
00083       next = vertices.begin();
00084     }
00085 
00086     // calculate the ray direction
00087     double rayXDiff = next->getX() - itr->getX();
00088     double rayYDiff = next->getY() - itr->getY();
00089 
00090     // test whether the test point is clockwise from the ray from the
00091     // current point to the next point
00092     double testXDiff = (testPnt.getX() - itr->getX());
00093     double testYDiff = (testPnt.getY() - itr->getY());
00094     if (testYDiff*rayXDiff - rayYDiff*testXDiff < -1.0e-14*scale*scale)
00095     {
00096       // if this condition is true (test would be < 0 in perfect arithmetic,
00097       // but want to guarantee the point is outside here), then the
00098       // test point is clockwise from the ray
00099       // (== 0 would mean collinear in perfect arithmetic)
00100       return false;
00101     }
00102   }
00103 
00104   return true;
00105 }
00106 
00107 bool AF2FreeZone::nearEqual(AF2Point2D const & pointAlpha,
00108     AF2Point2D const & pointBravo) const
00109 {
00110   double xDiff = pointAlpha.getX() - pointBravo.getX();
00111   double yDiff = pointAlpha.getY() - pointBravo.getY();
00112   return xDiff*xDiff + yDiff*yDiff < 2.0e-28*scale*scale;
00113 }
00114 
00115 bool AF2FreeZone::nearIntersects(AF2Point2D const & startPoint,
00116     AF2Point2D const & endPoint, bool const & containsBndry) const
00117 {
00118   // Remark: This method could execute faster if the class precomputed
00119   // and stored the coefficients needed for the comparisons.  That is
00120   // how netgen version 5.2 did it.
00121   // The edge endpoints x-coefficient would be -rayYDiff.
00122   // The edge endpoints y-coefficient would be rayXDiff.
00123   // The constant would be ( -rayXDiff * itr->getY() ) + . . .
00124   //     ( rayYDiff * itr->getX() )
00125 
00126   typedef std::list<AF2Point2D>::const_iterator ItrType;
00127 
00128   // check whether both endpoints lie outside of the free zone bounding box
00129   if ((startPoint.getX() < minX && endPoint.getX() < minX) ||
00130       (startPoint.getX() > maxX && endPoint.getX() > maxX) ||
00131       (startPoint.getY() < minY && endPoint.getY() < minY) ||
00132       (startPoint.getY() > maxY && endPoint.getY() > maxY))
00133   {
00134     return false;
00135   }
00136 
00137   // check whether both endpoints of the edge are clockwise (i.e., outside)
00138   // relative to one of the bounding line segments of the free zone
00139   for (ItrType itr = vertices.begin(); itr != vertices.end(); ++itr)
00140   {
00141     // get an iterator pointing to the point after the current point,
00142     // wrapping around if necessary
00143     ItrType next = itr;
00144     ++next;
00145     if (next == vertices.end())
00146     {
00147       next = vertices.begin();
00148     }
00149 
00150     // establish offsets that can be used to avoid detecting intersections
00151     // with the boundary as intersections
00152     double startXOffset = 0.0;
00153     double startYOffset = 0.0;
00154     double endXOffset = 0.0;
00155     double endYOffset = 0.0;
00156     bool offsetStart = false;
00157     bool offsetEnd = false;
00158 
00159     if (!containsBndry)
00160     {
00161       if (nearEqual(startPoint, *itr))
00162       {
00163         if (nearEqual(endPoint, *next))
00164         {
00165           // the query line segment matches a boundary line segment
00166           return false;
00167         }
00168         // the start point is approximately equal to a boundary point
00169         // so it should be offset in the direction of the end point
00170         offsetStart = true;
00171       }
00172       else if (nearEqual(startPoint, *next))
00173       {
00174         if (nearEqual(endPoint, *itr))
00175         {
00176           // the query line segment matches a boundary line segment
00177           return false;
00178         }
00179         // the start point is approximately equal to a boundary point
00180         // so it should be offset in the direction of the end point
00181         offsetStart = true;
00182       }
00183       else if (nearEqual(endPoint, *itr) || nearEqual(endPoint, *next))
00184       {
00185         offsetEnd = true;
00186       }
00187 
00188       if (offsetStart || offsetEnd)
00189       {
00190         double testSegXDiff = endPoint.getX() - startPoint.getX();
00191         double testSegYDiff = endPoint.getY() - startPoint.getY();
00192         double testSegLength = sqrt(testSegXDiff * testSegXDiff  +
00193             testSegYDiff * testSegYDiff);
00194         if (testSegLength == 0)
00195         {
00196           testSegLength = 1.0;
00197         }
00198         // The offset amount affects which segment directions are acceptable
00199         // near the borderline.  A factor of 2.0e-10 usually allows the
00200         // clockwise test to pass if the test segment makes an angle of more
00201         // than 1.0e-3 radians (in the appropriate direction) with the ray.
00202         if (offsetStart)
00203         {
00204           startXOffset = 2.0e-10 * scale * testSegXDiff / testSegLength;
00205           startYOffset = 2.0e-10 * scale * testSegYDiff / testSegLength;
00206         }
00207         if (offsetEnd)
00208         {
00209           endXOffset = -2.0e-10 * scale * testSegXDiff / testSegLength;
00210           endYOffset = -2.0e-10 * scale * testSegYDiff / testSegLength;
00211         }
00212       }
00213     }
00214 
00215     // calculate the ray direction
00216     double rayXDiff = next->getX() - itr->getX();
00217     double rayYDiff = next->getY() - itr->getY();
00218 
00219     // test whether the start point and end point are both clockwise from
00220     // the ray from the current point to the next point
00221     double startXDiff = (startPoint.getX() + startXOffset - itr->getX());
00222     double startYDiff = (startPoint.getY() + startYOffset - itr->getY());
00223     double endXDiff = (endPoint.getX() + endXOffset - itr->getX());
00224     double endYDiff = (endPoint.getY() + endYOffset - itr->getY());
00225     if ((startYDiff*rayXDiff - rayYDiff*startXDiff < -1.0e-14*scale*scale) &&
00226         (endYDiff*rayXDiff - rayYDiff*endXDiff < -1.0e-14*scale*scale))
00227     {
00228       // if these condition are true (test would be < 0 in perfect arithmetic,
00229       // but want to guarantee the point is outside here), then the
00230       // whole edge is clockwise from the ray
00231       // (== 0 would mean collinear in perfect arithmetic)
00232       return false;
00233     }
00234   }
00235 
00236   // Check whether the whole free zone is either clockwise (right) or
00237   // counterclockwise (left) of the test line segment
00238   bool allCW = true;
00239   bool allCCW = true;
00240   // First calculate the x and y difference of the test line segment
00241   double testSegXDiff = (endPoint.getX() - startPoint.getX());
00242   double testSegYDiff = (endPoint.getY() - startPoint.getY());
00243   for (ItrType itr = vertices.begin();
00244       itr != vertices.end() && (allCW || allCCW); ++itr)
00245   {
00246 
00247     // calculate the difference from the start point to the
00248     // free zone vertex
00249     double fzvXDiff = itr->getX() - startPoint.getX();
00250     double fzvYDiff = itr->getY() - startPoint.getY();
00251 
00252     // test whether the free zone vertex is clockwise from (right of)
00253     // the test line segment
00254     if (fzvYDiff*testSegXDiff - testSegYDiff*fzvXDiff > -1.0e-14*scale*scale)
00255     {
00256       // if this condition is true (test would be >= 0 in perfect arithmetic,
00257       // but want to guarantee the point is clockwise if the condition
00258       // is false), then the free zone vertex is (probably) not clockwise
00259       // from the test line segment.
00260       allCW = false;
00261     }
00262 
00263     // test whether the free zone vertex is counterclockwise from (left of)
00264     // the test line segment
00265     if (fzvYDiff*testSegXDiff - testSegYDiff*fzvXDiff < 1.0e-14*scale*scale)
00266     {
00267       // if this condition is true (test would be <= 0 in perfect arithmetic,
00268       // but want to guarantee the point is counterclockwise if the condition
00269       // is false), then the free zone vertex is (probably) not
00270       // counterclockwise from the test line segment.
00271       allCCW = false;
00272     }
00273   }
00274 
00275   return ( !allCW && !allCCW );
00276 }
00277 
00278 bool AF2FreeZone::isConvex() const
00279 {
00280   typedef std::list<AF2Point2D>::const_iterator ItrType;
00281 
00282   for (ItrType itr = vertices.begin(); itr != vertices.end(); ++itr)
00283   {
00284     // get an iterator pointing to the point after the current point,
00285     // wrapping around if necessary
00286     ItrType next = itr;
00287     ++next;
00288     if (next == vertices.end())
00289     {
00290       next = vertices.begin();
00291     }
00292 
00293     // calculate the ray direction
00294     double rayXDiff = next->getX() - itr->getX();
00295     double rayYDiff = next->getY() - itr->getY();
00296 
00297     // test whether all other points are counterclockwise
00298     for (ItrType testPnt = vertices.begin();
00299         testPnt != vertices.end(); ++testPnt)
00300     {
00301       // skip the test if the test point is the current point or the next point
00302       if ((*testPnt) == (*itr) || (*testPnt) == (*next))
00303       {
00304         continue;
00305       }
00306 
00307       // test that the test point is counterclockwise from the ray from the
00308       // current point to the next point
00309       double testXDiff = (testPnt->getX() - itr->getX());
00310       double testYDiff = (testPnt->getY() - itr->getY());
00311       if (testYDiff*rayXDiff - rayYDiff*testXDiff < -1.0e-14*scale*scale)
00312       {
00313         // if this condition is true (test would be < 0 in perfect arithmetic),
00314         // then the test point is clockwise from the ray
00315         // (== 0 would mean collinear in perfect arithmetic)
00316         return false;
00317       }
00318     }
00319   }
00320 
00321   return true;
00322 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines