ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 18 Jun 2012 07:16:12 +0200Forming a solid from a set of intersecting planeshttps://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/suppose I have equations of planes in 3D space in the form:
$A_ix + B_iy + C_iz = D_i$ for some values of i
**how can use sage to find the solid enclosed by them?** I need a plot of the surface of this solid. Actually I have many such sets of planes each forming a solid like this.(all concentric)
**For example:**
planes at x = -1, 1 in yz plane;
planes at y = -1, 1 in xz plane and planes at z = -1, 1 in xy plane intersect and form a unit cube centred at the origin.
ty ndomes!
But, its taking very long. I have like 12 or more equations (seen a maximum of 30 but can be more). I am now thinking of a different approach. We can easily make a function:
def inside(point , planes):
if the point is inside all the planes:
return 1
else:
return 0
so, now is there a good function in sagemath that can show all the region in space for which f(point) gives 1. (yea, its okay if the solids are filled as long as they're transparent enough. . Naively Taking lots of points and finding this solid(Like monte carlo) is very inefficient
Actually the program is taking very long time for these equations.
[-x - y - 1 == 0, -x - z - 1 == 0, -x + z - 1 == 0, -x + y - 1 == 0, -y- z - 1== 0, -y + z - 1 == 0, y - z - 1 == 0, y + z - 1 == 0, x - y - 1== 0, x - z - 1 == 0, x + z - 1 == 0, x + y - 1 == 0]
so, check it once. Thanks!Tue, 12 Jun 2012 11:26:42 +0200https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/Answer by ndomes for <p>suppose I have equations of planes in 3D space in the form:</p>
<p>$A_ix + B_iy + C_iz = D_i$ for some values of i</p>
<p><strong>how can use sage to find the solid enclosed by them?</strong> I need a plot of the surface of this solid. Actually I have many such sets of planes each forming a solid like this.(all concentric)</p>
<p><strong>For example:</strong></p>
<p>planes at x = -1, 1 in yz plane;
planes at y = -1, 1 in xz plane and planes at z = -1, 1 in xy plane intersect and form a unit cube centred at the origin.</p>
<p>ty ndomes!
But, its taking very long. I have like 12 or more equations (seen a maximum of 30 but can be more). I am now thinking of a different approach. We can easily make a function:</p>
<pre><code>def inside(point , planes):
if the point is inside all the planes:
return 1
else:
return 0
</code></pre>
<p>so, now is there a good function in sagemath that can show all the region in space for which f(point) gives 1. (yea, its okay if the solids are filled as long as they're transparent enough. . Naively Taking lots of points and finding this solid(Like monte carlo) is very inefficient</p>
<p>Actually the program is taking very long time for these equations. </p>
<pre><code>[-x - y - 1 == 0, -x - z - 1 == 0, -x + z - 1 == 0, -x + y - 1 == 0, -y- z - 1== 0, -y + z - 1 == 0, y - z - 1 == 0, y + z - 1 == 0, x - y - 1== 0, x - z - 1 == 0, x + z - 1 == 0, x + y - 1 == 0]
</code></pre>
<p>so, check it once. Thanks!</p>
https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?answer=13711#post-id-13711If you write an equation as list [A,B,C,D] instead of A*x + B*y + C*z == D,
you can easily use matrix operations to find the vertices.
This seems to be much faster.
equations =[[1,0,0,1],[1,0,0,-1],[0,1,0,1],[0,1,0,-1],[0,0,1,1],[0,0,1,-1]]
def get_matrix(L,K):
M = matrix(QQ,3,3)
M.set_row(0,L[K[0]][:3])
M.set_row(1,L[K[1]][:3])
M.set_row(2,L[K[2]][:3])
return M
def get_vector(L,K):
return vector(QQ,(L[K[0]][3],L[K[1]][3],L[K[2]][3]))
def find_vertices_2(eqns):
vertices = []
for C in combinations(range(len(eqns)),3):
M = get_matrix(eqns,C)
b = get_vector(eqns,C)
if M.is_invertible():
P = M.solve_right(b)
if not P in vertices:
vertices.append(P)
return vertices
I compared this function *find_vertices_2* to a function *find_vertices_1* based on solving systems of linear equations:
timeit('find_vertices_1(eqns)',number=20)
20 loops, best of 3: 1.3 s per loop
timeit('find_vertices_2(equations)',number=20)
20 loops, best of 3: 15.8 ms per loop
Sat, 16 Jun 2012 06:51:27 +0200https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?answer=13711#post-id-13711Comment by ebs for <p>If you write an equation as list [A,B,C,D] instead of A<em>x + B</em>y + C*z == D,
you can easily use matrix operations to find the vertices.
This seems to be much faster.</p>
<pre><code>equations =[[1,0,0,1],[1,0,0,-1],[0,1,0,1],[0,1,0,-1],[0,0,1,1],[0,0,1,-1]]
def get_matrix(L,K):
M = matrix(QQ,3,3)
M.set_row(0,L[K[0]][:3])
M.set_row(1,L[K[1]][:3])
M.set_row(2,L[K[2]][:3])
return M
def get_vector(L,K):
return vector(QQ,(L[K[0]][3],L[K[1]][3],L[K[2]][3]))
def find_vertices_2(eqns):
vertices = []
for C in combinations(range(len(eqns)),3):
M = get_matrix(eqns,C)
b = get_vector(eqns,C)
if M.is_invertible():
P = M.solve_right(b)
if not P in vertices:
vertices.append(P)
return vertices
</code></pre>
<p>I compared this function <em>find_vertices_2</em> to a function <em>find_vertices_1</em> based on solving systems of linear equations:</p>
<pre><code>timeit('find_vertices_1(eqns)',number=20)
</code></pre>
<p>20 loops, best of 3: 1.3 s per loop</p>
<pre><code>timeit('find_vertices_2(equations)',number=20)
</code></pre>
<p>20 loops, best of 3: 15.8 ms per loop</p>
https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19577#post-id-19577ok I"m marking this as solved because I think I can complete the code for sorting the vertices and joining them in orderMon, 18 Jun 2012 07:16:12 +0200https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19577#post-id-19577Comment by ndomes for <p>If you write an equation as list [A,B,C,D] instead of A<em>x + B</em>y + C*z == D,
you can easily use matrix operations to find the vertices.
This seems to be much faster.</p>
<pre><code>equations =[[1,0,0,1],[1,0,0,-1],[0,1,0,1],[0,1,0,-1],[0,0,1,1],[0,0,1,-1]]
def get_matrix(L,K):
M = matrix(QQ,3,3)
M.set_row(0,L[K[0]][:3])
M.set_row(1,L[K[1]][:3])
M.set_row(2,L[K[2]][:3])
return M
def get_vector(L,K):
return vector(QQ,(L[K[0]][3],L[K[1]][3],L[K[2]][3]))
def find_vertices_2(eqns):
vertices = []
for C in combinations(range(len(eqns)),3):
M = get_matrix(eqns,C)
b = get_vector(eqns,C)
if M.is_invertible():
P = M.solve_right(b)
if not P in vertices:
vertices.append(P)
return vertices
</code></pre>
<p>I compared this function <em>find_vertices_2</em> to a function <em>find_vertices_1</em> based on solving systems of linear equations:</p>
<pre><code>timeit('find_vertices_1(eqns)',number=20)
</code></pre>
<p>20 loops, best of 3: 1.3 s per loop</p>
<pre><code>timeit('find_vertices_2(equations)',number=20)
</code></pre>
<p>20 loops, best of 3: 15.8 ms per loop</p>
https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19586#post-id-19586I took for granted that all faces are convex. May be I misunderstood you. There is still your idea to order the vertices of each face by the respective rotation angle.Sat, 16 Jun 2012 09:10:10 +0200https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19586#post-id-19586Comment by ebs for <p>If you write an equation as list [A,B,C,D] instead of A<em>x + B</em>y + C*z == D,
you can easily use matrix operations to find the vertices.
This seems to be much faster.</p>
<pre><code>equations =[[1,0,0,1],[1,0,0,-1],[0,1,0,1],[0,1,0,-1],[0,0,1,1],[0,0,1,-1]]
def get_matrix(L,K):
M = matrix(QQ,3,3)
M.set_row(0,L[K[0]][:3])
M.set_row(1,L[K[1]][:3])
M.set_row(2,L[K[2]][:3])
return M
def get_vector(L,K):
return vector(QQ,(L[K[0]][3],L[K[1]][3],L[K[2]][3]))
def find_vertices_2(eqns):
vertices = []
for C in combinations(range(len(eqns)),3):
M = get_matrix(eqns,C)
b = get_vector(eqns,C)
if M.is_invertible():
P = M.solve_right(b)
if not P in vertices:
vertices.append(P)
return vertices
</code></pre>
<p>I compared this function <em>find_vertices_2</em> to a function <em>find_vertices_1</em> based on solving systems of linear equations:</p>
<pre><code>timeit('find_vertices_1(eqns)',number=20)
</code></pre>
<p>20 loops, best of 3: 1.3 s per loop</p>
<pre><code>timeit('find_vertices_2(equations)',number=20)
</code></pre>
<p>20 loops, best of 3: 15.8 ms per loop</p>
https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19587#post-id-19587Does your code go into an infinite loop if the polygon is concave? I think my planes can give rise to concave polygons. (Actually I am trying to get the brillouin zones of a 3d lattice). See this for example: http://www.wolfram.com/products/publicon/samples/default/HTMLFiles/DefaultSample_90.gifSat, 16 Jun 2012 08:13:51 +0200https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19587#post-id-19587Answer by ndomes for <p>suppose I have equations of planes in 3D space in the form:</p>
<p>$A_ix + B_iy + C_iz = D_i$ for some values of i</p>
<p><strong>how can use sage to find the solid enclosed by them?</strong> I need a plot of the surface of this solid. Actually I have many such sets of planes each forming a solid like this.(all concentric)</p>
<p><strong>For example:</strong></p>
<p>planes at x = -1, 1 in yz plane;
planes at y = -1, 1 in xz plane and planes at z = -1, 1 in xy plane intersect and form a unit cube centred at the origin.</p>
<p>ty ndomes!
But, its taking very long. I have like 12 or more equations (seen a maximum of 30 but can be more). I am now thinking of a different approach. We can easily make a function:</p>
<pre><code>def inside(point , planes):
if the point is inside all the planes:
return 1
else:
return 0
</code></pre>
<p>so, now is there a good function in sagemath that can show all the region in space for which f(point) gives 1. (yea, its okay if the solids are filled as long as they're transparent enough. . Naively Taking lots of points and finding this solid(Like monte carlo) is very inefficient</p>
<p>Actually the program is taking very long time for these equations. </p>
<pre><code>[-x - y - 1 == 0, -x - z - 1 == 0, -x + z - 1 == 0, -x + y - 1 == 0, -y- z - 1== 0, -y + z - 1 == 0, y - z - 1 == 0, y + z - 1 == 0, x - y - 1== 0, x - z - 1 == 0, x + z - 1 == 0, x + y - 1 == 0]
</code></pre>
<p>so, check it once. Thanks!</p>
https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?answer=13694#post-id-13694 var("x y z")
P = vector((x,y,z))
eqns = [ x==-1, x==1, y==-1, y==1, z==-1, z==1 ]
#eqns = [ z==0, x+y+z ==1, -x+y+z == 1, -x-y+z==1, x-y+z==1]
indizes = range(len(eqns))
vertices = []
for C in combinations(indizes,3):
lse = [eqns[i] for i in C]
sol = solve(lse,x,y,z,solution_dict=True)
if sol:
vertices.append(P.substitute(sol[0]))
#for V in vertices: print V
"draw the solid - start with vertices"
G = point3d(vertices,size=10,color='red')
" faces "
for eq in eqns:
poly = []
for V in vertices:
if(eq.substitute({x:V[0],y:V[1],z:V[2]})):
if not V in poly:
poly.append(V)
"""TODO poly should be convex"""
shuffle(poly)
# print poly
try:
G += polygon(poly,alpha=0.5,color=(random(),random(),random()))
except:
pass
G.show(aspect_ratio=1)Tue, 12 Jun 2012 14:54:00 +0200https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?answer=13694#post-id-13694Comment by ebs for <pre><code>var("x y z")
P = vector((x,y,z))
eqns = [ x==-1, x==1, y==-1, y==1, z==-1, z==1 ]
#eqns = [ z==0, x+y+z ==1, -x+y+z == 1, -x-y+z==1, x-y+z==1]
indizes = range(len(eqns))
vertices = []
for C in combinations(indizes,3):
lse = [eqns[i] for i in C]
sol = solve(lse,x,y,z,solution_dict=True)
if sol:
vertices.append(P.substitute(sol[0]))
#for V in vertices: print V
"draw the solid - start with vertices"
G = point3d(vertices,size=10,color='red')
" faces "
for eq in eqns:
poly = []
for V in vertices:
if(eq.substitute({x:V[0],y:V[1],z:V[2]})):
if not V in poly:
poly.append(V)
"""TODO poly should be convex"""
shuffle(poly)
# print poly
try:
G += polygon(poly,alpha=0.5,color=(random(),random(),random()))
except:
pass
G.show(aspect_ratio=1)
</code></pre>
https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19607#post-id-19607I think we can find all vertices on each face and sort them according to polar angle about the centroid of the face. Then we can plot each face separately.
Wed, 13 Jun 2012 00:44:36 +0200https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19607#post-id-19607Comment by ndomes for <pre><code>var("x y z")
P = vector((x,y,z))
eqns = [ x==-1, x==1, y==-1, y==1, z==-1, z==1 ]
#eqns = [ z==0, x+y+z ==1, -x+y+z == 1, -x-y+z==1, x-y+z==1]
indizes = range(len(eqns))
vertices = []
for C in combinations(indizes,3):
lse = [eqns[i] for i in C]
sol = solve(lse,x,y,z,solution_dict=True)
if sol:
vertices.append(P.substitute(sol[0]))
#for V in vertices: print V
"draw the solid - start with vertices"
G = point3d(vertices,size=10,color='red')
" faces "
for eq in eqns:
poly = []
for V in vertices:
if(eq.substitute({x:V[0],y:V[1],z:V[2]})):
if not V in poly:
poly.append(V)
"""TODO poly should be convex"""
shuffle(poly)
# print poly
try:
G += polygon(poly,alpha=0.5,color=(random(),random(),random()))
except:
pass
G.show(aspect_ratio=1)
</code></pre>
https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19592#post-id-19592I found another way to reorder the polygons.
See http://www.hr.shuttle.de:9000/home/pub/104/Thu, 14 Jun 2012 10:22:04 +0200https://ask.sagemath.org/question/9062/forming-a-solid-from-a-set-of-intersecting-planes/?comment=19592#post-id-19592