3D volume fails on tiny point sets #321
Comments
bqi343
changed the title
|
There are also many instances where no error is thrown but both the default library and cddlib give twice the volume of that returned by qhull. @testset "test_volume_tricky3" begin
points = [-1.0 1.0 1.0
-1.0 -1.0 1.0
1.0 0.0 -1.0
1.0 1.0 1.0
0.0 1.0 -1.0]
for lib in libs # default -> 5.333333333333333, qhull -> 2.6666666666666665 (qhull is correct)
vol = get_volume(points, lib, false)
@test vol ≈ 8.0 / 3
if !(vol ≈ 8.0 / 3)
println("failed $lib got $vol")
end
end
end
|
|
Hello guys,
With our integration tool, based on boundary triangulation, result is 2.0 = 8/4
The full code for example construction, visualization and computation is included.
The whole package will be fully rewritten and optimized after the summer.
Let me know if you have questions-
—alberto

… On Mar 26, 2023, at 7:14 PM, Benjamin Qi ***@***.***> wrote:
There are also many instances where no error is thrown but both the default library and cddlib give twice the volume of that returned by qhull.
@testset "test_volume_tricky3" begin
points = [-1.0 1.0 1.0
-1.0 -1.0 1.0
1.0 0.0 -1.0
1.0 1.0 1.0
0.0 1.0 -1.0]
for lib in libs # default -> 5.333333333333333, qhull -> 2.6666666666666665 (qhull is correct)
vol = get_volume(points, lib, false)
@test vol ≈ 8.0 / 3
if !(vol ≈ 8.0 / 3)
println("failed $lib got $vol")
end
end
end
—
Reply to this email directly, view it on GitHub <#321 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAG2COMFH7NMT3AK7ZTC26TW6B2QBANCNFSM6AAAAAAWIIHOD4>.
You are receiving this because you are subscribed to this thread.
|
|
I don't see any included visualization, and I'm pretty sure 8/4 is wrong considering that
import numpy as np
from scipy.spatial import ConvexHull
points = np.array([[-1.0, 1.0, 1.0],
[-1.0, -1.0, 1.0],
[1.0, 0.0, -1.0],
[1.0, 1.0, 1.0],
[0.0, 1.0, -1.0]])
print(ConvexHull(points).volume) # 2.6666666666666665
using Polyhedra: polyhedron, vrep, volume
using LinearAlgebra
points = [-1.0 1.0 1.0
-1.0 -1.0 1.0
1.0 0.0 -1.0
1.0 1.0 1.0
0.0 1.0 -1.0]
function volume_simplex(points)
A = Matrix{T}(undef, 3, 3)
for i in 1:3
A[i, :] = points[i+1, :] - points[1, :]
end
return abs(det(A)) / 6
end
p = polyhedron(vrep(points))
println(volume(p)) # 5.3..., wrong
vol_1 = volume_simplex(points[[1, 2, 4, 5], :])
vol_2 = volume_simplex(points[[2, 4, 5, 3], :])
println("$vol_1 $vol_2 $(vol_1 + vol_2)") # 1.3..., 1.3..., 2.6..., as expected
|
|
Your visualization:

