2010-08-11 1 views

ответ

27

Для Matlab:

point = [1,2,3]; 
normal = [1,1,2]; 

%# a plane is a*x+b*y+c*z+d=0 
%# [a,b,c] is the normal. Thus, we have to calculate 
%# d and we're set 
d = -point*normal'; %'# dot product for less typing 

%# create x,y 
[xx,yy]=ndgrid(1:10,1:10); 

%# calculate corresponding z 
z = (-normal(1)*xx - normal(2)*yy - d)/normal(3); 

%# plot the surface 
figure 
surf(xx,yy,z) 

enter image description here

Примечание: это решение работает только до тех пор, как обычно (3) не равен 0. Если плоскость параллельна оси, вы можете вращать размеры, чтобы сохранить тот же подход:

z = (-normal(3)*xx - normal(1)*yy - d)/normal(2); %% assuming normal(3)==0 and normal(2)~=0 

%% plot the surface 
figure 
surf(xx,yy,z) 

%% label the axis to avoid confusion 
xlabel('z') 
ylabel('x') 
zlabel('y') 
+0

О ничего себе, я никогда не знал, что даже была функция ndgrid. Здесь я прыгал через обручи с repmat и индексированием, чтобы создать их на лету все это время ха-ха. Благодаря! ** Редактирование: ** btw было бы z = -нормальным (1) * xx - нормальным (2) * yy-d; вместо? – Xzhsh

+0

@ Xzhsh: ой, да. Исправлена. – Jonas

+0

также делятся на нормальные (3);). На всякий случай кто-то еще смотрит на этот вопрос и путается – Xzhsh

45

для всех копий/pasters там, здесь похож код Python с использованием Matplotlib:

import numpy as np 
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D 

point = np.array([1, 2, 3]) 
normal = np.array([1, 1, 2]) 

# a plane is a*x+b*y+c*z+d=0 
# [a,b,c] is the normal. Thus, we have to calculate 
# d and we're set 
d = -point.dot(normal) 

# create x,y 
xx, yy = np.meshgrid(range(10), range(10)) 

# calculate corresponding z 
z = (-normal[0] * xx - normal[1] * yy - d) * 1. /normal[2] 

# plot the surface 
plt3d = plt.figure().gca(projection='3d') 
plt3d.plot_surface(xx, yy, z) 
plt.show() 

enter image description here

+1

Обратите внимание, что 'z' имеет тип' int' в исходном фрагменте, который создает поверхность wiggly. Я бы использовал 'z = (-normal [0] * xx - normal [1] * yy - d) * 1./normal [2]', чтобы преобразовать z в 'real'. – Falcon

+1

Большое спасибо Falcon, перед вашим комментарием я на самом деле думал, что это ограничение с matplotlib. Я попытался компенсировать с помощью 100 элементов -> диапазон (100), в то время как в примере Matlab использовалось только 10 -> 1:10. Я отредактировал свое решение соответствующим образом. –

+0

Если вы хотите сделать вывод более сопоставимым с примером @Jonas matlab, выполните следующие действия: a) замените 'range (10)' на 'np.arange (1,11)'. b) добавьте строку 'plt3d.azim = -135.0' перед' plt.show() '(поскольку Matlab и matplotlib, похоже, имеют разные по умолчанию вращения). c) Nitpicking: 'xlim ([0,10])' и 'ylim ([0, 10])'. Наконец, добавление меток оси помогло бы увидеть основное различие в первую очередь, поэтому я добавил бы «xlabel ('x')» и «ylabel ('y')» для ясности и, соответственно, для примера Matlab. – Joma

4

Для копировальных pasters желающих градиент на поверхности:

from mpl_toolkits.mplot3d import Axes3D 
from matplotlib import cm 
import numpy as np 
import matplotlib.pyplot as plt 

point = np.array([1, 2, 3]) 
normal = np.array([1, 1, 2]) 

# a plane is a*x+b*y+c*z+d=0 
# [a,b,c] is the normal. Thus, we have to calculate 
# d and we're set 
d = -point.dot(normal) 

# create x,y 
xx, yy = np.meshgrid(range(10), range(10)) 

# calculate corresponding z 
z = (-normal[0] * xx - normal[1] * yy - d) * 1./normal[2] 

# plot the surface 
plt3d = plt.figure().gca(projection='3d') 

Gx, Gy = np.gradient(xx * yy) # gradients with respect to x and y 
G = (Gx ** 2 + Gy ** 2) ** .5 # gradient magnitude 
N = G/G.max() # normalize 0..1 

