#ifndef CONTOURPY_BASE_IMPL_H #define CONTOURPY_BASE_IMPL_H #include "base.h" #include "converter.h" #include "util.h" #include 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 BaseContourGenerator::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(std::ceil((_nx-1.0) / _x_chunk_size))), _ny_chunks(static_cast(std::ceil((_ny-1.0) / _y_chunk_size))), _n_chunks(_nx_chunks*_ny_chunks), _corner_mask(corner_mask), _line_type(line_type), _fill_type(fill_type), _quad_as_tri(quad_as_tri), _z_interp(z_interp), _cache(new CacheItem[_n]), _filled(false), _lower_level(0.0), _upper_level(0.0), _identify_holes(false), _output_chunked(false), _direct_points(false), _direct_line_offsets(false), _direct_outer_offsets(false), _outer_offsets_into_points(false), _nan_separated(false), _return_list_count(0) { if (_x.ndim() != 2 || _y.ndim() != 2 || _z.ndim() != 2) throw std::invalid_argument("x, y and z must all be 2D arrays"); if (_x.shape(1) != _nx || _x.shape(0) != _ny || _y.shape(1) != _nx || _y.shape(0) != _ny) throw std::invalid_argument("x, y and z arrays must have the same shape"); if (_nx < 2 || _ny < 2) throw std::invalid_argument("x, y and z must all be at least 2x2 arrays"); if (mask.ndim() != 0) { // ndim == 0 if mask is not set, which is valid. if (mask.ndim() != 2) throw std::invalid_argument("mask array must be a 2D array"); if (mask.shape(1) != _nx || mask.shape(0) != _ny) throw std::invalid_argument( "If mask is set it must be a 2D array with the same shape as z"); } if (!supports_line_type(line_type)) throw std::invalid_argument("Unsupported LineType"); if (!supports_fill_type(fill_type)) throw std::invalid_argument("Unsupported FillType"); if (x_chunk_size < 0 || y_chunk_size < 0) throw std::invalid_argument("x_chunk_size and y_chunk_size cannot be negative"); if (_z_interp == ZInterp::Log) { const bool* mask_ptr = (mask.ndim() == 0 ? nullptr : mask.data()); for (index_t point = 0; point < _n; ++point) { if ( (mask_ptr == nullptr || !mask_ptr[point]) && _zptr[point] <= 0.0) throw std::invalid_argument("z values must be positive if using ZInterp.Log"); } } init_cache_grid(mask); } template BaseContourGenerator::~BaseContourGenerator() { delete [] _cache; } template typename BaseContourGenerator::ZLevel BaseContourGenerator::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 double BaseContourGenerator::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 void BaseContourGenerator::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 void BaseContourGenerator::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 void BaseContourGenerator::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 FillType BaseContourGenerator::default_fill_type() { FillType fill_type = FillType::OuterOffset; assert(supports_fill_type(fill_type)); return fill_type; } template LineType BaseContourGenerator::default_line_type() { LineType line_type = LineType::Separate; assert(supports_line_type(line_type)); return line_type; } template py::sequence BaseContourGenerator::filled(double lower_level, double upper_level) { if (lower_level >= upper_level) throw std::invalid_argument("upper_level must be larger than lower_level"); _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); _nan_separated = false; _return_list_count = (_fill_type == FillType::ChunkCombinedCodeOffset || _fill_type == FillType::ChunkCombinedOffsetOffset) ? 3 : 2; return march_wrapper(); } template index_t BaseContourGenerator::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 bool BaseContourGenerator::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 bool BaseContourGenerator::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 index_t BaseContourGenerator::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 py::tuple BaseContourGenerator::get_chunk_count() const { return py::make_tuple(_ny_chunks, _nx_chunks); } template void BaseContourGenerator::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 py::tuple BaseContourGenerator::get_chunk_size() const { return py::make_tuple(_y_chunk_size, _x_chunk_size); } template bool BaseContourGenerator::get_corner_mask() const { return _corner_mask; } template FillType BaseContourGenerator::get_fill_type() const { return _fill_type; } template index_t BaseContourGenerator::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 double BaseContourGenerator::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 LineType BaseContourGenerator::get_line_type() const { return _line_type; } template double BaseContourGenerator::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 double BaseContourGenerator::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 index_t BaseContourGenerator::get_n_chunks() const { return _n_chunks; } template void BaseContourGenerator::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 double BaseContourGenerator::get_point_x(index_t point) const { assert(point >= 0 && point < _n && "point index out of bounds"); return _xptr[point]; } template double BaseContourGenerator::get_point_y(index_t point) const { assert(point >= 0 && point < _n && "point index out of bounds"); return _yptr[point]; } template double BaseContourGenerator::get_point_z(index_t point) const { assert(point >= 0 && point < _n && "point index out of bounds"); return _zptr[point]; } template bool BaseContourGenerator::get_quad_as_tri() const { return _quad_as_tri; } template ZInterp BaseContourGenerator::get_z_interp() const { return _z_interp; } template bool BaseContourGenerator::has_direct_line_offsets() const { return _direct_line_offsets; } template bool BaseContourGenerator::has_direct_outer_offsets() const { return _direct_outer_offsets; } template bool BaseContourGenerator::has_direct_points() const { return _direct_points; } template void BaseContourGenerator::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 void BaseContourGenerator::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 void BaseContourGenerator::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 void BaseContourGenerator::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 bool BaseContourGenerator::is_filled() const { return _filled; } template bool BaseContourGenerator::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 bool BaseContourGenerator::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 bool BaseContourGenerator::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 void BaseContourGenerator::line(const Location& start_location, ChunkLocal& local) { // start_location.on_boundary indicates starts (and therefore also finishes) assert(is_quad_in_chunk(start_location.quad, local)); Location location = start_location; count_t point_count = 0; // Insert nan if required before start of new line. if (_nan_separated and local.pass > 0 && local.line_count > 0) { *local.points.current++ = Util::nan; *local.points.current++ = Util::nan; } // finished == true indicates closed line loop. bool finished = follow_interior(location, start_location, local, point_count); if (local.pass > 0) { assert(local.line_offsets.current == local.line_offsets.start + local.line_count); *local.line_offsets.current++ = local.total_point_count; } if (local.pass == 0 && !start_location.on_boundary && !finished) // An internal start that isn't a line loop is part of a line strip that starts on a // boundary and will be traced later. Do not count it as a valid start in pass 0 and remove // the first point or it will be duplicated by the correct boundary-started line later. point_count--; else local.line_count++; local.total_point_count += point_count; } template py::sequence BaseContourGenerator::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 || _line_type == LineType::ChunkCombinedNan) ? 1 : 2; _nan_separated = (_line_type == LineType::ChunkCombinedNan); if (_nan_separated) Util::ensure_nan_loaded(); return march_wrapper(); } template void BaseContourGenerator::march_chunk( ChunkLocal& local, std::vector& return_lists) { for (local.pass = 0; local.pass < 2; ++local.pass) { bool ignore_holes = (_identify_holes && local.pass == 1); index_t j_final_start = local.jstart; for (index_t j = local.jstart; j <= local.jend; ++j) { index_t quad = local.istart + j*_nx; if (NO_MORE_STARTS(quad)) break; if (NO_STARTS_IN_ROW(quad)) continue; // Want to count number of starts in this row, so store how many starts at start of row. auto prev_start_count = (_identify_holes ? local.line_count - local.hole_count : local.line_count); for (index_t i = local.istart; i <= local.iend; ++i, ++quad) { if (!ANY_START(quad)) continue; assert(EXISTS_ANY(quad)); if (_filled) { if (START_BOUNDARY_S(quad)) closed_line_wrapper(Location(quad, 1, _nx, Z_SW == 2, true), Outer, local); if (START_BOUNDARY_W(quad)) closed_line_wrapper(Location(quad, -_nx, 1, Z_NW == 2, true), Outer, local); if (START_CORNER(quad)) { switch (EXISTS_ANY_CORNER(quad)) { case MASK_EXISTS_NE_CORNER: closed_line_wrapper( Location(quad, -_nx+1, _nx+1, Z_NW == 2, true), Outer, local); break; case MASK_EXISTS_NW_CORNER: closed_line_wrapper( Location(quad, _nx+1, _nx-1, Z_SW == 2, true), Outer, local); break; case MASK_EXISTS_SE_CORNER: closed_line_wrapper( Location(quad, -_nx-1, -_nx+1, Z_NE == 2, true), Outer, local); break; default: assert(EXISTS_SW_CORNER(quad)); if (!ignore_holes) closed_line_wrapper( Location(quad, _nx-1, -_nx-1, false, true), Hole, local); break; } } if (START_N(quad)) closed_line_wrapper(Location(quad, -_nx, 1, Z_NW > 0, false), Outer, local); if (ignore_holes) continue; if (START_E(quad)) closed_line_wrapper(Location(quad, -1, -_nx, Z_NE > 0, false), Hole, local); if (START_HOLE_N(quad)) closed_line_wrapper(Location(quad, -1, -_nx, false, true), Hole, local); } else { // !_filled if (START_BOUNDARY_S(quad)) line(Location(quad, _nx, -1, false, true), local); if (START_BOUNDARY_W(quad)) line(Location(quad, 1, _nx, false, true), local); if (START_BOUNDARY_E(quad)) line(Location(quad, -1, -_nx, false, true), local); if (START_BOUNDARY_N(quad)) line(Location(quad, -_nx, 1, false, true), local); if (START_E(quad)) line(Location(quad, -1, -_nx, false, false), local); if (START_N(quad)) line(Location(quad, -_nx, 1, false, false), local); if (START_CORNER(quad)) { index_t forward, left; switch (EXISTS_ANY_CORNER(quad)) { case MASK_EXISTS_NE_CORNER: forward = _nx+1; left = _nx-1; break; case MASK_EXISTS_NW_CORNER: forward = _nx-1; left = -_nx-1; break; case MASK_EXISTS_SE_CORNER: forward = -_nx+1; left = _nx+1; break; default: assert(EXISTS_SW_CORNER(quad)); forward = -_nx-1; left = -_nx+1; break; } line(Location(quad, forward, left, false, true), local); } } // _filled } // i // Number of starts at end of row. auto start_count = (_identify_holes ? local.line_count - local.hole_count : local.line_count); if (start_count > prev_start_count) j_final_start = j; else _cache[local.istart + j*_nx] |= MASK_NO_STARTS_IN_ROW; } // j if (j_final_start < local.jend) _cache[local.istart + (j_final_start+1)*_nx] |= MASK_NO_MORE_STARTS; if (_nan_separated && local.line_count > 1) { // If _nan_separated, each line after the first has an extra nan to separate it from the // previous line's points. If we were returning line offsets to the caller then this // would need to occur in line() where the line_count is incremented. But as we are not // returning line offsets it is faster just to add the extra here all at once. local.total_point_count += local.line_count - 1; } if (local.pass == 0) { if (local.total_point_count == 0) { local.points.clear(); local.line_offsets.clear(); local.outer_offsets.clear(); break; // Do not need pass 1. } // Create arrays for points, line_offsets and optionally outer_offsets. Arrays may be // either C++ vectors or Python NumPy arrays. Want to group creation of the latter as // threaded code needs to lock creation of these to limit access to a single thread. if (_direct_points || _direct_line_offsets || _direct_outer_offsets) { typename Derived::Lock lock(static_cast(*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(*this)); for (auto& list : return_lists) list[local.chunk] = py::none(); } } else if (_filled) static_cast(this)->export_filled(local, return_lists); else static_cast(this)->export_lines(local, return_lists); } template py::sequence BaseContourGenerator::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 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(this)->march(return_lists); // Return to python objects. if (_return_list_count == 1) { assert(!_filled); if (_line_type == LineType::Separate) return return_lists[0]; else { assert(_line_type == LineType::ChunkCombinedNan); return py::make_tuple(return_lists[0]); } } else if (_return_list_count == 2) return py::make_tuple(return_lists[0], return_lists[1]); else { assert(_return_list_count == 3); return py::make_tuple(return_lists[0], return_lists[1], return_lists[2]); } } template void BaseContourGenerator::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 void BaseContourGenerator::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 bool BaseContourGenerator::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 bool BaseContourGenerator::supports_line_type(LineType line_type) { switch (line_type) { case LineType::Separate: case LineType::SeparateCode: case LineType::ChunkCombinedCode: case LineType::ChunkCombinedOffset: case LineType::ChunkCombinedNan: return true; default: return false; } } template void BaseContourGenerator::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 void BaseContourGenerator::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 BaseContourGenerator::ZLevel BaseContourGenerator::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