b2find_coinciding_points.H

//------------------------------------------------------------------------
// b2find_coinciding_points.H --
//
//     Find coinciding points.
//
// written by Thomas Ludwig
//
// Copyright (c) 2016
//     SMR Engineering & Development SA
//     2502 Bienne, Switzerland
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef _B2_FIND_COINCIDING_POINTS_H_
#define _B2_FIND_COINCIDING_POINTS_H_
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <vector>
 
namespace {
 
struct Signature {
    uint64_t d[3];
 
    bool operator==(const Signature& o) const {
        return d[0] == o.d[0] && d[1] == o.d[1] && d[2] == o.d[2];
    }
 
    bool operator!=(const Signature& o) const {
        return d[0] != o.d[0] || d[1] != o.d[1] || d[2] != o.d[2];
    }
 
    bool operator<(const Signature& o) const {
        return (
              d[0] < o.d[0]
              || (d[0] == o.d[0] && (d[1] < o.d[1] || (d[1] == o.d[1] && d[2] < o.d[2]))));
    }
};
 
struct Point {
    size_t index;
    Signature signature;
 
    bool operator<(const Point& o) const {
        return (signature < o.signature || (signature == o.signature && index < o.index));
    }
};
 
}  // namespace
 
namespace b2000 {
 
template <typename T = uint64_t>
inline void find_coinciding_points(
      const double delta, const bool delta_relative, const size_t n, const double* coor,
      std::vector<size_t>& coincide) {
    assert(delta > 0);
 
    coincide.clear();
 
    // Special case.
    if (n == 0) { return; }
 
    coincide.reserve(n);
    for (size_t i = 0; i != n; ++i) { coincide.push_back(i); }
 
    // Get axis-aligned bounding box.
    assert(n > 0);
    double bb_min[3] = {coor[0], coor[1], coor[2]};
    double bb_max[3] = {coor[0], coor[1], coor[2]};
    for (size_t i = 3; i != 3 * n; ++i) {
        size_t j = i % 3;
        if (coor[i] < bb_min[j]) { bb_min[j] = coor[i]; }
        if (coor[i] > bb_max[j]) { bb_max[j] = coor[i]; }
    }
 
    // Calculate the maximum distance in each direction by which two
    // "coinciding" points may be separated.
    double delta_abs = delta;
    if (delta_relative) {
        const double bb_extent[3] = {
              bb_max[0] - bb_min[0], bb_max[1] - bb_min[1], bb_max[2] - bb_min[2]};
 
        const double bb_extent_max = std::max(bb_extent[0], std::max(bb_extent[1], bb_extent[2]));
 
        delta_abs = bb_extent_max * delta;
    }
    const double delta_abs_2 = delta_abs * delta_abs;
 
    // The size of the boxes must be twice as large.
    const double box_size = 2 * delta_abs;
 
    std::vector<Point> points(n, Point());
 
    // First pass.
    double box_start_a[3] = {bb_min[0], bb_min[1], bb_min[2]};
    for (size_t i = 0; i != n; ++i) {
        Point& p = points[i];
        p.index = i;
        for (int j = 0; j != 3; ++j) {
            p.signature.d[j] = T((coor[3 * i + j] - box_start_a[j]) / box_size);
        }
    }
    std::sort(points.begin(), points.end());
    for (size_t i = 0; i + 1 < n;) {
        const size_t j = i++;
        for (; i < n; ++i) {
            if (points[i].signature != points[j].signature) { break; }
            const size_t ii = points[i].index;
            const size_t jj = points[j].index;
            const double* ci = &coor[3 * ii];
            const double* cj = &coor[3 * jj];
            const double d[3] = {ci[0] - cj[0], ci[1] - cj[1], ci[2] - cj[2]};
            if (coincide[ii] == ii && d[0] * d[0] + d[1] * d[1] + d[2] * d[2] < delta_abs_2) {
                coincide[ii] = coincide[jj];
            }
        }
    }
 
    // Second pass with offset boxes.
    double box_start_b[3] = {
          bb_min[0] + 0.5 * box_size, bb_min[1] + 0.5 * box_size, bb_min[2] + 0.5 * box_size};
    for (size_t i = 0; i != n; ++i) {
        Point& p = points[i];
        p.index = i;
        for (int j = 0; j != 3; ++j) {
            p.signature.d[j] = T((coor[3 * i + j] - box_start_b[j]) / box_size);
        }
    }
    std::sort(points.begin(), points.end());
    for (size_t i = 0; i + 1 < n;) {
        const size_t j = i++;
        for (; i < n; ++i) {
            if (points[i].signature != points[j].signature) { break; }
            const size_t ii = points[i].index;
            const size_t jj = points[j].index;
            const double* ci = &coor[3 * ii];
            const double* cj = &coor[3 * jj];
            const double d[3] = {ci[0] - cj[0], ci[1] - cj[1], ci[2] - cj[2]};
            if (coincide[ii] == ii && d[0] * d[0] + d[1] * d[1] + d[2] * d[2] < delta_abs_2) {
                coincide[ii] = coincide[jj];
            }
        }
    }
 
    // Resolve transitive relations due to the second pass.
    for (size_t i = 0; i != n; ++i) { coincide[i] = coincide[coincide[i]]; }
}
 
template <typename T = uint64_t>
inline void find_coinciding_points(
      const double delta, const bool delta_relative, const size_t n, const double* coor,
      std::vector<std::pair<size_t, size_t> >& coincide) {
    coincide.clear();
 
    std::vector<size_t> vcoincide;
    find_coinciding_points<T>(delta, delta_relative, n, coor, vcoincide);
 
    assert(vcoincide.size() == n);
    coincide.reserve(n);
    for (size_t i = 0; i != n; ++i) {
        coincide.push_back(std::pair<size_t, size_t>(vcoincide[i], i));
    }
    std::sort(coincide.begin(), coincide.end());
}
 
}  // namespace b2000
 
#endif /* _B2_FIND_COINCIDING_POINTS_H_ */