quickhull3d does not always work if all points are on a same plane

[Hermann-SW]

I worked on JSCAd application lattice_sphere_cmp and made use of new geom3.fromPoints() which is based on quickhull3d. I realized that if all 3d points are on a same plane, the result might be wrong as in this bug.jscad simple JSCAD example:

I created JSCAD issue #1347, but was asked today to get quickhull3d fixed instead:
jscad/OpenJSCAD.org#1347 (comment)

Just using quickhull3d alone does show the bug as well (point 1 should appear in some faces, but does not):
jscad/OpenJSCAD.org#1347 (comment)

pi@raspberrypi5:~/quickhull3d $ node bug.js
[ [ 3, 0, 2 ], [ 4, 3, 2 ], [ 4, 2, 0 ], [ 4, 0, 3 ] ]
pi@raspberrypi5:~/quickhull3d $ cat bug.js 
import qh from 'quickhull3d'

const runner = qh.default

const pts = [[-8,1,0],[-7,4,0],[-6,5,-2],[-7,0,-4],[-8,0,-1]]

const faces = runner(pts)

console.log(faces)
pi@raspberrypi5:~/quickhull3d $ 

I already implemented fixes in the mentioned github issue — they are part of my lattice_sphere_cmp already. Basically a test is done whether all points are on a same plane. If not, normal quickhull3d is used. Otherwise a point known to be not on plane is added, then quickhull3d works nicely on the not all points on plane input, and finally all faces containing the added point are removed. Finally the result is the single face left, together with its inverse (to guarantee that face is visible from both sides). First animation shows what happens if only the remaining face is used in geom3, the second animation if both faces are used:

A good argument for having two faces in the result is looking at a tetrahedron where one edge is shortened until length 0:

The top and bottom face will remain.

The fix I have as discussed has to deal with non-integer distance of plane from origin, so comparison does not always work. I introduced a tolerance value for the functions that made it work, but I do not really like that.

For my lattice_sphere_cmp application all points are from ℤ³, and I used cross product instead of normalized normal vector to determine whether all points are on a same plane. That is clean, but all points in ℤ³ is not the normal use case for quickhull3d.

This is just issue report, no complete quickhull3d solution yet.
The code implemented can be found in these two diffs as well:
Hermann-SW/lattice_sphere_cmp@7834f9a#diff-93e00aa7bb721b3cc9d54a52895f396ada5a8dc0093c4c223ef9d7c39a1c145c
Hermann-SW/lattice_sphere_cmp@78225b5#diff-93e00aa7bb721b3cc9d54a52895f396ada5a8dc0093c4c223ef9d7c39a1c145c
The whole discussion on the algorithm is in this discord thread:
https://discord.com/channels/775984095250612234/1248251413629370389

Initial fromPointsConvexPlaneℤ3(pts) code is here:
Hermann-SW/lattice_sphere_cmp@6461507#diff-93e00aa7bb721b3cc9d54a52895f396ada5a8dc0093c4c223ef9d7c39a1c145c

display="centroid!=normal +face" demonstrates the code discussed to display selected faces:

[mauriciopoppe]

Thanks for the bug report, a similar bug was opened a few years ago in #5 where the input only contains points in a plane, I remember I added this test:

quickhull3d/test/index.test.js

Lines 271 to 275 in 2012ad0

	test('predefined set of points #5', function () {
	const points = require('./issue5.json')
	const faces = qh(points)
	expect(isConvexHull(points, faces)).toBe(true)
	})

And verified that the output is a convex hull using:

quickhull3d/test/index.test.js

Lines 22 to 47 in 2012ad0

	function isConvexHull (points, faces) {
	let i, j
	const n = points.length
	let nError = 0
	for (i = 0; i < faces.length; i += 1) {
	const normal = getPlaneNormal(
	[],
	points[faces[i][0]],
	points[faces[i][1]],
	points[faces[i][2]]
	)
	const offset = vec3.dot(normal, points[faces[i][0]])
	for (j = 0; j < n; j += 1) {
	if (faces[i].indexOf(j) === -1) {
	const aboveFace = vec3.dot(points[j], normal) > offset + EPS
	if (aboveFace) {
	console.log('points', points)
	console.log('face %j with index %d', faces[i], j)
	console.log('%d should be less than %d', vec3.dot(points[j], normal), offset)
	}
	nError += Number(aboveFace)
	}
	}
	}
	return nError === 0
	}

