Machine perception of three-dimensional solids

解决问题：通过二维图像以及预存储的模型信息求出世界中的物体信息

$模型到图像变换：模型到物体变换：物体到图像变换：则有：目标：求出$
对于世界中的物体信息，假设在齐次坐标系中物体的某一点为，则物体到图像的变换矩阵：

其中：，描述相机在做中心投影时的缩放夹角，易得二维图像中的物体点为：

其中信息在压缩后无用，新的缩放比例因子调整为

人眼感知光线变化为非线性的，可用近似，感知图像边缘变化可以使用罗伯茨交叉算子（Roberts Cross operator）：刻画，通过交叉感知到所有方向的灰度变化，一种可能的复现代码：

# 罗伯茨交叉算子提取边缘

img_sqrt = np.sqrt(img_gray.astype(np.float32))

kx = np.array([[1, 0],
               [0, -1]], dtype=np.float32)
ky = np.array([[0, -1],
               [1, 0]], dtype=np.float32)

gx = cv2.filter2D(img_sqrt, cv2.CV_32F, kx)
gy = cv2.filter2D(img_sqrt, cv2.CV_32F, ky)

g = np.sqrt(gx**2 + gy**2)

g_norm = cv2.normalize(g, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)

cv2.imwrite('img_differentiated.png', g_norm)

提取特征点与方向：

程序首先在微分图像中，以 4x4 个像素为一个方块，寻找每个方块内的最大样本值。

如果该最大值超过了一个设定的较小阈值（作者建议将其设定为原始图像最大强度的十分之一左右），则将其视为一个候选点。

接着，程序以该候选点为中心，使用长度为 5 个样本、斜率分别为 0、1、0、-1 的四条线段，与周围的数据进行相关性拟合。

计算最佳拟合与最差拟合结果的比率。如果该比率大于第二个阈值（作者设定约为 3），程序就会记录下这个特征点及其算出的最佳方向。

# 特征点选取

T1 = np.max(g) * 0.1
T2 = 1.5

h, w = g_norm.shape

img_features = np.zeros((h, w), dtype=np.uint8)

mask_0 = np.zeros((5, 5), dtype=np.float32)
mask_0[2, :] = 1

mask_90 = np.zeros((5, 5), dtype=np.float32)
mask_90[:, 2] = 1

mask_45 = np.fliplr(np.eye(5, dtype=np.float32))
mask_135 = np.eye(5, dtype=np.float32)

masks = [mask_0, mask_45, mask_90, mask_135]

directions = [(2, 0), (2, -2), (0, 2), (2, 2)]
feature_points = []

for y in range(2, h - 6, 4):
    for x in range(2, w - 6, 4):
        block = g[y:y+4, x:x+4]
        min_val, max_val, min_loc, max_loc = cv2.minMaxLoc(block)

        if max_val > T1:
            gy_loc = y + max_loc[1]
            gx_loc = x + max_loc[0]

            roi = g_norm[gy_loc-2:gy_loc+3, gx_loc-2:gx_loc+3].astype(np.float32)

            correlations = []
            for mask in masks:
                corr = np.sum(roi * mask)
                correlations.append(corr)

            best_fit = max(correlations)
            worst_fit = min(correlations)

            ratio = best_fit / (worst_fit + 1e-5)

            if ratio > T2:
                best_idx = correlations.index(best_fit)
                feature_points.append((gx_loc, gy_loc, best_idx))

                dx, dy = directions[best_idx]
                pt1 = (gx_loc - dx, gy_loc - dy)
                pt2 = (gx_loc + dx, gy_loc + dy)
                cv2.line(img_features, pt1, pt2, 255, 1)

cv2.imwrite('img_features.png', img_features)

连接特征点：

程序会尝试在相邻的点之间建立连线。

两个点相连必须同时满足两个条件：第一，它们必须位于相互接触的 4x4 方块中；第二，这两点连线的实际方向，与这两个点各自在上一阶段记录的最佳方向相比，误差必须在范围内。

