```python
import numpy as np
import ipyvolume as ipv
import sympy as sy
from sympy.geometry import Plane as syPlane, Point3D as syPoint3D
import tqdm
```


```python
xyz = np.array((2 * np.random.random(1000) - 1, 2 * np.random.random(1000) - 1, 2 * np.random.random(1000) - 1))
```


```python
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
ipv.show()
```

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.


```python
x, y, z, d, e, f, a, b, c, theta = sy.symbols('x, y, z, d, e, f, a, b, c, theta')
```


```python
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))
```

To solve, we need to fill out the parameters we want. Let's try some simple things first.


```python
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)
```


```python
value = res[0].subs(x, 0.1).subs(z, 0.1)
```


```python
float(value)
```


```python
res0 = res[0]
```


```python
cone_eqn.subs
```

Ok, that seems reasonable actually!


```python
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)
```


```python
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone(res, x, z)
ipv.show()
```

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.


```python
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)
```


```python
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_2((0, 1, 0), (0, 0, 0), sy.pi/4)
ipv.show()
```


```python
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()
```

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 


```python
...
```

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)$.


```python
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)
```


```python
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_3(0.5, 0.5)
ipv.show()
```


```python
ipv.clear()
ipv.scatter(*xyz, marker='circle_2d')
plot_cone_3(1, 0.5)
ipv.show()
```


```python
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)
```


```python
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()
```


```python
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()
```


```python
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()
```

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.


```python
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)
```


```python
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()
```

Ok, that's it, now a thick cone...


```python
plot_cone = plot_cone_5
```


```python
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
```


```python
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
```


```python
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)
```


```python
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()
```

Nice!

# Filtering

This 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


```python
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})
```


```python
res
```

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.


```python
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
```


```python
def cone_opening_angle(height, radius):
    """
    Twice the opening angle is the maximum angle between directrices
    """
    return np.arctan(radius / height)
```


```python
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)))
```


```python
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
```

Ok, that should do as well, an all-in-one distance function, a bit complicated, but I see no way around that.


```python
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
```


```python
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))
```


```python
len(p_filtered)
```


```python
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()
```

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!

## Debugging

For debugging, I used some simple points like these:


```python
_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()
```

~~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).


```python
import ectopylasm.geometry
```


```python
ectopylasm.geometry.filter_points_cone(_xyz, 0.5, 0.5, 0.2)
```


```python
filter_points_cone(_xyz, 0.5, 0.5, 0.2)
```

Arrgh


```python
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()
```


```python
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)
```


```python
d_cone1, d_cone2
```


```python

```

748bc37f2c8968dc7b209d38382893d35b85e31c
ipynb
Jupyter Notebook
req2.2_filter_cone.ipynb
sundial-pointcloud-geometry/explore
f2aca1b787a5bcc8a187fefdf5048451fca1b2e1
MIT
req2.2_filter_cone.ipynb
sundial-pointcloud-geometry/explore
f2aca1b787a5bcc8a187fefdf5048451fca1b2e1
MIT
req2.2_filter_cone.ipynb
sundial-pointcloud-geometry/explore
f2aca1b787a5bcc8a187fefdf5048451fca1b2e1
MIT
Qwen/Qwen-72B
1. YES
2. YES