As you can see I'm using EPS of 1e-6.

---

I understand the motivation of outputing 2 faces if the input only contains coplanar points, as you mention this is not a normal input case. A similar problem happens in a lower dimension, e.g. find a triangle with 3 co-linear points, effectively there's no triangle and we only have a line segment.

If the input only contains coplanar points, we have a few options:

1. Crash the program because the input is not valid. This was the initial implementation but caused runtime errors so I went for option 2.
2. Provide an ordering of points and faces that doesn't break the guarantee that for each face all the other points are behind it (or in the same plane). This is the option that is currently implemented.
3. Create at least 2 opposite faces, keeping the guarantee that guarantee that for each face all the other points are behind it (or in the same plane). I'd argue that this is not a good output either with the triangle example above.

[z3dev]

Provide an ordering of points and faces that doesn't break the guarantee that for each face all the other points are behind it (or in the same plane). This is the option that is currently implemented.

Given the example, is this a problem of determining the orientation of the point to the face? It looks like the first point has been assigned to another face.

[Hermann-SW]

OK, no need for face to appear twice, that can be done by application using quickhull3d in postprocessing.

The 5 points of the example from description are in order, so complete face would be [0,1,2,3,4].
This script shows what quickhull3d does when adding a point outside the plane (origin here), and removing all faces containging new point (5). The result [ [ 4, 0, 1 ], [ 4, 1, 2 ], [ 4, 2, 3 ] ] is correct, all original five points 0..4 contained:

pi@raspberrypi5:~/quickhull3d $ node bug2.js
[[-8,1,0],[-7,4,0],[-6,5,-2],[-7,0,-4],[-8,0,-1]]
[[3,0,2],[4,3,2],[4,2,0],[4,0,3]]

[[-8,1,0],[-7,4,0],[-6,5,-2],[-7,0,-4],[-8,0,-1],[0,0,0]]
[[3,2,5],[1,0,5],[1,5,2],[4,0,1],[4,1,2],[4,2,3],[4,3,5],[4,5,0]]
[ [ 4, 0, 1 ], [ 4, 1, 2 ], [ 4, 2, 3 ] ]
pi@raspberrypi5:~/quickhull3d $ 
pi@raspberrypi5:~/quickhull3d $ cat bug2.js 
import qh from 'quickhull3d'

const runner = qh.default

const pts = [[-8,1,0],[-7,4,0],[-6,5,-2],[-7,0,-4],[-8,0,-1]]
const faces = runner(pts)
console.log(JSON.stringify(pts))
console.log(JSON.stringify(faces))
console.log()

const pts2 = [...pts,[0,0,0]]
const faces2 = runner(pts2)
console.log(JSON.stringify(pts2))
console.log(JSON.stringify(faces2))
const faces2b = faces2.reduce((a,p) =>
    (p.some((v)=>(v==pts2.length-1)))?a:a=[...a,p],[])
console.log(faces2b)
pi@raspberrypi5:~/quickhull3d $ 

So either Mauricio option (1) would be ok, reject input and let calling application deal with it (reposnse [] could indicate no faces, or coplanar points).
Or something else.

But not the current response [ [ 3, 0, 2 ], [ 4, 3, 2 ], [ 4, 2, 0 ], [ 4, 0, 3 ] ] because left top point 1 is in no triangle of response, so outside of convex hull contrary to definition of convex hull.

jscad app

Correct response:
[ [ 4, 0, 1 ], [ 4, 1, 2 ], [ 4, 2, 3 ] ]

[Hermann-SW]

Sorry, that response is not correct.
isPointInsideHull() is false for point [-6,5,2] for those faces.

