# How to plot a 3D plot of 3 inequalities?
**URL:** https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910
**Category:** General Usage
**Tags:** question, makie, polyhedra
**Created:** [February 2, 2023, 8:28am UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910 "2023-02-02T08:28:58Z")
**Posts on this page:** 11
**Page:** 1
### Author:  [@Kanishk](https://discourse.julialang.org/u/Kanishk)
#### Post date: [February 2, 2023, 8:28am UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/1 "2023-02-02T08:28:58Z")
Hello
I have 3 inequalities and each have 3 variables. Could you please explain how can I draw a combined one 3D plot with shaded area satisfying all the 3 inequalities? I have inequalities of the form as
```julia
f1(x,y,z) = a1*x + b1*y+ c1*z <= 0
f2(x,y,z) = a2*x + b2*y+ c2*z <= 0
f3(x,y,z) = a3*x + b3*y+ c3*z <= 0
```
Thank you
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [February 2, 2023, 10:33am UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/2 "2023-02-02T10:33:12Z")
A simple solution, but in the form of a 3d scatter plot is [given here](https://stackoverflow.com/questions/61732123/3d-plotting-region-fulfilling-inequalities-in-matplotlib).
---
### Author:  [@Kanishk](https://discourse.julialang.org/u/Kanishk)
#### Post date: [February 2, 2023, 11:12am UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/3 "2023-02-02T11:12:53Z")
Dear @rafael.guerra,
Thank you for providing a solution. But I am trying to plot using “IntervalArithmetic.jl” package.
I tried to plot using the code shown below-
```julia
f1(p, q, r) = p^2 + q^2 + r^2 - 4 #sphere
f2(p, q, r) = p^2 + q^2 + r^2 - 16 #sphere
# f3(p, q, r) = 2*p + 3*q
zrange = -5:.1:5
plt = plot(palette = palette(cgrad([:red,:green,:blue],length(zrange))),camera=(45,45),xlims=(-7,7), ylims=(-7,7), xlabel="x",ylabel="y",zlabel="z")
[plot!(plt,((p, q)->f1(p, q, r)) ≧ 0, zpos=r) for r=zrange]
[plot!(plt,((p, q)->f2(p, q, r)) ≦ 0, zpos=r) for r=zrange]
# [plot!(plt,((p, q)->f3(p, q, r)) ⩵ 0, zpos=r) for r=zrange]
plt
```
I am getting this plot-
But, don’t know why these additional square boxes and and lines are coming. I just want spherical region in the plot. Could you please help how can I remove everything other than the spherical region.
Thank you
---
### Author:  [@empet](https://discourse.julialang.org/u/empet)
#### Post date: [February 2, 2023, 9:09pm UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/4 "2023-02-02T21:09:59Z")
Your region of interest is delimited by three planes, passing all through the origin. Hence you should plot the three planes as parameterized surfaces. I don’t know how this could be plotted with Plots, but I explain with a PlotlyJS example.
Suppose that the three planes are defined by the linear functions:
```julia
f(x,y,z) =2x-y+3z
g(x,y,z)=-x+y+z
h(x,y,z) =x-3y-5z
```
```julia
using LinearAlgebra, PlotlyJS
import StaticArrays: SMatrix
function planeparam(normal::SMatrix{1, 3, T}; p=[0,0,0]) where T<:AbstractFloat
~iszero(normal) || error("normal vector is null vect")
M = nullspace(normal)
a, b = eachcol(M)
#parameterization of the plane π(p; normal=no) defined by the point p and the given normal
(u,v)->[a[k]*u+b[k]*v+p[k] for k in 1:3] # theoretically: (u, v)->[X(u, v), Y(u,v), Z(u,v)]
end
#Get the parameterization of each plane
F = planeparam(SMatrix{1,3}([2.0 -1.0 3]); )
G = planeparam(SMatrix{1,3}([-1.0 1 1.]); )
H = planeparam(SMatrix{1,3}([1.0 -3 -5.]); )
#define a meshgrid over a rectangular region in the parametric plane (u,v):
n=200
ul = vl= range(-3, 3, n)
u = ones(n) * ul'
v = vl *ones(n)'
fig= Plot(Layout(width=450, height=450, font_size=11,
coloraxis = attr(colorscale=colors.algae, reversescale=true,
showscale=false ),
scene_camera_eye=attr(x=1.5, y=1.5, z=1)))
#add the planes as parameterized surfaces to the figure fig
for K in [F, G, H]
X, Y, Z = eachrow(hcat(K.(u,v)...))
addtraces!(fig, surface( x=reshape(X, n, n), y=reshape(Y, n, n),
z=reshape(Z, n, n), opacity=0.75, coloraxis="coloraxis",
))
end
display(fig)
```
I selected three random planar surfaces, hence the plot is not quite convincing, but it can suggest an idea on how to proceed. Here is missing the condition \<=0. For linear functions it is very simple to get the region whree both three are negative, but it is complicated to shade that volume, because it is covered by the surfaces.
Note that depending on your linear functions, the volume where all three are negative can be bounded or unbounded.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [February 3, 2023, 9:47am UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/5 "2023-02-03T09:47:28Z")
This problem seems to correspond exactly to [Polyhedra.jl’s H-representation of halfspaces](https://juliapolyhedra.github.io/Polyhedra.jl/stable/representation/#H-representation).
We can use `Polyhedra.jl` and `Makie.jl` to plot it:
**NB:**
The result is not as expected, as per @sijo’s post below. Needs more work.
_Rome ne s’est pas faite en un jour…_
> **Polyhedra.jl and Makie.jl code**
>
> ```julia
> using Polyhedra
> using Makie, GLMakie
>
> const a1, b1, c1 = 0., 1., 1.
> const a2, b2, c2 = 1., 1., 1.
> const a3, b3, c3 = 1., -1., 1.
>
> f1(x,y,z) = a1*x + b1*y + c1*z
> f2(x,y,z) = a2*x + b2*y + c2*z
> f3(x,y,z) = a3*x + b3*y + c3*z
>
> h = HalfSpace([a1, b1, c1], 0) ∩ HalfSpace([a2, b2, c2], 0) ∩ HalfSpace([a3, b3, c3], 0)
> p = polyhedron(h)
> m = Polyhedra.Mesh(p)
>
> Makie.mesh(m, color=:blue)
>
> ```
---
### Author:  [@empet](https://discourse.julialang.org/u/empet)
#### Post date: [February 3, 2023, 12:30pm UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/6 "2023-02-03T12:30:46Z")
It’s a good idea to use Polyhedra.jl when the solution is a bounded volume, but if it unbounded it throws the error:
```julia
hred: Redundancy Polyhedra.UNKNOWN_REDUNDANCY
vrep: Nothing nothing
vred: Redundancy Polyhedra.UNKNOWN_REDUNDANCY
solver: Nothing nothing
```
But maybe for OP’s case it is bounded.
---
### Author:  [@sijo](https://discourse.julialang.org/u/sijo)
#### Post date: [February 3, 2023, 12:45pm UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/7 "2023-02-03T12:45:57Z")
I don’t know what Polyhedra.jl is doing in @rafael.guerra’s example but you can’t have a (nonzero) bounded volume delimited by three planes…
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [February 3, 2023, 1:21pm UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/8 "2023-02-03T13:21:07Z")
@sijo, good observation and I do not know the answer either.
Using the brute force scatter3d method, as per the Python post linked previously, we get an unbounded volume (computed over limited range here):
> **Plots.jl plotly() code**
>
> ```julia
> const a1, b1, c1 = 0., 1., 1.
> const a2, b2, c2 = 1., 1., 1.
> const a3, b3, c3 = 1., -1., 1.
>
> f1(x,y,z) = a1*x + b1*y + c1*z
> f2(x,y,z) = a2*x + b2*y + c2*z
> f3(x,y,z) = a3*x + b3*y + c3*z
>
> sup = -10:0.2:10
> X = Float32[]
> Y = Float32[]
> Z = Float32[]
> for x in sup, y in sup, z in sup
> if f1(x,y,z) ≤ 0 && f2(x,y,z) ≤ 0 && f3(x,y,z) ≤ 0
> push!(X, x)
> push!(Y, y)
> push!(Z, z)
> end
> end
>
> using Plots; plotly(size=(600,600))
> scatter3d(X,Y,Z, ms=1, msw=0.1, xlabel="x", ylabel="y", zlabel="z", c=:blues)
>
> ```
---
### Author:  [@Kanishk](https://discourse.julialang.org/u/Kanishk)
#### Post date: [February 3, 2023, 1:40pm UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/9 "2023-02-03T13:40:35Z")
Dear @rafael.guerra,
As I don’t know whether the actual volume will be bounded or unbounded, this type of solution (plot) seems useful to me. Could you please share julia code and julia packages used for this?
Thank you
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [February 3, 2023, 1:47pm UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/10 "2023-02-03T13:47:50Z")
Included it below the plot.
---
### Author:  [@empet](https://discourse.julialang.org/u/empet)
#### Post date: [February 4, 2023, 12:32pm UTC](https://discourse.julialang.org/t/how-to-plot-a-3d-plot-of-3-inequalities/93910/11 "2023-02-04T12:32:25Z")
Update: Such a system of linear inequalities can be solved with `MarchingCubes.jl` as follows:
1. define a vector valued function f(x,y, z) =[f1(x,y,z), f2(x,y,z), f3(x,y,z)],
and a cube with center at origin, \text{Cube}= [-a, a]^3
2. a grid on this cube and a function that gets 0 at the boundary of the region of interest, i.e. the volume V=\{(x,y,z)\:|\: f1(x,y, z)\<=0, f2(x,y, z)\<=0 , f3(x,y, z)\<=0 \}. Evaluate this function at the grid points, to get `fvals`
3. Call the marching cube method to retrieve a triangulation of the volume boundary and plot the corresponding mesh.
```julia
using MarchingCubes
f(x,y,z) = [y + z, x+y+z, x-y+z]
#f(x,y, z) = [2*x-4.5*y+3.2*z, -3x+2.1*y+1.0*z, x-3*y-5*z]
const Fl=Float64
#Define a cube of side 14
xl = yl = zl = collect(Float64,range(-7.0, 7.0, length=150))
x = [xx for zz in zl, yy in yl, xx in xl]
y = [yy for zz in zl, yy in yl, xx in xl]
z = [zz for zz in zl, yy in yl, xx in xl];
fvals = [maximum([f(xx, yy, zz)[k] for k=1:3]) for zz in zl, yy in yl, xx in xl]
#instantiate the structure MC (Marching Cube)
mc = MC(fvals, Int; x=xl, y=yl, z=zl)
m = march(mc, 0.0) #0.0 is the isovalue
```
1. Transform `mc.vertices` and `mc.triangles` to data types necessary to plot a 3d mesh in your favorite graphics library.
This is the plot I created with PlotlyJS.mesh3d: [https://github.com/empet/Datasets/blob/master/Images/sol-linineq1.png](https://github.com/empet/Datasets/blob/master/Images/sol-linineq1.png)