java - 如何用三角形绘制 3D 地形？

Question

我有一个代表一些 3D 地形的 3D 布尔数组。目前我可以通过在数组中由它的xy和z指定的位置绘制一个点来绘制它，它看起来像这样。

我不知道如何使用三角形来绘制它，所以它看起来像真实的地形。我也不想把每个都画成一个立方体。

是否有任何算法可以获得要绘制的点（请记住，为了提高效率，只应绘制陆地外部的点）？

Answer 1

Absolutely amazing question! I couldn't resist but to play with boxes, since boxes are sexy. It is actually fairly easy to produce boxes with the hidden faces omitted.

The following algorithm takes a list of 3D positions of true in your grid, which is easy to obtain by simply scanning through the grid and filling an array. Also with this data format you can store much larger grid, provided that the terrain is reasonably sparse. Up front, apologies for my spartan for-each loops, I just wanted to have the code in-place and avoid writing dozens of function objects or using lambdas. Also apologies for using C++ instead of Java, don't have javac at home and I'm not really good with it anyway. Here goes:

#include <vector>
#include <map>
#include <set>
#include <utility>
#include <assert.h>
#include <math.h>

struct Pos3 {
    int x, y, z;

    Pos3(int _x = 0, int _y = 0, int _z = 0);
    bool operator <(const Pos3 &other) const;
};

std::vector<int> index_buffer;
std::vector<float> vertex_buffer;

void onInitialize()
{
    const int N = 32;
    std::vector<Pos3> points;
    GeneratePoints(points, N);
    // input: bunch of points in NxNxN box (easy to get from a boolean array,
    // can have much larger terrains if stored like this)

    std::set<Pos3> point_set;
    point_set.insert(points.begin(), points.end());
    // put all the points to a set to be able to lookup neighbors (not needed with an array)

    std::vector<std::vector<int> > polygons;
    polygons.reserve(3 * points.size()); // guess
    std::map<Pos3, int> vertex_map;
    for(size_t i = 0, n = points.size(); i < n; ++ i) {
        Pos3 p = points[i], corners[8] = {
            p, Pos3(p.x + 1, p.y, p.z), Pos3(p.x + 1, p.y + 1, p.z), Pos3(p.x, p.y + 1, p.z),
            Pos3(p.x, p.y, p.z + 1), Pos3(p.x + 1, p.y, p.z + 1), Pos3(p.x + 1, p.y + 1, p.z + 1),
            Pos3(p.x, p.y + 1, p.z + 1)
        };
        // get corners of a cube

        static const int sides[][3 + 4] = {
            0, -1,  0,  4, 5, 1, 0,          1, 0, 0,   5, 6, 2, 1,
            0,  1,  0,  6, 7, 3, 2,         -1, 0, 0,   7, 4, 0, 3,
            0,  0, -1,  0, 1, 2, 3,          0, 0, 1,   7, 6, 5, 4
        };
        // directions and side quad indices

        for(int j = 0; j < 6; ++ j) {
            Pos3 n(p.x + sides[j][0], p.y + sides[j][1], p.z + sides[j][2]); // position of a neighbor
            if(point_set.find(n) != point_set.end())
                continue; // have a neighbor, will not triangulate this side

            polygons.resize(polygons.size() + 1);
            std::vector<int> &poly = polygons.back(); // or use emplace_back() in c++11
            poly.resize(4); // form quads
            for(int v = 0; v < 4; ++ v) {
                Pos3 vert = corners[sides[j][3 + v]];
                std::map<Pos3, int>::iterator it; // use map to reuse vertices
                if((it = vertex_map.find(vert)) == vertex_map.end())
                    vertex_map[vert] = poly[v] = vertex_map.size(); // new vertex
                else
                    poly[v] = (*it).second; // existing vertex
            }
        }
        // generate sides, skip invisible sides
        // note that this still triangulates cavities, would have to flood-fill
        // outside area and then set all that is not outside to opaque (did not
        // solve that as this is also a valid behavior)
    }

    vertex_buffer.resize(vertex_map.size() * 3);
    for(std::map<Pos3, int>::const_iterator it = vertex_map.begin(), e = vertex_map.end(); it != e; ++ it) {
        size_t i = (*it).second * 3;
        vertex_buffer[i + 0] = ((*it).first.x + .5f) / (N + 1) * 2 - 1;
        vertex_buffer[i + 1] = ((*it).first.y + .5f) / (N + 1) * 2 - 1;
        vertex_buffer[i + 2] = ((*it).first.z + .5f) / (N + 1) * 2 - 1;
    }
    // convert points from the discrete domain
    // to a unit 3D cube centered around the origin

    index_buffer.reserve(polygons.size() * 2 * 3); // approximate number of triangles
    for(size_t i = 0, n = polygons.size(); i < n; ++ i) {
        const std::vector<int> &poly = polygons[i];
        for(size_t j = 2, n = poly.size(); j < n; ++ j) {
            index_buffer.push_back(poly[0]);
            index_buffer.push_back(poly[j]);
            index_buffer.push_back(poly[j - 1]);
        }
    }
    // convert polygons (those are actually quads) to triangles
}