plt3d.plot_surface(xx, yy, z, rstride=1, cstride=1, 
        facecolors=cm.jet(N), 
        linewidth=0, antialiased=False, shade=False 
) 
plt.show() 

enter image description here

3

Приведенные выше ответы достаточно хороши. Следует упомянуть, что они используют тот же метод, который вычисляет значение z для данного (x, y). Обратный ход приходит к тому, что они пересекают плоскость, а плоскость в пространстве может меняться (только сохраняя ее проецирование одинаково). Например, вы не можете получить квадрат в 3D-пространстве (но искаженный).

Чтобы избежать этого, существует другой способ использования вращения. Если вы сначала генерируете данные в плоскости x-y (может быть любой формы), а затем поворачивайте ее на равную величину ([0 0 1] на ваш вектор), тогда вы получите то, что хотите. Просто запустите ниже код для справки.

point = [1,2,3]; 
normal = [1,2,2]; 
t=(0:10:360)'; 
circle0=[cosd(t) sind(t) zeros(length(t),1)]; 
r=vrrotvec2mat(vrrotvec([0 0 1],normal)); 
circle=circle0*r'+repmat(point,length(circle0),1); 
patch(circle(:,1),circle(:,2),circle(:,3),.5); 
axis square; grid on; 
%add line 
line=[point;point+normr(normal)] 
hold on;plot3(line(:,1),line(:,2),line(:,3),'LineWidth',5) 

Это получить круг в 3D: Resulting picture

0

уборщик пример Python, который также работает для хитрого $ Z, Y, Z $ ситуаций,

from mpl_toolkits.mplot3d import axes3d 
from matplotlib.patches import Circle, PathPatch 
import matplotlib.pyplot as plt 
from matplotlib.transforms import Affine2D 
from mpl_toolkits.mplot3d import art3d 
import numpy as np 

def plot_vector(fig, orig, v, color='blue'): 
    ax = fig.gca(projection='3d') 
    orig = np.array(orig); v=np.array(v) 
    ax.quiver(orig[0], orig[1], orig[2], v[0], v[1], v[2],color=color) 
    ax.set_xlim(0,10);ax.set_ylim(0,10);ax.set_zlim(0,10) 
    ax = fig.gca(projection='3d') 
    return fig 

def rotation_matrix(d): 
    sin_angle = np.linalg.norm(d) 
    if sin_angle == 0:return np.identity(3) 
    d /= sin_angle 
    eye = np.eye(3) 
    ddt = np.outer(d, d) 
    skew = np.array([[ 0, d[2], -d[1]], 
        [-d[2],  0, d[0]], 
        [d[1], -d[0], 0]], dtype=np.float64) 

    M = ddt + np.sqrt(1 - sin_angle**2) * (eye - ddt) + sin_angle * skew 
    return M 

def pathpatch_2d_to_3d(pathpatch, z, normal): 
    if type(normal) is str: #Translate strings to normal vectors 
     index = "xyz".index(normal) 
     normal = np.roll((1.0,0,0), index) 

    normal /= np.linalg.norm(normal) #Make sure the vector is normalised 
    path = pathpatch.get_path() #Get the path and the associated transform 
    trans = pathpatch.get_patch_transform() 

    path = trans.transform_path(path) #Apply the transform 

    pathpatch.__class__ = art3d.PathPatch3D #Change the class 
    pathpatch._code3d = path.codes #Copy the codes 
    pathpatch._facecolor3d = pathpatch.get_facecolor #Get the face color  

    verts = path.vertices #Get the vertices in 2D 

    d = np.cross(normal, (0, 0, 1)) #Obtain the rotation vector  
    M = rotation_matrix(d) #Get the rotation matrix 

    pathpatch._segment3d = np.array([np.dot(M, (x, y, 0)) + (0, 0, z) for x, y in verts]) 

def pathpatch_translate(pathpatch, delta): 
    pathpatch._segment3d += delta 

def plot_plane(ax, point, normal, size=10, color='y'):  
    p = Circle((0, 0), size, facecolor = color, alpha = .2) 
    ax.add_patch(p) 
    pathpatch_2d_to_3d(p, z=0, normal=normal) 
    pathpatch_translate(p, (point[0], point[1], point[2])) 


o = np.array([5,5,5]) 
v = np.array([3,3,3]) 
n = [0.5, 0.5, 0.5] 

from mpl_toolkits.mplot3d import Axes3D 
fig = plt.figure() 
ax = fig.gca(projection='3d') 
plot_plane(ax, o, n, size=3)  
ax.set_xlim(0,10);ax.set_ylim(0,10);ax.set_zlim(0,10) 
plt.show() 

enter image description here

Смежные вопросы