Jump to content
Sign in to follow this  
jkunz07

Help with svd shape matching using numpy

Recommended Posts

Hey, I'm trying to match the translation and rotation of another set of points using the svd function in numpy

 

The code I have right now is below:

 (Some of it I got from here http://nghiaho.com/?page_id=671)

import numpy, math

node = hou.pwd()
geo_A = node.geometry()
if node.inputs()[1] != None and len(node.inputs()) >= 2:
    geo_B = node.inputs()[1].geometry()

string_buffer = geo_A.pointFloatAttribValuesAsString("P", float_type=hou.numericData.Float32) #get array of all points attr "P" as string buffer
array_of_pos_A = numpy.fromstring(string_buffer, dtype='float32' ).reshape(-1, 3)   # reshape to make it 2dimentional (array of vectors)
string_buffer = geo_B.pointFloatAttribValuesAsString("P", float_type=hou.numericData.Float32)
array_of_pos_B = numpy.fromstring(string_buffer, dtype='float32' ).reshape(-1, 3)
N = array_of_pos_A.shape[0]; # total number of points

min_bbox_A = numpy.min(array_of_pos_A, axis=0)
max_bbox_A = numpy.max(array_of_pos_A, axis=0)
centroid_A = (min_bbox_A+max_bbox_A) / 2

min_bbox_B = numpy.min(array_of_pos_B, axis=0)
max_bbox_B = numpy.max(array_of_pos_B, axis=0)
centroid_B = (min_bbox_A+max_bbox_ / 2

AA = array_of_pos_A - numpy.tile(centroid_A, (N, 1))
BB = array_of_pos_B - numpy.tile(centroid_B, (N, 1))

H = numpy.dot( numpy.transpose(AA), BB)
U, S, Vt = numpy.linalg.svd(H)
R = Vt.T*U.T

# special reflection case
if numpy.linalg.det(R) < 0:
   Vt[2,:] *= -1
   R = Vt.T * U.T

t = -R*centroid_A.T + centroid_B.T


print 't :  '+str(t)+'\n\n'
print 'R :  '+str(R)+'\n\n'

geo_A.setPointFloatAttribValuesFromString("P", numpy.getbuffer(AA), float_type=hou.numericData.Float32)   # Put the values back to the geo point position attr

I can see 't' and 'R' matrices print out, but any attempts of actually transforming the first input geometry by them have resulted in broadcast errors.

Additionally, any attempts to output them as geometry attributes have failed as well. 

 

numpy_svd_test_v002.hip

 

Any help here would be greatly appreciated !

Share this post


Link to post
Share on other sites

that should do the trick ...

node = hou.pwd()
geo = node.geometry()
geoRef = hou.pwd().inputs()[1].geometry()

import numpy as np

posArray = np.asarray(geo.pointFloatAttribValues("P")).reshape(-1, 3)
posArrayRef = np.asarray(geoRef.pointFloatAttribValues("P")).reshape(-1, 3)

centroidRef = np.mean(posArrayRef, axis = 0)
centroid = np.mean(posArray, axis = 0)

posArrayRef -= centroidRef

covar = np.dot(posArray.T, posArrayRef)

u, s, v = np.linalg.svd(covar)

rotMat = np.dot(u, v)

mat = hou.hmath.buildTranslate(-centroid)
mat *= hou.Matrix4(hou.Matrix3(rotMat))
mat *= hou.hmath.buildTranslate(centroidRef)
geo.transform(mat)

hth. petz

  • Like 4

Share this post


Link to post
Share on other sites
On 27.9.2015 at 5:28 PM, petz said:

that should do the trick ...


node = hou.pwd()
geo = node.geometry()
geoRef = hou.pwd().inputs()[1].geometry()

import numpy as np

posArray = np.asarray(geo.pointFloatAttribValues("P")).reshape(-1, 3)
posArrayRef = np.asarray(geoRef.pointFloatAttribValues("P")).reshape(-1, 3)

centroidRef = np.mean(posArrayRef, axis = 0)
centroid = np.mean(posArray, axis = 0)

posArrayRef -= centroidRef

covar = np.dot(posArray.T, posArrayRef)

u, s, v = np.linalg.svd(covar)

rotMat = np.dot(u, v)

mat = hou.hmath.buildTranslate(-centroid)
mat *= hou.Matrix4(hou.Matrix3(rotMat))
mat *= hou.hmath.buildTranslate(centroidRef)
geo.transform(mat)

hth. petz

Very useful, thank you!

This does only work when the point count on both geometries is the same though. Any idea how to extend this for cases where the point count is not the same?

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now
Sign in to follow this  

×