Intersection of non-overlapping point clouds with Polyhedra.jl

General Usage
, , ,
1

Hi all,

I need to identify whether two point “clouds” of arbitrary dimension intersect. As a sanity check, I tested two non-overlapping point clouds to see if the results make sense. Strangely, the intersection is non-empty. Am I doing something incorrectly?

using Polyhedra, Random

# ranges from 0 to 1 in each dimension
points1 = rand(1000, 2)
# ranges from 10 to 11 in each dimension
points2 = rand(1000, 2) .+ 10

ph1 = polyhedron(vrep(points1))

ph2 = polyhedron(vrep(points2))

intersect(ph1, ph2)
2

On my machine the intersection is an empty polyhedron, using Polyhedra v0.6.17.

Which version of Polyhedra are you using? Maybe updating the package will fix the problem.

System info:

Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 11th Gen Intel(R) Core(TM) i7-11800H @ 2.30GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, tigerlake)
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 8
3

Hi Brian. Thanks for your reply. I am also using the same version:

Version

(@v1.7) pkg> st Polyhedra
      Status `~/.julia/environments/v1.7/Project.toml`
  [67491407] Polyhedra v0.6.17

julia> versioninfo()
Julia Version 1.7.1
Commit ac5cc99908 (2021-12-22 19:35 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-12.0.1 (ORCJIT, haswell)
Environment:
  JULIA_CMDSTAN_HOME = /home/dfish/cmdstan
  JULIA_NUM_THREADS = 4
  JULIA_EDITOR = code

Code

using Polyhedra, Random

Random.seed!(5448)

points1 = rand(1000, 2)

points2 = rand(1000, 2) .+ 10

ph1 = polyhedron(vrep(points1))

ph2 = polyhedron(vrep(points2))

intersect(ph1, ph2)

When I run the code above (with the RNG seed set this time), I get the following result:

Polyhedron DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}:
38-element iterator of HalfSpace{Float64, Vector{Float64}}:
 HalfSpace([-0.7462477519569687, -0.03845547811729231], -0.013480793553864966)
 HalfSpace([-0.631648303668281, -0.07925921671122552], -0.018694015873781493)
 HalfSpace([-0.6562091647821909, 0.2518118199560842], 0.21274256968080027)
 HalfSpace([-0.7183230037709791, 0.19596921461056738], 0.15986258590309105)
 HalfSpace([-0.003954717247456662, -0.23874782378908826], -0.0015132839275770697)
 HalfSpace([-0.0010517702816531934, -0.23877617777670027], -0.0009270320917026878)
 HalfSpace([-0.006522312712251365, 0.5979383834248235], 0.5958076163402416)
 HalfSpace([0.01792971747635548, 0.5615940039679551], 0.5749463593095993)
 HalfSpace([0.20507167125043874, 0.03585744534363632], 0.23837882028596058)
 HalfSpace([-0.07428990598987983, -0.2256683809307513], -0.006862438464505046)
 HalfSpace([0.1949577049237804, -0.025710282476602772], 0.19227225404479864)
 HalfSpace([0.21621597234019985, 0.0016964974125853165], 0.21629736074799183)
 HalfSpace([0.01574311843644082, -0.2265573697244256], 0.011312917854648953)
 HalfSpace([0.17023892872160043, -0.05653372519373371], 0.1650827541793411)
 HalfSpace([0.00013767154323721662, 0.5965297173714741], 0.5961925181340423)
 HalfSpace([-0.8213607519521141, -0.004166129056533232], -0.004672555440764067)
 HalfSpace([0.17893288766097307, -0.0457972346051876], 0.1745709824015048)

Are you saying that you are getting an empty array instead?

4

I’m guessing they do lazy evaluation of the polyhedron. If you do this:

julia> points(intersect(ph1,ph2))
0-element iterator of Vector{Float64}

the result looks like an empty polyhedron as opposed to the long list of points that you see returned from intersect(ph1,ph2).

I’m not an expert on this package though. Maybe someone who knows it better can verify if this is true.

1 Like
5

Thank you. As far as I can tell, using points as you suggested seems to be giving a more reasonable result.It looks like the points correspond to the boundary of the intersection. Nonetheless, I would be interested in additional input from someone who is more familiar with this package.