在完成所有符合条件的连接后，任何未能与其他点连接的孤立点都会被直接删除，从而过滤掉大部分由背景噪声产生的点。

# 连接特征点

T3 = 23

img_connected = np.zeros((h, w), dtype=np.uint8)

idx_to_angle = {0: 0, 1: 135, 2: 90, 3: 45}

def angle_diff(a1, a2):
    diff = abs(a1 - a2) % 180
    return min(diff, 180 - diff)

valid_lines = []
connected_indices = set()
n_pts = len(feature_points)

for i in range(n_pts):
    x1, y1, dir_idx1 = feature_points[i]
    ang1 = idx_to_angle[dir_idx1]

    bx1 = (x1 - 2) // 4
    by1 = (y1 - 2) // 4

    for j in range(i + 1, n_pts):
        x2, y2, dir_idx2 = feature_points[j]
        bx2 = (x2 - 2) // 4
        by2 = (y2 - 2) // 4

        if max(abs(bx1 - bx2), abs(by1 - by2)) <= 1:
            dx = x2 - x1
            dy = y2 - y1

            line_ang = math.degrees(math.atan2(dy, dx)) % 180
            ang2 = idx_to_angle[dir_idx2]

            if angle_diff(line_ang, ang1) <= T3 or angle_diff(line_ang, ang2) <= T3:
                valid_lines.append((x1, y1, x2, y2))
                connected_indices.add(i)
                connected_indices.add(j)

for x1, y1, x2, y2 in valid_lines:
    cv2.line(img_connected, (x1, y1), (x2, y2), 255, 1)

cv2.imwrite('img_connected.png', img_connected)

拓扑清理与降维：

清理三角形：寻找由三条线组成的闭合环（三角形），直接删除其中最长的一条边。

压缩四边形：寻找由四条线组成的闭合环，沿着其最短的对角线进行压缩（即直接将距离最近的两个非相邻顶点合并）。

去除短尾巴：寻找那些一端悬空（没有连接任何东西）、另一端却连接了多条线的小短线（spurs），并将它们删除。

# 降低复杂度

adj = {}
for x1, y1, x2, y2 in valid_lines:
    p1, p2 = (x1, y1), (x2, y2)
    if p1 not in adj:
        adj[p1] = set()
    if p2 not in adj:
        adj[p2] = set()
    adj[p1].add(p2)
    adj[p2].add(p1)

def dist(p1, p2):
    return math.hypot(p1[0] - p2[0], p1[1] - p2[1])

triangles = set()
for u in adj:
    for v in adj[u]:
        if u < v:
            common_neighbors = adj[u].intersection(adj[v])
            for w in common_neighbors:
                if v < w:
                    triangles.add((u, v, w))

edges_to_remove = set()
for u, v, w in triangles:
    d1, d2, d3 = dist(u, v), dist(v, w), dist(w, u)
    max_d = max(d1, d2, d3)
    if max_d == d1:
        edges_to_remove.add((u, v))
    elif max_d == d2:
        edges_to_remove.add((v, w))
    else:
        edges_to_remove.add((w, u))

for u, v in edges_to_remove:
    if v in adj.get(u, set()):
        adj[u].remove(v)
    if u in adj.get(v, set()):
        adj[v].remove(u)

quads = {}
for u in list(adj.keys()):
    if u not in adj: continue
    neighbors = list(adj[u])
    for i in range(len(neighbors)):
        for j in range(i + 1, len(neighbors)):
            v1, v2 = neighbors[i], neighbors[j]
            common = adj[v1].intersection(adj[v2]) - {u}
            for w in common:
                if w not in adj[u]:
                    quad_nodes = tuple(sorted([u, v1, w, v2]))
                    if quad_nodes not in quads:
                        quads[quad_nodes] = (u, v1, w, v2)

merges = []
for u, v1, w, v2 in quads.values():
    if dist(u, w) < dist(v1, v2):
        merges.append((u, w))
    else:
        merges.append((v1, v2))

