forked from alandefreitas/matplotplusplus
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathcontourc.cpp
More file actions
1971 lines (1751 loc) 路 79.4 KB
/
Copy pathcontourc.cpp
File metadata and controls
1971 lines (1751 loc) 路 79.4 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
//
// Adapted from
// https://github.com/matplotlib/matplotlib/blob/master/src/_contour.cpp
//
// This file contains liberal use of asserts to assist code development and
// debugging. Standard matplotlib builds disable asserts so they cause no
// performance reduction. To enable the asserts, you need to undefine the
// NDEBUG macro, which is achieved by adding the following
// undef_macros=['NDEBUG']
// to the appropriate make_extension call in setupext.py, and then rebuilding.
#define NO_IMPORT_ARRAY
// #include "mplutils.h"
// https://github.com/matplotlib/matplotlib/blob/master/src/mplutils.h
// https://github.com/matplotlib/matplotlib/blob/master/src/mplutils.cpp
#include <algorithm>
#include <cassert>
#include <matplot/axes_objects/contours.h>
#include <matplot/util/contourc.h>
namespace matplot {
// 'kind' codes.
#define MOVETO 1
#define LINETO 2
#define CLOSEPOLY 79
// Point indices from current quad index.
#define POINT_SW (quad)
#define POINT_SE (quad + 1)
#define POINT_NW (quad + _nx)
#define POINT_NE (quad + _nx + 1)
// CacheItem masks, only accessed directly to set. To read, use accessors
// detailed below. 1 and 2 refer to level indices (lower and upper).
#define MASK_Z_LEVEL 0x0003 // Combines the following two.
#define MASK_Z_LEVEL_1 0x0001 // z > lower_level.
#define MASK_Z_LEVEL_2 0x0002 // z > upper_level.
#define MASK_VISITED_1 0x0004 // Algorithm has visited this quad.
#define MASK_VISITED_2 0x0008
#define MASK_SADDLE_1 0x0010 // quad is a saddle quad.
#define MASK_SADDLE_2 0x0020
#define MASK_SADDLE_LEFT_1 0x0040 // Contours turn left at saddle quad.
#define MASK_SADDLE_LEFT_2 0x0080
#define MASK_SADDLE_START_SW_1 0x0100 // Next visit starts on S or W edge.
#define MASK_SADDLE_START_SW_2 0x0200
#define MASK_BOUNDARY_S 0x0400 // S edge of quad is a boundary.
#define MASK_BOUNDARY_W 0x0800 // W edge of quad is a boundary.
// EXISTS_QUAD bit is always used, but the 4 EXISTS_CORNER are only used if
// _corner_mask is true. Only one of EXISTS_QUAD or EXISTS_??_CORNER is ever
// set per quad, hence not using unique bits for each; care is needed when
// testing for these flags as they overlap.
#define MASK_EXISTS_QUAD 0x1000 // All of quad exists (is not masked).
#define MASK_EXISTS_SW_CORNER 0x2000 // SW corner exists, NE corner is masked.
#define MASK_EXISTS_SE_CORNER 0x3000
#define MASK_EXISTS_NW_CORNER 0x4000
#define MASK_EXISTS_NE_CORNER 0x5000
#define MASK_EXISTS 0x7000 // Combines all 5 EXISTS masks.
// The following are only needed for filled contours.
#define MASK_VISITED_S 0x10000 // Algorithm has visited S boundary.
#define MASK_VISITED_W 0x20000 // Algorithm has visited W boundary.
#define MASK_VISITED_CORNER 0x40000 // Algorithm has visited corner edge.
// Accessors for various CacheItem masks. li is shorthand for level_index.
#define Z_LEVEL(quad) (_cache[quad] & MASK_Z_LEVEL)
#define Z_NE Z_LEVEL(POINT_NE)
#define Z_NW Z_LEVEL(POINT_NW)
#define Z_SE Z_LEVEL(POINT_SE)
#define Z_SW Z_LEVEL(POINT_SW)
#define VISITED(quad, li) \
((_cache[quad] & (li == 1 ? MASK_VISITED_1 : MASK_VISITED_2)) != 0)
#define VISITED_S(quad) ((_cache[quad] & MASK_VISITED_S) != 0)
#define VISITED_W(quad) ((_cache[quad] & MASK_VISITED_W) != 0)
#define VISITED_CORNER(quad) ((_cache[quad] & MASK_VISITED_CORNER) != 0)
#define SADDLE(quad, li) \
((_cache[quad] & (li == 1 ? MASK_SADDLE_1 : MASK_SADDLE_2)) != 0)
#define SADDLE_LEFT(quad, li) \
((_cache[quad] & (li == 1 ? MASK_SADDLE_LEFT_1 : MASK_SADDLE_LEFT_2)) != 0)
#define SADDLE_START_SW(quad, li) \
((_cache[quad] & \
(li == 1 ? MASK_SADDLE_START_SW_1 : MASK_SADDLE_START_SW_2)) != 0)
#define BOUNDARY_S(quad) ((_cache[quad] & MASK_BOUNDARY_S) != 0)
#define BOUNDARY_W(quad) ((_cache[quad] & MASK_BOUNDARY_W) != 0)
#define BOUNDARY_N(quad) BOUNDARY_S(quad + _nx)
#define BOUNDARY_E(quad) BOUNDARY_W(quad + 1)
#define EXISTS_QUAD(quad) ((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_QUAD)
#define EXISTS_NONE(quad) ((_cache[quad] & MASK_EXISTS) == 0)
// The following are only used if _corner_mask is true.
#define EXISTS_SW_CORNER(quad) \
((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_SW_CORNER)
#define EXISTS_SE_CORNER(quad) \
((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_SE_CORNER)
#define EXISTS_NW_CORNER(quad) \
((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_NW_CORNER)
#define EXISTS_NE_CORNER(quad) \
((_cache[quad] & MASK_EXISTS) == MASK_EXISTS_NE_CORNER)
#define EXISTS_ANY_CORNER(quad) (!EXISTS_NONE(quad) && !EXISTS_QUAD(quad))
#define EXISTS_W_EDGE(quad) \
(EXISTS_QUAD(quad) || EXISTS_SW_CORNER(quad) || EXISTS_NW_CORNER(quad))
#define EXISTS_E_EDGE(quad) \
(EXISTS_QUAD(quad) || EXISTS_SE_CORNER(quad) || EXISTS_NE_CORNER(quad))
#define EXISTS_S_EDGE(quad) \
(EXISTS_QUAD(quad) || EXISTS_SW_CORNER(quad) || EXISTS_SE_CORNER(quad))
#define EXISTS_N_EDGE(quad) \
(EXISTS_QUAD(quad) || EXISTS_NW_CORNER(quad) || EXISTS_NE_CORNER(quad))
// Note that EXISTS_NE_CORNER(quad) is equivalent to BOUNDARY_SW(quad), etc.
QuadEdge::QuadEdge() : quad(-1), edge(Edge_None) {}
QuadEdge::QuadEdge(long quad_, Edge edge_) : quad(quad_), edge(edge_) {}
bool QuadEdge::operator<(const QuadEdge &other) const {
if (quad != other.quad)
return quad < other.quad;
else
return edge < other.edge;
}
bool QuadEdge::operator==(const QuadEdge &other) const {
return quad == other.quad && edge == other.edge;
}
bool QuadEdge::operator!=(const QuadEdge &other) const {
return !operator==(other);
}
std::ostream &operator<<(std::ostream &os, const QuadEdge &quad_edge) {
return os << quad_edge.quad << ' ' << quad_edge.edge;
}
XY::XY() {}
XY::XY(const double &x_, const double &y_) : x(x_), y(y_) {}
bool XY::operator==(const XY &other) const {
return x == other.x && y == other.y;
}
bool XY::operator!=(const XY &other) const {
return x != other.x || y != other.y;
}
XY XY::operator*(const double &multiplier) const {
return XY(x * multiplier, y * multiplier);
}
const XY &XY::operator+=(const XY &other) {
x += other.x;
y += other.y;
return *this;
}
const XY &XY::operator-=(const XY &other) {
x -= other.x;
y -= other.y;
return *this;
}
XY XY::operator+(const XY &other) const {
return XY(x + other.x, y + other.y);
}
XY XY::operator-(const XY &other) const {
return XY(x - other.x, y - other.y);
}
std::ostream &operator<<(std::ostream &os, const XY &xy) {
return os << '(' << xy.x << ' ' << xy.y << ')';
}
ContourLine::ContourLine(bool is_hole)
: std::vector<XY>(), _is_hole(is_hole), _parent(0) {}
void ContourLine::add_child(ContourLine *child) {
assert(!_is_hole && "Cannot add_child to a hole");
assert(child != 0 && "Null child ContourLine");
_children.push_back(child);
}
void ContourLine::clear_parent() {
assert(is_hole() && "Cannot clear parent of non-hole");
assert(_parent != 0 && "Null parent ContourLine");
_parent = 0;
}
const ContourLine::Children &ContourLine::get_children() const {
assert(!_is_hole && "Cannot get_children of a hole");
return _children;
}
const ContourLine *ContourLine::get_parent() const {
assert(_is_hole && "Cannot get_parent of a non-hole");
return _parent;
}
ContourLine *ContourLine::get_parent() {
assert(_is_hole && "Cannot get_parent of a non-hole");
return _parent;
}
bool ContourLine::is_hole() const { return _is_hole; }
void ContourLine::push_back(const XY &point) {
if (empty() || point != back())
std::vector<XY>::push_back(point);
}
void ContourLine::set_parent(ContourLine *parent) {
assert(_is_hole && "Cannot set parent of a non-hole");
assert(parent != 0 && "Null parent ContourLine");
_parent = parent;
}
void ContourLine::write() const {
std::cout << "ContourLine " << this << " of " << size() << " points:";
for (const_iterator it = begin(); it != end(); ++it)
std::cout << ' ' << *it;
if (is_hole())
std::cout << " hole, parent=" << get_parent();
else {
std::cout << " not hole";
if (!_children.empty()) {
std::cout << ", children=";
for (Children::const_iterator it = _children.begin();
it != _children.end(); ++it)
std::cout << *it << ' ';
}
}
std::cout << std::endl;
}
Contour::Contour() {}
Contour::~Contour() { delete_contour_lines(); }
void Contour::delete_contour_lines() {
for (iterator line_it = begin(); line_it != end(); ++line_it) {
delete *line_it;
*line_it = 0;
}
std::vector<ContourLine *>::clear();
}
void Contour::write() const {
std::cout << "Contour of " << size() << " lines." << std::endl;
for (const_iterator it = begin(); it != end(); ++it)
(*it)->write();
}
ParentCache::ParentCache(long nx, long x_chunk_points, long y_chunk_points)
: _nx(nx), _x_chunk_points(x_chunk_points),
_y_chunk_points(y_chunk_points),
_lines(0), // Initialised when first needed.
_istart(0), _jstart(0) {
assert(_x_chunk_points > 0 && _y_chunk_points > 0 &&
"Chunk sizes must be positive");
}
ContourLine *ParentCache::get_parent(long quad) {
long index = quad_to_index(quad);
ContourLine *parent = _lines[index];
while (parent == 0) {
index -= _x_chunk_points;
assert(index >= 0 && "Failed to find parent in chunk ParentCache");
parent = _lines[index];
}
assert(parent != 0 && "Failed to find parent in chunk ParentCache");
return parent;
}
long ParentCache::quad_to_index(long quad) const {
long i = quad % _nx;
long j = quad / _nx;
long index = (i - _istart) + (j - _jstart) * _x_chunk_points;
assert(i >= _istart && i < _istart + _x_chunk_points &&
"i-index outside chunk");
assert(j >= _jstart && j < _jstart + _y_chunk_points &&
"j-index outside chunk");
assert(index >= 0 && index < static_cast<long>(_lines.size()) &&
"ParentCache index outside chunk");
return index;
}
void ParentCache::set_chunk_starts(long istart, long jstart) {
assert(istart >= 0 && jstart >= 0 &&
"Chunk start indices cannot be negative");
_istart = istart;
_jstart = jstart;
if (_lines.empty())
_lines.resize(_x_chunk_points * _y_chunk_points, 0);
else
std::fill(_lines.begin(), _lines.end(), (ContourLine *)0);
}
void ParentCache::set_parent(long quad, ContourLine &contour_line) {
assert(!_lines.empty() &&
"Accessing ParentCache before it has been initialised");
long index = quad_to_index(quad);
if (_lines[index] == 0)
_lines[index] = (contour_line.is_hole() ? contour_line.get_parent()
: &contour_line);
}
QuadContourGenerator::QuadContourGenerator(const CoordinateArray &x,
const CoordinateArray &y,
const CoordinateArray &z,
bool corner_mask,
long chunk_size)
: _x(x), _y(y), _z(z), _nx(static_cast<long>(_x[0].size())),
_ny(static_cast<long>(_x.size())), _n(_nx * _ny),
_corner_mask(corner_mask),
_chunk_size(chunk_size > 0
? std::min(chunk_size, std::max(_nx, _ny) - 1)
: std::max(_nx, _ny) - 1),
_nxchunk(calc_chunk_count(_nx)), _nychunk(calc_chunk_count(_ny)),
_chunk_count(_nxchunk * _nychunk), _cache(std::vector<CacheItem>(_n)),
_parent_cache(_nx, chunk_size > 0 ? chunk_size + 1 : _nx,
chunk_size > 0 ? chunk_size + 1 : _ny) {
assert(!_x.empty() && !_y.empty() && !_z.empty() && "Empty array");
assert(_y.size() == _x.size() && _y[0].size() == _x[0].size() &&
"Different-sized y and x arrays");
assert(_z.size() == _x.size() && _z[0].size() == _x[0].size() &&
"Different-sized z and x arrays");
init_cache_grid();
}
void QuadContourGenerator::append_contour_line_to_vertices(
ContourLine &contour_line, vertices_list_type &vertices_list) const {
double x_diff = std::abs(_x[0][1] - _x[0][0]);
double y_diff = std::abs(_y[1][0] - _y[0][0]);
// Convert ContourLine to vertices_list
size_t i = 0;
size_t inserted = 0;
for (ContourLine::const_iterator point = contour_line.begin();
point != contour_line.end(); ++point, ++i) {
bool is_origin = point->x == 0. && point->y == 0.;
bool is_x_jump_to_origin =
is_origin && inserted != 0 &&
std::abs(vertices_list.first.back()) > 3 * x_diff;
bool is_y_jump_to_origin =
is_origin && inserted != 0 &&
std::abs(vertices_list.second.back()) > 3 * y_diff;
bool is_end_of_segment = is_x_jump_to_origin && is_y_jump_to_origin;
if (!is_end_of_segment) {
vertices_list.first.emplace_back(point->x);
vertices_list.second.emplace_back(point->y);
++inserted;
} else {
bool last_segment_is_not_empty =
inserted > 0 && std::isfinite(vertices_list.first.back()) &&
std::isfinite(vertices_list.first.back());
if (last_segment_is_not_empty) {
vertices_list.first.emplace_back(NaN);
vertices_list.second.emplace_back(NaN);
}
}
}
bool last_segment_is_not_empty =
inserted > 0 && std::isfinite(vertices_list.first.back()) &&
std::isfinite(vertices_list.first.back());
if (last_segment_is_not_empty) {
vertices_list.first.emplace_back(NaN);
vertices_list.second.emplace_back(NaN);
}
contour_line.clear();
}
void QuadContourGenerator::append_contour_to_vertices_and_codes(
Contour &contour, vertices_list_type &vertices_list,
codes_list_type &codes_list) const {
// Convert Contour to vertices and clear it.
// For each line in the contour
for (Contour::iterator line_it = contour.begin();
line_it != contour.end(); ++line_it) {
ContourLine &line = **line_it;
if (line.is_hole()) {
// If hole has already been converted to python its parent will
// be set to 0 and it can be deleted.
if (line.get_parent() != 0) {
delete *line_it;
*line_it = 0;
}
} else {
// Non-holes are converted to python together with their child
// holes so that they are rendered correctly.
ContourLine::const_iterator point;
ContourLine::Children::const_iterator children_it;
const ContourLine::Children &children = line.get_children();
// number of points in the line itself
size_t npoints = static_cast<size_t>(line.size() + 1);
// accumulate the number of points for each line child
for (children_it = children.begin();
children_it != children.end(); ++children_it) {
npoints += static_cast<size_t>((*children_it)->size() + 1);
}
auto is_valid_point = [](double x, double y) {
return std::isfinite(x) && std::isfinite(y) && x != 0. &&
y != 0.;
};
// for each point in this line
bool first = true;
for (point = line.begin(); point != line.end(); ++point) {
if (is_valid_point(point->x, point->y)) {
vertices_list.first.emplace_back(point->x);
vertices_list.second.emplace_back(point->y);
codes_list.emplace_back(static_cast<unsigned char>(
first ? MOVETO : LINETO));
first = false;
}
}
// Close polygon
// Emplace NaN instead of the first point again
// Because that's how gnuplot works.
// point = line.begin();
// vertices_list.first.emplace_back(point->x);
// vertices_list.second.emplace_back(point->y);
if (!first) {
vertices_list.first.emplace_back(NaN);
vertices_list.second.emplace_back(NaN);
codes_list.emplace_back(
static_cast<unsigned char>(CLOSEPOLY));
} else {
// We inserted no points
// ignore children (empty parents cannot have children)
for (children_it = children.begin();
children_it != children.end(); ++children_it) {
ContourLine &child = **children_it;
child.clear_parent();
}
// free parent
delete *line_it;
*line_it = 0;
// and go to next line
continue;
}
// Add the points from the children segments / holes
for (children_it = children.begin();
children_it != children.end(); ++children_it) {
ContourLine &child = **children_it;
bool never_emplaced_anything = true;
for (point = child.begin(); point != child.end(); ++point) {
if (is_valid_point(point->x, point->y)) {
vertices_list.first.emplace_back(point->x);
vertices_list.second.emplace_back(point->y);
codes_list.emplace_back(static_cast<unsigned char>(
(never_emplaced_anything ? MOVETO : LINETO)));
never_emplaced_anything = false;
}
}
// Close the polygon
// Emplace NaN instead of the first point again
// That's how gnuplot works
// point = child.begin();
// vertices_list.first.emplace_back(point->x);
// vertices_list.second.emplace_back(point->y);
if (!never_emplaced_anything) {
vertices_list.first.emplace_back(NaN);
vertices_list.second.emplace_back(NaN);
codes_list.emplace_back(
static_cast<unsigned char>(CLOSEPOLY));
}
// To indicate the parent can be deleted.
child.clear_parent();
}
// Include extra NaNs after all children
// This indicates the children are over
// and that the next polygon is another parent
vertices_list.first.emplace_back(NaN);
vertices_list.second.emplace_back(NaN);
codes_list.emplace_back(static_cast<unsigned char>(CLOSEPOLY));
delete *line_it;
*line_it = 0;
}
}
// Delete remaining contour lines.
contour.delete_contour_lines();
}
long QuadContourGenerator::calc_chunk_count(long point_count) const {
assert(point_count > 0 && "point count must be positive");
assert(_chunk_size > 0 && "Chunk size must be positive");
if (_chunk_size > 0) {
long count = (point_count - 1) / _chunk_size;
if (count * _chunk_size < point_count - 1)
++count;
assert(count >= 1 && "Invalid chunk count");
return count;
} else
return 1;
}
QuadContourGenerator::vertices_list_type
QuadContourGenerator::create_contour(const double &level) {
init_cache_levels(level, level);
vertices_list_type vertices_list;
// Lines that start and end on boundaries.
long ichunk, jchunk, istart, iend, jstart, jend;
for (long ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) {
get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart,
jend);
for (long j = jstart; j < jend; ++j) {
long quad_end = iend + j * _nx;
for (long quad = istart + j * _nx; quad < quad_end; ++quad) {
if (EXISTS_NONE(quad) || VISITED(quad, 1))
continue;
if (BOUNDARY_S(quad) && Z_SW >= 1 && Z_SE < 1 &&
start_line(vertices_list, quad, Edge_S, level))
continue;
if (BOUNDARY_W(quad) && Z_NW >= 1 && Z_SW < 1 &&
start_line(vertices_list, quad, Edge_W, level))
continue;
if (BOUNDARY_N(quad) && Z_NE >= 1 && Z_NW < 1 &&
start_line(vertices_list, quad, Edge_N, level))
continue;
if (BOUNDARY_E(quad) && Z_SE >= 1 && Z_NE < 1 &&
start_line(vertices_list, quad, Edge_E, level))
continue;
if (_corner_mask) {
// Equates to NE boundary.
if (EXISTS_SW_CORNER(quad) && Z_SE >= 1 && Z_NW < 1 &&
start_line(vertices_list, quad, Edge_NE, level))
continue;
// Equates to NW boundary.
if (EXISTS_SE_CORNER(quad) && Z_NE >= 1 && Z_SW < 1 &&
start_line(vertices_list, quad, Edge_NW, level))
continue;
// Equates to SE boundary.
if (EXISTS_NW_CORNER(quad) && Z_SW >= 1 && Z_NE < 1 &&
start_line(vertices_list, quad, Edge_SE, level))
continue;
// Equates to SW boundary.
if (EXISTS_NE_CORNER(quad) && Z_NW >= 1 && Z_SE < 1 &&
start_line(vertices_list, quad, Edge_SW, level))
continue;
}
}
}
}
// Internal loops.
ContourLine contour_line(false); // Reused for each contour line.
for (long ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) {
get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart,
jend);
for (long j = jstart; j < jend; ++j) {
long quad_end = iend + j * _nx;
for (long quad = istart + j * _nx; quad < quad_end; ++quad) {
if (EXISTS_NONE(quad) || VISITED(quad, 1))
continue;
Edge start_edge = get_start_edge(quad, 1);
if (start_edge == Edge_None)
continue;
QuadEdge quad_edge(quad, start_edge);
QuadEdge start_quad_edge(quad_edge);
// To obtain output identical to that produced by legacy
// code, sometimes need to ignore the first point and add it
// on the end instead.
bool ignore_first = (start_edge == Edge_N);
follow_interior(contour_line, quad_edge, 1, level,
!ignore_first, &start_quad_edge, 1, false);
if (ignore_first && !contour_line.empty())
contour_line.push_back(contour_line.front());
append_contour_line_to_vertices(contour_line,
vertices_list);
// Repeat if saddle point but not visited.
if (SADDLE(quad, 1) && !VISITED(quad, 1))
--quad;
}
}
}
return vertices_list;
}
std::pair<QuadContourGenerator::vertices_list_type,
QuadContourGenerator::codes_list_type>
QuadContourGenerator::create_filled_contour(const double &lower_level,
const double &upper_level) {
init_cache_levels(lower_level, upper_level);
Contour contour;
vertices_list_type vertices;
codes_list_type codes;
// for each chunk
long ichunk, jchunk, istart, iend, jstart, jend;
for (long ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) {
get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart,
jend);
_parent_cache.set_chunk_starts(istart, jstart);
// iterate quads
for (long j = jstart; j < jend; ++j) {
long quad_end = iend + j * _nx;
for (long quad = istart + j * _nx; quad < quad_end; ++quad) {
if (!EXISTS_NONE(quad))
single_quad_filled(contour, quad, lower_level,
upper_level);
}
}
// Clear VISITED_W and VISITED_S flags that are reused by later
// chunks.
if (jchunk < _nychunk - 1) {
long quad_end = iend + jend * _nx;
for (long quad = istart + jend * _nx; quad < quad_end; ++quad)
_cache[quad] &= ~MASK_VISITED_S;
}
if (ichunk < _nxchunk - 1) {
long quad_end = iend + jend * _nx;
for (long quad = iend + jstart * _nx; quad < quad_end;
quad += _nx)
_cache[quad] &= ~MASK_VISITED_W;
}
// Create python objects to return for this chunk.
append_contour_to_vertices_and_codes(contour, vertices, codes);
}
return std::make_pair(vertices, codes);
}
XY QuadContourGenerator::edge_interp(const QuadEdge &quad_edge,
const double &level) {
assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
"Quad index out of bounds");
assert(quad_edge.edge != Edge_None && "Invalid edge");
return interp(get_edge_point_index(quad_edge, true),
get_edge_point_index(quad_edge, false), level);
}
unsigned int QuadContourGenerator::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) {
assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
"Quad index out of bounds");
assert(quad_edge.edge != Edge_None && "Invalid edge");
assert(is_edge_a_boundary(quad_edge) && "Not a boundary edge");
assert((level_index == 1 || level_index == 2) &&
"level index must be 1 or 2");
assert(start_quad_edge.quad >= 0 && start_quad_edge.quad < _n &&
"Start quad index out of bounds");
assert(start_quad_edge.edge != Edge_None && "Invalid start edge");
// Only called for filled contours, so always updates _parent_cache.
unsigned int end_level = 0;
bool first_edge = true;
bool stop = false;
long &quad = quad_edge.quad;
while (true) {
// Levels of start and end points of quad_edge.
unsigned int start_level =
(first_edge ? Z_LEVEL(get_edge_point_index(quad_edge, true))
: end_level);
long end_point = get_edge_point_index(quad_edge, false);
end_level = Z_LEVEL(end_point);
if (level_index == 1) {
if (start_level <= level_index && end_level == 2) {
// Increasing z, switching levels from 1 to 2.
level_index = 2;
stop = true;
} else if (start_level >= 1 && end_level == 0) {
// Decreasing z, keeping same level.
stop = true;
}
} else { // level_index == 2
if (start_level <= level_index && end_level == 2) {
// Increasing z, keeping same level.
stop = true;
} else if (start_level >= 1 && end_level == 0) {
// Decreasing z, switching levels from 2 to 1.
level_index = 1;
stop = true;
}
}
if (!first_edge && !stop && quad_edge == start_quad_edge)
// Return if reached start point of contour line. Do this
// before checking/setting VISITED flags as will already have
// been visited.
break;
switch (quad_edge.edge) {
case Edge_E:
assert(!VISITED_W(quad + 1) && "Already visited");
_cache[quad + 1] |= MASK_VISITED_W;
break;
case Edge_N:
assert(!VISITED_S(quad + _nx) && "Already visited");
_cache[quad + _nx] |= MASK_VISITED_S;
break;
case Edge_W:
assert(!VISITED_W(quad) && "Already visited");
_cache[quad] |= MASK_VISITED_W;
break;
case Edge_S:
assert(!VISITED_S(quad) && "Already visited");
_cache[quad] |= MASK_VISITED_S;
break;
case Edge_NE:
case Edge_NW:
case Edge_SW:
case Edge_SE:
assert(!VISITED_CORNER(quad) && "Already visited");
_cache[quad] |= MASK_VISITED_CORNER;
break;
default:
assert(0 && "Invalid Edge");
break;
}
if (stop) {
// Exiting boundary to enter interior.
contour_line.push_back(edge_interp(
quad_edge, level_index == 1 ? lower_level : upper_level));
break;
}
move_to_next_boundary_edge(quad_edge);
// Just moved to new quad edge, so label parent of start of quad
// edge.
switch (quad_edge.edge) {
case Edge_W:
case Edge_SW:
case Edge_S:
case Edge_SE:
if (!EXISTS_SE_CORNER(quad))
_parent_cache.set_parent(quad, contour_line);
break;
case Edge_E:
case Edge_NE:
case Edge_N:
case Edge_NW:
if (!EXISTS_SW_CORNER(quad))
_parent_cache.set_parent(quad + 1, contour_line);
break;
default:
assert(0 && "Invalid edge");
break;
}
// Add point to contour.
contour_line.push_back(get_point_xy(end_point));
if (first_edge)
first_edge = false;
}
return level_index;
}
void QuadContourGenerator::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) {
assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
"Quad index out of bounds.");
assert(quad_edge.edge != Edge_None && "Invalid edge");
assert((level_index == 1 || level_index == 2) &&
"level index must be 1 or 2");
assert((start_quad_edge == 0 ||
(start_quad_edge->quad >= 0 && start_quad_edge->quad < _n)) &&
"Start quad index out of bounds.");
assert((start_quad_edge == 0 || start_quad_edge->edge != Edge_None) &&
"Invalid start edge");
assert((start_level_index == 1 || start_level_index == 2) &&
"start level index must be 1 or 2");
long &quad = quad_edge.quad;
Edge &edge = quad_edge.edge;
if (want_initial_point)
contour_line.push_back(edge_interp(quad_edge, level));
CacheItem visited_mask =
(level_index == 1 ? MASK_VISITED_1 : MASK_VISITED_2);
CacheItem saddle_mask =
(level_index == 1 ? MASK_SADDLE_1 : MASK_SADDLE_2);
Dir dir = Dir_Straight;
while (true) {
assert(!EXISTS_NONE(quad) && "Quad does not exist");
assert(!(_cache[quad] & visited_mask) && "Quad already visited");
// Determine direction to move to next quad. If the quad is already
// labelled as a saddle quad then the direction is easily read from
// the cache. Otherwise the direction is determined differently
// depending on whether the quad is a corner quad or not.
if (_cache[quad] & saddle_mask) {
// Already identified as a saddle quad, so direction is easy.
dir = (SADDLE_LEFT(quad, level_index) ? Dir_Left : Dir_Right);
_cache[quad] |= visited_mask;
} else if (EXISTS_ANY_CORNER(quad)) {
// Need z-level of point opposite the entry edge, as that
// determines whether contour turns left or right.
long point_opposite = -1;
switch (edge) {
case Edge_E:
point_opposite =
(EXISTS_SE_CORNER(quad) ? POINT_SW : POINT_NW);
break;
case Edge_N:
point_opposite =
(EXISTS_NW_CORNER(quad) ? POINT_SW : POINT_SE);
break;
case Edge_W:
point_opposite =
(EXISTS_SW_CORNER(quad) ? POINT_SE : POINT_NE);
break;
case Edge_S:
point_opposite =
(EXISTS_SW_CORNER(quad) ? POINT_NW : POINT_NE);
break;
case Edge_NE:
point_opposite = POINT_SW;
break;
case Edge_NW:
point_opposite = POINT_SE;
break;
case Edge_SW:
point_opposite = POINT_NE;
break;
case Edge_SE:
point_opposite = POINT_NW;
break;
default:
assert(0 && "Invalid edge");
break;
}
assert(point_opposite != -1 && "Failed to find opposite point");
// Lower-level polygons (level_index == 1) always have higher
// values to the left of the contour. Upper-level contours
// (level_index == 2) are reversed, which is what the fancy XOR
// does below.
if ((Z_LEVEL(point_opposite) >= level_index) ^
(level_index == 2))
dir = Dir_Right;
else
dir = Dir_Left;
_cache[quad] |= visited_mask;
} else {
// Calculate configuration of this quad.
long point_left = -1, point_right = -1;
switch (edge) {
case Edge_E:
point_left = POINT_SW;
point_right = POINT_NW;
break;
case Edge_N:
point_left = POINT_SE;
point_right = POINT_SW;
break;
case Edge_W:
point_left = POINT_NE;
point_right = POINT_SE;
break;
case Edge_S:
point_left = POINT_NW;
point_right = POINT_NE;
break;
default:
assert(0 && "Invalid edge");
break;
}
unsigned int config = (Z_LEVEL(point_left) >= level_index)
<< 1 |
(Z_LEVEL(point_right) >= level_index);
// Upper level (level_index == 2) polygons are reversed compared
// to lower level ones, i.e. higher values on the right rather
// than the left.
if (level_index == 2)
config = 3 - config;
// Calculate turn direction to move to next quad along contour
// line.
if (config == 1) {
// New saddle quad, set up cache bits for it.
double zmid =
0.25 * (get_point_z(POINT_SW) + get_point_z(POINT_SE) +
get_point_z(POINT_NW) + get_point_z(POINT_NE));
_cache[quad] |=
(level_index == 1 ? MASK_SADDLE_1 : MASK_SADDLE_2);
if ((zmid > level) ^ (level_index == 2)) {
dir = Dir_Right;
} else {
dir = Dir_Left;
_cache[quad] |= (level_index == 1 ? MASK_SADDLE_LEFT_1
: MASK_SADDLE_LEFT_2);
}
if (edge == Edge_N || edge == Edge_E) {
// Next visit to this quad must start on S or W.
_cache[quad] |=
(level_index == 1 ? MASK_SADDLE_START_SW_1
: MASK_SADDLE_START_SW_2);
}
} else {
// Normal (non-saddle) quad.
dir = (config == 0
? Dir_Left
: (config == 3 ? Dir_Right : Dir_Straight));
_cache[quad] |= visited_mask;
}
}
// Use dir to determine exit edge.
edge = get_exit_edge(quad_edge, dir);
if (set_parents) {
if (edge == Edge_E)
_parent_cache.set_parent(quad + 1, contour_line);
else if (edge == Edge_W)
_parent_cache.set_parent(quad, contour_line);
}
// Add new point to contour line.
contour_line.push_back(edge_interp(quad_edge, level));
// Stop if reached boundary.
if (is_edge_a_boundary(quad_edge))
break;
move_to_next_quad(quad_edge);
assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&
"Quad index out of bounds");
// Return if reached start point of contour line.
if (start_quad_edge != 0 && quad_edge == *start_quad_edge &&
level_index == start_level_index)
break;
}
}
void QuadContourGenerator::get_chunk_limits(long ijchunk, long &ichunk,
long &jchunk, long &istart,
long &iend, long &jstart,
long &jend) {
assert(ijchunk >= 0 && ijchunk < _chunk_count &&
"ijchunk out of bounds");
ichunk = ijchunk % _nxchunk;
jchunk = ijchunk / _nxchunk;
istart = ichunk * _chunk_size;
iend = (ichunk == _nxchunk - 1 ? _nx : (ichunk + 1) * _chunk_size);
jstart = jchunk * _chunk_size;
jend = (jchunk == _nychunk - 1 ? _ny : (jchunk + 1) * _chunk_size);
}
Edge QuadContourGenerator::get_corner_start_edge(
long quad, unsigned int level_index) const {
assert(quad >= 0 && quad < _n && "Quad index out of bounds");
assert((level_index == 1 || level_index == 2) &&
"level index must be 1 or 2");
assert(EXISTS_ANY_CORNER(quad) && "Quad is not a corner");
// Diagram for NE corner. Rotate for other corners.
//
// edge12