Efficiently find duplicated vectors in C++

I'm working on a R package dealing with 3D meshes. I have a function which constructs a parametric mesh but for certain parameterizations it can generate duplicated vertices (numeric vectors of length 3). I need to find the duplicated vertices in order to merge them into a single vertex. I implemented an algorithm with Rcpp to do that.

There's a function for that in the Rvcg package: vcgClean(mymesh, sel = 0). It is very fast, and as compared to it, my C++ algorithm is very slow. Similary, the R base function duplicated is highly faster.

Here is how my algorithm works. Say there are n vertices.

I take vertex 1, I compare it to vertices 2, 3, ..., n.
I take vertex 2, I compare it to vertices 3, 4, ..., n. Note that I could skip this comparison in case if this vertex has been marked as a duplicate of vertex 1 at the previous step. I don't do that. Anyway there's usually a small number of duplicates, so this is not the cause of the slowness.
And so on, vertex 3 is compared to vertices 4, 5, ..., n. Etc.

To compare two vertices, I test the near-equality of each coordinate:

bool test = nearEqual(v1(0), v2(0)) && nearEqual(v1(1), v2(1)) && nearEqual(v1(2), v2(2));

Before, I tested exact equality, this was slow as well. My nearEqual function is:

bool nearEqual(double l, double r) {
  return r == std::nextafter(l, r);
}

Why is my algorithm so slow? Is it due to the n-1 + n-2 + ... + 1 comparisons strategy? I don't see how I could do less comparisons. Or is it due to the nearEqual function? Or what? The vertices are stored in a 3 x n numeric matrix in Rcpp.

Edit

I've just tried another way to test near-equality of two columns but this is still slow:

const Rcpp::NumericMatrix::Column& v1 = Vertices.column(i);
const Rcpp::NumericMatrix::Column& v2 = Vertices.column(j);
bool test = Rcpp::max(Rcpp::abs(v1 - v2)) < 1e-16;

Solution

You could keep track of discovered vertex positions in a hash map and duplicated vertices in a set to only have to do one pass:

#include <Rcpp/Lightest>
#include <map>
#include <set>
#include <tuple>

typedef std::tuple<double, double, double> Point3;

// [[Rcpp::export()]]
Rcpp::List duplicated_vertices(Rcpp::NumericMatrix x) {
    std::map<Point3, Rcpp::IntegerVector> positions;
    std::set<Point3> duplicates;

    for (R_len_t i = 0; i < x.ncol(); ++i) {
        Point3 vertex = { x(0, i), x(1, i), x(2, i) };
        if (positions.count(vertex) == 0) {
            Rcpp::IntegerVector position = { i + 1 };
            positions.emplace(vertex, position);
        } else {
            duplicates.insert(vertex);
            positions.at(vertex).push_back(i + 1);
        }
    }

    Rcpp::List result;
    for (const Point3& vertex: duplicates) {
        result.push_back(positions.at(vertex));
    }

    return result;
}

On my machine I get a roughly 10-fold speedup compared to duplicated():

> set.seed(123)
> N <- 150
> M <- matrix(sample(c(1, 2, 3), N * 3, replace = TRUE), 3, N)
> 
> microbenchmark::microbenchmark(
+   Rcpp = dupes(M),
+   R = duplicated(t(M), fromLast = TRUE),
+   cpp_map = duplicated_vertices(M)
+ )
Unit: microseconds
    expr   min     lq    mean median    uq    max neval
    Rcpp 458.7 468.50 490.049  476.6 487.5 1092.7   100
       R 329.2 339.45 408.056  347.6 361.7 5103.3   100
 cpp_map  45.1  47.60  59.972   48.9  52.0  880.7   100