diff options
| author | maxim-yurchuk <maxim-yurchuk@yandex-team.com> | 2025-02-11 13:26:52 +0300 |
|---|---|---|
| committer | maxim-yurchuk <maxim-yurchuk@yandex-team.com> | 2025-02-11 13:57:59 +0300 |
| commit | f895bba65827952ed934b2b46f9a45e30a191fd2 (patch) | |
| tree | 03260c906d9ec41cdc03e2a496b15d407459cec0 /contrib/python/contourpy/src | |
| parent | 5f7060466f7b9707818c2091e1a25c14f33c3474 (diff) | |
| download | ydb-f895bba65827952ed934b2b46f9a45e30a191fd2.tar.gz | |
Remove deps on pandas
<https://github.com/ydb-platform/ydb/pull/14418>
<https://github.com/ydb-platform/ydb/pull/14419>
\-- аналогичные правки в gh
Хочу залить в обход синка, чтобы посмотреть удалится ли pandas в нашей gh репе через piglet
commit_hash:abca127aa37d4dbb94b07e1e18cdb8eb5b711860
Diffstat (limited to 'contrib/python/contourpy/src')
32 files changed, 0 insertions, 8441 deletions
diff --git a/contrib/python/contourpy/src/base.h b/contrib/python/contourpy/src/base.h deleted file mode 100644 index 6b4d99d17ea..00000000000 --- a/contrib/python/contourpy/src/base.h +++ /dev/null @@ -1,217 +0,0 @@ -// BaseContourGenerator class provides functionality common to multiple contouring algorithms. -// Uses the Curiously Recurring Template Pattern idiom whereby classes derived from this are -// declared as -// class Derived : public BaseContourGenerator<Derived> -// -// It provides static polymorphism at compile-time rather than the more common dynamic polymorphism -// at run-time using virtual functions. Hence it avoids the runtime overhead of calling virtual -// functions. - -#ifndef CONTOURPY_BASE_H -#define CONTOURPY_BASE_H - -#include "chunk_local.h" -#include "contour_generator.h" -#include "fill_type.h" -#include "line_type.h" -#include "outer_or_hole.h" -#include "z_interp.h" -#include <vector> - -namespace contourpy { - -template <typename Derived> -class BaseContourGenerator : public ContourGenerator -{ -public: - virtual ~BaseContourGenerator(); - - static FillType default_fill_type(); - static LineType default_line_type(); - - py::tuple filled(double lower_level, double upper_level) override; - - py::tuple get_chunk_count() const; // Return (y_chunk_count, x_chunk_count) - py::tuple get_chunk_size() const; // Return (y_chunk_size, x_chunk_size) - - bool get_corner_mask() const; - - FillType get_fill_type() const; - LineType get_line_type() const; - - bool get_quad_as_tri() const; - - ZInterp get_z_interp() const; - - py::sequence lines(double level) override; - - py::list multi_filled(const LevelArray levels) override; - py::list multi_lines(const LevelArray levels) override; - - static bool supports_fill_type(FillType fill_type); - static bool supports_line_type(LineType line_type); - - void write_cache() const; // For debug purposes only. - -protected: - BaseContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, bool corner_mask, LineType line_type, FillType fill_type, - bool quad_as_tri, ZInterp z_interp, index_t x_chunk_size, index_t y_chunk_size); - - typedef uint32_t CacheItem; - typedef CacheItem ZLevel; - - // C++11 scoped enum for direction of movement from one quad to the next. - enum class Direction - { - Left = 1, - Straight = 0, - Right = -1 - }; - - struct Location - { - Location(index_t quad_, index_t forward_, index_t left_, bool is_upper_, bool on_boundary_) - : quad(quad_), forward(forward_), left(left_), is_upper(is_upper_), - on_boundary(on_boundary_) - {} - - friend std::ostream &operator<<(std::ostream &os, const Location& location) - { - os << "quad=" << location.quad << " forward=" << location.forward << " left=" - << location.left << " is_upper=" << location.is_upper << " on_boundary=" - << location.on_boundary; - return os; - } - - index_t quad, forward, left; - bool is_upper, on_boundary; - }; - - // Calculate, set and return z-level at middle of quad. - ZLevel calc_and_set_middle_z_level(index_t quad); - - // Calculate and return z at middle of quad. - double calc_middle_z(index_t quad) const; - - void closed_line(const Location& start_location, OuterOrHole outer_or_hole, ChunkLocal& local); - - void closed_line_wrapper( - const Location& start_location, OuterOrHole outer_or_hole, ChunkLocal& local); - - // If point/line/hole counts not consistent, throw runtime error. - void check_consistent_counts(const ChunkLocal& local) const; - - index_t find_look_S(index_t look_N_quad) const; - - // Return true if finished (i.e. back to start quad, direction and upper). - bool follow_boundary( - Location& location, const Location& start_location, ChunkLocal& local, - count_t& point_count); - - // Return true if finished (i.e. back to start quad, direction and upper). - bool follow_interior( - Location& location, const Location& start_location, ChunkLocal& local, - count_t& point_count); - - index_t get_boundary_start_point(const Location& location) const; - - // These are quad chunk limits, not point chunk limits. - // chunk is index in range 0.._n_chunks-1. - void get_chunk_limits(index_t chunk, ChunkLocal& local) const; - - index_t get_interior_start_left_point( - const Location& location, bool& start_corner_diagonal) const; - - double get_interp_fraction(double z0, double z1, double level) const; - - double get_middle_x(index_t quad) const; - double get_middle_y(index_t quad) const; - - index_t get_n_chunks() const; - - void get_point_xy(index_t point, double*& points) const; - - double get_point_x(index_t point) const; - double get_point_y(index_t point) const; - double get_point_z(index_t point) const; - - bool has_direct_line_offsets() const; - bool has_direct_outer_offsets() const; - bool has_direct_points() const; - - void init_cache_grid(const MaskArray& mask); - - // Either for a single chunk, or the whole domain (all chunks) if local == nullptr. - void init_cache_levels_and_starts(const ChunkLocal* local = nullptr); - - // Increments local.points twice. - void interp(index_t point0, index_t point1, bool is_upper, double*& points) const; - - // Increments local.points twice. - void interp( - index_t point0, double x1, double y1, double z1, bool is_upper, double*& points) const; - - bool is_point_in_chunk(index_t point, const ChunkLocal& local) const; - - bool is_quad_in_bounds( - index_t quad, index_t istart, index_t iend, index_t jstart, index_t jend) const; - - bool is_quad_in_chunk(index_t quad, const ChunkLocal& local) const; - - void line(const Location& start_location, ChunkLocal& local); - - void march_chunk(ChunkLocal& local, std::vector<py::list>& return_lists); - - py::sequence march_wrapper(); - - void move_to_next_boundary_edge(index_t& quad, index_t& forward, index_t& left) const; - - // Common initialisation for filled/multi_filled and lines/multi_lines calls. - void pre_filled(); - void pre_lines(); - - void set_look_flags(index_t hole_start_quad); - - void write_cache_quad(index_t quad) const; - - ZLevel z_to_zlevel(double z_value) const; - - -private: - const CoordinateArray _x, _y, _z; - const double* _xptr; // For quick access to _x.data(). - const double* _yptr; - const double* _zptr; - const index_t _nx, _ny; // Number of points in each direction. - const index_t _n; // Total number of points (and quads). - const index_t _x_chunk_size, _y_chunk_size; - const index_t _nx_chunks, _ny_chunks; // Number of chunks in each direction. - const index_t _n_chunks; // Total number of chunks. - const bool _corner_mask; - const LineType _line_type; - const FillType _fill_type; - const bool _quad_as_tri; - const ZInterp _z_interp; - - CacheItem* _cache; - - // Current contouring operation. - bool _filled; - double _lower_level, _upper_level; - - // Current contouring operation, based on return type and filled or lines. - bool _identify_holes; - bool _output_chunked; // Implies empty chunks will have py::none(). - bool _direct_points; // Whether points array is written direct to Python. - bool _direct_line_offsets; // Whether line offsets array is written direct to Python. - bool _direct_outer_offsets; // Whether outer offsets array is written direct to Python. - bool _outer_offsets_into_points; // Otherwise into line offsets. Only used if _identify_holes. - bool _nan_separated; // Whether adjacent lines' points are separated by nans. - unsigned int _return_list_count; -}; - -} // namespace contourpy - -#endif // CONTOURPY_BASE_H diff --git a/contrib/python/contourpy/src/base_impl.h b/contrib/python/contourpy/src/base_impl.h deleted file mode 100644 index 80c9e51854a..00000000000 --- a/contrib/python/contourpy/src/base_impl.h +++ /dev/null @@ -1,2524 +0,0 @@ -#ifndef CONTOURPY_BASE_IMPL_H -#define CONTOURPY_BASE_IMPL_H - -#include "base.h" -#include "converter.h" -#include "util.h" -#include <iostream> - -namespace contourpy { - -// Point indices from current quad index. -#define POINT_NE (quad) -#define POINT_NW (quad-1) -#define POINT_SE (quad-_nx) -#define POINT_SW (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_1 (0x1 << 0) // z > lower_level. -#define MASK_Z_LEVEL_2 (0x1 << 1) // z > upper_level. -#define MASK_Z_LEVEL (MASK_Z_LEVEL_1 | MASK_Z_LEVEL_2) -#define MASK_MIDDLE_Z_LEVEL_1 (0x1 << 2) // middle z > lower_level -#define MASK_MIDDLE_Z_LEVEL_2 (0x1 << 3) // middle z > upper_level -#define MASK_MIDDLE (MASK_MIDDLE_Z_LEVEL_1 | MASK_MIDDLE_Z_LEVEL_2) -#define MASK_BOUNDARY_E (0x1 << 4) // E edge of quad is a boundary. -#define MASK_BOUNDARY_N (0x1 << 5) // N 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. -#define MASK_EXISTS_QUAD (0x1 << 6) // All of quad exists (is not masked). -#define MASK_EXISTS_NE_CORNER (0x1 << 7) // NE corner exists, SW corner is masked. -#define MASK_EXISTS_NW_CORNER (0x1 << 8) -#define MASK_EXISTS_SE_CORNER (0x1 << 9) -#define MASK_EXISTS_SW_CORNER (0x1 << 10) -#define MASK_EXISTS_ANY_CORNER (MASK_EXISTS_NE_CORNER | MASK_EXISTS_NW_CORNER | MASK_EXISTS_SE_CORNER | MASK_EXISTS_SW_CORNER) -#define MASK_EXISTS_ANY (MASK_EXISTS_QUAD | MASK_EXISTS_ANY_CORNER) -#define MASK_START_E (0x1 << 11) // E to N, filled and lines. -#define MASK_START_N (0x1 << 12) // N to E, filled and lines. -#define MASK_START_BOUNDARY_E (0x1 << 13) // Lines only. -#define MASK_START_BOUNDARY_N (0x1 << 14) // Lines only. -#define MASK_START_BOUNDARY_S (0x1 << 15) // Filled and lines. -#define MASK_START_BOUNDARY_W (0x1 << 16) // Filled and lines. -#define MASK_START_CORNER (0x1 << 18) // Filled and lines. -#define MASK_START_HOLE_N (0x1 << 17) // N boundary of EXISTS, E to W, filled only. -#define MASK_ANY_START (MASK_START_N | MASK_START_E | MASK_START_BOUNDARY_N | MASK_START_BOUNDARY_E | MASK_START_BOUNDARY_S | MASK_START_BOUNDARY_W | MASK_START_HOLE_N | MASK_START_CORNER) -#define MASK_LOOK_N (0x1 << 19) -#define MASK_LOOK_S (0x1 << 20) -#define MASK_NO_STARTS_IN_ROW (0x1 << 21) -#define MASK_NO_MORE_STARTS (0x1 << 22) - -// Accessors for various CacheItem masks. -#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 MIDDLE_Z_LEVEL(quad) ((_cache[quad] & MASK_MIDDLE) >> 2) -#define BOUNDARY_E(quad) (_cache[quad] & MASK_BOUNDARY_E) -#define BOUNDARY_N(quad) (_cache[quad] & MASK_BOUNDARY_N) -#define BOUNDARY_S(quad) (_cache[quad-_nx] & MASK_BOUNDARY_N) -#define BOUNDARY_W(quad) (_cache[quad-1] & MASK_BOUNDARY_E) -#define EXISTS_QUAD(quad) (_cache[quad] & MASK_EXISTS_QUAD) -#define EXISTS_NE_CORNER(quad) (_cache[quad] & MASK_EXISTS_NE_CORNER) -#define EXISTS_NW_CORNER(quad) (_cache[quad] & MASK_EXISTS_NW_CORNER) -#define EXISTS_SE_CORNER(quad) (_cache[quad] & MASK_EXISTS_SE_CORNER) -#define EXISTS_SW_CORNER(quad) (_cache[quad] & MASK_EXISTS_SW_CORNER) -#define EXISTS_ANY(quad) (_cache[quad] & MASK_EXISTS_ANY) -#define EXISTS_ANY_CORNER(quad) (_cache[quad] & MASK_EXISTS_ANY_CORNER) -#define EXISTS_E_EDGE(quad) (_cache[quad] & (MASK_EXISTS_QUAD | MASK_EXISTS_NE_CORNER | MASK_EXISTS_SE_CORNER)) -#define EXISTS_N_EDGE(quad) (_cache[quad] & (MASK_EXISTS_QUAD | MASK_EXISTS_NW_CORNER | MASK_EXISTS_NE_CORNER)) -#define EXISTS_S_EDGE(quad) (_cache[quad] & (MASK_EXISTS_QUAD | MASK_EXISTS_SW_CORNER | MASK_EXISTS_SE_CORNER)) -#define EXISTS_W_EDGE(quad) (_cache[quad] & (MASK_EXISTS_QUAD | MASK_EXISTS_NW_CORNER | MASK_EXISTS_SW_CORNER)) -// Note that EXISTS_NE_CORNER(quad) is equivalent to BOUNDARY_SW(quad), etc. -#define START_E(quad) (_cache[quad] & MASK_START_E) -#define START_N(quad) (_cache[quad] & MASK_START_N) -#define START_BOUNDARY_E(quad) (_cache[quad] & MASK_START_BOUNDARY_E) -#define START_BOUNDARY_N(quad) (_cache[quad] & MASK_START_BOUNDARY_N) -#define START_BOUNDARY_S(quad) (_cache[quad] & MASK_START_BOUNDARY_S) -#define START_BOUNDARY_W(quad) (_cache[quad] & MASK_START_BOUNDARY_W) -#define START_CORNER(quad) (_cache[quad] & MASK_START_CORNER) -#define START_HOLE_N(quad) (_cache[quad] & MASK_START_HOLE_N) -#define ANY_START(quad) ((_cache[quad] & MASK_ANY_START) != 0) -#define LOOK_N(quad) (_cache[quad] & MASK_LOOK_N) -#define LOOK_S(quad) (_cache[quad] & MASK_LOOK_S) -#define NO_STARTS_IN_ROW(quad) (_cache[quad] & MASK_NO_STARTS_IN_ROW) -#define NO_MORE_STARTS(quad) (_cache[quad] & MASK_NO_MORE_STARTS) -// Contour line/fill goes to the left or right of quad middle (quad_as_tri only). -#define LEFT_OF_MIDDLE(quad, is_upper) (MIDDLE_Z_LEVEL(quad) == (is_upper ? 2 : 0)) - - -template <typename Derived> -BaseContourGenerator<Derived>::BaseContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, bool corner_mask, LineType line_type, FillType fill_type, - bool quad_as_tri, ZInterp z_interp, index_t x_chunk_size, index_t y_chunk_size) - : _x(x), - _y(y), - _z(z), - _xptr(_x.data()), - _yptr(_y.data()), - _zptr(_z.data()), - _nx(_z.ndim() > 1 ? _z.shape(1) : 0), - _ny(_z.ndim() > 0 ? _z.shape(0) : 0), - _n(_nx*_ny), - _x_chunk_size(x_chunk_size > 0 ? std::min(x_chunk_size, _nx-1) : _nx-1), - _y_chunk_size(y_chunk_size > 0 ? std::min(y_chunk_size, _ny-1) : _ny-1), - _nx_chunks(static_cast<index_t>(std::ceil((_nx-1.0) / _x_chunk_size))), - _ny_chunks(static_cast<index_t>(std::ceil((_ny-1.0) / _y_chunk_size))), - _n_chunks(_nx_chunks*_ny_chunks), - _corner_mask(corner_mask), - _line_type(line_type), - _fill_type(fill_type), - _quad_as_tri(quad_as_tri), - _z_interp(z_interp), - _cache(new CacheItem[_n]), - _filled(false), - _lower_level(0.0), - _upper_level(0.0), - _identify_holes(false), - _output_chunked(false), - _direct_points(false), - _direct_line_offsets(false), - _direct_outer_offsets(false), - _outer_offsets_into_points(false), - _nan_separated(false), - _return_list_count(0) -{ - if (_x.ndim() != 2 || _y.ndim() != 2 || _z.ndim() != 2) - throw std::invalid_argument("x, y and z must all be 2D arrays"); - - if (_x.shape(1) != _nx || _x.shape(0) != _ny || - _y.shape(1) != _nx || _y.shape(0) != _ny) - throw std::invalid_argument("x, y and z arrays must have the same shape"); - - if (_nx < 2 || _ny < 2) - throw std::invalid_argument("x, y and z must all be at least 2x2 arrays"); - - if (mask.ndim() != 0) { // ndim == 0 if mask is not set, which is valid. - if (mask.ndim() != 2) - throw std::invalid_argument("mask array must be a 2D array"); - - if (mask.shape(1) != _nx || mask.shape(0) != _ny) - throw std::invalid_argument( - "If mask is set it must be a 2D array with the same shape as z"); - } - - if (!supports_line_type(line_type)) - throw std::invalid_argument("Unsupported LineType"); - - if (!supports_fill_type(fill_type)) - throw std::invalid_argument("Unsupported FillType"); - - if (x_chunk_size < 0 || y_chunk_size < 0) - throw std::invalid_argument("x_chunk_size and y_chunk_size cannot be negative"); - - if (_z_interp == ZInterp::Log) { - const bool* mask_ptr = (mask.ndim() == 0 ? nullptr : mask.data()); - for (index_t point = 0; point < _n; ++point) { - if ( (mask_ptr == nullptr || !mask_ptr[point]) && _zptr[point] <= 0.0) - throw std::invalid_argument("z values must be positive if using ZInterp.Log"); - } - } - - init_cache_grid(mask); -} - -template <typename Derived> -BaseContourGenerator<Derived>::~BaseContourGenerator() -{ - delete [] _cache; -} - -template <typename Derived> -typename BaseContourGenerator<Derived>::ZLevel - BaseContourGenerator<Derived>::calc_and_set_middle_z_level(index_t quad) -{ - ZLevel zlevel = z_to_zlevel(calc_middle_z(quad)); - _cache[quad] |= (zlevel << 2); - return zlevel; -} - -template <typename Derived> -double BaseContourGenerator<Derived>::calc_middle_z(index_t quad) const -{ - assert(quad >= 0 && quad < _n); - - switch (_z_interp) { - case ZInterp::Log: - return exp(0.25*(log(get_point_z(POINT_SW)) + - log(get_point_z(POINT_SE)) + - log(get_point_z(POINT_NW)) + - log(get_point_z(POINT_NE)))); - default: // ZInterp::Linear - return 0.25*(get_point_z(POINT_SW) + - get_point_z(POINT_SE) + - get_point_z(POINT_NW) + - get_point_z(POINT_NE)); - } -} - -template <typename Derived> -void BaseContourGenerator<Derived>::check_consistent_counts(const ChunkLocal& local) const -{ - if (local.total_point_count > 0) { - if (local.points.size != 2*local.total_point_count || - local.points.current != local.points.start + 2*local.total_point_count) { - throw std::runtime_error( - "Inconsistent total_point_count for chunk " + std::to_string(local.chunk) + - ". This may indicate a bug in ContourPy."); - } - } - else { - if (local.points.size != 0 || - local.points.start != nullptr || local.points.current != nullptr) { - throw std::runtime_error( - "Inconsistent zero total_point_count for chunk " + std::to_string(local.chunk) + - ". This may indicate a bug in ContourPy."); - } - } - - if (local.line_count > 0) { - if (local.line_offsets.size != local.line_count + 1 || - local.line_offsets.current == nullptr || - local.line_offsets.current != local.line_offsets.start + local.line_count + 1) { - throw std::runtime_error( - "Inconsistent line_count for chunk " + std::to_string(local.chunk) + - ". This may indicate a bug in ContourPy."); - } - } - else { - if (local.line_offsets.size != 0 || - local.line_offsets.start != nullptr || local.line_offsets.current != nullptr) { - throw std::runtime_error( - "Inconsistent zero line_count for chunk " + std::to_string(local.chunk) + - ". This may indicate a bug in ContourPy."); - } - } - - if (_identify_holes && local.line_count > 0) { - if (local.outer_offsets.size != local.line_count - local.hole_count + 1 || - local.outer_offsets.current == nullptr || - local.outer_offsets.current != local.outer_offsets.start + local.line_count - - local.hole_count + 1) { - throw std::runtime_error( - "Inconsistent hole_count for chunk " + std::to_string(local.chunk) + - ". This may indicate a bug in ContourPy."); - } - } - else { - if (local.outer_offsets.size != 0 || - local.outer_offsets.start != nullptr || local.outer_offsets.current != nullptr) { - throw std::runtime_error( - "Inconsistent zero hole_count for chunk " + std::to_string(local.chunk) + - ". This may indicate a bug in ContourPy."); - } - } -} - -template <typename Derived> -void BaseContourGenerator<Derived>::closed_line( - const Location& start_location, OuterOrHole outer_or_hole, ChunkLocal& local) -{ - assert(is_quad_in_chunk(start_location.quad, local)); - - Location location = start_location; - bool finished = false; - count_t point_count = 0; - - if (outer_or_hole == Hole && local.pass == 0 && _identify_holes) - set_look_flags(start_location.quad); - - while (!finished) { - if (location.on_boundary) - finished = follow_boundary(location, start_location, local, point_count); - else - finished = follow_interior(location, start_location, local, point_count); - location.on_boundary = !location.on_boundary; - } - - if (local.pass > 0) { - assert(local.line_offsets.current = local.line_offsets.start + local.line_count); - *local.line_offsets.current++ = local.total_point_count; - if (outer_or_hole == Outer && _identify_holes) { - assert(local.outer_offsets.current == - local.outer_offsets.start + local.line_count - local.hole_count); - if (_outer_offsets_into_points) - *local.outer_offsets.current++ = local.total_point_count; - else - *local.outer_offsets.current++ = local.line_count; - } - } - - local.total_point_count += point_count; - local.line_count++; - if (outer_or_hole == Hole) - local.hole_count++; -} - -template <typename Derived> -void BaseContourGenerator<Derived>::closed_line_wrapper( - const Location& start_location, OuterOrHole outer_or_hole, ChunkLocal& local) -{ - assert(is_quad_in_chunk(start_location.quad, local)); - - if (local.pass == 0 || !_identify_holes) { - closed_line(start_location, outer_or_hole, local); - } - else { - assert(outer_or_hole == Outer); - local.look_up_quads.clear(); - - closed_line(start_location, outer_or_hole, local); - - for (py::size_t i = 0; i < local.look_up_quads.size(); ++i) { - // Note that the collection can increase in size during this loop. - index_t quad = local.look_up_quads[i]; - - // Walk N to corresponding look S flag is reached. - quad = find_look_S(quad); - - // Only 3 possible types of hole start: START_E, START_HOLE_N or START_CORNER for SW - // corner. - if (START_E(quad)) { - closed_line(Location(quad, -1, -_nx, Z_NE > 0, false), Hole, local); - } - else if (START_HOLE_N(quad)) { - closed_line(Location(quad, -1, -_nx, false, true), Hole, local); - } - else { - assert(START_CORNER(quad) && EXISTS_SW_CORNER(quad)); - closed_line(Location(quad, _nx-1, -_nx-1, false, true), Hole, local); - } - } - } -} - -template <typename Derived> -FillType BaseContourGenerator<Derived>::default_fill_type() -{ - FillType fill_type = FillType::OuterOffset; - assert(supports_fill_type(fill_type)); - return fill_type; -} - -template <typename Derived> -LineType BaseContourGenerator<Derived>::default_line_type() -{ - LineType line_type = LineType::Separate; - assert(supports_line_type(line_type)); - return line_type; -} - -template <typename Derived> -py::tuple BaseContourGenerator<Derived>::filled(double lower_level, double upper_level) -{ - check_levels(lower_level, upper_level); - pre_filled(); - - _lower_level = lower_level; - _upper_level = upper_level; - return march_wrapper(); -} - -template <typename Derived> -index_t BaseContourGenerator<Derived>::find_look_S(index_t look_N_quad) const -{ - assert(_identify_holes); - - // Might need to be careful when looking in the same quad as the LOOK_UP. - index_t quad = look_N_quad; - - // look_S quad must have 1 of only 3 possible types of hole start (START_E, START_HOLE_N, - // START_CORNER for SW corner) but it may have other starts as well. - - // Start quad may be both a look_N and look_S quad. Only want to stop search here if look_S - // hole start is N of look_N. - - if (!LOOK_S(quad)) { - do - { - quad += _nx; - assert(quad >= 0 && quad < _n); - assert(EXISTS_ANY(quad)); - } while (!LOOK_S(quad)); - } - - return quad; -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::follow_boundary( - Location& location, const Location& start_location, ChunkLocal& local, count_t& point_count) -{ - // forward values for boundaries: - // -1 = N boundary, E to W. - // 1 = S boundary, W to E. - // -_nx = W boundary, N to S. - // _nx = E boundary, S to N. - // -_nx+1 = NE corner, NW to SE. - // _nx+1 = NW corner, SW to NE. - // -_nx-1 = SE corner, NE to SW. - // _nx-1 = SW corner, SE to NW. - - assert(is_quad_in_chunk(start_location.quad, local)); - assert(is_quad_in_chunk(location.quad, local)); - - // Local variables for faster access. - auto quad = location.quad; - auto forward = location.forward; - auto left = location.left; - auto start_quad = start_location.quad; - auto start_forward = start_location.forward; - auto start_left = start_location.left; - auto pass = local.pass; - double*& points = local.points.current; - - auto start_point = get_boundary_start_point(location); - auto end_point = start_point + forward; - - assert(is_point_in_chunk(start_point, local)); - assert(is_point_in_chunk(end_point, local)); - - auto start_z = Z_LEVEL(start_point); - auto end_z = Z_LEVEL(end_point); - - // Add new point, somewhere along start line. May be at start point of edge if this is a - // boundary start. - point_count++; - if (pass > 0) { - if (start_z == 1) - get_point_xy(start_point, points); - else // start_z != 1 - interp(start_point, end_point, location.is_upper, points); - } - - bool finished = false; - while (true) { - assert(is_quad_in_chunk(quad, local)); - - if (quad == start_quad && forward == start_forward && left == start_left) { - if (start_location.on_boundary && point_count > 1) { - // Polygon closed. - finished = true; - break; - } - } - else if (pass == 0) { - // Clear unwanted start locations. - if (left == _nx) { - if (START_BOUNDARY_S(quad)) { - assert(forward == 1); - _cache[quad] &= ~MASK_START_BOUNDARY_S; - } - } - else if (forward == -_nx) { - if (START_BOUNDARY_W(quad)) { - assert(left == 1); - _cache[quad] &= ~MASK_START_BOUNDARY_W; - } - } - else if (left == -_nx) { - if (START_HOLE_N(quad)) { - assert(forward == -1); - _cache[quad] &= ~MASK_START_HOLE_N; - } - } - else { - switch (EXISTS_ANY_CORNER(quad)) { - case MASK_EXISTS_NE_CORNER: - if (left == _nx+1) { - assert(forward == -_nx+1); - _cache[quad] &= ~MASK_START_CORNER; - } - break; - case MASK_EXISTS_NW_CORNER: - if (forward == _nx+1) { - assert(left == _nx-1); - _cache[quad] &= ~MASK_START_CORNER; - } - break; - case MASK_EXISTS_SE_CORNER: - if (forward == -_nx-1) { - assert(left == -_nx+1); - _cache[quad] &= ~MASK_START_CORNER; - } - break; - case MASK_EXISTS_SW_CORNER: - if (left == -_nx-1) { - assert(forward == _nx-1); - _cache[quad] &= ~MASK_START_CORNER; - } - break; - default: - // Not a corner. - break; - } - } - } - - // Check if need to leave boundary into interior. - if (end_z != 1) { - location.is_upper = (end_z == 2); // Leave via this level. - auto temp = forward; - forward = left; - left = -temp; - break; - } - - // Add end point. - point_count++; - if (pass > 0) { - get_point_xy(end_point, points); - - if (LOOK_N(quad) && _identify_holes && - (left == _nx || left == _nx+1 || forward == _nx+1)) { - assert(BOUNDARY_N(quad-_nx) || EXISTS_NE_CORNER(quad) || EXISTS_NW_CORNER(quad)); - local.look_up_quads.push_back(quad); - } - } - - move_to_next_boundary_edge(quad, forward, left); - - start_point = end_point; - start_z = end_z; - end_point = start_point + forward; - end_z = Z_LEVEL(end_point); - } - - location.quad = quad; - location.forward = forward; - location.left = left; - - return finished; -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::follow_interior( - Location& location, const Location& start_location, ChunkLocal& local, count_t& point_count) -{ - // Adds the start point in each quad visited, but not the end point unless closing the polygon. - // Only need to consider a single level of course. - assert(is_quad_in_chunk(start_location.quad, local)); - assert(is_quad_in_chunk(location.quad, local)); - - // Local variables for faster access. - auto quad = location.quad; - auto forward = location.forward; - auto left = location.left; - auto is_upper = location.is_upper; - auto start_quad = start_location.quad; - auto start_forward = start_location.forward; - auto start_left = start_location.left; - auto pass = local.pass; - double*& points = local.points.current; - - // left direction, and indices of points on entry edge. - bool start_corner_diagonal = false; - auto left_point = get_interior_start_left_point(location, start_corner_diagonal); - auto right_point = left_point - left; - bool want_look_N = _identify_holes && pass > 0; - - bool finished = false; // Whether finished line, i.e. returned to start. - while (true) { - assert(is_quad_in_chunk(quad, local)); - assert(is_point_in_chunk(left_point, local)); - assert(is_point_in_chunk(right_point, local)); - - if (pass > 0) - interp(left_point, right_point, is_upper, points); - point_count++; - - if (quad == start_quad && forward == start_forward && - left == start_left && is_upper == start_location.is_upper && - !start_location.on_boundary && point_count > 1) { - finished = true; // Polygon closed, exit immediately. - break; - } - - // Indices of the opposite points. - auto opposite_left_point = left_point + forward; - auto opposite_right_point = right_point + forward; - bool corner_opposite_is_right = false; // Only used for corners. - - if (start_corner_diagonal) { - // To avoid dealing with diagonal forward and left below, switch to direction 45 degrees - // to left, e.g. NW corner faces west using forward == -1. - corner_opposite_is_right = true; - switch (EXISTS_ANY_CORNER(quad)) { - case MASK_EXISTS_NW_CORNER: - forward = -1; - left = -_nx; - opposite_left_point = opposite_right_point = quad-1; - break; - case MASK_EXISTS_NE_CORNER: - forward = _nx; - left = -1; - opposite_left_point = opposite_right_point = quad; - break; - case MASK_EXISTS_SW_CORNER: - forward = -_nx; - left = 1; - opposite_left_point = opposite_right_point = quad-_nx-1; - break; - default: - assert(EXISTS_SE_CORNER(quad)); - forward = 1; - left = _nx; - opposite_left_point = opposite_right_point = quad-_nx; - break; - } - } - - // z-levels of the opposite points. - ZLevel z_opposite_left = Z_LEVEL(opposite_left_point); - ZLevel z_opposite_right = Z_LEVEL(opposite_right_point); - - Direction direction = Direction::Right; - ZLevel z_test = is_upper ? 2 : 0; - - if (EXISTS_QUAD(quad)) { - if (z_opposite_left == z_test) { - if (z_opposite_right == z_test || MIDDLE_Z_LEVEL(quad) == z_test) - direction = Direction::Left; - } - else if (z_opposite_right == z_test) - direction = Direction::Straight; - } - else if (start_corner_diagonal) { - direction = (z_opposite_left == z_test) ? Direction::Straight : Direction::Right; - } - else { - switch (EXISTS_ANY_CORNER(quad)) { - case MASK_EXISTS_NW_CORNER: - corner_opposite_is_right = (forward == -_nx); - break; - case MASK_EXISTS_NE_CORNER: - corner_opposite_is_right = (forward == -1); - break; - case MASK_EXISTS_SW_CORNER: - corner_opposite_is_right = (forward == 1); - break; - default: - assert(EXISTS_SE_CORNER(quad)); - corner_opposite_is_right = (forward == _nx); - break; - } - - if (corner_opposite_is_right) - direction = (z_opposite_right == z_test) ? Direction::Straight : Direction::Right; - else - direction = (z_opposite_left == z_test) ? Direction::Left : Direction::Straight; - } - - // Clear unwanted start locations. - if (pass == 0 && !(quad == start_quad && forward == start_forward && left == start_left)) { - if (START_E(quad) && forward == -1 && left == -_nx && direction == Direction::Right && - (is_upper ? Z_NE > 0 : Z_NE < 2)) { - _cache[quad] &= ~MASK_START_E; // E high if is_upper else low. - - if (!_filled && quad < start_location.quad) - // Already counted points from here onwards. - break; - } - else if (START_N(quad) && forward == -_nx && left == 1 && - direction == Direction::Left && (is_upper ? Z_NW > 0 : Z_NW < 2)) { - _cache[quad] &= ~MASK_START_N; // E high if is_upper else low. - - if (!_filled && quad < start_location.quad) - // Already counted points from here onwards. - break; - } - } - - // Extra quad_as_tri points. - if (_quad_as_tri && EXISTS_QUAD(quad)) { - if (pass == 0) { - switch (direction) { - case Direction::Left: - point_count += (LEFT_OF_MIDDLE(quad, is_upper) ? 1 : 3); - break; - case Direction::Right: - point_count += (LEFT_OF_MIDDLE(quad, is_upper) ? 3 : 1); - break; - case Direction::Straight: - point_count += 2; - break; - } - } - else { // pass == 1 - auto mid_x = get_middle_x(quad); - auto mid_y = get_middle_y(quad); - auto mid_z = calc_middle_z(quad); - - switch (direction) { - case Direction::Left: - if (LEFT_OF_MIDDLE(quad, is_upper)) { - interp(left_point, mid_x, mid_y, mid_z, is_upper, points); - point_count++; - } - else { - interp(right_point, mid_x, mid_y, mid_z, is_upper, points); - interp(opposite_right_point, mid_x, mid_y, mid_z, is_upper, points); - interp(opposite_left_point, mid_x, mid_y, mid_z, is_upper, points); - point_count += 3; - } - break; - case Direction::Right: - if (LEFT_OF_MIDDLE(quad, is_upper)) { - interp(left_point, mid_x, mid_y, mid_z, is_upper, points); - interp(opposite_left_point, mid_x, mid_y, mid_z, is_upper, points); - interp(opposite_right_point, mid_x, mid_y, mid_z, is_upper, points); - point_count += 3; - } - else { - interp(right_point, mid_x, mid_y, mid_z, is_upper, points); - point_count++; - } - break; - case Direction::Straight: - if (LEFT_OF_MIDDLE(quad, is_upper)) { - interp(left_point, mid_x, mid_y, mid_z, is_upper, points); - interp(opposite_left_point, mid_x, mid_y, mid_z, is_upper, points); - } - else { - interp(right_point, mid_x, mid_y, mid_z, is_upper, points); - interp(opposite_right_point, mid_x, mid_y, mid_z, is_upper, points); - } - point_count += 2; - break; - } - } - } - - bool reached_boundary = false; - - // Determine entry edge and left and right points of next quad. - // Do not update quad index yet. - switch (direction) { - case Direction::Left: { - auto temp = forward; - forward = left; - left = -temp; - // left_point unchanged. - right_point = opposite_left_point; - break; - } - case Direction::Right: { - auto temp = forward; - forward = -left; - left = temp; - left_point = opposite_right_point; - // right_point unchanged. - break; - } - case Direction::Straight: { - if (EXISTS_QUAD(quad)) { // Straight on in quad. - // forward and left stay the same. - left_point = opposite_left_point; - right_point = opposite_right_point; - } - else if (start_corner_diagonal) { // Straight on diagonal start corner. - // left point unchanged. - right_point = opposite_right_point; - } - else { // Straight on in a corner reaches boundary. - assert(EXISTS_ANY_CORNER(quad)); - reached_boundary = true; - - if (corner_opposite_is_right) { - // left_point unchanged. - right_point = opposite_right_point; - } - else { - left_point = opposite_left_point; - // right_point unchanged. - } - - // Set forward and left for correct exit along boundary. - switch (EXISTS_ANY_CORNER(quad)) { - case MASK_EXISTS_NW_CORNER: - forward = _nx+1; - left = _nx-1; - break; - case MASK_EXISTS_NE_CORNER: - forward = -_nx+1; - left = _nx+1; - break; - case MASK_EXISTS_SW_CORNER: - forward = _nx-1; - left = -_nx-1; - break; - default: - assert(EXISTS_SE_CORNER(quad)); - forward = -_nx-1; - left = -_nx+1; - break; - } - } - break; - } - } - - if (want_look_N && LOOK_N(quad) && forward == 1) { - // Only consider look_N if pass across E edge of this quad. - // Care needed if both look_N and look_S set in quad because this line corresponds to - // only one of them, so want to ignore the look_N if it is the other line otherwise it - // will be double counted. - if (!LOOK_S(quad) || (is_upper ? Z_NE < 2 : Z_NE > 0)) - local.look_up_quads.push_back(quad); - } - - // Check if reached NSEW boundary; already checked and noted if reached corner boundary. - if (!reached_boundary) { - if (forward > 0) - reached_boundary = (forward == 1 ? BOUNDARY_E(quad) : BOUNDARY_N(quad)); - else // forward < 0 - reached_boundary = (forward == -1 ? BOUNDARY_W(quad) : BOUNDARY_S(quad)); - - if (reached_boundary) { - auto temp = forward; - forward = left; - left = -temp; - } - } - - // If reached a boundary, return. - if (reached_boundary) { - if (!_filled) { - point_count++; - if (pass > 0) - interp(left_point, right_point, false, points); - } - break; - } - - quad += forward; - start_corner_diagonal = false; - } - - location.quad = quad; - location.forward = forward; - location.left = left; - location.is_upper = is_upper; - - return finished; -} - -template <typename Derived> -index_t BaseContourGenerator<Derived>::get_boundary_start_point(const Location& location) const -{ - auto quad = location.quad; - auto forward = location.forward; - auto left = location.left; - index_t start_point = -1; - - if (forward > 0) { - if (forward == _nx) { - assert(left == -1); - start_point = quad-_nx; - } - else if (left == _nx) { - assert(forward == 1); - start_point = quad-_nx-1; - } - else if (EXISTS_SW_CORNER(quad)) { - assert(forward == _nx-1 && left == -_nx-1); - start_point = quad-_nx; - } - else { - assert(EXISTS_NW_CORNER(quad) && forward == _nx+1 && left == _nx-1); - start_point = quad-_nx-1; - } - } - else { // forward < 0 - if (forward == -_nx) { - assert(left == 1); - start_point = quad-1; - } - else if (left == -_nx) { - assert(forward == -1); - start_point = quad; - } - else if (EXISTS_NE_CORNER(quad)) { - assert(forward == -_nx+1 && left == _nx+1); - start_point = quad-1; - } - else { - assert(EXISTS_SE_CORNER(quad) && forward == -_nx-1 && left == -_nx+1); - start_point = quad; - } - } - return start_point; -} - -template <typename Derived> -py::tuple BaseContourGenerator<Derived>::get_chunk_count() const -{ - return py::make_tuple(_ny_chunks, _nx_chunks); -} - -template <typename Derived> -void BaseContourGenerator<Derived>::get_chunk_limits(index_t chunk, ChunkLocal& local) const -{ - assert(chunk >= 0 && chunk < _n_chunks && "chunk index out of bounds"); - - local.chunk = chunk; - - auto ichunk = chunk % _nx_chunks; - auto jchunk = chunk / _nx_chunks; - - local.istart = ichunk*_x_chunk_size + 1; - local.iend = (ichunk < _nx_chunks-1 ? (ichunk+1)*_x_chunk_size : _nx-1); - - local.jstart = jchunk*_y_chunk_size + 1; - local.jend = (jchunk < _ny_chunks-1 ? (jchunk+1)*_y_chunk_size : _ny-1); -} - -template <typename Derived> -py::tuple BaseContourGenerator<Derived>::get_chunk_size() const -{ - return py::make_tuple(_y_chunk_size, _x_chunk_size); -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::get_corner_mask() const -{ - return _corner_mask; -} - -template <typename Derived> -FillType BaseContourGenerator<Derived>::get_fill_type() const -{ - return _fill_type; -} - -template <typename Derived> -index_t BaseContourGenerator<Derived>::get_interior_start_left_point( - const Location& location, bool& start_corner_diagonal) const -{ - auto quad = location.quad; - auto forward = location.forward; - auto left = location.left; - index_t left_point = -1; - - if (forward > 0) { - if (forward == _nx) { - assert(left == -1); - left_point = quad-_nx-1; - } - else if (left == _nx) { - assert(forward == 1); - left_point = quad-1; - } - else if (EXISTS_NW_CORNER(quad)) { - assert(forward == _nx-1 && left == -_nx-1); - left_point = quad-_nx-1; - start_corner_diagonal = true; - } - else { - assert(EXISTS_NE_CORNER(quad) && forward == _nx+1 && left == _nx-1); - left_point = quad-1; - start_corner_diagonal = true; - } - } - else { // forward < 0 - if (forward == -_nx) { - assert(left == 1); - left_point = quad; - } - else if (left == -_nx) { - assert(forward == -1); - left_point = quad-_nx; - } - else if (EXISTS_SW_CORNER(quad)) { - assert(forward == -_nx-1 && left == -_nx+1); - left_point = quad-_nx; - start_corner_diagonal = true; - } - else { - assert(EXISTS_SE_CORNER(quad) && forward == -_nx+1 && left == _nx+1); - left_point = quad; - start_corner_diagonal = true; - } - } - return left_point; -} - -template <typename Derived> -double BaseContourGenerator<Derived>::get_interp_fraction(double z0, double z1, double level) const -{ - switch (_z_interp) { - case ZInterp::Log: - // Equivalent to - // (log(z1) - log(level)) / (log(z1) - log(z0)) - // Same result obtained regardless of logarithm base. - return log(z1/level) / log(z1/z0); - default: // ZInterp::Linear - return (z1 - level) / (z1 - z0); - } -} - -template <typename Derived> -LineType BaseContourGenerator<Derived>::get_line_type() const -{ - return _line_type; -} - -template <typename Derived> -double BaseContourGenerator<Derived>::get_middle_x(index_t quad) const -{ - return 0.25*(get_point_x(POINT_SW) + get_point_x(POINT_SE) + - get_point_x(POINT_NW) + get_point_x(POINT_NE)); -} - -template <typename Derived> -double BaseContourGenerator<Derived>::get_middle_y(index_t quad) const -{ - return 0.25*(get_point_y(POINT_SW) + get_point_y(POINT_SE) + - get_point_y(POINT_NW) + get_point_y(POINT_NE)); -} - -template <typename Derived> -index_t BaseContourGenerator<Derived>::get_n_chunks() const -{ - return _n_chunks; -} - -template <typename Derived> -void BaseContourGenerator<Derived>::get_point_xy(index_t point, double*& points) const -{ - assert(point >= 0 && point < _n && "point index out of bounds"); - *points++ = _xptr[point]; - *points++ = _yptr[point]; -} - -template <typename Derived> -double BaseContourGenerator<Derived>::get_point_x(index_t point) const -{ - assert(point >= 0 && point < _n && "point index out of bounds"); - return _xptr[point]; -} - -template <typename Derived> -double BaseContourGenerator<Derived>::get_point_y(index_t point) const -{ - assert(point >= 0 && point < _n && "point index out of bounds"); - return _yptr[point]; -} - -template <typename Derived> -double BaseContourGenerator<Derived>::get_point_z(index_t point) const -{ - assert(point >= 0 && point < _n && "point index out of bounds"); - return _zptr[point]; -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::get_quad_as_tri() const -{ - return _quad_as_tri; -} - -template <typename Derived> -ZInterp BaseContourGenerator<Derived>::get_z_interp() const -{ - return _z_interp; -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::has_direct_line_offsets() const -{ - return _direct_line_offsets; -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::has_direct_outer_offsets() const -{ - return _direct_outer_offsets; -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::has_direct_points() const -{ - return _direct_points; -} - -template <typename Derived> -void BaseContourGenerator<Derived>::init_cache_grid(const MaskArray& mask) -{ - index_t i, j, quad; - if (mask.ndim() == 0) { - // No mask, easy to calculate quad existence and boundaries together. - for (j = 0, quad = 0; j < _ny; ++j) { - for (i = 0; i < _nx; ++i, ++quad) { - _cache[quad] = 0; - - if (i > 0 && j > 0) - _cache[quad] |= MASK_EXISTS_QUAD; - - if ((i % _x_chunk_size == 0 || i == _nx-1) && j > 0) - _cache[quad] |= MASK_BOUNDARY_E; - - if ((j % _y_chunk_size == 0 || j == _ny-1) && i > 0) - _cache[quad] |= MASK_BOUNDARY_N; - } - } - } - else { - // Could maybe speed this up and just have a single pass. - // Care would be needed with lookback of course. - const bool* mask_ptr = mask.data(); - - // Have mask so use two stages. - // Stage 1, determine if quads/corners exist. - quad = 0; - for (j = 0; j < _ny; ++j) { - for (i = 0; i < _nx; ++i, ++quad) { - _cache[quad] = 0; - - if (i > 0 && j > 0) { - unsigned int config = (mask_ptr[POINT_NW] << 3) | - (mask_ptr[POINT_NE] << 2) | - (mask_ptr[POINT_SW] << 1) | - (mask_ptr[POINT_SE] << 0); - if (_corner_mask) { - switch (config) { - case 0: _cache[quad] = MASK_EXISTS_QUAD; break; - case 1: _cache[quad] = MASK_EXISTS_NW_CORNER; break; - case 2: _cache[quad] = MASK_EXISTS_NE_CORNER; break; - case 4: _cache[quad] = MASK_EXISTS_SW_CORNER; break; - case 8: _cache[quad] = MASK_EXISTS_SE_CORNER; break; - default: - // Do nothing, quad is masked out. - break; - } - } - else if (config == 0) - _cache[quad] = MASK_EXISTS_QUAD; - } - } - } - - // Stage 2, calculate N and E boundaries. - quad = 0; - for (j = 0; j < _ny; ++j) { - bool j_chunk_boundary = j % _y_chunk_size == 0; - - for (i = 0; i < _nx; ++i, ++quad) { - bool i_chunk_boundary = i % _x_chunk_size == 0; - - if (_corner_mask) { - bool exists_E_edge = EXISTS_E_EDGE(quad); - bool E_exists_W_edge = (i < _nx-1 && EXISTS_W_EDGE(quad+1)); - bool exists_N_edge = EXISTS_N_EDGE(quad); - bool N_exists_S_edge = (j < _ny-1 && EXISTS_S_EDGE(quad+_nx)); - - if (exists_E_edge != E_exists_W_edge || - (i_chunk_boundary && exists_E_edge && E_exists_W_edge)) - _cache[quad] |= MASK_BOUNDARY_E; - - if (exists_N_edge != N_exists_S_edge || - (j_chunk_boundary && exists_N_edge && N_exists_S_edge)) - _cache[quad] |= MASK_BOUNDARY_N; - } - else { - bool E_exists_quad = (i < _nx-1 && EXISTS_QUAD(quad+1)); - bool N_exists_quad = (j < _ny-1 && EXISTS_QUAD(quad+_nx)); - bool exists = EXISTS_QUAD(quad); - - if (exists != E_exists_quad || (i_chunk_boundary && exists && E_exists_quad)) - _cache[quad] |= MASK_BOUNDARY_E; - - if (exists != N_exists_quad || (j_chunk_boundary && exists && N_exists_quad)) - _cache[quad] |= MASK_BOUNDARY_N; - } - } - } - } -} - -template <typename Derived> -void BaseContourGenerator<Derived>::init_cache_levels_and_starts(const ChunkLocal* local) -{ - bool ordered_chunks = (local == nullptr); - - // This function initialises the cache z-levels and starts for either a single chunk or the - // whole domain. If a single chunk, only the quads contained in the chunk are calculated and - // this includes the z-levels of the points that on the NE corners of those quads. In addition, - // chunks that are on the W (starting at i=1) also calculate the most westerly points (i=0), - // and similarly chunks that are on the S (starting at j=1) also calculate the most southerly - // points (j=0). Non W/S chunks do not do this as their neighboring chunks to the W/S are - // responsible for it. If ordered_chunks is true then those W/S points will already have had - // their cache items set so that their z-levels can be read from the cache as usual. But if - // ordered_chunks is false then we cannot rely upon those neighboring W/S points having their - // cache items already set and so must temporarily calculate those z-levels rather than reading - // the cache. - - constexpr CacheItem keep_mask = (MASK_EXISTS_ANY | MASK_BOUNDARY_N | MASK_BOUNDARY_E); - - index_t istart, iend, jstart, jend; // Loop indices. - index_t chunk_istart; // Actual start i-index of chunk. - - if (local != nullptr) { - chunk_istart = local->istart; - istart = chunk_istart > 1 ? chunk_istart : 0; - iend = local->iend; - jstart = local->jstart > 1 ? local->jstart : 0; - jend = local->jend; - } - else { - chunk_istart = 1; - istart = 0; - iend = _nx-1; - jstart = 0; - jend = _ny-1; - } - - index_t j_final_start = jstart - 1; - bool calc_W_z_level = (!ordered_chunks && istart == chunk_istart); - - for (index_t j = jstart; j <= jend; ++j) { - index_t quad = istart + j*_nx; - const double* z_ptr = _zptr + quad; - bool start_in_row = false; - bool calc_S_z_level = (!ordered_chunks && j == jstart); - - // z-level of NW point not needed if i == 0. - ZLevel z_nw = (istart == 0) ? 0 : (calc_W_z_level ? z_to_zlevel(*(z_ptr-1)) : Z_NW); - - // z-level of SW point not needed if i == 0 or j == 0. - ZLevel z_sw = (istart == 0 || j == 0) ? 0 : - ((calc_W_z_level || calc_S_z_level) ? z_to_zlevel(*(z_ptr-_nx-1)) : Z_SW); - - for (index_t i = istart; i <= iend; ++i, ++quad, ++z_ptr) { - // z-level of SE point not needed if j == 0. - ZLevel z_se = (j == 0) ? 0 : (calc_S_z_level ? z_to_zlevel(*(z_ptr-_nx)) : Z_SE); - - _cache[quad] &= keep_mask; - - // Calculate and cache z-level of NE point. - ZLevel z_ne = z_to_zlevel(*z_ptr); - _cache[quad] |= z_ne; - - switch (EXISTS_ANY(quad)) { - case MASK_EXISTS_QUAD: - if (_filled) { - switch ((z_nw << 6) | (z_ne << 4) | (z_sw << 2) | z_se) { // config - case 1: // 0001 - case 2: // 0002 - case 17: // 0101 - case 18: // 0102 - case 34: // 0202 - case 68: // 1010 - case 102: // 1212 - case 136: // 2020 - case 152: // 2120 - case 153: // 2121 - case 168: // 2220 - case 169: // 2221 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_S; - start_in_row = true; - } - break; - case 4: // 0010 - case 5: // 0011 - case 6: // 0012 - case 8: // 0020 - case 9: // 0021 - case 21: // 0111 - case 22: // 0112 - case 25: // 0121 - case 38: // 0212 - case 72: // 1020 - case 98: // 1202 - case 132: // 2010 - case 145: // 2101 - case 148: // 2110 - case 149: // 2111 - case 161: // 2201 - case 162: // 2202 - case 164: // 2210 - case 165: // 2211 - case 166: // 2212 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - start_in_row |= ANY_START(quad); - break; - case 10: // 0022 - case 26: // 0122 - case 42: // 0222 - case 64: // 1000 - case 106: // 1222 - case 128: // 2000 - case 144: // 2100 - case 160: // 2200 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_W(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_W; - start_in_row = true; - } - break; - case 16: // 0100 - case 154: // 2122 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - _cache[quad] |= MASK_START_N; - start_in_row = true; - break; - case 20: // 0110 - case 24: // 0120 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (MIDDLE_Z_LEVEL(quad) == 0) _cache[quad] |= MASK_START_N; - start_in_row |= ANY_START(quad); - break; - case 32: // 0200 - case 138: // 2022 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - _cache[quad] |= MASK_START_E; - _cache[quad] |= MASK_START_N; - start_in_row = true; - break; - case 33: // 0201 - case 69: // 1011 - case 70: // 1012 - case 100: // 1210 - case 101: // 1211 - case 137: // 2021 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - _cache[quad] |= MASK_START_E; - start_in_row = true; - break; - case 36: // 0210 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (MIDDLE_Z_LEVEL(quad) == 0) _cache[quad] |= MASK_START_N; - _cache[quad] |= MASK_START_E; - start_in_row = true; - break; - case 37: // 0211 - case 73: // 1021 - case 97: // 1201 - case 133: // 2011 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - _cache[quad] |= MASK_START_E; - start_in_row = true; - break; - case 40: // 0220 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (MIDDLE_Z_LEVEL(quad) < 2) _cache[quad] |= MASK_START_E; - if (MIDDLE_Z_LEVEL(quad) == 0) _cache[quad] |= MASK_START_N; - start_in_row |= ANY_START(quad); - break; - case 41: // 0221 - case 104: // 1220 - case 105: // 1221 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (MIDDLE_Z_LEVEL(quad) < 2) _cache[quad] |= MASK_START_E; - start_in_row |= ANY_START(quad); - break; - case 65: // 1001 - case 66: // 1002 - case 129: // 2001 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (MIDDLE_Z_LEVEL(quad) > 0) _cache[quad] |= MASK_START_E; - start_in_row |= ANY_START(quad); - break; - case 74: // 1022 - case 96: // 1200 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - _cache[quad] |= MASK_START_E; - start_in_row = true; - break; - case 80: // 1100 - case 90: // 1122 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (BOUNDARY_N(quad) && !START_HOLE_N(quad-1) && - j % _y_chunk_size > 0 && j != _ny-1 && i % _x_chunk_size > 1) - _cache[quad] |= MASK_START_HOLE_N; - start_in_row |= ANY_START(quad); - break; - case 81: // 1101 - case 82: // 1102 - case 88: // 1120 - case 89: // 1121 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (BOUNDARY_N(quad) && !START_HOLE_N(quad-1) && - j % _y_chunk_size > 0 && j != _ny-1 && i % _x_chunk_size > 1) - _cache[quad] |= MASK_START_HOLE_N; - start_in_row |= ANY_START(quad); - break; - case 84: // 1110 - case 85: // 1111 - case 86: // 1112 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_N(quad) && !START_HOLE_N(quad-1) && - j % _y_chunk_size > 0 && j != _ny-1 && i % _x_chunk_size > 1) - _cache[quad] |= MASK_START_HOLE_N; - start_in_row |= ANY_START(quad); - break; - case 130: // 2002 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (MIDDLE_Z_LEVEL(quad) > 0) _cache[quad] |= MASK_START_E; - if (MIDDLE_Z_LEVEL(quad) == 2) _cache[quad] |= MASK_START_N; - start_in_row |= ANY_START(quad); - break; - case 134: // 2012 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (MIDDLE_Z_LEVEL(quad) == 2) _cache[quad] |= MASK_START_N; - _cache[quad] |= MASK_START_E; - start_in_row = true; - break; - case 146: // 2102 - case 150: // 2112 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (MIDDLE_Z_LEVEL(quad) == 2) _cache[quad] |= MASK_START_N; - start_in_row |= ANY_START(quad); - break; - } - } - else { // !_filled quad - switch ((z_nw << 3) | (z_ne << 2) | (z_sw << 1) | z_se) { // config - case 1: // 0001 - case 3: // 0011 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_E(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_E; - start_in_row = true; - } - break; - case 2: // 0010 - case 10: // 1010 - case 14: // 1110 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_S(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_S; - start_in_row = true; - } - break; - case 4: // 0100 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_N(quad)) - _cache[quad] |= MASK_START_BOUNDARY_N; - else if (!BOUNDARY_E(quad)) - _cache[quad] |= MASK_START_N; - start_in_row |= ANY_START(quad); - break; - case 5: // 0101 - case 7: // 0111 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_N(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_N; - start_in_row = true; - } - break; - case 6: // 0110 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_N(quad)) - _cache[quad] |= MASK_START_BOUNDARY_N; - else if (!BOUNDARY_E(quad) && MIDDLE_Z_LEVEL(quad) == 0) - _cache[quad] |= MASK_START_N; - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - start_in_row |= ANY_START(quad); - break; - case 8: // 1000 - case 12: // 1100 - case 13: // 1101 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_W(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_W; - start_in_row = true; - } - break; - case 9: // 1001 - calc_and_set_middle_z_level(quad); - if (BOUNDARY_E(quad)) - _cache[quad] |= MASK_START_BOUNDARY_E; - else if (!BOUNDARY_N(quad) && MIDDLE_Z_LEVEL(quad) == 1) - _cache[quad] |= MASK_START_E; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - start_in_row |= ANY_START(quad); - break; - case 11: // 1011 - if (_quad_as_tri) calc_and_set_middle_z_level(quad); - if (BOUNDARY_E(quad)) - _cache[quad] |= MASK_START_BOUNDARY_E; - else if (!BOUNDARY_N(quad)) - _cache[quad] |= MASK_START_E; - start_in_row |= ANY_START(quad); - break; - } - } - break; - case MASK_EXISTS_NW_CORNER: - if (_filled) { - switch ((z_nw << 4) | (z_ne << 2) | z_sw) { // config - case 1: // 001 - case 5: // 011 - case 9: // 021 - case 10: // 022 - case 16: // 100 - case 17: // 101 - case 25: // 121 - case 26: // 122 - case 32: // 200 - case 33: // 201 - case 37: // 211 - case 41: // 221 - if (BOUNDARY_W(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_W; - start_in_row = true; - } - break; - case 2: // 002 - case 6: // 012 - case 18: // 102 - case 24: // 120 - case 36: // 210 - case 40: // 220 - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - case 4: // 010 - case 8: // 020 - case 34: // 202 - case 38: // 212 - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - case 20: // 110 - case 22: // 112 - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (BOUNDARY_N(quad) && !START_HOLE_N(quad-1) && - j % _y_chunk_size > 0 && j != _ny-1 && i % _x_chunk_size > 1) - _cache[quad] |= MASK_START_HOLE_N; - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - case 21: // 111 - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - if (BOUNDARY_N(quad) && !START_HOLE_N(quad-1) && - j % _y_chunk_size > 0 && j != _ny-1 && i % _x_chunk_size > 1) - _cache[quad] |= MASK_START_HOLE_N; - start_in_row |= ANY_START(quad); - break; - } - } - else { // !_filled NW corner. - switch ((z_nw << 2) | (z_ne << 1) | z_sw) { // config - case 1: // 001 - case 5: // 101 - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - case 2: // 010 - case 3: // 011 - if (BOUNDARY_N(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_N; - start_in_row = true; - } - break; - case 4: // 100 - case 6: // 110 - if (BOUNDARY_W(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_W; - start_in_row = true; - } - break; - } - } - break; - case MASK_EXISTS_NE_CORNER: - if (_filled) { - switch ((z_nw << 4) | (z_ne << 2) | z_se) { // config - case 1: // 001 - case 2: // 002 - case 5: // 011 - case 6: // 012 - case 10: // 022 - case 16: // 100 - case 26: // 122 - case 32: // 200 - case 36: // 210 - case 37: // 211 - case 40: // 220 - case 41: // 221 - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - case 4: // 010 - case 38: // 212 - _cache[quad] |= MASK_START_N; - start_in_row = true; - break; - case 8: // 020 - case 34: // 202 - _cache[quad] |= MASK_START_E; - _cache[quad] |= MASK_START_N; - start_in_row = true; - break; - case 9: // 021 - case 17: // 101 - case 18: // 102 - case 24: // 120 - case 25: // 121 - case 33: // 201 - _cache[quad] |= MASK_START_CORNER; - _cache[quad] |= MASK_START_E; - start_in_row = true; - break; - case 20: // 110 - case 21: // 111 - case 22: // 112 - if (BOUNDARY_N(quad) && !START_HOLE_N(quad-1) && - j % _y_chunk_size > 0 && j != _ny-1 && i % _x_chunk_size > 1) - _cache[quad] |= MASK_START_HOLE_N; - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - } - } - else { // !_filled NE corner. - switch ((z_nw << 2) | (z_ne << 1) | z_se) { // config - case 1: // 001 - if (BOUNDARY_E(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_E; - start_in_row = true; - } - break; - case 2: // 010 - if (BOUNDARY_N(quad)) - _cache[quad] |= MASK_START_BOUNDARY_N; - else if (!BOUNDARY_E(quad)) - _cache[quad] |= MASK_START_N; - start_in_row |= ANY_START(quad); - break; - case 3: // 011 - if (BOUNDARY_N(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_N; - start_in_row = true; - } - break; - case 4: // 100 - case 6: // 110 - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - case 5: // 101 - if (BOUNDARY_E(quad)) - _cache[quad] |= MASK_START_BOUNDARY_E; - else if (!BOUNDARY_N(quad)) - _cache[quad] |= MASK_START_E; - start_in_row |= ANY_START(quad); - break; - } - } - break; - case MASK_EXISTS_SW_CORNER: - if (_filled) { - switch ((z_nw << 4) | (z_sw << 2) | z_se) { // config - case 1: // 001 - case 2: // 002 - case 40: // 220 - case 41: // 221 - if (BOUNDARY_S(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_S; - start_in_row = true; - } - break; - case 4: // 010 - case 5: // 011 - case 6: // 012 - case 8: // 020 - case 9: // 021 - case 18: // 102 - case 24: // 120 - case 33: // 201 - case 34: // 202 - case 36: // 210 - case 37: // 211 - case 38: // 212 - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - start_in_row |= ANY_START(quad); - break; - case 10: // 022 - case 16: // 100 - case 26: // 122 - case 32: // 200 - if (BOUNDARY_W(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_W; - start_in_row = true; - } - break; - case 17: // 101 - case 25: // 121 - if (BOUNDARY_S(quad)) _cache[quad] |= MASK_START_BOUNDARY_S; - if (BOUNDARY_W(quad)) _cache[quad] |= MASK_START_BOUNDARY_W; - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - case 20: // 110 - case 21: // 111 - case 22: // 112 - if (BOUNDARY_S(quad)) - _cache[quad] |= MASK_START_BOUNDARY_S; - else - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - } - } - else { // !_filled SW corner. - switch ((z_nw << 2) | (z_sw << 1) | z_se) { // config - case 1: // 001 - case 3: // 011 - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - case 2: // 010 - case 6: // 110 - if (BOUNDARY_S(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_S; - start_in_row = true; - } - break; - case 4: // 100 - case 5: // 101 - if (BOUNDARY_W(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_W; - start_in_row = true; - } - break; - } - } - break; - case MASK_EXISTS_SE_CORNER: - if (_filled) { - switch ((z_ne << 4) | (z_sw << 2) | z_se) { // config - case 1: // 001 - case 2: // 002 - case 4: // 010 - case 5: // 011 - case 6: // 012 - case 8: // 020 - case 9: // 021 - case 17: // 101 - case 18: // 102 - case 20: // 110 - case 21: // 111 - case 22: // 112 - case 24: // 120 - case 25: // 121 - case 33: // 201 - case 34: // 202 - case 36: // 210 - case 37: // 211 - case 38: // 212 - case 40: // 220 - case 41: // 221 - if (BOUNDARY_S(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_S; - start_in_row = true; - } - break; - case 10: // 022 - case 16: // 100 - case 26: // 122 - case 32: // 200 - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - } - } - else { // !_filled SE corner. - switch ((z_ne << 2) | (z_sw << 1) | z_se) { // config - case 1: // 001 - case 3: // 011 - if (BOUNDARY_E(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_E; - start_in_row = true; - } - break; - case 2: // 010 - case 6: // 110 - if (BOUNDARY_S(quad)) { - _cache[quad] |= MASK_START_BOUNDARY_S; - start_in_row = true; - } - break; - case 4: // 100 - case 5: // 101 - _cache[quad] |= MASK_START_CORNER; - start_in_row = true; - break; - } - } - break; - } - - z_nw = z_ne; - z_sw = z_se; - } // i-loop. - - if (start_in_row) - j_final_start = j; - else if (j > 0) - _cache[chunk_istart + j*_nx] |= MASK_NO_STARTS_IN_ROW; - } // j-loop. - - if (j_final_start < jend) - _cache[chunk_istart + (j_final_start+1)*_nx] |= MASK_NO_MORE_STARTS; -} - -template <typename Derived> -void BaseContourGenerator<Derived>::interp( - index_t point0, index_t point1, bool is_upper, double*& points) const -{ - auto frac = get_interp_fraction( - get_point_z(point0), get_point_z(point1), is_upper ? _upper_level : _lower_level); - - assert(frac >= 0.0 && frac <= 1.0 && "Interp fraction out of bounds"); - - *points++ = get_point_x(point0)*frac + get_point_x(point1)*(1.0 - frac); - *points++ = get_point_y(point0)*frac + get_point_y(point1)*(1.0 - frac); -} - -template <typename Derived> -void BaseContourGenerator<Derived>::interp( - index_t point0, double x1, double y1, double z1, bool is_upper, double*& points) const -{ - auto frac = get_interp_fraction( - get_point_z(point0), z1, is_upper ? _upper_level : _lower_level); - - assert(frac >= 0.0 && frac <= 1.0 && "Interp fraction out of bounds"); - - *points++ = get_point_x(point0)*frac + x1*(1.0 - frac); - *points++ = get_point_y(point0)*frac + y1*(1.0 - frac); -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::is_point_in_chunk(index_t point, const ChunkLocal& local) const -{ - return is_quad_in_bounds(point, local.istart-1, local.iend, local.jstart-1, local.jend); -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::is_quad_in_bounds( - index_t quad, index_t istart, index_t iend, index_t jstart, index_t jend) const -{ - return (quad % _nx >= istart && quad % _nx <= iend && - quad / _nx >= jstart && quad / _nx <= jend); -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::is_quad_in_chunk(index_t quad, const ChunkLocal& local) const -{ - return is_quad_in_bounds(quad, local.istart, local.iend, local.jstart, local.jend); -} - -template <typename Derived> -void BaseContourGenerator<Derived>::line(const Location& start_location, ChunkLocal& local) -{ - // start_location.on_boundary indicates starts (and therefore also finishes) - - assert(is_quad_in_chunk(start_location.quad, local)); - - Location location = start_location; - count_t point_count = 0; - - // Insert nan if required before start of new line. - if (_nan_separated and local.pass > 0 && local.line_count > 0) { - *local.points.current++ = Util::nan; - *local.points.current++ = Util::nan; - } - - // finished == true indicates closed line loop. - bool finished = follow_interior(location, start_location, local, point_count); - - if (local.pass > 0) { - assert(local.line_offsets.current == local.line_offsets.start + local.line_count); - *local.line_offsets.current++ = local.total_point_count; - } - - if (local.pass == 0 && !start_location.on_boundary && !finished) - // An internal start that isn't a line loop is part of a line strip that starts on a - // boundary and will be traced later. Do not count it as a valid start in pass 0 and remove - // the first point or it will be duplicated by the correct boundary-started line later. - point_count--; - else - local.line_count++; - - local.total_point_count += point_count; -} - -template <typename Derived> -py::sequence BaseContourGenerator<Derived>::lines(double level) -{ - pre_lines(); - - _lower_level = _upper_level = level; - return march_wrapper(); -} - -template <typename Derived> -void BaseContourGenerator<Derived>::march_chunk( - ChunkLocal& local, std::vector<py::list>& return_lists) -{ - for (local.pass = 0; local.pass < 2; ++local.pass) { - bool ignore_holes = (_identify_holes && local.pass == 1); - - index_t j_final_start = local.jstart; - for (index_t j = local.jstart; j <= local.jend; ++j) { - index_t quad = local.istart + j*_nx; - - if (NO_MORE_STARTS(quad)) - break; - - if (NO_STARTS_IN_ROW(quad)) - continue; - - // Want to count number of starts in this row, so store how many starts at start of row. - auto prev_start_count = - (_identify_holes ? local.line_count - local.hole_count : local.line_count); - - for (index_t i = local.istart; i <= local.iend; ++i, ++quad) { - if (!ANY_START(quad)) - continue; - - assert(EXISTS_ANY(quad)); - - if (_filled) { - if (START_BOUNDARY_S(quad)) - closed_line_wrapper(Location(quad, 1, _nx, Z_SW == 2, true), Outer, local); - - if (START_BOUNDARY_W(quad)) - closed_line_wrapper(Location(quad, -_nx, 1, Z_NW == 2, true), Outer, local); - - if (START_CORNER(quad)) { - switch (EXISTS_ANY_CORNER(quad)) { - case MASK_EXISTS_NE_CORNER: - closed_line_wrapper( - Location(quad, -_nx+1, _nx+1, Z_NW == 2, true), Outer, local); - break; - case MASK_EXISTS_NW_CORNER: - closed_line_wrapper( - Location(quad, _nx+1, _nx-1, Z_SW == 2, true), Outer, local); - break; - case MASK_EXISTS_SE_CORNER: - closed_line_wrapper( - Location(quad, -_nx-1, -_nx+1, Z_NE == 2, true), Outer, local); - break; - default: - assert(EXISTS_SW_CORNER(quad)); - if (!ignore_holes) - closed_line_wrapper( - Location(quad, _nx-1, -_nx-1, false, true), Hole, local); - break; - } - } - - if (START_N(quad)) - closed_line_wrapper(Location(quad, -_nx, 1, Z_NW > 0, false), Outer, local); - - if (ignore_holes) - continue; - - if (START_E(quad)) - closed_line_wrapper(Location(quad, -1, -_nx, Z_NE > 0, false), Hole, local); - - if (START_HOLE_N(quad)) - closed_line_wrapper(Location(quad, -1, -_nx, false, true), Hole, local); - } - else { // !_filled - if (START_BOUNDARY_S(quad)) - line(Location(quad, _nx, -1, false, true), local); - - if (START_BOUNDARY_W(quad)) - line(Location(quad, 1, _nx, false, true), local); - - if (START_BOUNDARY_E(quad)) - line(Location(quad, -1, -_nx, false, true), local); - - if (START_BOUNDARY_N(quad)) - line(Location(quad, -_nx, 1, false, true), local); - - if (START_E(quad)) - line(Location(quad, -1, -_nx, false, false), local); - - if (START_N(quad)) - line(Location(quad, -_nx, 1, false, false), local); - - if (START_CORNER(quad)) { - index_t forward, left; - switch (EXISTS_ANY_CORNER(quad)) { - case MASK_EXISTS_NE_CORNER: - forward = _nx+1; - left = _nx-1; - break; - case MASK_EXISTS_NW_CORNER: - forward = _nx-1; - left = -_nx-1; - break; - case MASK_EXISTS_SE_CORNER: - forward = -_nx+1; - left = _nx+1; - break; - default: - assert(EXISTS_SW_CORNER(quad)); - forward = -_nx-1; - left = -_nx+1; - break; - } - line(Location(quad, forward, left, false, true), local); - } - } // _filled - } // i - - // Number of starts at end of row. - auto start_count = - (_identify_holes ? local.line_count - local.hole_count : local.line_count); - if (start_count > prev_start_count) - j_final_start = j; - else - _cache[local.istart + j*_nx] |= MASK_NO_STARTS_IN_ROW; - } // j - - if (j_final_start < local.jend) - _cache[local.istart + (j_final_start+1)*_nx] |= MASK_NO_MORE_STARTS; - - if (_nan_separated && local.line_count > 1) { - // If _nan_separated, each line after the first has an extra nan to separate it from the - // previous line's points. If we were returning line offsets to the caller then this - // would need to occur in line() where the line_count is incremented. But as we are not - // returning line offsets it is faster just to add the extra here all at once. - local.total_point_count += local.line_count - 1; - } - - if (local.pass == 0) { - if (local.total_point_count == 0) { - local.points.clear(); - local.line_offsets.clear(); - local.outer_offsets.clear(); - break; // Do not need pass 1. - } - - // Create arrays for points, line_offsets and optionally outer_offsets. Arrays may be - // either C++ vectors or Python NumPy arrays. Want to group creation of the latter as - // threaded code needs to lock creation of these to limit access to a single thread. - if (_direct_points || _direct_line_offsets || _direct_outer_offsets) { - typename Derived::Lock lock(static_cast<Derived&>(*this)); - - // Strictly speaking adding the NumPy arrays to return_lists does not need to be - // within the lock. - if (_direct_points) { - return_lists[0][local.chunk] = - local.points.create_python(local.total_point_count, 2); - } - if (_direct_line_offsets) { - return_lists[1][local.chunk] = - local.line_offsets.create_python(local.line_count + 1); - } - if (_direct_outer_offsets) { - return_lists[2][local.chunk] = - local.outer_offsets.create_python(local.line_count - local.hole_count + 1); - } - } - - if (!_direct_points) - local.points.create_cpp(2*local.total_point_count); - - if (!_direct_line_offsets) - local.line_offsets.create_cpp(local.line_count + 1); - - if (!_direct_outer_offsets) { - if (_identify_holes) - local.outer_offsets.create_cpp( local.line_count - local.hole_count + 1); - else - local.outer_offsets.clear(); - } - - // Reset counts for pass 1. - local.total_point_count = 0; - local.line_count = 0; - local.hole_count = 0; - } - } // pass - - // Set final line and outer offsets. - if (local.line_count > 0) { - *local.line_offsets.current++ = local.total_point_count; - - if (_identify_holes) { - if (_outer_offsets_into_points) - *local.outer_offsets.current++ = local.total_point_count; - else - *local.outer_offsets.current++ = local.line_count; - } - } - - // Throw exception if the two passes returned different number of points, lines, etc. - check_consistent_counts(local); - - if (local.total_point_count == 0) { - if (_output_chunked) { - typename Derived::Lock lock(static_cast<Derived&>(*this)); - for (auto& list : return_lists) - list[local.chunk] = py::none(); - } - } - else if (_filled) - static_cast<Derived*>(this)->export_filled(local, return_lists); - else - static_cast<Derived*>(this)->export_lines(local, return_lists); -} - -template <typename Derived> -py::sequence BaseContourGenerator<Derived>::march_wrapper() -{ - index_t list_len = _n_chunks; - if ((_filled && (_fill_type == FillType::OuterCode|| _fill_type == FillType::OuterOffset)) || - (!_filled && (_line_type == LineType::Separate || _line_type == LineType::SeparateCode))) - list_len = 0; - - // Prepare lists to return to python. - std::vector<py::list> return_lists; - return_lists.reserve(_return_list_count); - for (decltype(_return_list_count) i = 0; i < _return_list_count; ++i) - return_lists.emplace_back(list_len); - - static_cast<Derived*>(this)->march(return_lists); - - // Return to python objects. - if (_return_list_count == 1) { - assert(!_filled); - if (_line_type == LineType::Separate) - return return_lists[0]; - else { - assert(_line_type == LineType::ChunkCombinedNan); - return py::make_tuple(return_lists[0]); - } - } - else if (_return_list_count == 2) - return py::make_tuple(return_lists[0], return_lists[1]); - else { - assert(_return_list_count == 3); - return py::make_tuple(return_lists[0], return_lists[1], return_lists[2]); - } -} - -template <typename Derived> -void BaseContourGenerator<Derived>::move_to_next_boundary_edge( - index_t& quad, index_t& forward, index_t& left) const -{ - // edge == 0 for E edge (facing N), forward = +_nx - // 2 for S edge (facing E), forward = +1 - // 4 for W edge (facing S), forward = -_nx - // 6 for N edge (facing W), forward = -1 - // 1 for SE edge (NW corner) from SW facing NE, forward = +_nx+1 - // 3 for SW edge (NE corner) from NW facing SE, forward = -_nx+1 - // 5 for NW edge (SE corner) from NE facing SW, forward = -_nx-1 - // 7 for NE edge (SW corner) from SE facing NW, forward = +_nx-1 - int edge = 0; - - // Need index of quad that is the same as the end point, i.e. quad to SW of end point, as it is - // this point which we need to find the next available boundary of, looking clockwise. - if (forward > 0) { - if (forward == _nx) { - assert(left == -1); - // W edge facing N, no change to quad or edge. - } - else if (left == _nx) { - assert(forward == 1); - quad -= _nx; // S edge facing E. - edge = 2; - } - else if (EXISTS_SW_CORNER(quad)) { - assert(forward == _nx-1 && left == -_nx-1); - quad -= 1; - edge = 7; - } - else { - assert(EXISTS_NW_CORNER(quad) && forward == _nx+1 && _nx-1); - // quad unchanged. - edge = 1; - } - } - else { // forward < 0 - if (forward == -_nx) { - assert(left == 1); - quad -= _nx+1; // W edge facing S. - edge = 4; - } - else if (left == -_nx) { - assert(forward == -1); - quad -= 1; // N edge facing W. - edge = 6; - } - else if (EXISTS_NE_CORNER(quad)) { - assert(forward == -_nx+1 && left == _nx+1); - quad -= _nx; - edge = 3; - } - else { - assert(EXISTS_SE_CORNER(quad) && forward == -_nx-1 && left == -_nx+1); - quad -= _nx+1; - edge = 5; - } - } - - // If _corner_mask not set, only need to consider odd edge in loop below. - if (!_corner_mask) - ++edge; - - while (true) { - // Look at possible edges that leave NE point of quad. - // If something is wrong here or in the setup of the boundary flags, can end up with an - // infinite loop! - switch (edge) { - case 0: - // Is there an edge to follow towards SW? - if (EXISTS_SE_CORNER(quad)) { // Equivalent to BOUNDARY_NE. - // quad unchanged. - forward = -_nx-1; - left = -_nx+1; - return; - } - break; - case 1: - // Is there an edge to follow towards W? - if (BOUNDARY_N(quad)) { - // quad unchanged. - forward = -1; - left = -_nx; - return; - } - break; - case 2: - // Is there an edge to follow towards NW? - if (EXISTS_SW_CORNER(quad+_nx)) { // Equivalent to BOUNDARY_NE. - quad += _nx; - forward = _nx-1; - left = -_nx-1; - return; - } - break; - case 3: - // Is there an edge to follow towards N? - if (BOUNDARY_E(quad+_nx)) { // Really a BOUNDARY_W check. - quad += _nx; - forward = _nx; - left = -1; - return; - } - break; - case 4: - // Is there an edge to follow towards NE? - if (EXISTS_NW_CORNER(quad+_nx+1)) { // Equivalent to BOUNDARY_SE. - quad += _nx+1; - forward = _nx+1; - left = _nx-1; - return; - } - break; - case 5: - // Is there an edge to follow towards E? - if (BOUNDARY_N(quad+1)) { // Really a BOUNDARY_S check - quad += _nx+1; - forward = 1; - left = _nx; - return; - } - break; - case 6: - // Is there an edge to follow towards SE? - if (EXISTS_NE_CORNER(quad+1)) { // Equivalent to BOUNDARY_SW. - quad += 1; - forward = -_nx+1; - left = _nx+1; - return; - } - break; - case 7: - // Is there an edge to follow towards S? - if (BOUNDARY_E(quad)) { - quad += 1; - forward = -_nx; - left = 1; - return; - } - break; - default: - assert(0 && "Invalid edge index"); - break; - } - - edge = _corner_mask ? (edge + 1) % 8 : (edge + 2) % 8; - } -} - -template <typename Derived> -py::list BaseContourGenerator<Derived>::multi_filled(const LevelArray levels) -{ - check_levels(levels, true); - pre_filled(); - - auto levels_proxy = levels.unchecked<1>(); - auto n = levels_proxy.size(); - - py::list ret(n-1); - _lower_level = levels_proxy[0]; - for (decltype(n) i = 0; i < n-1; i++) { - _upper_level = levels_proxy[i+1]; - ret[i] = march_wrapper(); - _lower_level = _upper_level; - } - - return ret; -} - -template <typename Derived> -py::list BaseContourGenerator<Derived>::multi_lines(const LevelArray levels) -{ - check_levels(levels, false); - pre_lines(); - - auto levels_proxy = levels.unchecked<1>(); - auto n = levels_proxy.size(); - - py::list ret(n); - for (decltype(n) i = 0; i < n; i++) { - _lower_level = _upper_level = levels_proxy[i]; - ret[i] = march_wrapper(); - } - - return ret; -} - -template <typename Derived> -void BaseContourGenerator<Derived>::pre_filled() -{ - _filled = true; - - _identify_holes = !(_fill_type == FillType::ChunkCombinedCode || - _fill_type == FillType::ChunkCombinedOffset); - _output_chunked = !(_fill_type == FillType::OuterCode || _fill_type == FillType::OuterOffset); - _direct_points = _output_chunked; - _direct_line_offsets = (_fill_type == FillType::ChunkCombinedOffset|| - _fill_type == FillType::ChunkCombinedOffsetOffset); - _direct_outer_offsets = (_fill_type == FillType::ChunkCombinedCodeOffset || - _fill_type == FillType::ChunkCombinedOffsetOffset); - _outer_offsets_into_points = (_fill_type == FillType::ChunkCombinedCodeOffset); - _nan_separated = false; - _return_list_count = (_fill_type == FillType::ChunkCombinedCodeOffset || - _fill_type == FillType::ChunkCombinedOffsetOffset) ? 3 : 2; -} - -template <typename Derived> -void BaseContourGenerator<Derived>::pre_lines() -{ - _filled = false; - - _identify_holes = false; - _output_chunked = !(_line_type == LineType::Separate || _line_type == LineType::SeparateCode); - _direct_points = _output_chunked; - _direct_line_offsets = (_line_type == LineType::ChunkCombinedOffset); - _direct_outer_offsets = false; - _outer_offsets_into_points = false; - _return_list_count = (_line_type == LineType::Separate || - _line_type == LineType::ChunkCombinedNan) ? 1 : 2; - _nan_separated = (_line_type == LineType::ChunkCombinedNan); - - if (_nan_separated) - Util::ensure_nan_loaded(); -} - -template <typename Derived> -void BaseContourGenerator<Derived>::set_look_flags(index_t hole_start_quad) -{ - assert(_identify_holes); - - // The only possible hole starts are START_E (from E to N), START_HOLE_N (on N boundary, E to W) - // and START_CORNER for SW corner (on boundary, SE to NW). - assert(hole_start_quad >= 0 && hole_start_quad < _n); - assert(EXISTS_N_EDGE(hole_start_quad) || EXISTS_SW_CORNER(hole_start_quad)); - assert(!LOOK_S(hole_start_quad) && "Look S already set"); - - _cache[hole_start_quad] |= MASK_LOOK_S; - - // Walk S until find place to mark corresponding look N. - auto quad = hole_start_quad; - - while (true) { - assert(quad >= 0 && quad < _n); - assert(EXISTS_N_EDGE(quad) || (quad == hole_start_quad && EXISTS_SW_CORNER(quad))); - - if (BOUNDARY_S(quad) || EXISTS_NE_CORNER(quad) || EXISTS_NW_CORNER(quad) || Z_SE != 1) { - assert(!LOOK_N(quad) && "Look N already set"); - _cache[quad] |= MASK_LOOK_N; - break; - } - - quad -= _nx; - } -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::supports_fill_type(FillType fill_type) -{ - switch (fill_type) { - case FillType::OuterCode: - case FillType::OuterOffset: - case FillType::ChunkCombinedCode: - case FillType::ChunkCombinedOffset: - case FillType::ChunkCombinedCodeOffset: - case FillType::ChunkCombinedOffsetOffset: - return true; - default: - return false; - } -} - -template <typename Derived> -bool BaseContourGenerator<Derived>::supports_line_type(LineType line_type) -{ - switch (line_type) { - case LineType::Separate: - case LineType::SeparateCode: - case LineType::ChunkCombinedCode: - case LineType::ChunkCombinedOffset: - case LineType::ChunkCombinedNan: - return true; - default: - return false; - } -} - -template <typename Derived> -void BaseContourGenerator<Derived>::write_cache() const -{ - std::cout << "---------- Cache ----------" << std::endl; - index_t ny = _n / _nx; - for (index_t j = ny-1; j >= 0; --j) { - std::cout << "j=" << j << " "; - for (index_t i = 0; i < _nx; ++i) { - index_t quad = i + j*_nx; - write_cache_quad(quad); - } - std::cout << std::endl; - } - std::cout << " "; - for (index_t i = 0; i < _nx; ++i) - std::cout << "i=" << i << " "; - std::cout << std::endl; - std::cout << "---------------------------" << std::endl; -} - -template <typename Derived> -void BaseContourGenerator<Derived>::write_cache_quad(index_t quad) const -{ - assert(quad >= 0 && quad < _n && "quad index out of bounds"); - std::cout << (NO_MORE_STARTS(quad) ? 'x' : - (NO_STARTS_IN_ROW(quad) ? 'i' : '.')); - std::cout << (EXISTS_QUAD(quad) ? "Q_" : - (EXISTS_NW_CORNER(quad) ? "NW" : - (EXISTS_NE_CORNER(quad) ? "NE" : - (EXISTS_SW_CORNER(quad) ? "SW" : - (EXISTS_SE_CORNER(quad) ? "SE" : ".."))))); - std::cout << (BOUNDARY_N(quad) && BOUNDARY_E(quad) ? 'b' : ( - BOUNDARY_N(quad) ? 'n' : (BOUNDARY_E(quad) ? 'e' : '.'))); - std::cout << Z_LEVEL(quad); - std::cout << ((_cache[quad] & MASK_MIDDLE) >> 2); - std::cout << (START_BOUNDARY_S(quad) ? 's' : '.'); - std::cout << (START_BOUNDARY_W(quad) ? 'w' : '.'); - if (!_filled) { - std::cout << (START_BOUNDARY_E(quad) ? 'e' : '.'); - std::cout << (START_BOUNDARY_N(quad) ? 'n' : '.'); - } - std::cout << (START_E(quad) ? 'E' : '.'); - std::cout << (START_N(quad) ? 'N' : '.'); - if (_filled) - std::cout << (START_HOLE_N(quad) ? 'h' : '.'); - std::cout << (START_CORNER(quad) ? 'c' : '.'); - if (_filled) - std::cout << (LOOK_N(quad) && LOOK_S(quad) ? 'B' : - (LOOK_N(quad) ? '^' : (LOOK_S(quad) ? 'v' : '.'))); - std::cout << ' '; -} - -template <typename Derived> -typename BaseContourGenerator<Derived>::ZLevel BaseContourGenerator<Derived>::z_to_zlevel( - double z_value) const -{ - return (_filled && z_value > _upper_level) ? 2 : (z_value > _lower_level ? 1 : 0); -} - -} // namespace contourpy - -#endif // CONTOURPY_BASE_IMPL_H diff --git a/contrib/python/contourpy/src/chunk_local.cpp b/contrib/python/contourpy/src/chunk_local.cpp deleted file mode 100644 index 1c8833ab0b5..00000000000 --- a/contrib/python/contourpy/src/chunk_local.cpp +++ /dev/null @@ -1,59 +0,0 @@ -#include "chunk_local.h" -#include <iostream> -#include <limits> - -namespace contourpy { - -ChunkLocal::ChunkLocal() -{ - look_up_quads.reserve(100); - clear(); -} - -void ChunkLocal::clear() -{ - chunk = -1; - istart = iend = jstart = jend = -1; - pass = -1; - - total_point_count = 0; - line_count = 0; - hole_count = 0; - - points.clear(); - line_offsets.clear(); - outer_offsets.clear(); - - look_up_quads.clear(); -} - -std::ostream &operator<<(std::ostream &os, const ChunkLocal& local) -{ - os << "ChunkLocal:" - << " chunk=" << local.chunk - << " istart=" << local.istart - << " iend=" << local.iend - << " jstart=" << local.jstart - << " jend=" << local.jend - << " total_point_count=" << local.total_point_count - << " line_count=" << local.line_count - << " hole_count=" << local.hole_count; - - if (local.line_offsets.start != nullptr) { - os << " line_offsets="; - for (count_t i = 0; i < local.line_count + 1; ++i) { - os << local.line_offsets.start[i] << " "; - } - } - - if (local.outer_offsets.start != nullptr) { - os << " outer_offsets="; - for (count_t i = 0; i < local.line_count - local.hole_count + 1; ++i) { - os << local.outer_offsets.start[i] << " "; - } - } - - return os; -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/chunk_local.h b/contrib/python/contourpy/src/chunk_local.h deleted file mode 100644 index 0298ea4814b..00000000000 --- a/contrib/python/contourpy/src/chunk_local.h +++ /dev/null @@ -1,40 +0,0 @@ -#ifndef CONTOURPY_CHUNK_LOCAL_H -#define CONTOURPY_CHUNK_LOCAL_H - -#include "output_array.h" -#include <iosfwd> - -namespace contourpy { - -struct ChunkLocal -{ - ChunkLocal(); - - void clear(); - - friend std::ostream &operator<<(std::ostream &os, const ChunkLocal& local); - - - - index_t chunk; // Index in range 0 to _n_chunks-1. - index_t istart, iend, jstart, jend; // Chunk limits, inclusive. - int pass; - - // Data for whole pass. - count_t total_point_count; // Includes nan separators if used. - count_t line_count; // Count of all lines - count_t hole_count; // Count of holes only. - - // Output arrays that are initialised at the end of pass 0 and written to during pass 1. - OutputArray<double> points; - OutputArray<offset_t> line_offsets; // Into array of points. - OutputArray<offset_t> outer_offsets; // Into array of points or line offsets depending on - // fill_type. - - // Data for current outer. - std::vector<index_t> look_up_quads; // To find holes of current outer. -}; - -} // namespace contourpy - -#endif // CONTOURPY_CHUNK_LOCAL_H diff --git a/contrib/python/contourpy/src/common.h b/contrib/python/contourpy/src/common.h deleted file mode 100644 index 68c522d1b6e..00000000000 --- a/contrib/python/contourpy/src/common.h +++ /dev/null @@ -1,33 +0,0 @@ -#ifndef CONTOURPY_COMMON_H -#define CONTOURPY_COMMON_H - -#include <pybind11/pybind11.h> -#include <pybind11/numpy.h> - -namespace contourpy { - -namespace py = pybind11; - -// quad/point index type, the same as for numpy array shape/indices (== npy_intp). Also used for -// chunks. -typedef py::ssize_t index_t; - -// Count of points, lines and holes. -typedef py::size_t count_t; - -// Offsets into point arrays. -typedef uint32_t offset_t; - -// Input numpy array classes. -typedef py::array_t<double, py::array::c_style | py::array::forcecast> CoordinateArray; -typedef py::array_t<bool, py::array::c_style | py::array::forcecast> MaskArray; -typedef py::array_t<double> LevelArray; // Doesn't have to be contiguous. - -// Output numpy array classes. -typedef py::array_t<double> PointArray; -typedef py::array_t<uint8_t> CodeArray; -typedef py::array_t<offset_t> OffsetArray; - -} // namespace contourpy - -#endif // CONTOURPY_COMMON_H diff --git a/contrib/python/contourpy/src/contour_generator.cpp b/contrib/python/contourpy/src/contour_generator.cpp deleted file mode 100644 index 5ede6fa46ee..00000000000 --- a/contrib/python/contourpy/src/contour_generator.cpp +++ /dev/null @@ -1,79 +0,0 @@ -#include "contour_generator.h" -#include "util.h" - -namespace contourpy { - -void ContourGenerator::check_levels(const LevelArray& levels, bool filled) const -{ - if (levels.ndim() != 1) { - throw std::domain_error( - "Levels array must be 1D not " + std::to_string(levels.ndim()) + "D"); - } - - if (filled) { - auto n = levels.size(); - if (n < 2) { - throw std::invalid_argument( - "Levels array must have at least 2 elements, not " + std::to_string(n)); - } - - auto levels_proxy = levels.unchecked<1>(); - - // Check levels are not NaN. - for (decltype(n) i = 0; i < n; i++) { - if (Util::is_nan(levels_proxy[i])) - throw std::invalid_argument("Levels must not contain any NaN"); - } - - // Check levels are increasing. - auto lower_level = levels_proxy[0]; - for (decltype(n) i = 0; i < n-1; i++) { - auto upper_level = levels_proxy[i+1]; - if (lower_level >= upper_level) - throw std::invalid_argument("Levels must be increasing"); - lower_level = upper_level; - } - } -} - -void ContourGenerator::check_levels(double lower_level, double upper_level) const -{ - if (Util::is_nan(lower_level) || Util::is_nan(upper_level)) - throw std::invalid_argument("lower_level and upper_level cannot be NaN"); - if (lower_level >= upper_level) - throw std::invalid_argument("upper_level must be larger than lower_level"); -} - -py::list ContourGenerator::multi_filled(const LevelArray levels) -{ - check_levels(levels, true); - - auto levels_proxy = levels.unchecked<1>(); - auto n = levels_proxy.size(); - - py::list ret(n-1); - auto lower_level = levels_proxy[0]; - for (decltype(n) i = 0; i < n-1; i++) { - auto upper_level = levels_proxy[i+1]; - ret[i] = filled(lower_level, upper_level); - lower_level = upper_level; - } - - return ret; -} - -py::list ContourGenerator::multi_lines(const LevelArray levels) -{ - check_levels(levels, false); - - auto levels_proxy = levels.unchecked<1>(); - auto n = levels_proxy.size(); - - py::list ret(n); - for (decltype(n) i = 0; i < n; i++) - ret[i] = lines(levels_proxy[i]); - - return ret; -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/contour_generator.h b/contrib/python/contourpy/src/contour_generator.h deleted file mode 100644 index c4375800239..00000000000 --- a/contrib/python/contourpy/src/contour_generator.h +++ /dev/null @@ -1,36 +0,0 @@ -#ifndef CONTOURPY_CONTOUR_GENERATOR_H -#define CONTOURPY_CONTOUR_GENERATOR_H - -#include "common.h" - -namespace contourpy { - -class ContourGenerator -{ -public: - // Non-copyable and non-moveable. - ContourGenerator(const ContourGenerator& other) = delete; - ContourGenerator(const ContourGenerator&& other) = delete; - ContourGenerator& operator=(const ContourGenerator& other) = delete; - ContourGenerator& operator=(const ContourGenerator&& other) = delete; - - virtual ~ContourGenerator() = default; - - virtual py::tuple filled(double lower_level, double upper_level) = 0; - - virtual py::sequence lines(double level) = 0; - - virtual py::list multi_filled(const LevelArray levels); - virtual py::list multi_lines(const LevelArray levels); - -protected: - ContourGenerator() = default; - - // Check levels are acceptable, throw exception if not. - void check_levels(const LevelArray& levels, bool filled) const; - void check_levels(double lower_level, double upper_level) const; -}; - -} // namespace contourpy - -#endif // CONTOURPY_CONTOUR_GENERATOR_H diff --git a/contrib/python/contourpy/src/converter.cpp b/contrib/python/contourpy/src/converter.cpp deleted file mode 100644 index 6a800e2ec86..00000000000 --- a/contrib/python/contourpy/src/converter.cpp +++ /dev/null @@ -1,155 +0,0 @@ -#include "converter.h" -#include "mpl_kind_code.h" -#include <limits> - -namespace contourpy { - -void Converter::check_max_offset(count_t max_offset) -{ - if (max_offset > std::numeric_limits<OffsetArray::value_type>::max()) - throw std::range_error("Max offset too large to fit in np.uint32. Use smaller chunks."); -} - -CodeArray Converter::convert_codes( - count_t point_count, count_t cut_count, const offset_t* cut_start, offset_t subtract) -{ - assert(point_count > 0 && cut_count > 0); - assert(cut_start != nullptr); - - index_t codes_shape = static_cast<index_t>(point_count); - CodeArray py_codes(codes_shape); - convert_codes(point_count, cut_count, cut_start, subtract, py_codes.mutable_data()); - return py_codes; -} - -void Converter::convert_codes( - count_t point_count, count_t cut_count, const offset_t* cut_start, offset_t subtract, - CodeArray::value_type* codes) -{ - assert(point_count > 0 && cut_count > 0); - assert(cut_start != nullptr); - assert(codes != nullptr); - - std::fill(codes + 1, codes + point_count - 1, LINETO); - for (decltype(cut_count) i = 0; i < cut_count-1; ++i) { - codes[cut_start[i] - subtract] = MOVETO; - codes[cut_start[i+1] - 1 - subtract] = CLOSEPOLY; - } -} - -CodeArray Converter::convert_codes_check_closed( - count_t point_count, count_t cut_count, const offset_t* cut_start, const double* points) -{ - assert(point_count > 0 && cut_count > 0); - assert(cut_start != nullptr); - assert(points != nullptr); - - index_t codes_shape = static_cast<index_t>(point_count); - CodeArray codes(codes_shape); - convert_codes_check_closed(point_count, cut_count, cut_start, points, codes.mutable_data()); - return codes; -} - -void Converter::convert_codes_check_closed( - count_t point_count, count_t cut_count, const offset_t* cut_start, const double* points, - CodeArray::value_type* codes) -{ - assert(point_count > 0 && cut_count > 0); - assert(cut_start != nullptr); - assert(points != nullptr); - assert(codes != nullptr); - - std::fill(codes + 1, codes + point_count, LINETO); - for (decltype(cut_count) i = 0; i < cut_count-1; ++i) { - auto start = cut_start[i]; - auto end = cut_start[i+1]; - codes[start] = MOVETO; - bool closed = points[2*start] == points[2*end-2] && - points[2*start+1] == points[2*end-1]; - if (closed) - codes[end-1] = CLOSEPOLY; - } -} - -CodeArray Converter::convert_codes_check_closed_single( - count_t point_count, const double* points) -{ - assert(point_count > 0); - assert(points != nullptr); - - index_t codes_shape = static_cast<index_t>(point_count); - CodeArray py_codes(codes_shape); - convert_codes_check_closed_single(point_count, points, py_codes.mutable_data()); - return py_codes; -} - -void Converter::convert_codes_check_closed_single( - count_t point_count, const double* points, CodeArray::value_type* codes) -{ - assert(point_count > 0); - assert(points != nullptr); - assert(codes != nullptr); - - codes[0] = MOVETO; - auto start = points; - auto end = points + 2*point_count; - bool closed = *start == *(end-2) && *(start+1) == *(end-1); - if (closed) { - std::fill(codes + 1, codes + point_count - 1, LINETO); - codes[point_count-1] = CLOSEPOLY; - } - else - std::fill(codes + 1, codes + point_count, LINETO); -} - -OffsetArray Converter::convert_offsets( - count_t offset_count, const offset_t* start, offset_t subtract) -{ - assert(offset_count > 0); - assert(start != nullptr); - - index_t offsets_shape = static_cast<index_t>(offset_count); - OffsetArray py_offsets(offsets_shape); - convert_offsets(offset_count, start, subtract, py_offsets.mutable_data()); - return py_offsets; -} - -void Converter::convert_offsets( - count_t offset_count, const offset_t* start, offset_t subtract, - OffsetArray::value_type* offsets) -{ - assert(offset_count > 0); - assert(start != nullptr); - assert(offsets != nullptr); - - check_max_offset(*(start + offset_count - 1) - subtract); - - if (subtract == 0) - std::copy(start, start + offset_count, offsets); - else { - for (decltype(offset_count) i = 0; i < offset_count; ++i) - *offsets++ = start[i] - subtract; - } -} - -PointArray Converter::convert_points(count_t point_count, const double* start) -{ - assert(point_count > 0); - assert(start != nullptr); - - index_t points_shape[2] = {static_cast<index_t>(point_count), 2}; - PointArray py_points(points_shape); - convert_points(point_count, start, py_points.mutable_data()); - return py_points; -} - -void Converter::convert_points(count_t point_count, const double* start, double* points) -{ - assert(point_count > 0); - assert(start != nullptr); - assert(points != nullptr); - - std::copy(start, start + 2*point_count, points); -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/converter.h b/contrib/python/contourpy/src/converter.h deleted file mode 100644 index 5bf9ab5119d..00000000000 --- a/contrib/python/contourpy/src/converter.h +++ /dev/null @@ -1,63 +0,0 @@ -#ifndef CONTOURPY_CONVERTER_H -#define CONTOURPY_CONVERTER_H - -#include "common.h" - -namespace contourpy { - -// Conversion of C++ point/code/offset array to return to Python as a NumPy array. -// There are two versions of each function, the first creates, populates and returns a NumPy array -// whereas the second populates one that has already been created. The former are used in serial -// code and the latter in threaded code where the creation and manipulation of NumPy arrays needs to -// be threadlocked whereas the population of those arrays does not. -class Converter -{ -public: - // Create and populate codes array, - static CodeArray convert_codes( - count_t point_count, count_t cut_count, const offset_t* cut_start, offset_t subtract); - - // Populate codes array that has already been created. - static void convert_codes( - count_t point_count, count_t cut_count, const offset_t* cut_start, offset_t subtract, - CodeArray::value_type* codes); - - // Create and populate codes array, - static CodeArray convert_codes_check_closed( - count_t point_count, count_t cut_count, const offset_t* cut_start, const double* points); - - // Populate codes array that has already been created. - static void convert_codes_check_closed( - count_t point_count, count_t cut_count, const offset_t* cut_start, const double* points, - CodeArray::value_type* codes); - - // Create and populate codes array (single line loop/strip). - static CodeArray convert_codes_check_closed_single( - count_t point_count, const double* points); - - // Populate codes array that has already been created (single line loop/strip). - static void convert_codes_check_closed_single( - count_t point_count, const double* points, CodeArray::value_type* codes); - - // Create and populate offsets array, - static OffsetArray convert_offsets( - count_t offset_count, const offset_t* start, offset_t subtract); - - // Populate offsets array that has already been created. - static void convert_offsets( - count_t offset_count, const offset_t* start, offset_t subtract, - OffsetArray::value_type* offsets); - - // Create and populate points array, - static PointArray convert_points(count_t point_count, const double* start); - - // Populate points array that has already been created. - static void convert_points(count_t point_count, const double* start, double* points); - -private: - static void check_max_offset(count_t max_offset); -}; - -} // namespace contourpy - -#endif // CONTOURPY_CONVERTER_H diff --git a/contrib/python/contourpy/src/fill_type.cpp b/contrib/python/contourpy/src/fill_type.cpp deleted file mode 100644 index 79129e29092..00000000000 --- a/contrib/python/contourpy/src/fill_type.cpp +++ /dev/null @@ -1,31 +0,0 @@ -#include "fill_type.h" -#include <iostream> - -namespace contourpy { - -std::ostream &operator<<(std::ostream &os, const FillType& fill_type) -{ - switch (fill_type) { - case FillType::OuterCode: - os << "OuterCode"; - break; - case FillType::OuterOffset: - os << "OuterOffset"; - break; - case FillType::ChunkCombinedCode: - os << "ChunkCombinedCode"; - break; - case FillType::ChunkCombinedOffset: - os << "ChunkCombinedOffset"; - break; - case FillType::ChunkCombinedCodeOffset: - os << "ChunkCombinedCodeOffset"; - break; - case FillType::ChunkCombinedOffsetOffset: - os << "ChunkCombinedOffsetOffset"; - break; - } - return os; -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/fill_type.h b/contrib/python/contourpy/src/fill_type.h deleted file mode 100644 index 20d41aeb3b2..00000000000 --- a/contrib/python/contourpy/src/fill_type.h +++ /dev/null @@ -1,24 +0,0 @@ -#ifndef CONTOURPY_FILL_TYPE_H -#define CONTOURPY_FILL_TYPE_H - -#include <iosfwd> -#include <string> - -namespace contourpy { - -// C++11 scoped enum, must be fully qualified to use. -enum class FillType -{ - OuterCode = 201, - OuterOffset = 202, - ChunkCombinedCode = 203, - ChunkCombinedOffset = 204, - ChunkCombinedCodeOffset = 205, - ChunkCombinedOffsetOffset = 206, -}; - -std::ostream &operator<<(std::ostream &os, const FillType& fill_type); - -} // namespace contourpy - -#endif // CONTOURPY_FILL_TYPE_H diff --git a/contrib/python/contourpy/src/line_type.cpp b/contrib/python/contourpy/src/line_type.cpp deleted file mode 100644 index fbd97ea78ce..00000000000 --- a/contrib/python/contourpy/src/line_type.cpp +++ /dev/null @@ -1,28 +0,0 @@ -#include "line_type.h" -#include <iostream> - -namespace contourpy { - -std::ostream &operator<<(std::ostream &os, const LineType& line_type) -{ - switch (line_type) { - case LineType::Separate: - os << "Separate"; - break; - case LineType::SeparateCode: - os << "SeparateCode"; - break; - case LineType::ChunkCombinedCode: - os << "ChunkCombinedCode"; - break; - case LineType::ChunkCombinedOffset: - os << "ChunkCombinedOffset"; - break; - case LineType::ChunkCombinedNan: - os << "ChunkCombinedNan"; - break; - } - return os; -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/line_type.h b/contrib/python/contourpy/src/line_type.h deleted file mode 100644 index c1e5dda702e..00000000000 --- a/contrib/python/contourpy/src/line_type.h +++ /dev/null @@ -1,23 +0,0 @@ -#ifndef CONTOURPY_LINE_TYPE_H -#define CONTOURPY_LINE_TYPE_H - -#include <iosfwd> -#include <string> - -namespace contourpy { - -// C++11 scoped enum, must be fully qualified to use. -enum class LineType -{ - Separate = 101, - SeparateCode = 102, - ChunkCombinedCode = 103, - ChunkCombinedOffset = 104, - ChunkCombinedNan = 105, -}; - -std::ostream &operator<<(std::ostream &os, const LineType& line_type); - -} // namespace contourpy - -#endif // CONTOURPY_LINE_TYPE_H diff --git a/contrib/python/contourpy/src/mpl2005.cpp b/contrib/python/contourpy/src/mpl2005.cpp deleted file mode 100644 index d94432a44f0..00000000000 --- a/contrib/python/contourpy/src/mpl2005.cpp +++ /dev/null @@ -1,74 +0,0 @@ -#include "mpl2005.h" - -namespace contourpy { - -Mpl2005ContourGenerator::Mpl2005ContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, index_t x_chunk_size, index_t y_chunk_size) - : _x(x), - _y(y), - _z(z), - _site(cntr_new()) -{ - if (_x.ndim() != 2 || _y.ndim() != 2 || _z.ndim() != 2) - throw std::invalid_argument("x, y and z must all be 2D arrays"); - - auto nx = _z.shape(1); - auto ny = _z.shape(0); - - if (_x.shape(1) != nx || _x.shape(0) != ny || - _y.shape(1) != nx || _y.shape(0) != ny) - throw std::invalid_argument("x, y and z arrays must have the same shape"); - - if (nx < 2 || ny < 2) - throw std::invalid_argument("x, y and z must all be at least 2x2 arrays"); - - if (mask.ndim() != 0) { // ndim == 0 if mask is not set, which is valid. - if (mask.ndim() != 2) - throw std::invalid_argument("mask array must be a 2D array"); - - if (mask.shape(1) != nx || mask.shape(0) != ny) - throw std::invalid_argument( - "If mask is set it must be a 2D array with the same shape as z"); - } - - if (x_chunk_size < 0 || y_chunk_size < 0) - throw std::invalid_argument("x_chunk_size and y_chunk_size cannot be negative"); - - const bool* mask_data = (mask.ndim() > 0 ? mask.data() : nullptr); - - cntr_init( - _site, nx, ny, _x.data(), _y.data(), _z.data(), mask_data, x_chunk_size, y_chunk_size); -} - -Mpl2005ContourGenerator::~Mpl2005ContourGenerator() -{ - cntr_del(_site); -} - -py::tuple Mpl2005ContourGenerator::filled(double lower_level, double upper_level) -{ - check_levels(lower_level, upper_level); - double levels[2] = {lower_level, upper_level}; - return cntr_trace(_site, levels, 2); -} - -py::tuple Mpl2005ContourGenerator::get_chunk_count() const -{ - long nx_chunks = (long)(ceil((_site->imax-1.0) / _site->i_chunk_size)); - long ny_chunks = (long)(ceil((_site->jmax-1.0) / _site->j_chunk_size)); - return py::make_tuple(ny_chunks, nx_chunks); -} - -py::tuple Mpl2005ContourGenerator::get_chunk_size() const -{ - return py::make_tuple(_site->j_chunk_size, _site->i_chunk_size); -} - -py::sequence Mpl2005ContourGenerator::lines(double level) -{ - double levels[2] = {level, 0.0}; - return cntr_trace(_site, levels, 1); -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/mpl2005.h b/contrib/python/contourpy/src/mpl2005.h deleted file mode 100644 index 5801f5d17b3..00000000000 --- a/contrib/python/contourpy/src/mpl2005.h +++ /dev/null @@ -1,32 +0,0 @@ -#ifndef CONTOURPY_MPL_2005_H -#define CONTOURPY_MPL_2005_H - -#include "contour_generator.h" -#include "mpl2005_original.h" - -namespace contourpy { - -class Mpl2005ContourGenerator : public ContourGenerator -{ -public: - Mpl2005ContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, index_t x_chunk_size, index_t y_chunk_size); - - virtual ~Mpl2005ContourGenerator(); - - py::tuple filled(double lower_level, double upper_level) override; - - py::tuple get_chunk_count() const; // Return (y_chunk_count, x_chunk_count) - py::tuple get_chunk_size() const; // Return (y_chunk_size, x_chunk_size) - - py::sequence lines(double level) override; - -private: - CoordinateArray _x, _y, _z; - Csite *_site; -}; - -} // namespace contourpy - -#endif // CONTOURPY_MPL_2005_H diff --git a/contrib/python/contourpy/src/mpl2005_original.cpp b/contrib/python/contourpy/src/mpl2005_original.cpp deleted file mode 100644 index f2f0e8567c3..00000000000 --- a/contrib/python/contourpy/src/mpl2005_original.cpp +++ /dev/null @@ -1,1526 +0,0 @@ -/* - cntr.c - General purpose contour tracer for quadrilateral meshes. - Handles single level contours, or region between a pair of levels. - - The routines that do all the work, as well as the explanatory - comments, came from gcntr.c, part of the GIST package. The - original mpl interface was also based on GIST. The present - interface uses parts of the original, but places them in - the entirely different framework of a Python type. It - was written by following the Python "Extending and Embedding" - tutorial. - */ - -#include "mpl2005_original.h" -#include "mpl_kind_code.h" - -/* Note that all arrays in these routines are Fortran-style, - in the sense that the "i" index varies fastest; the dimensions - of the corresponding C array are z[jmax][imax] in the notation - used here. We can identify i and j with the x and y dimensions, - respectively. - -*/ - -/* What is a contour? - * - * Given a quadrilateral mesh (x,y), and values of a z at the points - * of that mesh, we seek a set of polylines connecting points at a - * particular value of z. Each point on such a contour curve lies - * on an edge of the mesh, at a point linearly interpolated to the - * contour level z0 between the given values of z at the endpoints - * of the edge. - * - * Identifying these points is easy. Figuring out how to connect them - * into a curve -- or possibly a set of disjoint curves -- is difficult. - * Each disjoint curve may be either a closed circuit, or it may begin - * and end on a mesh boundary. - * - * One of the problems with a quadrilateral mesh is that when the z - * values at one pair of diagonally opposite points lie below z0, and - * the values at the other diagonal pair of the same zone lie above z0, - * all four edges of the zone are cut, and there is an ambiguity in - * how we should connect the points. I call this a saddle zone. - * The problem is that two disjoint curves cut through a saddle zone - * (I reject the alternative of connecting the opposite points to make - * a single self-intersecting curve, since those make ugly contour plots - * -- I've tried it). The solution is to determine the z value of the - * centre of the zone, which is the mean of the z values of the four - * corner points. If the centre z is higher than the contour level of - * interest and you are moving along the line with higher values on the - * left, turn right to leave the saddle zone. If the centre z is lower - * than the contour level turn left. Whether the centre z is higher - * than the 1 or 2 contour levels is stored in the saddle array so that - * it does not need to be recalculated in subsequent passes. - * - * Another complicating factor is that there may be logical holes in - * the mesh -- zones which do not exist. We want our contours to stop - * if they hit the edge of such a zone, just as if they'd hit the edge - * of the whole mesh. The input region array addresses this issue. - * - * Yet another complication: We may want a list of closed polygons which - * outline the region between two contour levels z0 and z1. These may - * include sections of the mesh boundary (including edges of logical - * holes defined by the region array), in addition to sections of the - * contour curves at one or both levels. This introduces a huge - * topological problem -- if one of the closed contours (possibly - * including an interior logical hole in the mesh, but not any part of - * the boundary of the whole mesh) encloses a region which is not - * between z0 and z1, that curve must be connected by a slit (or "branch - * cut") to the enclosing curve, so that the list of disjoint polygons - * we return is each simply connected. - * - * Okay, one final stunning difficulty: For the two level case, no - * individual polygon should have more than a few thousand sides, since - * huge filled polygons place an inordinate load on rendering software, - * which needs an amount of scratch space proportional to the number - * of sides it needs to fill. So in the two level case, we want to - * chunk the mesh into rectangular pieces of no more than, say, 30x30 - * zones, which keeps each returned polygon to less than a few thousand - * sides (the worst case is very very bad -- you can easily write down - * a function and two level values which produce a polygon that cuts - * every edge of the mesh twice). - */ - -/* - * Here is the numbering scheme for points, edges, and zones in - * the mesh -- note that each ij corresponds to one point, one zone, - * one i-edge (i=constant edge) and one j-edge (j=constant edge): - * - * (ij-1)-------(ij)-------(ij) - * | | - * | | - * | | - * (ij-1) (ij) (ij) - * | | - * | | - * | | - * (ij-iX-1)----(ij-iX)----(ij-iX) - * - * At each point, the function value is either 0, 1, or 2, depending - * on whether it is below z0, between z0 and z1, or above z1. - * Each zone either exists (1) or not (0). - * From these three bits of data, all of the curve connectivity follows. - * - * The tracing algorithm is naturally edge-based: Either you are at a - * point where a level cuts an edge, ready to step across a zone to - * another edge, or you are drawing the edge itself, if it happens to - * be a boundary with at least one section between z0 and z1. - * - * In either case, the edge is a directed edge -- either the zone - * you are advancing into is to its left or right, or you are actually - * drawing it. I always trace curves keeping the region between z0 and - * z1 to the left of the curve. If I'm tracing a boundary, I'm always - * moving CCW (counter clockwise) around the zone that exists. And if - * I'm about to cross a zone, I'll make the direction of the edge I'm - * sitting on be such that the zone I'm crossing is to its left. - * - * I start tracing each curve near its lower left corner (mesh oriented - * as above), which is the first point I encounter scanning through the - * mesh in order. When I figure the 012 z values and zonal existence, - * I also mark the potential starting points: Each edge may harbor a - * potential starting point corresponding to either direction, so there - * are four start possibilities at each ij point. Only the following - * possibilities need to be marked as potential starting edges: - * - * +-+-+-+ - * | | | | - * A-0-C-+ One or both levels cut E and have z=1 above them, and - * | EZ| | 0A is cut and either 0C is cut or CD is cut. - * +-B-D-+ Or, one or both levels cut E and E is a boundary edge. - * | | | | (and Z exists) - * +-+-+-+ - * - * +-+-+-+ - * | | | | - * +-A-0-C One or both levels cut E and have z=1 below them, and - * | |ZE | 0A is cut and either 0C is cut or CD is cut. - * +-+-B-D Or, one or both levels cut E and E is a boundary edge. - * | | | | (and Z exists) - * +-+-+-+ - * - * +-+-+-+ - * | | | | - * +-+-+-+ E is a boundary edge, Z exists, at some point on E - * | |Z| | lies between the levels. - * +-+E+-+ - * | | | | - * +-+-+-+ - * - * +-+-+-+ - * | | | | - * +-+E+-+ E is a boundary edge, Z exists, at some point on E - * | |Z| | lies between the levels. - * +-+-+-+ - * | | | | - * +-+-+-+ - * - * During the first tracing pass, the start mark is erased whenever - * any non-starting edge is encountered, reducing the number of points - * that need to be considered for the second pass. The first pass - * makes the basic connectivity decisions. It figures out how many - * disjoint curves there will be, and identifies slits for the two level - * case or open contours for the single level case, and removes all but - * the actual start markers. A second tracing pass can perform the - * actual final trace. - */ - -/* ------------------------------------------------------------------------ */ - -namespace contourpy { - -void print_Csite(Csite *Csite) -{ - Cdata *data = Csite->data; - int i, j, ij; - int nd = Csite->imax * (Csite->jmax + 1) + 1; - printf("zlevels: %8.2lg %8.2lg\n", Csite->zlevel[0], Csite->zlevel[1]); - printf("edge %ld, left %ld, n %ld, count %ld, edge0 %ld, left0 %ld\n", - Csite->edge, Csite->left, Csite->n, Csite->count, - Csite->edge0, Csite->left0); - printf(" level0 %d, edge00 %ld\n", Csite->level0, Csite->edge00); - printf("%04x\n", data[nd-1]); - for (j = Csite->jmax; j >= 0; j--) - { - for (i=0; i < Csite->imax; i++) - { - ij = i + j * Csite->imax; - printf("%04x ", data[ij]); - } - printf("\n"); - } - printf("\n"); -} - - -/* the Cdata array consists of the following bits: - * Z_VALUE (2 bits) 0, 1, or 2 function value at point - * ZONE_EX 1 zone exists, 0 zone doesn't exist - * I_BNDY this i-edge (i=constant edge) is a mesh boundary - * J_BNDY this j-edge (i=constant edge) is a mesh boundary - * I0_START this i-edge is a start point into zone to left - * I1_START this i-edge is a start point into zone to right - * J0_START this j-edge is a start point into zone below - * J1_START this j-edge is a start point into zone above - * START_ROW next start point is in current row (accelerates 2nd pass) - * SLIT_UP marks this i-edge as the beginning of a slit upstroke - * SLIT_DN marks this i-edge as the beginning of a slit downstroke - * OPEN_END marks an i-edge start point whose other endpoint is - * on a boundary for the single level case - * ALL_DONE marks final start point - * SLIT_DN_VISITED this slit downstroke hasn't/has been visited in pass 2 - */ -#define Z_VALUE 0x0003 -#define ZONE_EX 0x0004 -#define I_BNDY 0x0008 -#define J_BNDY 0x0010 -#define I0_START 0x0020 -#define I1_START 0x0040 -#define J0_START 0x0080 -#define J1_START 0x0100 -#define START_ROW 0x0200 -#define SLIT_UP 0x0400 -#define SLIT_DN 0x0800 -#define OPEN_END 0x1000 -#define ALL_DONE 0x2000 -#define SLIT_DN_VISITED 0x4000 - -/* some helpful macros to find points relative to a given directed - * edge -- points are designated 0, 1, 2, 3 CCW around zone with 0 and - * 1 the endpoints of the current edge */ -#define FORWARD(left,ix) ((left)>0?((left)>1?1:-(ix)):((left)<-1?-1:(ix))) -#define POINT0(edge,fwd) ((edge)-((fwd)>0?fwd:0)) -#define POINT1(edge,fwd) ((edge)+((fwd)<0?fwd:0)) -#define IS_JEDGE(edge,left) ((left)>0?((left)>1?1:0):((left)<-1?1:0)) -#define ANY_START (I0_START|I1_START|J0_START|J1_START) -#define START_MARK(left) \ - ((left)>0?((left)>1?J1_START:I1_START):((left)<-1?J0_START:I0_START)) - -enum {kind_zone, kind_edge1, kind_edge2, - kind_slit_up, kind_slit_down, kind_start_slit=16}; - -/* Saddle zone array consists of the following bits: - * SADDLE_SET whether zone's saddle data has been set. - * SADDLE_GT0 whether z of centre of zone is higher than site->level[0]. - * SADDLE_GT1 whether z of centre of zone is higher than site->level[1]. - */ -#define SADDLE_SET 0x01 -#define SADDLE_GT0 0x02 -#define SADDLE_GT1 0x04 - -/* ------------------------------------------------------------------------ */ - -/* these actually mark points */ -static int zone_crosser (Csite * site, int level, int pass2); -static int edge_walker (Csite * site, int pass2); -static int slit_cutter (Csite * site, int up, int pass2); - -/* this calls the first three to trace the next disjoint curve - * -- return value is number of points on this curve, or - * 0 if there are no more curves this pass - * -(number of points) on first pass if: - * this is two level case, and the curve closed on a hole - * this is single level case, curve is open, and will start from - * a different point on the second pass - * -- in both cases, this curve will be combined with another - * on the second pass */ -static long curve_tracer (Csite * site, int pass2); - -/* this initializes the data array for curve_tracer */ -static void data_init (Csite * site); - -/* ------------------------------------------------------------------------ */ - -/* zone_crosser assumes you are sitting at a cut edge about to cross - * the current zone. It always marks the initial point, crosses at - * least one zone, and marks the final point. On non-boundary i-edges, - * it is responsible for removing start markers on the first pass. */ -static int -zone_crosser (Csite * site, int level, int pass2) -{ - Cdata * data = site->data; - long edge = site->edge; - long left = site->left; - long n = site->n; - long fwd = FORWARD (left, site->imax); - long p0, p1; - int jedge = IS_JEDGE (edge, left); - long edge0 = site->edge0; - long left0 = site->left0; - int level0 = site->level0 == level; - int two_levels = site->zlevel[1] > site->zlevel[0]; - Saddle* saddle = site->saddle; - - const double *x = pass2 ? site->x : 0; - const double *y = pass2 ? site->y : 0; - const double *z = site->z; - double zlevel = site->zlevel[level]; - double *xcp = pass2 ? site->xcp : 0; - double *ycp = pass2 ? site->ycp : 0; - short *kcp = pass2 ? site->kcp : 0; - - int z0, z1, z2, z3; - int done = 0; - int n_kind; - - if (level) - level = 2; - - for (;;) - { - n_kind = 0; - /* set edge endpoints */ - p0 = POINT0 (edge, fwd); - p1 = POINT1 (edge, fwd); - - /* always mark cut on current edge */ - if (pass2) - { - /* second pass actually computes and stores the point */ - double zcp = (zlevel - z[p0]) / (z[p1] - z[p0]); - xcp[n] = zcp * (x[p1] - x[p0]) + x[p0]; - ycp[n] = zcp * (y[p1] - y[p0]) + y[p0]; - kcp[n] = kind_zone; - n_kind = n; - } - if (!done && !jedge) - { - if (n) - { - /* if this is not the first point on the curve, and we're - * not done, and this is an i-edge, check several things */ - if (!two_levels && !pass2 && (data[edge] & OPEN_END)) - { - /* reached an OPEN_END mark, skip the n++ */ - done = 4; /* same return value 4 used below */ - break; - } - - /* check for curve closure -- if not, erase any start mark */ - if (edge == edge0 && left == left0) - { - /* may signal closure on a downstroke */ - if (level0) - done = (!pass2 && two_levels && left < 0) ? 5 : 3; - } - else if (!pass2) - { - Cdata start = - data[edge] & (fwd > 0 ? I0_START : I1_START); - if (start) - { - data[edge] &= ~start; - site->count--; - } - if (!two_levels) - { - start = data[edge] & (fwd > 0 ? I1_START : I0_START); - if (start) - { - data[edge] &= ~start; - site->count--; - } - } - } - } - } - n++; - if (done) - break; - - /* cross current zone to another cut edge */ - z0 = (data[p0] & Z_VALUE) != level; /* 1 if fill toward p0 */ - z1 = !z0; /* know level cuts edge */ - z2 = (data[p1 + left] & Z_VALUE) != level; - z3 = (data[p0 + left] & Z_VALUE) != level; - if (z0 == z2) - { - if (z1 == z3) - { - /* this is a saddle zone, determine whether to turn left or - * right depending on height of centre of zone relative to - * contour level. Set saddle[zone] if not already decided. */ - int turnRight; - long zone = edge + (left > 0 ? left : 0); - if (!(saddle[zone] & SADDLE_SET)) - { - double zcentre; - saddle[zone] = SADDLE_SET; - zcentre = (z[p0] + z[p0+left] + z[p1] + z[p1+left])/4.0; - if (zcentre > site->zlevel[0]) - saddle[zone] |= - (two_levels && zcentre > site->zlevel[1]) - ? SADDLE_GT0 | SADDLE_GT1 : SADDLE_GT0; - } - - turnRight = level == 2 ? (saddle[zone] & SADDLE_GT1) - : (saddle[zone] & SADDLE_GT0); - if (z1 ^ (level == 2)) - turnRight = !turnRight; - if (!turnRight) - goto bkwd; - } - /* bend forward (right along curve) */ - jedge = !jedge; - edge = p1 + (left > 0 ? left : 0); - { - long tmp = fwd; - fwd = -left; - left = tmp; - } - } - else if (z1 == z3) - { - bkwd: - /* bend backward (left along curve) */ - jedge = !jedge; - edge = p0 + (left > 0 ? left : 0); - { - long tmp = fwd; - fwd = left; - left = -tmp; - } - } - else - { - /* straight across to opposite edge */ - edge += left; - } - /* after crossing zone, edge/left/fwd is oriented CCW relative to - * the next zone, assuming we will step there */ - - /* now that we've taken a step, check for the downstroke - * of a slit on the second pass (upstroke checked above) - * -- taking step first avoids a race condition */ - if (pass2 && two_levels && !jedge) - { - if (left > 0) - { - if (data[edge] & SLIT_UP) - done = 6; - } - else - { - if (data[edge] & SLIT_DN) - done = 5; - } - } - - if (!done) - { - /* finally, check if we are on a boundary */ - if (data[edge] & (jedge ? J_BNDY : I_BNDY)) - { - done = two_levels ? 2 : 4; - /* flip back into the zone that exists */ - left = -left; - fwd = -fwd; - if (!pass2 && (edge != edge0 || left != left0)) - { - Cdata start = data[edge] & START_MARK (left); - if (start) - { - data[edge] &= ~start; - site->count--; - } - } - } - } - } - - site->edge = edge; - site->n = n; - site->left = left; - if (done <= 4) - { - return done; - } - if (pass2 && n_kind) - { - kcp[n_kind] += kind_start_slit; - } - return slit_cutter (site, done - 5, pass2); -} - -/* edge_walker assumes that the current edge is being drawn CCW - * around the current zone. Since only boundary edges are drawn - * and we always walk around with the filled region to the left, - * no edge is ever drawn CW. We attempt to advance to the next - * edge on this boundary, but if current second endpoint is not - * between the two contour levels, we exit back to zone_crosser. - * Note that we may wind up marking no points. - * -- edge_walker is never called for single level case */ -static int -edge_walker (Csite * site, int pass2) -{ - Cdata * data = site->data; - long edge = site->edge; - long left = site->left; - long n = site->n; - long fwd = FORWARD (left, site->imax); - long p0 = POINT0 (edge, fwd); - long p1 = POINT1 (edge, fwd); - int jedge = IS_JEDGE (edge, left); - long edge0 = site->edge0; - long left0 = site->left0; - int level0 = site->level0 == 2; - int marked; - int n_kind = 0; - - const double *x = pass2 ? site->x : 0; - const double *y = pass2 ? site->y : 0; - double *xcp = pass2 ? site->xcp : 0; - double *ycp = pass2 ? site->ycp : 0; - short *kcp = pass2 ? site->kcp : 0; - - int z0, z1, heads_up = 0; - - for (;;) - { - /* mark endpoint 0 only if value is 1 there, and this is a - * two level task */ - z0 = data[p0] & Z_VALUE; - z1 = data[p1] & Z_VALUE; - marked = 0; - n_kind = 0; - if (z0 == 1) - { - /* mark current boundary point */ - if (pass2) - { - xcp[n] = x[p0]; - ycp[n] = y[p0]; - kcp[n] = kind_edge1; - n_kind = n; - } - marked = 1; - } - else if (!n) - { - /* if this is the first point is not between the levels - * must do the job of the zone_crosser and mark the first cut here, - * so that it will be marked again by zone_crosser as it closes */ - if (pass2) - { - double zcp = site->zlevel[(z0 != 0)]; - zcp = (zcp - site->z[p0]) / (site->z[p1] - site->z[p0]); - xcp[n] = zcp * (x[p1] - x[p0]) + x[p0]; - ycp[n] = zcp * (y[p1] - y[p0]) + y[p0]; - kcp[n] = kind_edge2; - n_kind = n; - } - marked = 1; - } - if (n) - { - /* check for closure */ - if (level0 && edge == edge0 && left == left0) - { - site->edge = edge; - site->left = left; - site->n = n + marked; - /* if the curve is closing on a hole, need to make a downslit */ - if (fwd < 0 && !(data[edge] & (jedge ? J_BNDY : I_BNDY))) - { - if (n_kind) kcp[n_kind] += kind_start_slit; - return slit_cutter (site, 0, pass2); - } - if (fwd < 0 && level0 && left < 0) - { - /* remove J0_START from this boundary edge as boundary is - * included by the upwards slit from contour line below. */ - data[edge] &= ~J0_START; - if (n_kind) kcp[n_kind] += kind_start_slit; - return slit_cutter (site, 0, pass2); - } - return 3; - } - else if (pass2) - { - if (heads_up || (fwd < 0 && (data[edge] & SLIT_DN))) - { - if (!heads_up && !(data[edge] & SLIT_DN_VISITED)) - data[edge] |= SLIT_DN_VISITED; - else - { - site->edge = edge; - site->left = left; - site->n = n + marked; - if (n_kind) kcp[n_kind] += kind_start_slit; - return slit_cutter (site, heads_up, pass2); - } - } - } - else - { - /* if this is not first point, clear start mark for this edge */ - Cdata start = data[edge] & START_MARK (left); - if (start) - { - data[edge] &= ~start; - site->count--; - } - } - } - if (marked) - n++; - - /* if next endpoint not between levels, need to exit to zone_crosser */ - if (z1 != 1) - { - site->edge = edge; - site->left = left; - site->n = n; - return (z1 != 0); /* return level closest to p1 */ - } - - /* step to p1 and find next edge - * -- turn left if possible, else straight, else right - * -- check for upward slit beginning at same time */ - edge = p1 + (left > 0 ? left : 0); - if (pass2 && jedge && fwd > 0 && (data[edge] & SLIT_UP)) - { - jedge = !jedge; - heads_up = 1; - } - else if (data[edge] & (jedge ? I_BNDY : J_BNDY)) - { - long tmp = fwd; - fwd = left; - left = -tmp; - jedge = !jedge; - } - else - { - edge = p1 + (fwd > 0 ? fwd : 0); - if (pass2 && !jedge && fwd > 0 && (data[edge] & SLIT_UP)) - { - heads_up = 1; - } - else if (!(data[edge] & (jedge ? J_BNDY : I_BNDY))) - { - edge = p1 - (left < 0 ? left : 0); - jedge = !jedge; - { - long tmp = fwd; - fwd = -left; - left = tmp; - } - } - } - p0 = p1; - p1 = POINT1 (edge, fwd); - } -} - -/* -- slit_cutter is never called for single level case */ -static int -slit_cutter (Csite * site, int up, int pass2) -{ - Cdata * data = site->data; - long imax = site->imax; - long n = site->n; - - const double *x = pass2 ? site->x : 0; - const double *y = pass2 ? site->y : 0; - double *xcp = pass2 ? site->xcp : 0; - double *ycp = pass2 ? site->ycp : 0; - short *kcp = pass2 ? site->kcp : 0; - - if (up && pass2) - { - /* upward stroke of slit proceeds up left side of slit until - * it hits a boundary or a point not between the contour levels - * -- this never happens on the first pass */ - long p1 = site->edge; - int z1; - - for (;;) - { - z1 = data[p1] & Z_VALUE; - if (z1 != 1) - { - site->edge = p1; - site->left = -1; - site->n = n; - return (z1 != 0); - } - else if (data[p1] & J_BNDY) - { - /* this is very unusual case of closing on a mesh hole */ - site->edge = p1; - site->left = -imax; - site->n = n; - return 2; - } - xcp[n] = x[p1]; - ycp[n] = y[p1]; - kcp[n] = kind_slit_up; - n++; - p1 += imax; - } - - } - else - { - /* downward stroke proceeds down right side of slit until it - * hits a boundary or point not between the contour levels */ - long p0 = site->edge; - int z0; - /* at beginning of first pass, mark first i-edge with SLIT_DN */ - data[p0] |= SLIT_DN; - p0 -= imax; - for (;;) - { - z0 = data[p0] & Z_VALUE; - if (!pass2) - { - if (z0 != 1 || (data[p0] & I_BNDY) || (data[p0 + 1] & J_BNDY)) - { - /* at end of first pass, mark final i-edge with SLIT_UP */ - data[p0 + imax] |= SLIT_UP; - /* one extra count for splicing at outer curve */ - site->n = n + 1; - return 4; /* return same special value as for OPEN_END */ - } - } - else - { - if (z0 != 1) - { - site->edge = p0 + imax; - site->left = 1; - site->n = n; - return (z0 != 0); - } - else if (data[p0 + 1] & J_BNDY) - { - site->edge = p0 + 1; - site->left = imax; - site->n = n; - return 2; - } - else if (data[p0] & I_BNDY) - { - site->edge = p0; - site->left = 1; - site->n = n; - return 2; - } - } - if (pass2) - { - xcp[n] = x[p0]; - ycp[n] = y[p0]; - kcp[n] = kind_slit_down; - n++; - } - else - { - /* on first pass need to count for upstroke as well */ - n += 2; - } - p0 -= imax; - } - } -} - -/* ------------------------------------------------------------------------ */ - -/* curve_tracer finds the next starting point, then traces the curve, - * returning the number of points on this curve - * -- in a two level trace, the return value is negative on the - * first pass if the curve closed on a hole - * -- in a single level trace, the return value is negative on the - * first pass if the curve is an incomplete open curve - * -- a return value of 0 indicates no more curves */ -static long -curve_tracer (Csite * site, int pass2) -{ - Cdata * data = site->data; - long imax = site->imax; - long edge0 = site->edge0; - long left0 = site->left0; - long edge00 = site->edge00; - int two_levels = site->zlevel[1] > site->zlevel[0]; - int level, level0, mark_row; - long n; - - /* it is possible for a single i-edge to serve as two actual start - * points, one to the right and one to the left - * -- for the two level case, this happens on the first pass for - * a doubly cut edge, or on a chunking boundary - * -- for single level case, this is impossible, but a similar - * situation involving open curves is handled below - * a second two start possibility is when the edge0 zone does not - * exist and both the i-edge and j-edge boundaries are cut - * yet another possibility is three start points at a junction - * of chunk cuts - * -- sigh, several other rare possibilities, - * allow for general case, just go in order i1, i0, j1, j0 */ - int two_starts; - /* printf("curve_tracer pass %d\n", pass2); */ - /* print_Csite(site); */ - if (left0 == 1) - two_starts = data[edge0] & (I0_START | J1_START | J0_START); - else if (left0 == -1) - two_starts = data[edge0] & (J1_START | J0_START); - else if (left0 == imax) - two_starts = data[edge0] & J0_START; - else - two_starts = 0; - - if (pass2 || edge0 == 0) - { - /* zip up to row marked on first pass (or by data_init if edge0==0) - * -- but not for double start case */ - if (!two_starts) - { - /* final start point marked by ALL_DONE marker */ - int first = (edge0 == 0 && !pass2); - long e0 = edge0; - if (data[edge0] & ALL_DONE) - return 0; - while (!(data[edge0] & START_ROW)) - edge0 += imax; - if (e0 == edge0) - edge0++; /* two starts handled specially */ - if (first) - /* if this is the very first start point, we want to remove - * the START_ROW marker placed by data_init */ - data[edge0 - edge0 % imax] &= ~START_ROW; - } - - } - else - { - /* first pass ends when all potential start points visited */ - if (site->count <= 0) - { - /* place ALL_DONE marker for second pass */ - data[edge00] |= ALL_DONE; - /* reset initial site for second pass */ - site->edge0 = site->edge00 = site->left0 = 0; - return 0; - } - if (!two_starts) - edge0++; - } - - if (two_starts) - { - /* trace second curve with this start immediately */ - if (left0 == 1 && (data[edge0] & I0_START)) - { - left0 = -1; - level = (data[edge0] & I_BNDY) ? 2 : 0; - } - else if ((left0 == 1 || left0 == -1) && (data[edge0] & J1_START)) - { - left0 = imax; - level = 2; - } - else - { - left0 = -imax; - level = 2; - } - - } - else - { - /* usual case is to scan for next start marker - * -- on second pass, this is at most one row of mesh, but first - * pass hits nearly every point of the mesh, since it can't - * know in advance which potential start marks removed */ - while (!(data[edge0] & ANY_START)) - edge0++; - - if (data[edge0] & I1_START) - left0 = 1; - else if (data[edge0] & I0_START) - left0 = -1; - else if (data[edge0] & J1_START) - left0 = imax; - else /*data[edge0]&J0_START */ - left0 = -imax; - - if (data[edge0] & (I1_START | I0_START)) - level = (data[edge0] & I_BNDY) ? 2 : 0; - else - level = 2; - } - - /* this start marker will not be unmarked, but it has been visited */ - if (!pass2) - site->count--; - - /* if this curve starts on a non-boundary i-edge, we need to - * determine the level */ - if (!level && two_levels) - level = left0 > 0 ? - ((data[edge0 - imax] & Z_VALUE) != - 0) : ((data[edge0] & Z_VALUE) != 0); - - /* initialize site for this curve */ - site->edge = site->edge0 = edge0; - site->left = site->left0 = left0; - site->level0 = level0 = level; /* for open curve detection only */ - - /* single level case just uses zone_crosser */ - if (!two_levels) - level = 0; - - /* to generate the curve, alternate between zone_crosser and - * edge_walker until closure or first call to edge_walker in - * single level case */ - site->n = 0; - for (;;) - { - if (level < 2) - level = zone_crosser (site, level, pass2); - else if (level < 3) - level = edge_walker (site, pass2); - else - break; - } - n = site->n; - - /* single level case may have ended at a boundary rather than closing - * -- need to recognize this case here in order to place the - * OPEN_END mark for zone_crosser, remove this start marker, - * and be sure not to make a START_ROW mark for this case - * two level case may close with slit_cutter, in which case start - * must also be removed and no START_ROW mark made - * -- change sign of return n to inform caller */ - if (!pass2 && level > 3 && (two_levels || level0 == 0)) - { - if (!two_levels) - data[edge0] |= OPEN_END; - data[edge0] &= ~(left0 > 0 ? I1_START : I0_START); - mark_row = 0; /* do not mark START_ROW */ - n = -n; - } - else - { - if (two_levels) - mark_row = !two_starts; - else - mark_row = 1; - } - - /* on first pass, must apply START_ROW mark in column above previous - * start marker - * -- but skip if we just did second of two start case */ - if (!pass2 && mark_row) - { - data[edge0 - (edge0 - edge00) % imax] |= START_ROW; - site->edge00 = edge0; - } - - return n; -} - -/* ------------------------------------------------------------------------ */ - -static void -data_init (Csite * site) -{ - Cdata * data = site->data; - long imax = site->imax; - long jmax = site->jmax; - long ijmax = imax * jmax; - const double *z = site->z; - double zlev0 = site->zlevel[0]; - double zlev1 = site->zlevel[1]; - int two_levels = zlev1 > zlev0; - char *reg = site->reg; - long count = 0; - int started = 0; - int ibndy, jbndy, i_was_chunk; - - long ichunk, jchunk, i, j, ij; - long i_chunk_size = site->i_chunk_size; - long j_chunk_size = site->j_chunk_size; - - if (!two_levels) - { - /* Chunking not used for lines as start points are not correct. */ - i_chunk_size = imax - 1; - j_chunk_size = jmax - 1; - } - - /* do everything in a single pass through the data array to - * minimize cache faulting (z, reg, and data are potentially - * very large arrays) - * access to the z and reg arrays is strictly sequential, - * but we need two rows (+-imax) of the data array at a time */ - if (z[0] > zlev0) - data[0] = (two_levels && z[0] > zlev1) ? 2 : 1; - else - data[0] = 0; - jchunk = 0; - for (j = ij = 0; j < jmax; j++) - { - ichunk = i_was_chunk = 0; - for (i = 0; i < imax; i++, ij++) - { - /* transfer zonal existence from reg to data array - * -- get these for next row so we can figure existence of - * points and j-edges for this row */ - data[ij + imax + 1] = 0; - if (reg) - { - if (reg[ij + imax + 1] != 0) - data[ij + imax + 1] = ZONE_EX; - } - else - { - if (i < imax - 1 && j < jmax - 1) - data[ij + imax + 1] = ZONE_EX; - } - - /* translate z values to 0, 1, 2 flags */ - if (ij < imax) - data[ij + 1] = 0; - if (ij < ijmax - 1 && z[ij + 1] > zlev0) - data[ij + 1] |= (two_levels && z[ij + 1] > zlev1) ? 2 : 1; - - /* apply edge boundary marks */ - ibndy = i == ichunk - || (data[ij] & ZONE_EX) != (data[ij + 1] & ZONE_EX); - jbndy = j == jchunk - || (data[ij] & ZONE_EX) != (data[ij + imax] & ZONE_EX); - if (ibndy) - data[ij] |= I_BNDY; - if (jbndy) - data[ij] |= J_BNDY; - - /* apply i-edge start marks - * -- i-edges are only marked when actually cut - * -- no mark is necessary if one of the j-edges which share - * the lower endpoint is also cut - * -- no I0 mark necessary unless filled region below some cut, - * no I1 mark necessary unless filled region above some cut */ - if (j) - { - int v0 = (data[ij] & Z_VALUE); - int vb = (data[ij - imax] & Z_VALUE); - if (v0 != vb) - { /* i-edge is cut */ - if (ibndy) - { - if (data[ij] & ZONE_EX) - { - data[ij] |= I0_START; - count++; - } - if (data[ij + 1] & ZONE_EX) - { - data[ij] |= I1_START; - count++; - } - } - else - { - int va = (data[ij - 1] & Z_VALUE); - int vc = (data[ij + 1] & Z_VALUE); - int vd = (data[ij - imax + 1] & Z_VALUE); - if (v0 != 1 && va != v0 - && (vc != v0 || vd != v0) && (data[ij] & ZONE_EX)) - { - data[ij] |= I0_START; - count++; - } - if (vb != 1 && va == vb - && (vc == vb || vd == vb) - && (data[ij + 1] & ZONE_EX)) - { - data[ij] |= I1_START; - count++; - } - } - } - } - - /* apply j-edge start marks - * -- j-edges are only marked when they are boundaries - * -- all cut boundary edges marked - * -- for two level case, a few uncut edges must be marked - */ - if (i && jbndy) - { - int v0 = (data[ij] & Z_VALUE); - int vb = (data[ij - 1] & Z_VALUE); - if (v0 != vb) - { - if (data[ij] & ZONE_EX) - { - data[ij] |= J0_START; - count++; - } - if (data[ij + imax] & ZONE_EX) - { - data[ij] |= J1_START; - count++; - } - } - else if (two_levels && v0 == 1) - { - if (data[ij + imax] & ZONE_EX) - { - if (i_was_chunk || !(data[ij + imax - 1] & ZONE_EX)) - { - /* lower left is a drawn part of boundary */ - data[ij] |= J1_START; - count++; - } - } - else if (data[ij] & ZONE_EX) - { - if (data[ij + imax - 1] & ZONE_EX) - { - /* weird case of open hole at lower left */ - data[ij] |= J0_START; - count++; - } - } - } - } - - i_was_chunk = (i == ichunk); - if (i_was_chunk) - ichunk += i_chunk_size; - } - - if (j == jchunk) - jchunk += j_chunk_size; - - /* place first START_ROW marker */ - if (count && !started) - { - data[ij - imax] |= START_ROW; - started = 1; - } - } - - /* place immediate stop mark if nothing found */ - if (!count) - data[0] |= ALL_DONE; - else - for (i = 0; i < ijmax; ++i) site->saddle[i] = 0; - - /* initialize site */ - site->edge0 = site->edge00 = site->edge = 0; - site->left0 = site->left = 0; - site->n = 0; - site->count = count; -} - -/* ------------------------------------------------------------------------ - Original (slightly modified) core contour generation routines are above; - below are new routines for interfacing to mpl. - - ------------------------------------------------------------------------ */ - -/* Note: index order gets switched in the Python interface; - python Z[i,j] -> C z[j,i] - so if the array has shape Mi, Nj in python, - we have iMax = Nj, jMax = Mi in gcntr.c. - On the Python side: Ny, Nx = shape(z), - so in C, the x-dimension is the first index, the y-dimension - the second. -*/ - -/* reg should have the same dimensions as data, which - has an extra iMax + 1 points relative to Z. - It differs from mask in being the opposite (True - where a region exists, versus the mask, which is True - where a data point is bad), and in that it marks - zones, not points. All four zones sharing a bad - point must be marked as not existing. -*/ -static void -mask_zones (long iMax, long jMax, const bool *mask, char *reg) -{ - long i, j, ij; - long nreg = iMax * jMax + iMax + 1; - - for (ij = iMax+1; ij < iMax*jMax; ij++) - { - reg[ij] = 1; - } - - ij = 0; - for (j = 0; j < jMax; j++) - { - for (i = 0; i < iMax; i++, ij++) - { - if (i == 0 || j == 0) reg[ij] = 0; - if (mask[ij]) - { - reg[ij] = 0; - reg[ij + 1] = 0; - reg[ij + iMax] = 0; - reg[ij + iMax + 1] = 0; - } - } - } - for (; ij < nreg; ij++) - { - reg[ij] = 0; - } -} - -Csite * -cntr_new() -{ - Csite *site = new Csite; - if (site == nullptr) return nullptr; - site->data = nullptr; - site->reg = nullptr; - site->saddle = nullptr; - site->xcp = nullptr; - site->ycp = nullptr; - site->kcp = nullptr; - site->x = nullptr; - site->y = nullptr; - site->z = nullptr; - return site; -} - -void -cntr_init(Csite *site, long iMax, long jMax, const double *x, const double *y, - const double *z, const bool *mask, long i_chunk_size, long j_chunk_size) -{ - long ijmax = iMax * jMax; - long nreg = iMax * jMax + iMax + 1; - - site->imax = iMax; - site->jmax = jMax; - site->data = new Cdata[nreg]; - site->saddle = new Saddle[ijmax]; - if (mask != nullptr) - { - site->reg = new char[nreg]; - mask_zones(iMax, jMax, mask, site->reg); - } - /* I don't think we need to initialize site->data. */ - site->x = x; - site->y = y; - site->z = z; - site->xcp = nullptr; - site->ycp = nullptr; - site->kcp = nullptr; - - /* Store correct chunk sizes for filled contours. Chunking not used for - line contours. */ - if (i_chunk_size <= 0 || i_chunk_size > iMax - 1) - i_chunk_size = iMax - 1; - site->i_chunk_size = i_chunk_size; - - if (j_chunk_size <= 0 || j_chunk_size > jMax - 1) - j_chunk_size = jMax - 1; - site->j_chunk_size = j_chunk_size; -} - -void cntr_del(Csite *site) -{ - delete [] site->saddle; - delete [] site->reg; - delete [] site->data; - delete site; - site = nullptr; -} - -static int -reorder(double *xpp, double *ypp, short *kpp, double *xy, unsigned char *c, int npts, int nlevels) -{ - std::vector<int> subp; - int isp, nsp; - int iseg, nsegs; - int isegplus; - int i; - int k; - int started; - int maxnsegs = npts/2 + 1; - /* allocate maximum possible size--gross overkill */ - std::vector<int> i0(maxnsegs); - std::vector<int> i1(maxnsegs); - - /* Find the segments. */ - iseg = 0; - started = 0; - for (i=0; i<npts; i++) - { - if (started) - { - if ((kpp[i] >= kind_slit_up) || (i == npts-1)) - { - i1[iseg] = i; - started = 0; - iseg++; - if (iseg == maxnsegs) - { - k = -1; - return k; - } - } - } - else if ((kpp[i] < kind_slit_up) && (i < npts-1)) - { - i0[iseg] = i; - started = 1; - } - } - - nsegs = iseg; - - /* Find the subpaths as sets of connected segments. */ - - subp.resize(nsegs, false); - for (i=0; i<nsegs; i++) subp[i] = -1; - - nsp = 0; - for (iseg=0; iseg<nsegs; iseg++) - { - /* For each segment, if it is not closed, look ahead for - the next connected segment. - */ - double xend, yend; - xend = xpp[i1[iseg]]; - yend = ypp[i1[iseg]]; - if (subp[iseg] >= 0) continue; - subp[iseg] = nsp; - nsp++; - if (iseg == nsegs-1) continue; - for (isegplus = iseg+1; isegplus < nsegs; isegplus++) - { - if (subp[isegplus] >= 0) continue; - - if (xend == xpp[i0[isegplus]] && yend == ypp[i0[isegplus]]) - { - subp[isegplus] = subp[iseg]; - xend = xpp[i1[isegplus]]; - yend = ypp[i1[isegplus]]; - } - } - } - - /* Generate the verts and codes from the subpaths. */ - k = 0; - for (isp=0; isp<nsp; isp++) - { - int first = 1; - int kstart = k; - for (iseg=0; iseg<nsegs; iseg++) - { - int istart, iend; - if (subp[iseg] != isp) continue; - iend = i1[iseg]; - if (first) - { - istart = i0[iseg]; - } - else - { - istart = i0[iseg]+1; /* skip duplicate */ - } - for (i=istart; i<=iend; i++) - { - xy[2*k] = xpp[i]; - xy[2*k+1] = ypp[i]; - if (first) c[k] = MOVETO; - else c[k] = LINETO; - first = 0; - k++; - if (k > npts) /* should never happen */ - { - k = -1; - return k; - } - } - } - if (nlevels == 2 || - (xy[2*kstart] == xy[2*k-2] && xy[2*kstart+1] == xy[2*k-1])) - { - c[k-1] = CLOSEPOLY; - } - } - - return k; -} - -/* Build a list of XY 2-D arrays, shape (N,2), to which a list of path - code arrays is concatenated. -*/ -static py::tuple -build_cntr_list_v2(long *np, double *xp, double *yp, short *kp, - int nparts, long ntotal, int nlevels) -{ - int i; - long k; - py::ssize_t dims[2]; - py::ssize_t kdims[1]; - - py::list all_verts(nparts); - py::list all_codes(nparts); - - for (i=0, k=0; i < nparts; k+= np[i], i++) - { - double *xpp = xp+k; - double *ypp = yp+k; - short *kpp = kp+k; - int n; - - dims[0] = np[i]; - dims[1] = 2; - kdims[0] = np[i]; - PointArray xyv(dims); - CodeArray kv(kdims); - - n = reorder(xpp, ypp, kpp, xyv.mutable_data(), kv.mutable_data(), np[i], nlevels); - if (n == -1) - { - throw std::runtime_error("Error reordering vertices"); - } - - dims[0] = n; - xyv.resize(dims, false); - all_verts[i] = xyv; - - kdims[0] = n; - kv.resize(kdims, false); - all_codes[i] = kv; - } - - return py::make_tuple(all_verts, all_codes); -} - -/* cntr_trace is called once per contour level or level pair. - If nlevels is 1, a set of contour lines will be returned; if nlevels - is 2, the set of polygons bounded by the levels will be returned. - If points is True, the lines will be returned as a list of list - of points; otherwise, as a list of tuples of vectors. -*/ - -py::tuple -cntr_trace(Csite *site, double levels[], int nlevels) -{ - int iseg; - - long n; - long nparts = 0; - long ntotal = 0; - long ntotal2 = 0; - - site->zlevel[0] = levels[0]; - site->zlevel[1] = levels[0]; - if (nlevels == 2) - { - site->zlevel[1] = levels[1]; - } - site->n = site->count = 0; - data_init (site); - - /* make first pass to compute required sizes for second pass */ - for (;;) - { - n = curve_tracer (site, 0); - - if (!n) - break; - if (n > 0) - { - nparts++; - ntotal += n; - } - else - { - ntotal -= n; - } - } - std::vector<double> xp0(ntotal); - std::vector<double> yp0(ntotal); - std::vector<short> kp0(ntotal); - std::vector<long> nseg0(nparts); - - /* second pass */ - site->xcp = xp0.data(); - site->ycp = yp0.data(); - site->kcp = kp0.data(); - iseg = 0; - for (;;iseg++) - { - n = curve_tracer (site, 1); - if (ntotal2 + n > ntotal) - { - throw std::runtime_error("curve_tracer: ntotal2, pass 2 exceeds ntotal, pass 1"); - } - if (n == 0) - break; - if (n > 0) - { - /* could add array bounds checking */ - nseg0[iseg] = n; - site->xcp += n; - site->ycp += n; - site->kcp += n; - ntotal2 += n; - } - else - { - throw std::runtime_error("Negative n from curve_tracer in pass 2"); - } - } - - site->xcp = nullptr; - site->ycp = nullptr; - site->kcp = nullptr; - - return build_cntr_list_v2( - nseg0.data(), xp0.data(), yp0.data(), kp0.data(), nparts, ntotal, nlevels); -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/mpl2005_original.h b/contrib/python/contourpy/src/mpl2005_original.h deleted file mode 100644 index 93dba45ff9f..00000000000 --- a/contrib/python/contourpy/src/mpl2005_original.h +++ /dev/null @@ -1,56 +0,0 @@ -#ifndef CONTOURPY_MPL_2005_ORIGINAL_H -#define CONTOURPY_MPL_2005_ORIGINAL_H - -#include "common.h" - -namespace contourpy { - -/* the data about edges, zones, and points -- boundary or not, exists - * or not, z value 0, 1, or 2 -- is kept in a mesh sized data array */ -typedef short Cdata; - -/* information to decide on correct contour direction in saddle zones - * is stored in a mesh sized array. Only those entries corresponding - * to saddle zones have nonzero values in this array. */ -typedef char Saddle; - -/* here is the minimum structure required to tell where we are in the - * mesh sized data array */ -struct Csite -{ - long edge; /* ij of current edge */ - long left; /* +-1 or +-imax as the zone is to right, left, below, - * or above the edge */ - long imax; /* imax for the mesh */ - long jmax; /* jmax for the mesh */ - long n; /* number of points marked on this curve so far */ - long count; /* count of start markers visited */ - double zlevel[2]; /* contour levels, zlevel[1]<=zlevel[0] - * signals single level case */ - Saddle *saddle; /* saddle zone information for the mesh */ - char *reg; /* region array for the mesh (was int) */ - Cdata *data; /* added by EF */ - long edge0, left0; /* starting site on this curve for closure */ - int level0; /* starting level for closure */ - long edge00; /* site needing START_ROW mark */ - - /* making the actual marks requires a bunch of other stuff */ - const double *x, *y, *z; /* mesh coordinates and function values */ - double *xcp, *ycp; /* output contour points */ - short *kcp; /* kind of contour point */ - - long i_chunk_size, j_chunk_size; -}; - -Csite *cntr_new(); - -void cntr_init(Csite *site, long iMax, long jMax, const double *x, const double *y, - const double *z, const bool *mask, long i_chunk_size, long j_chunk_size); - -void cntr_del(Csite *site); - -py::tuple cntr_trace(Csite *site, double levels[], int nlevels); - -} // namespace contourpy - -#endif // CONTOURPY_MPL_2005_ORIGINAL_H diff --git a/contrib/python/contourpy/src/mpl2014.cpp b/contrib/python/contourpy/src/mpl2014.cpp deleted file mode 100644 index e1021572a7f..00000000000 --- a/contrib/python/contourpy/src/mpl2014.cpp +++ /dev/null @@ -1,1710 +0,0 @@ -// 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 setup.py, and then rebuilding. - -#include "mpl2014.h" -#include "mpl_kind_code.h" -#include <algorithm> - -namespace contourpy { -namespace mpl2014 { - -// 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(index_t quad_, Edge edge_) - : quad(quad_), edge(edge_) -{} - -bool QuadEdge::operator==(const QuadEdge& other) const -{ - return quad == other.quad && edge == other.edge; -} - -std::ostream& operator<<(std::ostream& os, const QuadEdge& quad_edge) -{ - return os << quad_edge.quad << ' ' << quad_edge.edge; -} - - - -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; -} - -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::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(index_t nx, index_t x_chunk_points, index_t 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(index_t quad) -{ - index_t index = index_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; -} - -index_t ParentCache::index_to_index(index_t quad) const -{ - index_t i = quad % _nx; - index_t j = quad / _nx; - index_t 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<index_t>(_lines.size()) && - "ParentCache index outside chunk"); - - return index; -} - -void ParentCache::set_chunk_starts(index_t istart, index_t 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(index_t quad, ContourLine& contour_line) -{ - assert(!_lines.empty() && "Accessing ParentCache before it has been initialised"); - index_t index = index_to_index(quad); - if (_lines[index] == 0) - _lines[index] = contour_line.is_hole() ? contour_line.get_parent() : &contour_line; -} - - - -Mpl2014ContourGenerator::Mpl2014ContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, bool corner_mask, index_t x_chunk_size, index_t y_chunk_size) - : _x(x), - _y(y), - _z(z), - _nx(_z.ndim() > 1 ? _z.shape(1) : 0), - _ny(_z.ndim() > 0 ? _z.shape(0) : 0), - _n(_nx*_ny), - _corner_mask(corner_mask), - _x_chunk_size(calc_chunk_size(_nx, x_chunk_size)), - _y_chunk_size(calc_chunk_size(_ny, y_chunk_size)), - _nxchunk(calc_chunk_count(_nx, _x_chunk_size)), - _nychunk(calc_chunk_count(_ny, _y_chunk_size)), - _chunk_count(_nxchunk*_nychunk), - _cache(new CacheItem[_n]), - _parent_cache(_nx, - _x_chunk_size > 0 ? _x_chunk_size+1 : _nx, - _y_chunk_size > 0 ? _y_chunk_size+1 : _ny) -{ - if (_x.ndim() != 2 || _y.ndim() != 2 || _z.ndim() != 2) - throw std::invalid_argument("x, y and z must all be 2D arrays"); - - if (_x.shape(1) != _nx || _x.shape(0) != _ny || - _y.shape(1) != _nx || _y.shape(0) != _ny) - throw std::invalid_argument("x, y and z arrays must have the same shape"); - - if (_nx < 2 || _ny < 2) - throw std::invalid_argument("x, y and z must all be at least 2x2 arrays"); - - if (mask.ndim() != 0) { // ndim == 0 if mask is not set, which is valid. - if (mask.ndim() != 2) - throw std::invalid_argument("mask array must be a 2D array"); - - if (mask.shape(1) != _nx || mask.shape(0) != _ny) - throw std::invalid_argument( - "If mask is set it must be a 2D array with the same shape as z"); - } - - if (x_chunk_size < 0 || y_chunk_size < 0) - throw std::invalid_argument("x_chunk_size and y_chunk_size cannot be negative"); - - init_cache_grid(mask); -} - -Mpl2014ContourGenerator::~Mpl2014ContourGenerator() -{ - delete [] _cache; -} - -void Mpl2014ContourGenerator::append_contour_line_to_vertices_and_codes( - ContourLine& contour_line, py::list& vertices_list, py::list& codes_list) const -{ - // Convert ContourLine to Python equivalent, and clear it for reuse. - // This function is called once for each line generated in create_contour(). - // A line is either a closed line loop (in which case the last point is - // identical to the first) or an open line strip. Two NumPy arrays are - // created for each line: - // vertices is a double array of shape (npoints, 2) containing the (x, y) - // coordinates of the points in the line - // codes is a uint8 array of shape (npoints,) containing the 'kind codes' - // which are defined in the Path class - // and they are appended to the Python lists vertices_list and codes_list - // respectively for return to the Python calling function. - - py::ssize_t npoints = static_cast<py::ssize_t>(contour_line.size()); - - py::ssize_t vertices_dims[2] = {npoints, 2}; - PointArray vertices(vertices_dims); - double* vertices_ptr = vertices.mutable_data(); - - py::ssize_t codes_dims[1] = {npoints}; - CodeArray codes(codes_dims); - unsigned char* codes_ptr = codes.mutable_data(); - - for (ContourLine::const_iterator point = contour_line.begin(); - point != contour_line.end(); ++point) { - *vertices_ptr++ = point->x; - *vertices_ptr++ = point->y; - *codes_ptr++ = (point == contour_line.begin() ? MOVETO : LINETO); - } - - // Closed line loop has identical first and last (x, y) points. - if (contour_line.size() > 1 && contour_line.front() == contour_line.back()) - *(codes_ptr-1) = CLOSEPOLY; - - vertices_list.append(vertices); - codes_list.append(codes); - - contour_line.clear(); -} - -void Mpl2014ContourGenerator::append_contour_to_vertices_and_codes( - Contour& contour, py::list& vertices_list, py::list& codes_list) const -{ - // Convert Contour to Python equivalent, and clear it for reuse. - // This function is called once for each polygon generated in - // create_filled_contour(). A polygon consists of an outer line loop - // (called the parent) and zero or more inner line loops or holes (called - // the children). Two NumPy arrays are created for each polygon: - // vertices is a double array of shape (npoints, 2) containing the (x, y) - // coordinates of the points in the polygon (parent plus children) - // codes is a uint8 array of shape (npoints,) containing the 'kind codes' - // which are defined in the Path class - // and they are appended to the Python lists vertices_list and codes_list - // respectively for return to the Python calling function. - - // Convert Contour to python equivalent, and clear it. - 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(); - py::ssize_t npoints = static_cast<py::ssize_t>(line.size() + 1); - for (children_it = children.begin(); children_it != children.end(); ++children_it) - // cppcheck-suppress useStlAlgorithm - npoints += static_cast<py::ssize_t>((*children_it)->size() + 1); - - py::ssize_t vertices_dims[2] = {npoints, 2}; - PointArray vertices(vertices_dims); - double* vertices_ptr = vertices.mutable_data(); - - py::ssize_t codes_dims[1] = {npoints}; - CodeArray codes(codes_dims); - unsigned char* codes_ptr = codes.mutable_data(); - - for (point = line.begin(); point != line.end(); ++point) { - *vertices_ptr++ = point->x; - *vertices_ptr++ = point->y; - *codes_ptr++ = (point == line.begin() ? MOVETO : LINETO); - } - point = line.begin(); - *vertices_ptr++ = point->x; - *vertices_ptr++ = point->y; - *codes_ptr++ = CLOSEPOLY; - - for (children_it = children.begin(); children_it != children.end(); ++children_it) { - ContourLine& child = **children_it; - for (point = child.begin(); point != child.end(); ++point) { - *vertices_ptr++ = point->x; - *vertices_ptr++ = point->y; - *codes_ptr++ = (point == child.begin() ? MOVETO : LINETO); - } - point = child.begin(); - *vertices_ptr++ = point->x; - *vertices_ptr++ = point->y; - *codes_ptr++ = CLOSEPOLY; - - child.clear_parent(); // To indicate it can be deleted. - } - - vertices_list.append(vertices); - codes_list.append(codes); - - delete *line_it; - *line_it = 0; - } - } - - // Delete remaining contour lines. - contour.delete_contour_lines(); -} - -index_t Mpl2014ContourGenerator::calc_chunk_count(index_t point_count, index_t chunk_size) -{ - // Accepts any values for point_count and chunk_size. - - if (chunk_size > 0 && point_count > 1) { - index_t 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; -} - -index_t Mpl2014ContourGenerator::calc_chunk_size(index_t point_count, index_t chunk_size) -{ - // Accepts any values for point_count and chunk_size. - index_t ret = point_count - 1; - if (chunk_size > 0) - ret = std::min(chunk_size, ret); - - return std::max<index_t>(ret, 1); -} - -void Mpl2014ContourGenerator::edge_interp( - const QuadEdge& quad_edge, const double& level, ContourLine& contour_line) -{ - assert(quad_edge.quad >= 0 && quad_edge.quad < _n && "Quad index out of bounds"); - assert(quad_edge.edge != Edge_None && "Invalid edge"); - interp(get_edge_point_index(quad_edge, true), - get_edge_point_index(quad_edge, false), - level, contour_line); -} - -py::tuple Mpl2014ContourGenerator::filled(double lower_level, double upper_level) -{ - check_levels(lower_level, upper_level); - - init_cache_levels(lower_level, upper_level); - - Contour contour; - py::list vertices, codes; - - index_t ichunk, jchunk, istart, iend, jstart, jend; - for (index_t ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) { - get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart, jend); - _parent_cache.set_chunk_starts(istart, jstart); - - for (index_t j = jstart; j < jend; ++j) { - index_t quad_end = iend + j*_nx; - for (index_t 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) { - index_t quad_end = iend + jend*_nx; - for (index_t quad = istart + jend*_nx; quad < quad_end; ++quad) - _cache[quad] &= ~MASK_VISITED_S; - } - - if (ichunk < _nxchunk-1) { - index_t quad_end = iend + jend*_nx; - for (index_t 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 py::make_tuple(vertices, codes); -} - -unsigned int Mpl2014ContourGenerator::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; - index_t& 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); - index_t 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. - edge_interp(quad_edge, level_index == 1 ? lower_level : upper_level, contour_line); - 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. - get_point_xy(end_point, contour_line); - - first_edge = false; - } - - return level_index; -} - -void Mpl2014ContourGenerator::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"); - - index_t& quad = quad_edge.quad; - Edge& edge = quad_edge.edge; - - if (want_initial_point) - edge_interp(quad_edge, level, contour_line); - - 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. - index_t 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. - index_t 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. - edge_interp(quad_edge, level, contour_line); - - // 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; - } -} - -py::tuple Mpl2014ContourGenerator::get_chunk_count() const -{ - return py::make_tuple(_nychunk, _nxchunk); -} - -void Mpl2014ContourGenerator::get_chunk_limits( - index_t ijchunk, index_t& ichunk, index_t& jchunk, index_t& istart, index_t& iend, - index_t& jstart, index_t& jend) -{ - assert(ijchunk >= 0 && ijchunk < _chunk_count && "ijchunk out of bounds"); - ichunk = ijchunk % _nxchunk; - jchunk = ijchunk / _nxchunk; - istart = ichunk*_x_chunk_size; - iend = (ichunk == _nxchunk-1 ? _nx : (ichunk+1)*_x_chunk_size); - jstart = jchunk*_y_chunk_size; - jend = (jchunk == _nychunk-1 ? _ny : (jchunk+1)*_y_chunk_size); -} - -py::tuple Mpl2014ContourGenerator::get_chunk_size() const -{ - return py::make_tuple(_y_chunk_size, _x_chunk_size); -} - -bool Mpl2014ContourGenerator::get_corner_mask() const -{ - return _corner_mask; -} - -Edge Mpl2014ContourGenerator::get_corner_start_edge( - index_t 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 - // point1 +---------+ point2 - // \ | - // \ | edge23 - // edge31 \ | - // \ | - // + point3 - // - index_t point1, point2, point3; - Edge edge12, edge23, edge31; - switch (_cache[quad] & MASK_EXISTS) { - case MASK_EXISTS_SW_CORNER: - point1 = POINT_SE; point2 = POINT_SW; point3 = POINT_NW; - edge12 = Edge_S; edge23 = Edge_W; edge31 = Edge_NE; - break; - case MASK_EXISTS_SE_CORNER: - point1 = POINT_NE; point2 = POINT_SE; point3 = POINT_SW; - edge12 = Edge_E; edge23 = Edge_S; edge31 = Edge_NW; - break; - case MASK_EXISTS_NW_CORNER: - point1 = POINT_SW; point2 = POINT_NW; point3 = POINT_NE; - edge12 = Edge_W; edge23 = Edge_N; edge31 = Edge_SE; - break; - case MASK_EXISTS_NE_CORNER: - point1 = POINT_NW; point2 = POINT_NE; point3 = POINT_SE; - edge12 = Edge_N; edge23 = Edge_E; edge31 = Edge_SW; - break; - default: - assert(0 && "Invalid EXISTS for quad"); - return Edge_None; - } - - unsigned int config = ((Z_LEVEL(point1) >= level_index) << 2) | - ((Z_LEVEL(point2) >= level_index) << 1) | - ((Z_LEVEL(point3) >= level_index) << 0); - - // 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 = 7 - config; - - switch (config) { - case 0: return Edge_None; - case 1: return edge23; - case 2: return edge12; - case 3: return edge12; - case 4: return edge31; - case 5: return edge23; - case 6: return edge31; - case 7: return Edge_None; - default: assert(0 && "Invalid config"); return Edge_None; - } -} - -index_t Mpl2014ContourGenerator::get_edge_point_index( - const QuadEdge& quad_edge, bool start) const -{ - assert(quad_edge.quad >= 0 && quad_edge.quad < _n && "Quad index out of bounds"); - assert(quad_edge.edge != Edge_None && "Invalid edge"); - - // Edges are ordered anticlockwise around their quad, as indicated by - // directions of arrows in diagrams below. - // Full quad NW corner (others similar) - // - // POINT_NW Edge_N POINT_NE POINT_NW Edge_N POINT_NE - // +----<-----+ +----<-----+ - // | | | / - // | | | quad / - // Edge_W V quad ^ Edge_E Edge_W V ^ - // | | | / Edge_SE - // | | | / - // +---->-----+ + - // POINT_SW Edge_S POINT_SE POINT_SW - // - const index_t& quad = quad_edge.quad; - switch (quad_edge.edge) { - case Edge_E: return (start ? POINT_SE : POINT_NE); - case Edge_N: return (start ? POINT_NE : POINT_NW); - case Edge_W: return (start ? POINT_NW : POINT_SW); - case Edge_S: return (start ? POINT_SW : POINT_SE); - case Edge_NE: return (start ? POINT_SE : POINT_NW); - case Edge_NW: return (start ? POINT_NE : POINT_SW); - case Edge_SW: return (start ? POINT_NW : POINT_SE); - case Edge_SE: return (start ? POINT_SW : POINT_NE); - default: assert(0 && "Invalid edge"); return 0; - } -} - -Edge Mpl2014ContourGenerator::get_exit_edge(const QuadEdge& quad_edge, Dir dir) const -{ - assert(quad_edge.quad >= 0 && quad_edge.quad < _n && "Quad index out of bounds"); - assert(quad_edge.edge != Edge_None && "Invalid edge"); - - const index_t& quad = quad_edge.quad; - const Edge& edge = quad_edge.edge; - if (EXISTS_ANY_CORNER(quad)) { - // Corner directions are always left or right. A corner is a triangle, - // entered via one edge so the other two edges are the left and right - // ones. - switch (edge) { - case Edge_E: - return (EXISTS_SE_CORNER(quad) - ? (dir == Dir_Left ? Edge_S : Edge_NW) - : (dir == Dir_Right ? Edge_N : Edge_SW)); - case Edge_N: - return (EXISTS_NW_CORNER(quad) - ? (dir == Dir_Right ? Edge_W : Edge_SE) - : (dir == Dir_Left ? Edge_E : Edge_SW)); - case Edge_W: - return (EXISTS_SW_CORNER(quad) - ? (dir == Dir_Right ? Edge_S : Edge_NE) - : (dir == Dir_Left ? Edge_N : Edge_SE)); - case Edge_S: - return (EXISTS_SW_CORNER(quad) - ? (dir == Dir_Left ? Edge_W : Edge_NE) - : (dir == Dir_Right ? Edge_E : Edge_NW)); - case Edge_NE: return (dir == Dir_Left ? Edge_S : Edge_W); - case Edge_NW: return (dir == Dir_Left ? Edge_E : Edge_S); - case Edge_SW: return (dir == Dir_Left ? Edge_N : Edge_E); - case Edge_SE: return (dir == Dir_Left ? Edge_W : Edge_N); - default: assert(0 && "Invalid edge"); return Edge_None; - } - } - else { - // A full quad has four edges, entered via one edge so that other three - // edges correspond to left, straight and right directions. - switch (edge) { - case Edge_E: - return (dir == Dir_Left ? Edge_S : (dir == Dir_Right ? Edge_N : Edge_W)); - case Edge_N: - return (dir == Dir_Left ? Edge_E : (dir == Dir_Right ? Edge_W : Edge_S)); - case Edge_W: - return (dir == Dir_Left ? Edge_N : (dir == Dir_Right ? Edge_S : Edge_E)); - case Edge_S: - return (dir == Dir_Left ? Edge_W : (dir == Dir_Right ? Edge_E : Edge_N)); - default: assert(0 && "Invalid edge"); return Edge_None; - } - } -} - -const double& Mpl2014ContourGenerator::get_point_x(index_t point) const -{ - assert(point >= 0 && point < _n && "Point index out of bounds."); - return _x.data()[point]; -} - -const double& Mpl2014ContourGenerator::get_point_y(index_t point) const -{ - assert(point >= 0 && point < _n && "Point index out of bounds."); - return _y.data()[point]; -} - -void Mpl2014ContourGenerator::get_point_xy( - index_t point, ContourLine& contour_line) const -{ - assert(point >= 0 && point < _n && "Point index out of bounds."); - contour_line.emplace_back(_x.data()[point], _y.data()[point]); -} - -const double& Mpl2014ContourGenerator::get_point_z(index_t point) const -{ - assert(point >= 0 && point < _n && "Point index out of bounds."); - return _z.data()[point]; -} - -Edge Mpl2014ContourGenerator::get_quad_start_edge( - index_t 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_QUAD(quad) && "Quad does not exist"); - - unsigned int config = ((Z_NW >= level_index) << 3) | - ((Z_NE >= level_index) << 2) | - ((Z_SW >= level_index) << 1) | - ((Z_SE >= level_index) << 0); - - // 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 = 15 - config; - - switch (config) { - case 0: return Edge_None; - case 1: return Edge_E; - case 2: return Edge_S; - case 3: return Edge_E; - case 4: return Edge_N; - case 5: return Edge_N; - case 6: - // If already identified as a saddle quad then the start edge is - // read from the cache. Otherwise return either valid start edge - // and the subsequent call to follow_interior() will correctly set - // up saddle bits in cache. - if (!SADDLE(quad,level_index) || SADDLE_START_SW(quad,level_index)) - return Edge_S; - else - return Edge_N; - case 7: return Edge_N; - case 8: return Edge_W; - case 9: - // See comment for 6 above. - if (!SADDLE(quad,level_index) || SADDLE_START_SW(quad,level_index)) - return Edge_W; - else - return Edge_E; - case 10: return Edge_S; - case 11: return Edge_E; - case 12: return Edge_W; - case 13: return Edge_W; - case 14: return Edge_S; - case 15: return Edge_None; - default: assert(0 && "Invalid config"); return Edge_None; - } -} - -Edge Mpl2014ContourGenerator::get_start_edge(index_t quad, unsigned int level_index) const -{ - if (EXISTS_ANY_CORNER(quad)) - return get_corner_start_edge(quad, level_index); - else - return get_quad_start_edge(quad, level_index); -} - -void Mpl2014ContourGenerator::init_cache_grid(const MaskArray& mask) -{ - index_t i, j, quad; - - if (mask.ndim() == 0) { - // No mask, easy to calculate quad existence and boundaries together. - quad = 0; - for (j = 0; j < _ny; ++j) { - for (i = 0; i < _nx; ++i, ++quad) { - _cache[quad] = 0; - - if (i < _nx-1 && j < _ny-1) - _cache[quad] |= MASK_EXISTS_QUAD; - - if ((i % _x_chunk_size == 0 || i == _nx-1) && j < _ny-1) - _cache[quad] |= MASK_BOUNDARY_W; - - if ((j % _y_chunk_size == 0 || j == _ny-1) && i < _nx-1) - _cache[quad] |= MASK_BOUNDARY_S; - } - } - } - else { - const bool* mask_ptr = mask.data(); - - // Have mask so use two stages. - // Stage 1, determine if quads/corners exist. - quad = 0; - for (j = 0; j < _ny; ++j) { - for (i = 0; i < _nx; ++i, ++quad) { - _cache[quad] = 0; - - if (i < _nx-1 && j < _ny-1) { - unsigned int config = (mask_ptr[POINT_NW] << 3) | - (mask_ptr[POINT_NE] << 2) | - (mask_ptr[POINT_SW] << 1) | - (mask_ptr[POINT_SE] << 0); - - if (_corner_mask) { - switch (config) { - case 0: _cache[quad] = MASK_EXISTS_QUAD; break; - case 1: _cache[quad] = MASK_EXISTS_NW_CORNER; break; - case 2: _cache[quad] = MASK_EXISTS_NE_CORNER; break; - case 4: _cache[quad] = MASK_EXISTS_SW_CORNER; break; - case 8: _cache[quad] = MASK_EXISTS_SE_CORNER; break; - default: - // Do nothing, quad is masked out. - break; - } - } - else if (config == 0) - _cache[quad] = MASK_EXISTS_QUAD; - } - } - } - - // Stage 2, calculate W and S boundaries. For each quad use boundary - // data already calculated for quads to W and S, so must iterate - // through quads in correct order (increasing i and j indices). - // Cannot use boundary data for quads to E and N as have not yet - // calculated it. - quad = 0; - for (j = 0; j < _ny; ++j) { - for (i = 0; i < _nx; ++i, ++quad) { - if (_corner_mask) { - bool W_exists_none = (i == 0 || EXISTS_NONE(quad-1)); - bool S_exists_none = (j == 0 || EXISTS_NONE(quad-_nx)); - bool W_exists_E_edge = (i > 0 && EXISTS_E_EDGE(quad-1)); - bool S_exists_N_edge = (j > 0 && EXISTS_N_EDGE(quad-_nx)); - - if ((EXISTS_W_EDGE(quad) && W_exists_none) || - (EXISTS_NONE(quad) && W_exists_E_edge) || - (i % _x_chunk_size == 0 && EXISTS_W_EDGE(quad) && W_exists_E_edge)) - _cache[quad] |= MASK_BOUNDARY_W; - - if ((EXISTS_S_EDGE(quad) && S_exists_none) || - (EXISTS_NONE(quad) && S_exists_N_edge) || - (j % _y_chunk_size == 0 && EXISTS_S_EDGE(quad) && S_exists_N_edge)) - _cache[quad] |= MASK_BOUNDARY_S; - } - else { - bool W_exists_quad = (i > 0 && EXISTS_QUAD(quad-1)); - bool S_exists_quad = (j > 0 && EXISTS_QUAD(quad-_nx)); - - if ((EXISTS_QUAD(quad) != W_exists_quad) || - (i % _x_chunk_size == 0 && EXISTS_QUAD(quad) && W_exists_quad)) - _cache[quad] |= MASK_BOUNDARY_W; - - if ((EXISTS_QUAD(quad) != S_exists_quad) || - (j % _y_chunk_size == 0 && EXISTS_QUAD(quad) && S_exists_quad)) - _cache[quad] |= MASK_BOUNDARY_S; - } - } - } - } -} - -void Mpl2014ContourGenerator::init_cache_levels( - const double& lower_level, const double& upper_level) -{ - assert(!(upper_level < lower_level) && "upper and lower levels are wrong way round"); - - bool two_levels = (lower_level != upper_level); - CacheItem keep_mask = - (_corner_mask ? MASK_EXISTS | MASK_BOUNDARY_S | MASK_BOUNDARY_W - : MASK_EXISTS_QUAD | MASK_BOUNDARY_S | MASK_BOUNDARY_W); - - if (two_levels) { - const double* z_ptr = _z.data(); - for (index_t quad = 0; quad < _n; ++quad, ++z_ptr) { - _cache[quad] &= keep_mask; - if (*z_ptr > upper_level) - _cache[quad] |= MASK_Z_LEVEL_2; - else if (*z_ptr > lower_level) - _cache[quad] |= MASK_Z_LEVEL_1; - } - } - else { - const double* z_ptr = _z.data(); - for (index_t quad = 0; quad < _n; ++quad, ++z_ptr) { - _cache[quad] &= keep_mask; - if (*z_ptr > lower_level) - _cache[quad] |= MASK_Z_LEVEL_1; - } - } -} - -void Mpl2014ContourGenerator::interp( - index_t point1, index_t point2, const double& level, ContourLine& contour_line) const -{ - assert(point1 >= 0 && point1 < _n && "Point index 1 out of bounds."); - assert(point2 >= 0 && point2 < _n && "Point index 2 out of bounds."); - assert(point1 != point2 && "Identical points"); - double fraction = (get_point_z(point2) - level) / (get_point_z(point2) - get_point_z(point1)); - contour_line.emplace_back( - get_point_x(point1)*fraction + get_point_x(point2)*(1.0 - fraction), - get_point_y(point1)*fraction + get_point_y(point2)*(1.0 - fraction)); -} - -bool Mpl2014ContourGenerator::is_edge_a_boundary( - const QuadEdge& quad_edge) const -{ - assert(quad_edge.quad >= 0 && quad_edge.quad < _n && "Quad index out of bounds"); - assert(quad_edge.edge != Edge_None && "Invalid edge"); - - switch (quad_edge.edge) { - case Edge_E: return BOUNDARY_E(quad_edge.quad); - case Edge_N: return BOUNDARY_N(quad_edge.quad); - case Edge_W: return BOUNDARY_W(quad_edge.quad); - case Edge_S: return BOUNDARY_S(quad_edge.quad); - case Edge_NE: return EXISTS_SW_CORNER(quad_edge.quad); - case Edge_NW: return EXISTS_SE_CORNER(quad_edge.quad); - case Edge_SW: return EXISTS_NE_CORNER(quad_edge.quad); - case Edge_SE: return EXISTS_NW_CORNER(quad_edge.quad); - default: assert(0 && "Invalid edge"); return true; - } -} - -py::sequence Mpl2014ContourGenerator::lines(double level) -{ - init_cache_levels(level, level); - - py::list vertices_list, codes_list; - - // Lines that start and end on boundaries. - index_t ichunk, jchunk, istart, iend, jstart, jend; - for (index_t ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) { - get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart, jend); - - for (index_t j = jstart; j < jend; ++j) { - index_t quad_end = iend + j*_nx; - for (index_t 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, codes_list, quad, Edge_S, level)) continue; - - if (BOUNDARY_W(quad) && Z_NW >= 1 && Z_SW < 1 && - start_line(vertices_list, codes_list, quad, Edge_W, level)) continue; - - if (BOUNDARY_N(quad) && Z_NE >= 1 && Z_NW < 1 && - start_line(vertices_list, codes_list, quad, Edge_N, level)) continue; - - if (BOUNDARY_E(quad) && Z_SE >= 1 && Z_NE < 1 && - start_line(vertices_list, codes_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, codes_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, codes_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, codes_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, codes_list, quad, Edge_SW, level)) continue; - } - } - } - } - - // Internal loops. - ContourLine contour_line(false); // Reused for each contour line. - for (index_t ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) { - get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart, jend); - - for (index_t j = jstart; j < jend; ++j) { - index_t quad_end = iend + j*_nx; - for (index_t 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_and_codes(contour_line, vertices_list, codes_list); - - // Repeat if saddle point but not visited. - if (SADDLE(quad,1) && !VISITED(quad,1)) - --quad; - } - } - } - - return py::make_tuple(vertices_list, codes_list); -} - -void Mpl2014ContourGenerator::move_to_next_boundary_edge( - QuadEdge& quad_edge) const -{ - assert(is_edge_a_boundary(quad_edge) && "QuadEdge is not a boundary"); - - index_t& quad = quad_edge.quad; - Edge& edge = quad_edge.edge; - - quad = get_edge_point_index(quad_edge, false); - - // quad is now such that POINT_SW is the end point of the quad_edge passed - // to this function. - - // To find the next boundary edge, first attempt to turn left 135 degrees - // and if that edge is a boundary then move to it. If not, attempt to turn - // left 90 degrees, then left 45 degrees, then straight on, etc, until can - // move. - // First determine which edge to attempt first. - int index = 0; - switch (edge) { - case Edge_E: index = 0; break; - case Edge_SE: index = 1; break; - case Edge_S: index = 2; break; - case Edge_SW: index = 3; break; - case Edge_W: index = 4; break; - case Edge_NW: index = 5; break; - case Edge_N: index = 6; break; - case Edge_NE: index = 7; break; - default: assert(0 && "Invalid edge"); break; - } - - // If _corner_mask not set, only need to consider odd index in loop below. - if (!_corner_mask) - ++index; - - // Try each edge in turn until a boundary is found. - int start_index = index; - do - { - switch (index) { - case 0: - if (EXISTS_SE_CORNER(quad-_nx-1)) { // Equivalent to BOUNDARY_NW - quad -= _nx+1; - edge = Edge_NW; - return; - } - break; - case 1: - if (BOUNDARY_N(quad-_nx-1)) { - quad -= _nx+1; - edge = Edge_N; - return; - } - break; - case 2: - if (EXISTS_SW_CORNER(quad-1)) { // Equivalent to BOUNDARY_NE - quad -= 1; - edge = Edge_NE; - return; - } - break; - case 3: - if (BOUNDARY_E(quad-1)) { - quad -= 1; - edge = Edge_E; - return; - } - break; - case 4: - if (EXISTS_NW_CORNER(quad)) { // Equivalent to BOUNDARY_SE - edge = Edge_SE; - return; - } - break; - case 5: - if (BOUNDARY_S(quad)) { - edge = Edge_S; - return; - } - break; - case 6: - if (EXISTS_NE_CORNER(quad-_nx)) { // Equivalent to BOUNDARY_SW - quad -= _nx; - edge = Edge_SW; - return; - } - break; - case 7: - if (BOUNDARY_W(quad-_nx)) { - quad -= _nx; - edge = Edge_W; - return; - } - break; - default: assert(0 && "Invalid index"); break; - } - - if (_corner_mask) - index = (index + 1) % 8; - else - index = (index + 2) % 8; - } while (index != start_index); - - assert(0 && "Failed to find next boundary edge"); -} - -void Mpl2014ContourGenerator::move_to_next_quad(QuadEdge& quad_edge) const -{ - assert(quad_edge.quad >= 0 && quad_edge.quad < _n && "Quad index out of bounds"); - assert(quad_edge.edge != Edge_None && "Invalid edge"); - - // Move from quad_edge.quad to the neighbouring quad in the direction - // specified by quad_edge.edge. - switch (quad_edge.edge) { - case Edge_E: quad_edge.quad += 1; quad_edge.edge = Edge_W; break; - case Edge_N: quad_edge.quad += _nx; quad_edge.edge = Edge_S; break; - case Edge_W: quad_edge.quad -= 1; quad_edge.edge = Edge_E; break; - case Edge_S: quad_edge.quad -= _nx; quad_edge.edge = Edge_N; break; - default: assert(0 && "Invalid edge"); break; - } -} - -void Mpl2014ContourGenerator::single_quad_filled( - Contour& contour, index_t quad, const double& lower_level, const double& upper_level) -{ - assert(quad >= 0 && quad < _n && "Quad index out of bounds"); - - // Order of checking is important here as can have different ContourLines - // from both lower and upper levels in the same quad. First check the S - // edge, then move up the quad to the N edge checking as required. - - // Possible starts from S boundary. - if (BOUNDARY_S(quad) && EXISTS_S_EDGE(quad)) { - - // Lower-level start from S boundary into interior. - if (!VISITED_S(quad) && Z_SW >= 1 && Z_SE == 0) - contour.push_back(start_filled( - quad, Edge_S, 1, NotHole, Interior, lower_level, upper_level)); - - // Upper-level start from S boundary into interior. - if (!VISITED_S(quad) && Z_SW < 2 && Z_SE == 2) - contour.push_back(start_filled( - quad, Edge_S, 2, NotHole, Interior, lower_level, upper_level)); - - // Lower-level start following S boundary from W to E. - if (!VISITED_S(quad) && Z_SW <= 1 && Z_SE == 1) - contour.push_back(start_filled( - quad, Edge_S, 1, NotHole, Boundary, lower_level, upper_level)); - - // Upper-level start following S boundary from W to E. - if (!VISITED_S(quad) && Z_SW == 2 && Z_SE == 1) - contour.push_back(start_filled( - quad, Edge_S, 2, NotHole, Boundary, lower_level, upper_level)); - } - - // Possible starts from W boundary. - if (BOUNDARY_W(quad) && EXISTS_W_EDGE(quad)) { - - // Lower-level start from W boundary into interior. - if (!VISITED_W(quad) && Z_NW >= 1 && Z_SW == 0) - contour.push_back(start_filled( - quad, Edge_W, 1, NotHole, Interior, lower_level, upper_level)); - - // Upper-level start from W boundary into interior. - if (!VISITED_W(quad) && Z_NW < 2 && Z_SW == 2) - contour.push_back(start_filled( - quad, Edge_W, 2, NotHole, Interior, lower_level, upper_level)); - - // Lower-level start following W boundary from N to S. - if (!VISITED_W(quad) && Z_NW <= 1 && Z_SW == 1) - contour.push_back(start_filled( - quad, Edge_W, 1, NotHole, Boundary, lower_level, upper_level)); - - // Upper-level start following W boundary from N to S. - if (!VISITED_W(quad) && Z_NW == 2 && Z_SW == 1) - contour.push_back(start_filled( - quad, Edge_W, 2, NotHole, Boundary, lower_level, upper_level)); - } - - // Possible starts from NE boundary. - if (EXISTS_SW_CORNER(quad)) { // i.e. BOUNDARY_NE - - // Lower-level start following NE boundary from SE to NW, hole. - if (!VISITED_CORNER(quad) && Z_NW == 1 && Z_SE == 1) - contour.push_back(start_filled( - quad, Edge_NE, 1, Hole, Boundary, lower_level, upper_level)); - } - // Possible starts from SE boundary. - else if (EXISTS_NW_CORNER(quad)) { // i.e. BOUNDARY_SE - - // Lower-level start from N to SE. - if (!VISITED(quad,1) && Z_NW == 0 && Z_SW == 0 && Z_NE >= 1) - contour.push_back(start_filled( - quad, Edge_N, 1, NotHole, Interior, lower_level, upper_level)); - - // Upper-level start from SE to N, hole. - if (!VISITED(quad,2) && Z_NW < 2 && Z_SW < 2 && Z_NE == 2) - contour.push_back(start_filled( - quad, Edge_SE, 2, Hole, Interior, lower_level, upper_level)); - - // Upper-level start from N to SE. - if (!VISITED(quad,2) && Z_NW == 2 && Z_SW == 2 && Z_NE < 2) - contour.push_back(start_filled( - quad, Edge_N, 2, NotHole, Interior, lower_level, upper_level)); - - // Lower-level start from SE to N, hole. - if (!VISITED(quad,1) && Z_NW >= 1 && Z_SW >= 1 && Z_NE == 0) - contour.push_back(start_filled( - quad, Edge_SE, 1, Hole, Interior, lower_level, upper_level)); - } - // Possible starts from NW boundary. - else if (EXISTS_SE_CORNER(quad)) { // i.e. BOUNDARY_NW - - // Lower-level start from NW to E. - if (!VISITED(quad,1) && Z_SW == 0 && Z_SE == 0 && Z_NE >= 1) - contour.push_back(start_filled( - quad, Edge_NW, 1, NotHole, Interior, lower_level, upper_level)); - - // Upper-level start from E to NW, hole. - if (!VISITED(quad,2) && Z_SW < 2 && Z_SE < 2 && Z_NE == 2) - contour.push_back(start_filled( - quad, Edge_E, 2, Hole, Interior, lower_level, upper_level)); - - // Upper-level start from NW to E. - if (!VISITED(quad,2) && Z_SW == 2 && Z_SE == 2 && Z_NE < 2) - contour.push_back(start_filled( - quad, Edge_NW, 2, NotHole, Interior, lower_level, upper_level)); - - // Lower-level start from E to NW, hole. - if (!VISITED(quad,1) && Z_SW >= 1 && Z_SE >= 1 && Z_NE == 0) - contour.push_back(start_filled( - quad, Edge_E, 1, Hole, Interior, lower_level, upper_level)); - } - // Possible starts from SW boundary. - else if (EXISTS_NE_CORNER(quad)) { // i.e. BOUNDARY_SW - - // Lower-level start from SW boundary into interior. - if (!VISITED_CORNER(quad) && Z_NW >= 1 && Z_SE == 0) - contour.push_back(start_filled( - quad, Edge_SW, 1, NotHole, Interior, lower_level, upper_level)); - - // Upper-level start from SW boundary into interior. - if (!VISITED_CORNER(quad) && Z_NW < 2 && Z_SE == 2) - contour.push_back(start_filled( - quad, Edge_SW, 2, NotHole, Interior, lower_level, upper_level)); - - // Lower-level start following SW boundary from NW to SE. - if (!VISITED_CORNER(quad) && Z_NW <= 1 && Z_SE == 1) - contour.push_back(start_filled( - quad, Edge_SW, 1, NotHole, Boundary, lower_level, upper_level)); - - // Upper-level start following SW boundary from NW to SE. - if (!VISITED_CORNER(quad) && Z_NW == 2 && Z_SE == 1) - contour.push_back(start_filled( - quad, Edge_SW, 2, NotHole, Boundary, lower_level, upper_level)); - } - - // A full (unmasked) quad can only have a start on the NE corner, i.e. from - // N to E (lower level) or E to N (upper level). Any other start will have - // already been created in a call to this function for a prior quad so we - // don't need to test for it again here. - // - // The situation is complicated by the possibility that the quad is a - // saddle quad, in which case a contour line starting on the N could leave - // by either the W or the E. We only need to consider those leaving E. - // - // A NE corner can also have a N to E or E to N start. - if (EXISTS_QUAD(quad) || EXISTS_NE_CORNER(quad)) { - - // Lower-level start from N to E. - if (!VISITED(quad,1) && Z_NW == 0 && Z_SE == 0 && Z_NE >= 1 && - (!SADDLE(quad,1) || SADDLE_LEFT(quad,1))) - contour.push_back(start_filled( - quad, Edge_N, 1, NotHole, Interior, lower_level, upper_level)); - - // Upper-level start from E to N, hole. - if (!VISITED(quad,2) && Z_NW < 2 && Z_SE < 2 && Z_NE == 2 && - (!SADDLE(quad,2) || !SADDLE_LEFT(quad,2))) - contour.push_back(start_filled( - quad, Edge_E, 2, Hole, Interior, lower_level, upper_level)); - - // Upper-level start from N to E. - if (!VISITED(quad,2) && Z_NW == 2 && Z_SE == 2 && Z_NE < 2 && - (!SADDLE(quad,2) || SADDLE_LEFT(quad,2))) - contour.push_back(start_filled( - quad, Edge_N, 2, NotHole, Interior, lower_level, upper_level)); - - // Lower-level start from E to N, hole. - if (!VISITED(quad,1) && Z_NW >= 1 && Z_SE >= 1 && Z_NE == 0 && - (!SADDLE(quad,1) || !SADDLE_LEFT(quad,1))) - contour.push_back(start_filled( - quad, Edge_E, 1, Hole, Interior, lower_level, upper_level)); - - // All possible contours passing through the interior of this quad - // should have already been created, so assert this. - assert((VISITED(quad,1) || get_start_edge(quad, 1) == Edge_None) && - "Found start of contour that should have already been created"); - assert((VISITED(quad,2) || get_start_edge(quad, 2) == Edge_None) && - "Found start of contour that should have already been created"); - } - - // Lower-level start following N boundary from E to W, hole. - // This is required for an internal masked region which is a hole in a - // surrounding contour line. - if (BOUNDARY_N(quad) && EXISTS_N_EDGE(quad) && - !VISITED_S(quad+_nx) && Z_NW == 1 && Z_NE == 1) - contour.push_back(start_filled( - quad, Edge_N, 1, Hole, Boundary, lower_level, upper_level)); -} - -ContourLine* Mpl2014ContourGenerator::start_filled( - index_t quad, Edge edge, unsigned int start_level_index, HoleOrNot hole_or_not, - BoundaryOrInterior boundary_or_interior, const double& lower_level, const double& upper_level) -{ - assert(quad >= 0 && quad < _n && "Quad index out of bounds"); - assert(edge != Edge_None && "Invalid edge"); - assert((start_level_index == 1 || start_level_index == 2) && - "start level index must be 1 or 2"); - - ContourLine* contour_line = new ContourLine(hole_or_not == Hole); - if (hole_or_not == Hole) { - // Find and set parent ContourLine. - ContourLine* parent = _parent_cache.get_parent(quad + 1); - assert(parent != 0 && "Failed to find parent ContourLine"); - contour_line->set_parent(parent); - parent->add_child(contour_line); - } - - QuadEdge quad_edge(quad, edge); - const QuadEdge start_quad_edge(quad_edge); - unsigned int level_index = start_level_index; - - // If starts on interior, can only finish on interior. - // If starts on boundary, can only finish on boundary. - - while (true) { - if (boundary_or_interior == Interior) { - double level = (level_index == 1 ? lower_level : upper_level); - follow_interior( - *contour_line, quad_edge, level_index, level, false, &start_quad_edge, - start_level_index, true); - } - else { - level_index = follow_boundary( - *contour_line, quad_edge, lower_level, upper_level, level_index, start_quad_edge); - } - - if (quad_edge == start_quad_edge && (boundary_or_interior == Boundary || - level_index == start_level_index)) - break; - - if (boundary_or_interior == Boundary) - boundary_or_interior = Interior; - else - boundary_or_interior = Boundary; - } - - return contour_line; -} - -bool Mpl2014ContourGenerator::start_line( - py::list& vertices_list, py::list& codes_list, index_t quad, Edge edge, const double& level) -{ - assert(is_edge_a_boundary(QuadEdge(quad, edge)) && "QuadEdge is not a boundary"); - - QuadEdge quad_edge(quad, edge); - ContourLine contour_line(false); - follow_interior(contour_line, quad_edge, 1, level, true, 0, 1, false); - - append_contour_line_to_vertices_and_codes(contour_line, vertices_list, codes_list); - - return VISITED(quad,1); -} - -/*void Mpl2014ContourGenerator::write_cache(bool grid_only) const -{ - std::cout << "-----------------------------------------------" << std::endl; - for (index_t quad = 0; quad < _n; ++quad) - write_cache_quad(quad, grid_only); - std::cout << "-----------------------------------------------" << std::endl; -} - -void Mpl2014ContourGenerator::write_cache_quad( - index_t quad, bool grid_only) const -{ - index_t j = quad / _nx; - index_t i = quad - j*_nx; - std::cout << quad << ": i=" << i << " j=" << j << " EXISTS=" << EXISTS_QUAD(quad); - if (_corner_mask) - std::cout << " CORNER=" << EXISTS_SW_CORNER(quad) << EXISTS_SE_CORNER(quad) - << EXISTS_NW_CORNER(quad) << EXISTS_NE_CORNER(quad); - std::cout << " BNDY=" << BOUNDARY_S(quad) << BOUNDARY_W(quad); - if (!grid_only) { - std::cout << " Z=" << Z_LEVEL(quad) - << " SAD=" << SADDLE(quad,1) << SADDLE(quad,2) - << " LEFT=" << SADDLE_LEFT(quad,1) << SADDLE_LEFT(quad,2) - << " NW=" << SADDLE_START_SW(quad,1) << SADDLE_START_SW(quad,2) - << " VIS=" << VISITED(quad,1) << VISITED(quad,2) - << VISITED_S(quad) << VISITED_W(quad) - << VISITED_CORNER(quad); - } - std::cout << std::endl; -}*/ - -} // namespace mpl2014 -} // namespace contourpy diff --git a/contrib/python/contourpy/src/mpl2014.h b/contrib/python/contourpy/src/mpl2014.h deleted file mode 100644 index 5872e8ddcff..00000000000 --- a/contrib/python/contourpy/src/mpl2014.h +++ /dev/null @@ -1,513 +0,0 @@ -/* - * Mpl2014ContourGenerator - * ----------------------- - * A ContourGenerator generates contours for scalar fields defined on - * quadrilateral grids. A single QuadContourGenerator object can create both - * line contours (at single levels) and filled contours (between pairs of - * levels) for the same field. - * - * A field to be contoured has nx, ny points in the x- and y-directions - * respectively. The quad grid is defined by x and y arrays of shape(ny, nx), - * and the field itself is the z array also of shape(ny, nx). There is an - * optional boolean mask; if it exists then it also has shape(ny, nx). The - * mask applies to grid points rather than quads. - * - * How quads are masked based on the point mask is determined by the boolean - * 'corner_mask' flag. If false then any quad that has one or more of its four - * corner points masked is itself masked. If true the behaviour is the same - * except that any quad which has exactly one of its four corner points masked - * has only the triangular corner (half of the quad) adjacent to that point - * masked; the opposite triangular corner has three unmasked points and is not - * masked. - * - * By default the entire domain of nx*ny points is contoured together which can - * result in some very long polygons. The alternative is to break up the - * domain into subdomains or 'chunks' of smaller size, each of which is - * independently contoured. The size of these chunks is controlled by the - * 'nchunk' (or 'chunk_size') parameter. Chunking not only results in shorter - * polygons but also requires slightly less RAM. It can result in rendering - * artifacts though, depending on backend, antialiased flag and alpha value. - * - * Notation - * -------- - * i and j are array indices in the x- and y-directions respectively. Although - * a single element of an array z can be accessed using z[j][i] or z(j,i), it - * is often convenient to use the single quad index z[quad], where - * quad = i + j*nx - * and hence - * i = quad % nx - * j = quad / nx - * - * Rather than referring to x- and y-directions, compass directions are used - * instead such that W, E, S, N refer to the -x, +x, -y, +y directions - * respectively. To move one quad to the E you would therefore add 1 to the - * quad index, to move one quad to the N you would add nx to the quad index. - * - * Cache - * ----- - * Lots of information that is reused during contouring is stored as single - * bits in a mesh-sized cache, indexed by quad. Each quad's cache entry stores - * information about the quad itself such as if it is masked, and about the - * point at the SW corner of the quad, and about the W and S edges. Hence - * information about each point and each edge is only stored once in the cache. - * - * Cache information is divided into two types: that which is constant over the - * lifetime of the QuadContourGenerator, and that which changes for each - * contouring operation. The former is all grid-specific information such - * as quad and corner masks, and which edges are boundaries, either between - * masked and non-masked regions or between adjacent chunks. The latter - * includes whether points lie above or below the current contour levels, plus - * some flags to indicate how the contouring is progressing. - * - * Line Contours - * ------------- - * A line contour connects points with the same z-value. Each point of such a - * contour occurs on an edge of the grid, at a point linearly interpolated to - * the contour z-level from the z-values at the end points of the edge. The - * direction of a line contour is such that higher values are to the left of - * the contour, so any edge that the contour passes through will have a left- - * hand end point with z > contour level and a right-hand end point with - * z <= contour level. - * - * Line contours are of two types. Firstly there are open line strips that - * start on a boundary, traverse the interior of the domain and end on a - * boundary. Secondly there are closed line loops that occur completely within - * the interior of the domain and do not touch a boundary. - * - * The ContourGenerator makes two sweeps through the grid to generate line - * contours for a particular level. In the first sweep it looks only for start - * points that occur on boundaries, and when it finds one it follows the - * contour through the interior until it finishes on another boundary edge. - * Each quad that is visited by the algorithm has a 'visited' flag set in the - * cache to indicate that the quad does not need to be visited again. In the - * second sweep all non-visited quads are checked to see if they contain part - * of an interior closed loop, and again each time one is found it is followed - * through the domain interior until it returns back to its start quad and is - * therefore completed. - * - * The situation is complicated by saddle quads that have two opposite corners - * with z >= contour level and the other two corners with z < contour level. - * These therefore contain two segments of a line contour, and the visited - * flags take account of this by only being set on the second visit. On the - * first visit a number of saddle flags are set in the cache to indicate which - * one of the two segments has been completed so far. - * - * Filled Contours - * --------------- - * Filled contours are produced between two contour levels and are always - * closed polygons. They can occur completely within the interior of the - * domain without touching a boundary, following either the lower or upper - * contour levels. Those on the lower level are exactly like interior line - * contours with higher values on the left. Those on the upper level are - * reversed such that higher values are on the right. - * - * Filled contours can also involve a boundary in which case they consist of - * one or more sections along a boundary and one or more sections through the - * interior. Interior sections can be on either level, and again those on the - * upper level have higher values on the right. Boundary sections can remain - * on either contour level or switch between the two. - * - * Once the start of a filled contour is found, the algorithm is similar to - * that for line contours in that it follows the contour to its end, which - * because filled contours are always closed polygons will be by returning - * back to the start. However, because two levels must be considered, each - * level has its own set of saddle and visited flags and indeed some extra - * visited flags for boundary edges. - * - * The major complication for filled contours is that some polygons can be - * holes (with points ordered clockwise) within other polygons (with points - * ordered anticlockwise). When it comes to rendering filled contours each - * non-hole polygon must be rendered along with its zero or more contained - * holes or the rendering will not be correct. The filled contour finding - * algorithm could progress pretty much as the line contour algorithm does, - * taking each polygon as it is found, but then at the end there would have to - * be an extra step to identify the parent non-hole polygon for each hole. - * This is not a particularly onerous task but it does not scale well and can - * easily dominate the execution time of the contour finding for even modest - * problems. It is much better to identity each hole's parent non-hole during - * the sweep algorithm. - * - * This requirement dictates the order that filled contours are identified. As - * the algorithm sweeps up through the grid, every time a polygon passes - * through a quad a ParentCache object is updated with the new possible parent. - * When a new hole polygon is started, the ParentCache is used to find the - * first possible parent in the same quad or to the S of it. Great care is - * needed each time a new quad is checked to see if a new polygon should be - * started, as a single quad can have multiple polygon starts, e.g. a quad - * could be a saddle quad for both lower and upper contour levels, meaning it - * has four contour line segments passing through it which could all be from - * different polygons. The S-most polygon must be started first, then the next - * S-most and so on until the N-most polygon is started in that quad. - */ -#ifndef CONTOURPY_MPL_2014_H -#define CONTOURPY_MPL_2014_H - -#include "contour_generator.h" -#include <list> -#include <iostream> -#include <vector> - -namespace contourpy { -namespace mpl2014 { - -// Edge of a quad including diagonal edges of masked quads if _corner_mask true. -typedef enum -{ - // Listing values here so easier to check for debug purposes. - Edge_None = -1, - Edge_E = 0, - Edge_N = 1, - Edge_W = 2, - Edge_S = 3, - // The following are only used if _corner_mask is true. - Edge_NE = 4, - Edge_NW = 5, - Edge_SW = 6, - Edge_SE = 7 -} Edge; - -// Combination of a quad and an edge of that quad. -// An invalid quad edge has quad of -1. -struct QuadEdge -{ - QuadEdge() = delete; - QuadEdge(index_t quad_, Edge edge_); - bool operator==(const QuadEdge& other) const; - friend std::ostream& operator<<(std::ostream& os, const QuadEdge& quad_edge); - - index_t quad; - Edge edge; -}; - -// 2D point with x,y coordinates. -struct XY -{ - XY() = delete; - XY(const double& x_, const double& y_); - bool operator==(const XY& other) const; - friend std::ostream& operator<<(std::ostream& os, const XY& xy); - - double x, y; -}; - -// A single line of a contour, which may be a closed line loop or an open line -// strip. Identical adjacent points are avoided using push_back(). -// A ContourLine is either a hole (points ordered clockwise) or it is not -// (points ordered anticlockwise). Each hole has a parent ContourLine that is -// not a hole; each non-hole contains zero or more child holes. A non-hole and -// its child holes must be rendered together to obtain the correct results. -class ContourLine : public std::vector<XY> -{ -public: - typedef std::list<ContourLine*> Children; - - explicit ContourLine(bool is_hole); - void add_child(ContourLine* child); - void clear_parent(); - const Children& get_children() const; - const ContourLine* get_parent() const; - ContourLine* get_parent(); - bool is_hole() const; - void set_parent(ContourLine* parent); - void write() const; - -private: - bool _is_hole; - ContourLine* _parent; // Only set if is_hole, not owned. - Children _children; // Only set if !is_hole, not owned. -}; - - -// A Contour is a collection of zero or more ContourLines. -class Contour : public std::vector<ContourLine*> -{ -public: - Contour(); - virtual ~Contour(); - void delete_contour_lines(); - void write() const; -}; - - -// Single chunk of ContourLine parents, indexed by quad. As a chunk's filled -// contours are created, the ParentCache is updated each time a ContourLine -// passes through each quad. When a new ContourLine is created, if it is a -// hole its parent ContourLine is read from the ParentCache by looking at the -// start quad, then each quad to the S in turn until a non-zero ContourLine is -// found. -class ParentCache -{ -public: - ParentCache(index_t nx, index_t x_chunk_points, index_t y_chunk_points); - ContourLine* get_parent(index_t quad); - void set_chunk_starts(index_t istart, index_t jstart); - void set_parent(index_t quad, ContourLine& contour_line); - -private: - index_t index_to_index(index_t quad) const; - - index_t _nx; - index_t _x_chunk_points, _y_chunk_points; // Number of points not quads. - std::vector<ContourLine*> _lines; // Not owned. - index_t _istart, _jstart; -}; - - -// See overview of algorithm at top of file. -class Mpl2014ContourGenerator : public ContourGenerator -{ -public: - // Constructor with optional mask. - // x, y, z: double arrays of shape (ny,nx). - // mask: boolean array, ether empty (if no mask), or of shape (ny,nx). - // corner_mask: flag for different masking behaviour. - // x_chunk_size: 0 for no chunking, or +ve integer for size of chunks that - // the x-direction is subdivided into. - // y_chunk_size: 0 for no chunking, or +ve integer for size of chunks that - // the y-direction is subdivided into. - Mpl2014ContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, bool corner_mask, index_t x_chunk_size, index_t y_chunk_size); - - // Destructor. - virtual ~Mpl2014ContourGenerator(); - - // Create and return polygons for a filled contour between the two - // specified levels. - py::tuple filled(double lower_level, double upper_level) override; - - py::tuple get_chunk_count() const; // Return (y_chunk_count, x_chunk_count) - py::tuple get_chunk_size() const; // Return (y_chunk_size, x_chunk_size) - - bool get_corner_mask() const; - - // Create and return polygons for a line (i.e. non-filled) contour at the - // specified level. - py::sequence lines(double level) override; - -private: - // Typedef for following either a boundary of the domain or the interior; - // clearer than using a boolean. - typedef enum - { - Boundary, - Interior - } BoundaryOrInterior; - - // Typedef for direction of movement from one quad to the next. - typedef enum - { - Dir_Right = -1, - Dir_Straight = 0, - Dir_Left = +1 - } Dir; - - // Typedef for a polygon being a hole or not; clearer than using a boolean. - typedef enum - { - NotHole, - Hole - } HoleOrNot; - - // Append a C++ ContourLine to the end of two python lists. Used for line - // contours where each ContourLine is converted to a separate numpy array - // of (x,y) points and a second numpy array of 'kinds' or 'codes'. - // Clears the ContourLine too. - void append_contour_line_to_vertices_and_codes( - ContourLine& contour_line, py::list& vertices_list, py::list& codes_list) const; - - // Append a C++ Contour to the end of two python lists. Used for filled - // contours where each non-hole ContourLine and its child holes are - // represented by a numpy array of (x,y) points and a second numpy array of - // 'kinds' or 'codes' that indicates where the points array is split into - // individual polygons. - // Clears the Contour too, freeing each ContourLine as soon as possible - // for minimum RAM usage. - void append_contour_to_vertices_and_codes( - Contour& contour, py::list& vertices_list, py::list& codes_list) const; - - // Return number of chunks that fit in the specified point_count. - static index_t calc_chunk_count(index_t point_count, index_t chunk_size); - - // Return actual chunk_size from specified point_count and requested chunk_size. - static index_t calc_chunk_size(index_t point_count, index_t chunk_size); - - // Append the point on the specified QuadEdge that intersects the specified - // level to the specified ContourLine. - void edge_interp(const QuadEdge& quad_edge, const double& level, ContourLine& contour_line); - - // Follow a contour along a boundary, appending points to the ContourLine - // as it progresses. Only called for filled contours. Stops when the - // contour leaves the boundary to move into the interior of the domain, or - // when the start_quad_edge is reached in which case the ContourLine is a - // completed closed loop. Always adds the end point of each boundary edge - // to the ContourLine, regardless of whether moving to another boundary - // edge or leaving the boundary into the interior. Never adds the start - // point of the first boundary edge to the ContourLine. - // contour_line: ContourLine to append points to. - // quad_edge: on entry the QuadEdge to start from, on exit the QuadEdge - // that is stopped on. - // lower_level: lower contour z-value. - // upper_level: upper contour z-value. - // level_index: level index started on (1 = lower, 2 = upper level). - // start_quad_edge: QuadEdge that the ContourLine started from, which is - // used to check if the ContourLine is finished. - // Returns the end level_index. - unsigned int 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); - - // Follow a contour across the interior of the domain, appending points to - // the ContourLine as it progresses. Called for both line and filled - // contours. Stops when the contour reaches a boundary or, if the - // start_quad_edge is specified, when quad_edge == start_quad_edge and - // level_index == start_level_index. Always adds the end point of each - // quad traversed to the ContourLine; only adds the start point of the - // first quad if want_initial_point flag is true. - // contour_line: ContourLine to append points to. - // quad_edge: on entry the QuadEdge to start from, on exit the QuadEdge - // that is stopped on. - // level_index: level index started on (1 = lower, 2 = upper level). - // level: contour z-value. - // want_initial_point: whether want to append the initial point to the - // ContourLine or not. - // start_quad_edge: the QuadEdge that the ContourLine started from to - // check if the ContourLine is finished, or 0 if no check should occur. - // start_level_index: the level_index that the ContourLine started from. - // set_parents: whether should set ParentCache as it progresses or not. - // This is true for filled contours, false for line contours. - void 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); - - // Return the index limits of a particular chunk. - void get_chunk_limits( - index_t ijchunk, index_t& ichunk, index_t& jchunk, index_t& istart, index_t& iend, - index_t& jstart, index_t& jend); - - // Check if a contour starts within the specified corner quad on the - // specified level_index, and if so return the start edge. Otherwise - // return Edge_None. - Edge get_corner_start_edge(index_t quad, unsigned int level_index) const; - - // Return index of point at start or end of specified QuadEdge, assuming - // anticlockwise ordering around non-masked quads. - index_t get_edge_point_index(const QuadEdge& quad_edge, bool start) const; - - // Return the edge to exit a quad from, given the specified entry quad_edge - // and direction to move in. - Edge get_exit_edge(const QuadEdge& quad_edge, Dir dir) const; - - const double& get_point_x(index_t point) const; - const double& get_point_y(index_t point) const; - - // Append the (x,y) coordinates of the specified point index to the - // specified ContourLine. - void get_point_xy(index_t point, ContourLine& contour_line) const; - - // Return the z-value of the specified point index. - const double& get_point_z(index_t point) const; - - // Check if a contour starts within the specified non-corner quad on the - // specified level_index, and if so return the start edge. Otherwise - // return Edge_None. - Edge get_quad_start_edge(index_t quad, unsigned int level_index) const; - - // Check if a contour starts within the specified quad, whether it is a - // corner or a full quad, and if so return the start edge. Otherwise - // return Edge_None. - Edge get_start_edge(index_t quad, unsigned int level_index) const; - - // Initialise the cache to contain grid information that is constant - // across the lifetime of this object, i.e. does not vary between calls to - // create_contour() and create_filled_contour(). - void init_cache_grid(const MaskArray& mask); - - // Initialise the cache with information that is specific to contouring the - // specified two levels. The levels are the same for contour lines, - // different for filled contours. - void init_cache_levels(const double& lower_level, const double& upper_level); - - // Append the (x,y) point at which the level intersects the line connecting - // the two specified point indices to the specified ContourLine. - void interp( - index_t point1, index_t point2, const double& level, ContourLine& contour_line) const; - - // Return true if the specified QuadEdge is a boundary, i.e. is either an - // edge between a masked and non-masked quad/corner or is a chunk boundary. - bool is_edge_a_boundary(const QuadEdge& quad_edge) const; - - // Follow a boundary from one QuadEdge to the next in an anticlockwise - // manner around the non-masked region. - void move_to_next_boundary_edge(QuadEdge& quad_edge) const; - - // Move from the quad specified by quad_edge.quad to the neighbouring quad - // by crossing the edge specified by quad_edge.edge. - void move_to_next_quad(QuadEdge& quad_edge) const; - - // Check for filled contours starting within the specified quad and - // complete any that are found, appending them to the specified Contour. - void single_quad_filled( - Contour& contour, index_t quad, const double& lower_level, const double& upper_level); - - // Start and complete a filled contour line. - // quad: index of quad to start ContourLine in. - // edge: edge of quad to start ContourLine from. - // start_level_index: the level_index that the ContourLine starts from. - // hole_or_not: whether the ContourLine is a hole or not. - // boundary_or_interior: whether the ContourLine starts on a boundary or - // the interior. - // lower_level: lower contour z-value. - // upper_level: upper contour z-value. - // Returns newly created ContourLine. - ContourLine* start_filled( - index_t quad, Edge edge, unsigned int start_level_index, HoleOrNot hole_or_not, - BoundaryOrInterior boundary_or_interior, const double& lower_level, - const double& upper_level); - - // Start and complete a line contour that both starts and end on a - // boundary, traversing the interior of the domain. - // vertices_list: Python list that the ContourLine should be appended to. - // codes_list: Python list that the kind codes should be appended to. - // quad: index of quad to start ContourLine in. - // edge: boundary edge to start ContourLine from. - // level: contour z-value. - // Returns true if the start quad does not need to be visited again, i.e. - // VISITED(quad,1). - bool start_line( - py::list& vertices_list, py::list& codes_list, index_t quad, Edge edge, - const double& level); - - // Debug function that writes the cache status to stdout. - //void write_cache(bool grid_only = false) const; - - // Debug function that writes that cache status for a single quad to - // stdout. - //void write_cache_quad(index_t quad, bool grid_only) const; - - - - // Note that mask is not stored as once it has been used to initialise the - // cache it is no longer needed. - CoordinateArray _x, _y, _z; - index_t _nx, _ny; // Number of points in each direction. - index_t _n; // Total number of points (and hence quads). - - bool _corner_mask; - index_t _x_chunk_size; // Number of quads per chunk (not points). - index_t _y_chunk_size; // Always > 0. - - index_t _nxchunk, _nychunk; // Number of chunks in each direction. - index_t _chunk_count; // Total number of chunks. - - typedef uint32_t CacheItem; - CacheItem* _cache; - - ParentCache _parent_cache; // On W quad sides. -}; - -} // namespace mpl2014 -} // namespace contourpy - -#endif // #define CONTOURPY_MPL_2014_H diff --git a/contrib/python/contourpy/src/mpl_kind_code.h b/contrib/python/contourpy/src/mpl_kind_code.h deleted file mode 100644 index 11209f8be86..00000000000 --- a/contrib/python/contourpy/src/mpl_kind_code.h +++ /dev/null @@ -1,17 +0,0 @@ -// Enum for kind codes used in Matplotlib Paths. - -#ifndef MPL_KIND_CODE -#define MPL_KIND_CODE - -namespace contourpy { - -typedef enum -{ - MOVETO = 1, - LINETO = 2, - CLOSEPOLY = 79 -} MplKindCode; - -} // namespace contourpy - -#endif // MPL_KIND_CODE diff --git a/contrib/python/contourpy/src/outer_or_hole.cpp b/contrib/python/contourpy/src/outer_or_hole.cpp deleted file mode 100644 index 47a91aeb7c1..00000000000 --- a/contrib/python/contourpy/src/outer_or_hole.cpp +++ /dev/null @@ -1,19 +0,0 @@ -#include "outer_or_hole.h" -#include <iostream> - -namespace contourpy { - -std::ostream &operator<<(std::ostream &os, const OuterOrHole& outer_or_hole) -{ - switch (outer_or_hole) { - case Outer: - os << "Outer"; - break; - case Hole: - os << "Hole"; - break; - } - return os; -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/outer_or_hole.h b/contrib/python/contourpy/src/outer_or_hole.h deleted file mode 100644 index 0089a0bf753..00000000000 --- a/contrib/python/contourpy/src/outer_or_hole.h +++ /dev/null @@ -1,18 +0,0 @@ -#ifndef CONTOURPY_OUTER_OR_HOLE_H -#define CONTOURPY_OUTER_OR_HOLE_H - -#include <iosfwd> - -namespace contourpy { - -typedef enum -{ - Outer, - Hole -} OuterOrHole; - -std::ostream &operator<<(std::ostream &os, const OuterOrHole& outer_or_hole); - -} // namespace contourpy - -#endif // CONTOURPY_OUTER_OR_HOLE_H diff --git a/contrib/python/contourpy/src/output_array.h b/contrib/python/contourpy/src/output_array.h deleted file mode 100644 index ce6be00d37c..00000000000 --- a/contrib/python/contourpy/src/output_array.h +++ /dev/null @@ -1,70 +0,0 @@ -#ifndef CONTOURPY_OUTPUT_ARRAY_H -#define CONTOURPY_OUTPUT_ARRAY_H - -#include "common.h" -#include <vector> - -namespace contourpy { - -// A reusable array that is output from C++ to Python. Depending on the chosen line or fill type, -// it can either be created as a NumPy array that will be directly returned to the Python caller, -// or as a C++ vector that will be further manipulated (such as split up) before being converted to -// NumPy array(s) for returning. BaseContourGenerator's marching does not care which form it is as -// it just writes values to either array using an incrementing pointer. -template <typename T> -class OutputArray -{ -public: - OutputArray() - : size(0), start(nullptr), current(nullptr) - {} - - void clear() - { - vector.clear(); - size = 0; - start = current = nullptr; - } - - void create_cpp(count_t new_size) - { - assert(new_size > 0); - size = new_size; - vector.resize(size); - start = current = vector.data(); - } - - py::array_t<T> create_python(count_t new_size) - { - assert(new_size > 0); - size = new_size; - py::array_t<T> py_array(size); - start = current = py_array.mutable_data(); - return py_array; - } - - py::array_t<T> create_python(count_t shape0, count_t shape1) - { - assert(shape0 > 0 && shape1 > 0); - size = shape0*shape1; - py::array_t<T> py_array({shape0, shape1}); - start = current = py_array.mutable_data(); - return py_array; - } - - // Non-copyable and non-moveable. - OutputArray(const OutputArray& other) = delete; - OutputArray(const OutputArray&& other) = delete; - OutputArray& operator=(const OutputArray& other) = delete; - OutputArray& operator=(const OutputArray&& other) = delete; - - - std::vector<T> vector; - count_t size; - T* start; // Start of array, whether C++ or Python. - T* current; // Where to write next value to before incrementing. -}; - -} // namespace contourpy - -#endif // CONTOURPY_OUTPUT_ARRAY_H diff --git a/contrib/python/contourpy/src/serial.cpp b/contrib/python/contourpy/src/serial.cpp deleted file mode 100644 index dd69471ef64..00000000000 --- a/contrib/python/contourpy/src/serial.cpp +++ /dev/null @@ -1,147 +0,0 @@ -#include "base_impl.h" -#include "converter.h" -#include "serial.h" - -namespace contourpy { - -SerialContourGenerator::SerialContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, bool corner_mask, LineType line_type, FillType fill_type, - bool quad_as_tri, ZInterp z_interp, index_t x_chunk_size, index_t y_chunk_size) - : BaseContourGenerator(x, y, z, mask, corner_mask, line_type, fill_type, quad_as_tri, z_interp, - x_chunk_size, y_chunk_size) -{} - -void SerialContourGenerator::export_filled( - const ChunkLocal& local, std::vector<py::list>& return_lists) -{ - assert(local.total_point_count > 0); - - switch (get_fill_type()) - { - case FillType::OuterCode: - case FillType::OuterOffset: { - assert(!has_direct_points() && !has_direct_line_offsets()); - auto outer_count = local.line_count - local.hole_count; - - for (decltype(outer_count) i = 0; i < outer_count; ++i) { - auto outer_start = local.outer_offsets.start[i]; - auto outer_end = local.outer_offsets.start[i+1]; - auto point_start = local.line_offsets.start[outer_start]; - auto point_end = local.line_offsets.start[outer_end]; - auto point_count = point_end - point_start; - assert(point_count > 2); - - return_lists[0].append(Converter::convert_points( - point_count, local.points.start + 2*point_start)); - - if (get_fill_type() == FillType::OuterCode) - return_lists[1].append(Converter::convert_codes( - point_count, outer_end - outer_start + 1, - local.line_offsets.start + outer_start, point_start)); - else - return_lists[1].append(Converter::convert_offsets( - outer_end - outer_start + 1, local.line_offsets.start + outer_start, - point_start)); - } - break; - } - case FillType::ChunkCombinedCode: - case FillType::ChunkCombinedCodeOffset: { - assert(has_direct_points() && !has_direct_line_offsets()); - - // return_lists[0][local_chunk] already contains combined points. - // If ChunkCombinedCodeOffset. return_lists[2][local.chunk] already contains outer - // offsets. - return_lists[1][local.chunk] = Converter::convert_codes( - local.total_point_count, local.line_count + 1, local.line_offsets.start, 0); - break; - } - case FillType::ChunkCombinedOffset: - case FillType::ChunkCombinedOffsetOffset: - assert(has_direct_points() && has_direct_line_offsets()); - if (get_fill_type() == FillType::ChunkCombinedOffsetOffset) { - assert(has_direct_outer_offsets()); - } - // return_lists[0][local_chunk] already contains combined points. - // return_lists[1][local.chunk] already contains line offsets. - // If ChunkCombinedOffsetOffset, return_lists[2][local.chunk] already contains - // outer offsets. - break; - } -} - -void SerialContourGenerator::export_lines( - const ChunkLocal& local, std::vector<py::list>& return_lists) -{ - assert(local.total_point_count > 0); - - switch (get_line_type()) - { - case LineType::Separate: - case LineType::SeparateCode: { - assert(!has_direct_points() && !has_direct_line_offsets()); - - bool separate_code = (get_line_type() == LineType::SeparateCode); - - for (decltype(local.line_count) i = 0; i < local.line_count; ++i) { - auto point_start = local.line_offsets.start[i]; - auto point_end = local.line_offsets.start[i+1]; - auto point_count = point_end - point_start; - assert(point_count > 1); - - return_lists[0].append(Converter::convert_points( - point_count, local.points.start + 2*point_start)); - - if (separate_code) { - return_lists[1].append( - Converter::convert_codes_check_closed_single( - point_count, local.points.start + 2*point_start)); - } - } - break; - } - case LineType::ChunkCombinedCode: { - assert(has_direct_points() && !has_direct_line_offsets()); - - // return_lists[0][local.chunk] already contains points. - return_lists[1][local.chunk] = Converter::convert_codes_check_closed( - local.total_point_count, local.line_count + 1, local.line_offsets.start, - local.points.start); - break; - } - case LineType::ChunkCombinedOffset: - assert(has_direct_points() && has_direct_line_offsets()); - // return_lists[0][local.chunk] already contains points. - // return_lists[1][local.chunk] already contains line offsets. - break; - case LineType::ChunkCombinedNan: - assert(has_direct_points()); - // return_lists[0][local.chunk] already contains points. - break; - } -} - -void SerialContourGenerator::march(std::vector<py::list>& return_lists) -{ - auto n_chunks = get_n_chunks(); - bool single_chunk = (n_chunks == 1); - - if (single_chunk) { - // Stage 1: If single chunk, initialise cache z-levels and starting locations for whole - // domain. - init_cache_levels_and_starts(); - } - - // Stage 2: Trace contours. - ChunkLocal local; - for (index_t chunk = 0; chunk < n_chunks; ++chunk) { - get_chunk_limits(chunk, local); - if (!single_chunk) - init_cache_levels_and_starts(&local); - march_chunk(local, return_lists); - local.clear(); - } -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/serial.h b/contrib/python/contourpy/src/serial.h deleted file mode 100644 index 0ee288f22e6..00000000000 --- a/contrib/python/contourpy/src/serial.h +++ /dev/null @@ -1,40 +0,0 @@ -#ifndef CONTOURPY_SERIAL_H -#define CONTOURPY_SERIAL_H - -#include "base.h" - -namespace contourpy { - -class SerialContourGenerator : public BaseContourGenerator<SerialContourGenerator> -{ -public: - SerialContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, bool corner_mask, LineType line_type, FillType fill_type, - bool quad_as_tri, ZInterp z_interp, index_t x_chunk_size, index_t y_chunk_size); - -private: - friend class BaseContourGenerator<SerialContourGenerator>; - - // Dummy Lock class which is the single-threaded version of ThreadedContourGenerator::Lock and - // does not do anything, allowing base class code to use Lock objects for both serial and - // multithreaded code. - class Lock - { - public: - explicit Lock(SerialContourGenerator& contour_generator) - {} - }; - - // Write points and offsets/codes to output numpy arrays. - void export_filled(const ChunkLocal& local, std::vector<py::list>& return_lists); - - // Write points and offsets/codes to output numpy arrays. - void export_lines(const ChunkLocal& local, std::vector<py::list>& return_lists); - - void march(std::vector<py::list>& return_lists); -}; - -} // namespace contourpy - -#endif // CONTOURPY_SERIAL_H diff --git a/contrib/python/contourpy/src/threaded.cpp b/contrib/python/contourpy/src/threaded.cpp deleted file mode 100644 index 9f6eec1beda..00000000000 --- a/contrib/python/contourpy/src/threaded.cpp +++ /dev/null @@ -1,316 +0,0 @@ -#include "base_impl.h" -#include "converter.h" -#include "threaded.h" -#include "util.h" -#include <thread> - -namespace contourpy { - -ThreadedContourGenerator::ThreadedContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, bool corner_mask, LineType line_type, FillType fill_type, - bool quad_as_tri, ZInterp z_interp, index_t x_chunk_size, index_t y_chunk_size, - index_t n_threads) - : BaseContourGenerator(x, y, z, mask, corner_mask, line_type, fill_type, quad_as_tri, z_interp, - x_chunk_size, y_chunk_size), - _n_threads(limit_n_threads(n_threads, get_n_chunks())), - _next_chunk(0) -{} - -void ThreadedContourGenerator::export_filled( - const ChunkLocal& local, std::vector<py::list>& return_lists) -{ - // Reimplementation of SerialContourGenerator::export_filled() to separate out the creation of - // numpy arrays (which requires a thread lock) from the population of those arrays (which does - // not). This minimises the time that the lock is used for. - - assert(local.total_point_count > 0); - - switch (get_fill_type()) - { - case FillType::OuterCode: - case FillType::OuterOffset: { - assert(!has_direct_points() && !has_direct_line_offsets()); - - auto outer_count = local.line_count - local.hole_count; - bool outer_code = (get_fill_type() == FillType::OuterCode); - std::vector<PointArray::value_type*> points_ptrs(outer_count); - std::vector<CodeArray::value_type*> codes_ptrs(outer_code ? outer_count: 0); - std::vector<OffsetArray::value_type*> offsets_ptrs(outer_code ? 0 : outer_count); - - { - Lock lock(*this); // cppcheck-suppress unreadVariable - for (decltype(outer_count) i = 0; i < outer_count; ++i) { - auto outer_start = local.outer_offsets.start[i]; - auto outer_end = local.outer_offsets.start[i+1]; - auto point_start = local.line_offsets.start[outer_start]; - auto point_end = local.line_offsets.start[outer_end]; - auto point_count = point_end - point_start; - assert(point_count > 2); - - index_t points_shape[2] = {static_cast<index_t>(point_count), 2}; - PointArray point_array(points_shape); - return_lists[0].append(point_array); - points_ptrs[i] = point_array.mutable_data(); - - if (outer_code) { - index_t codes_shape = static_cast<index_t>(point_count); - CodeArray code_array(codes_shape); - return_lists[1].append(code_array); - codes_ptrs[i] = code_array.mutable_data(); - } - else { - index_t offsets_shape = static_cast<index_t>(outer_end - outer_start + 1); - OffsetArray offset_array(offsets_shape); - return_lists[1].append(offset_array); - offsets_ptrs[i] = offset_array.mutable_data(); - } - } - } - - for (decltype(outer_count) i = 0; i < outer_count; ++i) { - auto outer_start = local.outer_offsets.start[i]; - auto outer_end = local.outer_offsets.start[i+1]; - auto point_start = local.line_offsets.start[outer_start]; - auto point_end = local.line_offsets.start[outer_end]; - auto point_count = point_end - point_start; - assert(point_count > 2); - - Converter::convert_points( - point_count, local.points.start + 2*point_start, points_ptrs[i]); - - if (outer_code) - Converter::convert_codes( - point_count, outer_end - outer_start + 1, - local.line_offsets.start + outer_start, point_start, codes_ptrs[i]); - else - Converter::convert_offsets( - outer_end - outer_start + 1, local.line_offsets.start + outer_start, - point_start, offsets_ptrs[i]); - } - break; - } - case FillType::ChunkCombinedCode: - case FillType::ChunkCombinedCodeOffset: { - assert(has_direct_points() && !has_direct_line_offsets()); - // return_lists[0][local_chunk] already contains combined points. - // If ChunkCombinedCodeOffset. return_lists[2][local.chunk] already contains outer - // offsets. - - index_t codes_shape = static_cast<index_t>(local.total_point_count); - CodeArray::value_type* codes_ptr = nullptr; - - { - Lock lock(*this); // cppcheck-suppress unreadVariable - CodeArray code_array(codes_shape); - return_lists[1][local.chunk] = code_array; - codes_ptr = code_array.mutable_data(); - } - - Converter::convert_codes( - local.total_point_count, local.line_count + 1, local.line_offsets.start, 0, - codes_ptr); - break; - } - case FillType::ChunkCombinedOffset: - case FillType::ChunkCombinedOffsetOffset: - assert(has_direct_points() && has_direct_line_offsets()); - if (get_fill_type() == FillType::ChunkCombinedOffsetOffset) { - assert(has_direct_outer_offsets()); - } - // return_lists[0][local_chunk] already contains combined points. - // return_lists[1][local.chunk] already contains line offsets. - // If ChunkCombinedOffsetOffset, return_lists[2][local.chunk] already contains - // outer offsets. - break; - } -} - -void ThreadedContourGenerator::export_lines( - const ChunkLocal& local, std::vector<py::list>& return_lists) -{ - // Reimplementation of SerialContourGenerator::export_lines() to separate out the creation of - // numpy arrays (which requires a thread lock) from the population of those arrays (which does - // not). This minimises the time that the lock is used for. - - assert(local.total_point_count > 0); - - switch (get_line_type()) - { - case LineType::Separate: - case LineType::SeparateCode: { - assert(!has_direct_points() && !has_direct_line_offsets()); - - bool separate_code = (get_line_type() == LineType::SeparateCode); - std::vector<PointArray::value_type*> points_ptrs(local.line_count); - std::vector<CodeArray::value_type*> codes_ptrs(separate_code ? local.line_count: 0); - - { - Lock lock(*this); // cppcheck-suppress unreadVariable - for (decltype(local.line_count) i = 0; i < local.line_count; ++i) { - auto point_start = local.line_offsets.start[i]; - auto point_end = local.line_offsets.start[i+1]; - auto point_count = point_end - point_start; - assert(point_count > 1); - - index_t points_shape[2] = {static_cast<index_t>(point_count), 2}; - PointArray point_array(points_shape); - - return_lists[0].append(point_array); - points_ptrs[i] = point_array.mutable_data(); - - if (separate_code) { - index_t codes_shape = static_cast<index_t>(point_count); - CodeArray code_array(codes_shape); - return_lists[1].append(code_array); - codes_ptrs[i] = code_array.mutable_data(); - } - } - } - - for (decltype(local.line_count) i = 0; i < local.line_count; ++i) { - auto point_start = local.line_offsets.start[i]; - auto point_end = local.line_offsets.start[i+1]; - auto point_count = point_end - point_start; - assert(point_count > 1); - - Converter::convert_points( - point_count, local.points.start + 2*point_start, points_ptrs[i]); - - if (separate_code) { - Converter::convert_codes_check_closed_single( - point_count, local.points.start + 2*point_start, codes_ptrs[i]); - } - } - break; - } - case LineType::ChunkCombinedCode: { - assert(has_direct_points() && !has_direct_line_offsets()); - // return_lists[0][local.chunk] already contains points. - - index_t codes_shape = static_cast<index_t>(local.total_point_count); - CodeArray::value_type* codes_ptr = nullptr; - - { - Lock lock(*this); // cppcheck-suppress unreadVariable - CodeArray code_array(codes_shape); - return_lists[1][local.chunk] = code_array; - codes_ptr = code_array.mutable_data(); - } - - Converter::convert_codes_check_closed( - local.total_point_count, local.line_count + 1, local.line_offsets.start, - local.points.start, codes_ptr); - break; - } - case LineType::ChunkCombinedOffset: - assert(has_direct_points() && has_direct_line_offsets()); - // return_lists[0][local.chunk] already contains points. - // return_lists[1][local.chunk] already contains line offsets. - break; - case LineType::ChunkCombinedNan: - assert(has_direct_points()); - // return_lists[0][local.chunk] already contains points. - break; - } -} - -index_t ThreadedContourGenerator::get_thread_count() const -{ - return _n_threads; -} - -index_t ThreadedContourGenerator::limit_n_threads(index_t n_threads, index_t n_chunks) -{ - index_t max_threads = std::max<index_t>(Util::get_max_threads(), 1); - if (n_threads == 0) - return std::min(max_threads, n_chunks); - else - return std::min({max_threads, n_chunks, n_threads}); -} - -void ThreadedContourGenerator::march(std::vector<py::list>& return_lists) -{ - // Each thread executes thread_function() which has two stages: - // 1) Initialise cache z-levels and starting locations - // 2) Trace contours - // Each stage is performed on a chunk by chunk basis. There is a barrier between the two stages - // to synchronise the threads so the cache setup is complete before being used by the trace. - _next_chunk = 0; // Next available chunk index. - _finished_count = 0; // Count of threads that have finished the cache init. - - // Main thread releases GIL for remainder of this function. - // It is temporarily reacquired as necessary within the scope of threaded Lock objects. - py::gil_scoped_release release; - - // Create (_n_threads-1) new worker threads. - std::vector<std::thread> threads; - threads.reserve(_n_threads-1); - for (index_t i = 0; i < _n_threads-1; ++i) - threads.emplace_back( - &ThreadedContourGenerator::thread_function, this, std::ref(return_lists)); - - thread_function(std::ref(return_lists)); // Main thread work. - - for (auto& thread : threads) - thread.join(); - assert(_next_chunk == 2*get_n_chunks()); - threads.clear(); -} - -void ThreadedContourGenerator::thread_function(std::vector<py::list>& return_lists) -{ - // Function that is executed by each of the threads. - // _next_chunk starts at zero and increases up to 2*_n_chunks. A thread in need of work reads - // _next_chunk and incremements it, then processes that chunk. For _next_chunk < _n_chunks this - // is stage 1 (init cache levels and starting locations) and for _next_chunk >= _n_chunks this - // is stage 2 (trace contours). There is a synchronisation barrier between the two stages so - // that the cache initialisation is complete before being used by the contour trace. - - auto n_chunks = get_n_chunks(); - index_t chunk; - ChunkLocal local; - - // Stage 1: Initialise cache z-levels and starting locations. - while (true) { - { - std::lock_guard<std::mutex> guard(_chunk_mutex); - if (_next_chunk < n_chunks) - chunk = _next_chunk++; - else - break; // No more work to do. - } - - get_chunk_limits(chunk, local); - init_cache_levels_and_starts(&local); - local.clear(); - } - - { - // Implementation of multithreaded barrier. Each thread increments the shared counter. - // Last thread to finish notifies the other threads that they can all continue. - std::unique_lock<std::mutex> lock(_chunk_mutex); - _finished_count++; - if (_finished_count == _n_threads) - _condition_variable.notify_all(); - else - _condition_variable.wait(lock); - } - - // Stage 2: Trace contours. - while (true) { - { - std::lock_guard<std::mutex> guard(_chunk_mutex); - if (_next_chunk < 2*n_chunks) - chunk = _next_chunk++ - n_chunks; - else - break; // No more work to do. - } - - get_chunk_limits(chunk, local); - march_chunk(local, return_lists); - local.clear(); - } -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/threaded.h b/contrib/python/contourpy/src/threaded.h deleted file mode 100644 index ff45c1c5446..00000000000 --- a/contrib/python/contourpy/src/threaded.h +++ /dev/null @@ -1,69 +0,0 @@ -#ifndef CONTOURPY_THREADED_H -#define CONTOURPY_THREADED_H - -#include "base.h" -#include <condition_variable> -#include <mutex> - -namespace contourpy { - -class ThreadedContourGenerator : public BaseContourGenerator<ThreadedContourGenerator> -{ -public: - ThreadedContourGenerator( - const CoordinateArray& x, const CoordinateArray& y, const CoordinateArray& z, - const MaskArray& mask, bool corner_mask, LineType line_type, FillType fill_type, - bool quad_as_tri, ZInterp z_interp, index_t x_chunk_size, index_t y_chunk_size, - index_t n_threads); - - index_t get_thread_count() const; - -private: - friend class BaseContourGenerator<ThreadedContourGenerator>; - - // Lock class is used to lock access to a single thread when creating/modifying Python objects. - // Also acquires the GIL for the duration of the lock. - class Lock - { - public: - explicit Lock(ThreadedContourGenerator& contour_generator) - : _lock(contour_generator._python_mutex) - {} - - // Non-copyable and non-moveable. - Lock(const Lock& other) = delete; - Lock(const Lock&& other) = delete; - Lock& operator=(const Lock& other) = delete; - Lock& operator=(const Lock&& other) = delete; - - private: - std::unique_lock<std::mutex> _lock; - py::gil_scoped_acquire _gil; - }; - - // Write points and offsets/codes to output numpy arrays. - void export_filled(const ChunkLocal& local, std::vector<py::list>& return_lists); - - // Write points and offsets/codes to output numpy arrays. - void export_lines(const ChunkLocal& local, std::vector<py::list>& return_lists); - - static index_t limit_n_threads(index_t n_threads, index_t n_chunks); - - void march(std::vector<py::list>& return_lists); - - void thread_function(std::vector<py::list>& return_lists); - - - - // Multithreading member variables. - index_t _n_threads; // Number of threads used. - index_t _next_chunk; // Next available chunk for thread to process. - index_t _finished_count; // Count of threads that have finished the cache init. - std::mutex _chunk_mutex; // Locks access to _next_chunk/_finished_count. - std::mutex _python_mutex; // Locks access to Python objects. - std::condition_variable _condition_variable; // Implements multithreaded barrier. -}; - -} // namespace contourpy - -#endif // CONTOURPY_THREADED_H diff --git a/contrib/python/contourpy/src/util.cpp b/contrib/python/contourpy/src/util.cpp deleted file mode 100644 index 82c1cc21d93..00000000000 --- a/contrib/python/contourpy/src/util.cpp +++ /dev/null @@ -1,30 +0,0 @@ -#include "util.h" -#include <cmath> -#include <thread> - -namespace contourpy { - -bool Util::_nan_loaded = false; - -double Util::nan = 0.0; - -void Util::ensure_nan_loaded() -{ - if (!_nan_loaded) { - auto numpy = py::module_::import("numpy"); - nan = numpy.attr("nan").cast<double>(); - _nan_loaded = true; - } -} - -index_t Util::get_max_threads() -{ - return static_cast<index_t>(std::thread::hardware_concurrency()); -} - -bool Util::is_nan(double value) -{ - return std::isnan(value); -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/util.h b/contrib/python/contourpy/src/util.h deleted file mode 100644 index dcafd131c15..00000000000 --- a/contrib/python/contourpy/src/util.h +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef CONTOURPY_UTIL_H -#define CONTOURPY_UTIL_H - -#include "common.h" - -namespace contourpy { - -class Util -{ -public: - static void ensure_nan_loaded(); - - static index_t get_max_threads(); - - static bool is_nan(double value); - - // This is the NaN used internally and returned to calling functions. The value is taken from - // numpy rather than the standard C++ approach so that it is guaranteed to work with - // numpy.isnan(). The value is actually the same for many platforms, but this approach - // guarantees it works for all platforms that numpy supports. - // - // ensure_nan_loaded() must be called before this value is read. - static double nan; - -private: - static bool _nan_loaded; -}; - -} // namespace contourpy - -#endif // CONTOURPY_UTIL_H diff --git a/contrib/python/contourpy/src/wrap.cpp b/contrib/python/contourpy/src/wrap.cpp deleted file mode 100644 index fede2280956..00000000000 --- a/contrib/python/contourpy/src/wrap.cpp +++ /dev/null @@ -1,419 +0,0 @@ -#include "base_impl.h" -#include "contour_generator.h" -#include "fill_type.h" -#include "line_type.h" -#include "mpl2005.h" -#include "mpl2014.h" -#include "serial.h" -#include "threaded.h" -#include "util.h" -#include "z_interp.h" - -#define STRINGIFY(x) #x -#define MACRO_STRINGIFY(x) STRINGIFY(x) - -namespace py = pybind11; -using namespace pybind11::literals; - -static contourpy::LineType mpl20xx_line_type = contourpy::LineType::SeparateCode; -static contourpy::FillType mpl20xx_fill_type = contourpy::FillType::OuterCode; - -PYBIND11_MODULE(_contourpy, m, py::mod_gil_not_used()) { - m.doc() = - "C++11 extension module wrapped using `pybind11`_.\n\n" - "It should not be necessary to access classes and functions in this extension module " - "directly. Instead, :func:`~contourpy.contour_generator` should be used to create " - ":class:`~.ContourGenerator` objects, and the enums " - "(:class:`~.FillType`, :class:`~.LineType` and " - ":class:`~.ZInterp`) and :func:`.max_threads` function are all available " - "in the :mod:`contourpy` module."; - - m.attr("__version__") = MACRO_STRINGIFY(CONTOURPY_VERSION); - - // asserts are enabled if NDEBUG is 0. - m.attr("NDEBUG") = -#ifdef NDEBUG - 1; -#else - 0; -#endif - - py::enum_<contourpy::FillType>(m, "FillType", - "Enum used for ``fill_type`` keyword argument in :func:`~.contour_generator`.\n\n" - "This controls the format of filled contour data returned from " - ":meth:`~.ContourGenerator.filled`.") - .value("OuterCode", contourpy::FillType::OuterCode) - .value("OuterOffset", contourpy::FillType::OuterOffset) - .value("ChunkCombinedCode", contourpy::FillType::ChunkCombinedCode) - .value("ChunkCombinedOffset", contourpy::FillType::ChunkCombinedOffset) - .value("ChunkCombinedCodeOffset", contourpy::FillType::ChunkCombinedCodeOffset) - .value("ChunkCombinedOffsetOffset", contourpy::FillType::ChunkCombinedOffsetOffset) - .export_values(); - - py::enum_<contourpy::LineType>(m, "LineType", - "Enum used for ``line_type`` keyword argument in :func:`~.contour_generator`.\n\n" - "This controls the format of contour line data returned from " - ":meth:`~.ContourGenerator.lines`.\n\n" - "``LineType.ChunkCombinedNan`` added in version 1.2.0") - .value("Separate", contourpy::LineType::Separate) - .value("SeparateCode", contourpy::LineType::SeparateCode) - .value("ChunkCombinedCode", contourpy::LineType::ChunkCombinedCode) - .value("ChunkCombinedOffset", contourpy::LineType::ChunkCombinedOffset) - .value("ChunkCombinedNan", contourpy::LineType::ChunkCombinedNan) - .export_values(); - - py::enum_<contourpy::ZInterp>(m, "ZInterp", - "Enum used for ``z_interp`` keyword argument in :func:`~.contour_generator`\n\n" - "This controls the interpolation used on ``z`` values to determine where contour lines " - "intersect the edges of grid quads, and ``z`` values at quad centres.") - .value("Linear", contourpy::ZInterp::Linear) - .value("Log", contourpy::ZInterp::Log) - .export_values(); - - m.def("max_threads", &contourpy::Util::get_max_threads, - "Return the maximum number of threads, obtained from " - "``std::thread::hardware_concurrency()``.\n\n" - "This is the number of threads used by a multithreaded ContourGenerator if the kwarg " - "``threads=0`` is passed to :func:`~.contour_generator`."); - - const char* chunk_count_doc = "Return tuple of (y, x) chunk counts."; - const char* chunk_size_doc = "Return tuple of (y, x) chunk sizes."; - const char* corner_mask_doc = "Return whether ``corner_mask`` is set or not."; - const char* create_contour_doc = - "Synonym for :meth:`~.ContourGenerator.lines` to provide backward compatibility " - "with Matplotlib."; - const char* create_filled_contour_doc = - "Synonym for :meth:`~.ContourGenerator.filled` to provide backward compatibility " - "with Matplotlib."; - const char* default_fill_type_doc = - "Return the default :class:`~.FillType` used by this algorithm."; - const char* default_line_type_doc = - "Return the default :class:`~.LineType` used by this algorithm."; - const char* fill_type_doc = "Return the :class:`~.FillType`."; - const char* filled_doc = - "Calculate and return filled contours between two levels.\n\n" - "Args:\n" - " lower_level (float): Lower z-level of the filled contours, cannot be ``np.nan``.\n" - " upper_level (float): Upper z-level of the filled contours, cannot be ``np.nan``.\n" - "Return:\n" - " Filled contour polygons as nested sequences of numpy arrays. The exact format is " - "determined by the :attr:`~.ContourGenerator.fill_type` used by the ``ContourGenerator`` " - "and the options are explained at :ref:`fill_type`.\n\n" - "Raises a ``ValueError`` if ``lower_level >= upper_level`` or if\n" - "``lower_level`` or ``upper_level`` are ``np.nan``.\n\n" - "To return filled contours below a ``level`` use ``filled(-np.inf, level)``.\n" - "To return filled contours above a ``level`` use ``filled(level, np.inf)``"; - const char* line_type_doc = "Return the :class:`~.LineType`."; - const char* lines_doc = - "Calculate and return contour lines at a particular level.\n\n" - "Args:\n" - " level (float): z-level to calculate contours at.\n\n" - "Return:\n" - " Contour lines (open line strips and closed line loops) as nested sequences of " - "numpy arrays. The exact format is determined by the :attr:`~.ContourGenerator.line_type` " - "used by the ``ContourGenerator`` and the options are explained at :ref:`line_type`.\n\n" - "``level`` may be ``np.nan``, ``np.inf`` or ``-np.inf``; they all return the same result " - "which is an empty line set."; - const char* multi_filled_doc = - "Calculate and return filled contours between multiple pairs of adjacent levels.\n\n" - "Args:\n" - " levels (array-like of floats): z-levels to calculate filled contours between. " - "There must be at least 2 levels, they cannot be NaN, and each level must be larger than " - "the previous level.\n\n" - "Return:\n" - " List of filled contours, one per pair of levels. The length of the returned list is " - "one less than ``len(levels)``. The format of returned contours is determined by the " - "``fill_type`` used by the ``ContourGenerator`` and the options are explained at " - ":ref:`fill_type`.\n\n" - "Raises a ``ValueError`` if ``levels`` are not increasing, or contain a NaN, or there are " - "fewer than 2 ``levels``.\n\n" - "Calling:\n\n" - ".. code-block:: python\n\n" - " ret = cont_gen.multi_filled(levels)\n\n" - "is equivalent to:\n\n" - ".. code-block:: python\n\n" - " ret = [cont_gen.filled(lower, upper) for lower, upper in zip(levels[:-1], levels[1:])]\n\n" - ".. versionadded:: 1.3.0"; - const char* multi_lines_doc = - "Calculate and return contour lines at multiple levels.\n\n" - "Args:\n" - " levels (array-like of floats): z-levels to calculate contours at.\n\n" - "Return:\n" - " List of contour lines, one per level. The format of returned lines is determined by " - "the ``line_type`` used by the ``ContourGenerator`` and the options are explained at " - ":ref:`line_type`.\n\n" - "Calling:\n\n" - ".. code-block:: python\n\n" - " ret = cont_gen.multi_lines(levels)\n\n" - "is equivalent to:\n\n" - ".. code-block:: python\n\n" - " ret = [cont_gen.lines(level) for level in levels]\n\n" - ".. versionadded:: 1.3.0"; - const char* quad_as_tri_doc = "Return whether ``quad_as_tri`` is set or not."; - const char* supports_corner_mask_doc = - "Return whether this algorithm supports ``corner_mask``."; - const char* supports_fill_type_doc = - "Return whether this algorithm supports a particular :class:`~.FillType`."; - const char* supports_line_type_doc = - "Return whether this algorithm supports a particular :class:`~.LineType`."; - const char* supports_quad_as_tri_doc = - "Return whether this algorithm supports ``quad_as_tri``."; - const char* supports_threads_doc = - "Return whether this algorithm supports the use of threads."; - const char* supports_z_interp_doc = - "Return whether this algorithm supports ``z_interp`` values other than ``ZInterp.Linear`` " - "which all support."; - const char* thread_count_doc = "Return the number of threads used."; - const char* z_interp_doc = "Return the ``ZInterp``."; - - py::class_<contourpy::ContourGenerator>(m, "ContourGenerator", - "Abstract base class for contour generator classes, defining the interface that they all " - "implement.") - .def("create_contour", &contourpy::ContourGenerator::lines, create_contour_doc, "level"_a) - .def("create_filled_contour", &contourpy::ContourGenerator::filled, - create_filled_contour_doc, "lower_level"_a, "upper_level"_a) - .def("filled", &contourpy::ContourGenerator::filled, filled_doc, - "lower_level"_a, "upper_level"_a) - .def("lines", &contourpy::ContourGenerator::lines, lines_doc, "level"_a) - .def("multi_filled", &contourpy::ContourGenerator::multi_filled, multi_filled_doc, - "levels"_a) - .def("multi_lines", &contourpy::ContourGenerator::multi_lines, multi_lines_doc, "levels"_a) - .def_property_readonly( - "chunk_count", [](py::object /* self */) {return py::make_tuple(1, 1);}, - chunk_count_doc) - .def_property_readonly( - "chunk_size", [](py::object /* self */) {return py::make_tuple(1, 1);}, chunk_size_doc) - .def_property_readonly( - "corner_mask", [](py::object /* self */) {return false;}, corner_mask_doc) - .def_property_readonly( - "fill_type", [](py::object /* self */) {return contourpy::FillType::OuterOffset;}, - fill_type_doc) - .def_property_readonly( - "line_type", [](py::object /* self */) {return contourpy::LineType::Separate;}, - line_type_doc) - .def_property_readonly( - "quad_as_tri", [](py::object /* self */) {return false;}, quad_as_tri_doc) - .def_property_readonly( - "thread_count", [](py::object /* self */) {return 1;}, thread_count_doc) - .def_property_readonly( - "z_interp", [](py::object /* self */) {return contourpy::ZInterp::Linear;}, - z_interp_doc) - .def_property_readonly_static( - "default_fill_type", - [](py::object /* self */) {return contourpy::FillType::OuterOffset;}, - default_fill_type_doc) - .def_property_readonly_static( - "default_line_type", [](py::object /* self */) {return contourpy::LineType::Separate;}, - default_line_type_doc) - .def_static("supports_corner_mask", []() {return false;}, supports_corner_mask_doc) - .def_static( - "supports_fill_type", [](contourpy::FillType /* fill_type */) {return false;}, - "fill_type"_a, supports_fill_type_doc) - .def_static( - "supports_line_type", [](contourpy::LineType /* line_type */) {return false;}, - "line_type"_a, supports_line_type_doc) - .def_static("supports_quad_as_tri", []() {return false;}, supports_quad_as_tri_doc) - .def_static("supports_threads", []() {return false;}, supports_threads_doc) - .def_static("supports_z_interp", []() {return false;}, supports_z_interp_doc); - - py::class_<contourpy::Mpl2005ContourGenerator, contourpy::ContourGenerator>( - m, "Mpl2005ContourGenerator", - "ContourGenerator corresponding to ``name=\"mpl2005\"``.\n\n" - "This is the original 2005 Matplotlib algorithm. " - "Does not support any of ``corner_mask``, ``quad_as_tri``, ``threads`` or ``z_interp``. " - "Only supports ``line_type=LineType.SeparateCode`` and ``fill_type=FillType.OuterCode``. " - "Only supports chunking for filled contours, not contour lines.\n\n" - ".. warning::\n" - " This algorithm is in ``contourpy`` for historic comparison. No new features or bug " - "fixes will be added to it, except for security-related bug fixes.") - .def(py::init<const contourpy::CoordinateArray&, - const contourpy::CoordinateArray&, - const contourpy::CoordinateArray&, - const contourpy::MaskArray&, - contourpy::index_t, - contourpy::index_t>(), - "x"_a, "y"_a, "z"_a, "mask"_a, py::kw_only(), "x_chunk_size"_a = 0, - "y_chunk_size"_a = 0) - .def_property_readonly( - "chunk_count", &contourpy::Mpl2005ContourGenerator::get_chunk_count, chunk_count_doc) - .def_property_readonly( - "chunk_size", &contourpy::Mpl2005ContourGenerator::get_chunk_size, chunk_size_doc) - .def_property_readonly( - "fill_type", [](py::object /* self */) {return mpl20xx_fill_type;}, fill_type_doc) - .def_property_readonly( - "line_type", [](py::object /* self */) {return mpl20xx_line_type;}, line_type_doc) - .def_property_readonly_static( - "default_fill_type", [](py::object /* self */) {return mpl20xx_fill_type;}, - default_fill_type_doc) - .def_property_readonly_static( - "default_line_type", [](py::object /* self */) {return mpl20xx_line_type;}, - default_line_type_doc) - .def_static( - "supports_fill_type", - [](contourpy::FillType fill_type) {return fill_type == mpl20xx_fill_type;}, - supports_fill_type_doc) - .def_static( - "supports_line_type", - [](contourpy::LineType line_type) {return line_type == mpl20xx_line_type;}, - supports_line_type_doc); - - py::class_<contourpy::mpl2014::Mpl2014ContourGenerator, contourpy::ContourGenerator>( - m, "Mpl2014ContourGenerator", - "ContourGenerator corresponding to ``name=\"mpl2014\"``.\n\n" - "This is the 2014 Matplotlib algorithm, a replacement of the original 2005 algorithm that " - "added ``corner_mask`` and made the code more maintainable. " - "Only supports ``corner_mask``, does not support ``quad_as_tri``, ``threads`` or " - "``z_interp``. \n" - "Only supports ``line_type=LineType.SeparateCode`` and " - "``fill_type=FillType.OuterCode``.\n\n" - ".. warning::\n" - " This algorithm is in ``contourpy`` for historic comparison. No new features or bug " - "fixes will be added to it, except for security-related bug fixes.") - .def(py::init<const contourpy::CoordinateArray&, - const contourpy::CoordinateArray&, - const contourpy::CoordinateArray&, - const contourpy::MaskArray&, - bool, - contourpy::index_t, - contourpy::index_t>(), - "x"_a, "y"_a, "z"_a, "mask"_a, py::kw_only(), "corner_mask"_a, "x_chunk_size"_a = 0, - "y_chunk_size"_a = 0) - .def_property_readonly( - "chunk_count", &contourpy::mpl2014::Mpl2014ContourGenerator::get_chunk_count, - chunk_count_doc) - .def_property_readonly( - "chunk_size", &contourpy::mpl2014::Mpl2014ContourGenerator::get_chunk_size, - chunk_size_doc) - .def_property_readonly( - "corner_mask", &contourpy::mpl2014::Mpl2014ContourGenerator::get_corner_mask, - corner_mask_doc) - .def_property_readonly( - "fill_type", [](py::object /* self */) {return mpl20xx_fill_type;}, fill_type_doc) - .def_property_readonly( - "line_type", [](py::object /* self */) {return mpl20xx_line_type;}, line_type_doc) - .def_property_readonly_static( - "default_fill_type", [](py::object /* self */) {return mpl20xx_fill_type;}, - default_fill_type_doc) - .def_property_readonly_static( - "default_line_type", [](py::object /* self */) {return mpl20xx_line_type;}, - default_line_type_doc) - .def_static("supports_corner_mask", []() {return true;}, supports_corner_mask_doc) - .def_static( - "supports_fill_type", - [](contourpy::FillType fill_type) {return fill_type == mpl20xx_fill_type;}, - supports_fill_type_doc) - .def_static( - "supports_line_type", - [](contourpy::LineType line_type) {return line_type == mpl20xx_line_type;}, - supports_line_type_doc); - - py::class_<contourpy::SerialContourGenerator, contourpy::ContourGenerator>( - m, "SerialContourGenerator", - "ContourGenerator corresponding to ``name=\"serial\"``, the default algorithm for " - "``contourpy``.\n\n" - "Supports ``corner_mask``, ``quad_as_tri`` and ``z_interp`` but not ``threads``. " - "Supports all options for ``line_type`` and ``fill_type``.") - .def(py::init<const contourpy::CoordinateArray&, - const contourpy::CoordinateArray&, - const contourpy::CoordinateArray&, - const contourpy::MaskArray&, - bool, - contourpy::LineType, - contourpy::FillType, - bool, - contourpy::ZInterp, - contourpy::index_t, - contourpy::index_t>(), - "x"_a, "y"_a, "z"_a, "mask"_a, py::kw_only(), "corner_mask"_a, "line_type"_a, - "fill_type"_a, "quad_as_tri"_a, "z_interp"_a, "x_chunk_size"_a = 0, - "y_chunk_size"_a = 0) - .def("_write_cache", &contourpy::SerialContourGenerator::write_cache) - .def_property_readonly( - "chunk_count", &contourpy::SerialContourGenerator::get_chunk_count, chunk_count_doc) - .def_property_readonly( - "chunk_size", &contourpy::SerialContourGenerator::get_chunk_size, chunk_size_doc) - .def_property_readonly( - "corner_mask", &contourpy::SerialContourGenerator::get_corner_mask, corner_mask_doc) - .def_property_readonly( - "fill_type", &contourpy::SerialContourGenerator::get_fill_type, fill_type_doc) - .def_property_readonly( - "line_type", &contourpy::SerialContourGenerator::get_line_type, line_type_doc) - .def_property_readonly( - "quad_as_tri", &contourpy::SerialContourGenerator::get_quad_as_tri, quad_as_tri_doc) - .def_property_readonly( - "z_interp", &contourpy::SerialContourGenerator::get_z_interp, z_interp_doc) - .def_property_readonly_static( - "default_fill_type", - [](py::object /* self */) {return contourpy::SerialContourGenerator::default_fill_type();}, - default_fill_type_doc) - .def_property_readonly_static("default_line_type", - [](py::object /* self */) {return contourpy::SerialContourGenerator::default_line_type();}, - default_line_type_doc) - .def_static("supports_corner_mask", []() {return true;}, supports_corner_mask_doc) - .def_static( - "supports_fill_type", &contourpy::SerialContourGenerator::supports_fill_type, - supports_fill_type_doc) - .def_static( - "supports_line_type", &contourpy::SerialContourGenerator::supports_line_type, - supports_line_type_doc) - .def_static("supports_quad_as_tri", []() {return true;}, supports_quad_as_tri_doc) - .def_static("supports_z_interp", []() {return true;}, supports_z_interp_doc); - - py::class_<contourpy::ThreadedContourGenerator, contourpy::ContourGenerator>( - m, "ThreadedContourGenerator", - "ContourGenerator corresponding to ``name=\"threaded\"``, the multithreaded version of " - ":class:`~.SerialContourGenerator`.\n\n" - "Supports ``corner_mask``, ``quad_as_tri`` and ``z_interp`` and ``threads``. " - "Supports all options for ``line_type`` and ``fill_type``.") - .def(py::init<const contourpy::CoordinateArray&, - const contourpy::CoordinateArray&, - const contourpy::CoordinateArray&, - const contourpy::MaskArray&, - bool, - contourpy::LineType, - contourpy::FillType, - bool, - contourpy::ZInterp, - contourpy::index_t, - contourpy::index_t, - contourpy::index_t>(), - "x"_a, "y"_a, "z"_a, "mask"_a, py::kw_only(), "corner_mask"_a, "line_type"_a, - "fill_type"_a, "quad_as_tri"_a, "z_interp"_a, "x_chunk_size"_a = 0, - "y_chunk_size"_a = 0, "thread_count"_a = 0) - .def("_write_cache", &contourpy::ThreadedContourGenerator::write_cache) - .def_property_readonly( - "chunk_count", &contourpy::ThreadedContourGenerator::get_chunk_count, chunk_count_doc) - .def_property_readonly( - "chunk_size", &contourpy::ThreadedContourGenerator::get_chunk_size, chunk_size_doc) - .def_property_readonly( - "corner_mask", &contourpy::ThreadedContourGenerator::get_corner_mask, corner_mask_doc) - .def_property_readonly( - "fill_type", &contourpy::ThreadedContourGenerator::get_fill_type, fill_type_doc) - .def_property_readonly( - "line_type", &contourpy::ThreadedContourGenerator::get_line_type, line_type_doc) - .def_property_readonly( - "quad_as_tri", &contourpy::ThreadedContourGenerator::get_quad_as_tri, quad_as_tri_doc) - .def_property_readonly( - "thread_count", &contourpy::ThreadedContourGenerator::get_thread_count, - thread_count_doc) - .def_property_readonly( - "z_interp", &contourpy::ThreadedContourGenerator::get_z_interp, z_interp_doc) - .def_property_readonly_static( - "default_fill_type", - [](py::object /* self */) {return contourpy::ThreadedContourGenerator::default_fill_type();}, - default_fill_type_doc) - .def_property_readonly_static( - "default_line_type", - [](py::object /* self */) {return contourpy::ThreadedContourGenerator::default_line_type();}, - default_line_type_doc) - .def_static("supports_corner_mask", []() {return true;}, supports_corner_mask_doc) - .def_static( - "supports_fill_type", &contourpy::ThreadedContourGenerator::supports_fill_type, - supports_fill_type_doc) - .def_static( - "supports_line_type", &contourpy::ThreadedContourGenerator::supports_line_type, - supports_line_type_doc) - .def_static("supports_quad_as_tri", []() {return true;}, supports_quad_as_tri_doc) - .def_static("supports_threads", []() {return true;}, supports_threads_doc) - .def_static("supports_z_interp", []() {return true;}, supports_z_interp_doc); -} diff --git a/contrib/python/contourpy/src/z_interp.cpp b/contrib/python/contourpy/src/z_interp.cpp deleted file mode 100644 index 92811ad0d02..00000000000 --- a/contrib/python/contourpy/src/z_interp.cpp +++ /dev/null @@ -1,19 +0,0 @@ -#include "z_interp.h" -#include <iostream> - -namespace contourpy { - -std::ostream &operator<<(std::ostream &os, const ZInterp& z_interp) -{ - switch (z_interp) { - case ZInterp::Linear: - os << "Linear"; - break; - case ZInterp::Log: - os << "Log"; - break; - } - return os; -} - -} // namespace contourpy diff --git a/contrib/python/contourpy/src/z_interp.h b/contrib/python/contourpy/src/z_interp.h deleted file mode 100644 index 76a7550be62..00000000000 --- a/contrib/python/contourpy/src/z_interp.h +++ /dev/null @@ -1,23 +0,0 @@ -#ifndef CONTOURPY_Z_INTERP_H -#define CONTOURPY_Z_INTERP_H - -#include <iosfwd> -#include <string> - -namespace contourpy { - -// Enum for type of interpolation used to find intersection of contour lines -// with grid cell edges. - -// C++11 scoped enum, must be fully qualified to use. -enum class ZInterp -{ - Linear = 1, - Log = 2 -}; - -std::ostream &operator<<(std::ostream &os, const ZInterp& z_interp); - -} // namespace contourpy - -#endif // CONTOURPY_ZINTERP_H |