ptsArray = ptsArray.reshape(len(ptsArray)/3 ,3) # reshape as array of vectors

ptsArray = ptsArray.reshape(len(ptsArray)//3 ,3) # reshape as array of vectors

import numpy as np

node = hou.pwd()
geo = node.geometry()

ptsString = geo.pointFloatAttribValuesAsString("P") # get pts data efficiently
ptsArray = np.frombuffer(ptsString, dtype=np.float32) # convert to np array
ptsArray = ptsArray.reshape(len(ptsArray)//3 ,3) # reshape as array of vectors

mean = np.mean(ptsArray, axis=0) # compute mean of pts positions
ptsArray = ptsArray - mean # transform pts to the center
mtx = np.cov(ptsArray.T) # compute covariation mtx

eigenvalues, eigenvectors = np.linalg.eig(mtx) # compute eigens

# sort eigenvalues
sort = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[sort]
eigenvectors = eigenvectors[:,sort]

e1 = eigenvectors[0]
    e2 = eigenvectors[1]
e3 = eigenvectors[2]

# compute bbox dimesions
orient = np.array( (e1, e2, e3) ).T

ee1 = np.dot((1,0,0), orient)
ee2 = np.dot((0,1,0), orient)
ee3 = np.dot((0,0,1), orient)

scale1 = np.dot(ptsArray, ee1)
scale2 = np.dot(ptsArray, ee2)
scale3 = np.dot(ptsArray, ee3)

scale1 = ( np.amax(scale1) - np.amin(scale1) ) / 2
scale2 = ( np.amax(scale2) - np.amin(scale2) ) / 2
scale3 = ( np.amax(scale3) - np.amin(scale3) ) / 2

scale = (scale1, scale2, scale3)

# output to attribs
geo.addAttrib(hou.attribType.Global, "eigenvector1", (0.,0.,0.))
geo.addAttrib(hou.attribType.Global, "eigenvector2", (0.,0.,0.))
geo.addAttrib(hou.attribType.Global, "eigenvector3", (0.,0.,0.))
geo.addAttrib(hou.attribType.Global, "scale", (0.,0.,0.))
geo.addAttrib(hou.attribType.Global, "mean", (0.,0.,0.))

geo.setGlobalAttribValue("eigenvector1", e1)
geo.setGlobalAttribValue("eigenvector2", e2)
geo.setGlobalAttribValue("eigenvector3", e3)
geo.setGlobalAttribValue("scale", scale)
geo.setGlobalAttribValue("mean", (float(mean[0]), float(mean[1]), float(mean[2])))

    e2 = eigenvectors[1]