Perhaps correct response is [ [ 0, 1, 2, 3, 4] ]?
All 5 points are "inside":

pi@raspberrypi5:~/quickhull3d $ node bug3.js
[[-8,1,0],[-7,4,0],[-6,5,-2],[-7,0,-4],[-8,0,-1]]
[ [ 0, 1, 2, 3, 4 ] ]
[ true, true, true, true, true ]
pi@raspberrypi5:~/quickhull3d $ cat bug3.js 
import qh from 'quickhull3d'

const isPointInsideHull = qh.isPointInsideHull

const pts = [[-8,1,0],[-7,4,0],[-6,5,-2],[-7,0,-4],[-8,0,-1]]
const faces = [ [ 0, 1, 2, 3, 4] ]

console.log(JSON.stringify(pts))
console.log(faces)
console.log(pts.map((p) => isPointInsideHull(p, pts, faces)))
pi@raspberrypi5:~/quickhull3d $ 

[mauriciopoppe]

Thanks for the graph, for the case and the expected output, I think this library can do a better job for this scenario so I have an algorithm that I haven't implemented yet.

• Before computing the initial tetrahedron check if the input is degenerate i.e. check if all the points are in the same plane
  • Pick 3 points and get a plane, project all the points towards the plane i.e. find the distance to the plane, if the projection distance is 0 for all the points then it means that all the points are in the same plane
• Translate the origin so that one of the 3 points becomes the origin, rotate the frame so that the plane is parallel to a plane formed by any 2 of the 3 axis e.g. any of the xy, xz, yz planes
• Discard the dimension where all the values are 0.
• Do a 2d convex hull
• Create a face where the indexes are the indexes of the convex hull points

[Hermann-SW]

Thanks Mauricio.
I tried to implement similar in JSCAD and tried to proces the generated geom2 with hull(), but that was not easy to get right.

