mikedh / trimesh

Python library for loading and using triangular meshes.
https://trimesh.org
MIT License
2.99k stars 577 forks source link

Section_multiplane issues, fully captured by an example #743

Open pmberkeley opened 4 years ago

pmberkeley commented 4 years ago

I have tried to define the issue a little better based on the example you included here: https://github.com/mikedh/trimesh/issues/688#issuecomment-572833398_

Here is a slight modification to your example:

import trimesh
import numpy as np

if __name__ == '__main__':

    m = trimesh.creation.box()

    sx = m.section_multiplane(plane_origin=[-0.5,0,0],
                             plane_normal=[1,0,0],
                             heights=np.linspace(0.0, 1.0, 10))

    sy = m.section_multiplane(plane_origin=[0,-0.5,0],
                         plane_normal=[0,1,0],
                         heights=np.linspace(0.0, 1.0, 10))

    sz = m.section_multiplane(plane_origin=[0,0,-0.5],
                         plane_normal=[0,0,1],
                         heights=np.linspace(0.0, 1.0, 10))

    sr = [m.section(plane_origin=[0,0,-0.5 +h],
                    plane_normal=[0,0,1])
      for h in np.linspace(0.0, 1.0, 10.0)]

If you check for the transforms to 3D on sz, sy, and sx, you see that they do not all correspond to what they should.

sx[0].metadata['to_3D']
sy[0].metadata['to_3D']
sz[0].metadata['to_3D']

the sz transform is as expected, but the sx and sy transforms have -1s in there that probably explain the axis inversion behavior I was witnessing. Example for sx:

array([[ 0. ,  0. ,  1. , -0.5],
       [ 0. ,  1. ,  0. ,  0. ],
       [-1. ,  0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  1. ]])

Then, if you use m.bounds[0] for the plane origin, section_multiplane will produce a None result for sections up until the positive quadrant, in addition to the inverted axis behavior. The section method does not have this behavior.

import trimesh
import numpy as np

if __name__ == '__main__':
    #trimesh.util.attach_to_log()

    m = trimesh.creation.box()

    #T = [[1,0,0,-m.bounds[0,0]],[0,1,0,-m.bounds[0,1]],[0,0,1,-m.bounds[0,2]],[0,0,0,1]]
    #m.vertices = trimesh.transform_points(m.vertices,T,translate=True)

    def mesh_sections(mesh,axis,step):
        plane_normal = [0,0,0]
        plane_normal[axis] = 1
        plane_origin = mesh.bounds[0]
        extents = [mesh.bounds[0,axis],mesh.bounds[1,axis]]
        levels = np.arange(*extents, step=step)
        sections = mesh.section_multiplane(plane_origin=plane_origin,plane_normal=plane_normal,heights=levels)
        return sections, levels

    mx_sections, mx_levels = mesh_sections(m,0,.1)
    my_sections, my_levels = mesh_sections(m,1,.1)
    mz_sections, mz_levels = mesh_sections(m,2,.1)

    srx = [m.section(plane_origin=m.bounds[0],
                    plane_normal=[1,0,0])
          for h in np.arange(m.bounds[0,0],m.bounds[1,0],.1)]

    sry = [m.section(plane_origin=m.bounds[0],
                    plane_normal=[0,1,0])
          for h in np.arange(m.bounds[0,1],m.bounds[1,1],.1)]

    srz = [m.section(plane_origin=m.bounds[0],
                    plane_normal=[0,0,1])
          for h in np.arange(m.bounds[0,2],m.bounds[1,2],.1)]

the results for mx_sections:

[None,
 None,
 None,
 None,
 None,
 <trimesh.path.path.Path2D at 0x1d2534926c8>,
 <trimesh.path.path.Path2D at 0x1d2531aea48>,
 <trimesh.path.path.Path2D at 0x1d2531ae288>,
 <trimesh.path.path.Path2D at 0x1d2534c79c8>,
 <trimesh.path.path.Path2D at 0x1d2534c7888>]

vs the results for srx:

