How to compute 3D Delaunay Triangulation ? (with or without library)

solub

Hi everyone,

I'm desperately and unsuccessfully trying to compute Delaunay triangulation from a 3D point cloud. It's a thing that many people have been trying to tackle for years (the same question has been asked over a dozen times on this forum: here, here, also here, here again...) but that has never found a proper solution.

To the best of my knowledge, there are 4 libraries capable of computing Delaunay triangulation:

• Mesh by Lee Byron
• iGeo by Satoru Sugihara
• Hemesh (WBlut)
• Toxiclibs by Karsten Schmidt

However:

• Mesh only works in 2 dimensions (+ is flawed)
• iGeo has not been ported yet to Processing 3
• Voronoi & DelaunayTriangulation classes from Toxiclibs do not support 3D

It seems therefore the only solution lies in the Hemesh library.

GOAL

I would like to transform the colors of a painting into a 3D mesh via Delaunay triangulation. It's a relatively widespread technique that media artist Quayola (among others) has been brilliantly using for years.

PROBLEMS (hard to explain in english)

• Instead of displaying a single 3D terrain (like in the pictures above) from a cloud of 3D points, the Hemesh library seems to compute multiple polyhedrons (like diamonds with cavities).
• There is a threshold (triangulation.getAlphaTriangles(threshold)) to limit the amount of surrounding points to connect to, but that very threshold prevents from connecting to neighboring points that are far.

As a result, this painting from John William Waterhouse:

(some dude getting hit on at a pool party)

is tranformed to this:

