# 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 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

• 4

Yes, it works !

##### 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?

## Create an account

Register a new account