[<trimesh.path.path.Path3D at 0x1d2534c70c8>,
 <trimesh.path.path.Path3D at 0x1d2531aca88>,
 <trimesh.path.path.Path3D at 0x1d2531ac1c8>,
 <trimesh.path.path.Path3D at 0x1d2531ace88>,
 <trimesh.path.path.Path3D at 0x1d2531ac648>,
 <trimesh.path.path.Path3D at 0x1d2531ac888>,
 <trimesh.path.path.Path3D at 0x1d2531a9388>,
 <trimesh.path.path.Path3D at 0x1d2531a9288>,
 <trimesh.path.path.Path3D at 0x1d2531a92c8>,
 <trimesh.path.path.Path3D at 0x1d2534c7088>]

I hope this helps clarify what the problems are with section_multiplane; I'm pretty sure they are there but I don't know what specifically is going on.

HakuG commented 4 years ago

Also getting some weird results due to section_multiplane compared to section.

using mesh.section seems to maintain coordinate space, while section_multiplane doesn't. Not sure if that's deliberate, but for my application, I need to know the original vertex directions.

pmberkeley commented 4 years ago

I have since found more problems with to_planar, and after doing more digging and troubleshooting, I'm noticing two things. When you put in the identity transformation matrix in this function

[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]

you get the error:

ValueError: points are not flat!

You have to put in a messed up transformation matrix to avoid this.

[[0,0,1,0],[0,1,0,0],[1,0,0,0],[0,0,0,1]]

And when you convert it to_2D or to_3D, it adds -0 into the matrix, which is supposed to be a problem with np.linalg.inv(), which is solved by np.linalg.pinv().

So I think something may be going on with both to_planar and mesh_multiplane, both places where np.linalg.inv() is used.

Something is also going on with the transformation matrices, but I haven't gotten to the bottom of that yet.

mikedh commented 4 years ago

Oh interesting... didn't know -0 was a possible issue. Yeah to_planar has a lot of questionable logic in there (also why the heck didn't I name it to_2D :). With regards to your ValueError, are the points flat? Some quick messing around:

In [7]: m = trimesh.creation.random_soup(10)

In [8]: m.outline()
Out[8]: <trimesh.Path3D(vertices.shape=(30, 3), len(entities)=10)>

In [9]: o = m.outline()

In [10]: o
Out[10]: <trimesh.Path3D(vertices.shape=(30, 3), len(entities)=10)>

In [13]: o.to_planar(to_2D=np.eye(4))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-13-f5fce6e80cc2> in <module>()
----> 1 o.to_planar(to_2D=np.eye(4))

/home/mikedh/trimesh/trimesh/path/path.py in to_planar(self, to_2D, normal, check)
    949             height = 0.0
    950             if check:
--> 951                 raise ValueError('points are not flat!')
    952         else:
    953             # if the points were planar store the height

ValueError: points are not flat!

In [15]: o.to_planar(to_2D=np.eye(4), check=False)
Out[15]: 
(<trimesh.Path2D(vertices.shape=(30, 2), len(entities)=10)>,
 array([[1., 0., 0., 0.],
        [0., 1., 0., 0.],
        [0., 0., 1., 0.],
        [0., 0., 0., 1.]]))

In [17]: trimesh.geometry.align_vectors([0,0,1], [0,0,1])
Out[17]: 
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

In [18]: trimesh.geometry.align_vectors([0,0,1], [0,0,1-1e-12])
Out[18]: 
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

In [19]: trimesh.geometry.align_vectors([0,0,1], [0,0,1+1e-12])
Out[19]: 
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.],
       [0., 0., 1., 0.],
       [0., 0., 0., 1.]])

In [20]: trimesh.geometry.align_vectors([0,0,-1], [0,0,1+1e-12])

Out[20]: 
array([[ 1.,  0.,  0.,  0.],
       [ 0., -1.,  0.,  0.],
       [ 0.,  0., -1.,  0.],
       [ 0.,  0.,  0.,  1.]])

In [21]: 

In [21]: 

In [21]: trimesh.geometry.plane_transform([0,102,0], [0,0,1])
Out[21]: 
array([[   1.,    0.,    0.,   -0.],
       [   0.,    1.,    0., -102.],
       [   0.,    0.,    1.,   -0.],
       [   0.,    0.,    0.,    1.]])
pmberkeley commented 4 years ago

