基于BowyerWatson的Delaunay三角化算法实现

实现效果如下图所示：

代码：

#include<iostream>
using namespace std;

#include "bowyer_watson.h"

#include <sstream>
#include <fstream>
#include <chrono>
#include <armadillo>
using namespace std;

//function to provide a first-order test of whether BowyerWatson::Triangle::contains() works correctly
void test()
{
    typedef BowyerWatson::Triangle Triangle;
    //create initial polyhedron
    arma::dvec3 u({ 0,0,1 });
    arma::dvec3 e({ 1,0,0 });
    arma::dvec3 n({ 0,1,0 });
    arma::dvec3 w({ -1,0,0 });
    arma::dvec3 s({ 0,-1,0 });
    arma::dvec3 d({ 0,0,-1 });
    Triangle* uen = new Triangle(u, e, n, 0);
    Triangle* unw = new Triangle(u, n, w, 0);
    Triangle* uws = new Triangle(u, w, s, 0);
    Triangle* use = new Triangle(u, s, e, 0);
    Triangle* des = new Triangle(d, e, s, 0);
    Triangle* dsw = new Triangle(d, s, w, 0);
    Triangle* dwn = new Triangle(d, w, n, 0);
    Triangle* dne = new Triangle(d, n, e, 0);
    uen->neighbors[0] = dne;    uen->neighbors[1] = unw;    uen->neighbors[2] = use;
    unw->neighbors[0] = dwn;    unw->neighbors[1] = uws;    unw->neighbors[2] = uen;
    uws->neighbors[0] = dsw;    uws->neighbors[1] = use;    uws->neighbors[2] = unw;
    use->neighbors[0] = des;    use->neighbors[1] = uen;    use->neighbors[2] = uws;
    des->neighbors[0] = use;    des->neighbors[1] = dsw;    des->neighbors[2] = dne;
    dsw->neighbors[0] = uws;    dsw->neighbors[1] = dwn;    dsw->neighbors[2] = des;
    dwn->neighbors[0] = unw;    dwn->neighbors[1] = dne;    dwn->neighbors[2] = dsw;
    dne->neighbors[0] = uen;    dne->neighbors[1] = des;    dne->neighbors[2] = dwn;
    Triangle* tris[8] = { uen, unw, use, uws, dne, dwn, des, dsw };
    //test point containment
    cout << "Should be ones:" << endl;
    for (unsigned int i = 0; i < 8; ++i)
    {
        arma::dvec3 to_test(
            {
                 1.0f * (tris[i]->vertices[0][0] + tris[i]->vertices[1][0] + tris[i]->vertices[2][0]) / 3.0f,
            1.0f * (tris[i]->vertices[0][1] + tris[i]->vertices[1][1] + tris[i]->vertices[2][1]) / 3.0f,
            1.0f * (tris[i]->vertices[0][2] + tris[i]->vertices[1][2] + tris[i]->vertices[2][2]) / 3.0f
            }
        );
        cout << (tris[i]->contains(to_test) ? '1' : '0');
    }
    cout << "\nShould be zeroes:" << endl;
    for (unsigned int i = 0; i < 8; ++i)
    {
        arma::dvec3 to_test(
            {
            -1.0f * (tris[i]->vertices[0][0] + tris[i]->vertices[1][0] + tris[i]->vertices[2][0]) / 3.0f,
            -1.0f * (tris[i]->vertices[0][1] + tris[i]->vertices[1][1] + tris[i]->vertices[2][1]) / 3.0f,
            -1.0f * (tris[i]->vertices[0][2] + tris[i]->vertices[1][2] + tris[i]->vertices[2][2]) / 3.0f
            }
        );
        cout << (tris[i]->contains(to_test) ? '1' : '0');
    }
    for (unsigned int j = 1; j < 8; ++j)
        for (unsigned int i = 0; i < 8; ++i)
        {
            arma::dvec3 to_test(
                {
                    1.0f * (tris[i]->vertices[0][0] + tris[i]->vertices[1][0] + tris[i]->vertices[2][0]) / 3.0f,
                1.0f * (tris[i]->vertices[0][1] + tris[i]->vertices[1][1] + tris[i]->vertices[2][1]) / 3.0f,
                1.0f * (tris[i]->vertices[0][2] + tris[i]->vertices[1][2] + tris[i]->vertices[2][2]) / 3.0f
                }
            );
            cout << (tris[(i + j) % 8]->contains(to_test) ? '1' : '0');
        }
    cout << endl;
}



