1 // This file contains liberal use of asserts to assist code development and
2 // debugging. Standard matplotlib builds disable asserts so they cause no
3 // performance reduction. To enable the asserts, you need to undefine the
4 // NDEBUG macro, which is achieved by adding the following
5 // undef_macros=['NDEBUG']
6 // to the appropriate make_extension call in setupext.py, and then rebuilding.
7 #define NO_IMPORT_ARRAY
8
9 #include "src/mplutils.h"
10 #include "src/_contour.h"
11 #include <algorithm>
12
13
14 // 'kind' codes.
15 #define MOVETO 1
16 #define LINETO 2
17 #define CLOSEPOLY 79
18
19 // Point indices from current quad index.
20 #define POINT_SW (quad)
21 #define POINT_SE (quad+1)
22 #define POINT_NW (quad+_nx)
23 #define POINT_NE (quad+_nx+1)
24
25 // CacheItem masks, only accessed directly to set. To read, use accessors
26 // detailed below. 1 and 2 refer to level indices (lower and upper).
27 #define MASK_Z_LEVEL 0x0003 // Combines the following two.
28 #define MASK_Z_LEVEL_1 0x0001 // z > lower_level.
29 #define MASK_Z_LEVEL_2 0x0002 // z > upper_level.
30 #define MASK_VISITED_1 0x0004 // Algorithm has visited this quad.
31 #define MASK_VISITED_2 0x0008
32 #define MASK_SADDLE_1 0x0010 // quad is a saddle quad.
33 #define MASK_SADDLE_2 0x0020
34 #define MASK_SADDLE_LEFT_1 0x0040 // Contours turn left at saddle quad.
35 #define MASK_SADDLE_LEFT_2 0x0080
36 #define MASK_SADDLE_START_SW_1 0x0100 // Next visit starts on S or W edge.
37 #define MASK_SADDLE_START_SW_2 0x0200
38 #define MASK_BOUNDARY_S 0x0400 // S edge of quad is a boundary.
39 #define MASK_BOUNDARY_W 0x0800 // W edge of quad is a boundary.
40 // EXISTS_QUAD bit is always used, but the 4 EXISTS_CORNER are only used if
41 // _corner_mask is true. Only one of EXISTS_QUAD or EXISTS_??_CORNER is ever
42 // set per quad, hence not using unique bits for each; care is needed when
43 // testing for these flags as they overlap.
44 #define MASK_EXISTS_QUAD 0x1000 // All of quad exists (is not masked).
45 #define MASK_EXISTS_SW_CORNER 0x2000 // SW corner exists, NE corner is masked.
46 #define MASK_EXISTS_SE_CORNER 0x3000
47 #define MASK_EXISTS_NW_CORNER 0x4000
48 #define MASK_EXISTS_NE_CORNER 0x5000
49 #define MASK_EXISTS 0x7000 // Combines all 5 EXISTS masks.
50
51 // The following are only needed for filled contours.
52 #define MASK_VISITED_S 0x10000 // Algorithm has visited S boundary.
53 #define MASK_VISITED_W 0x20000 // Algorithm has visited W boundary.
54 #define MASK_VISITED_CORNER 0x40000 // Algorithm has visited corner edge.
55
56
57 // Accessors for various CacheItem masks. li is shorthand for level_index.
58 #define Z_LEVEL(quad) (_cache[quad] & MASK_Z_LEVEL)
59 #define Z_NE Z_LEVEL(POINT_NE)
60 #define Z_NW Z_LEVEL(POINT_NW)
61 #define Z_SE Z_LEVEL(POINT_SE)
62 #define Z_SW Z_LEVEL(POINT_SW)
63 #define VISITED(quad,li) (_cache[quad] & (li==1 ? MASK_VISITED_1 : MASK_VISITED_2))
64 #define VISITED_S(quad) (_cache[quad] & MASK_VISITED_S)
65 #define VISITED_W(quad) (_cache[quad] & MASK_VISITED_W)
66 #define VISITED_CORNER(quad) (_cache[quad] & MASK_VISITED_CORNER)
67 #define SADDLE(quad,li) (_cache[quad] & (li==1 ? MASK_SADDLE_1 : MASK_SADDLE_2))
68 #define SADDLE_LEFT(quad,li) (_cache[quad] & (li==1 ? MASK_SADDLE_LEFT_1 : MASK_SADDLE_LEFT_2))
69 #define SADDLE_START_SW(quad,li) (_cache[quad] & (li==1 ? MASK_SADDLE_START_SW_1 : MASK_SADDLE_START_SW_2))
70 #define BOUNDARY_S(quad) (_cache[quad] & MASK_BOUNDARY_S)
71 #define BOUNDARY_W(quad) (_cache[quad] & MASK_BOUNDARY_W)
72 #define BOUNDARY_N(quad) BOUNDARY_S(quad+_nx)
73 #define BOUNDARY_E(quad) BOUNDARY_W(quad+1)
74 #define EXISTS_QUAD(quad) ((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_QUAD)
75 #define EXISTS_NONE(quad) ((_cache[quad] & MASK_EXISTS) == 0)
76 // The following are only used if _corner_mask is true.
77 #define EXISTS_SW_CORNER(quad) ((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_SW_CORNER)
78 #define EXISTS_SE_CORNER(quad) ((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_SE_CORNER)
79 #define EXISTS_NW_CORNER(quad) ((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_NW_CORNER)
80 #define EXISTS_NE_CORNER(quad) ((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_NE_CORNER)
81 #define EXISTS_ANY_CORNER(quad) (!EXISTS_NONE(quad) && !EXISTS_QUAD(quad))
82 #define EXISTS_W_EDGE(quad) (EXISTS_QUAD(quad) || EXISTS_SW_CORNER(quad) || EXISTS_NW_CORNER(quad))
83 #define EXISTS_E_EDGE(quad) (EXISTS_QUAD(quad) || EXISTS_SE_CORNER(quad) || EXISTS_NE_CORNER(quad))
84 #define EXISTS_S_EDGE(quad) (EXISTS_QUAD(quad) || EXISTS_SW_CORNER(quad) || EXISTS_SE_CORNER(quad))
85 #define EXISTS_N_EDGE(quad) (EXISTS_QUAD(quad) || EXISTS_NW_CORNER(quad) || EXISTS_NE_CORNER(quad))
86 // Note that EXISTS_NE_CORNER(quad) is equivalent to BOUNDARY_SW(quad), etc.
87
88
89
QuadEdge()90 QuadEdge::QuadEdge()
91 : quad(-1), edge(Edge_None)
92 {}
93
QuadEdge(long quad_,Edge edge_)94 QuadEdge::QuadEdge(long quad_, Edge edge_)
95 : quad(quad_), edge(edge_)
96 {}
97
operator <(const QuadEdge & other) const98 bool QuadEdge::operator<(const QuadEdge& other) const
99 {
100 if (quad != other.quad)
101 return quad < other.quad;
102 else
103 return edge < other.edge;
104 }
105
operator ==(const QuadEdge & other) const106 bool QuadEdge::operator==(const QuadEdge& other) const
107 {
108 return quad == other.quad && edge == other.edge;
109 }
110
operator !=(const QuadEdge & other) const111 bool QuadEdge::operator!=(const QuadEdge& other) const
112 {
113 return !operator==(other);
114 }
115
operator <<(std::ostream & os,const QuadEdge & quad_edge)116 std::ostream& operator<<(std::ostream& os, const QuadEdge& quad_edge)
117 {
118 return os << quad_edge.quad << ' ' << quad_edge.edge;
119 }
120
121
122
XY()123 XY::XY()
124 {}
125
XY(const double & x_,const double & y_)126 XY::XY(const double& x_, const double& y_)
127 : x(x_), y(y_)
128 {}
129
operator ==(const XY & other) const130 bool XY::operator==(const XY& other) const
131 {
132 return x == other.x && y == other.y;
133 }
134
operator !=(const XY & other) const135 bool XY::operator!=(const XY& other) const
136 {
137 return x != other.x || y != other.y;
138 }
139
operator *(const double & multiplier) const140 XY XY::operator*(const double& multiplier) const
141 {
142 return XY(x*multiplier, y*multiplier);
143 }
144
operator +=(const XY & other)145 const XY& XY::operator+=(const XY& other)
146 {
147 x += other.x;
148 y += other.y;
149 return *this;
150 }
151
operator -=(const XY & other)152 const XY& XY::operator-=(const XY& other)
153 {
154 x -= other.x;
155 y -= other.y;
156 return *this;
157 }
158
operator +(const XY & other) const159 XY XY::operator+(const XY& other) const
160 {
161 return XY(x + other.x, y + other.y);
162 }
163
operator -(const XY & other) const164 XY XY::operator-(const XY& other) const
165 {
166 return XY(x - other.x, y - other.y);
167 }
168
operator <<(std::ostream & os,const XY & xy)169 std::ostream& operator<<(std::ostream& os, const XY& xy)
170 {
171 return os << '(' << xy.x << ' ' << xy.y << ')';
172 }
173
174
175
ContourLine(bool is_hole)176 ContourLine::ContourLine(bool is_hole)
177 : std::vector<XY>(),
178 _is_hole(is_hole),
179 _parent(0)
180 {}
181
add_child(ContourLine * child)182 void ContourLine::add_child(ContourLine* child)
183 {
184 assert(!_is_hole && "Cannot add_child to a hole");
185 assert(child != 0 && "Null child ContourLine");
186 _children.push_back(child);
187 }
188
clear_parent()189 void ContourLine::clear_parent()
190 {
191 assert(is_hole() && "Cannot clear parent of non-hole");
192 assert(_parent != 0 && "Null parent ContourLine");
193 _parent = 0;
194 }
195
get_children() const196 const ContourLine::Children& ContourLine::get_children() const
197 {
198 assert(!_is_hole && "Cannot get_children of a hole");
199 return _children;
200 }
201
get_parent() const202 const ContourLine* ContourLine::get_parent() const
203 {
204 assert(_is_hole && "Cannot get_parent of a non-hole");
205 return _parent;
206 }
207
get_parent()208 ContourLine* ContourLine::get_parent()
209 {
210 assert(_is_hole && "Cannot get_parent of a non-hole");
211 return _parent;
212 }
213
is_hole() const214 bool ContourLine::is_hole() const
215 {
216 return _is_hole;
217 }
218
push_back(const XY & point)219 void ContourLine::push_back(const XY& point)
220 {
221 if (empty() || point != back())
222 std::vector<XY>::push_back(point);
223 }
224
set_parent(ContourLine * parent)225 void ContourLine::set_parent(ContourLine* parent)
226 {
227 assert(_is_hole && "Cannot set parent of a non-hole");
228 assert(parent != 0 && "Null parent ContourLine");
229 _parent = parent;
230 }
231
write() const232 void ContourLine::write() const
233 {
234 std::cout << "ContourLine " << this << " of " << size() << " points:";
235 for (const_iterator it = begin(); it != end(); ++it)
236 std::cout << ' ' << *it;
237 if (is_hole())
238 std::cout << " hole, parent=" << get_parent();
239 else {
240 std::cout << " not hole";
241 if (!_children.empty()) {
242 std::cout << ", children=";
243 for (Children::const_iterator it = _children.begin();
244 it != _children.end(); ++it)
245 std::cout << *it << ' ';
246 }
247 }
248 std::cout << std::endl;
249 }
250
251
252
Contour()253 Contour::Contour()
254 {}
255
~Contour()256 Contour::~Contour()
257 {
258 delete_contour_lines();
259 }
260
delete_contour_lines()261 void Contour::delete_contour_lines()
262 {
263 for (iterator line_it = begin(); line_it != end(); ++line_it) {
264 delete *line_it;
265 *line_it = 0;
266 }
267 std::vector<ContourLine*>::clear();
268 }
269
write() const270 void Contour::write() const
271 {
272 std::cout << "Contour of " << size() << " lines." << std::endl;
273 for (const_iterator it = begin(); it != end(); ++it)
274 (*it)->write();
275 }
276
277
278
ParentCache(long nx,long x_chunk_points,long y_chunk_points)279 ParentCache::ParentCache(long nx, long x_chunk_points, long y_chunk_points)
280 : _nx(nx),
281 _x_chunk_points(x_chunk_points),
282 _y_chunk_points(y_chunk_points),
283 _lines(0), // Initialised when first needed.
284 _istart(0),
285 _jstart(0)
286 {
287 assert(_x_chunk_points > 0 && _y_chunk_points > 0 &&
288 "Chunk sizes must be positive");
289 }
290
get_parent(long quad)291 ContourLine* ParentCache::get_parent(long quad)
292 {
293 long index = quad_to_index(quad);
294 ContourLine* parent = _lines[index];
295 while (parent == 0) {
296 index -= _x_chunk_points;
297 assert(index >= 0 && "Failed to find parent in chunk ParentCache");
298 parent = _lines[index];
299 }
300 assert(parent != 0 && "Failed to find parent in chunk ParentCache");
301 return parent;
302 }
303
quad_to_index(long quad) const304 long ParentCache::quad_to_index(long quad) const
305 {
306 long i = quad % _nx;
307 long j = quad / _nx;
308 long index = (i-_istart) + (j-_jstart)*_x_chunk_points;
309
310 assert(i >= _istart && i < _istart + _x_chunk_points &&
311 "i-index outside chunk");
312 assert(j >= _jstart && j < _jstart + _y_chunk_points &&
313 "j-index outside chunk");
314 assert(index >= 0 && index < static_cast<long>(_lines.size()) &&
315 "ParentCache index outside chunk");
316
317 return index;
318 }
319
set_chunk_starts(long istart,long jstart)320 void ParentCache::set_chunk_starts(long istart, long jstart)
321 {
322 assert(istart >= 0 && jstart >= 0 &&
323 "Chunk start indices cannot be negative");
324 _istart = istart;
325 _jstart = jstart;
326 if (_lines.empty())
327 _lines.resize(_x_chunk_points*_y_chunk_points, 0);
328 else
329 std::fill(_lines.begin(), _lines.end(), (ContourLine*)0);
330 }
331
set_parent(long quad,ContourLine & contour_line)332 void ParentCache::set_parent(long quad, ContourLine& contour_line)
333 {
334 assert(!_lines.empty() &&
335 "Accessing ParentCache before it has been initialised");
336 long index = quad_to_index(quad);
337 if (_lines[index] == 0)
338 _lines[index] = (contour_line.is_hole() ? contour_line.get_parent()
339 : &contour_line);
340 }
341
342
343
QuadContourGenerator(const CoordinateArray & x,const CoordinateArray & y,const CoordinateArray & z,const MaskArray & mask,bool corner_mask,long chunk_size)344 QuadContourGenerator::QuadContourGenerator(const CoordinateArray& x,
345 const CoordinateArray& y,
346 const CoordinateArray& z,
347 const MaskArray& mask,
348 bool corner_mask,
349 long chunk_size)
350 : _x(x),
351 _y(y),
352 _z(z),
353 _nx(static_cast<long>(_x.dim(1))),
354 _ny(static_cast<long>(_x.dim(0))),
355 _n(_nx*_ny),
356 _corner_mask(corner_mask),
357 _chunk_size(chunk_size > 0 ? std::min(chunk_size, std::max(_nx, _ny)-1)
358 : std::max(_nx, _ny)-1),
359 _nxchunk(calc_chunk_count(_nx)),
360 _nychunk(calc_chunk_count(_ny)),
361 _chunk_count(_nxchunk*_nychunk),
362 _cache(new CacheItem[_n]),
363 _parent_cache(_nx,
364 chunk_size > 0 ? chunk_size+1 : _nx,
365 chunk_size > 0 ? chunk_size+1 : _ny)
366 {
367 assert(!_x.empty() && !_y.empty() && !_z.empty() && "Empty array");
368 assert(_y.dim(0) == _x.dim(0) && _y.dim(1) == _x.dim(1) &&
369 "Different-sized y and x arrays");
370 assert(_z.dim(0) == _x.dim(0) && _z.dim(1) == _x.dim(1) &&
371 "Different-sized z and x arrays");
372 assert((mask.empty() ||
373 (mask.dim(0) == _x.dim(0) && mask.dim(1) == _x.dim(1))) &&
374 "Different-sized mask and x arrays");
375
376 init_cache_grid(mask);
377 }
378
~QuadContourGenerator()379 QuadContourGenerator::~QuadContourGenerator()
380 {
381 delete [] _cache;
382 }
383
append_contour_line_to_vertices(ContourLine & contour_line,PyObject * vertices_list) const384 void QuadContourGenerator::append_contour_line_to_vertices(
385 ContourLine& contour_line,
386 PyObject* vertices_list) const
387 {
388 assert(vertices_list != 0 && "Null python vertices_list");
389
390 // Convert ContourLine to python equivalent, and clear it.
391 npy_intp dims[2] = {static_cast<npy_intp>(contour_line.size()), 2};
392 numpy::array_view<double, 2> line(dims);
393 npy_intp i = 0;
394 for (ContourLine::const_iterator point = contour_line.begin();
395 point != contour_line.end(); ++point, ++i) {
396 line(i, 0) = point->x;
397 line(i, 1) = point->y;
398 }
399 if (PyList_Append(vertices_list, line.pyobj_steal())) {
400 Py_XDECREF(vertices_list);
401 throw std::runtime_error("Unable to add contour line to vertices_list");
402 }
403
404 contour_line.clear();
405 }
406
append_contour_to_vertices_and_codes(Contour & contour,PyObject * vertices_list,PyObject * codes_list) const407 void QuadContourGenerator::append_contour_to_vertices_and_codes(
408 Contour& contour,
409 PyObject* vertices_list,
410 PyObject* codes_list) const
411 {
412 assert(vertices_list != 0 && "Null python vertices_list");
413 assert(codes_list != 0 && "Null python codes_list");
414
415 // Convert Contour to python equivalent, and clear it.
416 for (Contour::iterator line_it = contour.begin(); line_it != contour.end();
417 ++line_it) {
418 ContourLine& line = **line_it;
419 if (line.is_hole()) {
420 // If hole has already been converted to python its parent will be
421 // set to 0 and it can be deleted.
422 if (line.get_parent() != 0) {
423 delete *line_it;
424 *line_it = 0;
425 }
426 }
427 else {
428 // Non-holes are converted to python together with their child
429 // holes so that they are rendered correctly.
430 ContourLine::const_iterator point;
431 ContourLine::Children::const_iterator children_it;
432
433 const ContourLine::Children& children = line.get_children();
434 npy_intp npoints = static_cast<npy_intp>(line.size() + 1);
435 for (children_it = children.begin(); children_it != children.end();
436 ++children_it)
437 npoints += static_cast<npy_intp>((*children_it)->size() + 1);
438
439 npy_intp vertices_dims[2] = {npoints, 2};
440 numpy::array_view<double, 2> vertices(vertices_dims);
441 double* vertices_ptr = vertices.data();
442
443 npy_intp codes_dims[1] = {npoints};
444 numpy::array_view<unsigned char, 1> codes(codes_dims);
445 unsigned char* codes_ptr = codes.data();
446
447 for (point = line.begin(); point != line.end(); ++point) {
448 *vertices_ptr++ = point->x;
449 *vertices_ptr++ = point->y;
450 *codes_ptr++ = (point == line.begin() ? MOVETO : LINETO);
451 }
452 point = line.begin();
453 *vertices_ptr++ = point->x;
454 *vertices_ptr++ = point->y;
455 *codes_ptr++ = CLOSEPOLY;
456
457 for (children_it = children.begin(); children_it != children.end();
458 ++children_it) {
459 ContourLine& child = **children_it;
460 for (point = child.begin(); point != child.end(); ++point) {
461 *vertices_ptr++ = point->x;
462 *vertices_ptr++ = point->y;
463 *codes_ptr++ = (point == child.begin() ? MOVETO : LINETO);
464 }
465 point = child.begin();
466 *vertices_ptr++ = point->x;
467 *vertices_ptr++ = point->y;
468 *codes_ptr++ = CLOSEPOLY;
469
470 child.clear_parent(); // To indicate it can be deleted.
471 }
472
473 if (PyList_Append(vertices_list, vertices.pyobj_steal()) ||
474 PyList_Append(codes_list, codes.pyobj_steal())) {
475 Py_XDECREF(vertices_list);
476 Py_XDECREF(codes_list);
477 contour.delete_contour_lines();
478 throw std::runtime_error("Unable to add contour line to vertices and codes lists");
479 }
480
481 delete *line_it;
482 *line_it = 0;
483 }
484 }
485
486 // Delete remaining contour lines.
487 contour.delete_contour_lines();
488 }
489
calc_chunk_count(long point_count) const490 long QuadContourGenerator::calc_chunk_count(long point_count) const
491 {
492 assert(point_count > 0 && "point count must be positive");
493 assert(_chunk_size > 0 && "Chunk size must be positive");
494
495 if (_chunk_size > 0) {
496 long count = (point_count-1) / _chunk_size;
497 if (count*_chunk_size < point_count-1)
498 ++count;
499
500 assert(count >= 1 && "Invalid chunk count");
501 return count;
502 }
503 else
504 return 1;
505 }
506
create_contour(const double & level)507 PyObject* QuadContourGenerator::create_contour(const double& level)
508 {
509 init_cache_levels(level, level);
510
511 PyObject* vertices_list = PyList_New(0);
512 if (vertices_list == 0)
513 throw std::runtime_error("Failed to create Python list");
514
515 // Lines that start and end on boundaries.
516 long ichunk, jchunk, istart, iend, jstart, jend;
517 for (long ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) {
518 get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart, jend);
519
520 for (long j = jstart; j < jend; ++j) {
521 long quad_end = iend + j*_nx;
522 for (long quad = istart + j*_nx; quad < quad_end; ++quad) {
523 if (EXISTS_NONE(quad) || VISITED(quad,1)) continue;
524
525 if (BOUNDARY_S(quad) && Z_SW >= 1 && Z_SE < 1 &&
526 start_line(vertices_list, quad, Edge_S, level)) continue;
527
528 if (BOUNDARY_W(quad) && Z_NW >= 1 && Z_SW < 1 &&
529 start_line(vertices_list, quad, Edge_W, level)) continue;
530
531 if (BOUNDARY_N(quad) && Z_NE >= 1 && Z_NW < 1 &&
532 start_line(vertices_list, quad, Edge_N, level)) continue;
533
534 if (BOUNDARY_E(quad) && Z_SE >= 1 && Z_NE < 1 &&
535 start_line(vertices_list, quad, Edge_E, level)) continue;
536
537 if (_corner_mask) {
538 // Equates to NE boundary.
539 if (EXISTS_SW_CORNER(quad) && Z_SE >= 1 && Z_NW < 1 &&
540 start_line(vertices_list, quad, Edge_NE, level)) continue;
541
542 // Equates to NW boundary.
543 if (EXISTS_SE_CORNER(quad) && Z_NE >= 1 && Z_SW < 1 &&
544 start_line(vertices_list, quad, Edge_NW, level)) continue;
545
546 // Equates to SE boundary.
547 if (EXISTS_NW_CORNER(quad) && Z_SW >= 1 && Z_NE < 1 &&
548 start_line(vertices_list, quad, Edge_SE, level)) continue;
549
550 // Equates to SW boundary.
551 if (EXISTS_NE_CORNER(quad) && Z_NW >= 1 && Z_SE < 1 &&
552 start_line(vertices_list, quad, Edge_SW, level)) continue;
553 }
554 }
555 }
556 }
557
558 // Internal loops.
559 ContourLine contour_line(false); // Reused for each contour line.
560 for (long ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) {
561 get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart, jend);
562
563 for (long j = jstart; j < jend; ++j) {
564 long quad_end = iend + j*_nx;
565 for (long quad = istart + j*_nx; quad < quad_end; ++quad) {
566 if (EXISTS_NONE(quad) || VISITED(quad,1))
567 continue;
568
569 Edge start_edge = get_start_edge(quad, 1);
570 if (start_edge == Edge_None)
571 continue;
572
573 QuadEdge quad_edge(quad, start_edge);
574 QuadEdge start_quad_edge(quad_edge);
575
576 // To obtain output identical to that produced by legacy code,
577 // sometimes need to ignore the first point and add it on the
578 // end instead.
579 bool ignore_first = (start_edge == Edge_N);
580 follow_interior(contour_line, quad_edge, 1, level,
581 !ignore_first, &start_quad_edge, 1, false);
582 if (ignore_first && !contour_line.empty())
583 contour_line.push_back(contour_line.front());
584 append_contour_line_to_vertices(contour_line, vertices_list);
585
586 // Repeat if saddle point but not visited.
587 if (SADDLE(quad,1) && !VISITED(quad,1))
588 --quad;
589 }
590 }
591 }
592
593 return vertices_list;
594 }
595
create_filled_contour(const double & lower_level,const double & upper_level)596 PyObject* QuadContourGenerator::create_filled_contour(const double& lower_level,
597 const double& upper_level)
598 {
599 init_cache_levels(lower_level, upper_level);
600
601 Contour contour;
602
603 PyObject* vertices = PyList_New(0);
604 if (vertices == 0)
605 throw std::runtime_error("Failed to create Python list");
606
607 PyObject* codes = PyList_New(0);
608 if (codes == 0) {
609 Py_XDECREF(vertices);
610 throw std::runtime_error("Failed to create Python list");
611 }
612
613 long ichunk, jchunk, istart, iend, jstart, jend;
614 for (long ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) {
615 get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart, jend);
616 _parent_cache.set_chunk_starts(istart, jstart);
617
618 for (long j = jstart; j < jend; ++j) {
619 long quad_end = iend + j*_nx;
620 for (long quad = istart + j*_nx; quad < quad_end; ++quad) {
621 if (!EXISTS_NONE(quad))
622 single_quad_filled(contour, quad, lower_level, upper_level);
623 }
624 }
625
626 // Clear VISITED_W and VISITED_S flags that are reused by later chunks.
627 if (jchunk < _nychunk-1) {
628 long quad_end = iend + jend*_nx;
629 for (long quad = istart + jend*_nx; quad < quad_end; ++quad)
630 _cache[quad] &= ~MASK_VISITED_S;
631 }
632
633 if (ichunk < _nxchunk-1) {
634 long quad_end = iend + jend*_nx;
635 for (long quad = iend + jstart*_nx; quad < quad_end; quad += _nx)
636 _cache[quad] &= ~MASK_VISITED_W;
637 }
638
639 // Create python objects to return for this chunk.
640 append_contour_to_vertices_and_codes(contour, vertices, codes);
641 }
642
643 PyObject* tuple = PyTuple_New(2);
644 if (tuple == 0) {
645 Py_XDECREF(vertices);
646 Py_XDECREF(codes);
647 throw std::runtime_error("Failed to create Python tuple");
648 }
649
650 // No error checking here as filling in a brand new pre-allocated tuple.
651 PyTuple_SET_ITEM(tuple, 0, vertices);
652 PyTuple_SET_ITEM(tuple, 1, codes);
653
654 return tuple;
655 }
656
edge_interp(const QuadEdge & quad_edge,const double & level)657 XY QuadContourGenerator::edge_interp(const QuadEdge& quad_edge,
658 const double& level)
659 {
660 assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
661 "Quad index out of bounds");
662 assert(quad_edge.edge != Edge_None && "Invalid edge");
663 return interp(get_edge_point_index(quad_edge, true),
664 get_edge_point_index(quad_edge, false),
665 level);
666 }
667
follow_boundary(ContourLine & contour_line,QuadEdge & quad_edge,const double & lower_level,const double & upper_level,unsigned int level_index,const QuadEdge & start_quad_edge)668 unsigned int QuadContourGenerator::follow_boundary(
669 ContourLine& contour_line,
670 QuadEdge& quad_edge,
671 const double& lower_level,
672 const double& upper_level,
673 unsigned int level_index,
674 const QuadEdge& start_quad_edge)
675 {
676 assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
677 "Quad index out of bounds");
678 assert(quad_edge.edge != Edge_None && "Invalid edge");
679 assert(is_edge_a_boundary(quad_edge) && "Not a boundary edge");
680 assert((level_index == 1 || level_index == 2) &&
681 "level index must be 1 or 2");
682 assert(start_quad_edge.quad >= 0 && start_quad_edge.quad < _n &&
683 "Start quad index out of bounds");
684 assert(start_quad_edge.edge != Edge_None && "Invalid start edge");
685
686 // Only called for filled contours, so always updates _parent_cache.
687 unsigned int end_level = 0;
688 bool first_edge = true;
689 bool stop = false;
690 long& quad = quad_edge.quad;
691
692 while (true) {
693 // Levels of start and end points of quad_edge.
694 unsigned int start_level =
695 (first_edge ? Z_LEVEL(get_edge_point_index(quad_edge, true))
696 : end_level);
697 long end_point = get_edge_point_index(quad_edge, false);
698 end_level = Z_LEVEL(end_point);
699
700 if (level_index == 1) {
701 if (start_level <= level_index && end_level == 2) {
702 // Increasing z, switching levels from 1 to 2.
703 level_index = 2;
704 stop = true;
705 }
706 else if (start_level >= 1 && end_level == 0) {
707 // Decreasing z, keeping same level.
708 stop = true;
709 }
710 }
711 else { // level_index == 2
712 if (start_level <= level_index && end_level == 2) {
713 // Increasing z, keeping same level.
714 stop = true;
715 }
716 else if (start_level >= 1 && end_level == 0) {
717 // Decreasing z, switching levels from 2 to 1.
718 level_index = 1;
719 stop = true;
720 }
721 }
722
723 if (!first_edge && !stop && quad_edge == start_quad_edge)
724 // Return if reached start point of contour line. Do this before
725 // checking/setting VISITED flags as will already have been
726 // visited.
727 break;
728
729 switch (quad_edge.edge) {
730 case Edge_E:
731 assert(!VISITED_W(quad+1) && "Already visited");
732 _cache[quad+1] |= MASK_VISITED_W;
733 break;
734 case Edge_N:
735 assert(!VISITED_S(quad+_nx) && "Already visited");
736 _cache[quad+_nx] |= MASK_VISITED_S;
737 break;
738 case Edge_W:
739 assert(!VISITED_W(quad) && "Already visited");
740 _cache[quad] |= MASK_VISITED_W;
741 break;
742 case Edge_S:
743 assert(!VISITED_S(quad) && "Already visited");
744 _cache[quad] |= MASK_VISITED_S;
745 break;
746 case Edge_NE:
747 case Edge_NW:
748 case Edge_SW:
749 case Edge_SE:
750 assert(!VISITED_CORNER(quad) && "Already visited");
751 _cache[quad] |= MASK_VISITED_CORNER;
752 break;
753 default:
754 assert(0 && "Invalid Edge");
755 break;
756 }
757
758 if (stop) {
759 // Exiting boundary to enter interior.
760 contour_line.push_back(edge_interp(quad_edge,
761 level_index == 1 ? lower_level
762 : upper_level));
763 break;
764 }
765
766 move_to_next_boundary_edge(quad_edge);
767
768 // Just moved to new quad edge, so label parent of start of quad edge.
769 switch (quad_edge.edge) {
770 case Edge_W:
771 case Edge_SW:
772 case Edge_S:
773 case Edge_SE:
774 if (!EXISTS_SE_CORNER(quad))
775 _parent_cache.set_parent(quad, contour_line);
776 break;
777 case Edge_E:
778 case Edge_NE:
779 case Edge_N:
780 case Edge_NW:
781 if (!EXISTS_SW_CORNER(quad))
782 _parent_cache.set_parent(quad + 1, contour_line);
783 break;
784 default:
785 assert(0 && "Invalid edge");
786 break;
787 }
788
789 // Add point to contour.
790 contour_line.push_back(get_point_xy(end_point));
791
792 if (first_edge)
793 first_edge = false;
794 }
795
796 return level_index;
797 }
798
follow_interior(ContourLine & contour_line,QuadEdge & quad_edge,unsigned int level_index,const double & level,bool want_initial_point,const QuadEdge * start_quad_edge,unsigned int start_level_index,bool set_parents)799 void QuadContourGenerator::follow_interior(ContourLine& contour_line,
800 QuadEdge& quad_edge,
801 unsigned int level_index,
802 const double& level,
803 bool want_initial_point,
804 const QuadEdge* start_quad_edge,
805 unsigned int start_level_index,
806 bool set_parents)
807 {
808 assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
809 "Quad index out of bounds.");
810 assert(quad_edge.edge != Edge_None && "Invalid edge");
811 assert((level_index == 1 || level_index == 2) &&
812 "level index must be 1 or 2");
813 assert((start_quad_edge == 0 ||
814 (start_quad_edge->quad >= 0 && start_quad_edge->quad < _n)) &&
815 "Start quad index out of bounds.");
816 assert((start_quad_edge == 0 || start_quad_edge->edge != Edge_None) &&
817 "Invalid start edge");
818 assert((start_level_index == 1 || start_level_index == 2) &&
819 "start level index must be 1 or 2");
820
821 long& quad = quad_edge.quad;
822 Edge& edge = quad_edge.edge;
823
824 if (want_initial_point)
825 contour_line.push_back(edge_interp(quad_edge, level));
826
827 CacheItem visited_mask = (level_index == 1 ? MASK_VISITED_1 : MASK_VISITED_2);
828 CacheItem saddle_mask = (level_index == 1 ? MASK_SADDLE_1 : MASK_SADDLE_2);
829 Dir dir = Dir_Straight;
830
831 while (true) {
832 assert(!EXISTS_NONE(quad) && "Quad does not exist");
833 assert(!(_cache[quad] & visited_mask) && "Quad already visited");
834
835 // Determine direction to move to next quad. If the quad is already
836 // labelled as a saddle quad then the direction is easily read from
837 // the cache. Otherwise the direction is determined differently
838 // depending on whether the quad is a corner quad or not.
839
840 if (_cache[quad] & saddle_mask) {
841 // Already identified as a saddle quad, so direction is easy.
842 dir = (SADDLE_LEFT(quad,level_index) ? Dir_Left : Dir_Right);
843 _cache[quad] |= visited_mask;
844 }
845 else if (EXISTS_ANY_CORNER(quad)) {
846 // Need z-level of point opposite the entry edge, as that
847 // determines whether contour turns left or right.
848 long point_opposite = -1;
849 switch (edge) {
850 case Edge_E:
851 point_opposite = (EXISTS_SE_CORNER(quad) ? POINT_SW
852 : POINT_NW);
853 break;
854 case Edge_N:
855 point_opposite = (EXISTS_NW_CORNER(quad) ? POINT_SW
856 : POINT_SE);
857 break;
858 case Edge_W:
859 point_opposite = (EXISTS_SW_CORNER(quad) ? POINT_SE
860 : POINT_NE);
861 break;
862 case Edge_S:
863 point_opposite = (EXISTS_SW_CORNER(quad) ? POINT_NW
864 : POINT_NE);
865 break;
866 case Edge_NE: point_opposite = POINT_SW; break;
867 case Edge_NW: point_opposite = POINT_SE; break;
868 case Edge_SW: point_opposite = POINT_NE; break;
869 case Edge_SE: point_opposite = POINT_NW; break;
870 default: assert(0 && "Invalid edge"); break;
871 }
872 assert(point_opposite != -1 && "Failed to find opposite point");
873
874 // Lower-level polygons (level_index == 1) always have higher
875 // values to the left of the contour. Upper-level contours
876 // (level_index == 2) are reversed, which is what the fancy XOR
877 // does below.
878 if ((Z_LEVEL(point_opposite) >= level_index) ^ (level_index == 2))
879 dir = Dir_Right;
880 else
881 dir = Dir_Left;
882 _cache[quad] |= visited_mask;
883 }
884 else {
885 // Calculate configuration of this quad.
886 long point_left = -1, point_right = -1;
887 switch (edge) {
888 case Edge_E: point_left = POINT_SW; point_right = POINT_NW; break;
889 case Edge_N: point_left = POINT_SE; point_right = POINT_SW; break;
890 case Edge_W: point_left = POINT_NE; point_right = POINT_SE; break;
891 case Edge_S: point_left = POINT_NW; point_right = POINT_NE; break;
892 default: assert(0 && "Invalid edge"); break;
893 }
894
895 unsigned int config = (Z_LEVEL(point_left) >= level_index) << 1 |
896 (Z_LEVEL(point_right) >= level_index);
897
898 // Upper level (level_index == 2) polygons are reversed compared to
899 // lower level ones, i.e. higher values on the right rather than
900 // the left.
901 if (level_index == 2)
902 config = 3 - config;
903
904 // Calculate turn direction to move to next quad along contour line.
905 if (config == 1) {
906 // New saddle quad, set up cache bits for it.
907 double zmid = 0.25*(get_point_z(POINT_SW) +
908 get_point_z(POINT_SE) +
909 get_point_z(POINT_NW) +
910 get_point_z(POINT_NE));
911 _cache[quad] |= (level_index == 1 ? MASK_SADDLE_1 : MASK_SADDLE_2);
912 if ((zmid > level) ^ (level_index == 2)) {
913 dir = Dir_Right;
914 }
915 else {
916 dir = Dir_Left;
917 _cache[quad] |= (level_index == 1 ? MASK_SADDLE_LEFT_1
918 : MASK_SADDLE_LEFT_2);
919 }
920 if (edge == Edge_N || edge == Edge_E) {
921 // Next visit to this quad must start on S or W.
922 _cache[quad] |= (level_index == 1 ? MASK_SADDLE_START_SW_1
923 : MASK_SADDLE_START_SW_2);
924 }
925 }
926 else {
927 // Normal (non-saddle) quad.
928 dir = (config == 0 ? Dir_Left
929 : (config == 3 ? Dir_Right : Dir_Straight));
930 _cache[quad] |= visited_mask;
931 }
932 }
933
934 // Use dir to determine exit edge.
935 edge = get_exit_edge(quad_edge, dir);
936
937 if (set_parents) {
938 if (edge == Edge_E)
939 _parent_cache.set_parent(quad+1, contour_line);
940 else if (edge == Edge_W)
941 _parent_cache.set_parent(quad, contour_line);
942 }
943
944 // Add new point to contour line.
945 contour_line.push_back(edge_interp(quad_edge, level));
946
947 // Stop if reached boundary.
948 if (is_edge_a_boundary(quad_edge))
949 break;
950
951 move_to_next_quad(quad_edge);
952 assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
953 "Quad index out of bounds");
954
955 // Return if reached start point of contour line.
956 if (start_quad_edge != 0 &&
957 quad_edge == *start_quad_edge &&
958 level_index == start_level_index)
959 break;
960 }
961 }
962
get_chunk_limits(long ijchunk,long & ichunk,long & jchunk,long & istart,long & iend,long & jstart,long & jend)963 void QuadContourGenerator::get_chunk_limits(long ijchunk,
964 long& ichunk,
965 long& jchunk,
966 long& istart,
967 long& iend,
968 long& jstart,
969 long& jend)
970 {
971 assert(ijchunk >= 0 && ijchunk < _chunk_count && "ijchunk out of bounds");
972 ichunk = ijchunk % _nxchunk;
973 jchunk = ijchunk / _nxchunk;
974 istart = ichunk*_chunk_size;
975 iend = (ichunk == _nxchunk-1 ? _nx : (ichunk+1)*_chunk_size);
976 jstart = jchunk*_chunk_size;
977 jend = (jchunk == _nychunk-1 ? _ny : (jchunk+1)*_chunk_size);
978 }
979
get_corner_start_edge(long quad,unsigned int level_index) const980 Edge QuadContourGenerator::get_corner_start_edge(long quad,
981 unsigned int level_index) const
982 {
983 assert(quad >= 0 && quad < _n && "Quad index out of bounds");
984 assert((level_index == 1 || level_index == 2) &&
985 "level index must be 1 or 2");
986 assert(EXISTS_ANY_CORNER(quad) && "Quad is not a corner");
987
988 // Diagram for NE corner. Rotate for other corners.
989 //
990 // edge12
991 // point1 +---------+ point2
992 // \ |
993 // \ | edge23
994 // edge31 \ |
995 // \ |
996 // + point3
997 //
998 long point1, point2, point3;
999 Edge edge12, edge23, edge31;
1000 switch (_cache[quad] & MASK_EXISTS) {
1001 case MASK_EXISTS_SW_CORNER:
1002 point1 = POINT_SE; point2 = POINT_SW; point3 = POINT_NW;
1003 edge12 = Edge_S; edge23 = Edge_W; edge31 = Edge_NE;
1004 break;
1005 case MASK_EXISTS_SE_CORNER:
1006 point1 = POINT_NE; point2 = POINT_SE; point3 = POINT_SW;
1007 edge12 = Edge_E; edge23 = Edge_S; edge31 = Edge_NW;
1008 break;
1009 case MASK_EXISTS_NW_CORNER:
1010 point1 = POINT_SW; point2 = POINT_NW; point3 = POINT_NE;
1011 edge12 = Edge_W; edge23 = Edge_N; edge31 = Edge_SE;
1012 break;
1013 case MASK_EXISTS_NE_CORNER:
1014 point1 = POINT_NW; point2 = POINT_NE; point3 = POINT_SE;
1015 edge12 = Edge_N; edge23 = Edge_E; edge31 = Edge_SW;
1016 break;
1017 default:
1018 assert(0 && "Invalid EXISTS for quad");
1019 return Edge_None;
1020 }
1021
1022 unsigned int config = (Z_LEVEL(point1) >= level_index) << 2 |
1023 (Z_LEVEL(point2) >= level_index) << 1 |
1024 (Z_LEVEL(point3) >= level_index);
1025
1026 // Upper level (level_index == 2) polygons are reversed compared to lower
1027 // level ones, i.e. higher values on the right rather than the left.
1028 if (level_index == 2)
1029 config = 7 - config;
1030
1031 switch (config) {
1032 case 0: return Edge_None;
1033 case 1: return edge23;
1034 case 2: return edge12;
1035 case 3: return edge12;
1036 case 4: return edge31;
1037 case 5: return edge23;
1038 case 6: return edge31;
1039 case 7: return Edge_None;
1040 default: assert(0 && "Invalid config"); return Edge_None;
1041 }
1042 }
1043
get_edge_point_index(const QuadEdge & quad_edge,bool start) const1044 long QuadContourGenerator::get_edge_point_index(const QuadEdge& quad_edge,
1045 bool start) const
1046 {
1047 assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
1048 "Quad index out of bounds");
1049 assert(quad_edge.edge != Edge_None && "Invalid edge");
1050
1051 // Edges are ordered anticlockwise around their quad, as indicated by
1052 // directions of arrows in diagrams below.
1053 // Full quad NW corner (others similar)
1054 //
1055 // POINT_NW Edge_N POINT_NE POINT_NW Edge_N POINT_NE
1056 // +----<-----+ +----<-----+
1057 // | | | /
1058 // | | | quad /
1059 // Edge_W V quad ^ Edge_E Edge_W V ^
1060 // | | | / Edge_SE
1061 // | | | /
1062 // +---->-----+ +
1063 // POINT_SW Edge_S POINT_SE POINT_SW
1064 //
1065 const long& quad = quad_edge.quad;
1066 switch (quad_edge.edge) {
1067 case Edge_E: return (start ? POINT_SE : POINT_NE);
1068 case Edge_N: return (start ? POINT_NE : POINT_NW);
1069 case Edge_W: return (start ? POINT_NW : POINT_SW);
1070 case Edge_S: return (start ? POINT_SW : POINT_SE);
1071 case Edge_NE: return (start ? POINT_SE : POINT_NW);
1072 case Edge_NW: return (start ? POINT_NE : POINT_SW);
1073 case Edge_SW: return (start ? POINT_NW : POINT_SE);
1074 case Edge_SE: return (start ? POINT_SW : POINT_NE);
1075 default: assert(0 && "Invalid edge"); return 0;
1076 }
1077 }
1078
get_exit_edge(const QuadEdge & quad_edge,Dir dir) const1079 Edge QuadContourGenerator::get_exit_edge(const QuadEdge& quad_edge,
1080 Dir dir) const
1081 {
1082 assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
1083 "Quad index out of bounds");
1084 assert(quad_edge.edge != Edge_None && "Invalid edge");
1085
1086 const long& quad = quad_edge.quad;
1087 const Edge& edge = quad_edge.edge;
1088 if (EXISTS_ANY_CORNER(quad)) {
1089 // Corner directions are always left or right. A corner is a triangle,
1090 // entered via one edge so the other two edges are the left and right
1091 // ones.
1092 switch (edge) {
1093 case Edge_E:
1094 return (EXISTS_SE_CORNER(quad)
1095 ? (dir == Dir_Left ? Edge_S : Edge_NW)
1096 : (dir == Dir_Right ? Edge_N : Edge_SW));
1097 case Edge_N:
1098 return (EXISTS_NW_CORNER(quad)
1099 ? (dir == Dir_Right ? Edge_W : Edge_SE)
1100 : (dir == Dir_Left ? Edge_E : Edge_SW));
1101 case Edge_W:
1102 return (EXISTS_SW_CORNER(quad)
1103 ? (dir == Dir_Right ? Edge_S : Edge_NE)
1104 : (dir == Dir_Left ? Edge_N : Edge_SE));
1105 case Edge_S:
1106 return (EXISTS_SW_CORNER(quad)
1107 ? (dir == Dir_Left ? Edge_W : Edge_NE)
1108 : (dir == Dir_Right ? Edge_E : Edge_NW));
1109 case Edge_NE: return (dir == Dir_Left ? Edge_S : Edge_W);
1110 case Edge_NW: return (dir == Dir_Left ? Edge_E : Edge_S);
1111 case Edge_SW: return (dir == Dir_Left ? Edge_N : Edge_E);
1112 case Edge_SE: return (dir == Dir_Left ? Edge_W : Edge_N);
1113 default: assert(0 && "Invalid edge"); return Edge_None;
1114 }
1115 }
1116 else {
1117 // A full quad has four edges, entered via one edge so that other three
1118 // edges correspond to left, straight and right directions.
1119 switch (edge) {
1120 case Edge_E:
1121 return (dir == Dir_Left ? Edge_S :
1122 (dir == Dir_Right ? Edge_N : Edge_W));
1123 case Edge_N:
1124 return (dir == Dir_Left ? Edge_E :
1125 (dir == Dir_Right ? Edge_W : Edge_S));
1126 case Edge_W:
1127 return (dir == Dir_Left ? Edge_N :
1128 (dir == Dir_Right ? Edge_S : Edge_E));
1129 case Edge_S:
1130 return (dir == Dir_Left ? Edge_W :
1131 (dir == Dir_Right ? Edge_E : Edge_N));
1132 default: assert(0 && "Invalid edge"); return Edge_None;
1133 }
1134 }
1135 }
1136
get_point_xy(long point) const1137 XY QuadContourGenerator::get_point_xy(long point) const
1138 {
1139 assert(point >= 0 && point < _n && "Point index out of bounds.");
1140 return XY(_x.data()[static_cast<npy_intp>(point)],
1141 _y.data()[static_cast<npy_intp>(point)]);
1142 }
1143
get_point_z(long point) const1144 const double& QuadContourGenerator::get_point_z(long point) const
1145 {
1146 assert(point >= 0 && point < _n && "Point index out of bounds.");
1147 return _z.data()[static_cast<npy_intp>(point)];
1148 }
1149
get_quad_start_edge(long quad,unsigned int level_index) const1150 Edge QuadContourGenerator::get_quad_start_edge(long quad,
1151 unsigned int level_index) const
1152 {
1153 assert(quad >= 0 && quad < _n && "Quad index out of bounds");
1154 assert((level_index == 1 || level_index == 2) &&
1155 "level index must be 1 or 2");
1156 assert(EXISTS_QUAD(quad) && "Quad does not exist");
1157
1158 unsigned int config = (Z_NW >= level_index) << 3 |
1159 (Z_NE >= level_index) << 2 |
1160 (Z_SW >= level_index) << 1 |
1161 (Z_SE >= level_index);
1162
1163 // Upper level (level_index == 2) polygons are reversed compared to lower
1164 // level ones, i.e. higher values on the right rather than the left.
1165 if (level_index == 2)
1166 config = 15 - config;
1167
1168 switch (config) {
1169 case 0: return Edge_None;
1170 case 1: return Edge_E;
1171 case 2: return Edge_S;
1172 case 3: return Edge_E;
1173 case 4: return Edge_N;
1174 case 5: return Edge_N;
1175 case 6:
1176 // If already identified as a saddle quad then the start edge is
1177 // read from the cache. Otherwise return either valid start edge
1178 // and the subsequent call to follow_interior() will correctly set
1179 // up saddle bits in cache.
1180 if (!SADDLE(quad,level_index) || SADDLE_START_SW(quad,level_index))
1181 return Edge_S;
1182 else
1183 return Edge_N;
1184 case 7: return Edge_N;
1185 case 8: return Edge_W;
1186 case 9:
1187 // See comment for 6 above.
1188 if (!SADDLE(quad,level_index) || SADDLE_START_SW(quad,level_index))
1189 return Edge_W;
1190 else
1191 return Edge_E;
1192 case 10: return Edge_S;
1193 case 11: return Edge_E;
1194 case 12: return Edge_W;
1195 case 13: return Edge_W;
1196 case 14: return Edge_S;
1197 case 15: return Edge_None;
1198 default: assert(0 && "Invalid config"); return Edge_None;
1199 }
1200 }
1201
get_start_edge(long quad,unsigned int level_index) const1202 Edge QuadContourGenerator::get_start_edge(long quad,
1203 unsigned int level_index) const
1204 {
1205 if (EXISTS_ANY_CORNER(quad))
1206 return get_corner_start_edge(quad, level_index);
1207 else
1208 return get_quad_start_edge(quad, level_index);
1209 }
1210
init_cache_grid(const MaskArray & mask)1211 void QuadContourGenerator::init_cache_grid(const MaskArray& mask)
1212 {
1213 long i, j, quad;
1214
1215 if (mask.empty()) {
1216 // No mask, easy to calculate quad existence and boundaries together.
1217 quad = 0;
1218 for (j = 0; j < _ny; ++j) {
1219 for (i = 0; i < _nx; ++i, ++quad) {
1220 _cache[quad] = 0;
1221
1222 if (i < _nx-1 && j < _ny-1)
1223 _cache[quad] |= MASK_EXISTS_QUAD;
1224
1225 if ((i % _chunk_size == 0 || i == _nx-1) && j < _ny-1)
1226 _cache[quad] |= MASK_BOUNDARY_W;
1227
1228 if ((j % _chunk_size == 0 || j == _ny-1) && i < _nx-1)
1229 _cache[quad] |= MASK_BOUNDARY_S;
1230 }
1231 }
1232 }
1233 else {
1234 // Casting avoids problem when sizeof(bool) != sizeof(npy_bool).
1235 const npy_bool* mask_ptr =
1236 reinterpret_cast<const npy_bool*>(mask.data());
1237
1238 // Have mask so use two stages.
1239 // Stage 1, determine if quads/corners exist.
1240 quad = 0;
1241 for (j = 0; j < _ny; ++j) {
1242 for (i = 0; i < _nx; ++i, ++quad) {
1243 _cache[quad] = 0;
1244
1245 if (i < _nx-1 && j < _ny-1) {
1246 unsigned int config = mask_ptr[POINT_NW] << 3 |
1247 mask_ptr[POINT_NE] << 2 |
1248 mask_ptr[POINT_SW] << 1 |
1249 mask_ptr[POINT_SE];
1250
1251 if (_corner_mask) {
1252 switch (config) {
1253 case 0: _cache[quad] = MASK_EXISTS_QUAD; break;
1254 case 1: _cache[quad] = MASK_EXISTS_NW_CORNER; break;
1255 case 2: _cache[quad] = MASK_EXISTS_NE_CORNER; break;
1256 case 4: _cache[quad] = MASK_EXISTS_SW_CORNER; break;
1257 case 8: _cache[quad] = MASK_EXISTS_SE_CORNER; break;
1258 default:
1259 // Do nothing, quad is masked out.
1260 break;
1261 }
1262 }
1263 else if (config == 0)
1264 _cache[quad] = MASK_EXISTS_QUAD;
1265 }
1266 }
1267 }
1268
1269 // Stage 2, calculate W and S boundaries. For each quad use boundary
1270 // data already calculated for quads to W and S, so must iterate
1271 // through quads in correct order (increasing i and j indices).
1272 // Cannot use boundary data for quads to E and N as have not yet
1273 // calculated it.
1274 quad = 0;
1275 for (j = 0; j < _ny; ++j) {
1276 for (i = 0; i < _nx; ++i, ++quad) {
1277 if (_corner_mask) {
1278 bool W_exists_none = (i == 0 || EXISTS_NONE(quad-1));
1279 bool S_exists_none = (j == 0 || EXISTS_NONE(quad-_nx));
1280 bool W_exists_E_edge = (i > 0 && EXISTS_E_EDGE(quad-1));
1281 bool S_exists_N_edge = (j > 0 && EXISTS_N_EDGE(quad-_nx));
1282
1283 if ((EXISTS_W_EDGE(quad) && W_exists_none) ||
1284 (EXISTS_NONE(quad) && W_exists_E_edge) ||
1285 (i % _chunk_size == 0 && EXISTS_W_EDGE(quad) &&
1286 W_exists_E_edge))
1287 _cache[quad] |= MASK_BOUNDARY_W;
1288
1289 if ((EXISTS_S_EDGE(quad) && S_exists_none) ||
1290 (EXISTS_NONE(quad) && S_exists_N_edge) ||
1291 (j % _chunk_size == 0 && EXISTS_S_EDGE(quad) &&
1292 S_exists_N_edge))
1293 _cache[quad] |= MASK_BOUNDARY_S;
1294 }
1295 else {
1296 bool W_exists_quad = (i > 0 && EXISTS_QUAD(quad-1));
1297 bool S_exists_quad = (j > 0 && EXISTS_QUAD(quad-_nx));
1298
1299 if ((EXISTS_QUAD(quad) != W_exists_quad) ||
1300 (i % _chunk_size == 0 && EXISTS_QUAD(quad) &&
1301 W_exists_quad))
1302 _cache[quad] |= MASK_BOUNDARY_W;
1303
1304 if ((EXISTS_QUAD(quad) != S_exists_quad) ||
1305 (j % _chunk_size == 0 && EXISTS_QUAD(quad) &&
1306 S_exists_quad))
1307 _cache[quad] |= MASK_BOUNDARY_S;
1308 }
1309 }
1310 }
1311 }
1312 }
1313
init_cache_levels(const double & lower_level,const double & upper_level)1314 void QuadContourGenerator::init_cache_levels(const double& lower_level,
1315 const double& upper_level)
1316 {
1317 assert(upper_level >= lower_level &&
1318 "upper and lower levels are wrong way round");
1319
1320 bool two_levels = (lower_level != upper_level);
1321 CacheItem keep_mask =
1322 (_corner_mask ? MASK_EXISTS | MASK_BOUNDARY_S | MASK_BOUNDARY_W
1323 : MASK_EXISTS_QUAD | MASK_BOUNDARY_S | MASK_BOUNDARY_W);
1324
1325 if (two_levels) {
1326 const double* z_ptr = _z.data();
1327 for (long quad = 0; quad < _n; ++quad, ++z_ptr) {
1328 _cache[quad] &= keep_mask;
1329 if (*z_ptr > upper_level)
1330 _cache[quad] |= MASK_Z_LEVEL_2;
1331 else if (*z_ptr > lower_level)
1332 _cache[quad] |= MASK_Z_LEVEL_1;
1333 }
1334 }
1335 else {
1336 const double* z_ptr = _z.data();
1337 for (long quad = 0; quad < _n; ++quad, ++z_ptr) {
1338 _cache[quad] &= keep_mask;
1339 if (*z_ptr > lower_level)
1340 _cache[quad] |= MASK_Z_LEVEL_1;
1341 }
1342 }
1343 }
1344
interp(long point1,long point2,const double & level) const1345 XY QuadContourGenerator::interp(
1346 long point1, long point2, const double& level) const
1347 {
1348 assert(point1 >= 0 && point1 < _n && "Point index 1 out of bounds.");
1349 assert(point2 >= 0 && point2 < _n && "Point index 2 out of bounds.");
1350 assert(point1 != point2 && "Identical points");
1351 double fraction = (get_point_z(point2) - level) /
1352 (get_point_z(point2) - get_point_z(point1));
1353 return get_point_xy(point1)*fraction + get_point_xy(point2)*(1.0 - fraction);
1354 }
1355
is_edge_a_boundary(const QuadEdge & quad_edge) const1356 bool QuadContourGenerator::is_edge_a_boundary(const QuadEdge& quad_edge) const
1357 {
1358 assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
1359 "Quad index out of bounds");
1360 assert(quad_edge.edge != Edge_None && "Invalid edge");
1361
1362 switch (quad_edge.edge) {
1363 case Edge_E: return BOUNDARY_E(quad_edge.quad);
1364 case Edge_N: return BOUNDARY_N(quad_edge.quad);
1365 case Edge_W: return BOUNDARY_W(quad_edge.quad);
1366 case Edge_S: return BOUNDARY_S(quad_edge.quad);
1367 case Edge_NE: return EXISTS_SW_CORNER(quad_edge.quad);
1368 case Edge_NW: return EXISTS_SE_CORNER(quad_edge.quad);
1369 case Edge_SW: return EXISTS_NE_CORNER(quad_edge.quad);
1370 case Edge_SE: return EXISTS_NW_CORNER(quad_edge.quad);
1371 default: assert(0 && "Invalid edge"); return true;
1372 }
1373 }
1374
move_to_next_boundary_edge(QuadEdge & quad_edge) const1375 void QuadContourGenerator::move_to_next_boundary_edge(QuadEdge& quad_edge) const
1376 {
1377 assert(is_edge_a_boundary(quad_edge) && "QuadEdge is not a boundary");
1378
1379 long& quad = quad_edge.quad;
1380 Edge& edge = quad_edge.edge;
1381
1382 quad = get_edge_point_index(quad_edge, false);
1383
1384 // quad is now such that POINT_SW is the end point of the quad_edge passed
1385 // to this function.
1386
1387 // To find the next boundary edge, first attempt to turn left 135 degrees
1388 // and if that edge is a boundary then move to it. If not, attempt to turn
1389 // left 90 degrees, then left 45 degrees, then straight on, etc, until can
1390 // move.
1391 // First determine which edge to attempt first.
1392 int index = 0;
1393 switch (edge) {
1394 case Edge_E: index = 0; break;
1395 case Edge_SE: index = 1; break;
1396 case Edge_S: index = 2; break;
1397 case Edge_SW: index = 3; break;
1398 case Edge_W: index = 4; break;
1399 case Edge_NW: index = 5; break;
1400 case Edge_N: index = 6; break;
1401 case Edge_NE: index = 7; break;
1402 default: assert(0 && "Invalid edge"); break;
1403 }
1404
1405 // If _corner_mask not set, only need to consider odd index in loop below.
1406 if (!_corner_mask)
1407 ++index;
1408
1409 // Try each edge in turn until a boundary is found.
1410 int start_index = index;
1411 do
1412 {
1413 switch (index) {
1414 case 0:
1415 if (EXISTS_SE_CORNER(quad-_nx-1)) { // Equivalent to BOUNDARY_NW
1416 quad -= _nx+1;
1417 edge = Edge_NW;
1418 return;
1419 }
1420 break;
1421 case 1:
1422 if (BOUNDARY_N(quad-_nx-1)) {
1423 quad -= _nx+1;
1424 edge = Edge_N;
1425 return;
1426 }
1427 break;
1428 case 2:
1429 if (EXISTS_SW_CORNER(quad-1)) { // Equivalent to BOUNDARY_NE
1430 quad -= 1;
1431 edge = Edge_NE;
1432 return;
1433 }
1434 break;
1435 case 3:
1436 if (BOUNDARY_E(quad-1)) {
1437 quad -= 1;
1438 edge = Edge_E;
1439 return;
1440 }
1441 break;
1442 case 4:
1443 if (EXISTS_NW_CORNER(quad)) { // Equivalent to BOUNDARY_SE
1444 edge = Edge_SE;
1445 return;
1446 }
1447 break;
1448 case 5:
1449 if (BOUNDARY_S(quad)) {
1450 edge = Edge_S;
1451 return;
1452 }
1453 break;
1454 case 6:
1455 if (EXISTS_NE_CORNER(quad-_nx)) { // Equivalent to BOUNDARY_SW
1456 quad -= _nx;
1457 edge = Edge_SW;
1458 return;
1459 }
1460 break;
1461 case 7:
1462 if (BOUNDARY_W(quad-_nx)) {
1463 quad -= _nx;
1464 edge = Edge_W;
1465 return;
1466 }
1467 break;
1468 default: assert(0 && "Invalid index"); break;
1469 }
1470
1471 if (_corner_mask)
1472 index = (index + 1) % 8;
1473 else
1474 index = (index + 2) % 8;
1475 } while (index != start_index);
1476
1477 assert(0 && "Failed to find next boundary edge");
1478 }
1479
move_to_next_quad(QuadEdge & quad_edge) const1480 void QuadContourGenerator::move_to_next_quad(QuadEdge& quad_edge) const
1481 {
1482 assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
1483 "Quad index out of bounds");
1484 assert(quad_edge.edge != Edge_None && "Invalid edge");
1485
1486 // Move from quad_edge.quad to the neighbouring quad in the direction
1487 // specified by quad_edge.edge.
1488 switch (quad_edge.edge) {
1489 case Edge_E: quad_edge.quad += 1; quad_edge.edge = Edge_W; break;
1490 case Edge_N: quad_edge.quad += _nx; quad_edge.edge = Edge_S; break;
1491 case Edge_W: quad_edge.quad -= 1; quad_edge.edge = Edge_E; break;
1492 case Edge_S: quad_edge.quad -= _nx; quad_edge.edge = Edge_N; break;
1493 default: assert(0 && "Invalid edge"); break;
1494 }
1495 }
1496
single_quad_filled(Contour & contour,long quad,const double & lower_level,const double & upper_level)1497 void QuadContourGenerator::single_quad_filled(Contour& contour,
1498 long quad,
1499 const double& lower_level,
1500 const double& upper_level)
1501 {
1502 assert(quad >= 0 && quad < _n && "Quad index out of bounds");
1503
1504 // Order of checking is important here as can have different ContourLines
1505 // from both lower and upper levels in the same quad. First check the S
1506 // edge, then move up the quad to the N edge checking as required.
1507
1508 // Possible starts from S boundary.
1509 if (BOUNDARY_S(quad) && EXISTS_S_EDGE(quad)) {
1510
1511 // Lower-level start from S boundary into interior.
1512 if (!VISITED_S(quad) && Z_SW >= 1 && Z_SE == 0)
1513 contour.push_back(start_filled(quad, Edge_S, 1, NotHole, Interior,
1514 lower_level, upper_level));
1515
1516 // Upper-level start from S boundary into interior.
1517 if (!VISITED_S(quad) && Z_SW < 2 && Z_SE == 2)
1518 contour.push_back(start_filled(quad, Edge_S, 2, NotHole, Interior,
1519 lower_level, upper_level));
1520
1521 // Lower-level start following S boundary from W to E.
1522 if (!VISITED_S(quad) && Z_SW <= 1 && Z_SE == 1)
1523 contour.push_back(start_filled(quad, Edge_S, 1, NotHole, Boundary,
1524 lower_level, upper_level));
1525
1526 // Upper-level start following S boundary from W to E.
1527 if (!VISITED_S(quad) && Z_SW == 2 && Z_SE == 1)
1528 contour.push_back(start_filled(quad, Edge_S, 2, NotHole, Boundary,
1529 lower_level, upper_level));
1530 }
1531
1532 // Possible starts from W boundary.
1533 if (BOUNDARY_W(quad) && EXISTS_W_EDGE(quad)) {
1534
1535 // Lower-level start from W boundary into interior.
1536 if (!VISITED_W(quad) && Z_NW >= 1 && Z_SW == 0)
1537 contour.push_back(start_filled(quad, Edge_W, 1, NotHole, Interior,
1538 lower_level, upper_level));
1539
1540 // Upper-level start from W boundary into interior.
1541 if (!VISITED_W(quad) && Z_NW < 2 && Z_SW == 2)
1542 contour.push_back(start_filled(quad, Edge_W, 2, NotHole, Interior,
1543 lower_level, upper_level));
1544
1545 // Lower-level start following W boundary from N to S.
1546 if (!VISITED_W(quad) && Z_NW <= 1 && Z_SW == 1)
1547 contour.push_back(start_filled(quad, Edge_W, 1, NotHole, Boundary,
1548 lower_level, upper_level));
1549
1550 // Upper-level start following W boundary from N to S.
1551 if (!VISITED_W(quad) && Z_NW == 2 && Z_SW == 1)
1552 contour.push_back(start_filled(quad, Edge_W, 2, NotHole, Boundary,
1553 lower_level, upper_level));
1554 }
1555
1556 // Possible starts from NE boundary.
1557 if (EXISTS_SW_CORNER(quad)) { // i.e. BOUNDARY_NE
1558
1559 // Lower-level start following NE boundary from SE to NW, hole.
1560 if (!VISITED_CORNER(quad) && Z_NW == 1 && Z_SE == 1)
1561 contour.push_back(start_filled(quad, Edge_NE, 1, Hole, Boundary,
1562 lower_level, upper_level));
1563 }
1564 // Possible starts from SE boundary.
1565 else if (EXISTS_NW_CORNER(quad)) { // i.e. BOUNDARY_SE
1566
1567 // Lower-level start from N to SE.
1568 if (!VISITED(quad,1) && Z_NW == 0 && Z_SW == 0 && Z_NE >= 1)
1569 contour.push_back(start_filled(quad, Edge_N, 1, NotHole, Interior,
1570 lower_level, upper_level));
1571
1572 // Upper-level start from SE to N, hole.
1573 if (!VISITED(quad,2) && Z_NW < 2 && Z_SW < 2 && Z_NE == 2)
1574 contour.push_back(start_filled(quad, Edge_SE, 2, Hole, Interior,
1575 lower_level, upper_level));
1576
1577 // Upper-level start from N to SE.
1578 if (!VISITED(quad,2) && Z_NW == 2 && Z_SW == 2 && Z_NE < 2)
1579 contour.push_back(start_filled(quad, Edge_N, 2, NotHole, Interior,
1580 lower_level, upper_level));
1581
1582 // Lower-level start from SE to N, hole.
1583 if (!VISITED(quad,1) && Z_NW >= 1 && Z_SW >= 1 && Z_NE == 0)
1584 contour.push_back(start_filled(quad, Edge_SE, 1, Hole, Interior,
1585 lower_level, upper_level));
1586 }
1587 // Possible starts from NW boundary.
1588 else if (EXISTS_SE_CORNER(quad)) { // i.e. BOUNDARY_NW
1589
1590 // Lower-level start from NW to E.
1591 if (!VISITED(quad,1) && Z_SW == 0 && Z_SE == 0 && Z_NE >= 1)
1592 contour.push_back(start_filled(quad, Edge_NW, 1, NotHole, Interior,
1593 lower_level, upper_level));
1594
1595 // Upper-level start from E to NW, hole.
1596 if (!VISITED(quad,2) && Z_SW < 2 && Z_SE < 2 && Z_NE == 2)
1597 contour.push_back(start_filled(quad, Edge_E, 2, Hole, Interior,
1598 lower_level, upper_level));
1599
1600 // Upper-level start from NW to E.
1601 if (!VISITED(quad,2) && Z_SW == 2 && Z_SE == 2 && Z_NE < 2)
1602 contour.push_back(start_filled(quad, Edge_NW, 2, NotHole, Interior,
1603 lower_level, upper_level));
1604
1605 // Lower-level start from E to NW, hole.
1606 if (!VISITED(quad,1) && Z_SW >= 1 && Z_SE >= 1 && Z_NE == 0)
1607 contour.push_back(start_filled(quad, Edge_E, 1, Hole, Interior,
1608 lower_level, upper_level));
1609 }
1610 // Possible starts from SW boundary.
1611 else if (EXISTS_NE_CORNER(quad)) { // i.e. BOUNDARY_SW
1612
1613 // Lower-level start from SW boundary into interior.
1614 if (!VISITED_CORNER(quad) && Z_NW >= 1 && Z_SE == 0)
1615 contour.push_back(start_filled(quad, Edge_SW, 1, NotHole, Interior,
1616 lower_level, upper_level));
1617
1618 // Upper-level start from SW boundary into interior.
1619 if (!VISITED_CORNER(quad) && Z_NW < 2 && Z_SE == 2)
1620 contour.push_back(start_filled(quad, Edge_SW, 2, NotHole, Interior,
1621 lower_level, upper_level));
1622
1623 // Lower-level start following SW boundary from NW to SE.
1624 if (!VISITED_CORNER(quad) && Z_NW <= 1 && Z_SE == 1)
1625 contour.push_back(start_filled(quad, Edge_SW, 1, NotHole, Boundary,
1626 lower_level, upper_level));
1627
1628 // Upper-level start following SW boundary from NW to SE.
1629 if (!VISITED_CORNER(quad) && Z_NW == 2 && Z_SE == 1)
1630 contour.push_back(start_filled(quad, Edge_SW, 2, NotHole, Boundary,
1631 lower_level, upper_level));
1632 }
1633
1634 // A full (unmasked) quad can only have a start on the NE corner, i.e. from
1635 // N to E (lower level) or E to N (upper level). Any other start will have
1636 // already been created in a call to this function for a prior quad so we
1637 // don't need to test for it again here.
1638 //
1639 // The situation is complicated by the possibility that the quad is a
1640 // saddle quad, in which case a contour line starting on the N could leave
1641 // by either the W or the E. We only need to consider those leaving E.
1642 //
1643 // A NE corner can also have a N to E or E to N start.
1644 if (EXISTS_QUAD(quad) || EXISTS_NE_CORNER(quad)) {
1645
1646 // Lower-level start from N to E.
1647 if (!VISITED(quad,1) && Z_NW == 0 && Z_SE == 0 && Z_NE >= 1 &&
1648 (!SADDLE(quad,1) || SADDLE_LEFT(quad,1)))
1649 contour.push_back(start_filled(quad, Edge_N, 1, NotHole, Interior,
1650 lower_level, upper_level));
1651
1652 // Upper-level start from E to N, hole.
1653 if (!VISITED(quad,2) && Z_NW < 2 && Z_SE < 2 && Z_NE == 2 &&
1654 (!SADDLE(quad,2) || !SADDLE_LEFT(quad,2)))
1655 contour.push_back(start_filled(quad, Edge_E, 2, Hole, Interior,
1656 lower_level, upper_level));
1657
1658 // Upper-level start from N to E.
1659 if (!VISITED(quad,2) && Z_NW == 2 && Z_SE == 2 && Z_NE < 2 &&
1660 (!SADDLE(quad,2) || SADDLE_LEFT(quad,2)))
1661 contour.push_back(start_filled(quad, Edge_N, 2, NotHole, Interior,
1662 lower_level, upper_level));
1663
1664 // Lower-level start from E to N, hole.
1665 if (!VISITED(quad,1) && Z_NW >= 1 && Z_SE >= 1 && Z_NE == 0 &&
1666 (!SADDLE(quad,1) || !SADDLE_LEFT(quad,1)))
1667 contour.push_back(start_filled(quad, Edge_E, 1, Hole, Interior,
1668 lower_level, upper_level));
1669
1670 // All possible contours passing through the interior of this quad
1671 // should have already been created, so assert this.
1672 assert((VISITED(quad,1) || get_start_edge(quad, 1) == Edge_None) &&
1673 "Found start of contour that should have already been created");
1674 assert((VISITED(quad,2) || get_start_edge(quad, 2) == Edge_None) &&
1675 "Found start of contour that should have already been created");
1676 }
1677
1678 // Lower-level start following N boundary from E to W, hole.
1679 // This is required for an internal masked region which is a hole in a
1680 // surrounding contour line.
1681 if (BOUNDARY_N(quad) && EXISTS_N_EDGE(quad) &&
1682 !VISITED_S(quad+_nx) && Z_NW == 1 && Z_NE == 1)
1683 contour.push_back(start_filled(quad, Edge_N, 1, Hole, Boundary,
1684 lower_level, upper_level));
1685 }
1686
start_filled(long quad,Edge edge,unsigned int start_level_index,HoleOrNot hole_or_not,BoundaryOrInterior boundary_or_interior,const double & lower_level,const double & upper_level)1687 ContourLine* QuadContourGenerator::start_filled(
1688 long quad,
1689 Edge edge,
1690 unsigned int start_level_index,
1691 HoleOrNot hole_or_not,
1692 BoundaryOrInterior boundary_or_interior,
1693 const double& lower_level,
1694 const double& upper_level)
1695 {
1696 assert(quad >= 0 && quad < _n && "Quad index out of bounds");
1697 assert(edge != Edge_None && "Invalid edge");
1698 assert((start_level_index == 1 || start_level_index == 2) &&
1699 "start level index must be 1 or 2");
1700
1701 ContourLine* contour_line = new ContourLine(hole_or_not == Hole);
1702 if (hole_or_not == Hole) {
1703 // Find and set parent ContourLine.
1704 ContourLine* parent = _parent_cache.get_parent(quad + 1);
1705 assert(parent != 0 && "Failed to find parent ContourLine");
1706 contour_line->set_parent(parent);
1707 parent->add_child(contour_line);
1708 }
1709
1710 QuadEdge quad_edge(quad, edge);
1711 const QuadEdge start_quad_edge(quad_edge);
1712 unsigned int level_index = start_level_index;
1713
1714 // If starts on interior, can only finish on interior.
1715 // If starts on boundary, can only finish on boundary.
1716
1717 while (true) {
1718 if (boundary_or_interior == Interior) {
1719 double level = (level_index == 1 ? lower_level : upper_level);
1720 follow_interior(*contour_line, quad_edge, level_index, level,
1721 false, &start_quad_edge, start_level_index, true);
1722 }
1723 else {
1724 level_index = follow_boundary(
1725 *contour_line, quad_edge, lower_level,
1726 upper_level, level_index, start_quad_edge);
1727 }
1728
1729 if (quad_edge == start_quad_edge && (boundary_or_interior == Boundary ||
1730 level_index == start_level_index))
1731 break;
1732
1733 if (boundary_or_interior == Boundary)
1734 boundary_or_interior = Interior;
1735 else
1736 boundary_or_interior = Boundary;
1737 }
1738
1739 return contour_line;
1740 }
1741
start_line(PyObject * vertices_list,long quad,Edge edge,const double & level)1742 bool QuadContourGenerator::start_line(
1743 PyObject* vertices_list, long quad, Edge edge, const double& level)
1744 {
1745 assert(vertices_list != 0 && "Null python vertices list");
1746 assert(is_edge_a_boundary(QuadEdge(quad, edge)) &&
1747 "QuadEdge is not a boundary");
1748
1749 QuadEdge quad_edge(quad, edge);
1750 ContourLine contour_line(false);
1751 follow_interior(contour_line, quad_edge, 1, level, true, 0, 1, false);
1752 append_contour_line_to_vertices(contour_line, vertices_list);
1753 return VISITED(quad,1);
1754 }
1755
write_cache(bool grid_only) const1756 void QuadContourGenerator::write_cache(bool grid_only) const
1757 {
1758 std::cout << "-----------------------------------------------" << std::endl;
1759 for (long quad = 0; quad < _n; ++quad)
1760 write_cache_quad(quad, grid_only);
1761 std::cout << "-----------------------------------------------" << std::endl;
1762 }
1763
write_cache_quad(long quad,bool grid_only) const1764 void QuadContourGenerator::write_cache_quad(long quad, bool grid_only) const
1765 {
1766 long j = quad / _nx;
1767 long i = quad - j*_nx;
1768 std::cout << quad << ": i=" << i << " j=" << j
1769 << " EXISTS=" << EXISTS_QUAD(quad);
1770 if (_corner_mask)
1771 std::cout << " CORNER=" << EXISTS_SW_CORNER(quad) << EXISTS_SE_CORNER(quad)
1772 << EXISTS_NW_CORNER(quad) << EXISTS_NE_CORNER(quad);
1773 std::cout << " BNDY=" << (BOUNDARY_S(quad)>0) << (BOUNDARY_W(quad)>0);
1774 if (!grid_only) {
1775 std::cout << " Z=" << Z_LEVEL(quad)
1776 << " SAD=" << (SADDLE(quad,1)>0) << (SADDLE(quad,2)>0)
1777 << " LEFT=" << (SADDLE_LEFT(quad,1)>0) << (SADDLE_LEFT(quad,2)>0)
1778 << " NW=" << (SADDLE_START_SW(quad,1)>0) << (SADDLE_START_SW(quad,2)>0)
1779 << " VIS=" << (VISITED(quad,1)>0) << (VISITED(quad,2)>0)
1780 << (VISITED_S(quad)>0) << (VISITED_W(quad)>0)
1781 << (VISITED_CORNER(quad)>0);
1782 }
1783 std::cout << std::endl;
1784 }
1785