The points are flat - they were produced by the section method as a work-around for the section_multiplane issues. I need the discrete (really just the properly-ordered) paths for regular-interval cross sections of my object, but the discrete method only works on 2D paths, and the section method only produces 3D paths! (The 2D paths produced by section_multiplane have the aforementioned issues with orientation).

The problem is a rotation problem, as far as I can tell. If I give it no transformation matrix, it works fine - it easily deduces a plane for the points (because they're created by the section method), but it shifts them away from their original location (which I need to know for my analysis). If I tell it the actual plane it should be in, or use the identity transformation matrix, it says they aren't in-plane. I also tried letting it find the plane and then just reversing the translation, but the translation is wrong! It seems to be swapping the signs and the translations of each axis that remains in the 2D path.

pmberkeley commented 4 years ago

I believe I've found it! The problems with the transformation matrices seem to be coming from geometry.plane_transform(origin, normal): https://github.com/pmberkeley/trimesh/blob/f0433cf6d886b1031d0d448daf5693df692d6f9f/trimesh/geometry.py#L14

The logic here assumes the normal is Z and the XY plane should be zeroed on the Z plane. This prevents people from using any other axes as the normal axis. I think an easy fix would be to have the desired/destination normal optionally specified for geometry.plane_transform and all methods that use it. Or, if a normal is specified by any parent methods, it could be assumed to be the desired normal. (This would fix the strange behavior with both section_multi-plane and to_planar, I believe - I'm about to test this out).

pmberkeley commented 4 years ago

I've been digging around in this for a while now, and I think I now (finally) understand what's going on. There ISN'T an issue with the code, it's just really, really confusing and under-explained about what it's giving you with mesh_multiplane and section_multiplane. It's stepping through and rotating each section into the XY plane, which also shifts it in an apparently (although not actually) unpredictable way.

I now think an added functions (which I may work out in a bit), plus documentation and examples, would help spare other people this confusion and related problems. For example, the expected result of a multiplane operation with the normal in the X direction would be for the points to map directly onto the YZ plane, but because they're rotated into the XY plane using an unexpected transformation matrix, you may end up with numbers that aren't expected (like negatives when everything is positive). This behavior is also leading to some errors in I think it was section_multiplane where it's recording "None" results for all negative mesh_multiplane operations.

I think I have a better idea now about how to clean this up, but have to get around to doing so. I'll pop back in if/when I make progress.

mikedh commented 4 years ago

Thanks for looking in to this in depth @pmberkeley! Hahah I'm gonna use "There ISN'T an issue with the code, it's just really, really confusing and under-explained about what it's giving you" as a product testimonial if we ever make trimesh t-shirts :laughing:

glitchland commented 4 years ago

Hrm, I just hit this issue too. Let me know how you get on with this @pmberkeley. Adjusting the mesh bounds a little and skipping 'None' types allowed me to work around the issue for now.

pmberkeley commented 4 years ago

Section_multiplane uses the origin as the start of the heights. So if your object goes from -5 to 5 on the x axis, you would set an origin at x=-5 (say, [-5,0,0], or whatever the lower bounds are), a normal of normal = [1,0,0], and then the heights array would be measured from that origin, NOT an array of the locations of the new origins. A heights array of heights = [0,2,4,6,8,10] would capture sections at [-5,-3,-1,1,3,5].

Then, the output of the section_multiplane will be 2D paths that have been rotated into the XY plane using the transformation matrices that are output by mesh_multiplane and cached. In order to get the section multiplane back into 3D space (from which you can get a projection onto some other plane, like the YZ plane, if that's what you want like I did), you have to use the to_3D method on each 2D section output by the section_multiplane method. If you don't put a specific transform into the method, it will pull the cached transforms to put it back in the original location. (@mikedh, is there some direct way to do this, or do you have to run the output of the section_multiplane method through a for loop like I did?)

This is my heights generator helper function; it's crude - it only works for major axis directions and it doesn't have an option for selecting n number of steps instead of setting step size. However, it should clarify which heights to choose for the methods to work.

def heights_generator(mesh,step,axis):
    """helper function to create heights appropriate for use in
    section_multiplane and mesh_multiplane

    Parameters
    _______________
    mesh : mesh object used in multiplane operations
    step : float, desired step size for sections
    axis: integer, 0, 1, or 2, for selecting axis of bounds for heights

    Returns
    ______________
    heights : array of heights suitable for use in multiplane methods
    """
    levels = np.arange(start=mesh.bounds[0][axis],
                       stop=mesh.bounds[1][axis] + step,
                       step=step)
    heights = levels + abs(mesh.bounds[0][axis])
    return heights
ranjanbsc commented 4 years ago

Hi @pmberkeley,

I am using trimesh and needed it for slicing meshes. I was trying to use this code given above (no transformations carried out):

srx = [m.section(plane_origin=m.bounds[0],
                plane_normal=[1,0,0])
      for h in np.arange(m.bounds[0,0],m.bounds[1,0],.1)]

After running the above code, I used "trimesh.Scene([srx,m]).show()" to visualize and check where the sections are created.

Every section is created only at one location, i.e., m.bounds[0] instead of at different h values (h is missing in the code "m.section(plane_origin=m.bounds[0],plane_normal=[1,0,0])").

Can you please check and confirm if this is the correct?

Ranjan

pmberkeley commented 4 years ago

@ranjanbsc , you haven't included enough context for me to be sure, but it looks like you're using m.section when you want m.section_multiplane. You're creating individual sections and then expecting Scene.show() to show them all, when they aren't all in the same object. (Or, you're expecting it to show at different levels, when the output of m.section isn't in 3D space. Not really sure which is happening).

ranjanbsc commented 4 years ago

Hi @pmberkeley ,

Yes, I need section multi-plane. I started out with section multi-plane. But in section multi-plane, we do not get 3D sections (I understand we can get back 3D sections using transformation matrix).

I was just trying to implement section tool itself in a loop and I started with this code. Now when I tried the above code (see below) in different orientation, the result was [None, None, None, None, None, None, None, None, None, None]:

sections_along_dirs_2 = [m.section(plane_origin=m.bounds[0], plane_normal=[0.5,0.9,0.1]) for h in np.arange(m.bounds[0,2],m.bounds[1,2],.1)]

So I started out with this and managed to get the sections in any direction. see below (i have added resulting sections and also image): [<trimesh.Path3D(vertices.shape=(7, 3), len(entities)=1)>, <trimesh.Path3D(vertices.shape=(7, 3), len(entities)=1)>, <trimesh.Path3D(vertices.shape=(7, 3), len(entities)=1)>, <trimesh.Path3D(vertices.shape=(8, 3), len(entities)=1)>, <trimesh.Path3D(vertices.shape=(8, 3), len(entities)=1)>, <trimesh.Path3D(vertices.shape=(8, 3), len(entities)=1)>, <trimesh.Path3D(vertices.shape=(8, 3), len(entities)=1)>, <trimesh.Path3D(vertices.shape=(9, 3), len(entities)=1)>, <trimesh.Path3D(vertices.shape=(7, 3), len(entities)=1)>, <trimesh.Path3D(vertices.shape=(7, 3), len(entities)=1)>]

image

I hope this provides the context.

pmberkeley commented 4 years ago

@ranjanbsc

You're trying to use trimesh to do a thing it isn't currently programmed to do. You can't show all the sections from a m.section call instead of a m.section_multiplane call, because it only "sees" one section at a time, I believe. And you're getting "None", because the way you set up the heights is wrong. As I said above, the heights for section_multiplane are supposed to be measured off the origin along the normal, not using absolute heights, as you did in your code. My heights-generator code is only meant for unit normal directions as an example for people, and for my own uses, but it's fairly trivial to create your own heights generator for whatever axis you want to look at: if you want a section every 1 unit, just do h = np.arange(0,n,1), where n is the total distance you need to find the slices for, starting from the origin you entered.

This is such a common source of confusion, and one that I took months to sort through myself, I think it might be worth it to have this at least clarified in the documentation, if it isn't already, or have the code changed so people have an option to put in the absolute heights, since this is the intuitive version of the algorithm for many, especially since m.section appears to work this way. But I'm not a maintainer of the package, nor do I have time right now to try to change the code, so we would have to turn to @mikedh or others for that.