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

问题描述

我想在尺寸为100x100x100的3d矩阵中包含一个圆锥体，其中1s代表圆锥体。然后，通过将它们加在一起并查看显示2的位置，将它们用于计算不同圆锥的相交。例如在2d中：

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,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,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,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,projection='3d')

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

cone1 = np.zeros((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,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

解决方法

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

恐怕分析方法太复杂了。但是出于绘图目的，您可以尝试简化一下。制作方程式-或使用顶点，轴方向，半径（如针对t参数的方程式）或通常使用quadric form进行参数化。

然后扫描x逐层扫描值。将电流x的值代入初始方程，得到两个未知数y和z的第二次幂的两个方程的系统。它可以解析地解决（使用scipy？），因此您最多可以获取4个解（或圆方程！）。绘制相应的点，然后转到下一个x。

geometry python