3D 体积与其主轴对齐

问题描述 投票：0回答：1

我正在尝试将 3D 

numpy
阵列中包含的矩形棱柱与阵列的主轴对齐。具体来说，我需要通过旋转整个阵列来完成这种对齐，而不是提取对象的坐标，然后旋转这些点（这是因为从背景中完美分割对象对于我的应用程序来说是不现实的）。

我有一种在非常特殊的情况下有效的方法，但考虑到该方法的普遍性，它似乎出人意料地敏感。简而言之，我的方法仅适用于非常特定的对象方向。

我正在寻求有关 (a) 我的方法有什么问题或 (b) 实现相同目标的不同方法的指导。

我的方法（受到这篇文章和这篇文章的启发）：

from scipy.ndimage import affine_transform
from skimage import measure

def find_orientation(prism):
    # Calculate second-order moments
    M = measure.moments_central(prism, order=2)

    # Constructing the covariance matrix from the central moments
    cov_matrix = np.array(
        [
            [M[2, 0, 0], M[1, 1, 0], M[1, 0, 1]],
            [M[1, 1, 0], M[0, 2, 0], M[0, 1, 1]],
            [M[1, 0, 1], M[0, 1, 1], M[0, 0, 2]],
        ]
    )

    # Compute eigenvectors from the covariance matrix
    eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)

    return eigenvectors


def rotate_to_align(prism):
    # Calculate orientation matrix using moments method
    orientation_matrix = find_orientation(prism)

    # Rotate the prism using the orientation matrix
    aligned_prism = affine_transform(prism, orientation_matrix.T, order=3)

    return aligned_prism


aligned_array = rotate_to_align(misaligned_array)

请参阅此笔记本了解完整详细信息（对于SO来说太长了）： https://github.com/petermattia/3d-alignment/blob/main/prism_alignment.ipynb

相关但不同的问题/资源：

谢谢！

image-processing linear-algebra image-rotation rotational-matrices scipy.ndimage
1个回答
0
投票

好的，这样你就有了体积，包含一个实心盒子，以某种方式定向。您正在寻找该物体的轴。

根据力矩计算轴似乎不稳定。

我建议使用对象每个面的表面法线。

您可以通过对体积的梯度进行聚类来获得这些值。这将为您提供六个稳定的表面法线簇。

配对相反的向量。合理地对它们进行排序，以减少意外的轮换。然后应用一些更多的线性代数来确保您实际上拥有彼此正交的向量。如果盒子碰巧有一点剪切或者只是由于有限算术引起的错误，则不必是这种情况。

这并不完美。插值/采样和有限算术是问题。

# imports
import numpy as np
from numpy.linalg import norm, inv
import cv2 as cv
import matplotlib.pyplot as plt
import scipy.ndimage as ndi

# first bunch of utility functions

def normalize(vec):
    vec = np.asarray(vec)
    return vec / norm(vec)

def translate3(tx=0, ty=0, tz=0):
    T = np.eye(4)
    T[:3, 3] = tx, ty, tz
    return T

def rotate3_axis(axis, radians=None, degrees=None):
    axis = np.asarray(axis, dtype=np.float64)

    assert (radians is not None) ^ (degrees is not None)

    if degrees is not None:
        radians = np.deg2rad(degrees)

    if radians is not None:
        axis = normalize(axis) * radians
    
    R, jac = cv.Rodrigues(axis)

    T = np.eye(4)
    T[:3, :3] = R
    return T

def rotate3_R(R):
    T = np.eye(4)
    T[:3, :3] = R
    return T

def rotate3(axis, **kwargs):
    if isinstance(axis, np.ndarray) and axis.ndim == 2:
        return rotate3_R(R=axis)
    else:
        return rotate3_axis(axis, **kwargs)
        

def scale3(s=1, sx=1, sy=1, sz=1):
    return np.diag([s*sx, s*sy, s*sz, 1.0])

