Reconstructing Shapes from Mean + Eigenvectors — How to Generate the Mesh?

[pstrimber]

Hi everyone,

I’m trying to reconstruct my shapes using the mean and the eigenvectors with the formula
reconstructed = mean_particles + score[0] * eigenvector[0] + score[1] * eigenvector[1] + …

My model contains 256 particles, so each eigenvector is 256 × 3 (X, Y, Z for every particle). The result reconstructed is therefore also a 256 × 3 array, which I take to be the particle positions of this new mesh.

What I’m unsure about:
How do I convert this particle matrix into a full mesh (OBJ or similar)?

I’d also like to confirm that:
score[i] is the value exported as “PCA component scores” for each image.
eigenvector[i] is the eigenvector exported in “Eigenvectors”, corresponding to mode i.

Any clarification or practical example would be greatly appreciated.

Thanks a lot!
Paula

[amorris]

To reconstruct a full mesh from a set of particles, you can use the MeshWarper to take a reference mesh + particles and warp it to a new set of particles.

You can follow this example:

https://github.com/SCIInstitute/ShapeWorks/blob/master/Testing/PythonTests/warp.py

def warpMesh():

  reference_mesh = Mesh(os.environ["DATA"] + "/mesh_warp/mesh_warp2.vtk")

  particles = ParticleSystem([os.environ["DATA"] + "/mesh_warp/mesh_warp2.particles"]).ShapeAsPointSet(0)

  warper = MeshWarper()

  warper.generateWarp(reference_mesh, particles)

  # Warp the mesh using the same particles as the reference --> warped mesh should be close to the reference mesh

  warped_mesh = warper.buildMesh(particles)

[pstrimber]

Hi Alan,

I hope you’re doing very well. Thank you so much for your previous response—it was extremely helpful. We have been able to make progress and reconstruct the images.
When reconstructing different images from our dataset using the code you provided and the formula reconstructed = mean_particles + score[0] * eigenvector[0] + score[1] * eigenvector[1] + …, we noticed something strange. The image reconstructed from the modes extracted from the model looked very different from the image reconstruction we see in the ShapeWorks desktop application. Moreover, it is noticeably different from the original image, which leads us to suspect that this Python-based reconstruction is not accurate, or that for some reason the modes exported from the model are corrupted.

Below I attach two examples. The first pair of images is the reconstruction using the code. The second pair of images is the reconstruction by ShapeWorks desktop. Finally, I’ve also attached the actual images after grooming. You can see that the reconstructions from ShapeWorks desktop are much closer to the real images.

Thank you in advance for your willingness and help, this is proving vital for continuing our research.

[amorris]

For the mesh warping, are you using the same reference template mesh in your code as Studio is using?

[pstrimber]

Hello Alan, thank you for your response.

So far we have tried two approaches. First, we exported the mean particles as shown in the first image. Using those particles as our reference mesh and the modes, we reconstructed the image with the formula we discussed earlier.

Then, as shown in the second image, we used the median particles as reference mesh—but again without success.

I’m attaching the results of using the mean and the median as reference meshes, respectively. This time I’ve attached photos of only a single bone, the one of the left in the previous question.

What do you think might be going wrong? Also, is it correct to use the exported PCA components directly as the modes, or should we apply a square root (or other transformation) to them? We read about PCA on the ShapeWorks website and we find this:

Lambda - This shows the PCA loading of current position of the slider. The middle of the slider, at the mean value, will be 0. The extent of lambda is defined by the number of standard deviations of the slider as described above. At standard deviation of 1.0, it will be the square root of the mode’s eigenvalue.

Should we raise the values of the PCA components score to the power of two?

FIRST OPTION:

SECOND OPTION:

STUDIO RECONSTRUCTION:

[pstrimber]

Hi, I’d also like to provide the code I’m using. It covers the process of reconstructing the new particles and then warping them to generate a VTK file.