There is only one issue with your approach:
... , if the projection distance is 0 for all the points ...
You need to Math.abs() compare to some tolerance instead.
And that depends on the max/min in the 3 dimensions.
A default tolerance of 1e-10 worked for me, but it is important to be able to raise that for inputs where max/min dimension values are very big numbers.
Below code is part of my lattice_sphere_cmp JSCAD app, it was used but currently is not (because I now use fromPointsConvexPlaneℤ3(pts) which allows to compare distance to 0 because all my application points are guaranteed to be in ℤ³ — I use cross procuct for "distance" computation instead of normalized normal, because cross product is guaranteed to have integer only entries in my scenario, and therefore all dot() results are integers). Below code is used by my fromPointsConvex(pts), which is workaround for geom3.fromPointsConvex() (issue #1347) that is based on quickhull3d:

const isPointOnPlane = (pla, pt, tol=1e-10) =>
  Math.abs(vec3.dot(pt, pla) - pla[3]) < tol

const arePointsOnPlane = (pla, pts, tol=1e-10) =>
  pts.every((p) => isPointOnPlane(pla, p, tol))

[mauriciopoppe]

This should be fixed in the alpha release 3.1.0-1 https://www.npmjs.com/package/quickhull3d/v/3.1.0-1, let me know how it goes and I'll create a stable release 3.1.0.

[Hermann-SW]

Hi,
I really want to test this.
I tried 3 methods and none worked :-(

1. From npm instructions, Cannot use import statement outside a module:

hermann@j4105:~$ npm install --save quickhull3d

up to date, audited 18 packages in 1s

found 0 vulnerabilities
hermann@j4105:~$ node bug.js 
(node:4689) Warning: To load an ES module, set "type": "module" in the package.json or use the .mjs extension.
(Use `node --trace-warnings ...` to show where the warning was created)
/home/hermann/bug.js:1
import qh from 'https://cdn.jsdelivr.net/npm/quickhull3d@3.1.0-1/+esm'
^^^^^^

SyntaxError: Cannot use import statement outside a module
    at wrapSafe (node:internal/modules/cjs/loader:1281:20)
    at Module._compile (node:internal/modules/cjs/loader:1321:27)
    at Module._extensions..js (node:internal/modules/cjs/loader:1416:10)
    at Module.load (node:internal/modules/cjs/loader:1208:32)
    at Module._load (node:internal/modules/cjs/loader:1024:12)
    at Function.executeUserEntryPoint [as runMain] (node:internal/modules/run_main:174:12)
    at node:internal/main/run_main_module:28:49

Node.js v20.13.1
hermann@j4105:~$ cat bug.js 
import qh from 'https://cdn.jsdelivr.net/npm/quickhull3d@3.1.0-1/+esm'

const runner = qh.default

const pts = [[-8,1,0],[-7,4,0],[-6,5,-2],[-7,0,-4],[-8,0,-1]]

const faces = runner(pts)

console.log(faces)
hermann@j4105:~$ 

[Hermann-SW]

2. in web developer tools:

import qh from 'https://cdn.jsdelivr.net/npm/quickhull3d@3.1.0-1/+esm'
Uncaught SyntaxError: import declarations may only appear at top level of a module debugger eval code:1:1

[Hermann-SW]

I did git pull, created package.json for module as learned from @z3dev
jscad/OpenJSCAD.org#1347 (comment)

Does not take your changes:

hermann@j4105:~/quickhull3d$ cat > package.json 
{
  "name": "new",
  "version": "1.0.0",
  "description": "a new project",
  "main": "index.js",
  "type": "module",
  "scripts": {
    "test": "echo \"Error: no test specified\" && exit 1"
  },
  "author": "",
  "license": "ISC",
  "dependencies": {
    "quickhull3d": "^2.1.0"
  }
}
hermann@j4105:~/quickhull3d$ npm install

removed 11 packages, changed 1 package, and audited 7 packages in 2s

found 0 vulnerabilities
hermann@j4105:~/quickhull3d$ node bug.js 
[ [ 3, 0, 2 ], [ 4, 3, 2 ], [ 4, 2, 0 ], [ 4, 0, 3 ] ]
hermann@j4105:~/quickhull3d$ 

[Hermann-SW]

Please let me know how I can undo whatever I did and test your current repo.
Can you please run bug.js locally on your computer? And show output?

[mauriciopoppe]

1. This module is now an ESmodule, in node please read https://nodejs.org/api/esm.html#enabling about how to enable it.
2. You can load it in a browser but not through devtools because it needs to be a file imported with <script type="module" src="main.js"></script>, in fact, this is what I'm using in gh-pages https://github.com/mauriciopoppe/quickhull3d/blob/gh-pages/index.html
3. The package is published as alpha so you need to explicitly use the version I sent above e.g. "quickhull3d": "3.1.0-1" and reinstall.

[Hermann-SW]

Hi,
thanks for the infos.
But I was not able to get anything to work but in browser.
Be aware that the "Minimal browser demo" output for your example changed
(now: [[2,1,0],[3,1,2],[3,0,1],[3,2,0]]):
https://www.npmjs.com/package/quickhull3d/v/3.1.0-1

This is above bug.js for the browser as bug.html:

<!DOCTYPE html>
<html>
<head><meta charset="UTF-8" /></head>

<body>
 <script type="module">
  import qh from 'https://cdn.jsdelivr.net/npm/quickhull3d@3.1.0-1/+esm'

  const pts = [[-8,1,0],[-7,4,0],[-6,5,-2],[-7,0,-4],[-8,0,-1]]

  const faces = qh(pts)

  alert(JSON.stringify(faces))
 </script>
</body>

</html>

And it correctly alerts now this:

[[2,1,0],[2,0,4],[2,4,3]]

Corresponding to this ASCII art diagram, similar to JSCAD screenshot above with bug output, but this time correct:

   1---0
  /   / \
2 ---/    \
| \        \
|   \ ------4
|          /
+---------3

[mauriciopoppe]

Thanks, I updated the output of the example in the README.

[z3dev]

@mauriciopoppe thanks a lot for the fix. And of course, this very nice library.

[mauriciopoppe]

Fix released in https://github.com/mauriciopoppe/quickhull3d/releases/tag/v3.1.0