Python生成简单3D管道
笔者最近研究了一下通过python生成简单3D管道的方法。主要通过空间坐标旋转平移的方法，得到管道起止点的圆截面（正多边形顶点——其法线向量与起止点向量相同），根据顶点构建管道面。并保存为obj文件。obj文件仅包含顶点v、组g、面f。
matplotlib,mpl_toolkits仅用于测试查看管道截面圆顶点的坐标生成情况。
代码如下：
import os, sys
from math import cos, sin, sqrt, atan2, pi, acos
from matplotlib import colors
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import psycopg2
import config
import pghelper
# 获取绕Z轴旋转矩阵
def GetRotateMatrix_AngleZ(az):
raz = np.array([[cos(az), -sin(az), 0.000, 0.000],
[sin(az), cos(az), 0.000, 0.000],
[0.0000, 0.0000, 1.0000, 0.0000],
[0.0000, 0.0000, 0.0000, 1.0000]])
return raz
# 获取绕X轴旋转矩阵
def GetRotateMatrix_AngleX(ax):
rax = np.array([[1.0000, 0.0000, 0.000, 0.000],
[0.000, cos(ax), -sin(ax), 0.0000],
[0.000, sin(ax), cos(ax), 0.000],
[0.0000, 0.0000, 0.0000, 1.0000]])
return rax
# 获取绕Y轴旋转矩阵
def GetRotateMatrix_AngleY(ay):
ray = np.array([[cos(ay), 0.0000, sin(ay), 0.0000],
[0.0000, 1.0000, 0.0000, 0.0000],
[-sin(ay), 0.0000, cos(ay), 0.0000],
[0.0000, 0.0000, 0.0000, 1.0000]])
return ray
# 获取旋转矩阵（对象坐标系坐标旋转到世界坐标系）
def GetRotationMatrix(ax, ay, az):
# 对象坐标系坐标旋转到世界坐标系，先旋转Y再X再Z,反之则先后循序相反
raz = GetRotateMatrix_AngleZ(az)
rax = GetRotateMatrix_AngleX(ax)
ray = GetRotateMatrix_AngleY(ay)
return np.dot(np.dot(ray, rax), raz)
# 世界坐标系到对象坐标系
# return np.dot(np.dot(raz, rax), ray)
# 生成空间直线起止点的圆柱体（圆或多边形顶点沿直线放样）底面和顶面多边形顶点points_start、points_end
def GetPipePoints(startpoint, endpoint, radius, num):
x1, y1, z1 = startpoint[0], startpoint[1], startpoint[2]
x2, y2, z2 = endpoint[0], endpoint[1], endpoint[2]
dx, dy, dz = x2 - x1, y2 - y1, z2 - z1
sxy = sqrt(dy**2 + dx**2)
az = -atan2(dy, dx)
ax = 0
ay = atan2(dz, sxy)
# print("sxy=", sxy)
# print("az=", az * 180 / pi)
# print("ay=", ay * 180 / pi)
rrM = GetRotationMatrix(ax, ay, az)
circlepoints = getcirclepoints(radius, num)
point3d = np.dot(circlepoints, rrM)
parr1 = np.array([x1, y1, z1, 1])
parr2 = np.array([x2, y2, z2, 1])
# print(point3d)
points_start = point3d + parr1
points_end = point3d + parr2
return points_start, points_end
def getcirclepoints(radius, num):
points = []
for i in range(num):
angle = i * pi * 2 / num
y = cos(angle) * radius
z = sin(angle) * radius
points.append([0, -y, z, 1])
return points
# 新坐标系（以空间直线(P1,P2，已知P1、P2在世界坐标系的三维坐标)为X'轴,过直线P1P2且垂直于世界OXY面的面内，垂直于P1P2的方向为Z'轴,垂直于X'和Z'轴的Y'轴构成新坐标系）
# 将新坐标系的多边形顶点转换为原世界坐标系的坐标。
def TestRotation():
startpoint = [50, 50, 50]
endpoint = [200.0, 100.0, 175.0]
radius = 20
num = 12
points_start, points_end = GetPipePoints(startpoint, endpoint, radius, num)
np.set_printoptions(suppress=True)
#定义坐标轴
fig = plt.figure()
axio = Axes3D(fig)
# 设置坐标轴标题和刻度
axio.set(xlabel='X', ylabel='Y', zlabel='Z')
axio.set_xlim3d(-50, 200)
axio.set_ylim3d(-50, 200)
axio.set_zlim3d(-50, 200)
for p in points_start:
x, y, z = p[0], p[1], p[2]
axio.scatter3D(x, y, z, cmap='Blues')
for p in points_end:
x, y, z = p[0], p[1], p[2]
print(x, y, z)
axio.scatter3D(x, y, z, cmap='Blues')
a = sqrt((startpoint[0] - endpoint[0])**2 +
(startpoint[1] - endpoint[1])**2 +
(startpoint[2] - endpoint[2])**2)
for i in range(len(points_start)):
# if i % 3 == 0: continue
pointp = points_start[i]
b = sqrt((pointp[0] - endpoint[0])**2 + (pointp[1] - endpoint[1])**2 +
(pointp[2] - endpoint[2])**2)
c = sqrt((pointp[0] - startpoint[0])**2 +
(pointp[1] - startpoint[1])**2 +
(pointp[2] - startpoint[2])**2)
cosA = (b**2 + c**2 - a**2) / (2 * b * c)
cosB = (a**2 + c**2 - b**2) / (2 * a * c)
cosC = (b**2 + a**2 - c**2) / (2 * b * a)
# print(f"cosA={cosA} cosB={cosB} cosC={cosC}")
a_angle = acos(cosA) * 180.0 / pi
b_angle = acos(cosB) * 180.0 / pi
c_angle = acos(cosC) * 180.0 / pi
print(
f"a={a} b={b} c={c} a_angle={a_angle} b_angle={b_angle} c_angle={c_angle}"
)
# print(f"a_angle={a_angle} b_angle={b_angle} c_angle={c_angle}")
linex = np.linspace(startpoint[0], endpoint[0], 100)
liney = np.linspace(startpoint[1], endpoint[1], 100)
linez = np.linspace(startpoint[2], endpoint[2], 100)
axio.plot3D(linex, liney, linez)
# XX = [point3d_2[0] for i in range(len(point3d_2))]
# YY = [point3d_2[1] for i in range(len(point3d_2))]
# ZZ = [point3d_2[2] for i in range(len(point3d_2))]
# axio.plot_surface(XX,YY,ZZ,color="r")
# print(XX)
# axio.legend()
# plt.zlabel('z')
plt.show()
def writepipeobj3d(pgconnectparam,file,facesnum):
pgconn = pghelper.connectdb(pgconnectparam)
pgcur = pgconn.cursor()
sql = '''select id,x1,y1,z1,x2,y2,z2 from js_line
where x1 is not null and
x2 is not null and
y1 is not null and
y1 is not null and
z1 is not null and
z1 is not null LIMIT 10 ;
'''
pgcur.execute(sql)
rows = pgcur.fetchall()
k = 0
num = facesnum
vertexes = []
Faces = []
for row in rows:
(id, x1, y1, z1, x2, y2, z2) = row
# 管道半径应该在数据row中获取
radius = 1.0
startpoint = [float(x1), float(y1), float(z1)]
endpoint = [float(x2), float(y2), float(z2)]
points_start, points_end = GetPipePoints(startpoint, endpoint,
radius, num)
k = len(vertexes)
for p in points_start:
x, y, z = p[0], p[1], p[2]
x = x - 3370000
y = y - 500000
vertexes.append([x, y, z])
file.write(f"v {y:.4f} {x:.4f} {z:.4f}\n")
for p in points_end:
x, y, z = p[0], p[1], p[2]
x = x - 3370000
y = y - 500000
vertexes.append([x, y, z])
file.write(f"v {y:.4f} {x:.4f} {z:.4f}\n")
file.write(f"g cube_{id}\n")
for i in range(k + 1, k + num):
# f = [i, i + 1, i + num, i + num + 1]
v1, v2, v3, v4 = i + 1, i, i + num, i + num + 1
Faces.append([v1, v2, v3, v4])
file.write(f"f {v1} {v2} {v3} {v4}\n")
v1, v2, v3, v4 = k + 1 , k + num, k + 2 * num, k + num+1
Faces.append([v1, v2, v3, v4])
file.write(f"f {v1} {v2} {v3} {v4}\n")
pgcur.close()
pgconn.close()
def main():
(dataconfigxmlfile, templatexmlfile, mdbfilename, pgconnectparam,
pipetables, surfaceobj) = config.readconfig()
file3d = "mypipe.obj"
open(file3d, "w")
with open(file3d, mode="a") as file:
writepipeobj3d(pgconnectparam,file,12)
print("管道生成完成！")
if __name__ == "__main__":
np.set_printoptions(suppress=True)
# TestRotation()
main()