Help with svd shape matching using numpy

[jkunz07]

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 !

[petz]

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

[jkunz07]

Yes, it works !

 

Thanks for your help

[RayFoundry]

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?