気泡の長軸・短軸の算出

import numpy as np
import vtk

def get_inertia_matrix(xyz):
    # Moment of inertia tensor
    x,y,z = xyz.T
    Ixx = np.sum(y**2 + z**2)
    Iyy = np.sum(x**2 + z**2)
    Izz = np.sum(x**2 + y**2)
    Ixy = -np.sum(x * y)
    Iyz = -np.sum(y * z)
    Ixz = -np.sum(x * z)
    I = np.array([[Ixx,Ixy,Ixz], [Ixy,Iyy,Iyz], [Ixz,Iyz,Izz]])
    return I


def get_axes(I):
    # Calculate the eigenvectors.
    eig = np.linalg.eig(I)

    # Sort by order of increasing eigenvalues
    eig_ids = np.argsort(eig[0])
    eig_vals = eig[0][eig_ids]
    eig_vecs = eig[1][eig_ids]
    return eig_vals, eig_vecs


input0 = inputs[0]

pnts = input0.GetPoints()
length = input0.GetLength()

# Calculate the center of gravity point and then
# move the point cloud of the bubble to the origin.
com = np.mean(pnts, axis=0)
print(f'com {com}')
xyz = pnts - com

# Prepared to find the intersection of the bubble and the axis
obbtree = vtk.vtkOBBTree()
obbtree.SetDataSet(self.GetInput())
obbtree.BuildLocator()

# Assuming the bubble to be ellipsoid-like, obtain the major and minor axes.
I = get_inertia_matrix(xyz)
eig_vals, eig_vecs = get_axes(I)
print('eig_vecs:')
print(eig_vecs)

axes = vtk.vtkAppendPolyData()
axis_ids = []
for i in range(3):
    # Determine the length of the axis by taking the 
    # intersection of the axial line and the bubble
    pnt = eig_vecs[i] * length + com

    intersects = vtk.vtkPoints()
    obbtree.IntersectWithLine(com, pnt, intersects, None)
    intersect = np.array(intersects.GetPoint(0))
    r = np.linalg.norm(intersect-com)
    print('r ', r)

    # Create an axis
    pnt1 = eig_vecs[i] * r
    pnt2 = -1.0 * pnt1
    line = vtk.vtkLineSource()
    line.SetPoint1(pnt1+com)
    line.SetPoint2(pnt2+com)
    line.Update()
    axis = line.GetOutput()

    axes.AddInputData(axis)
    axis_ids.append(i)

axes.Update()

output.DeepCopy(axes.GetOutput())
output.CellData.append(np.array(axis_ids), 'axis_id')