Snapping multiple meshes together and extract transformation matrices

[ttsesm]

Hi @marcomusy,

I have the following problem that I am trying to address and I am trying to figure out how possibly I could automate the whole procedure. Imagine that I have multiple pieces of a complete object which are randomly given as input (different orientation, position, etc) and then I would like to find a way to automatize (not perfectly) how to assemble them all together to the final object and extract the transformation matrices. So imagine that I have the following 5 pieces:

which if you put them together in the correct order you should get the following complete object:

Currently someone could do that manually by using a corresponding 3d analysis tool, e.g. Meshlab, CloudCompare, Blender, Meshmixer, etc... and as I did. However, this takes a lot of time especially if you plan to do it for multiple objects and moreover the result still might not be the best. Thus, I wanted to ask you from your experience if you know any tool that could help me on that or if you believe I could do something with vedo.

My idea would be to extract some kind of boundary/shape contours and try to apply some kind of shape fitting metric or something similar but I am not sure what that could be. I've found your discussion here about shape decomposition but I am not sure whether this could be related or not. I've tried to apply and test with different aligning algorithms but these are not working properly since they look for similar features that overlay each other while in this case I am looking for features that complement each other instead.

Any idea is welcome.

p.s. actually even an easier interactive mode, where I can select whether two edges should snap together would be helpful in regards to the current solution where I am trying to bring two pieces close together manually.

tombstone.zip

[marcomusy]

I think the problem is non-trivial.. it depends on how much effort you wish to invest on it..

1. if you already now more or less how the pieces must face together you may refine this initial alignment using vedo. The initial allignement can also be done in vedo by pressing a in the scene to drag/rotate different meshes (pressing a again toggles back to normal). Once this first alignment is obtained you can refine it with align(rigid=True) having selected a region of the mesh that you already confident it should match the facing opposite. This is me playing with manual alignment :)

vedo https://www.dropbox.com/s/b3pa0y4v10stxne/tombstone_scene.npz

2. you may look for cross correlation of one small region of interest to the other shape to identify where the matching occurs (this is easy to do in 2d, not sure if a 3d version of it exists). Or you cam compute a metric of how close the template face is to the target (with distanceTo()).
3. I would not go for shape decomposition - I have the impression it is not what you're looking for.

[ttsesm]

Thanks Marco for the feedback and your time.

I think the problem is non-trivial.. it depends on how much effort you wish to invest on it..

Indeed, actually it is still an open problem in the community (especially with 3D data). Well I plan to spend quite some time on it, since it is related with my next research project. However, I want to have these solved beforehand in order to have some ground truth to compare with.

1. if you already now more or less how the pieces must face together you may refine this initial alignment using `vedo`. The initial allignement can also be done in `vedo` by pressing `a` in the scene to drag/rotate different meshes (pressing `a` again toggles back to normal). Once this first alignment is obtained you can refine it with `align(rigid=True)` having selected a region of the mesh that you already confident it should match the facing opposite.
   This is me playing with manual alignment :)

vedo https://www.dropbox.com/s/b3pa0y4v10stxne/tombstone_scene.npz

Interesting I was not aware that I could do that with vedo. This is quite similar with what I can do in Meshlab already, though the manipulator in vedo seems much more trivial. Thanks for that ;-).

Though this solution would still need some manual intervention.

2. you may look for cross correlation of one small region of interest to the other shape to identify where the matching occurs (this is [easy to do in 2d](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.correlate2d.html), not sure if a 3d version of it exists). Or you cam compute a metric of how close the template face is to the target (with `distanceTo()`).

Yup, this seems a relative fair, first naive approach. For example I was thinking to extract the boundary lines (or big contours) of the mesh and try to do some shape fitting based on some distance metric, e.g. chamfer, hausdorf distance. I need to see though what the boundary/controur would give me since I would be interested to extract mainly the upper and bottom boundaries. I will try to work a bit on this and update.

3. I would not go for shape decomposition - I have the impression it is not what you're looking for.

I see, fair enough.

[marcomusy]

By pressing a you can toggle pan (shift-drag) rotate (ctrl-drag) any mesh, which is indeed quite convenient (I've used it a lot!).

I was thinking to extract the boundary lines (or big contours) of the mesh and try to do some shape fitting based on some distance metric

I would suggest not to :) beacuse the boundaries not always match. e.g. The above part does while the lower doesn't.. What I would try out is to pick a small region of the part that you already guess it's a reasonable match and work out a metric on that bit only.. (e.g. you may histogram the point to point distances and measure the peak at zero indicating a good match).