for n1, n2 in merges:
    if n1 in adj and n2 in adj:
        new_nodes = ((n1[0] + n2[0]) // 2, (n1[1] + n2[1]) // 2)
        new_neighbors = adj[n1].union(adj[n2]) - {n1, n2}
        del adj[n1]
        del adj[n2]
        if new_nodes not in adj:
            adj[new_nodes] = set()

        for neighbor in new_neighbors:
            if neighbor in adj:
                adj[neighbor].discard(n1)
                adj[neighbor].discard(n2)
                adj[neighbor].add(new_nodes)
                adj[new_nodes].add(neighbor)

spurs_to_remove = set()
for u in list(adj.keys()):
    if u in adj and len(adj[u]) == 1:
        v = list(adj[u])[0]
        if len(adj[v]) >= 3:
            spurs_to_remove.add((u, v))

for u, v in spurs_to_remove:
    if v in adj.get(u, set()):
        adj[u].remove(v)
        adj[v].remove(u)

h, w = g_norm.shape
img_reduced = np.zeros((h, w), dtype=np.uint8)
drawn_edges = set()
for u in adj:
    for v in adj[u]:
        edge = tuple(sorted([u, v]))
        if edge not in drawn_edges:
            cv2.line(img_reduced, u, v, 255, 1)
            drawn_edges.add(edge)

cv2.imwrite('img_reduced.png', img_reduced)

长直线拟合：

程序任选一小段线段作为起点，沿着一个方向不断将相邻的点纳入直线方程计算，直到遇到一个连接了其他两条以上线段的交叉点，或者直线的均方误差超过了设定的阈值。

当误差超过阈值时（通常意味着越过了真实的拐角点），程序会启动回退（backup）机制，逐个剔除刚才加入的点，直到原小线段与计算出的直线之间的夹角减小到原来的一半为止。

程序会在直线的另一个方向上重复上述前进和回退操作。

当所有长直线生成后，原来短线段的端点会被移动到计算出的长直线上。如果有几条长直线交汇于一点，该点的坐标会被重新计算为最长的两条非平行直线的交点。

# 初始直线拟合
chains = []
visited_edges = set()
keypoints = [node for node, neighbors in adj.items() if len(neighbors) != 2]

for u in keypoints:
    for v in adj[u]:
        edge = tuple(sorted([u, v]))
        if edge not in visited_edges:
            chain = [u, v]
            visited_edges.add(edge)

            prev, curr = u, v
            while len(adj[curr]) == 2:
                next_node = list(adj[curr] - {prev})[0]
                chain.append(next_node)
                visited_edges.add(tuple(sorted([curr, next_node])))
                prev, curr = curr, next_node

            chains.append(chain)

for u in list(adj.keys()):
    if len(adj[u]) == 2:
        for v in adj[u]:
            edge = tuple(sorted([u, v]))
            if edge not in visited_edges:
                chain = [u, v]
                visited_edges.add(edge)
                prev, curr = u, v
                while curr != u:
                    next_node = list(adj[curr] - {prev})[0]
                    chain.append(next_node)
                    visited_edges.add(tuple(sorted([curr, next_node])))
                    prev, curr = curr, next_node
                chains.append(chain)

h_img, w_img = g_norm.shape
img_fitted = np.zeros((h_img, w_img), dtype=np.uint8)

T4_error = 2.0

for chain in chains:
    pts = np.array(chain, dtype=np.int32)
    approx = cv2.approxPolyDP(pts, epsilon=T4_error, closed=False)
    approx = approx.reshape(-1, 2)

    for i in range(len(approx) - 1):
        pt1 = tuple(approx[i])
        pt2 = tuple(approx[i + 1])
        cv2.line(img_fitted, pt1, pt2, 255, 1)

cv2.imwrite('img_fitted.png', img_fitted)

最终补全与合并：

线条延伸补全：如果一条直线的长度足够长，程序会检查其两端沿延长线方向约 3 个光栅（像素）单位内的最近点。程序会在端点和这个新点之间与原微分图像进行相关性测试，如果平均值大于设定阈值，就会补上这段缺失的线。

超短线处理：对于非常短的线段，如果它的两个端点各自都只连接了一条其他线，程序就会将这两个端点合并；否则就直接删除这条短线。

最后，将整个线条结构再次送入均方误差直线拟合程序中处理一次，以消除在上述补全和合并过程中可能产生的新多余关节。

# 最终线条图
adj_H = {}
for chain in chains:
    pts = np.array(chain, dtype=np.int32)
    approx = cv2.approxPolyDP(pts, epsilon=T4_error, closed=False)
    approx = approx.reshape(-1, 2)

    for i in range(len(approx) - 1):
        u = (int(approx[i][0]), int(approx[i][1]))
        v = (int(approx[i + 1][0]), int(approx[i + 1][1]))
        if u == v: continue
        if u not in adj_H:
            adj_H[u] = set()
        if v not in adj_H:
            adj_H[v] = set()
        adj_H[u].add(v)
        adj_H[v].add(u)

merge_threshold = 4.0
nodes_H = list(adj_H.keys())
for i in range(len(nodes_H)):
    for j in range(i + 1, len(nodes_H)):
        n1, n2 = nodes_H[i], nodes_H[j]
        if n1 in adj_H and n2 in adj_H and n1 != n2:
            if dist(n1, n2) < merge_threshold:
                neighbors = list(adj_H[n2])
                del adj_H[n2]
                for neighbor in neighbors:
                    if neighbor != n1 and neighbor in adj_H:
                        adj_H[neighbor].discard(n2)
                        adj_H[neighbor].add(n1)
                        adj_H[n1].add(neighbor)

short_threshold = 4.0
spurs_to_delete = []
for u in list(adj_H.keys()):
    if u in adj_H and len(adj_H[u]) == 1:
        v = list(adj_H[u])[0]
        if dist(u, v) < short_threshold:
            spurs_to_delete.append((u, v))

for u, v in spurs_to_delete:
    if u in adj_H:
        adj_H[u].discard(v)
        if len(adj_H[u]) == 0: del adj_H[u]
    if v in adj_H and u in adj_H[v]:
        adj_H[v].discard(u)
        if len(adj_H[v]) == 0: del adj_H[v]

def angle_between(p1, p2, p3):
    v1 = (p1[0] - p2[0], p1[1] - p2[1])
    v2 = (p3[0] - p2[0], p3[1] - p2[1])
    dot = v1[0]*v2[0] + v1[1]*v2[1]
    mag1 = math.hypot(*v1)
    mag2 = math.hypot(*v2)
    if mag1 == 0 or mag2 == 0: return 0
    val = dot / (mag1 * mag2)
    val = max(-1.0, min(1.0, val))
    return math.degrees(math.acos(val))

parallel_threshold = 165

for curr in list(adj_H.keys()):
    if curr in adj_H and len(adj_H[curr]) == 2:
        n1, n2 = list(adj_H[curr])
        if n1 not in adj_H or n2 not in adj_H:
            continue
        ang = angle_between(n1, curr, n2)
        if ang > parallel_threshold:
            adj_H[n1].discard(curr)
            adj_H[n2].discard(curr)
            adj_H[n1].add(n2)
            adj_H[n2].add(n1)
            del adj_H[curr]

img_final = np.zeros((h_img, w_img), dtype=np.uint8)
drawn_H = set()
for u in adj_H:
    for v in adj_H[u]:
        edge = tuple(sorted([u, v]))
        if edge not in drawn_H:
            cv2.line(img_final, u, v, 255, 1)
            drawn_H.add(edge)

cv2.imwrite('img_final.png', img_final)

ConnectC

Machine perception of three-dimensional solids

Machine perception of three-dimensional solids

2026 BUAA 卓工保研经验

ConnectC

Comments | NOTHING