How to efficiently calculate a cone inside a matrix

Question

I want to have a cone inside a 3d Matrix of the size 100x100x100 where 1s represent the cone. These are then used to calculate the intersection of different cones by adding them together and seeing where the number 2 is displayed. example In 2d:

100000001

010000010

001000100

000101000

000010000

I currently have a code, which works but is very inefficient.This comes from the fact, that in order to get a specific thickness of the lateral surface I compute the cone multiple times. In the final program these cones can get quite big and as I need to calculate loads of them I have to wait a lot. I thought maybe someone can give me a tip or nod in a new direction how this could be solved more efficiently.

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)

The result looks something like this: Cone from the code

score 0 · Answer 1 · answered Oct 07 '20 at 10:22

This problem is not simple. Intersection of two arbitrary cones migth be empty, point, circle (when axes are collinear), crooked surface, two surfaces etc.

Analytical solution is too complex, I'm afraid. But for plotting purposes you can try to simplify it a bit. Make equations - or in parametric form using vertex point, axis direction, radius, (like your equations against t parameter), or in general quadric form.

Then scan x values layer by layer. Substitute current x value into initial equations and get a system of two equations of the second power for two unknowns y and z. It is solvable analytically (with scipy?), so you can get upto 4 solutions (or circle equation!). Plot corresponding points and go to the next x.

How to efficiently calculate a cone inside a matrix

1 Answers1