To get the boundaries use boundaries(), to fit a plane to a cloud of pts use fitPlane().

[ttsesm]

Hhmm :thinking: , boundaries() give me the following output for the pieces above: and none for the pieces below:

I tried to get only the largest region by mesh.boundaries().extractLargestRegion().c('red') but apparently something is messed up. Any idea how to workaround this?

[ttsesm]

Ok, I guess I will need to play with the angle or map it flat on one of the planes.

Btw, how do I get the Mesh() dimensions/size (I couldn't find any method returning this).

[ttsesm]

Btw, how do I get the Mesh() dimensions/size (I couldn't find any method returning this).

Ok got it, it is the bounds() method.

[marcomusy]

yes , you need to set a feature angle to identify boundaries.

bounds() gives the xyz bounding box of any vedo object.

[ttsesm]

yes , you need to set a feature angle to identify boundaries.

I do not think it will work, since my 3d models are watertight. Thus, this implies that the mesh outline would be an empty path. Actually, I played with different angles as well as transforming the meshes to be parallel to the plane but didn't manage to get any better result from the one shown above.

I found also this solution which seems interesting. Depending on the perspective/view angle you can retrieve the corresponding outline.

[marcomusy]

yes , you need to set a feature angle to identify boundaries.

I do not think it will work, since my 3d models are watertight. Thus, this implies that the mesh outline would be an empty path. Actually, I played with different angles as well as transforming the meshes to be parallel to the plane but didn't manage to get any better result from the one shown above.

I told you it would not work ... :)

I found also this solution which seems interesting. Depending on the perspective/view angle you can retrieve the corresponding outline.

such a thing is completely trivial in vedo check out: examples/basic/silhouette*.py But I don't think that it could turn out useful for your specific problem.

[ttsesm]

Hi Marco,

Thinking and experimenting about the current issue I wanted to test something else. However, for this I would need to create the missing part (kind of vessel) of each mesh object. How I could possibly do that with vedo? My idea is to extract the bounding box of the object as a separate mesh, then create the volume of both meshes, subtract the two volumes, and finally meshify the output. Would this work?

[marcomusy]

not sure if it would work ... but here you have a few ingredients that might be useful for your tasty recipe :) mesh.box() to return the mesh box as a new mesh mesh2Volume() convert a mesh it into a Volume where the foreground (exterior) voxels value is 1 and the background volume.isosurface() to get back the outer boundary

[ttsesm]

@marcomusy can I retrieve the mesh_file.printInfo() output as a dictionary or something?

[marcomusy]

no.. that is just a printout.. but you can access the same info from the object itself.

[ttsesm]

no.. that is just a printout.. but you can access the same info from the object itself.

I see, yes I know. I was just hopping that could retrieve them all together with a method like in printInfo()/print().

[ttsesm]

By pressing a you can toggle pan (shift-drag) rotate (ctrl-drag) any mesh, which is indeed quite convenient (I've used it a lot!).

Ciao Marco, is there a shortcut for locking pan/rotate only in one axis each time?

[marcomusy]

Hi! Not really .. vtk offers a number of standard ways to interact with the scene, you can go thourgh them by setting show(mode=0,1,2,3..) one that might be useful to you is maybe mode="image" or "2d".

[ttsesm]

Hi! Not really .. vtk offers a number of standard ways to interact with the scene, you can go thourgh them by setting show(mode=0,1,2,3..) one that might be useful to you is maybe mode="image" or "2d".

I see, thanks a lot. I will play with the modes options and check if any is useful.

[ttsesm]

To give a small update, recently I've came upon this manual tool http://vcg.isti.cnr.it/~pietroni/reassembly/index.html for snapping multiple pieces together. It is from the same guys that are running meshlab, though it runs quite nice they never did any further polishing for further usage optimization. In principle you just give some correspondence points on both two meshes that you want to merge and then they run a solver to find the optimal alignment. The more points you give the more accurate snapping you get. There is also a corresponding publication describing the details.

It would be interesting if something similar could be done with the interactive parts of vedo and if possibly could be further enriched with more constraint options. E.g. possibly to mark more than one point at a time with some kind of scribbles usage instead of individual points through a brush usage operation or something similar.

In any case it is a good reference, to be considered.

[marcomusy]

looks cool - thanks for the update :)

It would be interesting if something similar could be done with the interactive parts of vedo and if possibly could be further enriched with more constraint options. E.g. possibly to mark more than one point at a time with some kind of scribbles usage instead of individual points through a brush usage operation or something similar.

