`GEO::PeriodicDelaunay3d` returns empty results when computing Power Diagrams

[Canjia-Huang]

Hi Bruno,

I'm using GEO::PeriodicDelaunay3d to compute 3D Power Diagrams. However, I frequently encounter cases where the result is empty (delaunay->nb_cells() == 0), with no errors reported under debug mode. In contrast, GEO::RegularWeightedDelaunay3d works correctly on the same input.

I suspect this might be due to differences in how the two functions interpret or constrain the weight parameter. Could you clarify if GEO::PeriodicDelaunay3d has specific constraints on input weights that I might be missing?

Here is a minimal example that reproduces the issue. The code randomly samples points in [0,1]^3 and assigns each a random weight between 0 and 1:

#include <geogram/basic/command_line_args.h>
#include <geogram/delaunay/periodic_delaunay_3d.h>
#include <geogram/delaunay/delaunay_3d.h>

int main(){
    GEO::initialize(GEO::GEOGRAM_INSTALL_ALL);
    GEO::CmdLine::import_arg_group("algo");
    GEO::CmdLine::import_arg_group("standard");

    constexpr GEO::index_t points_nb = 1e5;

    std::vector<double> points;
    points.reserve(3*points_nb);
    for (GEO::index_t i = 0; i < points_nb; ++i) {
        points.push_back(GEO::Numeric::random_float32());
        points.push_back(GEO::Numeric::random_float32());
        points.push_back(GEO::Numeric::random_float32());
    }

    std::vector<double> weights;
    weights.reserve(points_nb);
    for (GEO::index_t i = 0; i < points_nb; ++i) {
        weights.push_back(GEO::Numeric::random_float32());
    }

    std::vector<double> weighted_points;
    weighted_points.reserve(4*points_nb);
    for (GEO::index_t i = 0; i < points_nb; ++i) {
        weighted_points.push_back(points[3*i]);
        weighted_points.push_back(points[3*i+1]);
        weighted_points.push_back(points[3*i+2]);
        weighted_points.push_back(weights[i]);
    }

    /* == PeriodicDelaunay3d */

    {
        const GEO::SmartPointer delaunay = new GEO::PeriodicDelaunay3d(false);
        delaunay->set_vertices(points_nb, points.data());
        delaunay->set_weights(weights.data());
        delaunay->compute();

        std::cout << "PeriodicDelaunay3d nb_cells: " << delaunay->nb_cells() << std::endl;
    }

    /* == RegularWeightedDelaunay3d */

    {
        const GEO::SmartPointer delaunay = new GEO::RegularWeightedDelaunay3d();
        delaunay->set_vertices(points_nb, weighted_points.data());

        std::cout << "RegularWeightedDelaunay3d nb_cells: " << delaunay->nb_cells() << std::endl;
    }
}

The output I observe is:

PeriodicDelaunay3d nb_cells: 0
RegularWeightedDelaunay3d nb_cells: 37478

Let me know if I can provide some other details.
Thanks in advance!

Best,
Canjia

[BrunoLevy]

Hello,
There is a difference in the way to specify weights for both Delaunay classes (which is a bit counter-intuitive and stupid, but it is like that for historical reasons). A refactoring is planned for when I have time (don't know when it will be). Until then, it works as follows:

• PeriodicDelaunay3d::set_weighs() uses the standard definition of weights.
• RegularWeightedDelaunay3d uses 4d points (as you already figured out), the fourth coordinate should be set to $\sqrt{W-w_i}$ where $W = max_i(w_i)$ and where $w_i$ is the weight associated with vertex $i$ (see the comment near the declaration of RegularWeigthedDelaunay3d). If you want to know, why is it so ? Because if you compute the 4D Voronoi diagram of the so-lifted points then intersect it with the 3D space, you obtain the power diagram defined by the weights $w_i$ (and it is the way I computed it in an earlier version of geogram).

However, this does not explain why you obtain zero triangle with PeriodicDelaunay3d (when I've time I'll try your example and investigate further).