axis_names  = ['x',  'y',  'z']
plane_names = ['yz', 'xz', 'xy']

def show(volume, gradient=False):
    vmin = volume.min()
    vmax = volume.max()

    fig, axes = plt.subplots(figsize=(10, 10), nrows=1, ncols=3)

    for k, axis in enumerate(axes):
        axis_slice = np.take(volume, indices=volume.shape[k] // 2, axis=k)
        axis.set_title(f"{axis_names[k]}-axis, {plane_names[k]}-plane")
        axis.set_xlabel(plane_names[k][1])
        axis.set_ylabel(plane_names[k][0])
        axis.imshow(axis_slice, cmap="gray", vmin=vmin, vmax=vmax)
    plt.tight_layout()
    plt.show()

# phantom

vsize = np.array([50, 50, 50])

# box with different sizes in each dimension
box_size = np.array([40, 20, 10], dtype=int)
box_p0 = (vsize - box_size) // 2
box_p1 = box_p0 + box_size

volume = np.zeros(vsize, dtype=np.float32)

(bx0, by0, bz0), (bx1, by1, bz1) = box_p0, box_p1
volume[bx0:bx1, by0:by1, bz0:bz1] = 1


Tt = translate3(*(vsize/2))
Tr = rotate3(axis=(0,0,1), degrees=30)
T = Tt @ Tr @ inv(Tt)

warped = ndi.affine_transform(
    input=volume,
    matrix=T,
    # offset ignored
    output_shape=vsize,
    order=1,
    mode='nearest',
)


# per-voxel gradient of the warped volume

def gradients(vol):
    gradx = np.gradient(vol, edge_order=2, axis=0)
    grady = np.gradient(vol, edge_order=2, axis=1)
    gradz = np.gradient(vol, edge_order=2, axis=2)
    return np.stack([gradx, grady, gradz], axis=-1)

grad = gradients(warped)
mag = norm(grad, axis=-1)


# cluster the gradient vectors
# seven clusters
# one per face, and background/zero

gradvectors = grad.reshape((-1, 3))

from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=7, random_state=0, n_init=10).fit(gradvectors)
labels = kmeans.labels_
clusters = kmeans.cluster_centers_

print("clusters (mean):")
print(clusters)

# median is not advisable. median is a single value, not a mean.
# smarter to discard outliers, then take mean
def inliers_of(target, vecs):
    dists = norm(vecs - target, axis=1)
    keep = (dists <= np.quantile(dists, 0.5))
    return vecs[keep]

clusters = np.array([
    inliers_of(cluster, gradvectors[labels == i]).mean(axis=0)
    for i, cluster in enumerate(clusters)
])

print("clusters (mean of 90%):")
print(clusters)

(hist, edges) = np.histogram(labels, bins=np.arange(kmeans.n_clusters+1))
print("counts:", hist)

#bg_label = np.argmax(hist) # most common
bg_label = np.argmin(norm(clusters, axis=1)) # closest to zero
print("background label:", bg_label)

face_vecs = np.delete(clusters, bg_label, axis=0)
#face_vecs /= norm(face_vecs, axis=1, keepdims=True)
print("foreground clusters:")
print(face_vecs)

clusters (mean):
[[-0.00001  0.00004  0.     ]
 [-0.00007  0.00034  0.49484]
 [-0.00007  0.00034 -0.49484]
 [ 0.29297  0.43958  0.     ]
 [-0.29142 -0.4349   0.     ]
 [-0.4351   0.27852  0.     ]
 [ 0.4278  -0.28183  0.     ]]
clusters (mean of 90%):
[[ 0.      -0.       0.     ]
 [ 0.       0.       0.5    ]
 [ 0.       0.      -0.5    ]
 [ 0.29585  0.45781  0.     ]
 [-0.28876 -0.4558   0.     ]
 [-0.44932  0.28824  0.     ]
 [ 0.4462  -0.27927  0.     ]]