(all the shapes are "closed" (P3 connected to P1) but it shouldn't be that way)

With a higher threshold:

(Between two layers of polyhedrons)

QUESTIONS

• How can I avoid this behavior (if possible) ?
• Do you know another library that could compute Delaunay triangulations in 3D space ?
• Is there another way to achieve what I'm trying to do here (3D terrain from point cloud) ?

Important parts of the code (it's in Python but PLEASE DO NOT MOVE THIS QUESTION TO PYTHON MODE):

setup()

threshold = 35
render = WB_Render(this)
for p in points:
        new_points = WB_Point(p.x, p.y, p.z)
        list.append(new_points)
triangulation = WB_Triangulate.alphaTriangulate3D(list)
triangles = triangulation.getAlphaTriangles(threshold)

draw()

noStroke()
for i in range(0, len(triangles), 3):
    render.drawTriangle(liste[triangles[i]], liste[triangles[i+1]], liste[triangles[i+2]])

Thank you

Chrisir

How can I avoid this behavior ?

What do you mean exactly? With a higher threshold you reached your goal already, that "connecting to neighboring points that are far" - thing.

What don't you like about this result?

solub

Oh man that's really hard to explain...

1/ On the example pictures from Quayola's work you can see triangles of various sizes (large ones and small ones) whereas on my sketch you can clearly see that the triangles are about the same size.

2/ On the example pictures, every point is connected to its neighbours. On my sketch they are not ! P3 is connected to P1 instead of being connected to P4. That's why in the last picture you can see 2 layers of polyhedrons. There shouldn't be any polyhedrons in the first place, just a 3D terrain.

In other words, there are polyhedrons because Hemesh automatically closes the shapes (P3 with P1). With Toxiclibs a volume is closed only when specified (volume.closeSides()). I don't know if it's possible to keep a volume "open" with Hemesh.

Chrisir

you could have an dynamic threshold in that it varies randomly or for different areas on the screen

yould also make it so that the connetion is done where the same colors are involved

koogs

Where does the 3d come into this? You're starting with a 2d image, you're finishing with a 2d image (the last one you posted). I don't get where the 3rd dimension comes from.

Can you not do a 2d transformation and then randomly / noisily offset the points in z to give you a terrain? (Like the last but one picture you posted)

solub

@koogs Indeed I didn't explain the process, sorry. After running a Sobel operator (for edge detection) x and y are placed to the edges. z is the brightness of the pixel. The last picture I posted is actually in 3d, it shows the cavity formed inside the painting because of the missing connections between P3s and P4s.

jeremydouglass

@solub,

Am I right in understanding that you want a perfect 2D division of a plane from a particular perspective, and then you want to height-map the points of that surface? This approach that is sometimes called 2.5D

If so, one approach:

1. use the Mesh library by Lee Byron to create a 2D Delaunay triangulation
2. in a second pass over the Delaunay point array sample the brightness of each point to generate your height map.
3. render each triangle from your point+height array. here is a recent discussion of using triangulate to recover the triangles from a Delaunay point array: https://forum.processing.org/two/discussion/23246/delaunay-triangles

solub

@jeremydouglass I realize my explanation was incomplete and somewhat inconsistent. No, I'm not looking for a pseudo-3D approach. The sketch is in true 3D with a dynamic camera in order to navigate through the significant points of a painting (projected on a z axis) or around the terrain generated from those points. If I recall correctly Nicolas Clavaud's library (Triangulate) doesn't work in full 3D and I also should mention that it is unable to read array lists in Python mode.

I have contacted Frederik Vanhoutte (W-Blut) on github regarding this 3d triangulation issue that I have, not sure he'll respond but be assured I'll post his answer here if he ever does.

In the meantime I'm considering alternative approaches, notably the Poisson surface reconstruction method. Maybe you have any info on that ?

jeremydouglass

@solub -- I think you may be misunderstanding my point -- whether your call it 2.5D or not, it still uses real 3D coordinates in a real OpenGL P3D view with a real dynamic camera etc -- it is a height-map projection of a 2D mesh, but every point has an xyz.

Triangulate (or any equivalent code, such as the do-it-yourself code that appears in that thread) does not need to work in full 3D -- you instead run it on a 2D triangulation generated by Mesh in order to recover the 2D triangles. You now have your x+y dimensions for a true 3D surface. You then use point sampling to generate your z heights. You now have an x, y, and z for each point, and you know how the points go together to form facets, so you can render them using vertex or PShape et cetera.

I also should mention that it is unable to read array lists in Python mode.

I'm not able to help you there, but the algorithm for retrieving the triangles is I believe in that thread / in the Triangulate source, so perhaps you could implement that in Python if you wish.

GoToLoop

... and I also should mention that it is unable to read array lists in Python mode.

If you're sure it's actually a Java's ArrayList issue:
https://Docs.Oracle.com/javase/10/docs/api/java/util/ArrayList.html

You may import it via: from java.util import ArrayList

So you have the real deal, rather than depending on Jython's automatic List/ArrayList conversion. :)>-

solub

@GoToLoop great tips, I didn't know it was possible to play with real Java's ArrayList in Python mode. I just got the library working ! @jeremydouglass I was wrong, Triangulate library DO work with 3D points. Thank you for the explanations.

Thank you both of you ! It's such a relief to finally see a script working...

jeremydouglass

@solub -- very glad to hear you got it working!

If you have a very simple example of the working approach you ended up taking to 3D Delaunay triangulation of an image in Processing.py, it might be nice of you to share an example here for future forum users.

solub

@jeremydouglass I was about to post a 3D example sketch yesterday but I ran into an issue. I'll post the final script as soon as I manage to solve the problem.

solub

To run the script you'll need Peasycam and Triangulate libraries (download) + the sobel fragment shader below from the PostFx library (only needed for edge detection). I first had my own sobel operator to detect the edges but found that @cansik's shader implementation was much more accurate and efficient.

Also please note that this sketch is in Python.

sobelFrag.glsl

#ifdef GL_ES
precision mediump float;
precision mediump int;
#endif

#define PROCESSING_TEXTURE_SHADER

uniform sampler2D texture;

varying vec4 vertColor;
varying vec4 vertTexCoord;

uniform vec2 resolution;

void main(void) {
  float x = 1.0 / resolution.x;
  float y = 1.0 / resolution.y;
  vec4 horizEdge = vec4( 0.0 );
  horizEdge -= texture2D( texture, vec2( vertTexCoord.x - x, vertTexCoord.y - y ) ) * 1.0;
  horizEdge -= texture2D( texture, vec2( vertTexCoord.x - x, vertTexCoord.y     ) ) * 2.0;
  horizEdge -= texture2D( texture, vec2( vertTexCoord.x - x, vertTexCoord.y + y ) ) * 1.0;
  horizEdge += texture2D( texture, vec2( vertTexCoord.x + x, vertTexCoord.y - y ) ) * 1.0;
  horizEdge += texture2D( texture, vec2( vertTexCoord.x + x, vertTexCoord.y     ) ) * 2.0;
  horizEdge += texture2D( texture, vec2( vertTexCoord.x + x, vertTexCoord.y + y ) ) * 1.0;
  vec4 vertEdge = vec4( 0.0 );
  vertEdge -= texture2D( texture, vec2( vertTexCoord.x - x, vertTexCoord.y - y ) ) * 1.0;
  vertEdge -= texture2D( texture, vec2( vertTexCoord.x    , vertTexCoord.y - y ) ) * 2.0;
  vertEdge -= texture2D( texture, vec2( vertTexCoord.x + x, vertTexCoord.y - y ) ) * 1.0;
  vertEdge += texture2D( texture, vec2( vertTexCoord.x - x, vertTexCoord.y + y ) ) * 1.0;
  vertEdge += texture2D( texture, vec2( vertTexCoord.x    , vertTexCoord.y + y ) ) * 2.0;
  vertEdge += texture2D( texture, vec2( vertTexCoord.x + x, vertTexCoord.y + y ) ) * 1.0;
  vec3 edge = sqrt((horizEdge.rgb * horizEdge.rgb) + (vertEdge.rgb * vertEdge.rgb));

  gl_FragColor = vec4(edge, texture2D(texture, vertTexCoord.xy).a);
}

sketch.pyde

add_library('peasycam')
add_library('triangulate')
from java.util import ArrayList

step, threshold = 5, 100
liste = ArrayList()

def setup():
    global img, triangles
    size(1800, 1000, P3D)
    smooth(8)

    cam = PeasyCam(this, 1600)
    sobel = loadShader('sobelFrag.glsl')
    img = loadImage("image_file_name.jpg")

    pg = createGraphics(img.width, img.height, P2D)
    pg.beginDraw()
    pg.image(img, 0, 0)
    pg.filter(GRAY)
    pg.filter(sobel)
    pg.loadPixels()
    for x in range(0, pg.width, step):
        for y in range(0, pg.height, step):
            i = x + y * pg.width
            col = img.pixels[i]
            b = pg.pixels[i] & 0xFF
            if b > threshold:
                liste.add(PVector(x, y , brightness(col) * .5))

    for e in range(0, pg.width, int(pg.width/20)):
        liste.add(PVector(e, 0, 0))
        liste.add(PVector(e, pg.height, 0))
    for e in range(0, pg.height, int(pg.height/20)):
        liste.add(PVector(0, e, 0))
        liste.add(PVector(pg.width, e, 0))

    pg.endDraw()
    triangles = Triangulate.triangulate(liste)

    beginShape(TRIANGLES)
    noStroke()

def draw():
    background(0)

    translate(-img.width/2, -img.height/2)
    for t in triangles:
        x = int((t.p1.x + t.p2.x + t.p3.x) / 3)
        y = int((t.p1.y + t.p2.y + t.p3.y) / 3)
        i = x + y * img.width
        col = img.pixels[i]
        fill(col)
        vertex(t.p1.x, t.p1.y, t.p1.z)
        vertex(t.p2.x, t.p2.y, t.p2.z)
        vertex(t.p3.x, t.p3.y, t.p3.z)
    endShape() 

Cheers !

jeremydouglass

Really beautiful work!

solub

Thank you Jeremy, I'm really glad to hear you like it.

I just updated the script to add points to the edges of the picture.

solub

I have an issue involving the script above so I hope you won't mind me asking a question on this thread again.

I'm trying to add a bloom effect (from the PostFX library...again) to the sketch but for some reason I have 2 canvas being displayed on the screen:

• one at the right location (center) but without the bloom effect

• the other one, smaller, at the bottom right corner with the bloom effect

Do you guys have an idea what I'm doing wrong here ?

add_library('peasycam')
add_library('triangulate')
from java.util import ArrayList

step, threshold = 5, 100
liste = ArrayList()

def setup():
    global img, triangles, fx
    size(1800, 1000, P3D)
    smooth(8)

    cam = PeasyCam(this, 1600)
    sobel = loadShader('sobelFrag.glsl')
    img = loadImage("image_file_name.jpg")
    fx = PostFX(this)

    pg = createGraphics(img.width, img.height, P2D)
    pg.beginDraw()
    pg.image(img, 0, 0)
    pg.filter(GRAY)
    pg.filter(sobel)
    pg.loadPixels()
    for x in range(0, pg.width, step):
        for y in range(0, pg.height, step):
            i = x + y * pg.width
            col = img.pixels[i]
            b = pg.pixels[i] & 0xFF
            if b > threshold:
                liste.add(PVector(x, y , brightness(col) * .5))

    for e in range(0, pg.width, int(pg.width/20)):
        liste.add(PVector(e, 0, 0))
        liste.add(PVector(e, pg.height, 0))
    for e in range(0, pg.height, int(pg.height/20)):
        liste.add(PVector(0, e, 0))
        liste.add(PVector(pg.width, e, 0))

    pg.endDraw()
    triangles = Triangulate.triangulate(liste)

    beginShape(TRIANGLES)
    noStroke()

def draw():
    background(0)

    translate(-img.width/2, -img.height/2)
    for t in triangles:
        x = int((t.p1.x + t.p2.x + t.p3.x) / 3)
        y = int((t.p1.y + t.p2.y + t.p3.y) / 3)
        i = x + y * img.width
        col = img.pixels[i]
        fill(col)
        vertex(t.p1.x, t.p1.y, t.p1.z)
        vertex(t.p2.x, t.p2.y, t.p2.z)
        vertex(t.p3.x, t.p3.y, t.p3.z)
    endShape() 

    fx.render().bloom(0.5, 20, 40).compose()

solub

So, I tried to draw the vertices on a canvas and then applying a bloom effect on this canvas but i still get strange artifacts.

in setup

PostFX fx;
PGraphics canvas;

void setup()
{
  size(500, 500, P2D);

  fx = new PostFX(this);  
  canvas = createGraphics(width, height, P3D);
}

in draw

void draw()
{
  canvas.beginDraw();
  // draw something onto the canvas
  canvas.endDraw();

  blendMode(BLEND);
  image(canvas, 0, 0);

  // add bloom filter
  blendMode(SCREEN);
  fx.render(canvas)
    .brightPass(0.5)
    .blur(20, 50)
    .compose();
}

Maybe @cansik have an idea what could possibly be wrong ?

cansik

@solub Can you show me / tell me what "strange artefacts" are?

The problem with your python sketch is, that you are translating the coordinate system, but are not using the matrix stack. So the image (output of the post fx) will be drawn at the wrong location.

pushMatrix()
translate(-img.width/2, -img.height/2)
// do your drawing stuff
popMatrix()

// do post fx stuff

solub

Hi @cansik, thank you for your answer.

I first though it was because of the missing pushMatrix/popMatrix as well but pushing the transformation matrix doesn't solve the problem.

When drawing on a canvas:

• I get 2 pictures at 2 different locations
• One with the bloom effect (left), the other (right) without it
• The first one doesn't rotate when moving the camera, the second do
• The 3D disappears on both pictures (flattened 3D delaunay)

sketch.pyde

add_library('peasycam')
add_library('triangulate')
add_library('PostFX')
from java.util import ArrayList

step, threshold = 5, 100
liste = ArrayList()

def setup():
    global img, triangles, fx, canvas
    size(1800, 1000, P3D)
    smooth(8)

    fx = PostFX(this)
    cam = PeasyCam(this, 1600)
    sobel = loadShader('sobelFrag.glsl')
    img = loadImage("test75.jpg")

    canvas = createGraphics(width, height, P3D)

    pg = createGraphics(img.width, img.height, P2D)
    pg.beginDraw()
    pg.image(img, 0, 0)
    pg.filter(GRAY)
    pg.filter(sobel)
    pg.loadPixels()
    for x in range(0, pg.width, step):
        for y in range(0, pg.height, step):
            i = x + y * pg.width
            col = img.pixels[i]
            b = pg.pixels[i] & 0xFF
            if b > threshold:
                liste.add(PVector(x, y , brightness(col) * .5))

    for e in range(0, pg.width, int(pg.width/20)):
        liste.add(PVector(e, 0, 0))
        liste.add(PVector(e, pg.height, 0))
    for e in range(0, pg.height, int(pg.height/20)):
        liste.add(PVector(0, e, 0))
        liste.add(PVector(pg.width, e, 0))

    pg.endDraw()
    triangles = Triangulate.triangulate(liste)


def draw():
    background(0)

    canvas.pushMatrix()
    canvas.translate(-img.width/2, -img.height/2)
    canvas.beginDraw()
    canvas.beginShape(TRIANGLES)
    canvas.noStroke()
    for t in triangles:
        x = int((t.p1.x + t.p2.x + t.p3.x) / 3)
        y = int((t.p1.y + t.p2.y + t.p3.y) / 3)
        i = x + y * img.width
        col = img.pixels[i]
        canvas.fill(col)
        canvas.vertex(t.p1.x, t.p1.y, t.p1.z)
        canvas.vertex(t.p2.x, t.p2.y, t.p2.z)
        canvas.vertex(t.p3.x, t.p3.y, t.p3.z)
    canvas.endShape() 
    canvas.endDraw()
    canvas.popMatrix()

    blendMode(BLEND)
    image(canvas,0,0)

    blendMode(SCREEN)

    fx.render(canvas).bloom(0.5, 20, 40).compose()

cansik

@solub Please show me the code, where you applied the bloom effect, with push and pop of the matrix stack.

And by the way, if you are using peasycam, you have to do the post processing in the HUD state (beginHUD, endHUD).

It maybe makes sense to look at the examples of the libraries first. There is an example with the peasycam in the postfx library, and there is an offscreen rendering example in the peasycam library.

Looking at them should solve all your problems.

solub

Please show me the code, where you applied the bloom effect, with push and pop of the matrix stack.

I'm not sure to understand, the code is right above your comment. Maybe you're refering to the first script (which is basically the same but without an off-screen graphics buffer)

And by the way, if you are using peasycam, you have to do the post processing in the HUD state (beginHUD, endHUD).

It did solve the problem, thank you very much cansik. Also let me thank you once again for all the work you've put in the PostFX library.