//

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

//

#include <cmath>

#include <sstream>

#include <numeric>

#include <regex>

#include <unordered_set>

#include <matplot/axes_objects/contours.h>

#include <matplot/axes_objects/histogram.h>

#include <matplot/core/axes.h>

#include <matplot/util/common.h>

namespace matplot {

contours::contours(class axes* parent, const vector_2d& X, const vector_2d& Y, const vector_2d& Z, const std::string& line_spec)

: axes_object(parent), X_data_(X), Y_data_(Y), Z_data_(Z), line_spec_(this, line_spec) {

initialize_preprocessed_data();

contour_generator_ = QuadContourGenerator(X_data_,Y_data_,Z_data_,_corner_mask,nchunk_);

}

contours::contours(class axes* parent, const vector_2d& Z, const std::string& line_spec)

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

initialize_preprocessed_data();

contour_generator_ = QuadContourGenerator(X_data_,Y_data_,Z_data_,_corner_mask,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

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

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

}

// 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 (!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_.empty();

}

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::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::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(const std::string& 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 (!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 (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 (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 !isfinite(x) || !isfinite(y);

};

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

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

continue;

}

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

if (!is_separator(lines_[i].first[j],lines_[i].second[j])) {

double x = lines_[i].first[j];

double y = lines_[i].second[j];

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

if (!line_spec_.user_color()) {

if (!filled_ || colormap_line_when_filled_) {

// plot palette colors

ss << " " << levels_[i];

}

}

ss << "\n";

} else {

// Skip the next nans separating parents and children

// to avoid extra useless empty lines in case there is

// more than one nan

while (is_separator(lines_[i].first[j+1],lines_[i].second[j+1])) {

++j;

}

// Include an empty line to indicate this polygon or line segment is over

ss << " \n";

continue;

}

}

ss << "e\n";

}

if (contour_text_) {

// We always check if a label is not too close to another label

constexpr double minimum_distance = 0.8;

auto too_close = [minimum_distance](double x1, double y1, double x2, double y2) {

return abs(x1 - x2) < minimum_distance && abs(y1 - y2) < minimum_distance;

};

// Instead of iterating each line, we iterate

// each line a little bit. The rationale is that

// we might put some labels in line i and there

// is no room left for more labels in line i+1

// with the minimum distance.

// The current index in each line

std::vector<size_t> point_index(lines_.size(), 0);

std::vector<size_t> n_points(lines_.size());

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

n_points[i] = lines_[i].first.size();

}

// Label positions for each line

std::vector<std::vector<double>> labels_xs(lines_.size());

std::vector<std::vector<double>> labels_ys(lines_.size());

std::vector<std::vector<double>> labels_degrees(lines_.size());

// Vectors with line direction for each line

std::vector<double> direction_us(lines_.size(),0.0);

std::vector<double> direction_vs(lines_.size(),0.0);

// while we have not iterated all points in all lines

while (point_index != n_points) {

// iterate each line until we put a label on that line or the line is over

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

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

continue;

}

// references to line directions

double& direction_u = direction_us[i];

double& direction_v = direction_vs[i];

// for each point in that line (starting from where we stopped)

for (; point_index[i] < lines_[i].first.size(); ++point_index[i]) {

size_t j = point_index[i];

if (is_separator(lines_[i].first[j], lines_[i].second[j])) {

// reset line direction

direction_u = 0.0;

direction_v = 0.0;

continue;

}

// check if not too close to another label

const double x = lines_[i].first[j];

const double y = lines_[i].second[j];

bool too_close_to_another_label = false;

for (size_t k = 0; !too_close_to_another_label && k < labels_xs.size(); ++k) {

for (size_t l = 0; !too_close_to_another_label && l < labels_xs[k].size(); ++l) {

if (too_close(x, y, labels_xs[k][l], labels_ys[k][l])) {

too_close_to_another_label = true;

}

}

}

// Update vector of line direction

constexpr double moving_average_reduction_factor = 0.1;

direction_u *= (1-moving_average_reduction_factor);

direction_v *= (1-moving_average_reduction_factor);

bool has_previous = j != 0 && isfinite(lines_[i].first[j-1]);

if (has_previous) {

direction_u += lines_[i].first[j] - lines_[i].first[j-1];

direction_v += lines_[i].second[j] - lines_[i].second[j-1];

} else {

bool has_next = j != lines_[i].first.size() - 1 && isfinite(lines_[i].first[j+1]);

if (has_next) {

direction_u += lines_[i].first[j + 1] - lines_[i].first[j];

direction_v += lines_[i].second[j + 1] - lines_[i].second[j];

}

}