//

// Created by Alan Freitas on 2020-07-07.

//

#include <algorithm>

#include <cmath>

#include <matplot/axes_objects/contours.h>

#include <matplot/axes_objects/histogram.h>

#include <matplot/core/axes_type.h>

#include <matplot/freestanding/plot.h>

#include <matplot/util/common.h>

#include <numeric>

#include <regex>

#include <sstream>

#include <unordered_set>

namespace matplot {

contours::contours(class axes_type *parent, const vector_2d &X,

const vector_2d &Y, const vector_2d &Z,

std::string_view line_spec)

: axes_object(parent), line_spec_(this, line_spec), X_data_(X),

Y_data_(Y), Z_data_(Z) {

initialize_preprocessed_data();

contour_generator_ =

QuadContourGenerator(X_data_, Y_data_, Z_data_, _corner_mask,

static_cast<long>(nchunk_));

}

contours::contours(class axes_type *parent, const vector_2d &Z,

std::string_view line_spec)

: axes_object(parent), line_spec_(this, line_spec), Z_data_(Z) {

initialize_preprocessed_data();

contour_generator_ =

QuadContourGenerator(X_data_, Y_data_, Z_data_, _corner_mask,

static_cast<long>(nchunk_));

}

std::string contours::set_variables_string() {

return " set style textbox opaque margins 0.5, 0.5 fc bgnd noborder "

"linewidth 1.0\n";

}

std::string contours::plot_string() {

make_sure_data_is_preprocessed();

// double zmax_ = zmax();

// double zmin_ = zmin();

auto [min_it, max_it] =

std::minmax_element(levels_.begin(), levels_.end());

double contour_min_level = *min_it;

double contour_max_level = *max_it;

std::stringstream ss;

if (filled_) {

auto [lower_levels, upper_levels] = get_lowers_and_uppers();

// Command for background filled curve

// The background polygon with the whole area has its level defined

// by the largest polygon. Whatever level is outside this largest

// polygon is the level for the "background". The largest polygon is

// always touching the background because because some polygon

// larger than it would always be touching it. This background

// polygon will appear if we define a level < z_min.

// We plot background only if one of the levels is below zmin

bool plot_background = lower_levels[0] < zmin_;

if (plot_background) {

auto &largest_segment_with_children = line_segments_[0];

auto &largest_segment =

std::get<0>(largest_segment_with_children);

size_t line_index = std::get<0>(largest_segment);

size_t segment_begin = std::get<1>(largest_segment);

size_t segment_end = std::get<2>(largest_segment);

bool parent_is_lower_level =

is_lower_level(line_index, segment_begin, segment_end);

auto previous_color = line_spec_.color();

bool previous_color_manual = line_spec_.user_color();

if (parent_is_lower_level) {

// background is 1 lower than parent on lower level

size_t level_index = line_index > 0 ? line_index - 1 : 0;

double background_z_level = lower_levels[level_index];

line_spec_.color(parent_->colormap_interpolation(

background_z_level, contour_min_level,

contour_max_level));

} else {

// background is 1 higher than parent on upper level

size_t level_index = line_index > 0 ? line_index - 1 : 0;

double background_z_level = upper_levels[level_index];

line_spec_.color(parent_->colormap_interpolation(

background_z_level, contour_min_level,

contour_max_level));

}

std::string ls =

" '-' with filledcurve " +

line_spec_.plot_string(

line_spec::style_to_plot::plot_line_only, false);

// filledcurves need to use the palette to initialize the

// colorbox

ls =

std::regex_replace(ls, std::regex(" linecolor rgb +[^ ]+ "),

" linecolor palette ");

ss << ls;

line_spec_.color(previous_color);

line_spec_.user_color(previous_color_manual);

}

// Other curves

for (size_t i = 0; i < line_segments_.size(); ++i) {

if (i != 0 || plot_background) {

ss << ",";

}

parent_and_children_type &parent_and_children_segment =

line_segments_[i];

line_segment_type &parent_segment =

std::get<0>(parent_and_children_segment);

size_t line_index = std::get<0>(parent_segment);

size_t segment_begin = std::get<1>(parent_segment);

size_t segment_end = std::get<2>(parent_segment);

bool is_ll =

is_lower_level(line_index, segment_begin, segment_end);

double segment_z_level =

is_ll ? lower_levels[line_index] : upper_levels[line_index];

auto previous_color = line_spec_.color();

bool previous_color_manual = line_spec_.user_color();

line_spec_.color(parent_->colormap_interpolation(

segment_z_level, contour_min_level, contour_max_level));

std::string ls =

" '-' with filledcurve " +

line_spec_.plot_string(

line_spec::style_to_plot::plot_line_only, false);

ls =

std::regex_replace(ls, std::regex(" linecolor rgb +[^ ]+ "),

" linecolor palette ");

ss << ls;

line_spec_.color(previous_color);

line_spec_.user_color(previous_color_manual);

// create filledcurves for the children

std::vector<line_segment_type> &children_segments =

std::get<1>(parent_and_children_segment);

for (size_t j = 0; j < children_segments.size(); ++j) {

ss << ",";

line_segment_type &child_segment = children_segments[j];

line_index = std::get<0>(child_segment);

segment_begin = std::get<1>(child_segment);

segment_end = std::get<2>(child_segment);

// upper level

// might need to fix that later:

// children is not always the upper level

// it might be the other way around

// we need to check which is clockwise before deciding on

// that

const bool is_ll2 =

is_lower_level(line_index, segment_begin, segment_end);

const double segment_z_level2 =

is_ll2 ? lower_levels[line_index]

: upper_levels[line_index];

line_spec_.color(parent_->colormap_interpolation(

segment_z_level2, contour_min_level,

contour_max_level));

std::string ls2 =

" '-' with filledcurve " +

line_spec_.plot_string(

line_spec::style_to_plot::plot_line_only, false);

ls2 = std::regex_replace(

ls2, std::regex(" linecolor rgb +[^ ]+ "),

" linecolor palette ");

ss << ls2;

}

line_spec_.color(previous_color);

line_spec_.user_color(previous_color_manual);

}

// check if Z has nans

bool z_has_nans = false;

for (size_t i = 0; !z_has_nans && i < Z_data_.size(); ++i) {

for (size_t j = 0; !z_has_nans && j < Z_data_[i].size(); ++j) {

if (!std::isfinite(Z_data_[i][j])) {

z_has_nans = true;

}

}

}

if (z_has_nans) {

// include another filled curve with background color

// this curve will hide everything that has nan in it

auto previous_color = line_spec_.color();

bool previous_color_manual = line_spec_.user_color();

line_spec_.color(parent_->color());

std::string ls =

", '-' with filledcurve " +

line_spec_.plot_string(

line_spec::style_to_plot::plot_line_only, false);

// this one does not use palletes. it's always the background

// color

ss << ls;

line_spec_.color(previous_color);

line_spec_.user_color(previous_color_manual);

}

}

// plot normal contour lines

bool first_non_empty_line = true;

for (size_t i = 0; i < lines_.size(); ++i) {

if (lines_[i].first.empty()) {

continue;

}

if (!first_non_empty_line || filled_) {

ss << ",";

} else {

first_non_empty_line = false;

}

auto previous_color = line_spec_.color();

bool previous_color_manual = line_spec_.user_color();

if (!previous_color_manual && filled_ &&

!colormap_line_when_filled_) {

line_spec_.color("black");

}

ss << " '-' ";

std::string ls = line_spec_.plot_string(

line_spec::style_to_plot::plot_line_only);

// we might need to use the palette to draw colors

if (!previous_color_manual) {

if (!filled_ || colormap_line_when_filled_) {

ls = std::regex_replace(

ls, std::regex(" linecolor rgb +[^ ]+ "),

" linecolor palette ");

}

}

ss << ls;

if (!previous_color_manual) {

line_spec_.color(previous_color);

line_spec_.user_color(previous_color_manual);

}

}

if (contour_text_) {

for (size_t i = 0; i < lines_.size(); ++i) {

if (lines_[i].first.empty()) {

continue;

}

ss << ", '-' with labels ";

if (!filled_ && iequals(font_weight_, "normal")) {

ss << " boxed ";

}

ss << " textcolor '" << to_string(font_color_) << "' ";

ss << " font '" << font() << "," << font_size() << "' ";

ss << " rotate variable ";

}

}

return ss.str();

}

void contours::make_sure_data_is_preprocessed() {

if (levels_.empty() || n_levels_ == 0 || layers_.empty()) {

process_contour_levels();

}

if (lines_.empty()) {

process_all_segs_and_all_kinds();

}

}

void contours::clear_preprocessed_data() {

layers_.clear();

if (!manual_n_levels_) {

n_levels_ = 0;

}

if (!manual_levels_) {

levels_.clear();

}

lines_.clear();

codes_.clear();

}

std::vector<double> contours::determine_contour_levels(double z_min,

double z_max,

size_t n_levels_,

extend_option ext) {

// generate levels with bin picker and remove levels out of range

// we start trying to generate n + 2 levels to account for the ones

// we will usually remove

size_t target_n_edges = n_levels_ + 2;

std::vector<double> levels = histogram::bin_picker(

z_min, z_max, target_n_edges, (z_max - z_min) / (target_n_edges));

auto remove_levels_out_of_range = [&](std::vector<double> &levels) {

while (!levels.empty() && levels.front() < z_min) {

levels.erase(levels.begin());

}

while (!levels.empty() && levels.back() > z_max) {

levels.pop_back();

}

};

remove_levels_out_of_range(levels);

// check how many levels we have left and if we need to include level =

// 0

bool zero_is_in_the_range = z_min <= 0.0 && z_max >= 0;

bool missing_zero =

zero_is_in_the_range &&

std::find(levels.begin(), levels.end(), 0.) == levels.end();

size_t excess =

levels.size() > n_levels_ ? levels.size() - n_levels_ : 0;

if (missing_zero || excess) {

// next target has one less because we will insert 0.0 in the range

target_n_edges -= missing_zero;

// next target has excess less

target_n_edges -= excess;

levels = histogram::bin_picker(z_min, z_max, target_n_edges,

(z_max - z_min) / (target_n_edges));

remove_levels_out_of_range(levels);

missing_zero =

zero_is_in_the_range &&

std::find(levels.begin(), levels.end(), 0.) == levels.end();

if (missing_zero) {

auto it = std::lower_bound(levels.begin(), levels.end(), 0.);

levels.insert(it, 0.);

}

}

// If we are going to extend later, we remove the extreme

// levels because we are going to include more extreme levels

// later.

if (ext == extend_option::min || ext == extend_option::both) {

levels.erase(levels.begin());

}

if (ext == extend_option::max || ext == extend_option::both) {

levels.pop_back();

}

if (levels.size() < 3) {

levels = histogram::bin_picker(z_min, z_max, target_n_edges,

(z_max - z_min) / (target_n_edges));

remove_levels_out_of_range(levels);

if (missing_zero) {

auto it = std::lower_bound(levels.begin(), levels.end(), 0.);

levels.insert(it, 0.);

}

}

return levels;

}

std::vector<double> contours::determine_contour_levels(const vector_2d &Z,

size_t n_levels_,

extend_option ext) {

double zmin_ = Z[0][0];

double zmax_ = Z[0][0];

for (const auto &row : Z) {

auto [row_min_it, row_max_it] =

std::minmax_element(row.begin(), row.end());

zmin_ = std::min(zmin_, *row_min_it);

zmax_ = std::max(zmax_, *row_max_it);

}

return determine_contour_levels(zmin_, zmax_, n_levels_, ext);

}

void contours::process_contour_levels() {

// Set automatic numer of levels and levels

if (levels_.empty()) {

if (n_levels_ == 0) {

n_levels_ = !filled_ ? 7 : 10;

}

// Something to do: Find a nicer algorithm (simple) for that later

// The same algorithm can be used to find ticks

levels_ =

determine_contour_levels(zmin(), zmax(), n_levels_, extend_);

} else {

if (n_levels_ == 0) {

n_levels_ = levels_.size();

}

}

// We need at least one level in the valid range

// If this is a filled contour, this is ok

if (!filled_) {

bool any_inside = false;

for (size_t i = 0; i < levels_.size(); ++i) {

if (levels_[i] > zmin() && levels_[i] < zmax()) {

any_inside = true;

break;

}

}

if (!any_inside) {

levels_ = {zmin()};

}

}

if (filled_) {

if ((extend_ == extend_option::neither && levels_.size() < 2) ||

levels_.empty()) {

throw std::logic_error(

"Filled contours require at least 2 levels");

}

}

if (levels_.size() > 1) {

for (size_t i = 0; i < levels_.size() - 1; ++i) {

if (levels_[i + 1] - levels_[i] <= 0.) {

throw std::logic_error("Contour levels must be increasing");

}

}

}

// Make a private _levels to include extended regions; we

// want to leave the original levels attribute unchanged.

// The plotting algorithm will use _levels only.

// (Colorbar needs this even for line contours.)

_levels = levels_;

// Extend minimum beyond zmin (for filled plots)

const bool log = parent_->z_axis().scale() == axis_type::axis_scale::log;

if (extend_ == extend_option::both || extend_ == extend_option::min) {

double lower = log ? 1e-250 : -1e250;

_levels.insert(_levels.begin(), lower);

}

// Extend maximum beyond zmax (for filled plots)

if (extend_ == extend_option::both || extend_ == extend_option::max) {

double upper = 1e250;

_levels.emplace_back(upper);

}

// Layers == levels if we are only plotting lines

if (!filled_) {

layers_ = levels_;

return;

}

// Layer values are mid-way between levels in screen space

if (log) {

// Avoid overflow by taking sqrt before multiplying.

layers_.resize(_levels.size() - 1);

for (size_t i = 0; i < _levels.size() - 1; ++i) {

layers_[i] = sqrt(_levels[i]) * sqrt(_levels[i + 1]);

}

} else {

layers_.resize(_levels.size() - 1);

for (size_t i = 0; i < _levels.size() - 1; ++i) {

layers_[i] = 0.5 * (_levels[i] + _levels[i + 1]);

}

}

}

double contours::zmin() { return zmin_; }

double contours::zmax() { return zmax_; }

std::pair<vector_1d, vector_1d> contours::get_lowers_and_uppers() {

vector_1d lowers(_levels.begin(), _levels.end() - 1);

if (lowers.empty()) {

return std::pair<vector_1d, vector_1d>{{}, {}};

}

if (zmin() == lowers[0]) {

// Include minimum values in lowest interval

// so we don't change levels_

if (parent_->z_axis().scale() == axis_type::axis_scale::log) {

lowers[0] = 0.99 * zmin();

} else {

lowers[0] -= 1;

}

}

vector_1d uppers(_levels.begin() + 1, _levels.end());

return std::make_pair(lowers, uppers);

}

std::string contours::legend_string(std::string_view title) {

auto [min_level_it, max_level_it] =

std::minmax_element(levels_.begin(), levels_.end());

double zmax = *max_level_it;

double zmin = *min_level_it;

std::stringstream ss;

for (size_t i = 0; i < lines_.size(); ++i) {

if (i != 0) {

ss << ",";

}

auto previous_color = line_spec_.color();

bool previous_color_manual = line_spec_.user_color();

if (!previous_color_manual) {

line_spec_.color(parent_->colormap_interpolation(

(zmax - zmin) - (levels_[i] - zmin), 0., zmax - zmin));

}

ss << " keyentry "

<< line_spec_.plot_string(

line_spec::style_to_plot::plot_line_only)

<< " title \"" << escape(title) << "\" ";

if (!previous_color_manual) {

line_spec_.color(previous_color);

line_spec_.user_color(previous_color_manual);

}

}

return ss.str();

}

/// If a line is lower (true) or upper level (false)

/// We have to know if it's lower or upper level

/// to decide its color.

bool contours::is_lower_level(size_t line_index, size_t segment_begin,

size_t segment_end) {

// The parent non-hole is not always the lower or upper level.

// That depends on whether the function is increasing

// or decreasing on that region.

// lower level <-> higher values on the left

// upper level <-> higher values on the right

// Find limits

double _xmax = xmax();

double _xmin = xmin();

double _ymax = ymax();

double _ymin = ymin();

// Take two points - outside the border when possible

double x1 = filled_lines_[line_index].first[segment_begin];

double x2 = filled_lines_[line_index].first[segment_begin + 1];

double y1 = filled_lines_[line_index].second[segment_begin];

double y2 = filled_lines_[line_index].second[segment_begin + 1];

auto is_on_border = [&]() {

return (x1 <= _xmin || x1 >= _xmax || x2 <= _xmin || x2 >= _xmax ||

y1 <= _ymin || y1 >= _ymax || y2 <= _ymin || y2 >= _ymax);

};

size_t sample_begin = segment_begin;

while (is_on_border() && sample_begin < segment_end - 1) {

++sample_begin;

x1 = filled_lines_[line_index].first[sample_begin];

x2 = filled_lines_[line_index].first[sample_begin + 1];

y1 = filled_lines_[line_index].second[sample_begin];

y2 = filled_lines_[line_index].second[sample_begin + 1];

}

double avg_x = 0.5 * (x1 + x2);

double avg_y = 0.5 * (y1 + y2);

bool x_is_increasing = x2 > x1;

bool y_is_increasing = y2 > y1;

// look for the grid position of (x > x1, y > y1) - NE

auto it_y =

std::find_if(Y_data_.begin(), Y_data_.end(),

[&](const auto &y_row) { return y_row[0] > avg_y; });

auto it_x = std::find_if(

X_data_[0].begin(), X_data_[0].end(),

[&](const double &x_row_value) { return x_row_value > avg_x; });

size_t n_row = it_y - Y_data_.begin();

size_t n_col = it_x - X_data_[0].begin();

// look at the left

// If x is increasing, the left is in the north

// - Do nothing because we are already at a position where y > avg_y

// If x is not increasing, the left is in the south

// - Try to reduce the n_row (our current grid position is NE)

if (!x_is_increasing && n_row > 0) {

n_row--;

}

// If y is increasing, the left is in the west

// - Try to reduce the n_col (our current grid position is NE)

// If y is not increasing, the left is in the east

// - Do nothing because we are already at a position where x > avg_x

if (y_is_increasing && n_col > 0) {

n_col--;

}

// Find the opposite grid position to compare whether higher

// values are on the left or on the right

// The logic is the same

size_t opposite_row = n_row;

size_t opposite_col = n_col;

if (x_is_increasing && opposite_row > 0) {

opposite_row--;

} else if (!x_is_increasing && (opposite_row < Y_data_.size() - 1)) {

opposite_row++;

}

if (y_is_increasing && opposite_col < X_data_[0].size() - 1) {

opposite_col++;

} else if (!y_is_increasing && opposite_col > 0) {

opposite_col--;

}

// if it increases

bool higher_values_on_left =

Z_data_[n_row][n_col] > Z_data_[opposite_row][opposite_col];

if (higher_values_on_left) {

// lower level

return true;

} else {

// else, upper level

return false;

}

}

/// Generate points to close a polygon if not closed yet

///

/// Some polygons might not be closed yet.

/// Large polygons that end on borders might not be going around

/// the borders properly as it should. That's fine for lines

/// but it's not fine for filled curves.

/// In these cases, we have to sanitize and continue clockwise or

/// anticlockwise around the borders until we reach the initial point.

///

/// If they start or end outside the borders, there is nothing

/// to do because we don't know the appropriate path to close

/// the polygon. For this reason, these interior areas are

/// already closed.

///

/// If they start and end at the borders, as we know it is

/// going clockwise, we know how to close it because we can

/// go around the borders. In these cases, the contourc

/// algorithm does not close the polygons for us.

std::pair<vector_1d, vector_1d>

contours::fill_border_jump(double start_x, double start_y, double end_x,

double end_y, double x_min, double x_max,

double y_min, double y_max, bool is_parent) {

// start_b = identify if starting border W,S,E,N

constexpr uint8_t NONE = 0;

constexpr uint8_t NORTH = 1;

constexpr uint8_t SOUTH = 2;

constexpr uint8_t WEST = 3;

constexpr uint8_t EAST = 4;

uint8_t xy1_border = NONE;

if (start_x <= x_min) {

xy1_border = WEST;

} else if (start_x >= x_max) {

xy1_border = EAST;

} else if (start_y <= y_min) {

xy1_border = SOUTH;

} else if (start_y >= y_max) {

xy1_border = NORTH;

}

if (xy1_border == NONE) {

return {};

}

// end_b = identify if closing border border W,S,E,N

uint8_t xy2_border = NONE;

if (end_x <= x_min) {

xy2_border = WEST;

} else if (end_x >= x_max) {

xy2_border = EAST;

} else if (end_y <= y_min) {

xy2_border = SOUTH;

} else if (end_y >= y_max) {

xy2_border = NORTH;

}

if (xy2_border == NONE) {

return {};

}

// parents go anticlockwise

// we need to know the direction we should use to fill the border

const bool clockwise = !is_parent;

std::pair<vector_1d, vector_1d> result;

if (clockwise) {

while (xy1_border != xy2_border) {

switch (xy1_border) {

case WEST:

// append NW

result.first.emplace_back(x_min);

result.second.emplace_back(y_max);

xy1_border = NORTH;

break;

case NORTH:

// append NE

result.first.emplace_back(x_max);

result.second.emplace_back(y_max);

xy1_border = EAST;

break;

case EAST:

// append SE

result.first.emplace_back(x_max);

result.second.emplace_back(y_min);

xy1_border = SOUTH;

break;

case SOUTH:

// append SW

result.first.emplace_back(x_min);

result.second.emplace_back(y_min);

xy1_border = WEST;

break;

default:

throw std::logic_error("Invalid direction");

}

}

} else {

// anticlockwise

while (xy1_border != xy2_border) {

switch (xy1_border) {

case WEST:

// append SW

result.first.emplace_back(x_min);

result.second.emplace_back(y_min);

xy1_border = SOUTH;

break;

case SOUTH:

// append SE

result.first.emplace_back(x_max);

result.second.emplace_back(y_min);

xy1_border = EAST;

break;

case EAST:

// append NE

result.first.emplace_back(x_max);

result.second.emplace_back(y_max);

xy1_border = NORTH;

break;

case NORTH:

// append NW

result.first.emplace_back(x_min);

result.second.emplace_back(y_max);

xy1_border = WEST;

break;

default:

throw std::logic_error("Invalid direction");

}

}

}

return result;

}

std::string contours::data_string() {

// If there is a jump from a border to the other, we need to complete

// the curve from one point to another to avoid filled curves that don't

// make sense.

double _xmax = xmax();

double _xmin = xmin();

double _ymax = ymax();

double _ymin = ymin();

auto is_border_jump = [&](double x1, double y1, double x2, double y2) {

constexpr uint8_t NONE = 0;

constexpr uint8_t NORTH = 1;

constexpr uint8_t SOUTH = 2;

constexpr uint8_t WEST = 3;

constexpr uint8_t EAST = 4;

uint8_t xy1_border = NONE;

if (x1 <= _xmin) {

xy1_border = WEST;

} else if (x1 >= _xmax) {

xy1_border = EAST;

} else if (y1 <= _ymin) {

xy1_border = SOUTH;

} else if (y1 >= _ymax) {

xy1_border = NORTH;

}

if (xy1_border != NONE) {

uint8_t xy2_border = NONE;

if (x2 <= _xmin) {

xy2_border = WEST;

} else if (x2 >= _xmax) {

xy2_border = EAST;

} else if (y2 <= _ymin) {

xy2_border = SOUTH;

} else if (y2 >= _ymax) {

xy2_border = NORTH;

}

if (xy2_border != NONE && xy2_border != xy1_border) {

return true;

}

}

return false;

};

auto [lower_levels, upper_levels] = get_lowers_and_uppers();

std::stringstream ss;

if (filled_) {

// Plot the line segments

// Create one filled curve for each segment of a contour.

bool plot_background = _levels[0] < zmin_;

if (plot_background) {

// find background polygon level

auto &largest_segment_with_children = line_segments_[0];

auto &largest_segment =

std::get<0>(largest_segment_with_children);

size_t line_index = std::get<0>(largest_segment);

size_t level_index = line_index > 0 ? line_index - 1 : 0;

size_t segment_begin = std::get<1>(largest_segment);

size_t segment_end = std::get<2>(largest_segment);

bool parent_is_lower_level =

is_lower_level(line_index, segment_begin, segment_end);

double background_z_level = parent_is_lower_level

? lower_levels[level_index]

: upper_levels[level_index];

// Plot background polygon

ss << " " << _xmin << " " << _ymin << " "

<< background_z_level << "\n";

ss << " " << _xmin << " " << _ymax << " "

<< background_z_level << "\n";

ss << " " << _xmax << " " << _ymax << " "

<< background_z_level << "\n";

ss << " " << _xmax << " " << _ymin << " "

<< background_z_level << "\n";

ss << " " << _xmin << " " << _ymin << " "

<< background_z_level << "\n";

ss << " e\n";

}

for (size_t i = 0; i < line_segments_.size(); ++i) {

// Send data for parent polygon

auto &parent_and_children = line_segments_[i];

auto &parent_segment = std::get<0>(parent_and_children);

size_t line_index = std::get<0>(parent_segment);

size_t begin_index = std::get<1>(parent_segment);

size_t end_index = std::get<2>(parent_segment);

for (size_t j = begin_index; j < end_index; ++j) {

double x = filled_lines_[line_index].first[j];

double y = filled_lines_[line_index].second[j];

// z = palette value

bool is_ll =

is_lower_level(line_index, begin_index, end_index);

double segment_z_level = is_ll ? lower_levels[line_index]

: upper_levels[line_index];

ss << " " << x << " " << y << " " << segment_z_level

<< "\n";

// work-around for edge cases

bool is_one_before_last = j == end_index - 2;

if (is_one_before_last) {

double next_x = filled_lines_[line_index].first[j + 1];

double next_y = filled_lines_[line_index].second[j + 1];

bool is_jump_to_border =

(next_x <= _xmin) || (next_x >= _xmax) ||

(next_y <= _ymin) || (next_y >= _ymax);

if (is_jump_to_border) {

// Parents are supposed to go anticlockwise

// Last move cannot be to the left

bool going_anticlockwise = next_x < x;

if (going_anticlockwise) {

// go to closest border and fill the path

std::array<double, 4> border_distances = {

(_xmax - x), (x - _xmin), (_ymax - y),

(y - _ymin)};

auto it =

std::min_element(border_distances.begin(),

border_distances.end());

size_t index_closest =

it - border_distances.begin();

if (index_closest == 0) {

x = _xmax;

} else if (index_closest == 1) {

x = _xmin;

} else if (index_closest == 2) {

y = _ymax;

} else if (index_closest == 3) {

y = _ymin;

}

ss << " " << x << " " << y << " "

<< segment_z_level << "\n";

auto [xs, ys] = fill_border_jump(

x, y, next_x, next_y, _xmin, _xmax, _ymin,

_ymax, true);

for (size_t k = 0; k < xs.size(); ++k) {

ss << " " << xs[k] << " " << ys[k]

<< " " << segment_z_level << "\n";

}

}

}

}

}

ss << " e\n";

// Send data for child/hole polygons

auto &child_segments = std::get<1>(parent_and_children);

for (size_t j = 0; j < child_segments.size(); ++j) {

auto &child_segment = child_segments[j];

size_t child_line_index = std::get<0>(child_segment);

size_t child_begin_index = std::get<1>(child_segment);

size_t child_end_index = std::get<2>(child_segment);

for (size_t k = child_begin_index; k < child_end_index;

++k) {

double x = filled_lines_[child_line_index].first[k];

double y = filled_lines_[child_line_index].second[k];

bool is_ll =

is_lower_level(child_line_index, child_begin_index,

child_end_index);

double segment_z_level =

is_ll ? lower_levels[child_line_index]

: upper_levels[child_line_index];

ss << " " << x << " " << y << " "

<< segment_z_level << "\n";

if (k != child_end_index - 1) {

double next_x =

filled_lines_[child_line_index].first[k + 1];

double next_y =

filled_lines_[child_line_index].second[k + 1];

if (is_border_jump(x, y, next_x, next_y)) {

auto [xs, ys] = fill_border_jump(

x, y, next_x, next_y, _xmin, _xmax, _ymin,

_ymax, true);

for (size_t l = 0; l < xs.size(); ++l) {

ss << " " << xs[l] << " " << ys[l]

<< " " << segment_z_level << "\n";

}

}

}

}

ss << " e\n";

}

}

// check if there are nans to hide

bool z_has_nans = false;

for (size_t i = 0; !z_has_nans && i < Z_data_.size(); ++i) {

for (size_t j = 0; !z_has_nans && j < Z_data_[i].size(); ++j) {

if (!std::isfinite(Z_data_[i][j])) {

z_has_nans = true;

}

}

}

if (z_has_nans) {

std::vector<size_t> nan_columns;

std::vector<size_t> nan_lines;

// for each line

for (size_t i = 0; i < Z_data_.size(); ++i) {

bool all_nan = true;

for (size_t j = 0; j < Z_data_[i].size(); ++j) {

if (std::isfinite(Z_data_[i][j])) {

all_nan = false;

break;

}

}

if (all_nan) {

nan_lines.emplace_back(i);

}

}

// for each column

for (size_t i = 0; i < Z_data_[0].size(); ++i) {

bool all_nan = true;

for (size_t j = 0; j < Z_data_.size(); ++j) {

if (std::isfinite(Z_data_[j][i])) {

all_nan = false;

break;

}

}

if (all_nan) {

nan_columns.emplace_back(i);

}

}

// draw polygons hiding these nans

for (size_t i = 0; i < nan_lines.size(); ++i) {

size_t first_line = nan_lines[i];

size_t end_line = nan_lines[i];

// find contiguous elements

for (size_t j = i + 1; j < nan_lines.size(); ++j) {

if (nan_lines[j] == end_line + 1) {

end_line = j;

} else {

break;

}

}

// plot a square hiding lines [first_line, end_line]

double ybegin = Y_data_[first_line][0];

double yend = end_line < Y_data_.size() - 1

? Y_data_[end_line + 1][0]

: _ymax;

ss << " " << _xmin << " " << ybegin << "\n";

ss << " " << _xmin << " " << yend << "\n";

ss << " " << _xmax << " " << yend << "\n";

ss << " " << _xmax << " " << ybegin << "\n";

ss << " " << _xmin << " " << ybegin << "\n";

// jump the lines we are plotting already

i = end_line;

}

for (size_t i = 0; i < nan_columns.size(); ++i) {

size_t first_col = nan_columns[i];

size_t end_col = nan_columns[i];

// find contiguous elements

for (size_t j = i + 1; j < nan_columns.size(); ++j) {

if (nan_columns[j] == end_col + 1) {

end_col = j;

} else {

break;

}

}

// plot a square hiding cols [first_col, end_col]

double xbegin = X_data_[0][first_col];

double xend = end_col < X_data_[0].size() - 1

? X_data_[0][end_col + 1]

: _xmax;

ss << " " << xbegin << " " << _ymin << "\n";

ss << " " << xbegin << " " << _ymax << "\n";

ss << " " << xend << " " << _ymax << "\n";

ss << " " << xend << " " << _ymin << "\n";

ss << " " << xbegin << " " << _ymin << "\n";

// jump the lines we are plotting already

i = end_col;

}

ss << " e\n";

}

}

// Plot the lines

// For the lines, we don't need to plot the segments separately

// One line is separated from the other with NaNs

auto is_separator = [](double x, double y) {

return !std::isfinite(x) || !std::isfinite(y);