counts: [120012   1546   1546    622    632    316    326]
background label: 0
foreground clusters:
[[ 0.       0.       0.5    ]
 [ 0.       0.      -0.5    ]
 [ 0.29585  0.45781  0.     ]
 [-0.28876 -0.4558   0.     ]
 [-0.44932  0.28824  0.     ]
 [ 0.4462  -0.27927  0.     ]]

# it should be three pairs of face vectors, where each pair is opposite faces
# check which vectors are a pair by checking the dot product

def find_opposite_faces(face_vecs):
    pairs = []

    for i, vec1 in enumerate(face_vecs):
        for j, vec2 in enumerate(face_vecs):
            if i >= j: continue
            dist = np.dot(normalize(vec1), normalize(vec2))
            if dist < -0.999:
                pairs.append((i, j))
                print(f"({i}, {j}) dist {dist}")

    return pairs

pairs = find_opposite_faces(face_vecs)
print(pairs)

(0, 1) dist -1.0
(2, 3) dist -0.9999593496322632
(4, 5) dist -0.9999381303787231
[(0, 1), (2, 3), (4, 5)]

# pick one vector per pair which should be pointing positively (dot with [1,1,1] is positive)

axis_vecs = np.array([
    face_vecs[i1] if np.dot(face_vecs[i1], [1,1,1]) >= 0 else face_vecs[i2]
    for (i1, i2) in pairs
])

print(axis_vecs)

[[ 0.       0.       0.5    ]
 [ 0.29585  0.45781  0.     ]
 [ 0.4462  -0.27927  0.     ]]

# vectors are normal
# build basis from the three vectors
# for least amount of "accidental" rotation, find closest to X, Y, Z axes
# 1. find the one closest to X axis, remove
# 2. find the one closest to Y axis, remove
# 3. remaining one taken for Z axis

def find_closest(vecs, target):
    dists = np.dot(target, vecs.T)
    return np.argmax(dists)

remaining = [0, 1, 2]

index = find_closest(axis_vecs[remaining], [1,0,0])
x_axis = axis_vecs[remaining[index]]
remaining.pop(index)

index = find_closest(axis_vecs[remaining], [0,1,0])
y_axis = axis_vecs[remaining[index]]
remaining.pop(index)

z_axis = axis_vecs[remaining[0]]

basis = np.array([x_axis, y_axis, z_axis]).T
print(basis)

[[ 0.4462   0.29585  0.     ]
 [-0.27927  0.45781  0.     ]
 [ 0.       0.       0.5    ]]

# this basis may not be orthogonal, so we need to adjust it

# Modified Gram-Schmidt process, for numerical stability
# https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process

def modified_gram_schmidt(basis):
    """
    basis: NxN matrix of column vectors
    """
    basis = np.asarray(basis, dtype=np.float64)
    basis = basis.copy()
    M,N = basis.shape

    for i in range(N):
        basis[:,i] /= norm(basis[:,i])

        for j in range(i+1, N):
            basis[:,j] -= np.dot(basis[:,j], basis[:,i]) * basis[:,i]

    return basis

mgs_basis = modified_gram_schmidt(basis)
print(mgs_basis)

[[ 0.84766  0.53054  0.     ]
 [-0.53054  0.84766  0.     ]
 [ 0.       0.       1.     ]]

# how close is this to the 30 degrees applied originally?
(rvec, jac) = cv.Rodrigues(mgs_basis)
print(np.rad2deg(norm(rvec)))

32.04207290905596

# now rotate the volume to align with the basis

Tt = translate3(*(vsize/2))
Tr = rotate3(mgs_basis)

T_unwarp = Tt @ Tr @ inv(Tt)

unwarped = ndi.affine_transform(
    input=warped,
    matrix=T_unwarp,
    # offset ignored
    output_shape=vsize,
    order=1,
    mode='nearest',
)

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.