Problem of type

[pierrotdu18]

Hi,
I've tried to use templates and smart pointers but it's like I don't even understand my own code, I'm completely lost and so sad because it's been now almost 1 week that my project hasn't compiled... (I had some nice 3D reconstructions before ^^)

I hope this is not irrespectful to give you my code like that but I hope you'll understand (I've discovered POO last week ^^)

I hope, as well, that my code is not ridiculous and that there are still some things that I've understood...

#ifndef INCLUDE_GEOMETRY
#define INCLUDE_GEOMETRY

#include <vector>
#include <set>
#include <limits>
#include <cmath>
#include <memory>

struct Coordinates;
struct Point;
struct Edge;
struct Triangle;

double distance(const Coordinates&, const Coordinates&);
double distance(const std::shared_ptr<Point>&, const std::shared_ptr<Point>&);
double distance(const std::shared_ptr<Point>&, const std::vector<std::shared_ptr<Point>>&);
shared_ptr<Point> farthest(const std::vector<Point>&, const std::vector<std::shared_ptr<Point>>&);
Edge get_edge(const std::shared_ptr<Point>&, const std::shared_ptr<Point>&);
Triangle get_triangle(const std::shared_ptr<Point>&, const std::shared_ptr<Point>&, const std::shared_ptr<Point>&);

struct Coordinates {
    double x;
    double y;
    double z;
};

struct SimplicialComplex;

struct Point {

};

template<double x, double y, double z>
struct TemplatedPoint : public Point {
    Coordinates coordinates {x, y, z};
    size_t index;
    TemplatedPoint(size_t index) {this->index = index; simplices = std::make_shared<SimplicialComplex>();}

    template<double ox, double oy, double oz>
    struct Compare {
        Coordinates origin {ox, oy, oz};
        bool operator()(TemplatedPoint* lhs, TemplatedPoint* rhs) {
            return distance(lhs->coordinates, this->coordinates) < distance(rhs->coordinates, this->coordinates);
        }
    };

    bool operator<(const TemplatedPoint& rhs) const {return this->index < rhs.index;}
    bool operator==(const TemplatedPoint& rhs) const {return this->index == rhs.index;}

    std::multiset<std::shared_ptr<TemplatedPoint>, Compare<x, y, z> > neighbourhood;

    std::shared_ptr<SimplicialComplex> simplices;
};

struct Edge {
    Edge(std::shared_ptr<Point>, std::shared_ptr<Point>);
    std::shared_ptr<Point> p1;
    std::shared_ptr<Point> p2;
};

struct Triangle {
    Triangle(std::shared_ptr<Point>, std::shared_ptr<Point>, std::shared_ptr<Point>);
    std::shared_ptr<Point> p1;
    std::shared_ptr<Point> p2;
    std::shared_ptr<Point> p3;
};

struct Edge_less {
    bool operator()(const Edge& e1, const Edge& e2) {
        return *e1.p1 < *e2.p1 || (*e1.p1 == *e2.p1 && *e1.p2 < *e2.p2);
    }
};

struct Triangle_less {
    bool operator()(Triangle t1, Triangle t2) {
        return *t1.p1 < *t2.p1 || (*t1.p1 == *t2.p1 && *t1.p2 < *t2.p2) || (*t1.p1 == *t2.p1 && *t1.p2 == *t2.p2 && *t1.p3 < *t2.p3);
    }
};

struct SimplicialComplex {
    std::set<Edge, Edge_less> edges;
    std::set<Triangle, Triangle_less> triangles;
};

#endif // INCLUDE_GEOMETRY 

I'm trying to initialise W like this :

vector<shared_ptr<Point> > W;
size_t i = 0;
while (!ifile.eof()) {
    double x, y, z;
    ifile >> x;
    ifile >> y;
    ifile >> z;
    auto w = make_shared(TemplatedPoint<x, y, z>(i));
    W.push_back(w);
    ++i;
}

[pierrotdu18]

Well I'm going to formulate some precise questions :

