iomgaa
00ebbffb28
- 第 1 步分屏后在两半扇形区内用边缘图归一化互相关估 (dx, dy)，再用逐帧多块模板匹配的中位数独立核对； 两者各分量之差超过门槛（15 像素）这一例第 1 步失败。删掉「两个视图左上角相减」的实现与测试。 - 门槛依据：34 例中走完第 1 步的 29 例，差值中位 2、95 分位 6.2、最大 13；参数与门槛在同一批视频上调出。 - scripts/ceus_align_review.py：只做第 1 步，出扇形区棋盘图与边缘叠加图，混入错开 15–30 像素的 阴性对照、随机编号；计分时任一对照未被判「没对齐」即作废。 - ceus_quant 读 review.csv：不通过的病例第 1 步失败；表里没有这段视频或平移对不上即停。 - 配准与紫圈颜色参数进分割缓存键。 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
278 lines
13 KiB
Python
278 lines
13 KiB
Python
"""第 1 步的两半配准：估出左半局部坐标到右半局部坐标的平移 `(dx, dy)`。
|
||
|
||
两半是同一个探头视野的两路成像（造影和灰阶），大小相同，但扇形在各自那一半里的位置不同：
|
||
造影那一半的扇形左边常有一段黑边，灰阶那一半的扇形贴着分屏线，而且同一帧宽的两段视频错位也
|
||
可能不同。所以每段视频单独配准，不写死、也不用两个视图左上角相减（任务书 02 第一阶段第 1 条）。
|
||
|
||
做法：均匀抽若干帧，两半各自求逐像素的时间平均灰度和时间标准差。
|
||
|
||
1. 扇形区：时间标准差高于阈值的像素（画面在动的地方；文字、色标、刻度这些静止的界面元素
|
||
标准差近 0），去掉顶部信息栏那几行，开运算断开细连接，取最大连通域的凸包。
|
||
2. 绿框、紫圈的颜色（`prompt_source`、`pan_source` 的 HSV 范围）在任一抽样帧出现过的像素
|
||
膨胀后用周围像素修补掉：它们只画在一半上，会被当成边缘或纹理去匹配。
|
||
3. 方法一（互相关，结果取它）：两半的边缘图（高斯平滑后的梯度幅值，只留扇形区和它的边）做
|
||
归一化互相关，在搜索范围内穷举平移取最大。扇形的轮廓和组织边界都参与。
|
||
4. 方法二（模板匹配，只用来核对）：在左半扇形内按网格定若干块，每一抽样帧上各切一次平滑后的
|
||
灰度，在右半搜索范围内、块中心落在右半扇形内的位置里找最相似的（归一化相关系数），全部
|
||
匹配的平移各分量取中位数。它用的是逐帧灰度纹理、局部块和中位数，不用边缘和扇形轮廓，和
|
||
方法一相互独立。
|
||
5. 两种估计各分量之差的最大值（切比雪夫距离：水平、竖直各自都不许错开超过门槛）超过门槛，
|
||
这一例第 1 步失败（`RegistrationError`）。
|
||
|
||
只做平移，不做缩放和旋转。
|
||
"""
|
||
|
||
from collections.abc import Sequence
|
||
from dataclasses import dataclass
|
||
from typing import Literal
|
||
|
||
import cv2
|
||
import numpy as np
|
||
import numpy.typing as npt
|
||
|
||
from cadx.ceus._config import PanSourceConfig, PromptSourceConfig, RegistrationConfig
|
||
from cadx.ceus._errors import CeusCaseError
|
||
from cadx.ceus._geometry import Frame, Mask, ViewRect, crop_frame
|
||
from cadx.ceus._mapping import Offset
|
||
from cadx.ceus._split_screen import sample_frames_for_split
|
||
|
||
REGISTRATION_STEP = 1
|
||
|
||
type _Gray = npt.NDArray[np.float32]
|
||
|
||
|
||
class RegistrationError(CeusCaseError):
|
||
"""两半配准失败：找不到扇形区、可用的模板块不够，或两种方法估出的平移不一致。"""
|
||
|
||
def __init__(self, reason: str) -> None:
|
||
super().__init__(REGISTRATION_STEP, reason)
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
class TemplatePatch:
|
||
"""方法二的一次匹配：模板块在左半局部坐标的左上角、第几张抽样帧、在右半找到的平移、相关系数。"""
|
||
|
||
x: int
|
||
y: int
|
||
frame: int
|
||
offset: Offset
|
||
score: float
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
class HalfRegistration:
|
||
"""左半局部坐标 → 右半局部坐标的平移：`右半局部 = 左半局部 + (dx, dy)`。
|
||
|
||
`correlation` 是方法一（边缘图互相关）的估计，采用它；`correlation_score` 是它的峰值；
|
||
`template` 是方法二（多块模板逐帧匹配取中位数）的估计，`patches` 是每次匹配的明细；
|
||
`disagreement_px` 是两者各分量之差的最大值。
|
||
"""
|
||
|
||
correlation: Offset
|
||
correlation_score: float
|
||
template: Offset
|
||
patches: tuple[TemplatePatch, ...]
|
||
disagreement_px: int
|
||
|
||
def contrast_to_target(self, contrast_side: Literal["left", "right"]) -> Offset:
|
||
"""造影那一半局部坐标 → 另一半局部坐标的平移。"""
|
||
dx, dy = self.correlation
|
||
if contrast_side == "left":
|
||
return (dx, dy)
|
||
return (-dx, -dy)
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
class _Half:
|
||
"""一半画面的配准素材：扇形区掩膜，修补绿框紫圈并平滑后的时间平均灰度和逐帧灰度。"""
|
||
|
||
fan: Mask
|
||
mean: _Gray
|
||
frames: list[_Gray]
|
||
|
||
|
||
def color_marks(
|
||
frames: Sequence[Frame], prompt_cfg: PromptSourceConfig, pan_cfg: PanSourceConfig, dilate: int
|
||
) -> Mask:
|
||
"""任一帧里落在绿框或紫圈 HSV 范围内的像素，膨胀 `dilate` 像素。"""
|
||
marks = np.zeros(frames[0].shape[:2], dtype=np.uint8)
|
||
for frame in frames:
|
||
hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
|
||
for cfg in (prompt_cfg, pan_cfg):
|
||
lower = np.asarray(cfg.hsv_lower, dtype=np.uint8)
|
||
upper = np.asarray(cfg.hsv_upper, dtype=np.uint8)
|
||
marks |= np.asarray(cv2.inRange(hsv, lower, upper), dtype=np.uint8)
|
||
if dilate > 0:
|
||
kernel = np.ones((2 * dilate + 1, 2 * dilate + 1), dtype=np.uint8)
|
||
marks = np.asarray(cv2.dilate(marks, kernel), dtype=np.uint8)
|
||
return marks > 0
|
||
|
||
|
||
def fan_mask(std: _Gray, cfg: RegistrationConfig) -> Mask:
|
||
"""扇形区：时间标准差高于阈值、去掉顶部信息栏、开运算后最大连通域的凸包。"""
|
||
moving = (std > cfg.fan_std_threshold).astype(np.uint8)
|
||
moving[: round(cfg.top_ignore_ratio * std.shape[0])] = 0
|
||
size = cfg.fan_open_kernel
|
||
moving = cv2.morphologyEx(moving, cv2.MORPH_OPEN, np.ones((size, size), dtype=np.uint8))
|
||
count, labels, stats, _ = cv2.connectedComponentsWithStats(moving, connectivity=8)
|
||
if count < 2:
|
||
raise RegistrationError("这一半画面里找不到扇形区（没有随时间变化的区域）")
|
||
largest = 1 + int(np.argmax(stats[1:, cv2.CC_STAT_AREA]))
|
||
component = np.asarray(labels, dtype=np.int32) == largest
|
||
points = np.argwhere(component)[:, ::-1].astype(np.int32) # (x, y)
|
||
hull = np.asarray(cv2.convexHull(points), dtype=np.int32)
|
||
filled = np.zeros(std.shape, dtype=np.uint8)
|
||
cv2.fillConvexPoly(filled, hull, 1)
|
||
return filled > 0
|
||
|
||
|
||
def _repair_smooth(gray: _Gray, marks: Mask, sigma: float) -> _Gray:
|
||
"""绿框紫圈像素用周围像素修补（Telea 法，半径 3），再高斯平滑。"""
|
||
repaired = cv2.inpaint(
|
||
np.clip(gray, 0, 255).astype(np.uint8), marks.astype(np.uint8), 3, cv2.INPAINT_TELEA
|
||
)
|
||
return np.asarray(cv2.GaussianBlur(repaired.astype(np.float32), (0, 0), sigma), np.float32)
|
||
|
||
|
||
def _prepare_half(gray: _Gray, view: ViewRect, marks: Mask, cfg: RegistrationConfig) -> _Half:
|
||
"""从整帧灰度栈 `(N, H, W)` 切出一半，求扇形区、平均灰度和逐帧灰度。"""
|
||
stack = gray[:, view.y1 : view.y2, view.x1 : view.x2]
|
||
local_marks = crop_frame(marks, view)
|
||
return _Half(
|
||
fan=fan_mask(stack.std(axis=0), cfg),
|
||
mean=_repair_smooth(stack.mean(axis=0), local_marks, cfg.blur_sigma),
|
||
frames=[_repair_smooth(frame, local_marks, cfg.blur_sigma) for frame in stack],
|
||
)
|
||
|
||
|
||
def _edge_map(half: _Half, sigma: float) -> _Gray:
|
||
"""时间平均灰度的梯度幅值，只留扇形区和它的边（扇形掩膜外扩 3σ，轮廓上的梯度也算进来）。"""
|
||
gx = cv2.Sobel(half.mean, cv2.CV_32F, 1, 0)
|
||
gy = cv2.Sobel(half.mean, cv2.CV_32F, 0, 1)
|
||
margin = int(np.ceil(3 * sigma))
|
||
kernel = np.ones((2 * margin + 1, 2 * margin + 1), dtype=np.uint8)
|
||
support = cv2.dilate(half.fan.astype(np.uint8), kernel) > 0
|
||
edges = np.hypot(gx, gy).astype(np.float32)
|
||
edges[~support] = 0
|
||
return edges
|
||
|
||
|
||
def _bbox(mask: Mask) -> tuple[int, int, int, int]:
|
||
"""掩膜的外接框 `(x1, y1, x2, y2)`，右下不含。"""
|
||
ys, xs = np.nonzero(mask)
|
||
return int(xs.min()), int(ys.min()), int(xs.max()) + 1, int(ys.max()) + 1
|
||
|
||
|
||
def _correlate_edges(
|
||
left: _Half, right: _Half, search: Offset, sigma: float
|
||
) -> tuple[Offset, float]:
|
||
"""方法一：左半边缘图（裁到扇形外接框）在补零后的右半边缘图上做归一化互相关。"""
|
||
sx, sy = search
|
||
left_edges = _edge_map(left, sigma)
|
||
right_edges = _edge_map(right, sigma)
|
||
x1, y1, x2, y2 = _bbox(left_edges > 0)
|
||
template = left_edges[y1:y2, x1:x2]
|
||
height, width = left_edges.shape
|
||
# 右半四周补零，使平移 dx ∈ [−sx, sx]、dy ∈ [−sy, sy] 都落在结果图里；右、下多补一个左半大小，
|
||
# 右半比左半窄几个像素时也够。
|
||
padded = cv2.copyMakeBorder(
|
||
right_edges, sy, sy + height, sx, sx + width, cv2.BORDER_CONSTANT, value=0.0
|
||
)
|
||
scores = cv2.matchTemplate(padded, template, cv2.TM_CCORR_NORMED)
|
||
# 结果图里位置 (u, v) 对应平移 dx = u − sx − x1、dy = v − sy − y1。
|
||
window = scores[y1 : y1 + 2 * sy + 1, x1 : x1 + 2 * sx + 1]
|
||
_, best, _, (u, v) = cv2.minMaxLoc(window)
|
||
return (u - sx, v - sy), float(best)
|
||
|
||
|
||
def _match_templates(
|
||
left: _Half, right: _Half, search: Offset, cfg: RegistrationConfig
|
||
) -> tuple[TemplatePatch, ...]:
|
||
"""方法二：左半扇形外接框上按网格切块，整块都在扇形内的才用；每块在每一抽样帧上各找一次。
|
||
|
||
逐帧匹配而不是在时间平均图上匹配：探头在视频里移动时，平均图上的组织纹理糊掉了。
|
||
"""
|
||
sx, sy = search
|
||
rows, cols = cfg.template_grid
|
||
side = round(cfg.template_size_ratio * left.fan.shape[1])
|
||
x1, y1, x2, y2 = _bbox(left.fan)
|
||
right_h, right_w = right.fan.shape
|
||
patches: list[TemplatePatch] = []
|
||
for row in range(rows):
|
||
for col in range(cols):
|
||
px = round(x1 + (col + 0.5) * (x2 - x1) / cols - side / 2)
|
||
py = round(y1 + (row + 0.5) * (y2 - y1) / rows - side / 2)
|
||
inside_view = px >= 0 and py >= 0 and px + side <= left.fan.shape[1]
|
||
inside_view = inside_view and py + side <= left.fan.shape[0]
|
||
if not inside_view or not left.fan[py : py + side, px : px + side].all():
|
||
continue
|
||
wx1, wy1 = max(0, px - sx), max(0, py - sy)
|
||
wx2, wy2 = min(right_w, px + side + sx), min(right_h, py + side + sy)
|
||
if wx2 - wx1 < side or wy2 - wy1 < side:
|
||
continue
|
||
# 只认块中心落在右半扇形内的位置：扇形外是纯黑，零方差窗口的相关系数 OpenCV 记为 1。
|
||
half = side // 2
|
||
centers = right.fan[
|
||
wy1 + half : wy2 - side + half + 1, wx1 + half : wx2 - side + half + 1
|
||
]
|
||
if not centers.any():
|
||
continue
|
||
for index, (left_frame, right_frame) in enumerate(
|
||
zip(left.frames, right.frames, strict=True)
|
||
):
|
||
template = left_frame[py : py + side, px : px + side]
|
||
scores = cv2.matchTemplate(
|
||
right_frame[wy1:wy2, wx1:wx2], template, cv2.TM_CCOEFF_NORMED
|
||
)
|
||
scores[~centers] = -1.0
|
||
_, best, _, (u, v) = cv2.minMaxLoc(scores)
|
||
offset = (wx1 + u - px, wy1 + v - py)
|
||
patches.append(TemplatePatch(px, py, index, offset, float(best)))
|
||
return tuple(patches)
|
||
|
||
|
||
def register_halves(
|
||
frames: Sequence[Frame],
|
||
left_view: ViewRect,
|
||
right_view: ViewRect,
|
||
cfg: RegistrationConfig,
|
||
prompt_cfg: PromptSourceConfig,
|
||
pan_cfg: PanSourceConfig,
|
||
) -> HalfRegistration:
|
||
"""两种方法各估一次左半 → 右半的平移，相差超过门槛抛 `RegistrationError`。"""
|
||
sampled = sample_frames_for_split(frames, cfg.sample_frames)
|
||
gray = np.stack([cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY) for frame in sampled])
|
||
gray = gray.astype(np.float32)
|
||
marks = color_marks(sampled, prompt_cfg, pan_cfg, cfg.mark_dilate)
|
||
left = _prepare_half(gray, left_view, marks, cfg)
|
||
right = _prepare_half(gray, right_view, marks, cfg)
|
||
search = (
|
||
round(cfg.search_x_ratio * left_view.width),
|
||
round(cfg.search_y_ratio * left_view.height),
|
||
)
|
||
|
||
correlation, correlation_score = _correlate_edges(left, right, search, cfg.blur_sigma)
|
||
patches = _match_templates(left, right, search, cfg)
|
||
blocks = len({(patch.x, patch.y) for patch in patches})
|
||
if blocks < cfg.template_min_patches:
|
||
raise RegistrationError(
|
||
f"左半扇形内可用的模板块只有 {blocks} 块，少于 {cfg.template_min_patches} 块"
|
||
)
|
||
template = (
|
||
round(float(np.median([patch.offset[0] for patch in patches]))),
|
||
round(float(np.median([patch.offset[1] for patch in patches]))),
|
||
)
|
||
disagreement = max(abs(correlation[0] - template[0]), abs(correlation[1] - template[1]))
|
||
if disagreement > cfg.max_disagreement_px:
|
||
raise RegistrationError(
|
||
f"两半配准不一致：互相关估计 (dx, dy) = {correlation}，模板匹配估计 {template}，"
|
||
f"相差 {disagreement} 像素，超过门槛 {cfg.max_disagreement_px}"
|
||
)
|
||
return HalfRegistration(
|
||
correlation=correlation,
|
||
correlation_score=correlation_score,
|
||
template=template,
|
||
patches=patches,
|
||
disagreement_px=disagreement,
|
||
)
|