3D volume fails on tiny point sets

[bqi343]

When using the default library (or CDDLib), an AssertionError is thrown when computing the volume of a 3D polyhedron where the same point is passed twice to vrep.

using Polyhedra, CDDLib, QHull

const libs = [Polyhedra.default_library(3, Float64), CDDLib.Library(:float), QHull.Library()]

function get_volume(point_set::Matrix{T}, lib) where {T}
    poly = polyhedron(vrep(point_set), lib)
    Polyhedra.volume(poly)
end

@testset "test_volume_tricky" begin
    points = [-1.0 0.0 1.0
        0.0 0.0 1.0
        0.0 0.0 0.0
        -1.0 -1.0 0.0
        -1.0 -1.0 0.0]
    for lib in libs
        # same error for default library or CDDLib
        # expected answer (0.16666666666666666) for qhull
        println(lib, " ", get_volume(points, lib))
    end
end

Full Stacktrace

test_volume_tricky: Error During Test at /Users/benq/.julia/packages/ReTest/WnRIG/src/ReTest.jl:517
  Got exception outside of a @test
  AssertionError: codim >= 0
  Stacktrace:
    [1] _triangulation(Δs::Vector{Vector{Polyhedra.Index{Float64, Vector{Float64}}}}, Δ::Vector{Polyhedra.Index{Float64, Vector{Float64}}}, v_idx::Vector{Polyhedra.Index{Float64, Vector{Float64}}}, h_idx::Vector{Polyhedra.Index{Float64, HalfSpace{Float64, Vector{Float64}}}}, incident_idx::Dict{Polyhedra.Index{Float64, HalfSpace{Float64, Vector{Float64}}}, Set{Polyhedra.Index{Float64, Vector{Float64}}}}, is_weak_adjacent::Dict{Tuple{Polyhedra.Index{Float64, HalfSpace{Float64, Vector{Float64}}}, Polyhedra.Index{Float64, HalfSpace{Float64, Vector{Float64}}}}, Bool}, codim::Int64)
      @ Polyhedra ~/.julia/packages/Polyhedra/Xhcfx/src/triangulation.jl:10
    [2] _triangulation(Δs::Vector{Vector{Polyhedra.Index{Float64, Vector{Float64}}}}, Δ::Vector{Polyhedra.Index{Float64, Vector{Float64}}}, v_idx::Vector{Polyhedra.Index{Float64, Vector{Float64}}}, h_idx::Vector{Polyhedra.Index{Float64, HalfSpace{Float64, Vector{Float64}}}}, incident_idx::Dict{Polyhedra.Index{Float64, HalfSpace{Float64, Vector{Float64}}}, Set{Polyhedra.Index{Float64, Vector{Float64}}}}, is_weak_adjacent::Dict{Tuple{Polyhedra.Index{Float64, HalfSpace{Float64, Vector{Float64}}}, Polyhedra.Index{Float64, HalfSpace{Float64, Vector{Float64}}}}, Bool}, codim::Int64) (repeats 4 times)
      @ Polyhedra ~/.julia/packages/Polyhedra/Xhcfx/src/triangulation.jl:34
    [3] triangulation_indices(p::DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}})
      @ Polyhedra ~/.julia/packages/Polyhedra/Xhcfx/src/triangulation.jl:46
    [4] triangulation
      @ ~/.julia/packages/Polyhedra/Xhcfx/src/triangulation.jl:50 [inlined]
    [5] volume(p::DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}})
      @ Polyhedra ~/.julia/packages/Polyhedra/Xhcfx/src/polyhedron.jl:100
    [6] get_volume(point_set::Matrix{Float64}, lib::DefaultLibrary{Float64})
      @ Main.InnerProductMaxTests ~/Dev/InnerProductMax/test/compare_chull.jl:9
    [7] macro expansion
      @ ~/Dev/InnerProductMax/test/compare_chull.jl:15 [inlined]
    [8] macro expansion
      @ ~/.julia/packages/ReTest/WnRIG/src/testset.jl:654 [inlined]
    [9] macro expansion
      @ ~/Dev/InnerProductMax/test/compare_chull.jl:13 [inlined]
   [10] top-level scope
      @ ~/.julia/packages/ReTest/WnRIG/src/ReTest.jl:517
   [11] eval
      @ ./boot.jl:368 [inlined]
   [12] (::ReTest.var"#80#98"{ReTest.Testset.Format})(mod::Module, ts::ReTest.TestsetExpr, pat::ReTest.And, chan::NamedTuple{(:out, :compute, :preview), Tuple{Distributed.RemoteChannel{Channel{Union{Nothing, ReTest.Testset.ReTestSet}}}, Channel{Nothing}, Nothing}})
      @ ReTest ~/.julia/packages/ReTest/WnRIG/src/ReTest.jl:1160
   [13] #invokelatest#2
      @ ./essentials.jl:729 [inlined]
   [14] invokelatest
      @ ./essentials.jl:726 [inlined]
   [15] #153
      @ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Distributed/src/remotecall.jl:425 [inlined]
   [16] run_work_thunk(thunk::Distributed.var"#153#154"{ReTest.var"#80#98"{ReTest.Testset.Format}, Tuple{Module, ReTest.TestsetExpr, ReTest.And, NamedTuple{(:out, :compute, :preview), Tuple{Distributed.RemoteChannel{Channel{Union{Nothing, ReTest.Testset.ReTestSet}}}, Channel{Nothing}, Nothing}}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, print_error::Bool)
      @ Distributed /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Distributed/src/process_messages.jl:70
   [17] remotecall_fetch(::Function, ::Distributed.LocalProcess, ::Module, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
      @ Distributed /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Distributed/src/remotecall.jl:450
   [18] remotecall_fetch(::Function, ::Distributed.LocalProcess, ::Module, ::Vararg{Any})
      @ Distributed /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Distributed/src/remotecall.jl:449
   [19] remotecall_fetch(::Function, ::Int64, ::Module, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
      @ Distributed /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Distributed/src/remotecall.jl:492
   [20] remotecall_fetch
      @ /Applications/Julia-1.8.app/Contents/Resources/julia/share/julia/stdlib/v1.8/Distributed/src/remotecall.jl:492 [inlined]
   [21] macro expansion
      @ ~/.julia/packages/ReTest/WnRIG/src/ReTest.jl:1157 [inlined]
   [22] (::ReTest.var"#78#96"{Int64, ReTest.Testset.Format, Nothing, ReTest.Testset.ReTestSet, Base.Threads.Atomic{Bool}, Channel{Nothing}, Distributed.RemoteChannel{Channel{Union{Nothing, ReTest.Testset.ReTestSet}}}, Vector{Any}, ReTest.And, Module})()
      @ ReTest ./task.jl:484

[bqi343]

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

[apaoluzzi]

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 https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/321#issuecomment-1484161652, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAG2COMFH7NMT3AK7ZTC26TW6B2QBANCNFSM6AAAAAAWIIHOD4. You are receiving this because you are subscribed to this thread.

[bqi343]

I don't see any included visualization, and I'm pretty sure 8/4 is wrong considering that

1. QHull agrees with 8/3, and presumably (?) it's been tested extensively

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

2. 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

[apaoluzzi]

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 https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/321#issuecomment-1485189902, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAG2COO4WECZVV3577CCD4LW6GONNANCNFSM6AAAAAAWIIHOD4. You are receiving this because you commented.

[bqi343]

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.

[apaoluzzi]

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 https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/321#issuecomment-1485567933, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAG2CONWS6M35KHONMZKUQLW6HGYHANCNFSM6AAAAAAWIIHOD4. You are receiving this because you commented.