Good Luck !!
Alberto _
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.8.5 (2023-01-08)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
julia> # polyhedron visualization and integration test
# ---------------------------------------------
using ViewerGL
julia> using LinearAlgebraicRepresentation
julia> GL = ViewerGL
ViewerGL
julia> Lar = LinearAlgebraicRepresentation
LinearAlgebraicRepresentation
julia> # your vertices
W = [-1.0 1.0 1.0
-1.0 -1.0 1.0
1.0 0.0 -1.0
1.0 1.0 1.0
0.0 1.0 -1.0]
5×3 Matrix{Float64}:
-1.0 1.0 1.0
-1.0 -1.0 1.0
1.0 0.0 -1.0
1.0 1.0 1.0
0.0 1.0 -1.0
julia> # Transform in our vertex repr (transpose)
W = convert(Matrix,W')
3×5 Matrix{Float64}:
-1.0 -1.0 1.0 1.0 0.0
1.0 -1.0 0.0 1.0 1.0
1.0 1.0 -1.0 1.0 -1.0
julia> # added two more triangle faces (your polyhedron is not linear)
# it has a face with four non-coplanar points (my last two triangles)
TW = [[1,2,4],[2,3,4],[3,5,4],[1,4,5],[1,3,2],[1,5,3]]
6-element Vector{Vector{Int64}}:
[1, 2, 4]
[2, 3, 4]
[3, 5, 4]
[1, 4, 5]
[1, 3, 2]
[1, 5, 3]
julia> # LAR representation of 8 unit cubes in positive octant
V,(VV,EV,FV,CV) = Lar.cuboidGrid([2,2,2],true)
([0.0 0.0 … 2.0 2.0; 0.0 0.0 … 2.0 2.0; 0.0 1.0 … 1.0 2.0], [[[1], [2], [3], [4], [5], [6], [7], [8], [9], [10] … [18], [19], [20], [21], [22], [23], [24], [25], [26], [27]], [[1, 2], [2, 3], [4, 5], [5, 6], [7, 8], [8, 9], [10, 11], [11, 12], [13, 14], [14, 15] … [9, 18], [10, 19], [11, 20], [12, 21], [13, 22], [14, 23], [15, 24], [16, 25], [17, 26], [18, 27]], [[1, 2, 4, 5], [2, 3, 5, 6], [4, 5, 7, 8], [5, 6, 8, 9], [10, 11, 13, 14], [11, 12, 14, 15], [13, 14, 16, 17], [14, 15, 17, 18], [19, 20, 22, 23], [20, 21, 23, 24] … [3, 6, 12, 15], [4, 7, 13, 16], [5, 8, 14, 17], [6, 9, 15, 18], [10, 13, 19, 22], [11, 14, 20, 23], [12, 15, 21, 24], [13, 16, 22, 25], [14, 17, 23, 26], [15, 18, 24, 27]], [[1, 2, 4, 5, 10, 11, 13, 14], [2, 3, 5, 6, 11, 12, 14, 15], [4, 5, 7, 8, 13, 14, 16, 17], [5, 6, 8, 9, 14, 15, 17, 18], [10, 11, 13, 14, 19, 20, 22, 23], [11, 12, 14, 15, 20, 21, 23, 24], [13, 14, 16, 17, 22, 23, 25, 26], [14, 15, 17, 18, 23, 24, 26, 27]]])
julia> # translate your polyhedral model to have the conic vertex in [2,2,2]'
W,TW = Lar.apply(Lar.t(1,1,1))((W,TW))
([0.0 0.0 … 2.0 1.0; 2.0 0.0 … 1.0 2.0; 2.0 2.0 … 0.0 0.0], [[1, 2, 3], [2, 4, 3], [4, 5, 3], [1, 3, 5], [1, 4, 2], [1, 5, 4]])
julia> # you may use the mouse to rotate/scale
GL.VIEW([
GL.GLGrid(V,EV,GL.COLORS[1],1.)
GL.GLGrid(W,TW,GL.COLORS[6],0.5)
GL.GLAxis(GL.Point3d(1,1,1),GL.Point3d(2,2,2))
]);
julia> # computation of 4x4 affine matrix of inertia and volume/surface using boundary triangulation, according to:
# C. Cattani and A. Paoluzzi: Boundary integration over linear polyhedra. Computer-Aided Design 22(2): 130-135 (1990) (doi:10.1016/0010-4485(90)90007-Y)
pol = (W,TW)
([0.0 0.0 … 2.0 1.0; 2.0 0.0 … 1.0 2.0; 2.0 2.0 … 0.0 0.0], [[1, 2, 3], [2, 4, 3], [4, 5, 3], [1, 3, 5], [1, 4, 2], [1, 5, 4]])
julia> vol = Lar.volume(pol)
2.0
julia>
… On Mar 27, 2023, at 4:16 PM, Benjamin Qi ***@***.***> wrote:
I don't see any included visualization, and I"m pretty sure 8/4 is wrong considering that
QHull agrees with 8/3:
import numpy as np
from scipy.spatial import ConvexHull
points = np.array([[-1.0, 1.0, 1.0],
[-1.0, -1.0, 1.0],
[1.0, 0.0, -1.0],
[1.0, 1.0, 1.0],
[0.0, 1.0, -1.0]])
print(ConvexHull(points).volume) # 2.6666666666666665
From manual inspection, the hull splits into the tetrahedrons labeled by $(1, 2, 4, 5)$ and $(2, 3, 4, 5)$, both of which have a volume of $8/6$, for a total of $8/3$.
using Polyhedra: polyhedron, vrep, volume
using LinearAlgebra
points = [-1.0 1.0 1.0
-1.0 -1.0 1.0
1.0 0.0 -1.0
1.0 1.0 1.0
0.0 1.0 -1.0]
function volume_simplex(points)
A = Matrix{T}(undef, 3, 3)
for i in 1:3
A[i, :] = points[i+1, :] - points[1, :]
end
return abs(det(A)) / 6
end
p = polyhedron(vrep(points))
println(volume(p)) # 5.3..., wrong
vol_1 = volume_simplex(points[[1, 2, 4, 5], :])
vol_2 = volume_simplex(points[[2, 4, 5, 3], :])
println("$vol_1 $vol_2 $(vol_1 + vol_2)") # 1.3..., 1.3..., 2.6..., as expected
<https://user-images.githubusercontent.com/21228024/227965164-656f2850-01fc-4df3-9d23-5907ce74f30b.png>
—
Reply to this email directly, view it on GitHub <#321 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAG2COO4WECZVV3577CCD4LW6GONNANCNFSM6AAAAAAWIIHOD4>.
You are receiving this because you commented.
|
|
I didn't run your code, but it looks like the last two faces in |
|
OK, with this triangulation:
julia> TW = [[1,2,4],[2,3,4],[3,5,4],[1,4,5],[2,1,5],[2,5,3]]
julia> vol = Lar.volume(pol)
2.666666666666666
you get the above.
… On Mar 27, 2023, at 7:44 PM, Benjamin Qi ***@***.***> wrote:
I didn't run your code, but it looks like the last two faces in TW should be 2,1,5 and 2,5,3.
—
Reply to this email directly, view it on GitHub <#321 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAG2CONWS6M35KHONMZKUQLW6HGYHANCNFSM6AAAAAAWIIHOD4>.
You are receiving this because you commented.
|
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
When using the default library (or
CDDLib), anAssertionErroris thrown when computing the volume of a 3D polyhedron where the same point is passed twice tovrep.Full Stacktrace
The text was updated successfully, but these errors were encountered: