如何有效地计算矩阵内的圆锥

Question

我想在大小为 100x100x100 的 3d 矩阵内有一个圆锥体，其中 1 代表圆锥体。然后通过将它们加在一起并查看数字 2 的显示位置来计算不同锥体的交集。 二维示例：

100000001

010000010

001000100

000101000

000010000

我目前有一个代码，它可以工作，但效率很低。这是因为，为了获得侧面的特定厚度，我多次计算圆锥体。在最后的程序中，这些锥体可能会变得很大，因为我需要计算它们的负载，所以我不得不等待很多。我想也许有人可以给我一个小费或向一个新的方向点头，如何更有效地解决这个问题。

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook
from scipy.linalg import norm
import matplotlib.pyplot as plt

def truncated_cone(p0, p1, R0, R1, color, cone_length):
    #
    #Based on https://stackoverflow.com/a/39823124/190597 (astrokeat)
    #
    # vector in direction of axis
    v = p1-p0
    # find magnitude of vector
    mag = norm(v)
    # unit vector in direction of axis
    v = v / mag
    # make some vector not in the same direction as v
    not_v = np.array([1, 1, 0])
    not_v = not_v/norm(not_v)
    if (v == not_v).all():
        not_v = np.array([0, 1, 0])
    # make vector perpendicular to v

    n1 = np.cross(v, not_v)
    # print n1,'\t',norm(n1)
    # normalize n1
    n1 /= norm(n1)
    # make unit vector perpendicular to v and n1
    n2 = np.cross(v, n1)
    # surface ranges over t from 0 to length of axis and 0 to 2*pi

    #depending on the cone n has to be high enough for an accurate cone
    circumference = int(2*np.pi*R1+1)
    pythagoras = int(np.sqrt(R1**2+cone_length**2)+1)
    if  circumference<pythagoras:
        n = pythagoras
    else:
        n = circumference

    t = np.linspace(0, mag, n)
    theta = np.linspace(0, 2 * np.pi, n)
    # use meshgrid to make 2d arrays
    t, theta = np.meshgrid(t, theta)

    R = np.linspace(R0, R1, n)
    # generate coordinates for surface
    X, Y, Z = [p0[i] + v[i] * t + R *
               np.sin(theta) * n1[i] + R * np.cos(theta) * n2[i] for i in [0, 1, 2]]

    # to fit in the coordinate system / matrices
    X = np.around(X)
    Y = np.around(Y)
    Z = np.around(Z)

    #Matrix which will contain the cone
    cone_array = np.zeros((100, 100, 100), dtype='int8')

    for i in np.arange(n):
        for j in np.arange(n):
            x = int(X[i][j])
            y = int(Y[i][j])
            z = int(Z[i][j])

            #stay inside the boundaries of the intersection matrix
            if x > 99 or y > 99 or z > 99 or x < 0 or y < 0 or z < 0: 
                pass
            else:
                cone_array[x][y][z] = 1

    return cone_array

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')

p0 = np.array([20,20,20])
p1 = np.array([50,50,50])
R0 = 0
R1 = 20

cone1 = np.zeros((100, 100, 100), dtype='int8')

#in order to get the thickness of the lateral surface the cone is calculated multiple times
for i in np.arange(10, 20+1, 1): 
    cone1 += truncated_cone(p0, p1, 0, i, 'blue', 52)

    
cone1_step = np.where(cone1>0)
ax.scatter(cone1_step[0], cone1_step[1], cone1_step[2], c='blue', zorder=1)

ax.set_xlim(0, 100)
ax.set_ylim(0, 100)
ax.set_zlim(0, 100)
        

结果如下所示： Cone from the code

Answer 1

这个问题并不简单。两个任意圆锥的交点可能是空的、点、圆（当轴共线时）、弯曲的表面、两个表面等。

恐怕解析解太复杂了。但出于绘图目的，您可以尝试将其简化一点。制作方程式 - 或使用顶点、轴方向、半径以参数形式（例如针对t 参数的方程式）或一般quadric form。

然后逐层扫描x值。将电流x 值代入初始方程，得到两个未知数y 和z 的二次方方程组。它可以解析求解（使用 scipy？），因此您最多可以获得 4 个解（或圆方程！）。画出对应点，跳转到下一个x。