diff options
| author | shumkovnd <shumkovnd@yandex-team.com> | 2023-11-10 14:39:34 +0300 |
|---|---|---|
| committer | shumkovnd <shumkovnd@yandex-team.com> | 2023-11-10 16:42:24 +0300 |
| commit | 77eb2d3fdcec5c978c64e025ced2764c57c00285 (patch) | |
| tree | c51edb0748ca8d4a08d7c7323312c27ba1a8b79a /contrib/python/contourpy/src | |
| parent | dd6d20cadb65582270ac23f4b3b14ae189704b9d (diff) | |
| download | ydb-77eb2d3fdcec5c978c64e025ced2764c57c00285.tar.gz | |
KIKIMR-19287: add task_stats_drawing script
Diffstat (limited to 'contrib/python/contourpy/src')
31 files changed, 8238 insertions, 0 deletions
diff --git a/contrib/python/contourpy/src/base.h b/contrib/python/contourpy/src/base.h new file mode 100644 index 0000000000..7faf9335ec --- /dev/null +++ b/contrib/python/contourpy/src/base.h @@ -0,0 +1,210 @@ +// 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: + ~BaseContourGenerator(); + + static FillType default_fill_type(); + static LineType default_line_type(); + + 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 filled(double lower_level, double upper_level); + py::sequence lines(double level); + + 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_filled() 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; + + 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. + 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 new file mode 100644 index 0000000000..61b9c46132 --- /dev/null +++ b/contrib/python/contourpy/src/base_impl.h @@ -0,0 +1,2454 @@ +#ifndef CONTOURPY_BASE_IMPL_H +#define CONTOURPY_BASE_IMPL_H + +#include "base.h" +#include "converter.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), + _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::sequence BaseContourGenerator<Derived>::filled(double lower_level, double upper_level) +{ + if (lower_level > upper_level) + throw std::invalid_argument("upper and lower levels are the wrong way round"); + + _filled = true; + _lower_level = lower_level; + _upper_level = upper_level; + + _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); + _return_list_count = (_fill_type == FillType::ChunkCombinedCodeOffset || + _fill_type == FillType::ChunkCombinedOffsetOffset) ? 3 : 2; + + 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_filled() const +{ + return _filled; +} + +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; + + // 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) +{ + _filled = false; + _lower_level = _upper_level = level; + + _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) ? 1 : 2; + + 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 (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 && _line_type == LineType::Separate); + return 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> +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: + 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 new file mode 100644 index 0000000000..1c8833ab0b --- /dev/null +++ b/contrib/python/contourpy/src/chunk_local.cpp @@ -0,0 +1,59 @@ +#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 new file mode 100644 index 0000000000..f9b68ea6be --- /dev/null +++ b/contrib/python/contourpy/src/chunk_local.h @@ -0,0 +1,40 @@ +#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; + 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 new file mode 100644 index 0000000000..cc4bc851a0 --- /dev/null +++ b/contrib/python/contourpy/src/common.h @@ -0,0 +1,32 @@ +#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; + +// 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.h b/contrib/python/contourpy/src/contour_generator.h new file mode 100644 index 0000000000..ce527b25f8 --- /dev/null +++ b/contrib/python/contourpy/src/contour_generator.h @@ -0,0 +1,23 @@ +#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; + +protected: + ContourGenerator() = default; +}; + +} // namespace contourpy + +#endif // CONTOURPY_CONTOUR_GENERATOR_H diff --git a/contrib/python/contourpy/src/converter.cpp b/contrib/python/contourpy/src/converter.cpp new file mode 100644 index 0000000000..1f693ceb1a --- /dev/null +++ b/contrib/python/contourpy/src/converter.cpp @@ -0,0 +1,155 @@ +#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 && subtract >= 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 && subtract >= 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 && subtract >= 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 && subtract >= 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 new file mode 100644 index 0000000000..5bf9ab5119 --- /dev/null +++ b/contrib/python/contourpy/src/converter.h @@ -0,0 +1,63 @@ +#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 new file mode 100644 index 0000000000..79129e2909 --- /dev/null +++ b/contrib/python/contourpy/src/fill_type.cpp @@ -0,0 +1,31 @@ +#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 new file mode 100644 index 0000000000..f77c765d28 --- /dev/null +++ b/contrib/python/contourpy/src/fill_type.h @@ -0,0 +1,24 @@ +#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 new file mode 100644 index 0000000000..5ef2d68dac --- /dev/null +++ b/contrib/python/contourpy/src/line_type.cpp @@ -0,0 +1,25 @@ +#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; + } + return os; +} + +} // namespace contourpy diff --git a/contrib/python/contourpy/src/line_type.h b/contrib/python/contourpy/src/line_type.h new file mode 100644 index 0000000000..afc6edf0c3 --- /dev/null +++ b/contrib/python/contourpy/src/line_type.h @@ -0,0 +1,22 @@ +#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, +}; + +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 new file mode 100644 index 0000000000..dfa5c54619 --- /dev/null +++ b/contrib/python/contourpy/src/mpl2005.cpp @@ -0,0 +1,76 @@ +#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(const double& lower_level, const double& upper_level) +{ + if (lower_level > upper_level) + throw std::invalid_argument("upper and lower levels are the wrong way round"); + + 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::tuple Mpl2005ContourGenerator::lines(const 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 new file mode 100644 index 0000000000..e112880765 --- /dev/null +++ b/contrib/python/contourpy/src/mpl2005.h @@ -0,0 +1,32 @@ +#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); + + ~Mpl2005ContourGenerator(); + + py::tuple filled(const double& lower_level, const double& upper_level); + + 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::tuple lines(const double& level); + +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 new file mode 100644 index 0000000000..f2f0e8567c --- /dev/null +++ b/contrib/python/contourpy/src/mpl2005_original.cpp @@ -0,0 +1,1526 @@ +/* + 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 new file mode 100644 index 0000000000..93dba45ff9 --- /dev/null +++ b/contrib/python/contourpy/src/mpl2005_original.h @@ -0,0 +1,56 @@ +#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 new file mode 100644 index 0000000000..99107ad9ea --- /dev/null +++ b/contrib/python/contourpy/src/mpl2014.cpp @@ -0,0 +1,1713 @@ +// 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) + 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( + const double& lower_level, const double& upper_level) +{ + if (lower_level > upper_level) + throw std::invalid_argument("upper and lower levels are the wrong way round"); + + 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); + + if (first_edge) + 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::tuple Mpl2014ContourGenerator::lines(const 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 new file mode 100644 index 0000000000..105ce23c1f --- /dev/null +++ b/contrib/python/contourpy/src/mpl2014.h @@ -0,0 +1,513 @@ +/* + * 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. + ~Mpl2014ContourGenerator(); + + // Create and return polygons for a filled contour between the two + // specified levels. + py::tuple filled(const double& lower_level, const double& upper_level); + + 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::tuple lines(const double& level); + +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 new file mode 100644 index 0000000000..11209f8be8 --- /dev/null +++ b/contrib/python/contourpy/src/mpl_kind_code.h @@ -0,0 +1,17 @@ +// 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 new file mode 100644 index 0000000000..47a91aeb7c --- /dev/null +++ b/contrib/python/contourpy/src/outer_or_hole.cpp @@ -0,0 +1,19 @@ +#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 new file mode 100644 index 0000000000..0089a0bf75 --- /dev/null +++ b/contrib/python/contourpy/src/outer_or_hole.h @@ -0,0 +1,18 @@ +#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 new file mode 100644 index 0000000000..ce6be00d37 --- /dev/null +++ b/contrib/python/contourpy/src/output_array.h @@ -0,0 +1,70 @@ +#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 new file mode 100644 index 0000000000..08d851ca34 --- /dev/null +++ b/contrib/python/contourpy/src/serial.cpp @@ -0,0 +1,143 @@ +#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; + } +} + +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 new file mode 100644 index 0000000000..0ee288f22e --- /dev/null +++ b/contrib/python/contourpy/src/serial.h @@ -0,0 +1,40 @@ +#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 new file mode 100644 index 0000000000..e27cd55d51 --- /dev/null +++ b/contrib/python/contourpy/src/threaded.cpp @@ -0,0 +1,312 @@ +#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; + } +} + +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 new file mode 100644 index 0000000000..ff45c1c544 --- /dev/null +++ b/contrib/python/contourpy/src/threaded.h @@ -0,0 +1,69 @@ +#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 new file mode 100644 index 0000000000..d41f0edef9 --- /dev/null +++ b/contrib/python/contourpy/src/util.cpp @@ -0,0 +1,11 @@ +#include "util.h" +#include <thread> + +namespace contourpy { + +index_t Util::get_max_threads() +{ + return static_cast<index_t>(std::thread::hardware_concurrency()); +} + +} // namespace contourpy diff --git a/contrib/python/contourpy/src/util.h b/contrib/python/contourpy/src/util.h new file mode 100644 index 0000000000..22dcaca7fd --- /dev/null +++ b/contrib/python/contourpy/src/util.h @@ -0,0 +1,16 @@ +#ifndef CONTOURPY_UTIL_H +#define CONTOURPY_UTIL_H + +#include "common.h" + +namespace contourpy { + +class Util +{ +public: + static index_t get_max_threads(); +}; + +} // namespace contourpy + +#endif // CONTOURPY_UTIL_H diff --git a/contrib/python/contourpy/src/wrap.cpp b/contrib/python/contourpy/src/wrap.cpp new file mode 100644 index 0000000000..2c79641990 --- /dev/null +++ b/contrib/python/contourpy/src/wrap.cpp @@ -0,0 +1,427 @@ +#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; + +static contourpy::LineType mpl20xx_line_type = contourpy::LineType::SeparateCode; +static contourpy::FillType mpl20xx_fill_type = contourpy::FillType::OuterCode; + +PYBIND11_MODULE(_contourpy, m) { + m.doc() = + "C++11 extension module wrapped using `pybind11`_.\n\n" + ".. note::\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:`~contourpy.ContourGenerator` objects, and the enums " + "(:class:`~contourpy.FillType`, :class:`~contourpy.LineType` and " + ":class:`~contourpy.ZInterp`) and :func:`contourpy.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:`~contourpy.contour_generator`.\n\n" + "This controls the format of filled contour data returned from " + ":meth:`~contourpy.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:`~contourpy.contour_generator`.\n\n" + "This controls the format of contour line data returned from " + ":meth:`~contourpy.ContourGenerator.lines`.") + .value("Separate", contourpy::LineType::Separate) + .value("SeparateCode", contourpy::LineType::SeparateCode) + .value("ChunkCombinedCode", contourpy::LineType::ChunkCombinedCode) + .value("ChunkCombinedOffset", contourpy::LineType::ChunkCombinedOffset) + .export_values(); + + py::enum_<contourpy::ZInterp>(m, "ZInterp", + "Enum used for ``z_interp`` keyword argument in :func:`~contourpy.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:`~contourpy.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 :func:`~contourpy.ContourGenerator.lines` to provide backward compatibility " + "with Matplotlib."; + const char* create_filled_contour_doc = + "Synonym for :func:`~contourpy.ContourGenerator.filled` to provide backward compatibility " + "with Matplotlib."; + const char* default_fill_type_doc = "Return the default ``FillType`` used by this algorithm."; + const char* default_line_type_doc = "Return the default ``LineType`` used by this algorithm."; + const char* fill_type_doc = "Return the ``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.\n" + " upper_level (float): Upper z-level of the filled contours.\n" + "Return:\n" + " Filled contour polygons as one or more sequences of numpy arrays. The exact format is " + "determined by the ``fill_type`` used by the ``ContourGenerator``."; + const char* line_type_doc = "Return the ``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 one or more sequences of " + "numpy arrays. The exact format is determined by the ``line_type`` used by the " + "``ContourGenerator``."; + 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 ``FillType``."; + const char* supports_line_type_doc = + "Return whether this algorithm supports a particular ``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", + [](py::object /* self */, double level) {return py::make_tuple();}, + py::arg("level"), create_contour_doc) + .def("create_filled_contour", + [](py::object /* self */, double lower_level, double upper_level) {return py::make_tuple();}, + py::arg("lower_level"), py::arg("upper_level"), create_filled_contour_doc) + .def("filled", + [](py::object /* self */, double lower_level, double upper_level) {return py::make_tuple();}, + py::arg("lower_level"), py::arg("upper_level"), filled_doc) + .def("lines", + [](py::object /* self */, double level) {return py::make_tuple();}, + py::arg("level"), lines_doc) + .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;}, + py::arg("fill_type"), supports_fill_type_doc) + .def_static( + "supports_line_type", [](contourpy::LineType /* line_type */) {return false;}, + py::arg("line_type"), 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>(), + py::arg("x"), + py::arg("y"), + py::arg("z"), + py::arg("mask"), + py::kw_only(), + py::arg("x_chunk_size") = 0, + py::arg("y_chunk_size") = 0) + .def("create_contour", &contourpy::Mpl2005ContourGenerator::lines, create_contour_doc) + .def("create_filled_contour", &contourpy::Mpl2005ContourGenerator::filled, + create_filled_contour_doc) + .def("filled", &contourpy::Mpl2005ContourGenerator::filled, filled_doc) + .def("lines", &contourpy::Mpl2005ContourGenerator::lines, lines_doc) + .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>(), + py::arg("x"), + py::arg("y"), + py::arg("z"), + py::arg("mask"), + py::kw_only(), + py::arg("corner_mask"), + py::arg("x_chunk_size") = 0, + py::arg("y_chunk_size") = 0) + .def("create_contour", &contourpy::mpl2014::Mpl2014ContourGenerator::lines, + create_contour_doc) + .def("create_filled_contour", &contourpy::mpl2014::Mpl2014ContourGenerator::filled, + create_filled_contour_doc) + .def("filled", &contourpy::mpl2014::Mpl2014ContourGenerator::filled, filled_doc) + .def("lines", &contourpy::mpl2014::Mpl2014ContourGenerator::lines, lines_doc) + .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>(), + py::arg("x"), + py::arg("y"), + py::arg("z"), + py::arg("mask"), + py::kw_only(), + py::arg("corner_mask"), + py::arg("line_type"), + py::arg("fill_type"), + py::arg("quad_as_tri"), + py::arg("z_interp"), + py::arg("x_chunk_size") = 0, + py::arg("y_chunk_size") = 0) + .def("_write_cache", &contourpy::SerialContourGenerator::write_cache) + .def("create_contour", &contourpy::SerialContourGenerator::lines, create_contour_doc) + .def("create_filled_contour", &contourpy::SerialContourGenerator::filled, + create_filled_contour_doc) + .def("filled", &contourpy::SerialContourGenerator::filled, filled_doc) + .def("lines", &contourpy::SerialContourGenerator::lines, lines_doc) + .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:`~contourpy._contourpy.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>(), + py::arg("x"), + py::arg("y"), + py::arg("z"), + py::arg("mask"), + py::kw_only(), + py::arg("corner_mask"), + py::arg("line_type"), + py::arg("fill_type"), + py::arg("quad_as_tri"), + py::arg("z_interp"), + py::arg("x_chunk_size") = 0, + py::arg("y_chunk_size") = 0, + py::arg("thread_count") = 0) + .def("_write_cache", &contourpy::ThreadedContourGenerator::write_cache) + .def("create_contour", &contourpy::ThreadedContourGenerator::lines, create_contour_doc) + .def("create_filled_contour", &contourpy::ThreadedContourGenerator::filled, + create_filled_contour_doc) + .def("filled", &contourpy::ThreadedContourGenerator::filled, filled_doc) + .def("lines", &contourpy::ThreadedContourGenerator::lines, lines_doc) + .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 new file mode 100644 index 0000000000..92811ad0d0 --- /dev/null +++ b/contrib/python/contourpy/src/z_interp.cpp @@ -0,0 +1,19 @@ +#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 new file mode 100644 index 0000000000..76a7550be6 --- /dev/null +++ b/contrib/python/contourpy/src/z_interp.h @@ -0,0 +1,23 @@ +#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 |