There is also some more code that generates normals (ommited for the sake of clarity), the output looks like this:

The shape is a Julia set generated on a discrete lattice, you might recognize the shape when you turn it arround.

This is actually pretty similar to what you would get by the Delaunay triangulation if you could easily remove your interior points. The generated shape is hollow. There can be some "bubbles" in the shape, in case the booleans also contain a bubble (does not occur with Julia). This is easily fixed by flood-filling the booleans in order to fill those up.

Next, we can apply Catmull-Clark subdivision in order to get a smoother mesh:

typedef std::map<std::pair<int, int>, std::pair<size_t, int> > EdgeMap;

static bool Get_EdgeID(size_t &eid, int a, int b, EdgeMap &edges)
{
    std::pair<int, int> e(std::min(a, b), std::max(a, b));
    EdgeMap::iterator it = edges.find(e);
    if(it == edges.end()) {
        edges[e] = std::make_pair(eid = edges.size(), 1); // id, count
        return true; // new edge
    } else {
        eid = (*it).second.first; // get id
        ++ (*it).second.second; // increase count
        return false; // no new edge
    }
}

void CatClark(std::vector<std::vector<int> > &src_quads, std::vector<float> &src_verts)
{
    const static float vpw[4] = {9.0f, 3.0f, 1.0f, 3.0f};
    const static float epw[4] = {3.0f, 3.0f, 1.0f, 1.0f};

    std::vector<std::vector<int> > dst_quads(src_quads.size() * 4, std::vector<int>(4)); // will produce quads
    std::vector<float> dst_verts(src_verts.size() + src_quads.size() * 3, 0); // alloc s¨pace for vertices

    EdgeMap edges;
    std::vector<int> face_valences(src_verts.size() / 3, 0);
    const size_t off_vp = src_quads.size(), off_ep = off_vp + src_verts.size() / 3;
    for(size_t j = 0; j < off_vp; ++ j) {
        assert(src_quads[j].size() == 4); // otherwise won't work
        size_t eid[4];
        for(int k = 0; k < 4; ++ k) {
            int quad[4];
            for(int i = 0; i < 4; ++ i)
                quad[i] = src_quads[j][(i + k) & 3]; // get the 4 vertices (but rotate them each k iteration)
            if(Get_EdgeID(eid[k], quad[0], quad[1], edges)) // create edges
                dst_verts.insert(dst_verts.end(), 3, .0f); // must add new vertex to accomodate subdivided edge point
            ++ face_valences[quad[0]]; // update face-valence
            for(int n = 0; n < 3; ++ n)
                dst_verts[j * 3 + n] += 0.25f * src_verts[quad[0] * 3 + n]; // increment face point 
            for(int i = 0; i < 4; ++ i) {
                for(int n = 0; n < 3; ++ n) {
                    dst_verts[(off_vp + quad[0]) * 3 + n] += vpw[i] * src_verts[quad[i] * 3 + n]; // incremente vertex point
                    dst_verts[(off_ep + eid[k]) * 3 + n] += epw[i] * src_verts[quad[i] * 3 + n]; // increment edge point
                }
            }
        }
        for(int k = 0; k < 4; ++ k) { // make child faces
            dst_quads[4 * j + k][0] = j;
            dst_quads[4 * j + k][4] = off_ep + eid[(3 + k) & 3]; 
            dst_quads[4 * j + k][5] = off_ep + eid[(0 + k) & 3]; 
            dst_quads[4 * j + k][6] = off_vp + src_quads[j][k];
        }
    }
    for(size_t j = 0, n = src_verts.size() / 3; j < n; ++ j) {
        for(int n = 0; n < 3; ++ n)
            dst_verts[(off_vp + j) * 3 + n] *= 0.0625f / float(face_valences[j]);
    }
    for(EdgeMap::const_iterator it = edges.begin(), e = edges.end(); it != e; ++ it) {
        size_t j = (*it).second.first;
        float rvalence = 0.1250f / float((*it).second.second);
        for(int n = 0; n < 3; ++ n)
            dst_verts[(off_ep + j) * 3 + n] *= rvalence;
    }
    dst_quads.swap(src_quads);
    dst_verts.swap(src_verts);
}

