: std::vector<XY>(),

_is_hole(is_hole),

_parent(0)

{

assert(!_is_hole && "Cannot add_child to a hole");

assert(child != 0 && "Null child ContourLine");

_children.push_back(child);

{

assert(is_hole() && "Cannot clear parent of non-hole");

assert(_parent != 0 && "Null parent ContourLine");

_parent = 0;

{

if (empty() || point != back())

std::vector<XY>::push_back(point);

{

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

assert(parent != 0 && "Null parent ContourLine");

_parent = parent;

{

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;

{

for (iterator line_it = begin(); line_it != end(); ++line_it) {

delete *line_it;

*line_it = 0;

}

std::vector<ContourLine*>::clear();

{

std::cout << "Contour of " << size() << " lines." << std::endl;

for (const_iterator it = begin(); it != end(); ++it)

(*it)->write();

: _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");

{

long i = quad % _nx;

long j = quad / _nx;

long index = (i-_istart) + (j-_jstart)*_x_chunk_points;

assert(i >= _istart && i < _istart + _x_chunk_points &&

"i-index outside chunk");

assert(j >= _jstart && j < _jstart + _y_chunk_points &&

"j-index outside chunk");

assert(index >= 0 && index < static_cast<long>(_lines.size()) &&

"ParentCache index outside chunk");

return index;

{

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

{

assert(!_lines.empty() &&

"Accessing ParentCache before it has been initialised");

long index = quad_to_index(quad);

if (_lines[index] == 0)

_lines[index] = (contour_line.is_hole() ? contour_line.get_parent()

: &contour_line);

const CoordinateArray& y,

const CoordinateArray& z,

const MaskArray& mask,

bool corner_mask,

long chunk_size)

: _x(x),

_y(y),

_z(z),

_nx(static_cast<long>(_x.dim(1))),

_ny(static_cast<long>(_x.dim(0))),

_n(_nx*_ny),

_corner_mask(corner_mask),

_chunk_size(chunk_size > 0 ? std::min(chunk_size, std::max(_nx, _ny)-1)

: std::max(_nx, _ny)-1),

_nxchunk(calc_chunk_count(_nx)),

_nychunk(calc_chunk_count(_ny)),

_chunk_count(_nxchunk*_nychunk),

_cache(new CacheItem[_n]),

_parent_cache(_nx,

chunk_size > 0 ? chunk_size+1 : _nx,

chunk_size > 0 ? chunk_size+1 : _ny)

{

assert(!_x.empty() && !_y.empty() && !_z.empty() && "Empty array");

assert(_y.dim(0) == _x.dim(0) && _y.dim(1) == _x.dim(1) &&

"Different-sized y and x arrays");

assert(_z.dim(0) == _x.dim(0) && _z.dim(1) == _x.dim(1) &&

"Different-sized z and x arrays");

assert((mask.empty() ||

(mask.dim(0) == _x.dim(0) && mask.dim(1) == _x.dim(1))) &&

"Different-sized mask and x arrays");

init_cache_grid(mask);

ContourLine& contour_line,

PyObject* vertices_list) const

{

assert(vertices_list != 0 && "Null python vertices_list");

// Convert ContourLine to python equivalent, and clear it.

npy_intp dims[2] = {static_cast<npy_intp>(contour_line.size()), 2};

numpy::array_view<double, 2> line(dims);

npy_intp i = 0;

for (ContourLine::const_iterator point = contour_line.begin();

point != contour_line.end(); ++point, ++i) {

line(i, 0) = point->x;

line(i, 1) = point->y;

}

if (PyList_Append(vertices_list, line.pyobj_steal())) {

Py_XDECREF(vertices_list);

throw std::runtime_error("Unable to add contour line to vertices_list");

}

contour_line.clear();

Contour& contour,

PyObject* vertices_list,

PyObject* codes_list) const

{

assert(vertices_list != 0 && "Null python vertices_list");

assert(codes_list != 0 && "Null python codes_list");

// 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();

npy_intp npoints = static_cast<npy_intp>(line.size() + 1);

for (children_it = children.begin(); children_it != children.end();

++children_it)

npoints += static_cast<npy_intp>((*children_it)->size() + 1);

npy_intp vertices_dims[2] = {npoints, 2};

numpy::array_view<double, 2> vertices(vertices_dims);

double* vertices_ptr = vertices.data();

npy_intp codes_dims[1] = {npoints};

numpy::array_view<unsigned char, 1> codes(codes_dims);

unsigned char* codes_ptr = codes.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.

}

if (PyList_Append(vertices_list, vertices.pyobj_steal()) ||

PyList_Append(codes_list, codes.pyobj_steal())) {

Py_XDECREF(vertices_list);

Py_XDECREF(codes_list);

contour.delete_contour_lines();

throw std::runtime_error("Unable to add contour line to vertices and codes lists");

}

delete *line_it;

*line_it = 0;

}

}

// Delete remaining contour lines.

contour.delete_contour_lines();

{

assert(point_count > 0 && "point count must be positive");

assert(_chunk_size > 0 && "Chunk size must be positive");

if (_chunk_size > 0) {

long count = (point_count-1) / _chunk_size;

if (count*_chunk_size < point_count-1)

++count;

assert(count >= 1 && "Invalid chunk count");

return count;

}

else

return 1;

const double& level)

{

assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&

"Quad index out of bounds");

assert(quad_edge.edge != Edge_None && "Invalid edge");

return interp(get_edge_point_index(quad_edge, true),

get_edge_point_index(quad_edge, false),

level);

ContourLine& contour_line,

QuadEdge& quad_edge,

const double& lower_level,

const double& upper_level,

unsigned int level_index,

const QuadEdge& start_quad_edge)

{

assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&

"Quad index out of bounds");

assert(quad_edge.edge != Edge_None && "Invalid edge");

assert(is_edge_a_boundary(quad_edge) && "Not a boundary edge");

assert((level_index == 1 || level_index == 2) &&

"level index must be 1 or 2");

assert(start_quad_edge.quad >= 0 && start_quad_edge.quad < _n &&

"Start quad index out of bounds");

assert(start_quad_edge.edge != Edge_None && "Invalid start edge");

// Only called for filled contours, so always updates _parent_cache.

unsigned int end_level = 0;

bool first_edge = true;

bool stop = false;

long& quad = quad_edge.quad;

while (true) {

// Levels of start and end points of quad_edge.

unsigned int start_level =

(first_edge ? Z_LEVEL(get_edge_point_index(quad_edge, true))

: end_level);

long end_point = get_edge_point_index(quad_edge, false);

end_level = Z_LEVEL(end_point);

if (level_index == 1) {

if (start_level <= level_index && end_level == 2) {

// Increasing z, switching levels from 1 to 2.

level_index = 2;

stop = true;

}

else if (start_level >= 1 && end_level == 0) {

// Decreasing z, keeping same level.

stop = true;

}

}

else { // level_index == 2

if (start_level <= level_index && end_level == 2) {

// Increasing z, keeping same level.

stop = true;

}

else if (start_level >= 1 && end_level == 0) {

// Decreasing z, switching levels from 2 to 1.

level_index = 1;

stop = true;

}

}

if (!first_edge && !stop && quad_edge == start_quad_edge)

// Return if reached start point of contour line. Do this before

// checking/setting VISITED flags as will already have been

// visited.

break;

switch (quad_edge.edge) {

case Edge_E:

assert(!VISITED_W(quad+1) && "Already visited");

_cache[quad+1] |= MASK_VISITED_W;

break;

case Edge_N:

assert(!VISITED_S(quad+_nx) && "Already visited");

_cache[quad+_nx] |= MASK_VISITED_S;

break;

case Edge_W:

assert(!VISITED_W(quad) && "Already visited");

_cache[quad] |= MASK_VISITED_W;

break;

case Edge_S:

assert(!VISITED_S(quad) && "Already visited");

_cache[quad] |= MASK_VISITED_S;

break;

case Edge_NE:

case Edge_NW:

case Edge_SW:

case Edge_SE:

assert(!VISITED_CORNER(quad) && "Already visited");

_cache[quad] |= MASK_VISITED_CORNER;

break;

default:

assert(0 && "Invalid Edge");

break;

}

if (stop) {

// Exiting boundary to enter interior.

contour_line.push_back(edge_interp(quad_edge,

level_index == 1 ? lower_level

: upper_level));

break;

}

move_to_next_boundary_edge(quad_edge);

// Just moved to new quad edge, so label parent of start of quad edge.

switch (quad_edge.edge) {

case Edge_W:

case Edge_SW:

case Edge_S:

case Edge_SE:

if (!EXISTS_SE_CORNER(quad))

_parent_cache.set_parent(quad, contour_line);

break;

case Edge_E:

case Edge_NE:

case Edge_N:

case Edge_NW:

if (!EXISTS_SW_CORNER(quad))

_parent_cache.set_parent(quad + 1, contour_line);

break;

default:

assert(0 && "Invalid edge");

break;

}

// Add point to contour.

contour_line.push_back(get_point_xy(end_point));

if (first_edge)

first_edge = false;

}

return level_index;

QuadEdge& quad_edge,

unsigned int level_index,

const double& level,

bool want_initial_point,

const QuadEdge* start_quad_edge,

unsigned int start_level_index,

bool set_parents)

{

assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&

"Quad index out of bounds.");

assert(quad_edge.edge != Edge_None && "Invalid edge");

assert((level_index == 1 || level_index == 2) &&

"level index must be 1 or 2");

assert((start_quad_edge == 0 ||

(start_quad_edge->quad >= 0 && start_quad_edge->quad < _n)) &&

"Start quad index out of bounds.");

assert((start_quad_edge == 0 || start_quad_edge->edge != Edge_None) &&

"Invalid start edge");

assert((start_level_index == 1 || start_level_index == 2) &&

"start level index must be 1 or 2");

long& quad = quad_edge.quad;

Edge& edge = quad_edge.edge;

if (want_initial_point)

contour_line.push_back(edge_interp(quad_edge, level));

CacheItem visited_mask = (level_index == 1 ? MASK_VISITED_1 : MASK_VISITED_2);

CacheItem saddle_mask = (level_index == 1 ? MASK_SADDLE_1 : MASK_SADDLE_2);

Dir dir = Dir_Straight;

while (true) {

assert(!EXISTS_NONE(quad) && "Quad does not exist");

assert(!(_cache[quad] & visited_mask) && "Quad already visited");

// Determine direction to move to next quad. If the quad is already

// labelled as a saddle quad then the direction is easily read from

// the cache. Otherwise the direction is determined differently

// depending on whether the quad is a corner quad or not.

if (_cache[quad] & saddle_mask) {

// Already identified as a saddle quad, so direction is easy.

dir = (SADDLE_LEFT(quad,level_index) ? Dir_Left : Dir_Right);

_cache[quad] |= visited_mask;

}

else if (EXISTS_ANY_CORNER(quad)) {

// Need z-level of point opposite the entry edge, as that

// determines whether contour turns left or right.

long point_opposite = -1;

switch (edge) {

case Edge_E:

point_opposite = (EXISTS_SE_CORNER(quad) ? POINT_SW

: POINT_NW);

break;

case Edge_N:

point_opposite = (EXISTS_NW_CORNER(quad) ? POINT_SW

: POINT_SE);

break;

case Edge_W:

point_opposite = (EXISTS_SW_CORNER(quad) ? POINT_SE

: POINT_NE);

break;

case Edge_S:

point_opposite = (EXISTS_SW_CORNER(quad) ? POINT_NW

: POINT_NE);

break;

case Edge_NE: point_opposite = POINT_SW; break;

case Edge_NW: point_opposite = POINT_SE; break;

case Edge_SW: point_opposite = POINT_NE; break;

case Edge_SE: point_opposite = POINT_NW; break;

default: assert(0 && "Invalid edge"); break;

}

assert(point_opposite != -1 && "Failed to find opposite point");

// Lower-level polygons (level_index == 1) always have higher

// values to the left of the contour. Upper-level contours

// (level_index == 2) are reversed, which is what the fancy XOR

// does below.

if ((Z_LEVEL(point_opposite) >= level_index) ^ (level_index == 2))

dir = Dir_Right;

else

dir = Dir_Left;

_cache[quad] |= visited_mask;

}

else {

// Calculate configuration of this quad.

long point_left = -1, point_right = -1;

switch (edge) {

case Edge_E: point_left = POINT_SW; point_right = POINT_NW; break;

case Edge_N: point_left = POINT_SE; point_right = POINT_SW; break;

case Edge_W: point_left = POINT_NE; point_right = POINT_SE; break;

case Edge_S: point_left = POINT_NW; point_right = POINT_NE; break;

default: assert(0 && "Invalid edge"); break;

}

unsigned int config = (Z_LEVEL(point_left) >= level_index) << 1 |

(Z_LEVEL(point_right) >= level_index);

// Upper level (level_index == 2) polygons are reversed compared to

// lower level ones, i.e. higher values on the right rather than

// the left.

if (level_index == 2)

config = 3 - config;

// Calculate turn direction to move to next quad along contour line.

if (config == 1) {

// New saddle quad, set up cache bits for it.

double zmid = 0.25*(get_point_z(POINT_SW) +

get_point_z(POINT_SE) +

get_point_z(POINT_NW) +

get_point_z(POINT_NE));

_cache[quad] |= (level_index == 1 ? MASK_SADDLE_1 : MASK_SADDLE_2);

if ((zmid > level) ^ (level_index == 2)) {

dir = Dir_Right;

}

else {

dir = Dir_Left;

_cache[quad] |= (level_index == 1 ? MASK_SADDLE_LEFT_1

: MASK_SADDLE_LEFT_2);

}

if (edge == Edge_N || edge == Edge_E) {

// Next visit to this quad must start on S or W.

_cache[quad] |= (level_index == 1 ? MASK_SADDLE_START_SW_1

: MASK_SADDLE_START_SW_2);

}

}

else {

// Normal (non-saddle) quad.

dir = (config == 0 ? Dir_Left

: (config == 3 ? Dir_Right : Dir_Straight));

_cache[quad] |= visited_mask;

}

}

// Use dir to determine exit edge.

edge = get_exit_edge(quad_edge, dir);

if (set_parents) {

if (edge == Edge_E)

_parent_cache.set_parent(quad+1, contour_line);

else if (edge == Edge_W)

_parent_cache.set_parent(quad, contour_line);

}

// Add new point to contour line.

contour_line.push_back(edge_interp(quad_edge, level));

// Stop if reached boundary.

if (is_edge_a_boundary(quad_edge))

break;

move_to_next_quad(quad_edge);

assert(quad_edge.quad >= 0 && quad_edge.quad < _n &&

"Quad index out of bounds");

// Return if reached start point of contour line.

if (start_quad_edge != 0 &&

quad_edge == *start_quad_edge &&

level_index == start_level_index)

break;

}

long& ichunk,

long& jchunk,

long& istart,

long& iend,

long& jstart,

long& jend)

{

assert(ijchunk >= 0 && ijchunk < _chunk_count && "ijchunk out of bounds");

ichunk = ijchunk % _nxchunk;

jchunk = ijchunk / _nxchunk;

istart = ichunk*_chunk_size;

iend = (ichunk == _nxchunk-1 ? _nx : (ichunk+1)*_chunk_size);

jstart = jchunk*_chunk_size;

jend = (jchunk == _nychunk-1 ? _ny : (jchunk+1)*_chunk_size);

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

//

long point1, point2, point3;

Edge edge12, edge23, edge31;

switch (_cache[quad] & MASK_EXISTS) {