#define NO_IMPORT_ARRAY

#include "mplutils.h"

#include "_contour.h"

#include <algorithm>

#define POINT_SW (quad)

#define POINT_SE (quad+1)

#define POINT_NW (quad+_nx)

#define POINT_NE (quad+_nx+1)

#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.

#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.

#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.

#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)

#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))

QuadEdge::QuadEdge()

: quad(-1), edge(Edge_None)

{}

QuadEdge::QuadEdge(long quad_, Edge edge_)

: quad(quad_), edge(edge_)

{}

bool QuadEdge::operator<(const QuadEdge& other) const

{

if (quad != other.quad)

return quad < other.quad;

else

return edge < other.edge;

}

bool QuadEdge::operator==(const QuadEdge& other) const

{

return quad == other.quad && edge == other.edge;

}

bool QuadEdge::operator!=(const QuadEdge& other) const

{

return !operator==(other);

}

std::ostream& operator<<(std::ostream& os, const QuadEdge& quad_edge)

{

return os << quad_edge.quad << ' ' << quad_edge.edge;

}

XY::XY()

{}

XY::XY(const double& x_, const double& y_)

: x(x_), y(y_)

{}

bool XY::operator==(const XY& other) const

{

return x == other.x && y == other.y;

}

bool XY::operator!=(const XY& other) const

{

return x != other.x || y != other.y;

}

XY XY::operator*(const double& multiplier) const

{

return XY(x*multiplier, y*multiplier);

}

const XY& XY::operator+=(const XY& other)

{

x += other.x;

y += other.y;

return *this;

}

const XY& XY::operator-=(const XY& other)

{

x -= other.x;

y -= other.y;

return *this;

}

XY XY::operator+(const XY& other) const

{

return XY(x + other.x, y + other.y);

}

XY XY::operator-(const XY& other) const

{

return XY(x - other.x, y - other.y);

}

std::ostream& operator<<(std::ostream& os, const XY& xy)

{

return os << '(' << xy.x << ' ' << xy.y << ')';

}

ContourLine::ContourLine(bool is_hole)

{}

void ContourLine::add_child(ContourLine* child)

}

void ContourLine::clear_parent()

}

const ContourLine::Children& ContourLine::get_children() const

{

assert(!_is_hole && "Cannot get_children of a hole");

return _children;

}

const ContourLine* ContourLine::get_parent() const

{

assert(_is_hole && "Cannot get_parent of a non-hole");

return _parent;

}

ContourLine* ContourLine::get_parent()

{

assert(_is_hole && "Cannot get_parent of a non-hole");

return _parent;

}

bool ContourLine::is_hole() const

{

return _is_hole;

}

void ContourLine::push_back(const XY& point)

}

void ContourLine::set_parent(ContourLine* parent)

}

void ContourLine::write() const

}

Contour::Contour()

{}

Contour::~Contour()

{

delete_contour_lines();

}

void Contour::delete_contour_lines()

}

void Contour::write() const

}

ParentCache::ParentCache(long nx, long x_chunk_points, long y_chunk_points)

}

ContourLine* ParentCache::get_parent(long quad)

{

long index = quad_to_index(quad);

ContourLine* parent = _lines[index];

while (parent == 0) {

index -= _x_chunk_points;

assert(index >= 0 && "Failed to find parent in chunk ParentCache");

parent = _lines[index];

}

assert(parent != 0 && "Failed to find parent in chunk ParentCache");

return parent;

}

long ParentCache::quad_to_index(long quad) const

}

void ParentCache::set_chunk_starts(long istart, long jstart)

}

void ParentCache::set_parent(long quad, ContourLine& contour_line)

}

QuadContourGenerator::QuadContourGenerator(const CoordinateArray& x,

}

QuadContourGenerator::~QuadContourGenerator()

{

delete [] _cache;

}

void QuadContourGenerator::append_contour_line_to_vertices_and_codes(

}

void QuadContourGenerator::append_contour_to_vertices_and_codes(

}

long QuadContourGenerator::calc_chunk_count(long point_count) const

}

PyObject* QuadContourGenerator::create_contour(const double& level)

