Branch data Line data Source code
1 : : #include "meshkit/AF2FreeZone.hpp"
2 : :
3 : : #include <algorithm>
4 : : #include <math.h>
5 : 264590 : AF2FreeZone::AF2FreeZone(
6 : : std::list<AF2Point2D> const & bndryPoints) :
7 [ + - ]: 264590 : vertices(bndryPoints.begin(), bndryPoints.end())
8 : : {
9 : : typedef std::list<AF2Point2D>::const_iterator ItrType;
10 : :
11 : : // determine the bounding box
12 : 264590 : minX = 0;
13 : 264590 : maxX = 0;
14 : 264590 : minY = 0;
15 : 264590 : maxY = 0;
16 : 264590 : bool isFirst = true;
17 [ + - ][ + - ]: 1617980 : for (ItrType itr = vertices.begin(); itr != vertices.end(); ++itr)
[ + - ][ + - ]
[ + + ]
18 : : {
19 [ + + ]: 1353390 : if (isFirst)
20 : : {
21 [ + - ][ + - ]: 264590 : minX = itr->getX();
22 [ + - ][ + - ]: 264590 : maxX = itr->getX();
23 [ + - ][ + - ]: 264590 : minY = itr->getY();
24 [ + - ][ + - ]: 264590 : maxY = itr->getY();
25 : 264590 : isFirst = false;
26 : : }
27 : : else
28 : : {
29 [ + - ][ + - ]: 1088800 : minX = std::min(minX, itr->getX());
[ + - ]
30 [ + - ][ + - ]: 1088800 : maxX = std::max(maxX, itr->getX());
[ + - ]
31 [ + - ][ + - ]: 1088800 : minY = std::min(minY, itr->getY());
[ + - ]
32 [ + - ][ + - ]: 1088800 : maxY = std::max(maxY, itr->getY());
[ + - ]
33 : : }
34 : : }
35 : :
36 : : // determine the scale
37 [ + - ]: 264590 : scale = std::max(maxX - minX, maxY - minY);
38 [ + - ]: 264590 : scale = std::max(scale, fabs(maxX));
39 [ + - ]: 264590 : scale = std::max(scale, fabs(minX));
40 [ + - ]: 264590 : scale = std::max(scale, fabs(maxY));
41 [ + - ]: 264590 : scale = std::max(scale, fabs(minY));
42 : 264590 : }
43 : :
44 : 6091599 : bool AF2FreeZone::nearContains(AF2Point2D const & testPnt,
45 : : bool const & containsBndry) const
46 : : {
47 : : // Remark: This method could execute faster if the class precomputed
48 : : // and stored the coefficients needed for the comparisons. That is
49 : : // how netgen version 5.2 did it.
50 : : // The test point x-coefficient would be -rayYDiff.
51 : : // The test point y-coefficient would be rayXDiff.
52 : : // The constant would be ( -rayXDiff * itr->getY() ) + . . .
53 : : // ( rayYDiff * itr->getX() )
54 : :
55 : : typedef std::list<AF2Point2D>::const_iterator ItrType;
56 : :
57 : : // check whether the test point lies outside of the free zone bounding box
58 [ + + + + ]: 17194657 : if (testPnt.getX() < minX || testPnt.getX() > maxX ||
[ + + ]
59 [ + + ][ + + ]: 11103058 : testPnt.getY() < minY || testPnt.getY() > maxY)
60 : : {
61 : 5429169 : return false;
62 : : }
63 : :
64 : : // check whether the test point is clockwise (i.e., outside) relative
65 : : // to one of the bounding line segments of the free zone
66 [ + - ][ + - ]: 1894690 : for (ItrType itr = vertices.begin(); itr != vertices.end(); ++itr)
[ + + ]
67 : : {
68 : : // if the parameter is defined such that the free zone should not
69 : : // be considered to contain points approximately equal to its boundary
70 : : // points, explicitly test whether the test point is approximately
71 : : // equal to the boundary point, and, if so, count it as clockwise
72 [ + + ][ + - ]: 1801018 : if (!containsBndry && nearEqual(testPnt, *itr))
[ + - ][ + + ]
[ + + ]
73 : : {
74 : 568758 : return false;
75 : : }
76 : :
77 : : // get an iterator pointing to the point after the current point,
78 : : // wrapping around if necessary
79 : 1281286 : ItrType next = itr;
80 [ + - ]: 1281286 : ++next;
81 [ + - ][ + + ]: 1281286 : if (next == vertices.end())
82 : : {
83 : 96393 : next = vertices.begin();
84 : : }
85 : :
86 : : // calculate the ray direction
87 [ + - ][ + - ]: 1281286 : double rayXDiff = next->getX() - itr->getX();
[ + - ][ + - ]
88 [ + - ][ + - ]: 1281286 : double rayYDiff = next->getY() - itr->getY();
[ + - ][ + - ]
89 : :
90 : : // test whether the test point is clockwise from the ray from the
91 : : // current point to the next point
92 [ + - ][ + - ]: 1281286 : double testXDiff = (testPnt.getX() - itr->getX());
[ + - ]
93 [ + - ][ + - ]: 1281286 : double testYDiff = (testPnt.getY() - itr->getY());
[ + - ]
94 [ + + ]: 1281286 : if (testYDiff*rayXDiff - rayYDiff*testXDiff < -1.0e-14*scale*scale)
95 : : {
96 : : // if this condition is true (test would be < 0 in perfect arithmetic,
97 : : // but want to guarantee the point is outside here), then the
98 : : // test point is clockwise from the ray
99 : : // (== 0 would mean collinear in perfect arithmetic)
100 : 49026 : return false;
101 : : }
102 : : }
103 : :
104 : 6091599 : return true;
105 : : }
106 : :
107 : 7305090 : bool AF2FreeZone::nearEqual(AF2Point2D const & pointAlpha,
108 : : AF2Point2D const & pointBravo) const
109 : : {
110 : 7305090 : double xDiff = pointAlpha.getX() - pointBravo.getX();
111 : 7305090 : double yDiff = pointAlpha.getY() - pointBravo.getY();
112 : 7305090 : return xDiff*xDiff + yDiff*yDiff < 2.0e-28*scale*scale;
113 : : }
114 : :
115 : 3682243 : bool AF2FreeZone::nearIntersects(AF2Point2D const & startPoint,
116 : : AF2Point2D const & endPoint, bool const & containsBndry) const
117 : : {
118 : : // Remark: This method could execute faster if the class precomputed
119 : : // and stored the coefficients needed for the comparisons. That is
120 : : // how netgen version 5.2 did it.
121 : : // The edge endpoints x-coefficient would be -rayYDiff.
122 : : // The edge endpoints y-coefficient would be rayXDiff.
123 : : // The constant would be ( -rayXDiff * itr->getY() ) + . . .
124 : : // ( rayYDiff * itr->getX() )
125 : :
126 : : typedef std::list<AF2Point2D>::const_iterator ItrType;
127 : :
128 : : // check whether both endpoints lie outside of the free zone bounding box
129 [ + + + + ]: 11191539 : if ((startPoint.getX() < minX && endPoint.getX() < minX) ||
[ + + ]
130 [ + + + + ]: 4743508 : (startPoint.getX() > maxX && endPoint.getX() > maxX) ||
131 [ + + ]: 8780244 : (startPoint.getY() < minY && endPoint.getY() < minY) ||
[ + + + + ]
132 [ + + ]: 1255767 : (startPoint.getY() > maxY && endPoint.getY() > maxY))
133 : : {
134 : 3082342 : return false;
135 : : }
136 : :
137 : : // check whether both endpoints of the edge are clockwise (i.e., outside)
138 : : // relative to one of the bounding line segments of the free zone
139 [ + - ][ + - ]: 1754648 : for (ItrType itr = vertices.begin(); itr != vertices.end(); ++itr)
[ + + ]
140 : : {
141 : : // get an iterator pointing to the point after the current point,
142 : : // wrapping around if necessary
143 : 1671345 : ItrType next = itr;
144 [ + - ]: 1671345 : ++next;
145 [ + - ][ + + ]: 1671345 : if (next == vertices.end())
146 : : {
147 : 185828 : next = vertices.begin();
148 : : }
149 : :
150 : : // establish offsets that can be used to avoid detecting intersections
151 : : // with the boundary as intersections
152 : 1671345 : double startXOffset = 0.0;
153 : 1671345 : double startYOffset = 0.0;
154 : 1671345 : double endXOffset = 0.0;
155 : 1671345 : double endYOffset = 0.0;
156 : 1671345 : bool offsetStart = false;
157 : 1671345 : bool offsetEnd = false;
158 : :
159 [ + + ]: 1671345 : if (!containsBndry)
160 : : {
161 [ + - ][ + - ]: 1671237 : if (nearEqual(startPoint, *itr))
[ + + ]
162 : : {
163 [ + - ][ + - ]: 375293 : if (nearEqual(endPoint, *next))
[ + + ]
164 : : {
165 : : // the query line segment matches a boundary line segment
166 : 516598 : return false;
167 : : }
168 : : // the start point is approximately equal to a boundary point
169 : : // so it should be offset in the direction of the end point
170 : 123239 : offsetStart = true;
171 : : }
172 [ + - ][ + - ]: 1295944 : else if (nearEqual(startPoint, *next))
[ + + ]
173 : : {
174 [ + - ][ + - ]: 253895 : if (nearEqual(endPoint, *itr))
[ + + ]
175 : : {
176 : : // the query line segment matches a boundary line segment
177 : 24 : return false;
178 : : }
179 : : // the start point is approximately equal to a boundary point
180 : : // so it should be offset in the direction of the end point
181 : 253871 : offsetStart = true;
182 : : }
183 [ + - ][ + - ]: 1042049 : else if (nearEqual(endPoint, *itr) || nearEqual(endPoint, *next))
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ]
184 : : {
185 : 334887 : offsetEnd = true;
186 : : }
187 : :
188 [ + + ][ + + ]: 1419159 : if (offsetStart || offsetEnd)
189 : : {
190 [ + - ][ + - ]: 711997 : double testSegXDiff = endPoint.getX() - startPoint.getX();
191 [ + - ][ + - ]: 711997 : double testSegYDiff = endPoint.getY() - startPoint.getY();
192 : 711997 : double testSegLength = sqrt(testSegXDiff * testSegXDiff +
193 : 711997 : testSegYDiff * testSegYDiff);
194 [ - + ]: 711997 : if (testSegLength == 0)
195 : : {
196 : 0 : testSegLength = 1.0;
197 : : }
198 : : // The offset amount affects which segment directions are acceptable
199 : : // near the borderline. A factor of 2.0e-10 usually allows the
200 : : // clockwise test to pass if the test segment makes an angle of more
201 : : // than 1.0e-3 radians (in the appropriate direction) with the ray.
202 [ + + ]: 711997 : if (offsetStart)
203 : : {
204 : 377110 : startXOffset = 2.0e-10 * scale * testSegXDiff / testSegLength;
205 : 377110 : startYOffset = 2.0e-10 * scale * testSegYDiff / testSegLength;
206 : : }
207 [ + + ]: 711997 : if (offsetEnd)
208 : : {
209 : 334887 : endXOffset = -2.0e-10 * scale * testSegXDiff / testSegLength;
210 : 1419159 : endYOffset = -2.0e-10 * scale * testSegYDiff / testSegLength;
211 : : }
212 : : }
213 : : }
214 : :
215 : : // calculate the ray direction
216 [ + - ][ + - ]: 1419267 : double rayXDiff = next->getX() - itr->getX();
[ + - ][ + - ]
217 [ + - ][ + - ]: 1419267 : double rayYDiff = next->getY() - itr->getY();
[ + - ][ + - ]
218 : :
219 : : // test whether the start point and end point are both clockwise from
220 : : // the ray from the current point to the next point
221 [ + - ][ + - ]: 1419267 : double startXDiff = (startPoint.getX() + startXOffset - itr->getX());
[ + - ]
222 [ + - ][ + - ]: 1419267 : double startYDiff = (startPoint.getY() + startYOffset - itr->getY());
[ + - ]
223 [ + - ][ + - ]: 1419267 : double endXDiff = (endPoint.getX() + endXOffset - itr->getX());
[ + - ]
224 [ + - ][ + - ]: 1419267 : double endYDiff = (endPoint.getY() + endYOffset - itr->getY());
[ + - ]
225 [ + + ][ + + ]: 1419267 : if ((startYDiff*rayXDiff - rayYDiff*startXDiff < -1.0e-14*scale*scale) &&
226 : 335526 : (endYDiff*rayXDiff - rayYDiff*endXDiff < -1.0e-14*scale*scale))
227 : : {
228 : : // if these condition are true (test would be < 0 in perfect arithmetic,
229 : : // but want to guarantee the point is outside here), then the
230 : : // whole edge is clockwise from the ray
231 : : // (== 0 would mean collinear in perfect arithmetic)
232 : 264520 : return false;
233 : : }
234 : : }
235 : :
236 : : // Check whether the whole free zone is either clockwise (right) or
237 : : // counterclockwise (left) of the test line segment
238 : 83303 : bool allCW = true;
239 : 83303 : bool allCCW = true;
240 : : // First calculate the x and y difference of the test line segment
241 : 83303 : double testSegXDiff = (endPoint.getX() - startPoint.getX());
242 : 83303 : double testSegYDiff = (endPoint.getY() - startPoint.getY());
243 [ + - ][ + + ]: 724662 : for (ItrType itr = vertices.begin();
244 [ + - ][ + + ]: 483108 : itr != vertices.end() && (allCW || allCCW); ++itr)
[ + + ][ + + ]
[ - + ][ # # ]
245 : : {
246 : :
247 : : // calculate the difference from the start point to the
248 : : // free zone vertex
249 [ + - ][ + - ]: 158251 : double fzvXDiff = itr->getX() - startPoint.getX();
[ + - ]
250 [ + - ][ + - ]: 158251 : double fzvYDiff = itr->getY() - startPoint.getY();
[ + - ]
251 : :
252 : : // test whether the free zone vertex is clockwise from (right of)
253 : : // the test line segment
254 [ + + ]: 158251 : if (fzvYDiff*testSegXDiff - testSegYDiff*fzvXDiff > -1.0e-14*scale*scale)
255 : : {
256 : : // if this condition is true (test would be >= 0 in perfect arithmetic,
257 : : // but want to guarantee the point is clockwise if the condition
258 : : // is false), then the free zone vertex is (probably) not clockwise
259 : : // from the test line segment.
260 : 144571 : allCW = false;
261 : : }
262 : :
263 : : // test whether the free zone vertex is counterclockwise from (left of)
264 : : // the test line segment
265 [ + + ]: 158251 : if (fzvYDiff*testSegXDiff - testSegYDiff*fzvXDiff < 1.0e-14*scale*scale)
266 : : {
267 : : // if this condition is true (test would be <= 0 in perfect arithmetic,
268 : : // but want to guarantee the point is counterclockwise if the condition
269 : : // is false), then the free zone vertex is (probably) not
270 : : // counterclockwise from the test line segment.
271 : 87952 : allCCW = false;
272 : : }
273 : : }
274 : :
275 [ + + ][ + + ]: 3682243 : return ( !allCW && !allCCW );
276 : : }
277 : :
278 : 264290 : bool AF2FreeZone::isConvex() const
279 : : {
280 : : typedef std::list<AF2Point2D>::const_iterator ItrType;
281 : :
282 [ + - ][ + - ]: 1616000 : for (ItrType itr = vertices.begin(); itr != vertices.end(); ++itr)
[ + + ]
283 : : {
284 : : // get an iterator pointing to the point after the current point,
285 : : // wrapping around if necessary
286 : 1351785 : ItrType next = itr;
287 [ + - ]: 1351785 : ++next;
288 [ + - ][ + + ]: 1351785 : if (next == vertices.end())
289 : : {
290 : 264215 : next = vertices.begin();
291 : : }
292 : :
293 : : // calculate the ray direction
294 [ + - ][ + - ]: 1351785 : double rayXDiff = next->getX() - itr->getX();
[ + - ][ + - ]
295 [ + - ][ + - ]: 1351785 : double rayYDiff = next->getY() - itr->getY();
[ + - ][ + - ]
296 : :
297 : : // test whether all other points are counterclockwise
298 [ + - + - ]: 16781182 : for (ItrType testPnt = vertices.begin();
[ + + ]
299 : 8390591 : testPnt != vertices.end(); ++testPnt)
300 : : {
301 : : // skip the test if the test point is the current point or the next point
302 [ + - ][ + - ]: 7038881 : if ((*testPnt) == (*itr) || (*testPnt) == (*next))
[ + - ][ + + ]
[ + - ][ + - ]
[ + - ][ + + ]
[ + + ]
303 : : {
304 : 2703562 : continue;
305 : : }
306 : :
307 : : // test that the test point is counterclockwise from the ray from the
308 : : // current point to the next point
309 [ + - ][ + - ]: 4335319 : double testXDiff = (testPnt->getX() - itr->getX());
[ + - ][ + - ]
310 [ + - ][ + - ]: 4335319 : double testYDiff = (testPnt->getY() - itr->getY());
[ + - ][ + - ]
311 [ + + ]: 4335319 : if (testYDiff*rayXDiff - rayYDiff*testXDiff < -1.0e-14*scale*scale)
312 : : {
313 : : // if this condition is true (test would be < 0 in perfect arithmetic),
314 : : // then the test point is clockwise from the ray
315 : : // (== 0 would mean collinear in perfect arithmetic)
316 : 75 : return false;
317 : : }
318 : : }
319 : : }
320 : :
321 : 264290 : return true;
322 : : }
|