import shapeworks as sw
import pyvista as pv
import numpy as np
import os

# --------- Creating file for new mesh particles ((256, 3)) ------------

# Mean particles file loading
particles = np.loadtxt(r'C:\Users\Paula\Desktop\ARCHIVOS INVESTIGACION\Shapeworks_modelos\Etiqueta_1\creacion_modelo_base\mean_particles_2.particles')

# Eigenvectors folder path
folder = r"C:\Users\Paula\Desktop\ARCHIVOS INVESTIGACION\Shapeworks_modelos\Etiqueta_1\creacion_modelo_base\eigenvectors"

# Quantity of modes to use for reconstruction
n_modes = 64 

# PCA Reconstructed image values (extracted from PCA component scores)
scores = np.array([-8.44069,-4.13383,-2.72596,-0.314619,-5.89917,-3.7897,4.76785,0.977776,-0.00548887,4.79821,-0.152873,-0.62478,0.15278,-0.266,-0.0415944,-0.6508,0.0303715,-1.36541,0.958535,0.468263,-0.675527,-0.59803,0.702887,-0.457504,-0.0244962,0.174829,-0.904587,1.49348,-0.110434,-0.26424,-0.85749,-0.175538,0.447883,-0.391096,1.21641,-1.2457,1.03426,-0.112783,0.179361,-0.821864,-0.848623,-0.55011,0.200428,0.00659176,0.430744,1.02361,-0.71389,0.295296,0.0974082,0.318745,-0.356043,0.026103,-0.0132103,0.217473,0.0640542,-0.123777,-0.322119,-0.512522,-0.369095,0.069574,-0.0316704,0.111362,-0.0652447,0], dtype=float)     # (n_modes,)

# Name for recreated image
name = "imagen1_reconstruida6"

# Path for reocntructed image particles
reconstructed_particles_path = r"C:\Users\Paula\Desktop\ARCHIVOS INVESTIGACION\Reconstrucciones\\"+name+".txt"

# Loading eigenvectors
eigenvecs = []
for i in range(n_modes):
    vec = np.loadtxt(f"{folder}\\eigenvector{i}.eval")
    vec = vec.reshape(256, 3)          # assures (256,3)
    eigenvecs.append(vec)
eigenvecs = np.stack(eigenvecs, axis=0)  # (n_modes, 256, 3)


# Applying formula: new mesh = media particles + sum(scores * eigenvectors)
delta = (scores[:, None, None] * eigenvecs).sum(axis=0)   # (256,3)
reconstructed = particles + delta      

# Estas particulas se usan en el warp para reconstruir la malla
np.savetxt(reconstructed_particles_path, reconstructed, fmt='%.6f', delimiter=' ')

# ------ Creating new vtk file with reconstructed mesh particles -----

DATA = r"C:\Users\Paula\Desktop\ARCHIVOS INVESTIGACION\Reconstrucciones"
os.environ["DATA"] = DATA

def warpMesh():

  #malla de referencia exportada desde Shapeworks
  reference_mesh = sw.Mesh(os.environ["DATA"] + "/mesh_warp/mean_mesh.vtk")

  # Particulas generadas desde el archivo reconstrucción imagenes
  particles = sw.ParticleSystem([reconstructed_particles_path]).ShapeAsPointSet(0)

  warper = sw.MeshWarper()

  warper.generateWarp(reference_mesh, particles)

  # Warp the mesh using the same particles as the reference --> warped mesh should be close to the reference mesh

  warped_mesh = warper.buildMesh(particles)

  warped_mesh.write(os.environ["DATA"] + "//" + name + ".vtk")

warpMesh()

print("done")

[amorris]

I would suggest removing the surface reconstruction from the equation. The mean shape particles and modes should match first. For example, if you have 512 particles and you are looking at +1SD of the first mode, when you perform that multiplication, then the particle positions should closely match from Studio and your Python. Have you confirmed this first?