Help with example CGAL 6.0 - 2D and Surface Function Interpolation

[DarkSun812]

Good morning everyone, today I tried to adapt a code that is in the CGAL examples (https://doc.cgal.org/latest/Interpolation/Interpolation_2interpolation_2_8cpp-example.html) but when trying to make the first interpolation, it does not print the result. I will leave the adaptation I made to the code to see if you can help me to detect the error I'm making.

// Compares the result of several interpolation methods

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Random.h>
#include <CGAL/squared_distance_2.h>

#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Interpolation_traits_2.h>
#include <CGAL/natural_neighbor_coordinates_2.h>
#include <CGAL/interpolation_functions.h>

#include <CGAL/point_generators_2.h>
#include <CGAL/algorithm.h>
#include <CGAL/Origin.h>

#include <cassert>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef K::FT                                         Coord_type;
typedef K::Vector_2                                   Vector;
typedef K::Point_2                                    Point;

typedef CGAL::Delaunay_triangulation_2<K>             Delaunay_triangulation;
typedef CGAL::Interpolation_traits_2<K>               Traits;

typedef std::vector< std::pair<Point, Coord_type> >   Coordinate_vector;
typedef std::map<Point, Coord_type, K::Less_xy_2>     Point_value_map;
typedef std::map<Point, Vector, K::Less_xy_2 >       Point_vector_map;


void readExcelCSV(std::vector<Point>& points, std::vector<Point>& pointsInterpolation, Point_value_map& values, Delaunay_triangulation& T) {
    std::ifstream file("datos.csv");
    std::ifstream filePoints("datosforInterpolation.csv");

    if (!file || !filePoints) {
        std::cerr << "Error in Open File." << std::endl;
        return;
    }


    std::string line;

    int cont = 0;
    while (std::getline(file, line)) {
        std::stringstream ss(line);
        std::string value;
        std::vector<double> row;
        
        for (int idx = 0; idx < 3; idx++)
        {
            std::getline(ss, value, ',');
            row.push_back(std::stod(value));
        }
       
        points.push_back(Point(row[0], row[1]));
        T.insert(points[cont]);
        Coord_type value2 = row[2];
        values.insert(std::make_pair(points[cont], value2));
        cont++;
    }
    cont = 0;
    while (std::getline(filePoints, line)) {
        std::stringstream ss(line);
        std::string value;
        std::vector<double> row;

        for (int idx = 0; idx < 3; idx++)
        {
            std::getline(ss, value, ',');
            row.push_back(std::stod(value));
        }

        pointsInterpolation.push_back(Point(row[0], row[1]));
        cont++;
    }

    filePoints.close();
    file.close();
    
};

int main()
{
    std::vector<Point> points;
    std::vector<Point> pointsInterpolation;
    Delaunay_triangulation T;

    Point_value_map values;
    Point_vector_map gradients;
    
    std::cout << "Go to read File csv..." << std::endl;
    readExcelCSV(points, pointsInterpolation,values, T);
    std::cout << "Read Excel and build delauny triangulation" << std::endl;
    
    //variables for statistics:
    std::pair<Coord_type, bool> res;
    Coord_type  l_total = Coord_type(0),
       l_value(l_total);
   

    //interpolation + error statistics
    for (int i = 0; i < pointsInterpolation.size(); i++) {
        std::cout << "point " << i << std::endl;
        Coord_type x(pointsInterpolation[i].x());
        Coord_type y(pointsInterpolation[i].y());

        
        std::cout << pointsInterpolation[i].x() << " " << pointsInterpolation[i].y()<< std::endl;
        //Coordinate_vector:
        std::vector< std::pair< Point, Coord_type > > coords;
        Coord_type norm =
            CGAL::natural_neighbor_coordinates_2(T, pointsInterpolation[i],
                std::back_inserter(coords)).second;

        assert(norm > 0);
        
        //linear interpolant:
        l_value = CGAL::linear_interpolation(coords.begin(), coords.end(),
            norm,
            CGAL::Data_access<Point_value_map>(values));
       
        std::cout << pointsInterpolation[i].x() << " " << pointsInterpolation[i].y() << " " << l_value << std::endl;
    }
    
    return EXIT_SUCCESS;
}

Here is an image of what the delauny triangle should look like and the points I want to find.

the csv file I want to enter

X	Y	Z
0	485	169
1	500	168
3	499	170
5	501	169
2	470	167.8
4	470	168.5
6	470	167.7

The csv file with interpolation points

X	Y	Z
0.5	490	?
2.5	480	?
4	489	?
5	489	?

I have no idea where I'm wrong, if you could help me to make this example work I would be very grateful, also if you could explain me what these lines are for since I adapted the code by logic since I saw how it worked, but I don't really know what I did.

[DarkSun812]

This is the failing line.

l_value = CGAL::linear_interpolation(coords.begin(), coords.end(),
    norm,
    CGAL::Data_access<Point_value_map>(values));

takes about 2 seconds and then the program closes.

[DarkSun812]

I know I just asked the question but I saw what the error was, the problem is that in the dataset I was using, the first value to interpolate did not return coordinates.

The following code before the failed line prevents the program from closing and continues with the other points.

if (coords.empty()) {
    std::cout << pointsInterpolation[i].x() << " " << pointsInterpolation[i].y() <<" Not valid coords" << std::endl;
    continue;
}

But this opens another question, how can I do to calculate those missing points, is there any way known to CGAL to handle those cases?

What do you recommend for these cases?