{

init_cache_levels(level, level);

PyObject* vertices_list = PyList_New(0);

if (vertices_list == 0)

throw std::runtime_error("Failed to create Python list");

PyObject* codes_list = PyList_New(0);

if (codes_list == 0) {

Py_XDECREF(vertices_list);

throw std::runtime_error("Failed to create Python list");

}

// Lines that start and end on boundaries.

long ichunk, jchunk, istart, iend, jstart, jend;

for (long ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) {

get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart, jend);

for (long j = jstart; j < jend; ++j) {

long quad_end = iend + j*_nx;

for (long quad = istart + j*_nx; quad < quad_end; ++quad) {

if (EXISTS_NONE(quad) || VISITED(quad,1)) continue;

if (BOUNDARY_S(quad) && Z_SW >= 1 && Z_SE < 1 &&

start_line(vertices_list, 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 (long ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) {

get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart, jend);

for (long j = jstart; j < jend; ++j) {

long quad_end = iend + j*_nx;

for (long quad = istart + j*_nx; quad < quad_end; ++quad) {

if (EXISTS_NONE(quad) || VISITED(quad,1))

continue;

Edge start_edge = get_start_edge(quad, 1);

if (start_edge == Edge_None)

continue;

QuadEdge quad_edge(quad, start_edge);

QuadEdge start_quad_edge(quad_edge);

// To obtain output identical to that produced by legacy code,

// sometimes need to ignore the first point and add it on the

// end instead.

bool ignore_first = (start_edge == Edge_N);

follow_interior(contour_line, quad_edge, 1, level,

!ignore_first, &start_quad_edge, 1, false);

if (ignore_first && !contour_line.empty())

contour_line.push_back(contour_line.front());

append_contour_line_to_vertices_and_codes(

contour_line, vertices_list, codes_list);

// Repeat if saddle point but not visited.

if (SADDLE(quad,1) && !VISITED(quad,1))

--quad;

}

}

}

PyObject* tuple = PyTuple_New(2);

if (tuple == 0) {

Py_XDECREF(vertices_list);

Py_XDECREF(codes_list);

throw std::runtime_error("Failed to create Python tuple");

}

// No error checking here as filling in a brand new pre-allocated tuple.

PyTuple_SET_ITEM(tuple, 0, vertices_list);

PyTuple_SET_ITEM(tuple, 1, codes_list);

return tuple;

}

PyObject* QuadContourGenerator::create_filled_contour(const double& lower_level,

const double& upper_level)

{

init_cache_levels(lower_level, upper_level);

Contour contour;

PyObject* vertices_list = PyList_New(0);

if (vertices_list == 0)

throw std::runtime_error("Failed to create Python list");

PyObject* codes_list = PyList_New(0);

if (codes_list == 0) {

Py_XDECREF(vertices_list);

throw std::runtime_error("Failed to create Python list");

}

long ichunk, jchunk, istart, iend, jstart, jend;

for (long ijchunk = 0; ijchunk < _chunk_count; ++ijchunk) {

get_chunk_limits(ijchunk, ichunk, jchunk, istart, iend, jstart, jend);

_parent_cache.set_chunk_starts(istart, jstart);

for (long j = jstart; j < jend; ++j) {

long quad_end = iend + j*_nx;

for (long quad = istart + j*_nx; quad < quad_end; ++quad) {

if (!EXISTS_NONE(quad))

single_quad_filled(contour, quad, lower_level, upper_level);

}

}

// Clear VISITED_W and VISITED_S flags that are reused by later chunks.

if (jchunk < _nychunk-1) {

long quad_end = iend + jend*_nx;

for (long quad = istart + jend*_nx; quad < quad_end; ++quad)

_cache[quad] &= ~MASK_VISITED_S;

}

if (ichunk < _nxchunk-1) {

long quad_end = iend + jend*_nx;

for (long quad = iend + jstart*_nx; quad < quad_end; quad += _nx)

_cache[quad] &= ~MASK_VISITED_W;

}

// Create python objects to return for this chunk.

append_contour_to_vertices_and_codes(contour, vertices_list, codes_list);

}

PyObject* tuple = PyTuple_New(2);

if (tuple == 0) {

Py_XDECREF(vertices_list);

Py_XDECREF(codes_list);

throw std::runtime_error("Failed to create Python tuple");

}

// No error checking here as filling in a brand new pre-allocated tuple.

PyTuple_SET_ITEM(tuple, 0, vertices_list);

PyTuple_SET_ITEM(tuple, 1, codes_list);

return tuple;

}

XY QuadContourGenerator::edge_interp(const QuadEdge& quad_edge,

}

unsigned int QuadContourGenerator::follow_boundary(

}

void QuadContourGenerator::follow_interior(ContourLine& contour_line,