This algorithm was adapted to work with STL containers from Iñigo 'iq' Quilez / rgba, "Tricks and techniques for rgba's past and future intros", Breakpoint, 2007.

This gives output quite similar to what you would get with marching cubes / marching tetrahedrons / marching triangles, except that it is always higher resolution than that of the original lattice (with the above methods you can easily change the triangulation resolution). A slightly different view of the same data:

Or without wireframe:

The complete source code, along with Visual Studio workspace and win32 binaries can be found at here. It uses GLUT and the old fixed-function pipeline to display the generated geometry (only for simplicity and portability, otherwise I'd have to include GLEW or the like). I really enjoyed toying around with your question, I hope you will like the output ...

If you want to use marching cubes, you can find many demos online, or check out this real-time water simulation I did a couple years ago.

Answer 2

如果您有描述地形外壳的点，则此软件包非常好（快速）用于计算该组点的 Delaunay 三角剖分：

http://www.cs.bgu.ac.il/~benmoshe/DT/

您可以稍后从三角剖分中绘制每个三角形。我做了几个或方法将所有内容转换回双精度以与 JOGL 一起使用，您可能会觉得有帮助：

public static ArrayList<Point_dt[]> DTtoTriListDT(Delaunay_Triangulation DT){
ArrayList<Point_dt[]> triangles = new ArrayList<Point_dt[]>();
Point_dt[] triangle = new Point_dt[3];

Iterator<Triangle_dt> surface = DT.trianglesIterator();
while(surface.hasNext()){
    Triangle_dt tri = surface.next();
    triangle[0] = tri.p1();
    triangle[1] = tri.p2();
    triangle[2] = tri.p3();
    triangles.add(triangle);
}

return triangles;}

和

public static ArrayList<double[][]> DTtoTriList(Delaunay_Triangulation DT){
ArrayList<Point_dt[]> trianglesdt = Algebra.DTtoTriListDT(DT);
ArrayList<double[][]> triangles = new ArrayList<double[][]>();
double[][] triangle = new double[3][3];

Iterator<Point_dt[]> surface = trianglesdt.iterator();
while(surface.hasNext()){
    Point_dt[] tri = surface.next();

    triangle[0][0] = tri[0].x();
    triangle[0][1] = tri[0].y();
    triangle[0][2] = tri[0].z();

    triangle[1][0] = tri[1].x();
    triangle[1][1] = tri[1].y();
    triangle[1][2] = tri[1].z();

    triangle[2][0] = tri[2].x();
    triangle[2][1] = tri[2].y();
    triangle[2][2] = tri[2].z();

    triangles.add(triangle);
}

return triangles;

}

我在代数课（FYI）上都有。使用最后一种方法，您将获得一个 ArrayList，其中每个条目上包含三个双精度的树集（每个集合包含每个点的三个坐标，每个条目是一个三角形）。