<jupyter_start><jupyter_text>Now, similarly to the [req2.1 notebook](req2.1_filter_plane.ipynb), we will want to filter and plot here a thick cone, i.e. basically two  the same cones where their surfaces are separatated by some thickness.The equation for a generically oriented and located regular (non-slanted) cone is given on [Wikipedia](https://en.wikipedia.org/wiki/Cone) as $F(u)=0$ where:$$F(u) = u \cdot d - |d| |u| \cos \theta$$where $u=(x,y,z)$, $d=(d,e,f)$ is the vector to which the cone axis is parallel and the aperture is $2\theta$. The vertex of this cone is at the origin, so to displace it to $(a,b,c)$ you need to subtract an additional vector $o = (a,b,c)$ from u:$$F(u) = (u - o) \cdot d - |d| |u - o| \cos \theta$$Alternatively, one can apparently use the equation$$F(u) = ((u - o) \cdot d)^2 - (d \cdot d) ((u - o) \cdot (u - o)) (\cos \theta)^2$$This is actually easier, because you can skip the square roots necessary for the norms.<jupyter_code>x, y, z, d, e, f, a, b, c, theta = sy.symbols('x, y, z, d, e, f, a, b, c, theta')
cone_eqn = ((x-a)*d + (y-b)*e + (z-c)*f)**2 \
           - sy.cos(theta)**2 * (d * d + e * e + f * f) \
           * ((x-a) * (x-a) + (y-b) * (y-b) + (z-c) * (z-c))<jupyter_output><empty_output><jupyter_text>To solve, we need to fill out the parameters we want. Let's try some simple things first.<jupyter_code>res = sy.solve(cone_eqn.subs(theta, sy.pi/4).subs(a, 0).subs(b, 0).subs(c, 0).subs(d, 0).subs(e, 1).subs(f, 0), y)
value = res[0].subs(x, 0.1).subs(z, 0.1)
float(value)
res0 = res[0]
cone_eqn.subs<jupyter_output><empty_output><jupyter_text>Ok, that seems reasonable actually!<jupyter_code>def plot_cone(res, x_sy, z_sy):
    """
    Draw a cone.
    """
    # get box limits in two dimensions
    x_lim = ipv.pylab.gcf().xlim
    z_lim = ipv.pylab.gcf().zlim
    x_steps = np.linspace(*x_lim, 10)
    z_steps = np.linspace(*z_lim, 10)
    x, z = np.meshgrid(x_steps, z_steps)

    # find corresponding y coordinates
    y1, y2 = [], []
    for xi, zi in zip(x.flatten(), z.flatten()):
        y1.append(float(res[0].evalf(subs={x_sy: xi, z_sy: zi})))
        y2.append(float(res[1].evalf(subs={x_sy: xi, z_sy: zi})))

    # plot
    ipv.plot_surface(x, np.array(y1), z)
    ipv.plot_surface(x, np.array(y2), z)
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone(res, x, z)
ipv.show()<jupyter_output><empty_output><jupyter_text>Nice! Could be better though.- The origin point must be included.- It looks a bit odd that the edges are cut off. It makes sense if you want to fill the whole box, but probably most people will prefer a circular ending. This we could add as an option.- Just pick one of the two solutions, probably the one that is actually in the same direction as the $d$ vector.Let's also try to make it a bit more general and try out some different examples to test other corner cases.<jupyter_code>def plot_cone_2(axis, origin, half_aperture):
    """
    Draw a cone.
    
    axis: three-component vector along which the cone axis lies
    origin: location of the apex of the cone
    half_aperture: the aperture of the cone is 2*half_aperture
    """
    x, y, z, d, e, f, a, b, c, theta = sy.symbols('x, y, z, d, e, f, a, b, c, theta')
    cone_eqn = ((x-a)*d + (y-b)*e + (z-c)*f)**2 \
               - sy.cos(theta)**2 * (d * d + e * e + f * f) \
               * ((x-a) * (x-a) + (y-b) * (y-b) + (z-c) * (z-c))
    
    res = sy.solve(cone_eqn.subs({d: axis[0], e: axis[1], f: axis[2],
                                  a: origin[0], b: origin[1], c: origin[2],
                                  theta: half_aperture}), y)
    
#     print(res)
    
    # get box limits in two dimensions
    x_lim = ipv.pylab.gcf().xlim
    z_lim = ipv.pylab.gcf().zlim
    x_steps = np.linspace(*x_lim, 20)
    z_steps = np.linspace(*z_lim, 20)
    x_array, z_array = np.meshgrid(x_steps, z_steps)
    
    y_lim = ipv.pylab.gcf().ylim
    
    # find corresponding y coordinates
    y1, y2 = [], []
    for xi, zi in zip(x_array.flatten(), z_array.flatten()):
        y1.append(float(res[0].evalf(subs={x: xi, z: zi})))
        y2.append(float(res[1].evalf(subs={x: xi, z: zi})))

    # plot
    ipv.plot_surface(x_array, np.array(y1), z_array)
    ipv.plot_surface(x_array, np.array(y2), z_array)
    
    ipv.ylim(*y_lim)
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_2((0, 1, 0), (0, 0, 0), sy.pi/4)
ipv.show()
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_2((0.25, 0.5, 0), (0, 0.5, 0), sy.pi/8)
ipv.show()<jupyter_output><empty_output><jupyter_text>Arrgh, what's going on here?Oh, wait, this must be what happens when the x and z parameters are, in fact, outside of where the cone is going to be.Is there a more reliable way to find points x, y, z that are actually on the cone surface instead of just randomly (or rather grid-wise) probing x and z and hoping there's a y point there? The trick is to find a parametric form, I think. Wikipedia has this form for a cone along the z-axis, so we could use that, but then we'd still have to displace it and rotate it.Maybe we can find our own parametric form though.One dimension that we would probably want is the same as the $u$ one of the Wikipedia parameterization, i.e. a coordinate that runs along the length of the axis up to its height $h$, i.e. $u \in [0,h]$. I think this means it should be that$$u^2 = x^2 + y^2 + z^2$$so that $u$ is the square-root of the sum of the three components. This indeed makes sense, because it has two solutions, a positive and a negative one, which technically the cone equations indeed must have. However, we only use the positive range.The remaining two parameters $s$ and $t$ must describe the remaining two<jupyter_code>...<jupyter_output><empty_output><jupyter_text>Ok... this is going to be too much work, and won't generalize nicely to other surfaces anyway. Let's try the easier approach of just rotating the simple parametric form. This will also automatically fix all three "niceness" issues mentioned above.Wikipedia's parametric description for a cone along the z-axis:$$F(s,t,u) = \left(u \tan s \cos t, u \tan s \sin t, u \right)$$where $s,t,u$ range over $[0,\theta)$, $[0,2\pi)$, and $[0,h]$, respectively. Actually, Wolfram Mathworld has a simpler parameterization with only two parameters:$$x = \frac{h-u}{h} r \cos\theta	\\y = \frac{h-u}{h} r \sin\theta	\\z = u$$for $u$ in $[0,h]$ and $\theta$ in $[0,2\pi)$.<jupyter_code>def plot_cone_3(height, radius, N_steps=20):
    """
    Draw a cone.
    
    height: height along the z-axis
    radius: radius of the circle
    """
    h, r, u, theta = sy.symbols('h, r, u, theta')
    x_eqn = (h - u) / h * r * sy.cos(theta)
    y_eqn = (h - u) / h * r * sy.sin(theta)
    z_eqn = u
            
    # get box limits in two dimensions
#     x_lim = ipv.pylab.gcf().xlim
#     z_lim = ipv.pylab.gcf().zlim
#     x_steps = np.linspace(*x_lim, 20)
#     z_steps = np.linspace(*z_lim, 20)
#     x_array, z_array = np.meshgrid(x_steps, z_steps)
    u_steps = np.linspace(0, height, N_steps)
    theta_steps = np.linspace(0, 2 * np.pi, N_steps)
    u_array, theta_array = np.meshgrid(u_steps, theta_steps)
    
    y_lim = ipv.pylab.gcf().ylim
    
    # find corresponding y coordinates
    x, y, z = [], [], []
    for ui, thetai in zip(u_array.flatten(), theta_array.flatten()):
        x.append(float(x_eqn.evalf(subs={h: height, r: radius, u: ui, theta: thetai})))
        y.append(float(y_eqn.evalf(subs={h: height, r: radius, u: ui, theta: thetai})))
        z.append(float(z_eqn.evalf(subs={h: height, r: radius, u: ui, theta: thetai})))

    # plot
    ipv.plot_surface(np.array(x).reshape(u_array.shape),
                     np.array(y).reshape(u_array.shape),
                     np.array(z).reshape(u_array.shape))
    
#     ipv.ylim(*y_lim)
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_3(0.5, 0.5)
ipv.show()
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_3(1, 0.5)
ipv.show()
def plot_cone_4(height, radius, rot_x=2*np.pi, rot_y=2*np.pi, N_steps=20):
    """
    Draw a cone.
    
    height: height along the z-axis
    radius: radius of the circle
    rot_y: rotation angle about the y axis (radians)
    rot_z: rotation angle about the z axis (radians)
    N_steps: number of steps in the parametric range used for drawing
    """
    h, r, u, theta = sy.symbols('h, r, u, theta')
    x_eqn = (h - u) / h * r * sy.cos(theta)
    y_eqn = (h - u) / h * r * sy.sin(theta)
    z_eqn = u

    x_rot_x = x_eqn
    y_rot_x = y_eqn * sy.cos(rot_x) + z_eqn * sy.sin(rot_x)
    z_rot_x = - y_eqn * sy.sin(rot_x) + z_eqn * sy.cos(rot_x)

    x_rot_y = x_rot_x * sy.cos(rot_y) + z_rot_x * sy.sin(rot_y)
    y_rot_y = y_rot_x
    z_rot_y = - x_rot_x * sy.sin(rot_y) + z_rot_x * sy.cos(rot_y)
            
    # get box limits in two dimensions
#     x_lim = ipv.pylab.gcf().xlim
#     z_lim = ipv.pylab.gcf().zlim
#     x_steps = np.linspace(*x_lim, 20)
#     z_steps = np.linspace(*z_lim, 20)
#     x_array, z_array = np.meshgrid(x_steps, z_steps)
    u_steps = np.linspace(0, height, N_steps)
    theta_steps = np.linspace(0, 2 * np.pi, N_steps)
    u_array, theta_array = np.meshgrid(u_steps, theta_steps)
    
    y_lim = ipv.pylab.gcf().ylim
    
    # find corresponding y coordinates
    x, y, z = [], [], []
    for ui, thetai in zip(u_array.flatten(), theta_array.flatten()):
        x.append(float(x_rot_y.evalf(subs={h: height, r: radius, u: ui, theta: thetai})))
        y.append(float(y_rot_y.evalf(subs={h: height, r: radius, u: ui, theta: thetai})))
        z.append(float(z_rot_y.evalf(subs={h: height, r: radius, u: ui, theta: thetai})))

    # plot
    ipv.plot_surface(np.array(x).reshape(u_array.shape),
                     np.array(y).reshape(u_array.shape),
                     np.array(z).reshape(u_array.shape))
    
#     ipv.ylim(*y_lim)
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_4(1, 0.5, rot_y=np.pi/4)#, rot_z=np.pi/2)
ipv.show()
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_4(1, 0.5, rot_x=np.pi/9)#, rot_z=np.pi/2)
ipv.show()
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_4(1, 0.5, rot_x=np.pi/9, rot_y=np.pi/4)#, rot_z=np.pi/2)
ipv.show()<jupyter_output><empty_output><jupyter_text>Ok, this we can work with, it rotates as expected. Probably, we should also add a version with a direction vector, but this will do for now. And we need a version in which you specify two positions: the base center and the apex position.Now just add a translation vector.<jupyter_code>def plot_cone_5(height, radius, rot_x=2*np.pi, rot_y=2*np.pi, base_pos=(0, 0, 0), N_steps=20, **kwargs):
    """
    Draw a cone.
    
    height: height along the z-axis
    radius: radius of the circle
    rot_y: rotation angle about the y axis (radians)
    rot_z: rotation angle about the z axis (radians)
    base_pos: translation of base of cone to this position, iterable of three numbers
    N_steps: number of steps in the parametric range used for drawing
    """
    h, r, u, theta = sy.symbols('h, r, u, theta')
    x_eqn = (h - u) / h * r * sy.cos(theta)
    y_eqn = (h - u) / h * r * sy.sin(theta)
    z_eqn = u

    x_rot_x = x_eqn
    y_rot_x = y_eqn * sy.cos(rot_x) - z_eqn * sy.sin(rot_x)
    z_rot_x = y_eqn * sy.sin(rot_x) + z_eqn * sy.cos(rot_x)

    x_rot_y = x_rot_x * sy.cos(rot_y) - z_rot_x * sy.sin(rot_y) + base_pos[0]
    y_rot_y = y_rot_x + base_pos[1]
    z_rot_y = x_rot_x * sy.sin(rot_y) + z_rot_x * sy.cos(rot_y) + base_pos[2]
            
    # get box limits in two dimensions
    u_steps = np.linspace(0, height, N_steps)
    theta_steps = np.linspace(0, 2 * np.pi, N_steps)
    u_array, theta_array = np.meshgrid(u_steps, theta_steps)
    
    y_lim = ipv.pylab.gcf().ylim
    
    # find corresponding y coordinates
    x, y, z = [], [], []
    for ui, thetai in zip(u_array.flatten(), theta_array.flatten()):
        x.append(float(x_rot_y.evalf(subs={h: height, r: radius, u: ui, theta: thetai})))
        y.append(float(y_rot_y.evalf(subs={h: height, r: radius, u: ui, theta: thetai})))
        z.append(float(z_rot_y.evalf(subs={h: height, r: radius, u: ui, theta: thetai})))

    # plot
    ipv.plot_surface(np.array(x).reshape(u_array.shape),
                     np.array(y).reshape(u_array.shape),
                     np.array(z).reshape(u_array.shape), **kwargs)
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_5(0.5, 0.5, rot_x=np.pi/9, rot_y=np.pi/4, base_pos=(0, 0.5, 0))#, rot_z=np.pi/2)
ipv.show()<jupyter_output><empty_output><jupyter_text>Ok, that's it, now a thick cone...<jupyter_code>plot_cone = plot_cone_5
def cone_axis_from_rotation(rot_x, rot_y):
    # z-unit vector (0, 0, 1) rotated twice
    cone_axis = (0, -np.sin(rot_x), np.cos(rot_x))  # rotation around x-axis
    cone_axis = np.array((-np.sin(rot_y) * cone_axis[2],
                          cone_axis[1],
                          np.cos(rot_y) * cone_axis[2]))  # around y
    return cone_axis
def thick_cone_base_positions(height, radius, thickness, rot_x, rot_y, base_pos):
    thickness = abs(thickness)
    base_distance = thickness / radius * height * np.sqrt(1 + radius**2 / height**2)  # trigonometry
    
    cone_axis = cone_axis_from_rotation(rot_x, rot_y)
    
    base_pos_1 = np.array(base_pos) - cone_axis * 0.5 * base_distance
    base_pos_2 = np.array(base_pos) + cone_axis * 0.5 * base_distance

    return base_pos_1, base_pos_2
def plot_thick_cone(height, radius, thickness,
                    rot_x=2*np.pi, rot_y=2*np.pi, base_pos=(0, 0, 0),
                    **kwargs):
    """
    Plot two cones separated by a distance `thickness`.
    
    Parameters: same as plot_cone, plus `thickness`.
    """
    base_pos_1, base_pos_2 = thick_cone_base_positions(height, radius, thickness, rot_x, rot_y, base_pos)
    plot_cone(height, radius, rot_x=rot_x, rot_y=rot_y, base_pos=base_pos_1, **kwargs)
    kwargs.pop('color', None)
    plot_cone(height, radius, rot_x=rot_x, rot_y=rot_y, base_pos=base_pos_2, color='blue', **kwargs)
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_thick_cone(0.5, 0.5, 0.2, rot_x=np.pi/9, rot_y=np.pi/4, base_pos=(0, 0.5, 0))
ipv.show()<jupyter_output><empty_output><jupyter_text>Nice!  FilteringThis is going to be a bit more difficult. We need a closest distance of each point to both surfaces. How to do this?Maybe we can try using Sympy to solve the system of equations describing the minimum distance of a point $P=(p_1, p_2, p_3)$ to a cone surface, like e.g. here: https://math.stackexchange.com/a/880971/258876<jupyter_code>h, r, u, theta, p1, p2, p3 = sy.symbols('h, r, u, theta, p1, p2, p3')
x_eqn = (h - u) / h * r * sy.cos(theta)
y_eqn = (h - u) / h * r * sy.sin(theta)
z_eqn = u

distance_eqn = (x_eqn - p1)**2 + (y_eqn - p2)**2 + (z_eqn - p3)**2

system = [sy.Derivative(distance_eqn, u),
          sy.Derivative(distance_eqn, theta)]

res = sy.solve(system, distance_eqn, subs={h: 0.5, r: 0.5, p1: 0, p2: 0, p3: 0})
res<jupyter_output><empty_output><jupyter_text>Not sure how to fix this.Another possibility may be to first filter by looking whether the angle of the vector from apex to $P$ with the vector along the cone axis is larger than $90^\circ + \theta/2$ where $\theta$ is the opening angle.<jupyter_code>def cone_apex_position(height, rot_x=2*np.pi, rot_y=2*np.pi, base_pos=(0, 0, 0)):
    cone_axis = cone_axis_from_rotation(rot_x, rot_y)
    return np.array(base_pos) + cone_axis * height
def cone_opening_angle(height, radius):
    """
    Twice the opening angle is the maximum angle between directrices
    """
    return np.arctan(radius / height)
def angle_between_two_vectors(a, b):
    return np.arccos(np.sum(a * b) / np.sqrt(np.sum(a**2)) / np.sqrt(np.sum(b**2)))
def point_distance_to_cone(P, height, radius,
                           rot_x=2*np.pi, rot_y=2*np.pi, base_pos=(0, 0, 0),
                           return_extra=False):
    """
    Check whether for a point P, the shortest path to the cone is
    perpendicular to the cone surface (and if so, return it). If
    not, it is either "above" the apex and the shortest path is simply
    the line straight to the apex, or it is "below" the base, and the
    shortest path is the shortest path to the directrix (the base
    circle).

    This function returns a second value depending on which of the
    three above cases is true for point P. If we're using the
    perpendicular, it is True, if we're above the apex it is False and
    if it is below the base, it is None.

    Extra values can be returned to be reused outside the function by
    setting return_extra to True.
    """
    cone_axis = cone_axis_from_rotation(rot_x, rot_y)
    apex_pos = cone_apex_position(height, rot_x=rot_x, rot_y=rot_y, base_pos=base_pos)
    point_apex_vec = np.array(P) - apex_pos
    point_apex_angle = np.pi - angle_between_two_vectors(cone_axis, point_apex_vec)
    opening_angle = cone_opening_angle(height, radius)

    # for the second conditional, we need the length of the component of the
    # difference vector between P and apex along the closest generatrix
    point_apex_generatrix_angle = point_apex_angle - opening_angle
    point_apex_distance = np.sqrt(np.sum(point_apex_vec**2))
    point_apex_generatrix_component = point_apex_distance * np.cos(point_apex_generatrix_angle)
    generatrix_length = np.sqrt(radius**2 + height**2)

    returnees = {}
    if return_extra:
        returnees['opening_angle'] = opening_angle
        returnees['point_apex_angle'] = point_apex_angle

    if point_apex_angle > opening_angle + np.pi / 2:
        # "above" the apex
        return point_apex_distance, False, returnees
    elif point_apex_generatrix_component > generatrix_length:
        # "below" the directrix
        # use cosine rule to find length of third side
        return np.sqrt(point_apex_distance**2 + generatrix_length**2
                       - 2 * point_apex_distance * generatrix_length
                       * np.cos(point_apex_generatrix_angle)), None, returnees
    else:
        # "perpendicular" to a generatrix
        return point_apex_distance * np.sin(point_apex_generatrix_angle), True, returnees<jupyter_output><empty_output><jupyter_text>Ok, that should do as well, an all-in-one distance function, a bit complicated, but I see no way around that.<jupyter_code>def filter_points_cone(points_xyz, height, radius, thickness,
                       rot_x=2*np.pi, rot_y=2*np.pi, base_pos=(0, 0, 0), verbose=False):
    base_pos_1, base_pos_2 = thick_cone_base_positions(height, radius, thickness, rot_x, rot_y, base_pos)
    
    p_filtered = []
    for ix, p_i in tqdm.tqdm(enumerate(points_xyz.T)):
        d_cone1, flag_cone1, vals1 = point_distance_to_cone(p_i, height, radius,
                                                            rot_x=rot_x, rot_y=rot_y,
                                                            base_pos=base_pos_1, return_extra=True)
        d_cone2, flag_cone2, _ = point_distance_to_cone(p_i, height, radius,
                                                        rot_x=rot_x, rot_y=rot_y,
                                                        base_pos=base_pos_2, return_extra=True)
        if flag_cone2 is False or flag_cone1 is None:
            # it is definitely outside of the cones' range
            pass
            if verbose: print(f"case 1: {p_i} was ignored")
        elif flag_cone1 is False:
            # the first condition is logically enclosed in the second, but the
            # first is faster and already covers a large part of the cases/volume:
            if d_cone1 <= thickness or \
               d_cone1 <= thickness / np.cos(vals1['point_apex_angle'] - vals1['opening_angle'] - np.pi/2):
                p_filtered.append(p_i)
                if verbose: print(f"case 2: {p_i} was added")
            else:
                pass
                if verbose: print(f"case 3: {p_i} was ignored")
        elif d_cone1 <= thickness and d_cone2 <= thickness:
            p_filtered.append(p_i)
            if verbose: print(f"case 4: {p_i} was added")
    return p_filtered
p_filtered = filter_points_cone(xyz, 0.5, 0.5, 0.2, rot_x=np.pi/9, rot_y=np.pi/4, base_pos=(0, 0.5, 0))
len(p_filtered)
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_thick_cone(0.5, 0.5, 0.2, rot_x=np.pi/9, rot_y=np.pi/4, base_pos=(0, 0.5, 0))
ipv.scatter(*np.array(p_filtered).T, marker='circle_2d', color='blue')
ipv.show()<jupyter_output><empty_output><jupyter_text>First try: That makes no sense...Second try (modified code, added `np.pi - ` in `point_distance_to_cone`): ok, better, at least the points are now in one of the cones, but there are still some weird things; false negatives (between the cones) and false positives (in the red cone).Third try (modified code, was using opening angle as if it was twice what it is, so now without divisions by 2): awesome!  DebuggingFor debugging, I used some simple points like these:<jupyter_code>_xyz = np.array([[0, 0, 0.5], # inside blue, above red
                 [0, 0, 0.7], # above blue (i.e. above both), within "thickness" from blue
                 [0, 0, 1],   # far above both
                 [0, 0, 0.3], # inside red
                 [0, 0, 0.1], # inside red, "below" blue
                 [0, 0, -0.1], # below both
                ]).T

ipv.clear()
ipv.scatter(*_xyz, color="black")
p_filtered = filter_points_cone(_xyz, 0.5, 0.5, 0.2, verbose=True)
plot_thick_cone(0.5, 0.5, 0.2)
ipv.scatter(*np.array(p_filtered).T, color='green')
ipv.show()<jupyter_output><empty_output><jupyter_text>~~Ok, so the ones above are correct and the one inside is correct, but the ones inside the red cone are wrong. Let's go through the wrong cases one by one.~~~~Why is `[0, 0, 0.3]` added? It's in case 4, so `d_cone1 <= thickness and d_cone2 <= thickness` is True. What are the values?~~Before, the ones in the red cone were also added, which is what I'm debugging below (but makes no sense anymore now, i.e. now it makes sense, because it's no longer wrong).<jupyter_code>import ectopylasm.geometry
ectopylasm.geometry.filter_points_cone(_xyz, 0.5, 0.5, 0.2)
filter_points_cone(_xyz, 0.5, 0.5, 0.2)<jupyter_output><empty_output><jupyter_text>Arrgh<jupyter_code>p_filtered = ectopylasm.geometry.filter_points_cone(xyz, 0.5, 0.5, 0.2, rot_x=np.pi/9, rot_y=np.pi/4, base_pos=(0, 0.5, 0))
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_thick_cone(0.5, 0.5, 0.2, rot_x=np.pi/9, rot_y=np.pi/4, base_pos=(0, 0.5, 0))
ipv.scatter(*np.array(p_filtered).T, marker='circle_2d', color='blue')
ipv.show()
height = 0.5
radius = 0.5
thickness = 0.2
rot_x = 0
rot_y = 0
base_pos = (0, 0, 0)

p_i = [0, 0, 0.3]

base_pos_1, base_pos_2 = thick_cone_base_positions(height, radius, thickness, rot_x, rot_y, base_pos)

d_cone1, flag_cone1, vals1 = point_distance_to_cone(p_i, height, radius,
                                                    rot_x=rot_x, rot_y=rot_y,
                                                    base_pos=base_pos_1, return_extra=True)
d_cone2, flag_cone2, _ = point_distance_to_cone(p_i, height, radius,
                                                rot_x=rot_x, rot_y=rot_y,
                                                base_pos=base_pos_2, return_extra=True)
d_cone1, d_cone2<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13276654214436, 'fraction_numerical': 0.018721284726695007, 'mean_word_length': 3.0752405949256345, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 59, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '296741', 'license': 'MIT', 'n_tokens_mistral': 8977, 'n_tokens_neox': 8613, 'n_words': 2404, 'path': 'req2.2_filter_cone.ipynb', 'repo_name': 'sundial-pointcloud-geometry/explore', 'size': '23289'}