//arguments: (pointcount = 1000, relaxations = 5, filename = bwoutput.py, seed = 666)
int main(int argc, char* argv[])
{
    unsigned int pointcount = 20;
    unsigned int relaxations = 5;
    std::string output_filename("bwoutput.py");
    unsigned int seed = 666;

    istringstream ss;
    bool input_erroneous(false);
    bool display_help(false);
    switch (argc)
    {
    case(5):
        ss.clear();
        ss.str(argv[4]);
        if (!(ss >> seed))
        {
            cout << "Received invalid input \"" << argv[4] << "\" as argument 4; this should be an integer specifying the seed to use for the random number generator. Use argument \"--help\" to see help menu." << endl;
            input_erroneous = true;
        }
    case(4):
        output_filename = argv[3];
        if (output_filename.substr(output_filename.length() - 3) != ".py")
            output_filename += ".py";
    case(3):
        ss.clear();
        ss.str(argv[2]);
        if (!(ss >> relaxations))
        {
            cout << "Received invalid input \"" << argv[2] << "\" as argument 2; this should be a nonnegative integer specifying the number of times to perform modified Lloyd relaxation. Use argument \"--help\" to see help menu." << endl;
            input_erroneous = true;
        }
    case(2):
        if (strcmp(argv[1], "--help") == 0 || strcmp(argv[1], "-help") == 0 || strcmp(argv[1], "help") == 0 || strcmp(argv[1], "\"--help\"") == 0)
        {
            display_help = true;
            input_erroneous = true;
            break;
        }
        ss.clear();
        ss.str(argv[1]);
        if (!(ss >> pointcount) || pointcount == 0)
        {
            cout << "Received invalid input \"" << argv[1] << "\" as argument 1; this should be a positive integer specifying the number of points to generate. Use argument \"--help\" to see help menu." << endl;
            input_erroneous = true;
        }
    case(1):
        break;
    default:
        display_help = true; //the OS provides argv[0], so this line should never be reached
        input_erroneous = true;
    }
    if (display_help)
    {
        cout << "\tBowyer-Watson Algorithm: Generates random points on the unit sphere, constructs a triangulation thereof, and creates a python file for viewing the resuling mesh.\n";
        cout << "\tTo view this menu, pass \"--help\" as the first argument.\n";
        cout << "\tArguments expected, in order, and their defaults:\n";
        cout << "\t\tpointcount = 1000: a positive integer specifying the number of points to generate.\n";
        cout << "\t\trelaxations = 5: a nonnegative integer specifying the number of times to relax the points (modified Lloyd relaxation).\n";
        cout << "\t\tfilename = bwoutput.py: the filename for the output python file. If it does not end in \".py\", that extension will be appended to it.\n";
        cout << "\t\tseed = 666: a nonnegative integer specifying the seed to be used for the random number generator used in generating points.";
        cout << endl;
    }
    if (input_erroneous)
        return(0);

    srand(seed);
    chrono::high_resolution_clock::time_point starttime;
    chrono::high_resolution_clock::time_point endtime;
    chrono::high_resolution_clock::time_point meshendtime;

    BowyerWatson bw;
    cout << "\tBeginning spherical Bowyer-Watson test with " << pointcount << " points & " << relaxations << " iterations of the relaxation algorithm..." << endl;
    starttime = chrono::high_resolution_clock::now();

    std::vector<std::pair<arma::dvec3, unsigned int> > point_vector;
    bw.generate_points(pointcount, point_vector);
    BowyerWatson::Triangle* results = bw.perform(bw.create_initial_triangles(point_vector));
    for (unsigned int relaxations_performed = 0; relaxations_performed < relaxations; ++relaxations_performed)
    {
        bw.relax_points(results, point_vector);
        results = bw.perform(bw.create_initial_triangles(point_vector));
    }
    endtime = chrono::high_resolution_clock::now();

    cout << "\tRunning time: " << chrono::duration_cast<chrono::milliseconds>(endtime - starttime).count() << " milliseconds." << endl;
    cout << "\tCompleted calculations. Transforming to mesh..." << endl;
    BowyerWatson::Mesh mesh = bw.get_mesh(results, true);
    meshendtime = chrono::high_resolution_clock::now();
    cout << "\tMesh running time: " << chrono::duration_cast<chrono::milliseconds>(meshendtime - endtime).count() << " milliseconds." << endl;
    cout << "\tw00t, you got a mesh! Writing to file..." << endl;

    ///
    ofstream file;
    file.open(output_filename);
    file << "from mpl_toolkits.mplot3d import Axes3D\n";
    file << "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n";
    file << "import matplotlib.pyplot as plt\n";
    file << "fig = plt.figure()\n";
    file << "ax = Axes3D(fig)\n";
    for (unsigned int i = 0; i < mesh.triangles.size(); i += 3)
    {
        file << "ax.add_collection3d(Poly3DCollection([list(zip([" <<
            mesh.vertices[mesh.triangles[i]] << "," << mesh.vertices[mesh.triangles[i + 1]] << "," << mesh.vertices[mesh.triangles[i + 2]]
            << "],[" <<
            mesh.vertices[mesh.triangles[i] + 1] << "," << mesh.vertices[mesh.triangles[i + 1] + 1] << "," << mesh.vertices[mesh.triangles[i + 2] + 1]
            << "],[" <<
            mesh.vertices[mesh.triangles[i] + 2] << "," << mesh.vertices[mesh.triangles[i + 1] + 2] << "," << mesh.vertices[mesh.triangles[i + 2] + 2]
            << "]))], facecolors='w', edgecolors='b'))\n";
    }
    file << "ax.set_xlim(-1.4, 1.4)\n";
    file << "ax.set_ylim(-1.4, 1.4)\n";
    file << "ax.set_zlim(-1.4, 1.4)\n";
    file << "plt.show()\n";
    file.close();
    ///

    cout << "\tBowyer-Watson test complete! Wrote to: " << output_filename << endl;
    return(0);
}