Yes! you can "brush" the surface to determine the areas that should match (without even considering the order of the points picked by the brush). See e.g. examples/basic/mousehover2.py

I think the key to the solution is to write down a clever metric - not the trivial L2 - that favors the matching of small areas where the matching is indeed already almost perfect and discard outliers, without trying to average out over larger regions. Maybe (i'm thinking aloud) you may even draw the vector differences of the surfaces in a different R^3 parameter space and use tools like removeOutliers() in this space..

[ttsesm]

I think the key to the solution is to write down a clever metric - not the trivial L2 - that favors the matching of small areas where the matching is indeed already almost perfect and discard outliers, without trying to average out over larger regions. Maybe (i'm thinking aloud) you may even draw the vector differences of the surfaces in a different R^3 parameter space and use tools like removeOutliers() in this space..

Actually in the paper they describe an enerrgy minimization metric which could be adapted. Soon I will have a better look to understand the implementation since it seems interesting and simple to be adapted. After all, I want to implement an automatic solution as well and this looks a nice starting point.

[ttsesm]

Hi @marcomusy, I found some time to have a look on the aforementioned paper above. As I pointed to you also earlier In the section 3.2 they mention about the metric that they use in order to minimize the energy and extract the transformation matrices. They also mention about different constraint pipelines but in the paper only the adjacency one is implemented. As you can see in the provided formulation

they use this k stiffness factor as they call it. However, it is not clear to me how this factor prevents let's say the two pieces to overlap each other instead of snapping them together. Thus, I would like to ask you whether you get a better understanding on how this is implemented. I would be interested to re-implement their solver.

[marcomusy]

I'll think about it and let you know if I get any ideas :)

[marcomusy]

some semi random thoughts:

• I would forget about a completely automated process at least at the very first iteration..
• assume for the user the need to mark some areas that are meant to overlap and then explore the parameter space with a metric (what they in fact do in the paper..)
• a possible metric... why not playing with the normals too , in order to avoid mesh penetration one can add a repulsive term:

  np and nq are the normals (in purple, I draw them inverted for surface point q).

probably E won't work exactly as I wrote but it's to give a gist of what it may look like.

• another possible metric could be made by looking at the histogram of the (signed) distances of one mesh to the other, maybe on a subset of points (so it's fast).
• i would create a small exploratory application where you can paint-brush some areas and drag the surfaces to see how your metric behaves.. you can do this easily in vedo and have the advantage that you will almost certainly reuse it in your final code.

hope it helps!

[ttsesm]

Thanks Marco, yes using the normals was in my intentions as well. Let me process a bit your thoughts and your suggestions and I will try to create a small demo so we can play with.

[ttsesm]

@marcomusy if q and p are my points (or set of points) in the formulation of E above what are the A, B and C?

[marcomusy]

simple numeric constants, for A and B only their ratio matters so you can treat them as dimensionless , not so for C which should be scaled by the object units

[ttsesm]

Ciao Marco,

Ok, I've created a small script where I am extracting some correspondences from two pieces. I do not care about the method for the moment (I am doing it manually at the moment with the closestPoint() method but this can be done in different ways as you also pointed out) since I just want to test the solver. The snippet is the following:

from vedo import *
import numpy as np

def find_correspondences(meshes):

    m1_pts = meshes[0].points()
    m2_pts = meshes[1].points()

    correspondences = []

    for i, p in enumerate(m1_pts):
        iclos = meshes[1].closestPoint(p, returnPointId=True)
        correspondences.append([i, iclos])

    dist_thres = 0.05
    correspondences = np.asarray(correspondences).reshape(-1, 2)
    correspondences = correspondences[np.where(mag(m1_pts[correspondences[:,0]] - m2_pts[correspondences[:,1]]) < dist_thres), :]

    return correspondences.squeeze()

def transform_mesh(m):
    cm = m.centerOfMass()
    m.shift(-cm)
    elli = pcaEllipsoid(m, pvalue=0.5)

    ax1 = versor(elli.axis1)
    ax2 = versor(elli.axis2)
    ax3 = versor(elli.axis3)

    T = np.array([ax1, ax2, ax3])  # the transposed matrix is already the inverse

    return m.applyTransform(T, reset=True)

def main():
    files = ['./cake_part01.xyz', './cake_part02.xyz']

    mesh_files = []
    mesh_models = []

    for i, file in enumerate(files):
        m = Mesh(file).color(i)
        mesh_files.append(m)
        mesh = transform_mesh(m.clone())
        mesh_models.append(mesh)

    corrs = find_correspondences(mesh_files)

    pts1 = Points(mesh_files[0].points()[list(corrs[:,0])], r=5, c='blue')
    pts2 = Points(mesh_files[1].points()[list(corrs[:, 1])], r=5, c='red')

    pts1_ = Points(mesh_models[0].points()[list(corrs[:, 0])], r=5, c='blue')
    pts2_ = Points(mesh_models[1].points()[list(corrs[:, 1])], r=5, c='red')

    show(mesh_files, pts1, pts2, "Original", at=0, N=2, axes=1, sharecam=False)
    show(mesh_models, pts1_, pts2_, "Testing", at=1, interactive=True, sharecam=False).close()

    return 0

if __name__ == "__main__":
    main()

So now I have my matching points from the two pieces and I want to go from the "Testing" window back to "Original" let's say. In first place before including normals etc. I would like to test minimization metric from the paper. Any suggestion how I could it given the matching points? I mean I do not get how to obtain the rigid transformations R and T for each set of points minimize the distance and apply the stiffness factor.

cake_pieces.zip

[marcomusy]

Sorry for the late reply... Easter holidays :) Why do you work with pointclouds instead of Meshes? Do you have normals? In the case of the tombstone I imagined as an initial alignment something like:

"""Press 1 or 2 to select mesh to paint
Press spacebar to release selection
Press Z to align"""
from vedo import settings, Mesh, Text2D, Points, Plotter

settings.enableDefaultKeyboardCallbacks = False

def func(event):  # callback function
    p = event.picked3d
    if p is None:
        return
    if event.keyPressed not in ['1','2']:
        return
    i = int(event.keyPressed)-1
    msh = meshes[i]
    if event.actor != msh:
        return
    ids = msh.closestPoint(p, N=150, returnPointId=True)
    match[i].extend(ids)
    pids = list(set(match[i])) # make them unique
    plt.remove("mypoints")
    pts = Points(msh.points()[pids], r=6)
    pts.name = "mypoints"
    plt.add(pts)
    # print(len(match[0]), len(match[1]))

def key(event):
    if event.keyPressed=="1":
        txt.text("Selected mesh #1")
        plt.render()
    elif event.keyPressed=="2":
        txt.text("Selected mesh #2")
        plt.render()
    elif event.keyPressed=="Z":
        ids0 = list(set(match[0])) # make them unique
        ids1 = list(set(match[1]))
        pts0 = meshes_pts[0][ids0]
        pts1 = meshes_pts[1][ids1]
        p1 = Points(pts1).alignTo(Points(pts0), rigid=True)
        msh1.applyTransform(p1.transform, reset=True)
        plt.remove("mypoints").render()
        # compute here some metrics
        # ...
    elif event.keyPressed=="q":
        # plt.interactor.ExitCallback() # add in case you get a seg-fault
        plt.close()
    else:
        # disable picking
        txt.text(__doc__)
        plt.render()

msh0 = Mesh('data/tombstone/Tombstone1_low.obj').c(0).print()
msh1 = Mesh('data/tombstone/Tombstone2_low.vtk').c(1)
meshes = [msh0, msh1]
meshes_pts = [msh0.points(), msh1.points()]
match = [[],[]]

txt = Text2D(__doc__, bg='yellow', font='Calco')

plt = Plotter(axes=1)
plt.addCallback('mouse hover', func)
plt.addCallback('key press', key)
plt.show(msh0, msh1, txt)
plt.close()

Tombstone.zip

Next steps:

• implement E as above and see how it behaves when shfiting the matching points around (or some mesh region close to those)
• minimize E

Hope this helps!

[ttsesm]

Thanks Marco, no worries for the late response. I will have a look and check on the code and update you. It will definitely help.

Point clouds or meshes should be the same, it shouldn't make any difference. Yes, normals will be available and plan to use.

[marcomusy]

Any progress on this? I'll close the issue by now, let me know if you need any help.

[ttsesm]

Any progress on this? I'll close the issue by now, let me know if you need any help.

Sure, no worries. I can write here in any case that I need any help or to update you. For the moment I am working on a data driven approach since anything else that I've tried was not possible to generalize easily. Thus, I am checking on some data driven solutions and possibly adapt them for my use cases.

In case you are interested and you have some time, you can have some look in these two works https://neural-shape-mating.github.io/ and https://breaking-bad-dataset.github.io/ they look promising.