In the code above, TemplatedPoint is a kind of Point, but even when I use a std::vector<std::shared_ptr<Point>> and when I put some pointers to objets of type TemplatedPoint in the vector, then I can't use their attributes. I think this is because of the vector converts my ptrs into Point ptrs. But as I use those elements very often in the code, I can't convert it with static_cast all the time (and I don't even know if it would work).
So, how could I do ?

Moreover you talked about double& for the non-type template, I'm very interrested because my points have plenty of decimals and it would be very annoying to make a tool to convert double to int, int to double (fixed points etc.)

I hope you'll have time to answer me.
Joyeux noël à vous aussi

[gunnerfunner]

For incomplete types such as SimplicalComplex, use (smart) pointers

... and you've gone and changed everything to smart pointers, making life unnecessarily difficult for yourself!

If you'd overloaded the < operator for Edge and Triangle as I'd suggested in my last post yesterday with a general example then you'd not have needed Edge_less and Triangle_less

You don't need Point and TemplatePoint, just Point will suffice and see it's implementation below

my points have plenty of decimals ...

that's perfectly fine, both Coordinates and Point objects have double data members, the only restriction you have is that the 'Origin' from which all distances are measured have integer co-ordinates

The following program (a) compiles successfully, (b) reads the file accurately into Coordinates objects, (c) shows how you can convert Coordinates objects into Point objects, (d) or how you can set up the neighborhood multi-set around any Origin and (e) illustrates how to access the SimplicalComplex members of any Coordinates object

Please read everything here carefully first and foremost. I think you're very close to your goal, just don't panic or lose concentration. Good luck!

#include <iostream>
#include <memory>
#include <set>
#include <cmath>
#include <fstream>
#include <algorithm>

struct Coordinates
{
    double _x;
    double _y;
    double _z;

    Coordinates() {}
    Coordinates (double x, double y, double z) : _x(x), _y(y), _z(z){};
};
std::istream& operator >> (std::istream& is, Coordinates& rhs)
{
    is >> rhs._x >> rhs._y >> rhs._z;
    return is;
}
std::ostream& operator << (std::ostream& os, const Coordinates& p)
{
    os << "(" << p._x << ", " << p._y << ", " << p._z << ")";
    return os;
}
double distance(const Coordinates& lhs, const Coordinates& rhs);

template <typename T, int a, int b, int c>
struct SimplicalComplex;

template <typename T, int a, int b, int c>
struct Point
{
    double _P_x;
    double _P_y;
    double _P_z;
    std::unique_ptr<SimplicalComplex<T, a, b, c>> _simplices;

    Point(){};
    Point (double x, double y, double z) : _P_x(x), _P_y(y), _P_z(z){}
    template <typename T1, int a1, int b1, int c1>
    struct DistanceCompare
    {
        T _origin = T(a1, b1, c1);

        bool operator()(const Coordinates& lhs, const Coordinates& rhs)const
        {
            return distance(lhs, _origin) < distance (rhs, _origin);
        }
    };
    std::multiset<Coordinates, DistanceCompare<T, a, b, c>> neighborhood;

    bool operator < (const Point<T, a, b, c>& rhs)const
    {
        if(_P_x < rhs._P_x){return *this;}
        else if ((_P_x == rhs._P_x) && (_P_y < rhs._P_y)){return *this;}
        else if ((_P_x == rhs._P_x) && (_P_y == rhs._P_y) && (_P_z < rhs._P_y)){return *this;}
    }
    bool operator == (const Point<T, a, b, c>& rhs)const
    {
        if((_P_x == rhs._P_x) && (_P_y == rhs._P_y) && (_P_z == rhs._P_y)) {return true;}
        else {return false;}
    }
    void fromCoordinates (const Coordinates& rhs)
    {
        _P_x = rhs._x; _P_y = rhs._y; _P_z = rhs._z;
    }
};
template <typename T, int a, int b, int c>
std::ostream& operator << (std::ostream& os, const Point<T, a, b, c>& rhs)
{
    os << "(" << rhs._P_x << ", " << rhs._P_y << ", " << rhs._P_z << ")";
    return os;
}
template <typename T, int a, int b, int c>
struct Edge
{
    Point<T, a, b, c> _less;
    Point<T, a, b, c> _more;

    bool operator < (const Edge<T, a, b, c>& rhs) const
    {
        if(_less < rhs._less){return *this;}
        else if((_less == rhs._less) && (_more < rhs._more)){return * this;}
    }
};
template <typename T, int a, int b, int c>
struct Triangle
{
    Point <T, a, b, c> _low;
    Point <T, a, b, c> _mid;
    Point <T, a, b, c> _high;

     bool operator < (const Triangle<T, a, b, c>& rhs) const
    {
        if(_low < rhs._low){return *this;}
        else if((_low == rhs._low) && (_mid < rhs._mid)){return * this;}
        else if ((_low == rhs._low) && (_mid == rhs._mid) && (_high < rhs._high)){return *this;}
    }
};
template <typename T, int a, int b, int c>
struct SimplicalComplex
{
    std::set<Edge<T, a, b, c>> _edges;
    std::set<Triangle<T, a, b, c>> _triangles;
};

int main()
{
    using Origin = Point<Coordinates, 0, 0, 0>;//to save some typing
    Origin p(0, 0, 0);
    Coordinates p1(0, 0, 0);
    std::vector<Coordinates> vec;
    std::ifstream infile("D:\\input.txt");
    while(infile)
    {
        Coordinates temp;
        infile >> temp;//read the file into Coordinates, can convert to type Origin later if reqd (see below)
        if(infile)
        {
            vec.push_back(std::move(temp));
        }
    }
    Origin temp;
    std::vector<Origin> neighbors;//you can convert the Coordinates to Point<Coordinates, 0, 0, 0> if you wish ...
    for (auto& elem: vec)
    {
        temp.fromCoordinates(elem);
        neighbors.push_back(std::move(temp));
    }
    for (auto& elem: vec)
    {
        p.neighborhood.insert(elem);// ... or insert the Coordinates directly into the neighborhood of p
    }

    for (auto& elem: p.neighborhood)
    {
        std::cout << elem << '\n';
    }
    //to access the SimplicalComplex of p, use:
    //p._simplices->SimplicalComplex<Point, 0, 0, 0>::_edges //
     //p._simplices->SimplicalComplex<Point, 0, 0, 0>::_triangles //
}
double distance(const Coordinates& lhs, const Coordinates& rhs)
{
    return std::pow(lhs._x - rhs._x, 2) + std::pow(lhs._y - rhs._y, 2) + std::pow(lhs._z - rhs._z, 2);
}

//tested with the following data file: http://pastie.org/10985735 

ps: why looping on eof() is not good:http://stackoverflow.com/questions/5605125/why-is-iostreameof-inside-a-loop-condition-considered-wrong

[pierrotdu18]

First I want to thank you a lot for your reply.

Then, I'm fearing that there are some things that I didn't explain very clearly, and that you understood wrong (I hope I'm wrong and I juste have to read the code more carefuly, but to be sure....).

The fact that a point has an origin is not essential, it is even useless. At first I put a "mayor" but it was certainly just a bad idea for making possible to order the neighbours by distance to the point.
Let's recap : a point has just (for real) those information to carry :
- its coordinates (x, y, z) which are of type double
- references or pointers to its neighbours, which are in a std::multiset and must be order by distance from the point itself (there is no orgin but the point itself actually)
- some simplicial (which are "witnessed" by the point, see the article if you want to learn more), so some edges and triangles

When I'm reading your code I feel like there's a misunderstanding somewhere (I hope there isn't any)...
Indeed, I think that despite all that I havn't really understood why I need templates, as the parameters a, b, c are actually exactly the coordinates of the point...

Also, another thing, I'd like to have the minimum number of points to allocate, so I'd like to initialize W and then use references or pointers to those points (for the simplexes for instance, or the neighbours). It is important that when a function modifies a point, the whole program is aware of the modification

Thanks for the link !

[gunnerfunner]

... the parameters a, b, c are actually exactly the coordinates of the point...

This actually makes things much simpler and you can work with doubles throughout, no need for any integers or templates (as you rightly guessed). The correct choice of container in this case would be std::multimap (not map for the same reason as not set earlier). I'll show you a stylized example of how to set up the neighborhood for any Point:

#include <iostream>
#include <map>
#include <cmath>
#include <iomanip>

struct Point
{
    double _x;
    double _y;
    double _z;

    Point() {}
    Point (double x, double y, double z) : _x(x), _y(y), _z(z){};
    std::multimap<double, Point> neighborhood;
};
std::ostream& operator << (std::ostream& os, const Point& p)
{
    os << "(" << p._x << ", " << p._y << ", " << p._z << ")";
}
double distance(const Point& lhs, const Point& rhs);
int main()
{
    Point origin(0, 0, 0);
    Point a(11, 21, -3), b(12, 5, 6), c(-12, 5, 6), d(2, 3, 4), e(13, 6, 1), f(1, 2, 2);

    origin.neighborhood.insert(std::pair<double, Point>(distance(origin, a), a));
    origin.neighborhood.insert(std::pair<double, Point>(distance(origin, b), b));
    origin.neighborhood.insert(std::pair<double, Point>(distance(origin, c), c));
    origin.neighborhood.insert(std::pair<double, Point>(distance(origin, d), e));
    origin.neighborhood.insert(std::pair<double, Point>(distance(origin, e), e));
    origin.neighborhood.insert(std::pair<double, Point>(distance(origin, f), f));

    std::cout << std::setprecision(4) << std::left ;
    for (auto& elem : origin.neighborhood)
    {
        std::cout << elem.first << "\t" << elem.second << '\n';
    }

}
double distance(const Point& lhs, const Point& rhs)
{
    return std::pow((std::pow(lhs._x - rhs._x, 2) + std::pow(lhs._y - rhs._y, 2) + std::pow(lhs._z - rhs._z, 2)), 0.5);
}

Output

3       (1, 2, 2)
5.385   (13, 6, 1)
14.32   (12, 5, 6)
14.32   (-12, 5, 6)
14.35   (13, 6, 1)
23.9    (11, 21, -3)

Process returned 0 (0x0)   execution time : 0.109 s
Press any key to continue.

It is important that when a function modifies a point, the whole program is aware of the modification

That would automatically happen if you delete the original element and insert the new (modified) element into the neighborhood, however it would not be possible to edit the value type (here Point) stored in the multiset directly because changing the Point would also change the key value which is used to store the elements in the container

[pierrotdu18]

Wow thank you so much, I knew there was a misunderstanding.
Now I have a last question : I don't really see how the second part of your last reply answers my question, I must have not understood everything of course.
Let's consider a point p.
Let's say that there are some other points : a, b, c, d...
If a and b belong to p's neighbours, is it possible to do something like that so the points are not copied but without using pointers ?
p.neighbourhood.insert(reference of b)

It is important that every information a point carries about other points should be references or pointers (a neighbour is somewhere a volatile pointer to another point which is a neighbour of p in the "real world")

[pierrotdu18]

I'm so relieved, my code is compiling.
However, the execution of the program is really quickly interrupted...
I'm looking for a reason why this should be wrong, so far without success.
I've updated my code here : https://github.com/pierrotdu18/TIPE/blob/master/Reconstruction/Algorithmes/3D/geometry.h

If you want I've also updated every other file in the parent directory.

Well, with all the help you've given me I will cite you in the "thanks section" of my written report !
(and send you this report if you'd like to)

[gunnerfunner]

There are some other things I've spotted that might be troubling your program:

- there is no need to overload distance(), the second distance function seems to be calculating the nearest Point and so it should be called nearest (just as you have a farthest). Then consider this statement from .cpp, line 14:
double dl = distance(p, l);
p and l are pointer variables and so, without dereferencing, you are calculating the 'distance' between their memory addresses. Same problem with the Edge ctor (.cpp, line 65) - just comparing memory addresses w/o de-referencing

It'd also be better to define the distance with Coordinates objects to be consistent with actual geometry, in which case the 'distance' b/w 2 const Point* objects, p1 and p2, would simply be:
distance(p1->coordinates, p2->coordinates)

Now what is going on with the farthest() function, .cpp line 26:
double d = distance(p, L); where L is a vector of pointers to Point??

In case you're trying to get the nearest and farthest Point objects from any given Point object this info is readily available in the neighborhood multimap of the Point object, you just have to dereference the begin() and end() iterators and pick off the second element of the <double, Point> pair, there is no need for additional functions

The Triangle ctor (.cpp, lines 35 - 62) can be simplified with a custom sort function:

struct PointerSort
{
    bool operator < (Point* lhs, Point* rhs)
    {
        return lhs ->index < rhs -> index;
    }
};

Triangle::Triangle(Point* i, Point* j,  Point* k)
{
    std::vector<Point*> ptr_vec;
    ptr_vec.push_back(i);
    ptr_vec.push_back(j);
    ptr_vec.push_back(k);
    std::sort(ptr_vec.begin(), ptr_vec.end(), PointerSort());
    //then you can assign the pointers as you wish
}

operator < overload for Edge, Triangle - check the if, else-if syntax I sent you yesterday, I am not sure your current code is working without any conditional statements. Also, since you're using pointers, please make double-sure that you're dereferencing them accurately throughout unlike in the case of distance() and Edge above. I'm not checking all the formulae as I think you know sufficient C++ to do this yourself

I cannot find the reconstruction() function (main(), line 52) declared/defined anywhere;

main(), line 42: still looping on eof() despite @cire's warnings and my alternative code and link why such loop is bad

main(), line 48: W[i]->simplices = make_unique<SimplicialComplex>();
I can see what you're trying to do here, viz. assign to the simplices member of each Point object pointed by each cell of W but think of the data-type of simplices which is SimplicalComplex which in turn has a set<Edge>, edges, and set<Triangle> triangles. So it seems that each Point object has its corresponding set<Edge> and set<Triangle> and you need to define these first, just as we did with the neighborhood, and then make the simplices member

Let's say that there are some other points : a, b, c, d...
If a and b belong to p's neighbours, is it possible to do something like that so the points are not copied

Yes, you can use move or forwarding semantics to do this while a, b are being moved to p's neighborhood. Referring to my previous post, use this instead:

//#include<utility> for both
origin.neighborhood.insert(std::make_pair(distance(origin, a), std::move(a)));
//or
origin.neighborhood.insert(std::make_pair(distance(origin, a), std::forward<Point>(a)));

a is passed to distance by reference and then moved or forwarded to make_pair and inserted into origin's neighborhood (you can look up the difference online). However be aware that after this a is left in a valid but unspecified state, which means that it is possible to destroy and to assign to the variable, but nothing else

[pierrotdu18]

For the distance everything is normal :
distance(Point*, Point*) : I dereference the points in the function and evaluate with the coordinates
distance(Point* vector<Point*>) : distance from between the point and the set of points (in maths we use the word "distance" for this function that's why I also used this word)
farthest(W, L) : returns the farthest point of W\L from L

I've done every other thing you've proposed :)

A question : when I do this :

vector<Point> W;

    size_t i = 0;
    while (ifile) {
        double x, y, z;
        ifile >> x;
        ifile >> y;
        ifile >> z;
        Point w(i, x, y, z);
        W.push_back(move(w));
        ++i;
    }

every point of W has the same adress ! How to avoid that ?

[keskiverto]

every point of W has the same adress

How do you check that and does it matter?


PS. I do have some trouble to comprehend.

You do have 3D-points as input. Each point should be unique.
Somehow, you do determine edges that exist between some points.
Somehow, you do determine triangles that are defined by some points.
Either necessary for the above, or totally unrelated, is the need to find points that are near each other.
Perhaps you want to draw images?

Am I close?

[pierrotdu18]

To check that I make another vector which contains every address of the other vector's elements and I print them.
It does matter because I need to have a vector containing all the points of the input file, not just one point!

You're very close : from an input cloud of points, which are supposed to be sampled over an underlying shape, I try to reconstruct the shape with a good extraction of a good triangulation (and then I print my reconstruction in 3D with the software geomview)

You can have a look to the research article I linked in a previous post

[keskiverto]

It seems that the input data is essentially:
std::vector<Coordinates> W; // "witnesses"
Furthermore, is the W const within the algorithm?

The L ("landmarks") seems to be some subset of W.
Are the xyz-coordinates of any point changed at any time?
Do you create any new points, rather than always use members of W?

If a point from W is added to L, will it be considered any more on further iterations that add to L?


IMHO the lists of triangles, edges and whatnot are not a property of the "witnesses". They are separate concepts.

To check that I make another vector which contains every address of the other vector's elements and I print them.

Can you show that code?

[gunnerfunner]

Github still shows your old code but, anyways, I have some general queries:

farthest(W, L) : returns the farthest point of W\L from L

L is vector<Point>, so what does it mean to say 'farthest point ... from L'? Also, W used to be vector<Point>* previously I recall, have you changed all pointers to objects?

The other thing that I am wondering about is how Edge and Triangle are constructed? Are p1, p2 (for Edge) and p1, p2, p3 (for Triangle) coming from the Coordinates objects read from the file or different?. If they are coming from the file on what basis are the ctor arguments chosen?

The reason I'm asking is that it seems you have a circularity currently: Point will not be a complete type until we have information about Point::simplices which will need information about SimplicalComplex::edges and SimplicalComplex::triangles of which the former will need information about Edge::p1, Edge::p2 (which are of type Point*) and the latter information about Triangle::p1, Triangle::p2, Triangle::p3 (which are also type Point*). I am not sure even using pointers removes this circularity, it might have if SimplicalComplex did not depend on Point*.

So are you sure you need Point* as ctor parameters for Edge and Triangle or can they be constructed with Coordinates*? Of course that would also imply that the < operator overloads for these two classes have to be changed

[pierrotdu18]

@keskiverto

What the article doesn't mention very explicitly is how the witness complex is really updated, and what data structure is advised for the witness complex.

I have a class Point which carries some information : the index of the point (to make other things easier but it doesn't have any mathematical meaning), it's neighbours (pointers to other points, there are at most max(\nu_i) neighbours to take into account), and the simplexes which are witnessed by this point.
I chose so add these to the class Point because on one hand, it makes other things in the code easier, and on the other hand, it is my way of describing the witness complex (it makes things very nice to store it like this for its update).

I have the vector W, which contain points of type Point and not only coordinates. Every other point (L, neighbours, the simplexes) are actually pointers to some points of W. I don't create any point either modify x, y, z. The only things I do are : add a Point pointer to L at each iteration of the main loop, modify the each point's neighbours and at the same time update the list of simplices which are witnessed by this point.

I also use some other vectors but they only contain pointers to points of W.

My code was for (auto w : W) ofile << &w << endl;

[keskiverto]

L is vector<Point>, so what does it mean to say 'farthest point ... from L'?

Based on Github implementation file:

We have two sets of points, U and L.
For each member u of U we compute distance( u, L ).
That essentially creates a set D of distances:
Di = distance( Ui, L )

max_element( D ) leads to Ui -- the "farthest".


What was the distance( u, L )? Nearest:
Di = distance( u, Li )
return *min_element( D )

[gunnerfunner]

So it is the sum of the distances of each element of W (not U) from each element in L which gives some measure of aggregate distance from L and then take the maximum? The code is not written like that but at least conceptually it makes sense

[pierrotdu18]

I've updated me GitHub : https://github.com/pierrotdu18/TIPE/tree/master/Reconstruction/Algorithmes/3D

@gunnerfunner :
distance from x to y : d(x, y)
distance from a point x to a finite set of points d(x, L) : minimum of d(x, y) where y discribes the set of point

Those two definitions are very common in maths. This one is not, but I need it for my algorithm :
farthest(W, L) : x in W such as d(x, L) >= d(y, L) for every point y of W, in other words, x in W such that it maximizes the distance to L (sorry if this is not correct English)

To answer your last question : what could be possible is that Edge and Triangle contain not pointers to the points, but the points' indexes (size_t) : indeed it would simplify the problem

[keskiverto]

Not a sum, but lets play some:

double distance( const std::vector<std::vector<double>> & DM,
                 size_t row,
                 size_t col )
{
  if ( col < row ) std::swap( row, col );
  return DM[row][col];
}


double distance( const std::vector<std::vector<double>> & DM,
                 size_t u,
                 const std::vector<size_t> & L )
{
  double winner = ...
  for ( auto v : L ) {
    double d = distance( DM, u, v );
    if ( d < winner ) ...
  }
  return winner;
}


size_t farthest( const std::vector<std::vector<double>> & DM,
                 const std::vector<size_t> & U,
                 const std::vector<size_t> & L )
{
  double wd = -1.0;
  size_t winner = 0;
  for ( auto u : U ) {
    double d = distance( DM, u, L );
    if ( wd < d ) {
      wd = d;
      winner = u;
    }
  }
  return winner;
}


void foo( const std::vector<Coordinates> & W )
{
  const size_t N = W.size();

  // cache squared distances
  std::vector<std::vector<double>> DM( N, std::vector::<double>(N) );
  for ( size_t row = 0; row < N; ++row ) {
    const auto & R = W[row];
    for ( size_t col = row + 1; col < N; ++col ) {
      DM[row][col] = ( R - W[col] ).dot();
    }
  }

  std::vector<size_t> L;
  L.reserve( N );
  L.emplace_back( 0 ); // start with the first witness

  std::vector<size_t> U( N-1 );
  std::iota( U.begin(), U.end(), 1 ); // all, but 0

  // repeat the following as much as needed to create L // while( ! U.empty() )
  size_t u = farthest( DM, U, L );
  // move from U to L
  L.emplace_back( U[u] );
  U.erase( U.begin() + u ); // not efficient 

I might have some errors there. The main points are:

1. Coordinates are not copied, nor pointers. The referring is via indices.

2. The index of a "point" is implicit; its index within the vector W.

3. The distances (or actually just squared distances) are computed only once and the "farthest" uses those precomputed values.

[pierrotdu18]

Was this message intended for me ?

I've updated my GitHub with the new implementation of SimplicialComplex

But I do still have the problem of W which contains the same point 62312 times (it's the number of points of my test example)

EDIT:

That's ok for W I've solved the problem.
I've upated my code and there's still something wrong. It's ok for i = 1, 2 and for i = 3, it bugs when it reachs the x-th point where x is exactly 1883 (why??)

EDIT:

I've found the problem and the program is executed without error.
However the result is totally unsatisfying