feat(ceus,segment): 任务书 02 第二阶段——分割和后半段算法
- ceus 迁入 polygons、dark_region、reference_region、tic_fitting、tic、anchor_frames、 perfusion_metrics、trajectory、tic_curves、tic_plot、overlay，旧测试照搬（改 import 与配置读取）。 - 对外两段接口：prepare_segmentation（第 1–2 步）与 quantify_case（第 4–9 步）。 - 第 3 条：删掉 D1 取医生标注时刻的模式，D1 只从拟合曲线算。 - 第 4 条：F1/F2 取视频第 20−X、60−X 秒（X 为视频第 0 秒时已推注的秒数），落在视频外记未定位。 - 绿框融合权重、紫圈颜色阈值照搬；关键帧容差窗口的出处写进 ceus.yaml 注释。 - segment：同一进程调用 MedSAM2，设置写进 ceus.yaml 的 segment 段，全部必填；gpu 测试真跑通过。 - 依赖加 scipy、matplotlib、medsam2（可编辑安装，不编译 CUDA 扩展）；uv 只解析 linux x86_64。 - 导入契约：只有 segment 引用 torch、sam2；segment 不引用其他 cadx 模块。 Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
+142
-3
@@ -1,10 +1,23 @@
|
||||
"""ceus：CEUS 定量的全部算法和参数，纯计算，不引用其他 cadx 模块。
|
||||
|
||||
目前有第 1、2 步（视频解码、分屏、找绿框、选提示帧、判断造影画面在哪一半）和两半之间的坐标换算；
|
||||
第 3 步以后随任务书 02 第二阶段加入。
|
||||
某一例不满足算法前提时抛 `CeusCaseError`（带失败的步号和原因）。
|
||||
对外两段接口：`prepare_segmentation` 做第 1–2 步（解码、分屏、定造影那一半、选提示帧），交出分割的
|
||||
输入；`quantify_case` 拿分割结果做第 4–9 步（轮廓、暗区、背景圆、TIC、拟合、关键帧、灌注指标）。
|
||||
第 3 步分割由 segment 做，ceus 不引用它。某一例不满足算法前提时抛 `CeusCaseError`
|
||||
（带失败的步号和原因）。
|
||||
写产物用：`pack_trajectory`、`pack_tic_curves`、`render_tic_png` + `tic_plot_markers`、
|
||||
`render_case_overlay`、`keyframe_images`、`quant_summary`。
|
||||
"""
|
||||
|
||||
from cadx.ceus._anchor_frames import (
|
||||
ANCHOR_SLOTS,
|
||||
STATUS_ANCHOR_UNDEFINED,
|
||||
STATUS_NO_FRAME_IN_WINDOW,
|
||||
STATUS_NO_ORDERED_FRAME,
|
||||
STATUS_NOT_COVERED,
|
||||
STATUS_OK,
|
||||
AnchorFrameLocation,
|
||||
locate_anchor_frames,
|
||||
)
|
||||
from cadx.ceus._config import (
|
||||
AnchorArrivalConfig,
|
||||
AnchorFramesConfig,
|
||||
@@ -25,6 +38,15 @@ from cadx.ceus._config import (
|
||||
ViewsConfig,
|
||||
)
|
||||
from cadx.ceus._contrast_side import ContrastSide, ContrastSideError, locate_contrast_side
|
||||
from cadx.ceus._dark_region import (
|
||||
apply_dark_morphology,
|
||||
dark_candidate_mask,
|
||||
extract_dark_region_masks,
|
||||
extract_dark_region_polygons,
|
||||
fill_binary_holes,
|
||||
score_dark_candidate,
|
||||
split_dark_components,
|
||||
)
|
||||
from cadx.ceus._errors import CeusCaseError
|
||||
from cadx.ceus._geometry import (
|
||||
ViewRect,
|
||||
@@ -49,7 +71,33 @@ from cadx.ceus._green_prompt import (
|
||||
prompt_color_mask,
|
||||
)
|
||||
from cadx.ceus._mapping import map_source_mask_to_target, target_polygons_to_source, view_offset
|
||||
from cadx.ceus._overlay import OverlayConfig, OverlayLayers, draw_polygon_overlay, render_overlay
|
||||
from cadx.ceus._perfusion_metrics import (
|
||||
LognormalPerfusionMetrics,
|
||||
compute_lognormal_perfusion_metrics,
|
||||
)
|
||||
from cadx.ceus._pipeline import (
|
||||
LESION_CURVES,
|
||||
QuantResult,
|
||||
SegmentationInput,
|
||||
keyframe_images,
|
||||
prepare_segmentation,
|
||||
quant_summary,
|
||||
quantify_case,
|
||||
render_case_overlay,
|
||||
tic_plot_markers,
|
||||
)
|
||||
from cadx.ceus._polygons import (
|
||||
NoValidContourError,
|
||||
align_polygon_to_previous,
|
||||
apply_prompt_frame_position_constraints,
|
||||
fuse_prompt_prior_with_medsam,
|
||||
largest_contour,
|
||||
mask_center,
|
||||
masks_to_smoothed_polygons,
|
||||
)
|
||||
from cadx.ceus._prompt_frames import PromptFrame, select_prompt_frames
|
||||
from cadx.ceus._reference_region import select_reference_region
|
||||
from cadx.ceus._split_screen import (
|
||||
SplitDetectionError,
|
||||
detect_split_x,
|
||||
@@ -57,12 +105,54 @@ from cadx.ceus._split_screen import (
|
||||
split_views,
|
||||
view_rect,
|
||||
)
|
||||
from cadx.ceus._tic import (
|
||||
CURVE_NAMES,
|
||||
BaselineResult,
|
||||
TicResult,
|
||||
compute_tic_for_video,
|
||||
detect_baseline,
|
||||
mean_or_nan,
|
||||
)
|
||||
from cadx.ceus._tic_curves import (
|
||||
CURVES_SCHEMA_VERSION,
|
||||
TicCurves,
|
||||
TicCurveSeries,
|
||||
from_tic_result,
|
||||
pack_tic_curves,
|
||||
unpack_tic_curves,
|
||||
)
|
||||
from cadx.ceus._tic_fitting import (
|
||||
LognormalFitResult,
|
||||
LognormalParameters,
|
||||
evaluate_lognormal,
|
||||
evaluate_lognormal_derivative,
|
||||
fit_lognormal_curve,
|
||||
)
|
||||
from cadx.ceus._tic_plot import TicPlotMarker, render_tic_png
|
||||
from cadx.ceus._trajectory import (
|
||||
POLYGON_NUM_POINTS,
|
||||
Trajectory,
|
||||
pack_trajectory,
|
||||
unpack_trajectory,
|
||||
)
|
||||
from cadx.ceus._video import VideoFile
|
||||
|
||||
__all__ = [
|
||||
"ANCHOR_SLOTS",
|
||||
"CURVES_SCHEMA_VERSION",
|
||||
"CURVE_NAMES",
|
||||
"LESION_CURVES",
|
||||
"POLYGON_NUM_POINTS",
|
||||
"STATUS_ANCHOR_UNDEFINED",
|
||||
"STATUS_NOT_COVERED",
|
||||
"STATUS_NO_FRAME_IN_WINDOW",
|
||||
"STATUS_NO_ORDERED_FRAME",
|
||||
"STATUS_OK",
|
||||
"AnchorArrivalConfig",
|
||||
"AnchorFrameLocation",
|
||||
"AnchorFramesConfig",
|
||||
"AnchorLatePhaseConfig",
|
||||
"BaselineResult",
|
||||
"CeusCaseError",
|
||||
"CeusConfig",
|
||||
"ColorPrompt",
|
||||
@@ -72,42 +162,91 @@ __all__ = [
|
||||
"InteriorDarkRegionConfig",
|
||||
"KeyframePositionConstraint",
|
||||
"KeyframesConfig",
|
||||
"LognormalFitResult",
|
||||
"LognormalParameters",
|
||||
"LognormalPerfusionMetrics",
|
||||
"NoValidContourError",
|
||||
"OverlayConfig",
|
||||
"OverlayLayers",
|
||||
"PanSourceConfig",
|
||||
"PromptFrame",
|
||||
"PromptObservation",
|
||||
"PromptSourceConfig",
|
||||
"QuantResult",
|
||||
"ReferenceRegionConfig",
|
||||
"RoiShapeConfig",
|
||||
"SegmentationInput",
|
||||
"SmoothingConfig",
|
||||
"SplitDetectionError",
|
||||
"TicAnalysisConfig",
|
||||
"TicBaselineConfig",
|
||||
"TicCurveSeries",
|
||||
"TicCurves",
|
||||
"TicFittingConfig",
|
||||
"TicPlotMarker",
|
||||
"TicResult",
|
||||
"Trajectory",
|
||||
"VideoFile",
|
||||
"ViewRect",
|
||||
"ViewsConfig",
|
||||
"align_polygon_to_previous",
|
||||
"apply_dark_morphology",
|
||||
"apply_prompt_frame_position_constraints",
|
||||
"collect_prompt_observations",
|
||||
"compute_lognormal_perfusion_metrics",
|
||||
"compute_tic_for_video",
|
||||
"crop_frame",
|
||||
"dark_candidate_mask",
|
||||
"detect_baseline",
|
||||
"detect_split_x",
|
||||
"draw_polygon_overlay",
|
||||
"evaluate_lognormal",
|
||||
"evaluate_lognormal_derivative",
|
||||
"extract_color_prompt",
|
||||
"extract_dark_region_masks",
|
||||
"extract_dark_region_polygons",
|
||||
"fill_binary_holes",
|
||||
"fit_lognormal_curve",
|
||||
"from_tic_result",
|
||||
"fuse_prompt_prior_with_medsam",
|
||||
"keyframe_images",
|
||||
"largest_contour",
|
||||
"local_rect_within",
|
||||
"locate_anchor_frames",
|
||||
"locate_contrast_side",
|
||||
"map_source_mask_to_target",
|
||||
"mask_center",
|
||||
"mask_iou",
|
||||
"masks_to_smoothed_polygons",
|
||||
"mean_or_nan",
|
||||
"odd_window",
|
||||
"optional_positive_int",
|
||||
"pack_tic_curves",
|
||||
"pack_trajectory",
|
||||
"polygon_area",
|
||||
"polygon_center",
|
||||
"polygon_to_mask",
|
||||
"polygons_to_mask",
|
||||
"prepare_segmentation",
|
||||
"prompt_color_mask",
|
||||
"quant_summary",
|
||||
"quantify_case",
|
||||
"rect_from_values",
|
||||
"render_case_overlay",
|
||||
"render_overlay",
|
||||
"render_tic_png",
|
||||
"resample_closed_contour",
|
||||
"sample_frames_for_split",
|
||||
"score_dark_candidate",
|
||||
"select_prompt_frames",
|
||||
"select_reference_region",
|
||||
"split_dark_components",
|
||||
"split_views",
|
||||
"target_polygons_to_source",
|
||||
"tic_plot_markers",
|
||||
"triangular_smooth",
|
||||
"unpack_tic_curves",
|
||||
"unpack_trajectory",
|
||||
"view_offset",
|
||||
"view_rect",
|
||||
]
|
||||
|
||||
@@ -0,0 +1,419 @@
|
||||
"""第 8 步：把关键帧时刻映射到真实、有效且不重复的视频帧，共 7 个槽位。
|
||||
|
||||
- B0、D1、D2、D3、D4：从 `roi_perf`（病灶去暗区）曲线的 Lognormal 拟合上取。B0 是到达时刻再往前
|
||||
退 `baseline_frame_backoff_s` 秒，只在观测到的基线段里、到达之前取帧；D1 是拟合曲线上的到达
|
||||
时刻（design 0002 决定 6：不用医生标注的造影剂到达时间）；D2、D3、D4 是洗入拐点、峰值、
|
||||
洗出拐点。病灶两条曲线（`roi_all`、`roi_perf`）的拟合没有都过质量门时，这 5 个槽位全部未定位。
|
||||
D1–D4 在各自的容差窗口里找帧，并保持时间先后（不会出现 D2 的帧早于 D1 的帧）。
|
||||
- F1、F2：推注后第 `f1_anchor_s`、`late_phase.preferred_time_s` 秒。`injection_offset_s`（下称 X）
|
||||
是超声报告表的「推注开始时间s」，含义是视频第 0 秒时已经推注了 X 秒（design 0002 决定 5），
|
||||
所以推注后第 T 秒在视频第 T − X 秒。落在视频有效帧范围之外（早于第一帧或晚于最后一帧），
|
||||
这一槽位记为未定位，不拿别的帧补。
|
||||
|
||||
每个槽位交出帧号、视频时间、状态和未定位的原因；个别槽位未定位不算这一例失败。
|
||||
|
||||
照搬旧项目 `anchor_frames.py`，两处改动（任务书 02「四处要改的地方」第 3、4 条）：删掉 D1 取医生
|
||||
标注时刻的 `excel_video` 模式；推注时间的符号按上面的含义（旧 `excel_video` 模式按「推注发生在
|
||||
视频第 X 秒」算成了 T + X）。
|
||||
"""
|
||||
|
||||
import math
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._config import AnchorFramesConfig
|
||||
from cadx.ceus._geometry import AnyArray
|
||||
from cadx.ceus._tic_fitting import LognormalFitResult
|
||||
|
||||
ANCHOR_SLOTS = ("B0", "D1", "D2", "D3", "D4", "F1", "F2")
|
||||
|
||||
STATUS_OK = "ok"
|
||||
STATUS_ANCHOR_UNDEFINED = "anchor_undefined"
|
||||
STATUS_NO_FRAME_IN_WINDOW = "no_frame_in_window"
|
||||
STATUS_NO_ORDERED_FRAME = "no_ordered_frame"
|
||||
STATUS_NOT_COVERED = "not_covered"
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class AnchorFrameLocation:
|
||||
"""一个关键帧槽位的定位结果。
|
||||
|
||||
- `frame_index`、`frame_time_video_s`：定位到的帧号和它的视频时间，未定位为 None；
|
||||
- `time_s`：这一帧在推注时间轴上的时刻（视频时间 + X），没有推注时间或未定位为 None；
|
||||
- `status`：`ok` 或未定位的类别（`anchor_undefined`、`not_covered`、`no_frame_in_window`、
|
||||
`no_ordered_frame`）；`reason`：未定位的原因（中文），定位到时为 None；
|
||||
- `anchor_time_video_s`：要找的那个时刻（视频时间），没有时为 None；
|
||||
- `time_error_s`：帧时间减锚点时刻；`source_curve`：时刻来自哪条曲线或推注时钟。
|
||||
"""
|
||||
|
||||
slot: str
|
||||
frame_index: int | None
|
||||
time_s: float | None
|
||||
status: str
|
||||
reason: str | None
|
||||
anchor_time_video_s: float | None
|
||||
frame_time_video_s: float | None
|
||||
time_error_s: float | None
|
||||
source_curve: str | None
|
||||
|
||||
|
||||
def _finite_offset(injection_offset_s: float | None) -> float | None:
|
||||
"""推注时间 X：None 表示没有；有就必须是有限非负秒数，否则报错。"""
|
||||
if injection_offset_s is None:
|
||||
return None
|
||||
parsed = float(injection_offset_s)
|
||||
if not math.isfinite(parsed) or parsed < 0.0:
|
||||
raise ValueError(f"injection_offset_s 必须是有限非负秒数，收到 {injection_offset_s}")
|
||||
return parsed
|
||||
|
||||
|
||||
def _validate_axis(frame_indices: AnyArray, time_video_s: AnyArray) -> tuple[AnyArray, AnyArray]:
|
||||
"""帧号与视频时间必须等长、非空、严格递增。"""
|
||||
indices = np.asarray(frame_indices, dtype=np.int32)
|
||||
times = np.asarray(time_video_s, dtype=np.float64)
|
||||
if indices.ndim != 1 or times.shape != indices.shape or indices.size == 0:
|
||||
raise ValueError("frame_indices 与 time_video_s 必须是等长非空一维数组")
|
||||
if np.any(np.diff(indices) <= 0):
|
||||
raise ValueError("frame_indices 必须严格递增")
|
||||
if not np.isfinite(times).all() or np.any(np.diff(times) <= 0.0):
|
||||
raise ValueError("time_video_s 必须有限且严格递增")
|
||||
return indices, times
|
||||
|
||||
|
||||
def _nearest_frame(
|
||||
anchor_time_s: float,
|
||||
window_before_s: float,
|
||||
window_after_s: float,
|
||||
frame_indices: AnyArray,
|
||||
time_video_s: AnyArray,
|
||||
eligible_mask: AnyArray,
|
||||
) -> int | None:
|
||||
"""窗口 `[锚点 − before, 锚点 + after]` 内离锚点最近的可用帧；距离相同取更早的。"""
|
||||
eligible = (
|
||||
np.asarray(eligible_mask, dtype=bool)
|
||||
& (time_video_s >= anchor_time_s - window_before_s)
|
||||
& (time_video_s <= anchor_time_s + window_after_s)
|
||||
)
|
||||
positions = np.flatnonzero(eligible)
|
||||
if positions.size == 0:
|
||||
return None
|
||||
distances = np.abs(time_video_s[positions] - anchor_time_s)
|
||||
order = np.lexsort((time_video_s[positions], distances))
|
||||
return int(frame_indices[positions[int(order[0])]])
|
||||
|
||||
|
||||
type _State = tuple[int, float, tuple[int | None, ...]]
|
||||
|
||||
|
||||
def _ordered_anchor_frames(
|
||||
anchor_specs: tuple[tuple[str, float, float, float], ...],
|
||||
frame_indices: AnyArray,
|
||||
time_video_s: AnyArray,
|
||||
eligible_mask: AnyArray,
|
||||
) -> dict[str, int]:
|
||||
"""给 D1–D4 选帧：先保留最多的槽位，再让严格递增映射的总时间误差最小（动态规划）。"""
|
||||
eligible = np.asarray(eligible_mask, dtype=bool)
|
||||
states: dict[int, _State] = {-1: (0, 0.0, ())}
|
||||
for _slot, anchor_time_s, before_s, after_s in anchor_specs:
|
||||
positions = np.flatnonzero(
|
||||
eligible
|
||||
& (time_video_s >= anchor_time_s - before_s)
|
||||
& (time_video_s <= anchor_time_s + after_s)
|
||||
)
|
||||
next_states: dict[int, _State] = {}
|
||||
for previous, (count, cost, path) in states.items():
|
||||
skipped: _State = (count, cost, (*path, None))
|
||||
incumbent = next_states.get(previous)
|
||||
if incumbent is None or _state_key(skipped) < _state_key(incumbent):
|
||||
next_states[previous] = skipped
|
||||
for position_value in positions.tolist():
|
||||
position = int(position_value)
|
||||
if position <= previous:
|
||||
continue
|
||||
assigned: _State = (
|
||||
count + 1,
|
||||
cost + abs(float(time_video_s[position]) - anchor_time_s),
|
||||
(*path, position),
|
||||
)
|
||||
incumbent = next_states.get(position)
|
||||
if incumbent is None or _state_key(assigned) < _state_key(incumbent):
|
||||
next_states[position] = assigned
|
||||
states = next_states
|
||||
|
||||
_count, _cost, best_path = min(states.values(), key=_state_key)
|
||||
return {
|
||||
spec[0]: int(frame_indices[position])
|
||||
for spec, position in zip(anchor_specs, best_path, strict=True)
|
||||
if position is not None
|
||||
}
|
||||
|
||||
|
||||
def _state_key(state: _State) -> tuple[int, float, tuple[int, ...]]:
|
||||
"""动态规划的比较键：槽位多者优先，其次总误差小者，再按路径定序。"""
|
||||
count, cost, path = state
|
||||
sortable_path = tuple(-1 if position is None else position for position in path)
|
||||
return -count, cost, sortable_path
|
||||
|
||||
|
||||
def _video_time_for_frame(
|
||||
frame_index: int | None,
|
||||
frame_indices: AnyArray,
|
||||
time_video_s: AnyArray,
|
||||
) -> float | None:
|
||||
"""帧号对应的视频时间；帧号必须在时间轴里恰好出现一次。"""
|
||||
if frame_index is None:
|
||||
return None
|
||||
positions = np.flatnonzero(frame_indices == frame_index)
|
||||
if positions.size != 1:
|
||||
raise ValueError(f"帧号 {frame_index} 在时间轴中不是唯一值")
|
||||
return float(time_video_s[int(positions[0])])
|
||||
|
||||
|
||||
def _fit_eligible(fit: LognormalFitResult) -> AnyArray:
|
||||
"""拟合用到的有效且不重复的样本。"""
|
||||
return fit.valid_sample_mask & ~fit.duplicate_sample_mask
|
||||
|
||||
|
||||
def _span_text(start_s: float, end_s: float) -> str:
|
||||
return f"[{start_s:.2f}, {end_s:.2f}] 秒"
|
||||
|
||||
|
||||
def locate_anchor_frames(
|
||||
*,
|
||||
perf_fit: LognormalFitResult,
|
||||
perf_baseline_end_index: int,
|
||||
lesion_fit_ready: bool,
|
||||
frame_indices: AnyArray,
|
||||
time_video_s: AnyArray,
|
||||
injection_offset_s: float | None,
|
||||
cfg: AnchorFramesConfig,
|
||||
) -> dict[str, AnchorFrameLocation]:
|
||||
"""定位 7 个槽位，按 `ANCHOR_SLOTS` 的顺序交出。
|
||||
|
||||
`lesion_fit_ready` 是病灶两条曲线的拟合是否都过了质量门；`injection_offset_s` 是推注时间 X
|
||||
（秒，视频第 0 秒时已推注的时长），None 时 F1、F2 未定位。
|
||||
"""
|
||||
indices, times = _validate_axis(frame_indices, time_video_s)
|
||||
if perf_fit.valid_sample_mask.shape != times.shape:
|
||||
raise ValueError("perf_fit 的样本掩膜长度与时间轴不一致")
|
||||
injection_axis_offset_s = _finite_offset(injection_offset_s)
|
||||
arrival_time = perf_fit.d1_time_s
|
||||
locations: dict[str, AnchorFrameLocation] = {}
|
||||
|
||||
def record(
|
||||
slot: str,
|
||||
frame_index: int | None,
|
||||
status: str,
|
||||
reason: str | None,
|
||||
*,
|
||||
anchor_time_video_s: float | None,
|
||||
source_curve: str | None,
|
||||
) -> None:
|
||||
frame_time_video_s = _video_time_for_frame(frame_index, indices, times)
|
||||
locations[slot] = AnchorFrameLocation(
|
||||
slot=slot,
|
||||
frame_index=frame_index,
|
||||
time_s=(
|
||||
None
|
||||
if frame_time_video_s is None or injection_axis_offset_s is None
|
||||
else frame_time_video_s + injection_axis_offset_s
|
||||
),
|
||||
status=status,
|
||||
reason=None if status == STATUS_OK else reason,
|
||||
anchor_time_video_s=anchor_time_video_s,
|
||||
frame_time_video_s=frame_time_video_s,
|
||||
time_error_s=(
|
||||
None
|
||||
if frame_time_video_s is None or anchor_time_video_s is None
|
||||
else frame_time_video_s - anchor_time_video_s
|
||||
),
|
||||
source_curve=source_curve,
|
||||
)
|
||||
|
||||
model_ready = lesion_fit_ready and perf_fit.success
|
||||
if not lesion_fit_ready:
|
||||
fit_failure = "病灶两条曲线（roi_all、roi_perf）的 Lognormal 拟合没有都过质量门"
|
||||
elif not perf_fit.success:
|
||||
fit_failure = f"roi_perf 曲线的 Lognormal 拟合没过质量门（{perf_fit.failure_reason}）"
|
||||
else:
|
||||
fit_failure = ""
|
||||
|
||||
perf_eligible = _fit_eligible(perf_fit)
|
||||
|
||||
b0_ready = arrival_time is not None and perf_fit.baseline_source == "observed" and model_ready
|
||||
if b0_ready and arrival_time is not None:
|
||||
b0_time = float(arrival_time) - cfg.baseline_frame_backoff_s
|
||||
frame = _nearest_frame(
|
||||
b0_time,
|
||||
cfg.b0_window_s,
|
||||
cfg.b0_window_s,
|
||||
indices,
|
||||
times,
|
||||
perf_eligible
|
||||
& (np.arange(times.size) < perf_baseline_end_index)
|
||||
& (times < float(arrival_time)),
|
||||
)
|
||||
record(
|
||||
"B0",
|
||||
frame,
|
||||
STATUS_OK if frame is not None else STATUS_NO_FRAME_IN_WINDOW,
|
||||
f"B0 锚点 {b0_time:.2f} 秒的容差窗口 "
|
||||
f"{_span_text(b0_time - cfg.b0_window_s, b0_time + cfg.b0_window_s)} 内，"
|
||||
"观测到的基线段里、到达时刻之前没有可用帧",
|
||||
anchor_time_video_s=b0_time,
|
||||
source_curve="roi_perf",
|
||||
)
|
||||
else:
|
||||
if not model_ready:
|
||||
b0_reason = fit_failure
|
||||
elif arrival_time is None:
|
||||
b0_reason = "拟合曲线升不过到达门槛，没有到达时刻"
|
||||
else:
|
||||
b0_reason = "视频开头没有观测到稳定的增强前基线（基线是估计的），B0 不定义"
|
||||
record(
|
||||
"B0",
|
||||
None,
|
||||
STATUS_ANCHOR_UNDEFINED,
|
||||
b0_reason,
|
||||
anchor_time_video_s=None,
|
||||
source_curve="roi_perf",
|
||||
)
|
||||
|
||||
raw_specs = (
|
||||
("D1", arrival_time, 0.0, cfg.d1_window_after_s),
|
||||
("D2", perf_fit.d2_time_s, cfg.d2_window_s, cfg.d2_window_s),
|
||||
("D3", perf_fit.d3_time_s, cfg.d3_window_s, cfg.d3_window_s),
|
||||
("D4", perf_fit.d4_time_s, cfg.d4_window_s, cfg.d4_window_s),
|
||||
)
|
||||
eligible_times = times[perf_eligible]
|
||||
video_start_s = float(eligible_times[0]) if eligible_times.size else math.inf
|
||||
video_end_s = float(eligible_times[-1]) if eligible_times.size else -math.inf
|
||||
covered_specs = tuple(
|
||||
(slot, float(anchor_time), before, after)
|
||||
for slot, anchor_time, before, after in raw_specs
|
||||
if model_ready
|
||||
and anchor_time is not None
|
||||
and video_start_s <= float(anchor_time) <= video_end_s
|
||||
)
|
||||
dynamic_eligible = perf_eligible
|
||||
if model_ready:
|
||||
dynamic_eligible &= times < perf_fit.fit_end_s
|
||||
ordered = _ordered_anchor_frames(covered_specs, indices, times, dynamic_eligible)
|
||||
for slot, anchor_time, before, after in raw_specs:
|
||||
if not model_ready:
|
||||
record(
|
||||
slot,
|
||||
None,
|
||||
STATUS_ANCHOR_UNDEFINED,
|
||||
fit_failure,
|
||||
anchor_time_video_s=anchor_time,
|
||||
source_curve="roi_perf",
|
||||
)
|
||||
elif anchor_time is None:
|
||||
record(
|
||||
slot,
|
||||
None,
|
||||
STATUS_NOT_COVERED if slot == "D4" else STATUS_ANCHOR_UNDEFINED,
|
||||
"峰后观测不够长，或洗出拐点落在观测范围之外，拟合曲线上没有 D4"
|
||||
if slot == "D4"
|
||||
else f"拟合曲线上没有 {slot} 时刻（升不过到达门槛）",
|
||||
anchor_time_video_s=None,
|
||||
source_curve="roi_perf",
|
||||
)
|
||||
elif not video_start_s <= float(anchor_time) <= video_end_s:
|
||||
record(
|
||||
slot,
|
||||
None,
|
||||
STATUS_NOT_COVERED,
|
||||
f"{slot} 时刻 {float(anchor_time):.2f} 秒落在视频有效帧范围 "
|
||||
f"{_span_text(video_start_s, video_end_s)} 之外",
|
||||
anchor_time_video_s=float(anchor_time),
|
||||
source_curve="roi_perf",
|
||||
)
|
||||
elif slot in ordered:
|
||||
record(
|
||||
slot,
|
||||
ordered[slot],
|
||||
STATUS_OK,
|
||||
None,
|
||||
anchor_time_video_s=float(anchor_time),
|
||||
source_curve="roi_perf",
|
||||
)
|
||||
else:
|
||||
window = _span_text(float(anchor_time) - before, float(anchor_time) + after)
|
||||
has_candidate = bool(
|
||||
np.any(
|
||||
dynamic_eligible
|
||||
& (times >= float(anchor_time) - before)
|
||||
& (times <= float(anchor_time) + after)
|
||||
)
|
||||
)
|
||||
record(
|
||||
slot,
|
||||
None,
|
||||
STATUS_NO_ORDERED_FRAME if has_candidate else STATUS_NO_FRAME_IN_WINDOW,
|
||||
f"{slot} 容差窗口 {window} 内的候选帧都被时间更早的槽位占用，无法保持先后顺序"
|
||||
if has_candidate
|
||||
else f"{slot} 容差窗口 {window} 内没有可用帧（拟合区间内、不重复）",
|
||||
anchor_time_video_s=float(anchor_time),
|
||||
source_curve="roi_perf",
|
||||
)
|
||||
|
||||
globally_eligible = ~perf_fit.duplicate_sample_mask
|
||||
for slot, injection_time_s, tolerance_s in (
|
||||
("F1", cfg.f1_anchor_s, None),
|
||||
("F2", cfg.late_phase.preferred_time_s, cfg.late_phase.tolerance_s),
|
||||
):
|
||||
if injection_axis_offset_s is None:
|
||||
record(
|
||||
slot,
|
||||
None,
|
||||
STATUS_NOT_COVERED,
|
||||
"没有推注时间，推注后时刻换算不到视频时间",
|
||||
anchor_time_video_s=None,
|
||||
source_curve="injection_clock",
|
||||
)
|
||||
continue
|
||||
video_time_s = injection_time_s - injection_axis_offset_s
|
||||
fixed_eligible_times = times[globally_eligible]
|
||||
if (
|
||||
fixed_eligible_times.size == 0
|
||||
or video_time_s < float(fixed_eligible_times[0])
|
||||
or video_time_s > float(fixed_eligible_times[-1])
|
||||
):
|
||||
span = (
|
||||
"没有不重复的帧"
|
||||
if fixed_eligible_times.size == 0
|
||||
else f"视频有效帧范围 "
|
||||
f"{_span_text(float(fixed_eligible_times[0]), float(fixed_eligible_times[-1]))}"
|
||||
)
|
||||
record(
|
||||
slot,
|
||||
None,
|
||||
STATUS_NOT_COVERED,
|
||||
f"推注后第 {injection_time_s:g} 秒在视频第 {video_time_s:.2f} 秒"
|
||||
f"（推注时间 X={injection_axis_offset_s:g} 秒），落在{span}之外",
|
||||
anchor_time_video_s=video_time_s,
|
||||
source_curve="injection_clock",
|
||||
)
|
||||
continue
|
||||
window = float("inf") if tolerance_s is None else tolerance_s
|
||||
frame = _nearest_frame(
|
||||
video_time_s,
|
||||
window,
|
||||
window,
|
||||
indices,
|
||||
times,
|
||||
globally_eligible,
|
||||
)
|
||||
record(
|
||||
slot,
|
||||
frame,
|
||||
STATUS_OK if frame is not None else STATUS_NO_FRAME_IN_WINDOW,
|
||||
f"视频第 {video_time_s:.2f} 秒 ±{window:g} 秒内没有不重复的帧",
|
||||
anchor_time_video_s=video_time_s,
|
||||
source_curve="injection_clock",
|
||||
)
|
||||
|
||||
return {slot: locations[slot] for slot in ANCHOR_SLOTS}
|
||||
@@ -211,6 +211,10 @@ class InteriorDarkRegionConfig(BaseModel):
|
||||
threshold_offset: int
|
||||
min_area: int
|
||||
min_area_ratio: float
|
||||
max_regions: int | None
|
||||
max_contours_per_region: int | None
|
||||
preferred_area_ratio_range: list[float] | None
|
||||
prompt_exclude_dilate: int
|
||||
preserve_holes: bool
|
||||
fill_holes: bool
|
||||
morphology_kernel: int
|
||||
@@ -235,6 +239,25 @@ class InteriorDarkRegionConfig(BaseModel):
|
||||
raise ValueError(f"interior_dark_region.adaptive_block_size 必须 ≥3，收到 {validated}")
|
||||
return validated
|
||||
|
||||
@field_validator("max_regions", "max_contours_per_region")
|
||||
@classmethod
|
||||
def _optional_count_positive(cls, value: int | None) -> int | None:
|
||||
"""截断个数为 null（不限）或正整数；0 会把暗区全部截掉而不报错。"""
|
||||
if value is not None and value < 1:
|
||||
raise ValueError(f"interior_dark_region 的截断个数必须为 null 或 ≥1，收到 {value}")
|
||||
return value
|
||||
|
||||
@field_validator("preferred_area_ratio_range")
|
||||
@classmethod
|
||||
def _area_range_pair(cls, value: list[float] | None) -> list[float] | None:
|
||||
"""偏好面积比区间为 null 或 [lo, hi]，0 ≤ lo ≤ hi ≤ 1。"""
|
||||
if value is not None and not (len(value) == 2 and 0.0 <= value[0] <= value[1] <= 1.0):
|
||||
raise ValueError(
|
||||
"interior_dark_region.preferred_area_ratio_range 必须为 null 或 "
|
||||
f"0≤lo≤hi≤1 的 [lo, hi]，收到 {value}"
|
||||
)
|
||||
return value
|
||||
|
||||
@field_validator("percentile")
|
||||
@classmethod
|
||||
def _percentile_within_bounds(cls, value: float) -> float:
|
||||
|
||||
@@ -0,0 +1,670 @@
|
||||
"""第 5 步：在病灶 ROI 内找暗区（不进造影剂、一直发黑的囊液部分）。
|
||||
|
||||
哪种阈值法在当前帧最可信事先定不了，所以 `threshold_method="auto_ensemble"` 时对
|
||||
`candidate_methods` 里的每种阈值法各出一份暗块候选，用统一评分（明暗对比 + 面积分
|
||||
− 碎片惩罚 + 与上一帧的时间一致性）择优。
|
||||
|
||||
- `dark_candidate_mask`：按方法在 ROI 灰度内出候选暗掩膜（otsu / percentile / absolute /
|
||||
adaptive / kmeans / clahe_otsu）；`threshold_offset` 只对前三种生效；末尾过形态学；
|
||||
- `apply_dark_morphology`：闭、开运算各一次（`morphology_kernel` 方核，≤1 跳过）；
|
||||
`fill_holes` 且不 `preserve_holes` 时填洞；最后与采样掩膜取交；
|
||||
- `split_dark_components`：连通域拆分，面积门槛 `max(min_area, ROI 面积 × min_area_ratio)`；
|
||||
`preserve_holes` 时把洞挖回；按面积降序，`max_regions` 截断；
|
||||
- `score_dark_candidate`：评分（三个常量 64.0 / 4.0 / 2.0）；
|
||||
- `extract_dark_region_masks`：单帧——剔绿框色 → 各候选法择优 → 拆分；
|
||||
- `extract_dark_region_polygons`：逐帧编排，`enabled=False` 时每帧空列表（仍保留帧对齐）；
|
||||
否则各暗块取显著轮廓、等弧长重采样到 `contour_points` 点。
|
||||
|
||||
输入是造影那一半画面（`analysis_view: source`）的逐帧 BGR 与 ROI 多边形，输出同一坐标系的暗区
|
||||
多边形。照搬旧项目 `dark_region.py`，去掉产品层；旧代码用 `getattr(..., 默认值)` 读的四个旧 schema
|
||||
里没有的键（`max_regions`、`max_contours_per_region`、`preferred_area_ratio_range`、
|
||||
`prompt_exclude_dilate`）改成配置字段、值照旧；剔绿用的 HSV 阈值只取 `prompt_source`，
|
||||
旧代码在没传时回退的写死阈值删掉。
|
||||
"""
|
||||
|
||||
from collections.abc import Iterable, Sequence
|
||||
from typing import Any, Protocol, cast
|
||||
|
||||
import cv2
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._geometry import (
|
||||
AnyArray,
|
||||
Mask,
|
||||
mask_iou,
|
||||
polygon_to_mask,
|
||||
resample_closed_contour,
|
||||
)
|
||||
|
||||
# 评分公式的三个常量（照参考组件 P:1296 / 1308 / 1310）。
|
||||
_CONTRAST_NORMALIZER = 64.0 # 明暗对比归一化除数
|
||||
_AREA_SLOPE = 4.0 # preferred_area_ratio_range 之外面积分的线性衰减斜率
|
||||
_AREA_SCORE_GAIN = 2.0 # 没有 preferred_area_ratio_range 时面积分的系数
|
||||
|
||||
_EMPTY_SCORE = -1e9 # 空候选的哨兵分：择优时永远不会被选中
|
||||
|
||||
|
||||
class _DarkRegionConfigLike(Protocol):
|
||||
"""暗区配置（`InteriorDarkRegionConfig` 符合）。"""
|
||||
|
||||
@property
|
||||
def enabled(self) -> bool: ...
|
||||
@property
|
||||
def exclude_prompt_color(self) -> bool: ...
|
||||
@property
|
||||
def prompt_exclude_dilate(self) -> int: ...
|
||||
@property
|
||||
def threshold_method(self) -> str: ...
|
||||
@property
|
||||
def candidate_methods(self) -> Sequence[str]: ...
|
||||
@property
|
||||
def threshold_offset(self) -> int: ...
|
||||
@property
|
||||
def min_area(self) -> int: ...
|
||||
@property
|
||||
def min_area_ratio(self) -> float: ...
|
||||
@property
|
||||
def max_regions(self) -> int | None: ...
|
||||
@property
|
||||
def max_contours_per_region(self) -> int | None: ...
|
||||
@property
|
||||
def preferred_area_ratio_range(self) -> list[float] | None: ...
|
||||
@property
|
||||
def preserve_holes(self) -> bool: ...
|
||||
@property
|
||||
def fill_holes(self) -> bool: ...
|
||||
@property
|
||||
def morphology_kernel(self) -> int: ...
|
||||
@property
|
||||
def fragment_penalty_weight(self) -> float: ...
|
||||
@property
|
||||
def temporal_consistency_weight(self) -> float: ...
|
||||
@property
|
||||
def contour_points(self) -> int: ...
|
||||
@property
|
||||
def percentile(self) -> float: ...
|
||||
@property
|
||||
def threshold(self) -> int: ...
|
||||
@property
|
||||
def adaptive_block_size(self) -> int: ...
|
||||
@property
|
||||
def adaptive_c(self) -> float: ...
|
||||
@property
|
||||
def kmeans_max_iter(self) -> int: ...
|
||||
@property
|
||||
def kmeans_attempts(self) -> int: ...
|
||||
@property
|
||||
def clahe_clip_limit(self) -> float: ...
|
||||
@property
|
||||
def clahe_tile_grid_size(self) -> int: ...
|
||||
|
||||
|
||||
class _PromptHsvConfigLike(Protocol):
|
||||
"""剔绿框色用的 HSV 阈值（`PromptSourceConfig` 符合）。"""
|
||||
|
||||
@property
|
||||
def hsv_lower(self) -> list[int]: ...
|
||||
@property
|
||||
def hsv_upper(self) -> list[int]: ...
|
||||
|
||||
|
||||
def fill_binary_holes(mask: AnyArray) -> AnyArray:
|
||||
"""对二值 mask 填充内部空洞（参考 P:1112-1119）。
|
||||
|
||||
参数：
|
||||
mask: ``(h, w)`` 的二值数组（bool 或非零即真）。
|
||||
|
||||
返回：
|
||||
填洞后的 ``dtype=bool`` mask（外接背景之外的所有像素均置 True）。
|
||||
|
||||
关键实现：
|
||||
从角点 ``(0, 0)`` 泛洪标记"可从外部到达的背景"，取反即得"被前景包围的洞"，
|
||||
与原前景按位或合并。要求 ``(0, 0)`` 像素为背景（暗区不会触达画面角点，成立）。
|
||||
"""
|
||||
filled = np.asarray(mask, dtype=np.uint8) * np.uint8(255)
|
||||
flood = filled.copy()
|
||||
height, width = flood.shape
|
||||
flood_mask = np.zeros((height + 2, width + 2), dtype=np.uint8)
|
||||
cv2.floodFill(flood, flood_mask, (0, 0), 255)
|
||||
holes = np.asarray(cv2.bitwise_not(flood), dtype=np.uint8)
|
||||
return (filled | holes) > 0
|
||||
|
||||
|
||||
def apply_dark_morphology(
|
||||
dark: AnyArray,
|
||||
sample_mask: AnyArray,
|
||||
dark_cfg: _DarkRegionConfigLike,
|
||||
) -> AnyArray:
|
||||
"""对候选暗 mask 做形态学闭 / 开运算并约束回 sample_mask（参考 P:1122-1136）。
|
||||
|
||||
参数：
|
||||
dark: 候选暗 mask（bool）。
|
||||
sample_mask: 分析采样掩膜（ROI 内、已剔绿的有效区）；结果与其取交。
|
||||
dark_cfg: 暗区配置（取 ``morphology_kernel`` / ``fill_holes`` / ``preserve_holes``）。
|
||||
|
||||
返回：
|
||||
形态学处理后、且**仅落在 sample_mask 内**的 ``dtype=bool`` mask。
|
||||
|
||||
关键实现：
|
||||
方核大小 ``> 1`` 时先 CLOSE（补小缺口）再 OPEN（去毛刺）各一次；
|
||||
仅当 ``fill_holes=True 且 preserve_holes=False`` 才填洞（保洞优先于填洞）；
|
||||
末尾与 sample_mask 取交，杜绝形态学膨胀外溢 ROI。
|
||||
"""
|
||||
kernel_size = int(dark_cfg.morphology_kernel)
|
||||
if kernel_size > 1:
|
||||
kernel = np.ones((kernel_size, kernel_size), dtype=np.uint8)
|
||||
dark_u8 = dark.astype(np.uint8) * 255
|
||||
dark_u8 = cv2.morphologyEx(dark_u8, cv2.MORPH_CLOSE, kernel, iterations=1)
|
||||
dark_u8 = cv2.morphologyEx(dark_u8, cv2.MORPH_OPEN, kernel, iterations=1)
|
||||
dark = dark_u8 > 0
|
||||
if bool(dark_cfg.fill_holes) and not bool(dark_cfg.preserve_holes):
|
||||
dark = fill_binary_holes(dark)
|
||||
return dark & sample_mask
|
||||
|
||||
|
||||
def _otsu_threshold(values: AnyArray) -> float:
|
||||
"""对一维灰度样本求 Otsu 阈值（参考 P:1150-1155 / 1206-1211）。"""
|
||||
threshold, _ = cv2.threshold(values.reshape(-1, 1), 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
|
||||
return float(threshold)
|
||||
|
||||
|
||||
def _candidate_otsu(gray: AnyArray, roi_values: AnyArray) -> AnyArray:
|
||||
"""Otsu 法候选：ROI 灰度 Otsu 阈值，``gray <= 阈值`` 为暗（P:1149-1156）。"""
|
||||
return gray <= _otsu_threshold(roi_values)
|
||||
|
||||
|
||||
def _candidate_percentile(
|
||||
gray: AnyArray, roi_values: AnyArray, dark_cfg: _DarkRegionConfigLike
|
||||
) -> AnyArray:
|
||||
"""分位法候选：ROI 灰度的 ``percentile`` 分位作阈值（P:1157-1159）。"""
|
||||
threshold = float(np.percentile(roi_values, float(dark_cfg.percentile)))
|
||||
return gray <= threshold
|
||||
|
||||
|
||||
def _candidate_absolute(gray: AnyArray, dark_cfg: _DarkRegionConfigLike) -> AnyArray:
|
||||
"""绝对阈法候选：固定灰度阈 ``threshold``（P:1160-1162）。"""
|
||||
return gray <= float(dark_cfg.threshold)
|
||||
|
||||
|
||||
def _candidate_adaptive(gray: AnyArray, dark_cfg: _DarkRegionConfigLike) -> AnyArray:
|
||||
"""自适应阈法候选：高斯邻域自适应 + ``THRESH_BINARY_INV``（P:1163-1176）。
|
||||
|
||||
block_size 强制奇数且 ``>= 3``（偶数加 1）；inverse 二值化使暗部为前景。
|
||||
"""
|
||||
block_size = int(dark_cfg.adaptive_block_size)
|
||||
if block_size % 2 == 0:
|
||||
block_size += 1
|
||||
block_size = max(3, block_size)
|
||||
adaptive = cv2.adaptiveThreshold(
|
||||
gray,
|
||||
255,
|
||||
cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
|
||||
cv2.THRESH_BINARY_INV,
|
||||
block_size,
|
||||
float(dark_cfg.adaptive_c),
|
||||
)
|
||||
return adaptive > 0
|
||||
|
||||
|
||||
def _candidate_kmeans(
|
||||
sample_mask: AnyArray, roi_values: AnyArray, dark_cfg: _DarkRegionConfigLike
|
||||
) -> AnyArray:
|
||||
"""K-means 二聚类候选：ROI 灰度聚成两簇，低均值簇为暗（P:1177-1196）。
|
||||
|
||||
EPS=0.5 硬编码（与参考一致）；样本不足 2 个时返回全 False。结果仅写回 sample_mask 内。
|
||||
"""
|
||||
values = roi_values.reshape(-1, 1).astype(np.float32)
|
||||
if len(values) < 2:
|
||||
return np.zeros_like(sample_mask, dtype=bool)
|
||||
criteria = (
|
||||
cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER,
|
||||
int(dark_cfg.kmeans_max_iter),
|
||||
0.5,
|
||||
)
|
||||
# cv2 的类型存根不接受 bestLabels=None（运行时接受，表示不给初始标签）。
|
||||
kmeans = cast(Any, cv2.kmeans)
|
||||
_, labels, centers = cast(
|
||||
tuple[float, AnyArray, AnyArray],
|
||||
kmeans(
|
||||
values,
|
||||
2,
|
||||
None,
|
||||
criteria,
|
||||
int(dark_cfg.kmeans_attempts),
|
||||
cv2.KMEANS_PP_CENTERS,
|
||||
),
|
||||
)
|
||||
dark_label = int(np.argmin(centers.reshape(-1)))
|
||||
dark = np.zeros_like(sample_mask, dtype=bool)
|
||||
dark[sample_mask] = labels.reshape(-1) == dark_label
|
||||
return dark
|
||||
|
||||
|
||||
def _candidate_clahe_otsu(
|
||||
gray: AnyArray, sample_mask: AnyArray, dark_cfg: _DarkRegionConfigLike
|
||||
) -> AnyArray:
|
||||
"""CLAHE 增强后 Otsu 候选：先局部直方图均衡再在 ROI 内 Otsu（P:1197-1212）。"""
|
||||
clahe = cv2.createCLAHE(
|
||||
clipLimit=float(dark_cfg.clahe_clip_limit),
|
||||
tileGridSize=(
|
||||
int(dark_cfg.clahe_tile_grid_size),
|
||||
int(dark_cfg.clahe_tile_grid_size),
|
||||
),
|
||||
)
|
||||
enhanced = clahe.apply(gray)
|
||||
threshold = _otsu_threshold(enhanced[sample_mask])
|
||||
return enhanced <= threshold
|
||||
|
||||
|
||||
def _apply_threshold_offset(
|
||||
dark: AnyArray,
|
||||
gray: AnyArray,
|
||||
sample_mask: AnyArray,
|
||||
method: str,
|
||||
dark_cfg: _DarkRegionConfigLike,
|
||||
) -> AnyArray:
|
||||
"""对 otsu / percentile / absolute 三法施加 ``threshold_offset``（参考 P:1216-1233）。
|
||||
|
||||
其余方法（adaptive / kmeans / clahe_otsu）不支持偏移，原样返回。偏移为 0 时也原样返回。
|
||||
重新以"基准阈值 + 偏移"二值化（不是在已二值化结果上平移），与参考严格等价。
|
||||
"""
|
||||
threshold_offset = float(dark_cfg.threshold_offset)
|
||||
if threshold_offset == 0 or method not in {"otsu", "percentile", "absolute"}:
|
||||
return dark
|
||||
base_values = gray[sample_mask]
|
||||
if method == "otsu":
|
||||
return gray <= (_otsu_threshold(base_values) + threshold_offset)
|
||||
if method == "percentile":
|
||||
base = float(np.percentile(base_values, float(dark_cfg.percentile)))
|
||||
return gray <= (base + threshold_offset)
|
||||
# absolute
|
||||
return gray <= (float(dark_cfg.threshold) + threshold_offset)
|
||||
|
||||
|
||||
def dark_candidate_mask(
|
||||
gray: AnyArray,
|
||||
sample_mask: AnyArray,
|
||||
method: str,
|
||||
dark_cfg: _DarkRegionConfigLike,
|
||||
) -> AnyArray:
|
||||
"""按指定阈值法在 ROI 灰度内出候选暗 mask（参考 P:1139-1235）。
|
||||
|
||||
参数：
|
||||
gray: ``(h, w)`` 单通道 uint8 灰度图（整帧）。
|
||||
sample_mask: 分析采样掩膜（ROI 内、已剔绿）；阈值统计只在其内取灰度。
|
||||
method: 阈值法名（otsu / percentile / absolute / adaptive / kmeans / clahe_otsu）。
|
||||
dark_cfg: 暗区配置（提供各法子参与 ``threshold_offset``）。
|
||||
|
||||
返回：
|
||||
过形态学并约束回 sample_mask 后的 ``dtype=bool`` 暗 mask；
|
||||
sample_mask 为空时返回全 False。
|
||||
|
||||
关键实现：
|
||||
各法分派到独立 helper（控复杂度）；除 adaptive/clahe 在全图取阈外，统计量
|
||||
均取自 ``gray[sample_mask]``。``threshold_offset`` 仅对 otsu/percentile/absolute
|
||||
生效（见 ``_apply_threshold_offset``）。末尾 ``(dark & sample_mask)`` 再过
|
||||
``apply_dark_morphology``（其内会再次与 sample_mask 取交）。未知方法直接报错（不兜底）。
|
||||
"""
|
||||
roi_values = gray[sample_mask]
|
||||
if len(roi_values) == 0:
|
||||
return np.zeros_like(sample_mask, dtype=bool)
|
||||
|
||||
if method == "otsu":
|
||||
dark = _candidate_otsu(gray, roi_values)
|
||||
elif method == "percentile":
|
||||
dark = _candidate_percentile(gray, roi_values, dark_cfg)
|
||||
elif method == "absolute":
|
||||
dark = _candidate_absolute(gray, dark_cfg)
|
||||
elif method == "adaptive":
|
||||
dark = _candidate_adaptive(gray, dark_cfg)
|
||||
elif method == "kmeans":
|
||||
dark = _candidate_kmeans(sample_mask, roi_values, dark_cfg)
|
||||
elif method == "clahe_otsu":
|
||||
dark = _candidate_clahe_otsu(gray, sample_mask, dark_cfg)
|
||||
else:
|
||||
raise RuntimeError(f"不支持的 interior_dark_region.threshold_method: {method}")
|
||||
|
||||
dark = _apply_threshold_offset(dark, gray, sample_mask, method, dark_cfg)
|
||||
return apply_dark_morphology(dark & sample_mask, sample_mask, dark_cfg)
|
||||
|
||||
|
||||
def split_dark_components(
|
||||
dark: AnyArray,
|
||||
roi_mask: AnyArray,
|
||||
dark_cfg: _DarkRegionConfigLike,
|
||||
) -> list[AnyArray]:
|
||||
"""把候选暗 mask 拆成达标连通域列表（参考 P:1238-1274）。
|
||||
|
||||
参数：
|
||||
dark: 候选暗 mask（bool）。
|
||||
roi_mask: ROI 掩膜（**未剔绿的原始 ROI**，用于面积阈与组件裁剪）。
|
||||
dark_cfg: 暗区配置（``min_area`` / ``min_area_ratio`` / ``preserve_holes`` /
|
||||
``max_regions``）。
|
||||
|
||||
返回：
|
||||
按面积降序排列的连通域 ``dtype=bool`` mask 列表；
|
||||
每个面积须 ``>= max(min_area, ROI面积×min_area_ratio)``；``max_regions`` 截断。
|
||||
|
||||
关键实现：
|
||||
组件阈 ``max(min_area, roi_mask.sum()×min_area_ratio)``（``min_area_ratio>0`` 才生效）。
|
||||
``findContours(RETR_CCOMP)`` 取外 / 内两层；仅遍历外轮廓（``hierarchy[3]<0``），
|
||||
实心填充后——若 ``preserve_holes``——再把其子轮廓（洞）挖回 0，最后与 roi_mask 取交。
|
||||
``max_regions`` 为 None 时不限。
|
||||
"""
|
||||
min_area = float(dark_cfg.min_area)
|
||||
min_area_ratio = float(dark_cfg.min_area_ratio)
|
||||
if min_area_ratio > 0:
|
||||
min_area = max(min_area, float(roi_mask.sum()) * min_area_ratio)
|
||||
|
||||
contours, hierarchy = cv2.findContours(
|
||||
dark.astype(np.uint8), cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE
|
||||
)
|
||||
# 没有轮廓时 cv2 交回的 hierarchy 是 None。
|
||||
if not contours:
|
||||
return []
|
||||
hierarchy = np.asarray(hierarchy)[0]
|
||||
preserve_holes = bool(dark_cfg.preserve_holes)
|
||||
components: list[tuple[int, AnyArray]] = []
|
||||
for idx, _contour in enumerate(contours):
|
||||
if hierarchy[idx][3] >= 0:
|
||||
continue
|
||||
out = np.zeros_like(dark, dtype=np.uint8)
|
||||
cv2.drawContours(out, contours, idx, 255, thickness=-1)
|
||||
if preserve_holes:
|
||||
child = int(hierarchy[idx][2])
|
||||
while child >= 0:
|
||||
cv2.drawContours(out, contours, child, 0, thickness=-1)
|
||||
child = int(hierarchy[child][0])
|
||||
component = (out > 0) & roi_mask
|
||||
area = int(component.sum())
|
||||
if area >= min_area:
|
||||
components.append((area, component))
|
||||
|
||||
components.sort(key=lambda item: item[0], reverse=True)
|
||||
max_regions = dark_cfg.max_regions
|
||||
selected = components if max_regions is None else components[:max_regions]
|
||||
return [component for _, component in selected]
|
||||
|
||||
|
||||
def _area_score(area_ratio: float, dark_cfg: _DarkRegionConfigLike) -> float:
|
||||
"""面积分项（参考 P:1298-1310，含斜率 4.0 / 系数 2.0 两个硬编码常量）。
|
||||
|
||||
无 ``preferred_area_ratio_range`` 时 ``min(r, 1−r) × 2.0``（偏好"占一半"的暗块）；
|
||||
给定区间则命中得 1.0，越界按 ``max(0, 1 − 与最近边距×4.0)`` 线性衰减。
|
||||
"""
|
||||
preferred_range = dark_cfg.preferred_area_ratio_range
|
||||
if preferred_range is None:
|
||||
return min(area_ratio, 1.0 - area_ratio) * _AREA_SCORE_GAIN
|
||||
preferred_low = float(preferred_range[0])
|
||||
preferred_high = float(preferred_range[1])
|
||||
if preferred_low <= area_ratio <= preferred_high:
|
||||
return 1.0
|
||||
nearest = preferred_low if area_ratio < preferred_low else preferred_high
|
||||
return max(0.0, 1.0 - abs(area_ratio - nearest) * _AREA_SLOPE)
|
||||
|
||||
|
||||
def score_dark_candidate(
|
||||
gray: AnyArray,
|
||||
sample_mask: AnyArray,
|
||||
component_masks: list[AnyArray],
|
||||
previous_mask: AnyArray | None,
|
||||
dark_cfg: _DarkRegionConfigLike,
|
||||
) -> float:
|
||||
"""对一份候选暗块集打分，用于 auto_ensemble 择优（参考 P:1277-1319）。
|
||||
|
||||
参数：
|
||||
gray: 整帧灰度图。
|
||||
sample_mask: 分析采样掩膜（ROI 内、已剔绿）。
|
||||
component_masks: 候选暗块 mask 列表（``split_dark_components`` 的输出）。
|
||||
previous_mask: 上一帧选中暗部并集（无则 None），用于时间一致性。
|
||||
dark_cfg: 暗区配置（``fragment_penalty_weight`` / ``temporal_consistency_weight``
|
||||
/ 可选 ``preferred_area_ratio_range``）。
|
||||
|
||||
返回：
|
||||
``对比分 + 面积分 + 时间一致性×权重 − 碎片惩罚``；空候选 / 暗或亮像素为空
|
||||
时返回哨兵 ``-1e9``（择优时永不胜出）。
|
||||
|
||||
关键实现（四项逐行核对参考）：
|
||||
- 对比分：``max(0, (median(亮) − median(暗)) / 64.0)``（除数 64.0 硬编码）。
|
||||
- 面积分：见 ``_area_score``（斜率 4.0 / 系数 2.0）。
|
||||
- 碎片惩罚：``fragment_penalty_weight × max(0, 块数 − 1)``（块数 ≤1 不罚）。
|
||||
- 时间一致性：仅当 ``previous_mask`` 同形且非空，取 ``IoU(union, previous_mask)``
|
||||
再乘 ``temporal_consistency_weight``。
|
||||
"""
|
||||
if not component_masks:
|
||||
return _EMPTY_SCORE
|
||||
union = np.zeros_like(sample_mask, dtype=bool)
|
||||
for component in component_masks:
|
||||
union |= component
|
||||
dark_values = gray[union & sample_mask]
|
||||
bright_values = gray[(~union) & sample_mask]
|
||||
if len(dark_values) == 0 or len(bright_values) == 0:
|
||||
return _EMPTY_SCORE
|
||||
|
||||
contrast_score = max(
|
||||
0.0,
|
||||
(float(np.median(bright_values)) - float(np.median(dark_values))) / _CONTRAST_NORMALIZER,
|
||||
)
|
||||
area_ratio = float(union.sum()) / max(1.0, float(sample_mask.sum()))
|
||||
area_score = _area_score(area_ratio, dark_cfg)
|
||||
|
||||
fragment_penalty = float(dark_cfg.fragment_penalty_weight) * max(0, len(component_masks) - 1)
|
||||
temporal_score = 0.0
|
||||
if previous_mask is not None and previous_mask.shape == union.shape and previous_mask.any():
|
||||
temporal_score = mask_iou(union, previous_mask)
|
||||
temporal_weight = float(dark_cfg.temporal_consistency_weight)
|
||||
return contrast_score + area_score + temporal_score * temporal_weight - fragment_penalty
|
||||
|
||||
|
||||
def _prompt_exclusion_mask(
|
||||
target_bgr: AnyArray,
|
||||
dark_cfg: _DarkRegionConfigLike,
|
||||
prompt_cfg: _PromptHsvConfigLike,
|
||||
) -> Mask:
|
||||
"""要从分析中剔除的绿框色像素（参考 P:1334-1345）。
|
||||
|
||||
HSV 阈值取 `prompt_cfg`（与找绿框共用同一组阈值）；`prompt_exclude_dilate` 大于 1 时按方核
|
||||
膨胀一次。
|
||||
"""
|
||||
lower = np.asarray(prompt_cfg.hsv_lower, dtype=np.uint8)
|
||||
upper = np.asarray(prompt_cfg.hsv_upper, dtype=np.uint8)
|
||||
hsv = cv2.cvtColor(target_bgr, cv2.COLOR_BGR2HSV)
|
||||
prompt_pixels: Mask = np.asarray(cv2.inRange(hsv, lower, upper)) > 0
|
||||
kernel_size = int(dark_cfg.prompt_exclude_dilate)
|
||||
if kernel_size > 1:
|
||||
kernel = np.ones((kernel_size, kernel_size), dtype=np.uint8)
|
||||
prompt_pixels = cv2.dilate(prompt_pixels.astype(np.uint8) * 255, kernel, iterations=1) > 0
|
||||
return prompt_pixels
|
||||
|
||||
|
||||
def extract_dark_region_masks(
|
||||
target_bgr: AnyArray,
|
||||
roi_polygon: AnyArray,
|
||||
dark_cfg: _DarkRegionConfigLike,
|
||||
previous_mask: AnyArray | None,
|
||||
prompt_cfg: _PromptHsvConfigLike,
|
||||
) -> tuple[list[AnyArray], AnyArray | None]:
|
||||
"""单帧暗区提取：剔绿 → 候选法（auto_ensemble 择优）→ 拆分（参考 P:1322-1373）。
|
||||
|
||||
参数：
|
||||
target_bgr: 分析坐标系下的整帧 BGR uint8 图（``analysis_view='source'`` 即左视图局部）。
|
||||
roi_polygon: 该帧 ROI 多边形 ``(N, 2)``（分析坐标系）。
|
||||
dark_cfg: 暗区配置。
|
||||
previous_mask: 上一帧选中暗部并集（供时间一致性评分），首帧传 None。
|
||||
prompt_cfg: 剔绿用 HSV 阈值（``exclude_prompt_color=True`` 时使用）。
|
||||
|
||||
返回：
|
||||
``(选中暗块 mask 列表, 暗块并集 mask)``；无有效 ROI / 无暗块时返回 ``([], None)``。
|
||||
|
||||
关键实现：
|
||||
``sample_mask`` 由 ROI mask 剔除绿框像素得到（评分统计在其内），但
|
||||
``split_dark_components`` 的面积阈与裁剪仍用**原始 roi_mask**（与参考 P:1363 一致）。
|
||||
``threshold_method='auto_ensemble'`` 时遍历 ``candidate_methods`` 各出候选 → 拆分 →
|
||||
``score_dark_candidate`` 取最高分者；否则单法。返回的并集供下一帧 previous_mask 用。
|
||||
"""
|
||||
roi_mask = polygon_to_mask(roi_polygon, target_bgr.shape[0], target_bgr.shape[1])
|
||||
if not roi_mask.any():
|
||||
return [], None
|
||||
|
||||
gray = cv2.cvtColor(target_bgr, cv2.COLOR_BGR2GRAY)
|
||||
sample_mask = roi_mask.copy()
|
||||
if bool(dark_cfg.exclude_prompt_color):
|
||||
sample_mask &= ~_prompt_exclusion_mask(target_bgr, dark_cfg, prompt_cfg)
|
||||
|
||||
if len(gray[sample_mask]) == 0:
|
||||
return [], None
|
||||
|
||||
method = dark_cfg.threshold_method
|
||||
candidate_methods = list(dark_cfg.candidate_methods) if method == "auto_ensemble" else [method]
|
||||
components, union = _select_best_candidate(
|
||||
gray, sample_mask, roi_mask, candidate_methods, previous_mask, dark_cfg
|
||||
)
|
||||
if not components:
|
||||
return [], None
|
||||
return components, union
|
||||
|
||||
|
||||
def _select_best_candidate(
|
||||
gray: AnyArray,
|
||||
sample_mask: AnyArray,
|
||||
roi_mask: AnyArray,
|
||||
candidate_methods: list[str],
|
||||
previous_mask: AnyArray | None,
|
||||
dark_cfg: _DarkRegionConfigLike,
|
||||
) -> tuple[list[AnyArray], AnyArray | None]:
|
||||
"""对每个候选阈值法出暗块 + 评分，取最高分者（参考 P:1360-1373）。
|
||||
|
||||
参数：
|
||||
gray: 整帧灰度图。
|
||||
sample_mask: 分析采样掩膜（ROI 内、已剔绿）。
|
||||
roi_mask: 原始 ROI 掩膜（供 ``split_dark_components`` 面积阈与裁剪）。
|
||||
candidate_methods: 待评估的阈值法名列表（auto_ensemble 即全部，否则单元素）。
|
||||
previous_mask: 上一帧选中暗部并集（时间一致性评分用）。
|
||||
dark_cfg: 暗区配置。
|
||||
|
||||
返回：
|
||||
``(最优候选的暗块列表, 其并集 mask)``；无任何候选（列表空）时并集为 None。
|
||||
|
||||
关键实现：
|
||||
逐候选 ``dark_candidate_mask → split_dark_components → score_dark_candidate``，
|
||||
``score`` 严格大于现任最优才更换（首个候选无条件入选），并随即缓存其并集，
|
||||
避免对落选候选白算并集。空候选的评分哨兵 ``-1e9`` 保证其永不胜出。
|
||||
"""
|
||||
best: tuple[float, list[AnyArray], AnyArray] | None = None
|
||||
for candidate_method in candidate_methods:
|
||||
dark = dark_candidate_mask(gray, sample_mask, candidate_method, dark_cfg)
|
||||
components = split_dark_components(dark, roi_mask, dark_cfg)
|
||||
score = score_dark_candidate(gray, sample_mask, components, previous_mask, dark_cfg)
|
||||
if best is None or score > best[0]:
|
||||
union = np.zeros_like(sample_mask, dtype=bool)
|
||||
for component in components:
|
||||
union |= component
|
||||
best = (score, components, union)
|
||||
if best is None:
|
||||
return [], None
|
||||
return best[1], best[2]
|
||||
|
||||
|
||||
def extract_dark_region_polygons(
|
||||
analysis_frames: Iterable[AnyArray],
|
||||
roi_polygons: Sequence[AnyArray],
|
||||
dark_cfg: _DarkRegionConfigLike,
|
||||
prompt_cfg: _PromptHsvConfigLike,
|
||||
) -> list[list[AnyArray]]:
|
||||
"""逐帧暗区多边形提取编排（参考 P:1376-1466，按本项目可测性裁剪视图映射）。
|
||||
|
||||
参数：
|
||||
analysis_frames: 分析坐标系下的逐帧整帧 BGR uint8 可迭代对象（``analysis_view='source'``
|
||||
即造影那一半的局部坐标；调用方已解码并切好）。**逐帧消费、不留**，
|
||||
所以可以直接喂一个解码器的迭代器。
|
||||
roi_polygons: 与帧序列等长的逐帧 ROI 多边形（同一坐标系）。
|
||||
dark_cfg: 暗区配置。
|
||||
prompt_cfg: 剔绿用 HSV 阈值（透传给 ``extract_dark_region_masks``）。
|
||||
|
||||
返回：
|
||||
与帧序列**等长**的逐帧暗部多边形列表（外层帧对齐，内层变长 0~N 个）：
|
||||
每个多边形为 ``(contour_points, 2)`` float32（等弧长重采样）。
|
||||
``enabled=False`` 时每帧返回空列表 ``[]``（**仍保留帧对齐**，N 帧进 N 帧出）。
|
||||
|
||||
关键实现：
|
||||
逐帧 ``extract_dark_region_masks`` 取选中暗块，其并集作为下一帧时间一致性的
|
||||
``previous_mask`` 串起来；对每个暗块用 ``significant_contours``（轮廓面积阈
|
||||
``min_area×0.1``，与参考 P:1437 一致）取显著外轮廓，再
|
||||
``resample_closed_contour`` 重采样到 ``contour_points`` 点，重采样退化（None）则跳过。
|
||||
``max_contours_per_region`` 为 None 时不限。
|
||||
校验帧数与多边形数一致，越界即配置缺陷直接报错（不静默对齐）——流式之后这条校验
|
||||
由 ``zip(..., strict=True)`` 执行，因为帧数要遍历完才知道。
|
||||
"""
|
||||
if not bool(dark_cfg.enabled):
|
||||
# 关掉也要保持帧对齐（N 帧进 N 帧出），所以仍旧把帧走一遍，按帧数产空列表。
|
||||
return [[] for _, _ in zip(analysis_frames, roi_polygons, strict=True)]
|
||||
|
||||
num_points = int(dark_cfg.contour_points)
|
||||
contour_area_floor = float(dark_cfg.min_area) * 0.1
|
||||
max_contours = dark_cfg.max_contours_per_region
|
||||
|
||||
dark_polygons_by_frame: list[list[AnyArray]] = []
|
||||
previous_mask: AnyArray | None = None
|
||||
for frame, polygon in zip(analysis_frames, roi_polygons, strict=True):
|
||||
dark_masks, selected_union = extract_dark_region_masks(
|
||||
frame, polygon, dark_cfg, previous_mask, prompt_cfg
|
||||
)
|
||||
previous_mask = selected_union
|
||||
frame_polygons: list[AnyArray] = []
|
||||
for dark_mask in dark_masks:
|
||||
for contour in _significant_contours(dark_mask, contour_area_floor, max_contours):
|
||||
poly = resample_closed_contour(contour, num_points)
|
||||
if poly is not None:
|
||||
frame_polygons.append(poly)
|
||||
dark_polygons_by_frame.append(frame_polygons)
|
||||
|
||||
return dark_polygons_by_frame
|
||||
|
||||
|
||||
def _significant_contours(
|
||||
mask: AnyArray,
|
||||
min_area: float,
|
||||
max_regions: int | None,
|
||||
) -> list[AnyArray]:
|
||||
"""取 mask 的显著外轮廓（面积降序，``min_area`` 过滤，可截断；参考 P:665-690）。
|
||||
|
||||
参数：
|
||||
mask: 单个暗块的 ``dtype=bool`` mask。
|
||||
min_area: 轮廓面积下限（``cv2.contourArea``）。
|
||||
max_regions: 保留的外轮廓数上限（None 不限）。
|
||||
|
||||
返回：
|
||||
显著轮廓点列表（每个 ``(M, 2)`` float32，OpenCV ``CHAIN_APPROX_NONE`` 原始稠密点）。
|
||||
含外轮廓与其下的子轮廓（洞）——外轮廓阈 ``min_area``、子轮廓阈 ``min_area×0.1``。
|
||||
|
||||
关键实现（严格对照参考 ``significant_contours`` P:665-690）：
|
||||
``RETR_CCOMP`` 取外 / 内两层；外轮廓（``hierarchy[3]<0``）面积 ``>= min_area``，
|
||||
按面积降序、``max_regions`` 截断；对每个入选外轮廓，沿 ``hierarchy`` 的
|
||||
first_child（``[idx][2]``）+ next_sibling（``[child][0]``）遍历其洞，面积
|
||||
``>= min_area×0.1`` 的子轮廓一并追加。``split_dark_components`` 已按
|
||||
``preserve_holes`` 决定 mask 是否含洞，故洞是否出现由 mask 自然决定。
|
||||
"""
|
||||
contours, hierarchy = cv2.findContours(
|
||||
mask.astype(np.uint8), cv2.RETR_CCOMP, cv2.CHAIN_APPROX_NONE
|
||||
)
|
||||
# 没有轮廓时 cv2 交回的 hierarchy 是 None。
|
||||
if not contours:
|
||||
return []
|
||||
hierarchy = np.asarray(hierarchy)[0]
|
||||
external = [
|
||||
(cv2.contourArea(contour), idx, contour)
|
||||
for idx, contour in enumerate(contours)
|
||||
if hierarchy[idx][3] < 0 and cv2.contourArea(contour) >= min_area
|
||||
]
|
||||
external.sort(key=lambda item: item[0], reverse=True)
|
||||
selected_external = external if max_regions is None else external[:max_regions]
|
||||
selected: list[AnyArray] = []
|
||||
for _, idx, contour in selected_external:
|
||||
selected.append(np.asarray(contour).reshape(-1, 2).astype(np.float32))
|
||||
child = int(hierarchy[idx][2])
|
||||
while child >= 0:
|
||||
child_contour = np.asarray(contours[child])
|
||||
if cv2.contourArea(child_contour) >= min_area * 0.1:
|
||||
selected.append(child_contour.reshape(-1, 2).astype(np.float32))
|
||||
child = int(hierarchy[child][0])
|
||||
return selected
|
||||
@@ -0,0 +1,322 @@
|
||||
"""overlay.mp4：把病灶轮廓、暗区、背景圆叠回整幅画面逐帧渲染，供人目视检查。
|
||||
|
||||
各视图局部坐标的多边形加上所属视图的左上角偏移，画到整帧坐标上。颜色（BGR）和线宽照参考组件：
|
||||
造影那一半（source）的 ROI 黄 (0,255,255)、分割那一半（target）的 ROI 绿 (0,255,0)、暗区橙红
|
||||
(255,80,20)；背景圆画青 (255,255,0) 而不照参考组件画紫，理由见 `_BACKGROUND_CIRCLE_CYAN`。
|
||||
ROI 半透明填充 alpha=0.18；线条是黑色外描边加彩色内线。
|
||||
|
||||
`render_overlay` 接收已解码的逐帧整帧图像（可以是迭代器）和逐帧图层，返回 mp4 字节：
|
||||
`cv2.VideoWriter` 不能直写内存，所以写临时文件、读回、在 finally 里删掉。fps 由调用方给，非正报错，
|
||||
不回退 25。
|
||||
|
||||
照搬旧项目 `overlay.py`，去掉产品层；`OverlayConfig`、`OverlayLayers` 的字段默认值去掉，调用方全部
|
||||
显式给出。
|
||||
"""
|
||||
|
||||
import tempfile
|
||||
from collections.abc import Iterable, Sequence
|
||||
from dataclasses import dataclass
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import cv2
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._geometry import AnyArray, ViewRect
|
||||
|
||||
# 各层颜色（BGR，照旧项目 findings §3.10 / 参考组件 P:1469-1544）。
|
||||
_OUTLINE_BLACK = (0, 0, 0)
|
||||
_SOURCE_ROI_YELLOW = (0, 255, 255)
|
||||
_TARGET_ROI_GREEN = (0, 255, 0)
|
||||
_DARK_ORANGE_RED = (255, 80, 20)
|
||||
_BACKGROUND_CIRCLE_CYAN = (255, 255, 0)
|
||||
"""自动挑的那个背景圆画成青色。**不能是紫**——紫色是医生画在同一画面上的 ``ROI_pan``
|
||||
轮廓（邻近正常胰腺实质，也就是病变旁边那块没病的胰腺腺体组织）的颜色，而这个圆是程序
|
||||
自己找的一小块干净背景，两者是不同的东西（旧项目 design/0030 决策四）。参考实现
|
||||
把它画成品红 (255,0,255)，而紫色轮廓已经在算了（``tic.py`` 的 ``CURVE_NAMES`` 四条里就有
|
||||
``roi_pan``），两圈都是紫的话叠加在视频上人眼分不出哪个是参照区、哪个是医生画的——
|
||||
而**看 overlay 是第一步唯一的人工验收手段**，分不出就等于没验。
|
||||
青色和绿框、紫色轮廓在色相上都隔得开。"""
|
||||
# ROI 半透明填充透明度（P:1500 硬编码 0.18）。
|
||||
_FILL_ALPHA = 0.18
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class OverlayConfig:
|
||||
"""各画层开关与填充：画不画造影那一半的 ROI（黄）、分割那一半的 ROI（绿）、分割那一半的暗区、
|
||||
造影那一半的暗区、背景圆（青），以及要不要对分割那一半的 ROI 做 alpha=0.18 半透明填充。"""
|
||||
|
||||
draw_source_roi: bool
|
||||
draw_target_polygon: bool
|
||||
draw_dark_region: bool
|
||||
draw_source_dark_region: bool
|
||||
draw_reference_region: bool
|
||||
fill_polygon: bool
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class OverlayLayers:
|
||||
"""一帧要画的各层几何（多边形都是各自视图的局部坐标）。
|
||||
|
||||
- `target_polygon`、`target_rect`：分割那一半的 ROI 和它在整帧中的矩形；
|
||||
- `dark_polygons`：分割那一半的暗区多边形（可为空列表）；
|
||||
- `source_roi_polygon`、`source_rect`：造影那一半的 ROI 和矩形（不画时可为 None）；
|
||||
- `source_dark_polygons`：造影那一半的暗区多边形（可为空列表）；
|
||||
- `reference_region`：背景圆元数据（`enabled`、`center_xy`、`radius_px`，圆心是造影那一半的
|
||||
局部坐标），没有时为 None。
|
||||
"""
|
||||
|
||||
target_polygon: AnyArray
|
||||
target_rect: ViewRect
|
||||
dark_polygons: list[AnyArray]
|
||||
source_roi_polygon: AnyArray | None
|
||||
source_rect: ViewRect | None
|
||||
source_dark_polygons: list[AnyArray]
|
||||
reference_region: dict[str, Any] | None
|
||||
|
||||
|
||||
def _to_global_points(polygon: AnyArray, rect: ViewRect) -> AnyArray:
|
||||
"""把视图局部坐标多边形抬到整帧坐标，并量化为 cv2 折线点 ``(M, 1, 2)`` int32。
|
||||
|
||||
参数:
|
||||
polygon: 形如 ``(N, 2)`` 的局部坐标顶点。
|
||||
rect: 该多边形所属视图在整帧中的矩形（提供左上角偏移）。
|
||||
|
||||
返回:
|
||||
``cv2.polylines`` / ``cv2.fillPoly`` 直接可用的 ``(N, 1, 2)`` int32 点阵。
|
||||
|
||||
关键实现:
|
||||
加 ``rect.x1`` / ``rect.y1`` 偏移后四舍五入取整（参考 P:1487-1489），先 copy
|
||||
避免就地污染入参多边形。
|
||||
"""
|
||||
global_poly = polygon.copy().astype(np.float64)
|
||||
global_poly[:, 0] += rect.x1
|
||||
global_poly[:, 1] += rect.y1
|
||||
return np.round(global_poly).astype(np.int32).reshape(-1, 1, 2)
|
||||
|
||||
|
||||
def _draw_outlined_polygon(
|
||||
frame: AnyArray,
|
||||
points: AnyArray,
|
||||
color: tuple[int, int, int],
|
||||
outline_thickness: int,
|
||||
color_thickness: int,
|
||||
) -> None:
|
||||
"""就地画"黑色外描边 + 彩色内线"的闭合多边形（参考 P:1490-1491 等的双描边手法）。
|
||||
|
||||
参数:
|
||||
frame: 待绘制的整帧（就地修改）。
|
||||
points: ``(N, 1, 2)`` int32 折线点阵（整帧坐标）。
|
||||
color: 彩色内线 BGR。
|
||||
outline_thickness: 黑色外描边线宽（> 彩色线宽，作描边效果）。
|
||||
color_thickness: 彩色内线线宽。
|
||||
"""
|
||||
cv2.polylines(
|
||||
frame,
|
||||
[points],
|
||||
isClosed=True,
|
||||
color=_OUTLINE_BLACK,
|
||||
thickness=outline_thickness,
|
||||
)
|
||||
cv2.polylines(frame, [points], isClosed=True, color=color, thickness=color_thickness)
|
||||
|
||||
|
||||
def _fill_target_polygon(frame: AnyArray, points: AnyArray) -> AnyArray:
|
||||
"""对 target ROI 内部做 alpha=0.18 绿色半透明填充（参考 P:1497-1500）。
|
||||
|
||||
参数:
|
||||
frame: 当前帧。
|
||||
points: target ROI 整帧坐标折线点阵。
|
||||
|
||||
返回:
|
||||
混合后的新帧（``addWeighted(fill, 0.18, frame, 0.82, 0)``）。
|
||||
|
||||
关键实现:
|
||||
在帧副本上实心填充绿色后与原帧按 0.18 / 0.82 加权混合，得到半透明效果。
|
||||
"""
|
||||
fill = frame.copy()
|
||||
cv2.fillPoly(fill, [points], _TARGET_ROI_GREEN)
|
||||
return cv2.addWeighted(fill, _FILL_ALPHA, frame, 1.0 - _FILL_ALPHA, 0)
|
||||
|
||||
|
||||
def _draw_source_roi_layer(frame: AnyArray, layers: OverlayLayers) -> None:
|
||||
"""画左 / source ROI 轮廓（黄，线宽 5→2，参考 P:1481-1491）。需 source_rect。"""
|
||||
if layers.source_roi_polygon is None or layers.source_rect is None:
|
||||
return
|
||||
points = _to_global_points(layers.source_roi_polygon, layers.source_rect)
|
||||
_draw_outlined_polygon(frame, points, _SOURCE_ROI_YELLOW, 5, 2)
|
||||
|
||||
|
||||
def _draw_dark_layer(
|
||||
frame: AnyArray,
|
||||
polygons: Sequence[AnyArray],
|
||||
rect: ViewRect,
|
||||
outline_thickness: int,
|
||||
color_thickness: int,
|
||||
) -> None:
|
||||
"""画一组暗区轮廓（橙红，参考 P:1504-1526）；逐多边形抬到整帧再双描边。
|
||||
|
||||
参数:
|
||||
frame: 当前帧（就地修改）。
|
||||
polygons: 暗区局部坐标多边形列表。
|
||||
rect: 这些暗区所属视图矩形（target 或 source）。
|
||||
outline_thickness: 黑色外描边线宽（target 6 / source 5，参考 P:1510/1525）。
|
||||
color_thickness: 橙红内线线宽（target 3 / source 2，参考 P:1511/1526）。
|
||||
"""
|
||||
for polygon in polygons:
|
||||
points = _to_global_points(polygon, rect)
|
||||
_draw_outlined_polygon(frame, points, _DARK_ORANGE_RED, outline_thickness, color_thickness)
|
||||
|
||||
|
||||
def _draw_reference_layer(frame: AnyArray, layers: OverlayLayers) -> None:
|
||||
"""画 source 视图那个自动挑的背景圆（青，参考 P:1528-1543 画的是紫）。
|
||||
|
||||
关键实现:
|
||||
仅当 ``reference_region.enabled`` 且 ``source_rect`` / ``center_xy`` /
|
||||
``radius_px`` 齐备时绘制；圆心由 source 局部坐标加 ``source_rect`` 偏移抬到整帧。
|
||||
黑色外描边线宽 6、青色内线线宽 3（参考 P:1542-1543）。
|
||||
"""
|
||||
region = layers.reference_region
|
||||
if region is None or layers.source_rect is None or not region.get("enabled", False):
|
||||
return
|
||||
center = region.get("center_xy")
|
||||
radius = region.get("radius_px")
|
||||
if center is None or radius is None:
|
||||
return
|
||||
center_global = (
|
||||
round(float(center[0]) + layers.source_rect.x1),
|
||||
round(float(center[1]) + layers.source_rect.y1),
|
||||
)
|
||||
radius_int = round(float(radius))
|
||||
cv2.circle(frame, center_global, radius_int, color=_OUTLINE_BLACK, thickness=6)
|
||||
cv2.circle(frame, center_global, radius_int, color=_BACKGROUND_CIRCLE_CYAN, thickness=3)
|
||||
|
||||
|
||||
def draw_polygon_overlay(
|
||||
frame: AnyArray,
|
||||
layers: OverlayLayers,
|
||||
cfg: OverlayConfig,
|
||||
) -> AnyArray:
|
||||
"""在整帧坐标系上叠加 ROI / 暗区 / reference 各层，返回画好的新帧（不改入参）。
|
||||
|
||||
参数:
|
||||
frame: 整帧 BGR 图像 ``(H, W, 3)``。
|
||||
layers: 单帧各层几何（局部坐标多边形 + 视图矩形）。
|
||||
cfg: 各画层开关与填充配置。
|
||||
|
||||
返回:
|
||||
叠加后的整帧副本（尺寸 / dtype 与入参一致）。
|
||||
|
||||
关键实现:
|
||||
先 copy 帧避免污染入参；按参考 P:1469-1544 的层序绘制——source ROI（黄）→
|
||||
target ROI 填充 + 轮廓（绿）→ target 暗区（橙红）→ source 暗区（橙红）→
|
||||
背景圆（青）。各层的局部→整帧坐标变换在 ``_to_global_points`` 内统一完成。
|
||||
填充层（``fill_polygon``）须先于轮廓绘制，否则填充会盖住轮廓。
|
||||
"""
|
||||
out = frame.copy()
|
||||
|
||||
# Phase 1: source ROI 轮廓（黄）。
|
||||
if cfg.draw_source_roi:
|
||||
_draw_source_roi_layer(out, layers)
|
||||
|
||||
# Phase 2: target ROI——先半透明填充，再画绿色轮廓（主产物）。
|
||||
target_points = _to_global_points(layers.target_polygon, layers.target_rect)
|
||||
if cfg.fill_polygon:
|
||||
out = _fill_target_polygon(out, target_points)
|
||||
if cfg.draw_target_polygon:
|
||||
_draw_outlined_polygon(out, target_points, _TARGET_ROI_GREEN, 7, 3)
|
||||
|
||||
# Phase 3: target 暗区轮廓（橙红）。
|
||||
if cfg.draw_dark_region and layers.dark_polygons:
|
||||
_draw_dark_layer(out, layers.dark_polygons, layers.target_rect, 6, 3)
|
||||
|
||||
# Phase 4: source 暗区轮廓（橙红）。
|
||||
if (
|
||||
cfg.draw_source_dark_region
|
||||
and layers.source_dark_polygons
|
||||
and layers.source_rect is not None
|
||||
):
|
||||
_draw_dark_layer(out, layers.source_dark_polygons, layers.source_rect, 5, 2)
|
||||
|
||||
# Phase 5: 背景圆（青）。
|
||||
if cfg.draw_reference_region:
|
||||
_draw_reference_layer(out, layers)
|
||||
|
||||
return out
|
||||
|
||||
|
||||
def render_overlay(
|
||||
frames: Iterable[AnyArray],
|
||||
per_frame_layers: Sequence[OverlayLayers],
|
||||
fps: float,
|
||||
cfg: OverlayConfig,
|
||||
*,
|
||||
frame_size: tuple[int, int],
|
||||
) -> bytes:
|
||||
"""逐帧渲染 overlay 并编码为 mp4 字节（fourcc mp4v、原帧尺寸、调用方 fps）。
|
||||
|
||||
参数:
|
||||
frames: 已解码的逐帧整帧 BGR 图像，非空且尺寸一致。**逐帧消费、不留**，
|
||||
所以可以直接喂一个解码器的迭代器。
|
||||
per_frame_layers: 与 ``frames`` 一一对应的逐帧图层；长度须等于帧数。
|
||||
fps: 视频帧率，**必须 > 0**（调用方已探测）；非正直接报错，不静默回退 25.0。
|
||||
cfg: 各画层开关与填充配置。
|
||||
frame_size: 输出视频的 ``(width, height)``。**由调用方给，不再从首帧读**——
|
||||
编码器要在写第一帧之前就建好，而流式之后「先看一眼首帧再决定尺寸」意味着
|
||||
要么把首帧留住、要么把迭代器消费掉一帧。调用方本来就从探测拿到了这个尺寸。
|
||||
|
||||
返回:
|
||||
mp4 视频字节（mp4v 编码），由调用方写文件。
|
||||
|
||||
异常:
|
||||
ValueError: ``frames`` 为空、``frames`` 与 ``per_frame_layers`` 长度不一致
|
||||
（由 ``zip`` 的 ``strict`` 执行——流式之后帧数要遍历完才知道）、
|
||||
``fps`` 非正，或 ``frame_size`` 非正。
|
||||
RuntimeError: ``cv2.VideoWriter`` 无法打开（编码器不可用 / 尺寸非法）。
|
||||
|
||||
关键实现:
|
||||
``cv2.VideoWriter`` 不能直写内存，故写到临时 ``.mp4`` 文件、读回字节、``finally``
|
||||
``unlink`` 删除（仿 ``roi_inference/frames.py`` 临时文件纪律，确保异常路径也清理）；
|
||||
本模块不落本地、不碰对象存储。逐帧 ``draw_polygon_overlay`` 保持帧尺寸不变，
|
||||
故无需逐帧重校尺寸。
|
||||
"""
|
||||
# Phase 1: 防御校验（不静默兜底，缺陷直接报错）。
|
||||
if float(fps) <= 0.0:
|
||||
raise ValueError(f"fps 必须为正（调用方应已探测），得到 {fps}。")
|
||||
width, height = frame_size
|
||||
if width <= 0 or height <= 0:
|
||||
raise ValueError(f"frame_size 必须为正，得到 {frame_size}。")
|
||||
if not per_frame_layers:
|
||||
raise ValueError("per_frame_layers 不能为空：一帧都不写产不出视频。")
|
||||
|
||||
# Phase 2: 写临时 mp4，结束（含异常）后必删。
|
||||
with tempfile.NamedTemporaryFile(suffix=".mp4", delete=False) as handle:
|
||||
temp_path = Path(handle.name)
|
||||
try:
|
||||
writer = cv2.VideoWriter(
|
||||
str(temp_path),
|
||||
cv2.VideoWriter.fourcc(*"mp4v"),
|
||||
float(fps),
|
||||
(width, height),
|
||||
)
|
||||
if not writer.isOpened():
|
||||
raise RuntimeError(f"cv2.VideoWriter 无法打开（尺寸 {width}x{height}、fps {fps}）。")
|
||||
try:
|
||||
# strict：帧数要遍历完才知道，长度不一致在这里当场抛，不静默按较短的截断。
|
||||
for index, (frame, layers) in enumerate(zip(frames, per_frame_layers, strict=True)):
|
||||
if index == 0 and frame.shape[:2] != (height, width):
|
||||
# **尺寸不符的帧会被 VideoWriter 静默丢掉**：实测写 6 帧读回 3 帧,
|
||||
# 只在 stderr 留一行警告，不抛异常也不改返回值。结果是 overlay 视频
|
||||
# 短一截甚至是空的，而它照样有字节、照样被写成文件。
|
||||
# 尺寸以前由首帧决定所以结构上不可能不符，现在由调用方给，这条得自己守。
|
||||
raise ValueError(
|
||||
f"帧尺寸 {frame.shape[:2]} 与 frame_size {(width, height)} 不符；"
|
||||
"尺寸不符的帧会被编码器静默丢掉"
|
||||
)
|
||||
writer.write(draw_polygon_overlay(frame, layers, cfg))
|
||||
finally:
|
||||
writer.release()
|
||||
|
||||
return temp_path.read_bytes()
|
||||
finally:
|
||||
temp_path.unlink(missing_ok=True)
|
||||
@@ -0,0 +1,151 @@
|
||||
"""第 9 步：从过了质量门的 Lognormal 拟合曲线算 11 个灌注指标（口径仿照商用造影定量软件 VueBOX）。
|
||||
|
||||
指标在 MP4 显示灰度域上算（不是线性化的回声功率），时间以视频开头为零点，面积用「拟合曲线减
|
||||
基线」。依赖观测范围的指标（如洗出段面积）在观测不够时为 None，不外推。
|
||||
|
||||
- PE：峰值增强（峰值减基线）；TTP：达峰时间（D3）；WiR：最大洗入斜率（D2 处导数）；
|
||||
- TI：洗入切线与基线的交点；RT：上升时间 TTP − TI；WiAUC：TI 到 TTP 的面积；WiPI：WiAUC / RT；
|
||||
- WoR：最大洗出斜率（D4 处导数）；TO：洗出切线与基线的交点；FT：下降时间 TO − TTP；
|
||||
- WoAUC：TTP 到 TO 的面积；WiWoAUC：两段面积之和；mTT：TI 到 TO 之间的平均渡越时刻；mTTl：mTT − TI。
|
||||
|
||||
照搬旧项目 `perfusion_metrics.py`，去掉产品层。
|
||||
"""
|
||||
|
||||
import math
|
||||
from dataclasses import dataclass
|
||||
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._tic_fitting import (
|
||||
LognormalFitResult,
|
||||
LognormalParameters,
|
||||
evaluate_lognormal,
|
||||
evaluate_lognormal_derivative,
|
||||
)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class LognormalPerfusionMetrics:
|
||||
"""11 项指标及三个支撑时刻（TI、mTT、TO），单位见字段名后缀；观测不够的为 None。"""
|
||||
|
||||
pe_au: float
|
||||
wi_auc_au_s: float | None
|
||||
rt_s: float | None
|
||||
mttl_s: float | None
|
||||
ttp_s: float
|
||||
wir_au_per_s: float
|
||||
wipi_au: float | None
|
||||
wo_auc_au_s: float | None
|
||||
wiwo_auc_au_s: float | None
|
||||
ft_s: float | None
|
||||
wor_au_per_s: float | None
|
||||
ti_s: float | None
|
||||
mtt_s: float | None
|
||||
to_s: float | None
|
||||
|
||||
|
||||
def _normal_cdf(value: float) -> float:
|
||||
return 0.5 * (1.0 + math.erf(value / math.sqrt(2.0)))
|
||||
|
||||
|
||||
def _elapsed_cdf(time_s: float, parameters: LognormalParameters) -> float:
|
||||
elapsed_s = time_s - parameters.arrival_offset_s
|
||||
if elapsed_s <= 0.0:
|
||||
return 0.0
|
||||
return _normal_cdf((math.log(elapsed_s) - parameters.mu) / parameters.sigma)
|
||||
|
||||
|
||||
def _elapsed_first_moment_cdf(time_s: float, parameters: LognormalParameters) -> float:
|
||||
elapsed_s = time_s - parameters.arrival_offset_s
|
||||
if elapsed_s <= 0.0:
|
||||
return 0.0
|
||||
standardized = (math.log(elapsed_s) - parameters.mu - parameters.sigma**2) / parameters.sigma
|
||||
return math.exp(parameters.mu + 0.5 * parameters.sigma**2) * _normal_cdf(standardized)
|
||||
|
||||
|
||||
def _pulse_auc(start_s: float, end_s: float, parameters: LognormalParameters) -> float:
|
||||
if end_s <= start_s:
|
||||
raise ValueError("Lognormal 面积终点必须晚于起点")
|
||||
return parameters.area * (_elapsed_cdf(end_s, parameters) - _elapsed_cdf(start_s, parameters))
|
||||
|
||||
|
||||
def _pulse_first_moment(start_s: float, end_s: float, parameters: LognormalParameters) -> float:
|
||||
probability = _elapsed_cdf(end_s, parameters) - _elapsed_cdf(start_s, parameters)
|
||||
elapsed_moment = _elapsed_first_moment_cdf(end_s, parameters) - _elapsed_first_moment_cdf(
|
||||
start_s, parameters
|
||||
)
|
||||
return parameters.area * (parameters.arrival_offset_s * probability + elapsed_moment)
|
||||
|
||||
|
||||
def _observed_time(value: float, lower: float, upper: float) -> float | None:
|
||||
return value if lower <= value <= upper else None
|
||||
|
||||
|
||||
def compute_lognormal_perfusion_metrics(
|
||||
fit: LognormalFitResult,
|
||||
) -> LognormalPerfusionMetrics | None:
|
||||
"""拟合没过质量门时返回 None；过了就算 11 项指标，观测范围不够的那几项为 None。"""
|
||||
if not fit.success:
|
||||
return None
|
||||
parameters = fit.parameters
|
||||
if parameters is None or fit.d2_time_s is None or fit.d3_time_s is None:
|
||||
raise ValueError("成功 Lognormal 拟合缺少参数或洗入/峰值时刻")
|
||||
|
||||
d2_time_s = float(fit.d2_time_s)
|
||||
ttp_s = float(fit.d3_time_s)
|
||||
peak_value = float(evaluate_lognormal(np.asarray([ttp_s]), parameters)[0])
|
||||
pe_au = peak_value - parameters.baseline
|
||||
wir_au_per_s = evaluate_lognormal_derivative(d2_time_s, parameters)
|
||||
if not math.isfinite(pe_au) or pe_au <= 0.0:
|
||||
raise ValueError("成功 Lognormal 拟合的峰值增强必须为有限正数")
|
||||
if not math.isfinite(wir_au_per_s) or wir_au_per_s <= 0.0:
|
||||
raise ValueError("成功 Lognormal 拟合的最大洗入斜率必须为有限正数")
|
||||
|
||||
d2_value = float(evaluate_lognormal(np.asarray([d2_time_s]), parameters)[0])
|
||||
ti_candidate_s = d2_time_s - (d2_value - parameters.baseline) / wir_au_per_s
|
||||
ti_s = _observed_time(ti_candidate_s, fit.fit_start_s, d2_time_s)
|
||||
rt_s = None if ti_s is None else ttp_s - ti_s
|
||||
wi_auc_au_s = None if ti_s is None else _pulse_auc(ti_s, ttp_s, parameters)
|
||||
wipi_au = None if wi_auc_au_s is None or rt_s is None else wi_auc_au_s / rt_s
|
||||
|
||||
wor_au_per_s: float | None = None
|
||||
to_s: float | None = None
|
||||
ft_s: float | None = None
|
||||
wo_auc_au_s: float | None = None
|
||||
if fit.d4_time_s is not None:
|
||||
d4_time_s = float(fit.d4_time_s)
|
||||
wor_au_per_s = evaluate_lognormal_derivative(d4_time_s, parameters)
|
||||
if not math.isfinite(wor_au_per_s) or wor_au_per_s >= 0.0:
|
||||
raise ValueError("成功 Lognormal 拟合的最大洗出斜率必须为有限负数")
|
||||
d4_value = float(evaluate_lognormal(np.asarray([d4_time_s]), parameters)[0])
|
||||
to_candidate_s = d4_time_s - (d4_value - parameters.baseline) / wor_au_per_s
|
||||
to_s = _observed_time(to_candidate_s, d4_time_s, fit.fit_end_s)
|
||||
if to_s is not None:
|
||||
ft_s = to_s - ttp_s
|
||||
wo_auc_au_s = _pulse_auc(ttp_s, to_s, parameters)
|
||||
|
||||
wiwo_auc_au_s = (
|
||||
None if wi_auc_au_s is None or wo_auc_au_s is None else wi_auc_au_s + wo_auc_au_s
|
||||
)
|
||||
mtt_s: float | None = None
|
||||
mttl_s: float | None = None
|
||||
if ti_s is not None and to_s is not None and wiwo_auc_au_s is not None:
|
||||
mtt_s = _pulse_first_moment(ti_s, to_s, parameters) / wiwo_auc_au_s
|
||||
mttl_s = mtt_s - ti_s
|
||||
|
||||
return LognormalPerfusionMetrics(
|
||||
pe_au=pe_au,
|
||||
wi_auc_au_s=wi_auc_au_s,
|
||||
rt_s=rt_s,
|
||||
mttl_s=mttl_s,
|
||||
ttp_s=ttp_s,
|
||||
wir_au_per_s=wir_au_per_s,
|
||||
wipi_au=wipi_au,
|
||||
wo_auc_au_s=wo_auc_au_s,
|
||||
wiwo_auc_au_s=wiwo_auc_au_s,
|
||||
ft_s=ft_s,
|
||||
wor_au_per_s=wor_au_per_s,
|
||||
ti_s=ti_s,
|
||||
mtt_s=mtt_s,
|
||||
to_s=to_s,
|
||||
)
|
||||
@@ -0,0 +1,480 @@
|
||||
"""ceus 对外的两段接口，中间由 `scripts/ceus_quant.py` 调 segment 做第 3 步分割。
|
||||
|
||||
- `prepare_segmentation`（第 1–2 步）：读视频、找分屏线、判断造影画面在哪一半、选提示帧，交出
|
||||
要交给分割的那一半（灰阶那一半，`target_view`）的逐帧图像和提示帧掩膜（`target_view` 局部坐标）。
|
||||
- `quantify_case`（第 4–9 步）：拿分割得到的逐帧掩膜（`target_view` 局部坐标）整理成平滑轮廓并与
|
||||
绿框融合、挪回造影那一半，找暗区、选背景圆、逐帧采 TIC、拟合、定位 7 个关键帧、算灌注指标。
|
||||
|
||||
ceus 不引用 segment（`结构.md` 第三部分），所以分割夹在两段中间。编排顺序照旧 debug 脚本
|
||||
`t1_local_debug.py`（第 1–6 步）和 `tic_lognormal_replay.py`（第 6–9 步）：旧 replay 读回轨迹后
|
||||
重新从视频采了一遍 TIC，两遍输入相同，这里只采一遍。
|
||||
|
||||
另有给第三阶段脚本写产物用的几个函数：`tic_plot_markers`（tic.png 上 D1–D4 的竖线）、
|
||||
`render_case_overlay`（overlay.mp4 的字节）、`keyframe_images`（定位到的关键帧整幅原图）、
|
||||
`quant_summary`（summary.json 里拟合、指标、关键帧那几段）。
|
||||
"""
|
||||
|
||||
import math
|
||||
from collections.abc import Iterator
|
||||
from dataclasses import dataclass
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._anchor_frames import (
|
||||
ANCHOR_SLOTS,
|
||||
STATUS_OK,
|
||||
AnchorFrameLocation,
|
||||
locate_anchor_frames,
|
||||
)
|
||||
from cadx.ceus._config import CeusConfig
|
||||
from cadx.ceus._contrast_side import ContrastSide, locate_contrast_side
|
||||
from cadx.ceus._dark_region import extract_dark_region_polygons
|
||||
from cadx.ceus._geometry import AnyArray, Frame, Mask, crop_frame
|
||||
from cadx.ceus._mapping import map_source_mask_to_target, target_polygons_to_source
|
||||
from cadx.ceus._overlay import OverlayConfig, OverlayLayers, render_overlay
|
||||
from cadx.ceus._perfusion_metrics import (
|
||||
LognormalPerfusionMetrics,
|
||||
compute_lognormal_perfusion_metrics,
|
||||
)
|
||||
from cadx.ceus._polygons import (
|
||||
apply_prompt_frame_position_constraints,
|
||||
fuse_prompt_prior_with_medsam,
|
||||
masks_to_smoothed_polygons,
|
||||
)
|
||||
from cadx.ceus._reference_region import select_reference_region
|
||||
from cadx.ceus._split_screen import split_views
|
||||
from cadx.ceus._tic import CURVE_NAMES, TicResult, compute_tic_for_video
|
||||
from cadx.ceus._tic_curves import TicCurves, from_tic_result
|
||||
from cadx.ceus._tic_fitting import LognormalFitResult, evaluate_lognormal
|
||||
from cadx.ceus._tic_plot import TicPlotMarker
|
||||
from cadx.ceus._trajectory import Trajectory
|
||||
from cadx.ceus._video import VideoFile
|
||||
|
||||
# 灌注指标只对病灶两条曲线算（同旧 replay 的 `_LESION_CURVES`）。
|
||||
LESION_CURVES = ("roi_all", "roi_perf")
|
||||
|
||||
# tic.png 上 D1–D4 竖线的颜色（BGR，同旧 replay）。
|
||||
_ANCHOR_COLORS = {
|
||||
"D1": (30, 30, 220),
|
||||
"D2": (40, 150, 40),
|
||||
"D3": (180, 60, 160),
|
||||
"D4": (180, 130, 30),
|
||||
}
|
||||
|
||||
# overlay.mp4 的画层（同旧 debug 脚本 `render_overlay_bytes`）：两半都画 ROI，暗区只画在造影那一半。
|
||||
_OVERLAY_CONFIG = OverlayConfig(
|
||||
draw_source_roi=True,
|
||||
draw_target_polygon=True,
|
||||
draw_dark_region=False,
|
||||
draw_source_dark_region=True,
|
||||
draw_reference_region=True,
|
||||
fill_polygon=False,
|
||||
)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class SegmentationInput:
|
||||
"""第一段（第 1–2 步）的结果。
|
||||
|
||||
- `video_path`、`fps`、`frame_count`、`frame_size`（宽, 高）、`frame_times_s`：视频本身；
|
||||
- `split_x`：分屏线的列；`contrast`：哪一半是造影画面、两个视图矩形（整帧坐标）、提示帧
|
||||
（造影那一半局部坐标的绿框掩膜）；
|
||||
- `target_frames`：要交给分割的那一半（`contrast.target_view`，灰阶画面）的逐帧 BGR 图像，
|
||||
与帧号一一对应；
|
||||
- `prompt_masks`：提示帧帧号 → 这一帧绿框挪到 `target_view` 局部坐标后的布尔掩膜；
|
||||
- `first_contrast_frame`：造影那一半的第一帧（选背景圆用）。
|
||||
"""
|
||||
|
||||
video_path: Path
|
||||
fps: float
|
||||
frame_count: int
|
||||
frame_size: tuple[int, int]
|
||||
frame_times_s: AnyArray
|
||||
split_x: int
|
||||
contrast: ContrastSide
|
||||
target_frames: list[Frame]
|
||||
prompt_masks: dict[int, Mask]
|
||||
first_contrast_frame: Frame
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class QuantResult:
|
||||
"""第二段（第 4–9 步）的结果。
|
||||
|
||||
- `trajectory`：逐帧轮廓、暗区、背景圆和几何（写 trajectory.npz 用 `pack_trajectory`）；
|
||||
- `tic`：逐帧采样和四条曲线的基线、拟合；`curves`：同一份数据的 npz 载体（`pack_tic_curves`）；
|
||||
- `injection_offset_s`：推注时间 X（秒，视频第 0 秒时已推注的时长）；
|
||||
- `lesion_fit_ready`：病灶两条曲线的拟合是否都过了质量门；
|
||||
- `anchors`：7 个关键帧槽位的定位（按 `ANCHOR_SLOTS` 顺序）；
|
||||
- `perfusion`：`roi_all`、`roi_perf` 各自的灌注指标，拟合没过质量门为 None；
|
||||
- `pan_detected_fraction`：紫圈检出率（非重复帧里认出紫圈的帧占比）。
|
||||
"""
|
||||
|
||||
trajectory: Trajectory
|
||||
tic: TicResult
|
||||
curves: TicCurves
|
||||
injection_offset_s: float
|
||||
lesion_fit_ready: bool
|
||||
anchors: dict[str, AnchorFrameLocation]
|
||||
perfusion: dict[str, LognormalPerfusionMetrics | None]
|
||||
pan_detected_fraction: float
|
||||
|
||||
|
||||
def prepare_segmentation(video_path: Path, cfg: CeusConfig) -> SegmentationInput:
|
||||
"""第 1–2 步：解码、分屏、找绿框定造影那一半、选提示帧，交出分割的输入。
|
||||
|
||||
找不到可信分屏线抛 `SplitDetectionError`（第 1 步），两半都有或都没有绿框抛 `ContrastSideError`
|
||||
（第 2 步），两者都是 `CeusCaseError`。
|
||||
"""
|
||||
video = VideoFile(video_path)
|
||||
fps, frame_count, width, height = video.probe()
|
||||
frame_times_s = video.frame_times_s()
|
||||
full_frames = [frame for _, frame in video.iter_frames()]
|
||||
if len(full_frames) != frame_count:
|
||||
raise RuntimeError(
|
||||
f"解码帧数 {len(full_frames)} 与容器记录的帧数 {frame_count} 不一致：{video_path}"
|
||||
)
|
||||
|
||||
split_x, left_view, right_view = split_views(full_frames, cfg.views)
|
||||
contrast = locate_contrast_side(
|
||||
full_frames, left_view, right_view, cfg.prompt_source, cfg.keyframes
|
||||
)
|
||||
target_frames = [crop_frame(frame, contrast.target_view).copy() for frame in full_frames]
|
||||
first_contrast_frame = crop_frame(full_frames[0], contrast.contrast_view).copy()
|
||||
prompt_masks = {
|
||||
prompt.frame_index: map_source_mask_to_target(
|
||||
prompt.mask, contrast.contrast_view, contrast.target_view
|
||||
)
|
||||
for prompt in contrast.prompt_frames
|
||||
}
|
||||
return SegmentationInput(
|
||||
video_path=video_path,
|
||||
fps=fps,
|
||||
frame_count=frame_count,
|
||||
frame_size=(width, height),
|
||||
frame_times_s=frame_times_s,
|
||||
split_x=split_x,
|
||||
contrast=contrast,
|
||||
target_frames=target_frames,
|
||||
prompt_masks=prompt_masks,
|
||||
first_contrast_frame=first_contrast_frame,
|
||||
)
|
||||
|
||||
|
||||
def _contrast_frames(prepared: SegmentationInput) -> Iterator[Frame]:
|
||||
"""重新解码，逐帧交出造影那一半（不把整段留在内存里）。"""
|
||||
for _, frame in VideoFile(prepared.video_path).iter_frames():
|
||||
yield crop_frame(frame, prepared.contrast.contrast_view)
|
||||
|
||||
|
||||
def _validate_masks(prepared: SegmentationInput, target_masks: AnyArray) -> None:
|
||||
"""分割结果必须是 `(帧数, target 高, target 宽)` 的布尔数组。"""
|
||||
view = prepared.contrast.target_view
|
||||
expected = (prepared.frame_count, view.height, view.width)
|
||||
if target_masks.dtype != np.bool_ or target_masks.shape != expected:
|
||||
raise ValueError(
|
||||
f"分割掩膜应为形状 {expected} 的布尔数组，"
|
||||
f"收到 {target_masks.shape} {target_masks.dtype}"
|
||||
)
|
||||
|
||||
|
||||
def quantify_case(
|
||||
prepared: SegmentationInput,
|
||||
target_masks: AnyArray,
|
||||
injection_offset_s: float,
|
||||
cfg: CeusConfig,
|
||||
) -> QuantResult:
|
||||
"""第 4–9 步：轮廓、暗区、背景圆、TIC、拟合、关键帧、灌注指标。
|
||||
|
||||
`target_masks` 是分割在 `target_view` 局部坐标上的逐帧掩膜 `(N, H, W)` 布尔数组；
|
||||
`injection_offset_s` 是超声报告表的「推注开始时间s」X。分割掩膜一帧有效轮廓都没有时抛
|
||||
`NoValidContourError`（第 4 步，`CeusCaseError`）。拟合没过质量门不算失败：D1–D4 未定位、
|
||||
灌注指标为 None。
|
||||
"""
|
||||
_validate_masks(prepared, target_masks)
|
||||
contrast = prepared.contrast
|
||||
frame_count = prepared.frame_count
|
||||
prompt_frame_dicts: list[dict[str, Any]] = [
|
||||
{"frame_index": prompt.frame_index, "mask": prompt.mask}
|
||||
for prompt in contrast.prompt_frames
|
||||
]
|
||||
|
||||
# 第 4 步：掩膜 → 平滑多边形 → 位置约束 → 与绿框融合 → 挪回造影那一半。
|
||||
masks_by_frame = {index: target_masks[index] for index in range(frame_count)}
|
||||
target_polygons = masks_to_smoothed_polygons(masks_by_frame, frame_count, cfg.smoothing)
|
||||
target_polygons = apply_prompt_frame_position_constraints(
|
||||
target_polygons,
|
||||
prompt_frame_dicts,
|
||||
contrast.contrast_view,
|
||||
contrast.target_view,
|
||||
cfg.constraints.keyframe_position,
|
||||
)
|
||||
target_polygons = fuse_prompt_prior_with_medsam(
|
||||
target_polygons,
|
||||
prompt_frame_dicts,
|
||||
contrast.contrast_view,
|
||||
contrast.target_view,
|
||||
cfg.roi_shape,
|
||||
cfg.smoothing,
|
||||
)
|
||||
source_polygons = target_polygons_to_source(
|
||||
target_polygons, contrast.contrast_view, contrast.target_view
|
||||
)
|
||||
|
||||
# 第 5 步：暗区（在造影那一半上找）。
|
||||
dark_polygons = extract_dark_region_polygons(
|
||||
_contrast_frames(prepared),
|
||||
source_polygons,
|
||||
cfg.interior_dark_region,
|
||||
cfg.prompt_source,
|
||||
)
|
||||
|
||||
# 第 6 步：背景圆、四条 TIC；第 7 步拟合在 compute_tic_for_video 里。
|
||||
source_polygon_array = np.asarray(source_polygons, dtype=np.float32)
|
||||
reference_mask, reference_meta = select_reference_region(
|
||||
prepared.first_contrast_frame,
|
||||
source_polygon_array,
|
||||
cfg.prompt_source,
|
||||
cfg.pan_source,
|
||||
cfg.tic_analysis,
|
||||
)
|
||||
tic = compute_tic_for_video(
|
||||
_contrast_frames(prepared),
|
||||
source_polygons,
|
||||
dark_polygons,
|
||||
reference_mask,
|
||||
prepared.fps,
|
||||
cfg.prompt_source,
|
||||
cfg.pan_source,
|
||||
cfg.tic_analysis,
|
||||
prepared.frame_times_s,
|
||||
)
|
||||
|
||||
# 第 8 步：关键帧。
|
||||
lesion_fit_ready = all(tic.fits[name].success for name in LESION_CURVES)
|
||||
frame_indices = np.arange(frame_count, dtype=np.int32)
|
||||
anchors = locate_anchor_frames(
|
||||
perf_fit=tic.fits["roi_perf"],
|
||||
perf_baseline_end_index=tic.baselines["roi_perf"].frame_end,
|
||||
lesion_fit_ready=lesion_fit_ready,
|
||||
frame_indices=frame_indices,
|
||||
time_video_s=prepared.frame_times_s,
|
||||
injection_offset_s=injection_offset_s,
|
||||
cfg=cfg.tic_analysis.anchor_frames,
|
||||
)
|
||||
|
||||
# 第 9 步：灌注指标。
|
||||
perfusion = {
|
||||
name: compute_lognormal_perfusion_metrics(tic.fits[name]) for name in LESION_CURVES
|
||||
}
|
||||
|
||||
trajectory = Trajectory(
|
||||
frame_indices=frame_indices,
|
||||
fps=float(prepared.fps),
|
||||
target_polygons=np.asarray(target_polygons, dtype=np.float32),
|
||||
source_polygons=source_polygon_array,
|
||||
source_view_rect=contrast.contrast_view,
|
||||
target_view_rect=contrast.target_view,
|
||||
source_effective_rect=contrast.contrast_view,
|
||||
target_effective_rect=contrast.target_view,
|
||||
source_polygon_coordinate_space="source_view_local",
|
||||
target_polygon_coordinate_space="target_effective_view_local",
|
||||
reference_mask=reference_mask.astype(np.uint8),
|
||||
reference_region=reference_meta,
|
||||
keyframe_indices=np.asarray(
|
||||
[prompt.frame_index for prompt in contrast.prompt_frames], dtype=np.int32
|
||||
),
|
||||
dark_polygons=dark_polygons,
|
||||
calibration={
|
||||
"split_x": prepared.split_x,
|
||||
"source": "auto_detected",
|
||||
"contrast_side": contrast.side,
|
||||
},
|
||||
)
|
||||
curves = from_tic_result(
|
||||
tic,
|
||||
fps=prepared.fps,
|
||||
injection_offset_s=injection_offset_s,
|
||||
config_version=cfg.config_version,
|
||||
)
|
||||
return QuantResult(
|
||||
trajectory=trajectory,
|
||||
tic=tic,
|
||||
curves=curves,
|
||||
injection_offset_s=injection_offset_s,
|
||||
lesion_fit_ready=lesion_fit_ready,
|
||||
anchors=anchors,
|
||||
perfusion=perfusion,
|
||||
pan_detected_fraction=tic.pan_detected_fraction,
|
||||
)
|
||||
|
||||
|
||||
def tic_plot_markers(anchors: dict[str, AnchorFrameLocation]) -> tuple[TicPlotMarker, ...]:
|
||||
"""tic.png 上 D1–D4 的竖线：有锚点时刻的画在锚点时刻，没定位到帧的在标签后注明状态。"""
|
||||
markers: list[TicPlotMarker] = []
|
||||
for slot in ("D1", "D2", "D3", "D4"):
|
||||
anchor = anchors[slot]
|
||||
if anchor.anchor_time_video_s is None:
|
||||
continue
|
||||
status = "" if anchor.status == STATUS_OK else f" [{anchor.status}]"
|
||||
markers.append(
|
||||
TicPlotMarker(f"{slot} model{status}", anchor.anchor_time_video_s, _ANCHOR_COLORS[slot])
|
||||
)
|
||||
return tuple(markers)
|
||||
|
||||
|
||||
def render_case_overlay(prepared: SegmentationInput, result: QuantResult) -> bytes:
|
||||
"""overlay.mp4 的字节：重新解码整段视频，两半都画 ROI，造影那一半画暗区和背景圆。"""
|
||||
trajectory = result.trajectory
|
||||
layers = [
|
||||
OverlayLayers(
|
||||
target_polygon=trajectory.target_polygons[index],
|
||||
target_rect=trajectory.target_view_rect,
|
||||
dark_polygons=[],
|
||||
source_roi_polygon=trajectory.source_polygons[index],
|
||||
source_rect=trajectory.source_view_rect,
|
||||
source_dark_polygons=list(trajectory.dark_polygons[index]),
|
||||
reference_region=trajectory.reference_region,
|
||||
)
|
||||
for index in range(prepared.frame_count)
|
||||
]
|
||||
frames = (frame for _, frame in VideoFile(prepared.video_path).iter_frames())
|
||||
return render_overlay(
|
||||
frames, layers, prepared.fps, _OVERLAY_CONFIG, frame_size=prepared.frame_size
|
||||
)
|
||||
|
||||
|
||||
def keyframe_images(
|
||||
prepared: SegmentationInput, anchors: dict[str, AnchorFrameLocation]
|
||||
) -> dict[str, Frame]:
|
||||
"""定位到的关键帧的整幅原始画面（左右两半都在，不裁剪），槽位 → BGR 图像；未定位的槽位没有。"""
|
||||
wanted: dict[int, list[str]] = {}
|
||||
for slot, anchor in anchors.items():
|
||||
if anchor.frame_index is not None:
|
||||
wanted.setdefault(anchor.frame_index, []).append(slot)
|
||||
images: dict[str, Frame] = {}
|
||||
for index, frame in VideoFile(prepared.video_path).iter_frames():
|
||||
for slot in wanted.get(index, []):
|
||||
images[slot] = frame.copy()
|
||||
missing = sorted(slot for slots in wanted.values() for slot in slots if slot not in images)
|
||||
if missing:
|
||||
raise RuntimeError(f"重新解码时没读到关键帧 {missing}：{prepared.video_path}")
|
||||
return {slot: images[slot] for slot in ANCHOR_SLOTS if slot in images}
|
||||
|
||||
|
||||
def _finite_or_none(value: float | None) -> float | None:
|
||||
"""非有限数（NaN、inf）在 JSON 里写 null。"""
|
||||
if value is None or not math.isfinite(float(value)):
|
||||
return None
|
||||
return float(value)
|
||||
|
||||
|
||||
def _fit_summary(fit: LognormalFitResult) -> dict[str, Any]:
|
||||
"""一条曲线的拟合摘要（同旧 replay 的 `_fit_summary`）。"""
|
||||
parameters = fit.parameters
|
||||
anchor_times = {
|
||||
"D2": _finite_or_none(fit.d2_time_s),
|
||||
"D3": _finite_or_none(fit.d3_time_s),
|
||||
"D4": _finite_or_none(fit.d4_time_s),
|
||||
}
|
||||
anchor_values = {
|
||||
slot: None
|
||||
if parameters is None or time_s is None
|
||||
else float(evaluate_lognormal(np.asarray([time_s]), parameters)[0])
|
||||
for slot, time_s in anchor_times.items()
|
||||
}
|
||||
return {
|
||||
"success": fit.success,
|
||||
"failure_reason": fit.failure_reason or None,
|
||||
"curve_model": fit.curve_model,
|
||||
"baseline_source": fit.baseline_source,
|
||||
"valid_samples": int(fit.valid_sample_mask.sum()),
|
||||
"duplicate_samples": int(fit.duplicate_sample_mask.sum()),
|
||||
"parameters": None
|
||||
if parameters is None
|
||||
else {
|
||||
"area": parameters.area,
|
||||
"mu": parameters.mu,
|
||||
"sigma": parameters.sigma,
|
||||
"arrival_offset_s": parameters.arrival_offset_s,
|
||||
"baseline": parameters.baseline,
|
||||
},
|
||||
"quality": {
|
||||
"r2": _finite_or_none(fit.r2),
|
||||
"nrmse": _finite_or_none(fit.nrmse),
|
||||
"residual_noise_ratio": _finite_or_none(fit.residual_noise_ratio),
|
||||
"multi_peak_flag": fit.multi_peak_flag,
|
||||
},
|
||||
"model_arrival_time_video_s": _finite_or_none(fit.d1_time_s),
|
||||
"model_landmark_times_video_s": anchor_times,
|
||||
"model_landmark_intensities": anchor_values,
|
||||
}
|
||||
|
||||
|
||||
def _perfusion_summary(metrics: LognormalPerfusionMetrics | None) -> dict[str, Any] | None:
|
||||
"""灌注指标摘要（同旧 replay 的 `_perfusion_metrics_summary`），拟合没过质量门为 None。"""
|
||||
if metrics is None:
|
||||
return None
|
||||
return {
|
||||
"signal_domain": "db_display",
|
||||
"time_origin": "video_start",
|
||||
"area_curve": "fitted_minus_baseline",
|
||||
"values": {
|
||||
"PE_au": metrics.pe_au,
|
||||
"WiAUC_au_s": metrics.wi_auc_au_s,
|
||||
"RT_s": metrics.rt_s,
|
||||
"mTTl_s": metrics.mttl_s,
|
||||
"TTP_s": metrics.ttp_s,
|
||||
"WiR_au_per_s": metrics.wir_au_per_s,
|
||||
"WiPI_au": metrics.wipi_au,
|
||||
"WoAUC_au_s": metrics.wo_auc_au_s,
|
||||
"WiWoAUC_au_s": metrics.wiwo_auc_au_s,
|
||||
"FT_s": metrics.ft_s,
|
||||
"WoR_au_per_s": metrics.wor_au_per_s,
|
||||
},
|
||||
"supporting_times_s": {
|
||||
"TI": metrics.ti_s,
|
||||
"mTT": metrics.mtt_s,
|
||||
"TO": metrics.to_s,
|
||||
},
|
||||
}
|
||||
|
||||
|
||||
def _anchor_summary(anchor: AnchorFrameLocation) -> dict[str, Any]:
|
||||
"""一个关键帧槽位的定位：帧号、视频时间、推注后时刻、状态、未定位原因。"""
|
||||
return {
|
||||
"frame_index": anchor.frame_index,
|
||||
"frame_time_video_s": anchor.frame_time_video_s,
|
||||
"time_since_injection_s": anchor.time_s,
|
||||
"status": anchor.status,
|
||||
"reason": anchor.reason,
|
||||
"anchor_time_video_s": anchor.anchor_time_video_s,
|
||||
"time_error_s": anchor.time_error_s,
|
||||
"source_curve": anchor.source_curve,
|
||||
}
|
||||
|
||||
|
||||
def quant_summary(result: QuantResult) -> dict[str, Any]:
|
||||
"""一例的定量摘要（可直接写成 JSON）：拟合、灌注指标、紫圈检出率、7 个关键帧的定位。
|
||||
|
||||
不含视频号、路径等病例标识，由脚本按需另加。
|
||||
"""
|
||||
tic = result.tic
|
||||
dynamic_slots = ("B0", "D1", "D2", "D3", "D4")
|
||||
return {
|
||||
"injection_offset_s": result.injection_offset_s,
|
||||
"reference_valid": tic.reference_valid,
|
||||
"pan_detected_fraction": result.pan_detected_fraction,
|
||||
"early_truncation_flag": tic.early_truncation_flag,
|
||||
"lesion_fit_ready": result.lesion_fit_ready,
|
||||
"fits": {name: _fit_summary(tic.fits[name]) for name in CURVE_NAMES},
|
||||
"perfusion_metrics": {
|
||||
name: _perfusion_summary(result.perfusion[name]) for name in LESION_CURVES
|
||||
},
|
||||
"anchor_frames": {slot: _anchor_summary(item) for slot, item in result.anchors.items()},
|
||||
"available_dynamic_slots": [
|
||||
slot for slot in dynamic_slots if result.anchors[slot].status == STATUS_OK
|
||||
],
|
||||
}
|
||||
@@ -0,0 +1,357 @@
|
||||
"""第 4 步：逐帧分割掩膜 → 平滑多边形 → 位置约束 → 与绿框融合。
|
||||
|
||||
分割在灰度那一半（target_view）上做，交来的是这一半局部坐标的逐帧二值掩膜。这里把它们变成
|
||||
顶点数统一的平滑多边形轨迹，再用提示帧上的绿框（挪到 target_view 后）施加两类约束：
|
||||
|
||||
- `masks_to_smoothed_polygons`：逐帧取最大外轮廓 → 等弧长重采样到 `contour_points` 个顶点 →
|
||||
与上一帧对齐起点和朝向 → 空掩膜帧由相邻帧沿时间线性插值补齐（N 帧进 N 帧出，不剔除帧）→
|
||||
沿时间轴三角平滑；
|
||||
- `apply_prompt_frame_position_constraints`：提示帧绿框质心之间线性插值出逐帧中心轨迹，把每帧
|
||||
多边形按 `strength` 整体平移过去（不改形状和面积）；
|
||||
- `fuse_prompt_prior_with_medsam`：`roi_shape.source` 是 `prompt_prior` 时，由提示帧绿框推出逐帧
|
||||
先验多边形，逐顶点混合 `多边形 = 分割 × (1 − s) + 先验 × s`（提示帧用 `keyframe_shape_strength`，
|
||||
其余帧用 `propagation_shape_strength`），再把面积收紧到先验面积的 `[min, max]` 倍
|
||||
（绿框在最终轮廓里占多少，见 `待定问题/01`）。
|
||||
|
||||
照搬旧项目 `polygons.py`，去掉产品层；坐标换算改用两个视图矩形左上角之差（`_mapping.py`），
|
||||
不再需要旧的有效区矩形和映射锚点。
|
||||
"""
|
||||
|
||||
from typing import Any, Protocol
|
||||
|
||||
import cv2
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._errors import CeusCaseError
|
||||
from cadx.ceus._geometry import (
|
||||
AnyArray,
|
||||
Points,
|
||||
ViewRect,
|
||||
polygon_area,
|
||||
polygon_center,
|
||||
resample_closed_contour,
|
||||
triangular_smooth,
|
||||
)
|
||||
from cadx.ceus._mapping import map_source_mask_to_target
|
||||
|
||||
|
||||
class _SmoothingConfigLike(Protocol):
|
||||
"""平滑配置（`SmoothingConfig` 符合）。"""
|
||||
|
||||
@property
|
||||
def contour_points(self) -> int: ...
|
||||
@property
|
||||
def temporal_window(self) -> int: ...
|
||||
|
||||
|
||||
class _RoiShapeConfigLike(Protocol):
|
||||
"""绿框融合配置（`RoiShapeConfig` 符合）。"""
|
||||
|
||||
@property
|
||||
def source(self) -> str: ...
|
||||
@property
|
||||
def keyframe_shape_strength(self) -> float: ...
|
||||
@property
|
||||
def propagation_shape_strength(self) -> float: ...
|
||||
@property
|
||||
def min_area_ratio_to_prompt(self) -> float: ...
|
||||
@property
|
||||
def max_area_ratio_to_prompt(self) -> float: ...
|
||||
|
||||
|
||||
class _KeyframePositionConstraintLike(Protocol):
|
||||
"""提示帧位置约束配置（`KeyframePositionConstraint` 符合）。"""
|
||||
|
||||
@property
|
||||
def enabled(self) -> bool: ...
|
||||
@property
|
||||
def mode(self) -> str: ...
|
||||
@property
|
||||
def strength(self) -> float: ...
|
||||
|
||||
|
||||
_SUPPORTED_CONSTRAINT_MODE = "center_interpolation"
|
||||
|
||||
POLYGON_STEP = 4
|
||||
|
||||
|
||||
class NoValidContourError(CeusCaseError, RuntimeError):
|
||||
"""这一例的掩膜序列（分割结果或挪过来的绿框）一帧有效轮廓都没有，第 4 步做不下去。
|
||||
|
||||
同时是 `RuntimeError`：旧测试按 `RuntimeError` 断言，语义不变。
|
||||
"""
|
||||
|
||||
def __init__(self, reason: str) -> None:
|
||||
super().__init__(POLYGON_STEP, reason)
|
||||
|
||||
|
||||
def largest_contour(mask: AnyArray) -> Points | None:
|
||||
"""二值掩膜面积最大的外轮廓顶点 `(M, 2)` float32；没有轮廓返回 None。
|
||||
|
||||
取全部轮廓点（`CHAIN_APPROX_NONE`），方便后面等弧长重采样。
|
||||
"""
|
||||
contours, _ = cv2.findContours(mask.astype(np.uint8), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
|
||||
if not contours:
|
||||
return None
|
||||
best = max(contours, key=cv2.contourArea)
|
||||
return np.asarray(best).reshape(-1, 2).astype(np.float32)
|
||||
|
||||
|
||||
def mask_center(mask: AnyArray) -> Points:
|
||||
"""二值掩膜前景像素的质心 `[cx, cy]`；掩膜全空时返回 (0, 0)。"""
|
||||
ys, xs = np.where(mask)
|
||||
if len(xs) == 0:
|
||||
return np.array([0.0, 0.0], dtype=np.float32)
|
||||
return np.array([float(xs.mean()), float(ys.mean())], dtype=np.float32)
|
||||
|
||||
|
||||
def align_polygon_to_previous(poly: AnyArray, previous: AnyArray | None) -> AnyArray:
|
||||
"""把多边形的起点和朝向对齐到上一帧，让逐顶点的时间序列连续。
|
||||
|
||||
枚举正序、逆序两种朝向 × 全部循环移位，取与上一帧逐顶点平方差最小者；不对齐的话，时间
|
||||
平滑会把不同位置的顶点混在一起。`previous` 为 None（首帧）时原样返回。
|
||||
"""
|
||||
if previous is None:
|
||||
return poly
|
||||
candidates = [poly, poly[::-1].copy()]
|
||||
best: tuple[float, AnyArray] | None = None
|
||||
for candidate in candidates:
|
||||
for shift in range(len(candidate)):
|
||||
shifted = np.roll(candidate, shift, axis=0)
|
||||
score = float(((shifted - previous) ** 2).sum())
|
||||
if best is None or score < best[0]:
|
||||
best = (score, shifted)
|
||||
if best is None:
|
||||
raise ValueError("多边形没有顶点，无法对齐。")
|
||||
return best[1]
|
||||
|
||||
|
||||
def _interpolate_polygons(polygons: list[AnyArray | None], num_points: int) -> AnyArray:
|
||||
"""含 None 的多边形序列沿帧轴线性插值，补齐成 `(N, num_points, 2)`。
|
||||
|
||||
一帧有效轮廓都没有时无从插值，报 `RuntimeError`。
|
||||
"""
|
||||
valid_indices = np.asarray(
|
||||
[idx for idx, poly in enumerate(polygons) if poly is not None], dtype=np.int32
|
||||
)
|
||||
if len(valid_indices) == 0:
|
||||
raise NoValidContourError("逐帧 mask 未产出任何有效 ROI 轮廓，无从时间插值。")
|
||||
valid_values = np.asarray(
|
||||
[poly.reshape(-1) for poly in polygons if poly is not None], dtype=np.float32
|
||||
)
|
||||
frame_indices = np.arange(len(polygons), dtype=np.float32)
|
||||
interpolated = np.zeros((len(polygons), num_points * 2), dtype=np.float32)
|
||||
for col in range(num_points * 2):
|
||||
interpolated[:, col] = np.interp(
|
||||
frame_indices, valid_indices.astype(np.float32), valid_values[:, col]
|
||||
)
|
||||
return interpolated.reshape(len(polygons), num_points, 2)
|
||||
|
||||
|
||||
def _interpolate_centers(centers_by_frame: dict[int, AnyArray], frame_count: int) -> AnyArray:
|
||||
"""稀疏的提示帧中心沿帧轴线性插值成逐帧中心轨迹 `(N, 2)`，两端保持端值。"""
|
||||
frame_ids = np.asarray(sorted(centers_by_frame), dtype=np.float32)
|
||||
centers = np.asarray(
|
||||
[centers_by_frame[int(frame_id)] for frame_id in frame_ids.tolist()], dtype=np.float32
|
||||
)
|
||||
all_frames = np.arange(frame_count, dtype=np.float32)
|
||||
interpolated = np.zeros((frame_count, 2), dtype=np.float32)
|
||||
for col in range(2):
|
||||
interpolated[:, col] = np.interp(all_frames, frame_ids, centers[:, col])
|
||||
return interpolated
|
||||
|
||||
|
||||
def _clamp_polygon_to_view(polygon: AnyArray, target_rect: ViewRect) -> AnyArray:
|
||||
"""顶点裁到视图边界 `[0, w-1] × [0, h-1]`（返回副本）。"""
|
||||
out = polygon.copy()
|
||||
out[:, 0] = np.clip(out[:, 0], 0, target_rect.width - 1)
|
||||
out[:, 1] = np.clip(out[:, 1], 0, target_rect.height - 1)
|
||||
return out
|
||||
|
||||
|
||||
def _scale_polygon_about_center(polygon: AnyArray, center: AnyArray, scale: float) -> AnyArray:
|
||||
"""绕给定中心等比缩放多边形。"""
|
||||
return center[None, :] + (polygon - center[None, :]) * float(scale)
|
||||
|
||||
|
||||
def masks_to_smoothed_polygons(
|
||||
masks_by_frame: dict[int, AnyArray],
|
||||
frame_count: int,
|
||||
smoothing_cfg: _SmoothingConfigLike,
|
||||
) -> list[AnyArray]:
|
||||
"""逐帧掩膜 → 顶点统一、沿时间平滑的多边形序列，长度 `frame_count`。
|
||||
|
||||
缺帧或空掩膜视为退化帧：保留位置，由时间插值补齐。一帧有效轮廓都没有时报 `NoValidContourError`。
|
||||
"""
|
||||
num_points = int(smoothing_cfg.contour_points)
|
||||
temporal_window = int(smoothing_cfg.temporal_window)
|
||||
|
||||
aligned: list[AnyArray | None] = []
|
||||
previous: AnyArray | None = None
|
||||
for frame_idx in range(frame_count):
|
||||
mask = masks_by_frame.get(frame_idx)
|
||||
poly: AnyArray | None = None
|
||||
if mask is not None:
|
||||
poly = resample_closed_contour(largest_contour(mask), num_points)
|
||||
if poly is not None:
|
||||
poly = align_polygon_to_previous(poly, previous)
|
||||
previous = poly
|
||||
aligned.append(poly)
|
||||
|
||||
interpolated = _interpolate_polygons(aligned, num_points)
|
||||
smoothed = triangular_smooth(interpolated, temporal_window)
|
||||
return [smoothed[idx] for idx in range(frame_count)]
|
||||
|
||||
|
||||
def _prompt_frame_target_centers(
|
||||
prompt_frames: list[dict[str, Any]],
|
||||
source_view_rect: ViewRect,
|
||||
target_view_rect: ViewRect,
|
||||
) -> dict[int, AnyArray]:
|
||||
"""各提示帧的绿框掩膜挪到 target_view 后的质心，`帧号 → (cx, cy)`。"""
|
||||
centers: dict[int, AnyArray] = {}
|
||||
for prompt_frame in prompt_frames:
|
||||
target_mask = map_source_mask_to_target(
|
||||
prompt_frame["mask"], source_view_rect, target_view_rect
|
||||
)
|
||||
centers[int(prompt_frame["frame_index"])] = mask_center(target_mask)
|
||||
return centers
|
||||
|
||||
|
||||
def apply_prompt_frame_position_constraints(
|
||||
polygons: list[AnyArray],
|
||||
prompt_frames: list[dict[str, Any]],
|
||||
source_view_rect: ViewRect,
|
||||
target_view_rect: ViewRect,
|
||||
constraint_cfg: _KeyframePositionConstraintLike,
|
||||
) -> list[AnyArray]:
|
||||
"""按提示帧中心插值轨迹整体平移逐帧多边形（不改形状和面积），再裁到 target_view 内。
|
||||
|
||||
`prompt_frames` 每项含 `frame_index` 和造影那一半局部坐标的绿框 `mask`。每帧位移
|
||||
= (插值中心 − 当前质心) × strength，strength 截到 [0, 1]。`enabled=False` 时原对象返回；
|
||||
`mode` 不是 `center_interpolation` 报 `RuntimeError`。
|
||||
"""
|
||||
if not constraint_cfg.enabled:
|
||||
return polygons
|
||||
|
||||
if constraint_cfg.mode != _SUPPORTED_CONSTRAINT_MODE:
|
||||
raise RuntimeError(
|
||||
f"constraints.keyframe_position.mode 仅支持 '{_SUPPORTED_CONSTRAINT_MODE}'，"
|
||||
f"得到 '{constraint_cfg.mode}'。"
|
||||
)
|
||||
|
||||
strength = float(np.clip(constraint_cfg.strength, 0.0, 1.0))
|
||||
anchor_centers = _prompt_frame_target_centers(prompt_frames, source_view_rect, target_view_rect)
|
||||
center_track = _interpolate_centers(anchor_centers, len(polygons))
|
||||
|
||||
corrected: list[AnyArray] = []
|
||||
for frame_idx, polygon in enumerate(polygons):
|
||||
shift = (center_track[frame_idx] - polygon_center(polygon)) * strength
|
||||
corrected.append(_clamp_polygon_to_view(polygon + shift[None, :], target_view_rect))
|
||||
return corrected
|
||||
|
||||
|
||||
def _prompt_prior_polygons(
|
||||
prompt_frames: list[dict[str, Any]],
|
||||
frame_count: int,
|
||||
source_view_rect: ViewRect,
|
||||
target_view_rect: ViewRect,
|
||||
smoothing_cfg: _SmoothingConfigLike,
|
||||
) -> list[AnyArray]:
|
||||
"""由提示帧绿框推出逐帧先验多边形（target_view 局部坐标）。
|
||||
|
||||
每个提示帧的绿框挪到 target_view，取最大外轮廓重采样、对齐上一帧，放在它的帧号上；
|
||||
其余帧沿时间插值、平滑；最后提示帧那几帧用原始绿框多边形覆盖（不被平滑改动）。
|
||||
"""
|
||||
num_points = int(smoothing_cfg.contour_points)
|
||||
temporal_window = int(smoothing_cfg.temporal_window)
|
||||
|
||||
sparse: list[AnyArray | None] = [None] * frame_count
|
||||
raw_key_polygons: dict[int, AnyArray] = {}
|
||||
previous: AnyArray | None = None
|
||||
for prompt_frame in prompt_frames:
|
||||
target_mask = map_source_mask_to_target(
|
||||
prompt_frame["mask"], source_view_rect, target_view_rect
|
||||
)
|
||||
resampled = resample_closed_contour(largest_contour(target_mask), num_points)
|
||||
if resampled is None:
|
||||
continue
|
||||
poly = align_polygon_to_previous(resampled, previous)
|
||||
previous = poly
|
||||
frame_idx = int(prompt_frame["frame_index"])
|
||||
sparse[frame_idx] = poly
|
||||
raw_key_polygons[frame_idx] = poly
|
||||
|
||||
interpolated = _interpolate_polygons(sparse, num_points)
|
||||
if temporal_window > 1:
|
||||
interpolated = triangular_smooth(interpolated, temporal_window)
|
||||
for frame_idx, polygon in raw_key_polygons.items():
|
||||
interpolated[frame_idx] = polygon
|
||||
return [interpolated[idx] for idx in range(frame_count)]
|
||||
|
||||
|
||||
def _tighten_area_to_prior_band(
|
||||
polygon: AnyArray,
|
||||
prior_poly: AnyArray,
|
||||
min_area_ratio: float,
|
||||
max_area_ratio: float,
|
||||
) -> AnyArray:
|
||||
"""把多边形面积收紧到先验面积的 `[min, max]` 倍之内（绕先验质心按 sqrt(目标/当前) 缩放）。
|
||||
|
||||
先验面积或当前面积不超过 1 像素时不收紧（避免对退化面积做伪缩放）；已在带内原样返回。
|
||||
"""
|
||||
prior_area = polygon_area(prior_poly)
|
||||
current_area = polygon_area(polygon)
|
||||
if prior_area <= 1.0 or current_area <= 1.0:
|
||||
return polygon
|
||||
min_area = prior_area * min_area_ratio
|
||||
max_area = prior_area * max_area_ratio
|
||||
target_area = float(np.clip(current_area, min_area, max_area))
|
||||
if np.isclose(target_area, current_area):
|
||||
return polygon
|
||||
scale = float(np.sqrt(target_area / current_area))
|
||||
return _scale_polygon_about_center(polygon, polygon_center(prior_poly), scale)
|
||||
|
||||
|
||||
def fuse_prompt_prior_with_medsam(
|
||||
medsam_polygons: list[AnyArray],
|
||||
prompt_frames: list[dict[str, Any]],
|
||||
source_view_rect: ViewRect,
|
||||
target_view_rect: ViewRect,
|
||||
shape_cfg: _RoiShapeConfigLike,
|
||||
smoothing_cfg: _SmoothingConfigLike,
|
||||
) -> list[AnyArray]:
|
||||
"""分割多边形与绿框先验逐顶点融合并收紧面积，裁到 target_view 内。
|
||||
|
||||
`shape_cfg.source` 不是 `prompt_prior` 时原对象返回（不融合）。每帧先把分割多边形对齐到先验的
|
||||
起点和朝向，再按权重 s（截到 [0, 1]）混合。
|
||||
"""
|
||||
if shape_cfg.source != "prompt_prior":
|
||||
return medsam_polygons
|
||||
|
||||
prior_polygons = _prompt_prior_polygons(
|
||||
prompt_frames,
|
||||
len(medsam_polygons),
|
||||
source_view_rect,
|
||||
target_view_rect,
|
||||
smoothing_cfg,
|
||||
)
|
||||
prompt_frame_indices = {int(prompt_frame["frame_index"]) for prompt_frame in prompt_frames}
|
||||
prompt_frame_strength = float(shape_cfg.keyframe_shape_strength)
|
||||
propagation_strength = float(shape_cfg.propagation_shape_strength)
|
||||
min_area_ratio = float(shape_cfg.min_area_ratio_to_prompt)
|
||||
max_area_ratio = float(shape_cfg.max_area_ratio_to_prompt)
|
||||
|
||||
fused: list[AnyArray] = []
|
||||
for frame_idx, (medsam_poly, prior_poly) in enumerate(
|
||||
zip(medsam_polygons, prior_polygons, strict=True)
|
||||
):
|
||||
raw_strength = (
|
||||
prompt_frame_strength if frame_idx in prompt_frame_indices else propagation_strength
|
||||
)
|
||||
strength = float(np.clip(raw_strength, 0.0, 1.0))
|
||||
aligned_medsam = align_polygon_to_previous(medsam_poly, prior_poly)
|
||||
polygon = aligned_medsam * (1.0 - strength) + prior_poly * strength
|
||||
polygon = _tighten_area_to_prior_band(polygon, prior_poly, min_area_ratio, max_area_ratio)
|
||||
fused.append(_clamp_polygon_to_view(polygon, target_view_rect))
|
||||
return fused
|
||||
@@ -0,0 +1,408 @@
|
||||
"""第 6 步之前：在造影那一半的第一帧上选背景圆（`background_circle` 曲线的采样区）。
|
||||
|
||||
背景圆是 ROI 之外一小块干净背景，用来看整机增益漂移：病灶有了 ROI 曲线还不够，还要一块同帧
|
||||
的背景作对照，才能区分「造影剂真增强」和「全局亮度漂移」。它不是紫圈（`roi_pan`）：紫圈是医生
|
||||
画的正常胰腺对照，由 `_tic.py` 逐帧直接读颜色；背景圆是程序自动挑的，不喂任何口径字段。
|
||||
|
||||
`select_reference_region` 按模式分派：
|
||||
- `enabled=False`：空掩膜 + `{"enabled": False}`；
|
||||
- `manual`：以配置 `center_xy`（局部像素坐标）和 `radius_px` 直接画圆；
|
||||
- `auto_near_roi`：从 ROI 质心沿 8 个偏好方向（左下优先）× 若干搜索半径放候选圆心，避开 ROI
|
||||
膨胀 `min_distance_px` 的禁区，候选圆内「有效背景」（非黑、不是绿框色或紫圈色）占比须 ≥
|
||||
`min_valid_fraction`，按「有效占比 − 离 ROI 远近惩罚」择优；找不到时返回空掩膜和失败原因，
|
||||
背景曲线保持 NaN。
|
||||
|
||||
搜索半径倍率和偏好方向是本模块的常量，不在配置里。照搬旧项目 `reference_region.py`，去掉产品层。
|
||||
"""
|
||||
|
||||
from typing import Any, Protocol, cast
|
||||
|
||||
import cv2
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._config import TicAnalysisConfig
|
||||
from cadx.ceus._geometry import AnyArray, polygon_to_mask
|
||||
from cadx.ceus._green_prompt import prompt_color_mask
|
||||
|
||||
# auto 模式搜索半径（None 时）相对 ROI 等效尺度的倍率（参考 P:1058）。
|
||||
_SEARCH_RADIUS_FACTORS = (0.65, 0.85, 1.05, 1.25, 1.5)
|
||||
|
||||
# 8 个偏好方向（左下优先，参考 P:1062）：左下 / 左 / 下 / 右下 / 右 / 左上 / 上 / 右上。
|
||||
_PREFERRED_DIRECTIONS = (
|
||||
(-1.0, 1.0),
|
||||
(-1.0, 0.0),
|
||||
(0.0, 1.0),
|
||||
(1.0, 1.0),
|
||||
(1.0, 0.0),
|
||||
(-1.0, -1.0),
|
||||
(0.0, -1.0),
|
||||
(1.0, -1.0),
|
||||
)
|
||||
|
||||
# 评分中"离 ROI 距离"惩罚系数（参考 P:1091，硬编码不可改）。
|
||||
_DISTANCE_PENALTY_WEIGHT = 0.05
|
||||
|
||||
|
||||
class _ReferenceRegionConfigLike(Protocol):
|
||||
"""背景圆配置（`ReferenceRegionConfig` 符合）。"""
|
||||
|
||||
@property
|
||||
def enabled(self) -> bool: ...
|
||||
@property
|
||||
def mode(self) -> str: ...
|
||||
@property
|
||||
def radius_px(self) -> int: ...
|
||||
@property
|
||||
def min_distance_px(self) -> int: ...
|
||||
@property
|
||||
def nonblack_threshold(self) -> int: ...
|
||||
@property
|
||||
def min_valid_fraction(self) -> float: ...
|
||||
@property
|
||||
def center_xy(self) -> list[float] | None: ...
|
||||
|
||||
|
||||
class _PromptHsvConfigLike(Protocol):
|
||||
"""一组颜色的 HSV 阈值（`PromptSourceConfig`、`PanSourceConfig` 都符合）。"""
|
||||
|
||||
@property
|
||||
def hsv_lower(self) -> list[int]: ...
|
||||
@property
|
||||
def hsv_upper(self) -> list[int]: ...
|
||||
|
||||
|
||||
def _circle_mask(shape: tuple[int, int], center_xy: AnyArray, radius: float) -> AnyArray:
|
||||
"""以 ``center_xy`` 为圆心、``radius`` 为半径栅格化实心圆 mask（参考 P:975-979）。
|
||||
|
||||
参数：
|
||||
shape: 目标 mask 形状 ``(height, width)``。
|
||||
center_xy: 圆心 ``[cx, cy]`` 像素坐标（可为浮点）。
|
||||
radius: 圆半径（像素）。
|
||||
|
||||
返回：
|
||||
``dtype=bool`` 的 ``(height, width)`` mask，圆内（含边界，``<= r²``）为 ``True``。
|
||||
|
||||
关键实现：
|
||||
用 ``np.ogrid`` 广播算每像素到圆心的平方距离，与 ``radius²`` 比较——闭区间
|
||||
（``<=``）与参考一致，避免边界像素被漏。
|
||||
"""
|
||||
height, width = shape
|
||||
yy, xx = np.ogrid[:height, :width]
|
||||
cx, cy = float(center_xy[0]), float(center_xy[1])
|
||||
return (xx - cx) ** 2 + (yy - cy) ** 2 <= float(radius) ** 2
|
||||
|
||||
|
||||
def _polygon_union_mask(polygons: AnyArray, height: int, width: int) -> AnyArray:
|
||||
"""逐帧 ROI 多边形并集 mask（参考 P:982-986）。
|
||||
|
||||
参数：
|
||||
polygons: 逐帧 ROI 多边形序列（每帧 ``(N, 2)``）。
|
||||
height: mask 高度。
|
||||
width: mask 宽度。
|
||||
|
||||
返回：
|
||||
所有帧 ROI 多边形按位或的 ``dtype=bool`` 并集 mask——禁区取并集而非单帧，
|
||||
保证参考圆在整段视频内都不与任何帧的 ROI 重叠。
|
||||
"""
|
||||
union = np.zeros((height, width), dtype=bool)
|
||||
for polygon in polygons:
|
||||
union |= polygon_to_mask(polygon, height, width)
|
||||
return union
|
||||
|
||||
|
||||
def _forbidden_roi_mask(roi_union: AnyArray, min_distance_px: int) -> AnyArray:
|
||||
"""把 ROI 并集按 ``min_distance_px`` 膨胀为参考圆禁区（参考 P:1040-1045）。
|
||||
|
||||
参数：
|
||||
roi_union: ROI 并集 mask（bool）。
|
||||
min_distance_px: 参考圆与 ROI 的最小间距（像素）；``<= 0`` 时不膨胀。
|
||||
|
||||
返回：
|
||||
膨胀后的 ``dtype=bool`` 禁区 mask——参考圆不得触及其内任何像素。
|
||||
|
||||
关键实现：
|
||||
膨胀核边长 ``2·min_distance_px + 1``（方核），与参考一致；``<= 0`` 直接返回原并集
|
||||
（禁区即 ROI 本身）。
|
||||
"""
|
||||
if min_distance_px <= 0:
|
||||
return roi_union
|
||||
kernel = np.ones((min_distance_px * 2 + 1, min_distance_px * 2 + 1), dtype=np.uint8)
|
||||
return cv2.dilate(roi_union.astype(np.uint8) * 255, kernel, iterations=1) > 0
|
||||
|
||||
|
||||
def _search_radii(roi_scale: float) -> list[float]:
|
||||
"""求 auto 模式的候选搜索半径列表（参考 P:1056-1058）。
|
||||
|
||||
参数：
|
||||
roi_scale: ROI 等效尺度（bbox 长 / 宽较大者，下限 1.0）。
|
||||
|
||||
返回：
|
||||
候选圆心到 ROI 质心的距离列表：``roi_scale × _SEARCH_RADIUS_FACTORS``。
|
||||
|
||||
关键实现：
|
||||
倍率不从配置读——``ReferenceRegionConfig`` 里没有这一项，而它是 ``extra="forbid"``，
|
||||
往 yaml 里写 ``search_radii_px`` 会在加载时被拒。要做成可配就得先进 schema。
|
||||
"""
|
||||
return [roi_scale * factor for factor in _SEARCH_RADIUS_FACTORS]
|
||||
|
||||
|
||||
def _preferred_directions() -> AnyArray:
|
||||
"""auto 模式的偏好方向数组（参考 P:1059-1065）。
|
||||
|
||||
返回：
|
||||
``(8, 2)`` float32 方向数组（左下优先排序）——评分相同则靠前方向先入选。
|
||||
|
||||
关键实现：
|
||||
方向表不从配置读，同 ``_search_radii``：``preferred_directions`` 不在 schema 里，
|
||||
写进 yaml 会被 ``extra="forbid"`` 拒掉。
|
||||
"""
|
||||
return np.asarray(_PREFERRED_DIRECTIONS, dtype=np.float32)
|
||||
|
||||
|
||||
def _candidate_centers(
|
||||
roi_center: AnyArray,
|
||||
search_radii: list[float],
|
||||
directions: AnyArray,
|
||||
) -> list[AnyArray]:
|
||||
"""生成 auto 模式的候选圆心列表（参考 P:1068-1073 的外两层循环）。
|
||||
|
||||
参数：
|
||||
roi_center: ROI 质心 ``[cx, cy]``（float32）。
|
||||
search_radii: 候选距离列表。
|
||||
directions: ``(K, 2)`` 偏好方向数组。
|
||||
|
||||
返回：
|
||||
候选圆心 ``[cx, cy]`` 列表——半径 × 方向（单位化）双重笛卡尔积；零向量跳过。
|
||||
外层遍历距离、内层遍历方向，保持与参考一致的先小半径、左下优先的枚举顺序。
|
||||
"""
|
||||
centers: list[AnyArray] = []
|
||||
for distance in search_radii:
|
||||
for row in range(directions.shape[0]):
|
||||
direction = directions[row]
|
||||
norm = float(cast(float, np.linalg.norm(direction)))
|
||||
if norm <= 1e-6:
|
||||
continue
|
||||
centers.append(roi_center + direction / norm * float(distance))
|
||||
return centers
|
||||
|
||||
|
||||
def _circle_within_bounds(center: AnyArray, radius: float, height: int, width: int) -> bool:
|
||||
"""判断以 ``center`` 为心、``radius`` 为半径的圆是否完整落在画面内（参考 P:1074-1079）。"""
|
||||
return not (
|
||||
center[0] - radius < 0
|
||||
or center[1] - radius < 0
|
||||
or center[0] + radius >= width
|
||||
or center[1] + radius >= height
|
||||
)
|
||||
|
||||
|
||||
def _valid_fraction(mask: AnyArray, valid_background: AnyArray) -> float:
|
||||
"""候选圆内"有效背景"像素占比（参考 P:1085）。
|
||||
|
||||
参数：
|
||||
mask: 候选圆 mask（bool）。
|
||||
valid_background: 有效背景 mask（非黑、非医生画的线、非禁区，bool）。
|
||||
|
||||
返回：
|
||||
``(圆 ∩ 有效背景) / 圆面积``；圆面积为 0 时返回 0.0（不会被选中）。
|
||||
"""
|
||||
area = int(mask.sum())
|
||||
if area == 0:
|
||||
return 0.0
|
||||
return float((mask & valid_background).sum()) / float(area)
|
||||
|
||||
|
||||
def _score_candidate(
|
||||
center: AnyArray,
|
||||
valid_fraction: float,
|
||||
first_polygon: AnyArray,
|
||||
roi_scale: float,
|
||||
) -> float:
|
||||
"""候选圆评分：有效占比越高越好，离 ROI 越远略减分（参考 P:1088-1091）。
|
||||
|
||||
参数：
|
||||
center: 候选圆心 ``[cx, cy]``。
|
||||
valid_fraction: 该候选圆内有效背景占比。
|
||||
first_polygon: 第一帧 ROI 多边形（``(N, 2)``），用于算到 ROI 的最近点距离。
|
||||
roi_scale: ROI 等效尺度（归一化距离用，下限 1.0）。
|
||||
|
||||
返回：
|
||||
``valid_fraction − 到ROI最近距离 / roi_scale × 0.05``——倾向"贴近 ROI 的干净背景"
|
||||
（同帧、就近，光照 / 探头条件最一致），但有效占比仍是主导项。
|
||||
"""
|
||||
distance_to_roi = float(np.sqrt(((center[None, :] - first_polygon) ** 2).sum(axis=1)).min())
|
||||
return valid_fraction - distance_to_roi / max(roi_scale, 1.0) * _DISTANCE_PENALTY_WEIGHT
|
||||
|
||||
|
||||
def _roi_scale(first_polygon: AnyArray) -> float:
|
||||
"""ROI 等效尺度：bbox 长 / 宽中较大者，下限 1.0（参考 P:1053-1055）。"""
|
||||
bbox_min = first_polygon.min(axis=0)
|
||||
bbox_max = first_polygon.max(axis=0)
|
||||
return max(float(bbox_max[0] - bbox_min[0]), float(bbox_max[1] - bbox_min[1]), 1.0)
|
||||
|
||||
|
||||
def _manual_reference(
|
||||
shape: tuple[int, int], ref_cfg: _ReferenceRegionConfigLike, radius: float
|
||||
) -> tuple[AnyArray, dict[str, Any]]:
|
||||
"""manual 模式：以配置 ``center_xy`` + ``radius`` 直接出圆 + metadata（参考 P:1021-1033）。
|
||||
|
||||
参数：
|
||||
shape: 第一帧 ``(height, width)``。
|
||||
ref_cfg: 参考区配置（取 ``center_xy``）。
|
||||
radius: 圆半径（像素）。
|
||||
|
||||
返回：
|
||||
``(圆 mask, metadata dict)``；metadata 含 enabled/mode/shape/center_xy/radius_px/area_px。
|
||||
|
||||
关键实现：
|
||||
``center_xy`` 形状必须为 ``(2,)``，否则视为配置缺陷直接报错（不静默兜底）。schema
|
||||
的 ``model_validator`` 已强制 manual 模式 center_xy 非空，这里再校验维度做纵深防御。
|
||||
"""
|
||||
center = np.asarray(ref_cfg.center_xy, dtype=np.float32)
|
||||
if center.shape != (2,):
|
||||
raise RuntimeError("tic_analysis.reference_region.center_xy 必须是 [x, y]。")
|
||||
mask = _circle_mask(shape, center, radius)
|
||||
return mask, {
|
||||
"enabled": True,
|
||||
"mode": "manual",
|
||||
"shape": "circle",
|
||||
"center_xy": [float(center[0]), float(center[1])],
|
||||
"radius_px": radius,
|
||||
"area_px": int(mask.sum()),
|
||||
}
|
||||
|
||||
|
||||
def _auto_reference(
|
||||
first_source_bgr: AnyArray,
|
||||
source_polygons: AnyArray,
|
||||
ref_cfg: _ReferenceRegionConfigLike,
|
||||
prompt_cfg: _PromptHsvConfigLike,
|
||||
pan_cfg: _PromptHsvConfigLike,
|
||||
tic_cfg: TicAnalysisConfig,
|
||||
radius: float,
|
||||
) -> tuple[AnyArray, dict[str, Any]]:
|
||||
"""auto_near_roi 模式：ROI 外搜索一块有效背景圆择优（参考 P:1038-1109）。
|
||||
|
||||
参数：
|
||||
first_source_bgr: 第一帧 source 视图 BGR uint8 ``(h, w, 3)``（已是 source 局部坐标）。
|
||||
source_polygons: 逐帧 ROI 多边形序列（禁区取并集，避开整段 ROI）。
|
||||
ref_cfg: 参考区配置。
|
||||
prompt_cfg: 剔绿 HSV 阈值（透传给 ``prompt_color_mask``）。
|
||||
pan_cfg: 剔紫 HSV 阈值（同一路传法，见「关键实现」为什么紫色也要排）。
|
||||
tic_cfg: TIC 配置（剔绿开关 / 膨胀核，透传给 ``prompt_color_mask``）。
|
||||
radius: 圆半径（像素）。
|
||||
|
||||
返回：
|
||||
``(选中圆 mask, metadata dict)``；metadata 额外含 ``valid_fraction_first_frame``。
|
||||
|
||||
关键实现：
|
||||
有效背景 = 非黑（灰度 > ``nonblack_threshold``）且不是医生画上去的线（绿框色与
|
||||
紫色轮廓色都排掉）且非禁区。逐候选圆需满足：完整在画面内、不触禁区、有效占比
|
||||
≥ ``min_valid_fraction``，再按 ``_score_candidate`` 取最高分。无合规候选返回空掩膜和
|
||||
``failure_reason=no_valid_candidate``；背景曲线据此保持 NaN，不阻断病灶拟合。
|
||||
|
||||
**紫色和绿色一起排。** 紫色轮廓是医生画的邻近正常胰腺实质参照区，TIC 采样时把
|
||||
这些像素从掩膜里扣掉（``tic._sample_frame_row`` 的排除色 = 绿 ∪ 紫）。这里若只
|
||||
排绿，挑中的背景圆可能压在那圈紫色上，``valid_background_circle_px`` 被悄悄削掉
|
||||
一块，而没有任何一处会报出来——紫色的灰度不低，``nonblack_threshold`` 也拦不住它。
|
||||
"""
|
||||
height, width = first_source_bgr.shape[:2]
|
||||
roi_union = _polygon_union_mask(source_polygons, height, width)
|
||||
forbidden_roi = _forbidden_roi_mask(roi_union, int(ref_cfg.min_distance_px))
|
||||
|
||||
ink_mask = prompt_color_mask(first_source_bgr, prompt_cfg, tic_cfg) | prompt_color_mask(
|
||||
first_source_bgr, pan_cfg, tic_cfg
|
||||
)
|
||||
gray = cv2.cvtColor(first_source_bgr, cv2.COLOR_BGR2GRAY)
|
||||
nonblack = gray > int(ref_cfg.nonblack_threshold)
|
||||
valid_background = nonblack & ~ink_mask & ~forbidden_roi
|
||||
|
||||
first_polygon = source_polygons[0]
|
||||
roi_center = first_polygon.mean(axis=0).astype(np.float32)
|
||||
roi_scale = _roi_scale(first_polygon)
|
||||
search_radii = _search_radii(roi_scale)
|
||||
directions = _preferred_directions()
|
||||
min_valid_fraction = float(ref_cfg.min_valid_fraction)
|
||||
|
||||
best: tuple[float, AnyArray, AnyArray, float] | None = None
|
||||
for candidate_center in _candidate_centers(roi_center, search_radii, directions):
|
||||
if not _circle_within_bounds(candidate_center, radius, height, width):
|
||||
continue
|
||||
mask = _circle_mask((height, width), candidate_center, radius)
|
||||
if int(mask.sum()) == 0 or np.any(mask & forbidden_roi):
|
||||
continue
|
||||
valid_fraction = _valid_fraction(mask, valid_background)
|
||||
if valid_fraction < min_valid_fraction:
|
||||
continue
|
||||
score = _score_candidate(candidate_center, valid_fraction, first_polygon, roi_scale)
|
||||
if best is None or score > best[0]:
|
||||
best = (score, candidate_center, mask, valid_fraction)
|
||||
|
||||
if best is None:
|
||||
return np.zeros((height, width), dtype=bool), {
|
||||
"enabled": True,
|
||||
"mode": "auto_near_roi",
|
||||
"shape": "circle",
|
||||
"available": False,
|
||||
"failure_reason": "no_valid_candidate",
|
||||
}
|
||||
|
||||
_, selected_center, selected_mask, selected_fraction = best
|
||||
return selected_mask, {
|
||||
"enabled": True,
|
||||
"mode": "auto_near_roi",
|
||||
"shape": "circle",
|
||||
"center_xy": [float(selected_center[0]), float(selected_center[1])],
|
||||
"radius_px": radius,
|
||||
"area_px": int(selected_mask.sum()),
|
||||
"valid_fraction_first_frame": float(selected_fraction),
|
||||
}
|
||||
|
||||
|
||||
def select_reference_region(
|
||||
first_source_bgr: AnyArray,
|
||||
source_polygons: AnyArray,
|
||||
prompt_cfg: _PromptHsvConfigLike,
|
||||
pan_cfg: _PromptHsvConfigLike,
|
||||
tic_cfg: TicAnalysisConfig,
|
||||
) -> tuple[AnyArray, dict[str, Any]]:
|
||||
"""在第一帧 source 视图上选定 TIC 参考区圆，出 mask + 元数据（参考 P:1007-1109）。
|
||||
|
||||
参数：
|
||||
first_source_bgr: 第一帧 source 视图 BGR uint8 ``(h, w, 3)``（已是 source 局部坐标）。
|
||||
source_polygons: 逐帧 ROI 多边形序列 ``(T, N, 2)``（参考区只看第一帧定位，禁区取整段并集）。
|
||||
prompt_cfg: 绿框 HSV 阈值配置（剔绿用，透传 ``prompt_color_mask``）。
|
||||
pan_cfg: 紫色轮廓 HSV 阈值配置（剔紫用，同一路传法）。
|
||||
tic_cfg: TIC 分析配置；从中取 ``reference_region`` 子配置与剔绿开关 / 膨胀核。
|
||||
|
||||
返回：
|
||||
``(reference mask, metadata dict)``——
|
||||
mask：``dtype=bool`` 的 ``(h, w)`` 圆 mask（TIC 据此采样 ``background_circle`` 曲线）；
|
||||
metadata：写进 trajectory.npz 的 ``reference_region_json``。各模式的键集：
|
||||
- enabled=False → ``{"enabled": False}``（空 mask）；
|
||||
- manual → ``{enabled, mode, shape, center_xy, radius_px, area_px}``；
|
||||
- auto_near_roi → 额外含 ``valid_fraction_first_frame``。
|
||||
``shape`` 恒为 ``"circle"``（出自旧项目 findings §2B：输出死键，不从配置读）。
|
||||
|
||||
异常：
|
||||
RuntimeError: ``mode`` 非法，或 manual 的 center_xy 维度非法。
|
||||
|
||||
关键实现：
|
||||
三模式分派——关闭直接空返；manual 直接画圆；auto 走 ``_auto_reference`` 搜索择优。
|
||||
参考区配置经 ``tic_cfg.reference_region`` 取出，不在本模块写默认值。
|
||||
"""
|
||||
ref_cfg = tic_cfg.reference_region
|
||||
|
||||
if not bool(ref_cfg.enabled):
|
||||
empty = np.zeros(first_source_bgr.shape[:2], dtype=bool)
|
||||
return empty, {"enabled": False}
|
||||
|
||||
radius = float(ref_cfg.radius_px)
|
||||
mode = ref_cfg.mode
|
||||
if mode == "manual":
|
||||
return _manual_reference(first_source_bgr.shape[:2], ref_cfg, radius)
|
||||
if mode != "auto_near_roi":
|
||||
raise RuntimeError(f"不支持的 tic_analysis.reference_region.mode: {mode}")
|
||||
return _auto_reference(
|
||||
first_source_bgr, source_polygons, ref_cfg, prompt_cfg, pan_cfg, tic_cfg, radius
|
||||
)
|
||||
@@ -0,0 +1,438 @@
|
||||
"""第 6 步：TIC 四条曲线的逐帧采样、增强前基线，以及（调第 7 步）Lognormal 拟合。
|
||||
|
||||
TIC（时间-强度曲线）是病灶亮度随时间的变化。在造影那一半画面（局部坐标）上逐帧采四条
|
||||
灰度均值曲线：
|
||||
|
||||
- `roi_all`：病灶 ROI 内（剔除绿框色与紫圈色后）的灰度均值；
|
||||
- `roi_perf`：在 `roi_all` 基础上再剔除暗区后的灰度均值（perf 指灌注，见 `_tic_curves.py`）；
|
||||
- `roi_pan`：医生画的紫圈内部（一块正常胰腺，作对照）；
|
||||
- `background_circle`：程序自动挑的一小块背景圆，只用来看整机增益漂移，不是 `roi_pan`。
|
||||
|
||||
掩膜构成：
|
||||
|
||||
exclude = 绿框色 ∪ 紫圈色
|
||||
valid_roi_all = roi_mask & ~exclude
|
||||
valid_roi_perf = valid_roi_all & ~dark_union
|
||||
valid_pan = 紫圈填充区 & ~exclude
|
||||
valid_background = background_circle_mask & ~exclude
|
||||
|
||||
四条原始曲线采好后，每条各自识别增强前基线（`detect_baseline`），再交给 `_tic_fitting` 拟合。
|
||||
曲线最前面不插「时间 0、值 0」的合成点：那个点不是量出来的，会在它与首个真实帧之间造出
|
||||
一段最陡的假上升，把「变亮最快」的时刻钉在视频开头。
|
||||
|
||||
照搬旧项目 `tic.py`，去掉产品层；另加 `pan_detected_fraction`（紫圈检出率，任务书 02 要报）。
|
||||
"""
|
||||
|
||||
from collections.abc import Iterable, Sequence
|
||||
from dataclasses import dataclass
|
||||
from typing import Any
|
||||
|
||||
import cv2
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._config import (
|
||||
PanSourceConfig,
|
||||
PromptSourceConfig,
|
||||
TicAnalysisConfig,
|
||||
TicBaselineConfig,
|
||||
)
|
||||
from cadx.ceus._geometry import AnyArray, Frame, Mask, polygon_to_mask, polygons_to_mask
|
||||
from cadx.ceus._green_prompt import extract_color_prompt, prompt_color_mask
|
||||
from cadx.ceus._tic_fitting import LognormalFitResult, fit_lognormal_curve
|
||||
|
||||
__all__ = [
|
||||
"CURVE_NAMES",
|
||||
"BaselineResult",
|
||||
"TicResult",
|
||||
"compute_tic_for_video",
|
||||
"detect_baseline",
|
||||
"mean_or_nan",
|
||||
]
|
||||
|
||||
# 四条曲线的名字；逐帧行键按 `{name}_raw` / `{name}_processed` 等派生。
|
||||
CURVE_NAMES = ("roi_all", "roi_perf", "roi_pan", "background_circle")
|
||||
|
||||
# 临床口径认得的三条。背景圆不在里面：它不喂任何口径字段，基线被截短的质控位也不由它置位。
|
||||
_SPEC_CURVE_NAMES = ("roi_all", "roi_perf", "roi_pan")
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class BaselineResult:
|
||||
"""一条曲线的增强前基线区间（还没被造影剂点亮的那一段）。
|
||||
|
||||
`mean`、`sd` 是区间内的稳健中心和噪声（中位数、1.4826×MAD），`valid=False` 时为 None；
|
||||
`frame_start` 含、`frame_end` 不含，无效时都是 -1。无效时用 None 而不用 0.0：0.0 能一路
|
||||
算下去把下游悄悄算错，None 会在第一个用它的地方报错。
|
||||
"""
|
||||
|
||||
mean: float | None
|
||||
sd: float | None
|
||||
frame_start: int
|
||||
frame_end: int
|
||||
valid: bool
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class TicResult:
|
||||
"""TIC 采样与拟合的结果。
|
||||
|
||||
- `rows`：逐帧一行（只有真实帧），键见 `compute_tic_for_video`；
|
||||
- `parameters`：定量参数，这里恒为 None（灌注指标由 `_perfusion_metrics` 从拟合算）；
|
||||
- `baselines`、`fits`：曲线名 → 基线 / 拟合，键集恒为 `CURVE_NAMES`；
|
||||
- `early_truncation_flag`：口径三条曲线里有一条采到过值、却凑不出足够长的基线；
|
||||
- `reference_valid`：紫圈检出率达标，`roi_pan` 可用；为 False 时 `roi_pan` 整条 NaN；
|
||||
- `pan_detected_fraction`：非重复帧里真正采到紫圈像素的帧占比（紫圈检出率）。
|
||||
"""
|
||||
|
||||
rows: list[dict[str, Any]]
|
||||
parameters: dict[str, Any] | None
|
||||
baselines: dict[str, BaselineResult]
|
||||
fits: dict[str, LognormalFitResult]
|
||||
early_truncation_flag: bool
|
||||
reference_valid: bool
|
||||
pan_detected_fraction: float
|
||||
|
||||
|
||||
def mean_or_nan(gray: AnyArray, mask: AnyArray) -> float:
|
||||
"""掩膜内灰度均值；掩膜全空时返回 NaN（这一帧这条曲线没采到，不当作亮度 0）。"""
|
||||
values = gray[mask]
|
||||
if len(values) == 0:
|
||||
return float("nan")
|
||||
return float(values.mean())
|
||||
|
||||
|
||||
_INVALID_BASELINE = BaselineResult(mean=None, sd=None, frame_start=-1, frame_end=-1, valid=False)
|
||||
|
||||
|
||||
def _robust_center_noise(values: AnyArray) -> tuple[float, float]:
|
||||
"""中位数与 1.4826×MAD（正态下等于标准差的稳健估计）。"""
|
||||
center = float(np.median(values))
|
||||
noise = 1.4826 * float(np.median(np.abs(values - center)))
|
||||
return center, noise
|
||||
|
||||
|
||||
def _window_is_stable(window: AnyArray, sd_multiplier: float) -> bool | None:
|
||||
"""候选窗对半分，后半中位数有没有明显高过前半。
|
||||
|
||||
稳定返回 True、已在上涨返回 False；任一半没有有效采样、判不出来时返回 None。
|
||||
"""
|
||||
half = len(window) // 2
|
||||
first = window[:half]
|
||||
second = window[half:]
|
||||
first = first[np.isfinite(first)]
|
||||
second = second[np.isfinite(second)]
|
||||
if first.size == 0 or second.size == 0:
|
||||
return None
|
||||
first_center, first_noise = _robust_center_noise(first)
|
||||
second_center = float(np.median(second))
|
||||
first_differences = np.diff(first)
|
||||
if first_differences.size:
|
||||
difference_center = float(np.median(first_differences))
|
||||
difference_noise = 1.4826 * float(np.median(np.abs(first_differences - difference_center)))
|
||||
else:
|
||||
difference_noise = first_noise
|
||||
return bool(
|
||||
second_center - first_center
|
||||
<= sd_multiplier * max(difference_noise, float(np.finfo(np.float64).eps))
|
||||
)
|
||||
|
||||
|
||||
def detect_baseline(
|
||||
values: AnyArray,
|
||||
time_s: AnyArray,
|
||||
cfg: TicBaselineConfig,
|
||||
sample_mask: AnyArray | None,
|
||||
) -> BaselineResult:
|
||||
"""在曲线开头按真实时间找增强前基线。
|
||||
|
||||
候选窗从视频开头覆盖 `preferred_duration_s` 秒；对半分后后半中位数比前半高出超过
|
||||
`stability_sd_multiplier × 前半差分噪声`，就从右侧逐个样本缩窗；有效样本覆盖的时长不足
|
||||
`minimum_duration_s` 时基线无效。NaN 帧和 `sample_mask` 为 False 的帧不算基线样本
|
||||
（`sample_mask` 为 None 表示都算）。
|
||||
|
||||
这一步只出基线均值和噪声，不出到达时刻和峰值：到达门槛含峰高，要先找峰，所以到达时刻
|
||||
在拟合那一步回头定。
|
||||
"""
|
||||
series = np.asarray(values, dtype=np.float64)
|
||||
times = np.asarray(time_s, dtype=np.float64)
|
||||
mask = (
|
||||
np.ones(series.shape, dtype=bool)
|
||||
if sample_mask is None
|
||||
else np.asarray(sample_mask, dtype=bool)
|
||||
)
|
||||
if series.ndim != 1 or times.shape != series.shape or mask.shape != series.shape:
|
||||
raise ValueError("values、time_s 与 sample_mask 必须是等长一维数组。")
|
||||
if series.size == 0 or not np.isfinite(times).all() or np.any(np.diff(times) <= 0):
|
||||
raise ValueError("time_s 必须非空、有限且严格递增。")
|
||||
|
||||
window_end_s = float(times[0]) + cfg.preferred_duration_s
|
||||
end = int(np.searchsorted(times, window_end_s, side="left"))
|
||||
end = min(max(end, 1), int(series.size))
|
||||
cadence = float(np.median(np.diff(times))) if series.size > 1 else 0.0
|
||||
|
||||
while end > 0:
|
||||
window_mask = np.isfinite(series[:end]) & mask[:end]
|
||||
finite_indices = np.flatnonzero(window_mask)
|
||||
if finite_indices.size < 2:
|
||||
return _INVALID_BASELINE
|
||||
observed_duration = float(times[finite_indices[-1]] - times[finite_indices[0]]) + cadence
|
||||
if observed_duration + 1e-9 < cfg.minimum_duration_s:
|
||||
return _INVALID_BASELINE
|
||||
window = np.where(window_mask, series[:end], np.nan)
|
||||
stable = _window_is_stable(window, cfg.stability_sd_multiplier)
|
||||
if stable is None:
|
||||
return _INVALID_BASELINE
|
||||
if stable:
|
||||
finite = series[:end][window_mask]
|
||||
center, noise = _robust_center_noise(finite)
|
||||
return BaselineResult(
|
||||
mean=center,
|
||||
sd=noise,
|
||||
frame_start=0,
|
||||
frame_end=end,
|
||||
valid=True,
|
||||
)
|
||||
end -= 1
|
||||
return _INVALID_BASELINE
|
||||
|
||||
|
||||
def _sample_frame_row(
|
||||
frame_index: int,
|
||||
time_s: float,
|
||||
source_bgr: Frame,
|
||||
roi_polygon: AnyArray,
|
||||
frame_dark_polygons: Sequence[AnyArray],
|
||||
background_circle_mask: Mask,
|
||||
prompt_cfg: PromptSourceConfig,
|
||||
pan_cfg: PanSourceConfig,
|
||||
tic_cfg: TicAnalysisConfig,
|
||||
) -> tuple[dict[str, Any], AnyArray]:
|
||||
"""一帧的四掩膜采样，返回这一帧的原始行和判断相邻重复帧用的 ROI 灰度签名。
|
||||
|
||||
紫圈的线本身是医生画的，不是组织，所以排除色是「绿 ∪ 紫」。`roi_pan` 取紫圈实心填充
|
||||
再减排除色（轮廓线本身不算参照区内部），这一帧检不到紫圈时 `pan_area_px=0`、
|
||||
`roi_pan_raw` 为 NaN；紫圈不走分割、不做传播，画面上每帧都有，直接读。
|
||||
|
||||
`pan_area_px` 是紫圈填充面积（未扣排除色），`valid_pan_px` 是扣完后真正参与统计的像素数；
|
||||
检出率数后者。暗区多边形一次奇偶填充成掩膜（洞不计回暗区），再收口到 `valid_roi_all`；
|
||||
收口后暗区面积为 0 时暗块计数也记 0。
|
||||
"""
|
||||
height, width = source_bgr.shape[:2]
|
||||
gray = np.asarray(cv2.cvtColor(source_bgr, cv2.COLOR_BGR2GRAY))
|
||||
exclude_mask = prompt_color_mask(source_bgr, prompt_cfg, tic_cfg) | prompt_color_mask(
|
||||
source_bgr, pan_cfg, tic_cfg
|
||||
)
|
||||
|
||||
roi_mask = polygon_to_mask(roi_polygon, height, width)
|
||||
valid_roi_all_mask = roi_mask & ~exclude_mask
|
||||
valid_background_mask = background_circle_mask & ~exclude_mask
|
||||
|
||||
dark_union = polygons_to_mask(frame_dark_polygons, height, width)
|
||||
dark_union &= valid_roi_all_mask
|
||||
valid_roi_perf_mask = valid_roi_all_mask & ~dark_union
|
||||
|
||||
pan_prompt = extract_color_prompt(source_bgr, pan_cfg)
|
||||
if pan_prompt is None:
|
||||
pan_mask = np.zeros((height, width), dtype=bool)
|
||||
pan_area_px = 0
|
||||
else:
|
||||
pan_mask = pan_prompt.mask & ~exclude_mask
|
||||
pan_area_px = int(pan_prompt.mask.sum())
|
||||
|
||||
dark_area_px = int(dark_union.sum())
|
||||
dark_component_count = len(frame_dark_polygons) if dark_area_px > 0 else 0
|
||||
|
||||
x_coord, y_coord, box_width, box_height = cv2.boundingRect(
|
||||
np.asarray(roi_polygon, dtype=np.float32)
|
||||
)
|
||||
roi_gray = gray[y_coord : y_coord + box_height, x_coord : x_coord + box_width]
|
||||
roi_crop_mask = roi_mask[y_coord : y_coord + box_height, x_coord : x_coord + box_width]
|
||||
signature_gray = np.asarray(cv2.resize(roi_gray, (32, 32), interpolation=cv2.INTER_AREA))
|
||||
signature_mask = np.asarray(
|
||||
cv2.resize(roi_crop_mask.astype(np.uint8), (32, 32), interpolation=cv2.INTER_NEAREST)
|
||||
).astype(bool)
|
||||
signature = np.where(signature_mask, signature_gray, 0).astype(np.float32)
|
||||
|
||||
row: dict[str, Any] = {
|
||||
"frame_index": frame_index,
|
||||
"time_sec": float(time_s),
|
||||
"roi_all_raw": mean_or_nan(gray, valid_roi_all_mask),
|
||||
"roi_perf_raw": mean_or_nan(gray, valid_roi_perf_mask),
|
||||
"roi_pan_raw": mean_or_nan(gray, pan_mask),
|
||||
"background_circle_raw": mean_or_nan(gray, valid_background_mask),
|
||||
"roi_area_px": int(roi_mask.sum()),
|
||||
"dark_area_px": dark_area_px,
|
||||
"pan_area_px": pan_area_px,
|
||||
"background_circle_area_px": int(background_circle_mask.sum()),
|
||||
"valid_roi_all_px": int(valid_roi_all_mask.sum()),
|
||||
"valid_roi_perf_px": int(valid_roi_perf_mask.sum()),
|
||||
"valid_pan_px": int(pan_mask.sum()),
|
||||
"valid_background_circle_px": int(valid_background_mask.sum()),
|
||||
"dark_component_count": dark_component_count,
|
||||
}
|
||||
return row, signature
|
||||
|
||||
|
||||
def _attach_baselines_and_fits(
|
||||
rows: list[dict[str, Any]], tic_cfg: TicAnalysisConfig
|
||||
) -> tuple[dict[str, BaselineResult], dict[str, LognormalFitResult], bool]:
|
||||
"""逐条曲线识别基线并拟合，把候选、拟合、正式曲线原地写回逐帧行。
|
||||
|
||||
每条曲线各算各的基线窗：胰腺比病灶早增强，共用一个窗会把 `roi_pan` 的窗拖进已在上涨的区间。
|
||||
`roi_perf` 只用「去暗区面积占比与上一帧相比变化不超过 `maximum_perf_area_step_fraction`」
|
||||
的帧（`perf_mask_stable`），暗区掩膜跳变的帧不进基线和拟合。
|
||||
|
||||
`early_truncation_flag` 只看口径三条：某条采到过实测值、却凑不出足够长的基线就置位；
|
||||
整条 NaN 的曲线不算（那是没检出，由 `reference_valid` 管）；背景圆不参与。
|
||||
"""
|
||||
baselines: dict[str, BaselineResult] = {}
|
||||
fits: dict[str, LognormalFitResult] = {}
|
||||
early_truncation_flag = False
|
||||
times = np.asarray([row["time_sec"] for row in rows], dtype=np.float64)
|
||||
duplicates = np.asarray([row["duplicate_frame"] for row in rows], dtype=bool)
|
||||
roi_areas = np.asarray([row["roi_area_px"] for row in rows], dtype=np.float64)
|
||||
perf_areas = np.asarray([row["valid_roi_perf_px"] for row in rows], dtype=np.float64)
|
||||
perf_fractions = np.divide(
|
||||
perf_areas,
|
||||
roi_areas,
|
||||
out=np.zeros_like(perf_areas),
|
||||
where=roi_areas > 0.0,
|
||||
)
|
||||
perf_mask_stable = roi_areas > 0.0
|
||||
if len(rows) > 1:
|
||||
stable_steps = (
|
||||
np.abs(np.diff(perf_fractions)) <= tic_cfg.fitting.maximum_perf_area_step_fraction
|
||||
)
|
||||
perf_mask_stable[1:] &= stable_steps
|
||||
for index, row in enumerate(rows):
|
||||
row["perf_mask_stable"] = bool(perf_mask_stable[index])
|
||||
for name in CURVE_NAMES:
|
||||
raw = np.asarray([row[f"{name}_raw"] for row in rows], dtype=np.float64)
|
||||
curve_mask = perf_mask_stable if name == "roi_perf" else np.ones(len(rows), dtype=bool)
|
||||
baseline = detect_baseline(raw, times, tic_cfg.baseline, ~duplicates & curve_mask)
|
||||
baselines[name] = baseline
|
||||
if name in _SPEC_CURVE_NAMES and not baseline.valid and bool(np.isfinite(raw).any()):
|
||||
early_truncation_flag = True
|
||||
fit = fit_lognormal_curve(
|
||||
times,
|
||||
raw,
|
||||
baseline=baseline.mean,
|
||||
baseline_noise=baseline.sd,
|
||||
baseline_end_index=baseline.frame_end,
|
||||
duplicate_mask=duplicates,
|
||||
sample_mask=curve_mask,
|
||||
cfg=tic_cfg.fitting,
|
||||
anchor_cfg=tic_cfg.anchor_frames,
|
||||
)
|
||||
fits[name] = fit
|
||||
for idx, row in enumerate(rows):
|
||||
row[f"{name}_candidate"] = float(fit.candidate[idx])
|
||||
row[f"{name}_fitted"] = float(fit.fitted[idx])
|
||||
row[f"{name}_processed"] = float(fit.processed[idx])
|
||||
return baselines, fits, early_truncation_flag
|
||||
|
||||
|
||||
def compute_tic_for_video(
|
||||
source_frames: Iterable[Frame],
|
||||
roi_polygons: Sequence[AnyArray],
|
||||
dark_polygons_per_frame: Sequence[Sequence[AnyArray]],
|
||||
reference_mask: AnyArray,
|
||||
fps: float,
|
||||
prompt_cfg: PromptSourceConfig,
|
||||
pan_cfg: PanSourceConfig,
|
||||
tic_cfg: TicAnalysisConfig,
|
||||
frame_times_s: Sequence[float] | AnyArray | None,
|
||||
) -> TicResult:
|
||||
"""逐帧采四掩膜灰度 → 标记重复帧 → 基线与 Lognormal 拟合，交出 `TicResult`。
|
||||
|
||||
- `source_frames`：造影那一半的逐帧 BGR（局部坐标），逐帧消费、不留，可以直接喂解码迭代器；
|
||||
- `roi_polygons`、`dark_polygons_per_frame`：与帧等长的逐帧 ROI 多边形和暗区多边形列表；
|
||||
- `reference_mask`：背景圆掩膜（第一帧定位、整段沿用）；
|
||||
- `frame_times_s`：逐帧显示时间戳（秒）；为 None 时按 `fps` 构造等间隔时间轴（只给合成测试用）。
|
||||
|
||||
`rows` 每行键：`frame_index`、`time_sec`、`duplicate_frame`、`perf_mask_stable`，四条曲线各自的
|
||||
`_raw`、`_candidate`、`_fitted`、`_processed`，以及各掩膜面积和有效像素数、暗块计数。
|
||||
|
||||
紫圈检出率 = 非重复帧里 `valid_pan_px > 0` 的帧占比（数的是真正参与统计的像素，不是填充面积：
|
||||
一整块实心紫色会让填充面积过门槛而扣完排除色一个像素不剩）。达到
|
||||
`pan_cfg.min_detected_frame_fraction` 才算紫圈可用；不达标时 `roi_pan` 整条记 NaN、基线无效，
|
||||
不拿背景圆去填（两者不是同一种东西）。
|
||||
"""
|
||||
if len(roi_polygons) != len(dark_polygons_per_frame):
|
||||
raise ValueError(
|
||||
"roi_polygons 与 dark_polygons_per_frame 长度必须一致，得到 "
|
||||
f"{len(roi_polygons)} / {len(dark_polygons_per_frame)}。"
|
||||
)
|
||||
if float(fps) <= 0.0:
|
||||
raise ValueError(f"fps 必须为正，得到 {fps}。")
|
||||
if tic_cfg.brightness_metric != "grayscale_mean":
|
||||
raise ValueError(
|
||||
f"仅支持 brightness_metric=grayscale_mean，得到 {tic_cfg.brightness_metric}。"
|
||||
)
|
||||
background_circle_mask = reference_mask.astype(bool)
|
||||
if frame_times_s is None:
|
||||
times = np.arange(len(roi_polygons), dtype=np.float64) / float(fps)
|
||||
else:
|
||||
times = np.asarray(frame_times_s, dtype=np.float64)
|
||||
if times.shape != (len(roi_polygons),):
|
||||
raise ValueError(
|
||||
"frame_times_s 必须与 roi_polygons 等长，得到 "
|
||||
f"{times.shape} / {len(roi_polygons)}。"
|
||||
)
|
||||
if not np.isfinite(times).all() or np.any(np.diff(times) <= 0.0):
|
||||
raise ValueError("frame_times_s 必须有限且严格递增。")
|
||||
|
||||
# 帧是流式的，长度要遍历完才知道：zip(strict=True) 保证三个序列等长，多了少了当场报错。
|
||||
rows: list[dict[str, Any]] = []
|
||||
previous_unique_signature: AnyArray | None = None
|
||||
for frame_index, (frame, roi_polygon, dark_polygons) in enumerate(
|
||||
zip(source_frames, roi_polygons, dark_polygons_per_frame, strict=True)
|
||||
):
|
||||
if frame_index == 0 and background_circle_mask.shape != frame.shape[:2]:
|
||||
raise ValueError(
|
||||
f"reference_mask 形状 {background_circle_mask.shape} 与 source 帧 "
|
||||
f"{frame.shape[:2]} 不符。"
|
||||
)
|
||||
row, signature = _sample_frame_row(
|
||||
frame_index,
|
||||
float(times[frame_index]),
|
||||
frame,
|
||||
roi_polygon,
|
||||
dark_polygons,
|
||||
background_circle_mask,
|
||||
prompt_cfg,
|
||||
pan_cfg,
|
||||
tic_cfg,
|
||||
)
|
||||
duplicate = False
|
||||
if tic_cfg.fitting.detect_duplicate_frames and previous_unique_signature is not None:
|
||||
difference = float(np.mean(np.abs(signature - previous_unique_signature)))
|
||||
duplicate = difference <= tic_cfg.fitting.duplicate_mean_absdiff_threshold
|
||||
row["duplicate_frame"] = duplicate
|
||||
rows.append(row)
|
||||
if not duplicate:
|
||||
previous_unique_signature = signature
|
||||
if not rows:
|
||||
raise ValueError("source_frames 不能为空。")
|
||||
|
||||
unique_rows = [row for row in rows if not row["duplicate_frame"]]
|
||||
detected_frames = sum(1 for row in unique_rows if row["valid_pan_px"] > 0)
|
||||
pan_detected_fraction = detected_frames / len(unique_rows) if unique_rows else 0.0
|
||||
reference_valid = bool(unique_rows) and (
|
||||
pan_detected_fraction >= pan_cfg.min_detected_frame_fraction
|
||||
)
|
||||
if not reference_valid:
|
||||
for row in rows:
|
||||
row["roi_pan_raw"] = float("nan")
|
||||
|
||||
baselines, fits, early_truncation_flag = _attach_baselines_and_fits(rows, tic_cfg)
|
||||
|
||||
return TicResult(
|
||||
rows=rows,
|
||||
parameters=None,
|
||||
baselines=baselines,
|
||||
fits=fits,
|
||||
early_truncation_flag=early_truncation_flag,
|
||||
reference_valid=reference_valid,
|
||||
pan_detected_fraction=pan_detected_fraction,
|
||||
)
|
||||
@@ -0,0 +1,436 @@
|
||||
"""tic_curves.npz 的打包和解包：四条曲线连同基线、拟合、面积、有效像素数（格式同旧项目）。
|
||||
|
||||
曲线数值单独存一份，问「第 37 秒那条曲线是多少」不必重跑整段视频。
|
||||
|
||||
- 解包用 `np.load(..., allow_pickle=False)`，坏字节直接报错；
|
||||
- 打包前校验各逐帧数组长度一致；
|
||||
- 字符串存 0 维 str 数组，不写 object 数组；
|
||||
- NaN 原样写入：它是「这一帧没采到」的编码，本模块不用 `np.nan_to_num`，填 0 会把
|
||||
「没采到」变成「亮度是 0」而不报错。
|
||||
|
||||
照搬旧项目 `tic_curves.py`，去掉产品层。
|
||||
"""
|
||||
|
||||
import io
|
||||
from dataclasses import dataclass
|
||||
from typing import Any, cast
|
||||
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._geometry import AnyArray
|
||||
from cadx.ceus._tic import CURVE_NAMES, TicResult
|
||||
|
||||
# 这个文件自己的布局版本。键集、dtype、形状任一处变了就要换它。
|
||||
CURVES_SCHEMA_VERSION = "3"
|
||||
|
||||
# 逐帧整数列的键名（面积 / 有效像素数 / 暗块计数），行键与工件键同名。
|
||||
_COUNT_KEYS = (
|
||||
"roi_area_px",
|
||||
"dark_area_px",
|
||||
"pan_area_px",
|
||||
"background_circle_area_px",
|
||||
"valid_roi_all_px",
|
||||
"valid_roi_perf_px",
|
||||
"valid_pan_px",
|
||||
"valid_background_circle_px",
|
||||
"dark_component_count",
|
||||
)
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class TicCurveSeries:
|
||||
"""一条曲线的原始样本、候选曲线、Lognormal 结果与拟合质量。
|
||||
|
||||
``fitted`` 保留优化器得到但可能未过质量门的曲线，便于复核失败；``processed`` 只在
|
||||
质量门通过时等于 ``fitted``，失败时整条为 NaN。这样候选或失败拟合不会冒充正式曲线。
|
||||
"""
|
||||
|
||||
raw: AnyArray
|
||||
candidate: AnyArray
|
||||
fitted: AnyArray
|
||||
processed: AnyArray
|
||||
valid_sample_mask: AnyArray
|
||||
duplicate_sample_mask: AnyArray
|
||||
baseline_mean: float
|
||||
baseline_sd: float
|
||||
baseline_frames: tuple[int, int]
|
||||
baseline_valid: bool
|
||||
fit_success: bool
|
||||
fit_failure_reason: str
|
||||
curve_model: str
|
||||
smoothing_method: str
|
||||
baseline_source: str
|
||||
fit_area: float
|
||||
fit_mu: float
|
||||
fit_sigma: float
|
||||
fit_arrival_offset_s: float
|
||||
fit_baseline: float
|
||||
fit_r2: float
|
||||
fit_nrmse: float
|
||||
residual_noise_ratio: float
|
||||
fit_interval_s: tuple[float, float]
|
||||
multi_peak_flag: bool
|
||||
anchor_times_s: tuple[float, float, float, float]
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class TicCurves:
|
||||
"""TIC 曲线工件的内存表示（npz pack/unpack 的双向载体）。
|
||||
|
||||
属性：
|
||||
frame_indices: ``(T,)`` int32，逐帧下标。``T`` 是真实帧数，**不含任何合成点**。
|
||||
time_video_s: ``(T,)`` float32，来自 MP4 presentation timestamp。
|
||||
injection_offset_s: 推注偏移（秒）；取不到时为 ``float("nan")``，**不是 0**
|
||||
——填 0 等于断言「视频和推注同时开始」，而那是我们没有的信息。
|
||||
time_origin: ``"injection"`` 或 ``"video_start"``，说明绝对时间轴以谁为零点。
|
||||
fps: 视频帧率。
|
||||
curves_schema_version: 这个文件自己的布局版本。
|
||||
curves: 曲线名 → ``TicCurveSeries``，键集恒为 ``tic.CURVE_NAMES``。
|
||||
roi_area_px / dark_area_px / pan_area_px / background_circle_area_px:
|
||||
``(T,)`` int32，逐帧各掩膜栅格面积。
|
||||
valid_roi_all_px / valid_roi_perf_px / valid_pan_px /
|
||||
valid_background_circle_px: ``(T,)`` int32，逐帧真正参与统计的像素数。
|
||||
dark_component_count: ``(T,)`` int32，逐帧暗块个数。
|
||||
reference_valid: 紫色轮廓检出率达标，``roi_pan`` 那条可用。
|
||||
early_truncation_flag: 有曲线没能凑出足够长的基线区间。
|
||||
|
||||
**``roi_perf`` 这个名字断言了「暗区 = 无灌注区」这个等号。** 「灌注」指血液连同
|
||||
造影剂真正流进一块组织，没有灌注的地方（囊液、坏死、出血）不会被点亮；``_perf``
|
||||
后缀在规范里指「去掉无灌注区之后剩下的可灌注组织」。而我们量到的是**暗区**——
|
||||
病变轮廓内部灰度明显偏低的那一块，它不必然等于无灌注区，**伪影**（超声成像自身
|
||||
产生的假象，比如前方有气体或钙化挡住声波、后方整片变暗的声影）也会让有灌注的组织
|
||||
显得暗。**这个等号是工科方在 2026-08-12 会上提出的，临床方带着「肯定会有被伪影干扰」
|
||||
这个保留意见暂时认下**（旧项目会议转录 475 与 477 行；旧 design/0030 决策七据此写成）。
|
||||
2026-08-16 项目负责人拍板先按这个等号做、暂不为伪影做任何区分，将来它被推翻时，要回头
|
||||
重算的正是所有带 ``_perf`` 后缀的量。
|
||||
|
||||
**三样东西这里不存**，各有各的理由：
|
||||
|
||||
``S(t) = raw − B`` 不存。``B`` 已经在文件里了，再存一遍就是同一个数放两份，
|
||||
迟早有一份被改、另一份没改。
|
||||
|
||||
``t_injection = time_video_s + injection_offset_s`` 不存。同上，两个加数都在。
|
||||
|
||||
无灌注比例不存。临床方原规范 §6 给了公式 ``(A_all − A_perf) / A_all``，但没给它的时间下标
|
||||
——取哪个时点还挂在临床方名下。逐帧的分子分母（``roi_area_px`` 与 ``dark_area_px``）
|
||||
都在文件里，时点定下来之后是一次派生，不必重跑视频。
|
||||
"""
|
||||
|
||||
frame_indices: AnyArray
|
||||
time_video_s: AnyArray
|
||||
injection_offset_s: float
|
||||
time_origin: str
|
||||
fps: float
|
||||
config_version: str
|
||||
algorithm_version: str
|
||||
curves_schema_version: str
|
||||
curves: dict[str, TicCurveSeries]
|
||||
roi_area_px: AnyArray
|
||||
dark_area_px: AnyArray
|
||||
pan_area_px: AnyArray
|
||||
background_circle_area_px: AnyArray
|
||||
valid_roi_all_px: AnyArray
|
||||
valid_roi_perf_px: AnyArray
|
||||
valid_pan_px: AnyArray
|
||||
valid_background_circle_px: AnyArray
|
||||
dark_component_count: AnyArray
|
||||
perf_mask_stable: AnyArray
|
||||
reference_valid: bool
|
||||
early_truncation_flag: bool
|
||||
|
||||
|
||||
def _nan_if_none(value: float | None) -> float:
|
||||
"""缺失的时刻在 npz 里记 NaN（float 数组没有 null）。"""
|
||||
return float("nan") if value is None else float(value)
|
||||
|
||||
|
||||
def from_tic_result(
|
||||
result: TicResult,
|
||||
fps: float,
|
||||
injection_offset_s: float | None,
|
||||
config_version: str,
|
||||
) -> TicCurves:
|
||||
"""把 ``TicResult`` 转成 ``TicCurves``，调用方不必知道行键名。
|
||||
|
||||
参数：
|
||||
result: ``tic.compute_tic_for_video`` 的产物。
|
||||
fps: 视频帧率（>0）。
|
||||
injection_offset_s: 推注偏移（秒）；源资料里没有就传 ``None``。
|
||||
config_version: 生成本工件的科研配置版本。
|
||||
|
||||
返回：
|
||||
``TicCurves``。
|
||||
|
||||
异常：
|
||||
ValueError: ``result.rows`` 为空，或 ``fps`` 非正。
|
||||
|
||||
``time_origin`` 由 ``injection_offset_s`` 派生，不单独作为入参：两者只有一种自洽的
|
||||
搭配（有偏移才谈得上以推注为零点），拆成两个参数只是多一个能填错的地方。
|
||||
"""
|
||||
if float(fps) <= 0.0:
|
||||
raise ValueError(f"fps 必须为正，得到 {fps}。")
|
||||
rows = result.rows
|
||||
if not rows:
|
||||
raise ValueError("TicResult.rows 不能为空。")
|
||||
if not config_version:
|
||||
raise ValueError("config_version 不能为空。")
|
||||
algorithm_versions = {result.fits[name].algorithm_version for name in CURVE_NAMES}
|
||||
if len(algorithm_versions) != 1:
|
||||
raise ValueError(f"四条曲线的 algorithm_version 不一致：{algorithm_versions}")
|
||||
|
||||
curves: dict[str, TicCurveSeries] = {}
|
||||
for name in CURVE_NAMES:
|
||||
baseline = result.baselines[name]
|
||||
fit = result.fits[name]
|
||||
parameters = fit.parameters
|
||||
curves[name] = TicCurveSeries(
|
||||
raw=np.asarray([row[f"{name}_raw"] for row in rows], dtype=np.float32),
|
||||
candidate=np.asarray([row[f"{name}_candidate"] for row in rows], dtype=np.float32),
|
||||
fitted=np.asarray([row[f"{name}_fitted"] for row in rows], dtype=np.float32),
|
||||
processed=np.asarray([row[f"{name}_processed"] for row in rows], dtype=np.float32),
|
||||
valid_sample_mask=fit.valid_sample_mask.astype(bool),
|
||||
duplicate_sample_mask=fit.duplicate_sample_mask.astype(bool),
|
||||
baseline_mean=float("nan") if baseline.mean is None else baseline.mean,
|
||||
baseline_sd=float("nan") if baseline.sd is None else baseline.sd,
|
||||
baseline_frames=(baseline.frame_start, baseline.frame_end),
|
||||
baseline_valid=baseline.valid,
|
||||
fit_success=fit.success,
|
||||
fit_failure_reason=fit.failure_reason,
|
||||
curve_model=fit.curve_model,
|
||||
smoothing_method=fit.smoothing_method,
|
||||
baseline_source=fit.baseline_source,
|
||||
fit_area=float("nan") if parameters is None else parameters.area,
|
||||
fit_mu=float("nan") if parameters is None else parameters.mu,
|
||||
fit_sigma=float("nan") if parameters is None else parameters.sigma,
|
||||
fit_arrival_offset_s=(
|
||||
float("nan") if parameters is None else parameters.arrival_offset_s
|
||||
),
|
||||
fit_baseline=float("nan") if parameters is None else parameters.baseline,
|
||||
fit_r2=fit.r2,
|
||||
fit_nrmse=fit.nrmse,
|
||||
residual_noise_ratio=fit.residual_noise_ratio,
|
||||
fit_interval_s=(fit.fit_start_s, fit.fit_end_s),
|
||||
multi_peak_flag=fit.multi_peak_flag,
|
||||
anchor_times_s=(
|
||||
_nan_if_none(fit.d1_time_s),
|
||||
_nan_if_none(fit.d2_time_s),
|
||||
_nan_if_none(fit.d3_time_s),
|
||||
_nan_if_none(fit.d4_time_s),
|
||||
),
|
||||
)
|
||||
|
||||
counts = {key: np.asarray([row[key] for row in rows], dtype=np.int32) for key in _COUNT_KEYS}
|
||||
frame_indices = np.asarray([row["frame_index"] for row in rows], dtype=np.int32)
|
||||
return TicCurves(
|
||||
frame_indices=frame_indices,
|
||||
time_video_s=np.asarray([row["time_sec"] for row in rows], dtype=np.float32),
|
||||
injection_offset_s=(
|
||||
float("nan") if injection_offset_s is None else float(injection_offset_s)
|
||||
),
|
||||
time_origin="video_start" if injection_offset_s is None else "injection",
|
||||
fps=float(fps),
|
||||
config_version=config_version,
|
||||
algorithm_version=algorithm_versions.pop(),
|
||||
curves_schema_version=CURVES_SCHEMA_VERSION,
|
||||
curves=curves,
|
||||
perf_mask_stable=np.asarray([row["perf_mask_stable"] for row in rows], dtype=bool),
|
||||
reference_valid=result.reference_valid,
|
||||
early_truncation_flag=result.early_truncation_flag,
|
||||
**counts,
|
||||
)
|
||||
|
||||
|
||||
def _validate_lengths(curves: TicCurves) -> int:
|
||||
"""校验逐帧列长度一致，返回帧数 T；任何不一致直接 raise。"""
|
||||
if curves.frame_indices.ndim != 1:
|
||||
raise ValueError(f"frame_indices 必须一维，得到 {curves.frame_indices.shape}。")
|
||||
num_frames = int(curves.frame_indices.shape[0])
|
||||
curve_names = set(curves.curves)
|
||||
expected_names = set(CURVE_NAMES)
|
||||
if curve_names != expected_names:
|
||||
raise ValueError(f"曲线键集必须为 {sorted(expected_names)}，得到 {sorted(curve_names)}。")
|
||||
if curves.curves_schema_version != CURVES_SCHEMA_VERSION:
|
||||
raise ValueError(
|
||||
f"只支持 curves schema {CURVES_SCHEMA_VERSION}，得到 {curves.curves_schema_version}。"
|
||||
)
|
||||
if not curves.config_version or not curves.algorithm_version:
|
||||
raise ValueError("config_version 与 algorithm_version 不能为空。")
|
||||
if num_frames == 0 or np.any(np.diff(curves.frame_indices) <= 0):
|
||||
raise ValueError("frame_indices 必须非空且严格递增。")
|
||||
if not np.isfinite(curves.time_video_s).all() or np.any(np.diff(curves.time_video_s) <= 0.0):
|
||||
raise ValueError("time_video_s 必须有限且严格递增。")
|
||||
|
||||
named: list[tuple[str, AnyArray]] = [("time_video_s", curves.time_video_s)]
|
||||
named += [(key, getattr(curves, key)) for key in _COUNT_KEYS]
|
||||
named.append(("perf_mask_stable", curves.perf_mask_stable))
|
||||
for name in CURVE_NAMES:
|
||||
series = curves.curves[name]
|
||||
named.append((f"{name}_raw", series.raw))
|
||||
named.append((f"{name}_candidate", series.candidate))
|
||||
named.append((f"{name}_fitted", series.fitted))
|
||||
named.append((f"{name}_processed", series.processed))
|
||||
named.append((f"{name}_valid_sample_mask", series.valid_sample_mask))
|
||||
named.append((f"{name}_duplicate_sample_mask", series.duplicate_sample_mask))
|
||||
for key, array in named:
|
||||
if array.shape != (num_frames,):
|
||||
raise ValueError(f"{key} 必须形如 ({num_frames},)，得到 {array.shape}。")
|
||||
|
||||
return num_frames
|
||||
|
||||
|
||||
def pack_tic_curves(curves: TicCurves) -> bytes:
|
||||
"""把 ``TicCurves`` 校验并打包为压缩 npz 字节（产物工件）。
|
||||
|
||||
参数：
|
||||
curves: 待打包的曲线载体。
|
||||
|
||||
返回：
|
||||
压缩 npz 的原始字节。
|
||||
|
||||
异常：
|
||||
ValueError: 逐帧列长度不一致，或曲线名缺失。
|
||||
|
||||
关键实现：
|
||||
先 ``_validate_lengths`` 钉死长度（不一致即 raise，不静默对齐）；字符串键存为
|
||||
0 维 str 数组，标量存 ``(1,)`` 数组，布尔存 ``int8``，全程不写入任何 object
|
||||
数组（禁 pickle）。NaN 原样写入 float32——它是「这一帧没采到」的编码，
|
||||
任何形式的填零都会把这个事实抹掉。
|
||||
"""
|
||||
_validate_lengths(curves)
|
||||
|
||||
payload: dict[str, AnyArray] = {
|
||||
"frame_indices": curves.frame_indices.astype(np.int32),
|
||||
"time_video_s": curves.time_video_s.astype(np.float32),
|
||||
"injection_offset_s": np.asarray([curves.injection_offset_s], dtype=np.float32),
|
||||
"time_origin": np.array(curves.time_origin),
|
||||
"fps": np.asarray([curves.fps], dtype=np.float32),
|
||||
"config_version": np.array(curves.config_version),
|
||||
"algorithm_version": np.array(curves.algorithm_version),
|
||||
"curves_schema_version": np.array(curves.curves_schema_version),
|
||||
"reference_valid": np.asarray([int(curves.reference_valid)], dtype=np.int8),
|
||||
"early_truncation_flag": np.asarray([int(curves.early_truncation_flag)], dtype=np.int8),
|
||||
"perf_mask_stable": curves.perf_mask_stable.astype(np.int8),
|
||||
}
|
||||
for key in _COUNT_KEYS:
|
||||
payload[key] = np.asarray(getattr(curves, key), dtype=np.int32)
|
||||
for name in CURVE_NAMES:
|
||||
series = curves.curves[name]
|
||||
payload[f"{name}_raw"] = series.raw.astype(np.float32)
|
||||
payload[f"{name}_candidate"] = series.candidate.astype(np.float32)
|
||||
payload[f"{name}_fitted"] = series.fitted.astype(np.float32)
|
||||
payload[f"{name}_processed"] = series.processed.astype(np.float32)
|
||||
payload[f"{name}_valid_sample_mask"] = series.valid_sample_mask.astype(np.int8)
|
||||
payload[f"{name}_duplicate_sample_mask"] = series.duplicate_sample_mask.astype(np.int8)
|
||||
payload[f"{name}_baseline_mean"] = np.asarray([series.baseline_mean], dtype=np.float32)
|
||||
payload[f"{name}_baseline_sd"] = np.asarray([series.baseline_sd], dtype=np.float32)
|
||||
payload[f"{name}_baseline_frames"] = np.asarray(series.baseline_frames, dtype=np.int32)
|
||||
payload[f"{name}_baseline_valid"] = np.asarray([int(series.baseline_valid)], dtype=np.int8)
|
||||
payload[f"{name}_fit_success"] = np.asarray([int(series.fit_success)], dtype=np.int8)
|
||||
payload[f"{name}_fit_failure_reason"] = np.array(series.fit_failure_reason)
|
||||
payload[f"{name}_curve_model"] = np.array(series.curve_model)
|
||||
payload[f"{name}_smoothing_method"] = np.array(series.smoothing_method)
|
||||
payload[f"{name}_baseline_source"] = np.array(series.baseline_source)
|
||||
payload[f"{name}_fit_parameters"] = np.asarray(
|
||||
[
|
||||
series.fit_area,
|
||||
series.fit_mu,
|
||||
series.fit_sigma,
|
||||
series.fit_arrival_offset_s,
|
||||
series.fit_baseline,
|
||||
],
|
||||
dtype=np.float64,
|
||||
)
|
||||
payload[f"{name}_fit_quality"] = np.asarray(
|
||||
[series.fit_r2, series.fit_nrmse, series.residual_noise_ratio],
|
||||
dtype=np.float64,
|
||||
)
|
||||
payload[f"{name}_fit_interval_s"] = np.asarray(series.fit_interval_s, dtype=np.float64)
|
||||
payload[f"{name}_multi_peak_flag"] = np.asarray(
|
||||
[int(series.multi_peak_flag)], dtype=np.int8
|
||||
)
|
||||
payload[f"{name}_anchor_times_s"] = np.asarray(series.anchor_times_s, dtype=np.float64)
|
||||
|
||||
buffer = io.BytesIO()
|
||||
# 存根把 **payload 也当成可能填进 allow_pickle 等关键字参数，这里只有数组键。
|
||||
cast(Any, np.savez_compressed)(buffer, **payload)
|
||||
return buffer.getvalue()
|
||||
|
||||
|
||||
def unpack_tic_curves(raw: bytes) -> TicCurves:
|
||||
"""把 npz 字节解包还原为 ``TicCurves``（禁 pickle，坏字节 raise）。
|
||||
|
||||
参数：
|
||||
raw: ``pack_tic_curves`` 产出的 npz 原始字节。
|
||||
|
||||
返回：
|
||||
还原后的 ``TicCurves``，dtype 与 pack 时一致。
|
||||
|
||||
异常：
|
||||
ValueError / OSError / KeyError: 非法 npz、含 object pickle 成员、或缺键。
|
||||
|
||||
关键实现：
|
||||
``np.load(..., allow_pickle=False)`` 硬禁 pickle——任何含 object pickle 的恶意
|
||||
字节都会在访问该成员时抛错而非被反序列化执行。NaN 原样取回，不做任何替换。
|
||||
"""
|
||||
with np.load(io.BytesIO(raw), allow_pickle=False) as archive:
|
||||
schema_version = str(archive["curves_schema_version"].item())
|
||||
if schema_version != CURVES_SCHEMA_VERSION:
|
||||
raise ValueError(
|
||||
f"只支持 curves schema {CURVES_SCHEMA_VERSION}，得到 {schema_version}。"
|
||||
)
|
||||
curves: dict[str, TicCurveSeries] = {}
|
||||
for name in CURVE_NAMES:
|
||||
frames = archive[f"{name}_baseline_frames"].astype(np.int32)
|
||||
fit_parameters = archive[f"{name}_fit_parameters"].astype(np.float64)
|
||||
fit_quality = archive[f"{name}_fit_quality"].astype(np.float64)
|
||||
fit_interval = archive[f"{name}_fit_interval_s"].astype(np.float64)
|
||||
anchor_times = archive[f"{name}_anchor_times_s"].astype(np.float64)
|
||||
curves[name] = TicCurveSeries(
|
||||
raw=archive[f"{name}_raw"].astype(np.float32),
|
||||
candidate=archive[f"{name}_candidate"].astype(np.float32),
|
||||
fitted=archive[f"{name}_fitted"].astype(np.float32),
|
||||
processed=archive[f"{name}_processed"].astype(np.float32),
|
||||
valid_sample_mask=archive[f"{name}_valid_sample_mask"].astype(bool),
|
||||
duplicate_sample_mask=archive[f"{name}_duplicate_sample_mask"].astype(bool),
|
||||
baseline_mean=float(archive[f"{name}_baseline_mean"][0]),
|
||||
baseline_sd=float(archive[f"{name}_baseline_sd"][0]),
|
||||
baseline_frames=(int(frames[0]), int(frames[1])),
|
||||
baseline_valid=bool(archive[f"{name}_baseline_valid"][0]),
|
||||
fit_success=bool(archive[f"{name}_fit_success"][0]),
|
||||
fit_failure_reason=str(archive[f"{name}_fit_failure_reason"].item()),
|
||||
curve_model=str(archive[f"{name}_curve_model"].item()),
|
||||
smoothing_method=str(archive[f"{name}_smoothing_method"].item()),
|
||||
baseline_source=str(archive[f"{name}_baseline_source"].item()),
|
||||
fit_area=float(fit_parameters[0]),
|
||||
fit_mu=float(fit_parameters[1]),
|
||||
fit_sigma=float(fit_parameters[2]),
|
||||
fit_arrival_offset_s=float(fit_parameters[3]),
|
||||
fit_baseline=float(fit_parameters[4]),
|
||||
fit_r2=float(fit_quality[0]),
|
||||
fit_nrmse=float(fit_quality[1]),
|
||||
residual_noise_ratio=float(fit_quality[2]),
|
||||
fit_interval_s=(float(fit_interval[0]), float(fit_interval[1])),
|
||||
multi_peak_flag=bool(archive[f"{name}_multi_peak_flag"][0]),
|
||||
anchor_times_s=(
|
||||
float(anchor_times[0]),
|
||||
float(anchor_times[1]),
|
||||
float(anchor_times[2]),
|
||||
float(anchor_times[3]),
|
||||
),
|
||||
)
|
||||
counts = {key: archive[key].astype(np.int32) for key in _COUNT_KEYS}
|
||||
return TicCurves(
|
||||
frame_indices=archive["frame_indices"].astype(np.int32),
|
||||
time_video_s=archive["time_video_s"].astype(np.float32),
|
||||
injection_offset_s=float(archive["injection_offset_s"][0]),
|
||||
time_origin=str(archive["time_origin"].item()),
|
||||
fps=float(archive["fps"][0]),
|
||||
config_version=str(archive["config_version"].item()),
|
||||
algorithm_version=str(archive["algorithm_version"].item()),
|
||||
curves_schema_version=schema_version,
|
||||
curves=curves,
|
||||
perf_mask_stable=archive["perf_mask_stable"].astype(bool),
|
||||
reference_valid=bool(archive["reference_valid"][0]),
|
||||
early_truncation_flag=bool(archive["early_truncation_flag"][0]),
|
||||
**counts,
|
||||
)
|
||||
@@ -0,0 +1,928 @@
|
||||
"""第 7 步：TIC 的约束 Lognormal 拟合，以及拟合曲线上的连续时刻 D1–D4。
|
||||
|
||||
正式模型直接拟合清洗后的原始灰度样本。二阶 Savitzky-Golay 平滑只用来找候选峰、圈出第一遍
|
||||
团注（造影剂第一次流过）的时间段和给优化器初值；质量门不通过时，不把候选曲线顶替成正式的
|
||||
`processed` 曲线（`processed` 整条为 NaN），免得一条没拟合好的曲线冒充正式结果。
|
||||
|
||||
照搬旧项目 `tic_fitting.py`，去掉产品层；旧 `_empty_result` 与 `fit_lognormal_curve` 的参数
|
||||
默认值改成调用处显式传入。
|
||||
"""
|
||||
|
||||
import itertools
|
||||
import math
|
||||
from dataclasses import dataclass
|
||||
from typing import Any, cast
|
||||
|
||||
import numpy as np
|
||||
from scipy.optimize import ( # pyright: ignore[reportMissingTypeStubs]
|
||||
brentq, # pyright: ignore[reportUnknownVariableType]
|
||||
least_squares, # pyright: ignore[reportUnknownVariableType]
|
||||
)
|
||||
from scipy.signal import ( # pyright: ignore[reportMissingTypeStubs]
|
||||
find_peaks, # pyright: ignore[reportUnknownVariableType]
|
||||
savgol_filter, # pyright: ignore[reportUnknownVariableType]
|
||||
)
|
||||
|
||||
from cadx.ceus._config import AnchorFramesConfig, TicFittingConfig
|
||||
from cadx.ceus._geometry import AnyArray
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class LognormalParameters:
|
||||
"""Lognormal 团注曲线的五个参数：面积、对数均值、对数标准差、到达偏移（秒）、基线。"""
|
||||
|
||||
area: float
|
||||
mu: float
|
||||
sigma: float
|
||||
arrival_offset_s: float
|
||||
baseline: float
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class LognormalFitResult:
|
||||
"""一条曲线的拟合结果。
|
||||
|
||||
`candidate` 是候选平滑曲线（只作检测用）；`fitted` 是优化器给出的曲线（可能没过质量门）；
|
||||
`processed` 只在 `success` 时等于 `fitted`，否则整条 NaN。`d1_time_s`–`d4_time_s` 是拟合
|
||||
曲线上的到达、洗入拐点、峰值、洗出拐点时刻（视频时间，秒），失败时为 None。
|
||||
"""
|
||||
|
||||
success: bool
|
||||
failure_reason: str
|
||||
algorithm_version: str
|
||||
curve_model: str
|
||||
smoothing_method: str
|
||||
baseline_source: str
|
||||
candidate: AnyArray
|
||||
fitted: AnyArray
|
||||
processed: AnyArray
|
||||
valid_sample_mask: AnyArray
|
||||
duplicate_sample_mask: AnyArray
|
||||
parameters: LognormalParameters | None
|
||||
r2: float
|
||||
nrmse: float
|
||||
residual_noise_ratio: float
|
||||
fit_start_s: float
|
||||
fit_end_s: float
|
||||
multi_peak_flag: bool
|
||||
d1_time_s: float | None
|
||||
d2_time_s: float | None
|
||||
d3_time_s: float | None
|
||||
d4_time_s: float | None
|
||||
|
||||
|
||||
def evaluate_lognormal(time_s: AnyArray | float, parameters: LognormalParameters) -> AnyArray:
|
||||
"""带固定基线与到达偏移的 Lognormal 团注曲线在给定时刻的值。"""
|
||||
times = np.asarray(time_s, dtype=np.float64)
|
||||
values = np.full(times.shape, parameters.baseline, dtype=np.float64)
|
||||
shifted = times - parameters.arrival_offset_s
|
||||
positive = shifted > 0.0
|
||||
if not positive.any():
|
||||
return values
|
||||
x_values = shifted[positive]
|
||||
exponent = -((np.log(x_values) - parameters.mu) ** 2 / (2.0 * parameters.sigma**2))
|
||||
values[positive] += (
|
||||
parameters.area
|
||||
/ (x_values * parameters.sigma * math.sqrt(2.0 * math.pi))
|
||||
* np.exp(exponent)
|
||||
)
|
||||
return values
|
||||
|
||||
|
||||
def _empty_result(
|
||||
size: int,
|
||||
cfg: TicFittingConfig,
|
||||
reason: str,
|
||||
*,
|
||||
candidate: AnyArray | None,
|
||||
fitted: AnyArray | None,
|
||||
valid_mask: AnyArray | None,
|
||||
duplicate_mask: AnyArray | None,
|
||||
parameters: LognormalParameters | None,
|
||||
r2: float,
|
||||
nrmse: float,
|
||||
residual_noise_ratio: float,
|
||||
fit_start_s: float,
|
||||
fit_end_s: float,
|
||||
multi_peak_flag: bool,
|
||||
baseline_source: str,
|
||||
) -> LognormalFitResult:
|
||||
"""拟合被拒时的结果：`success=False`、`processed` 全 NaN、四个时刻为 None。"""
|
||||
nan_curve = np.full(size, np.nan, dtype=np.float64)
|
||||
return LognormalFitResult(
|
||||
success=False,
|
||||
failure_reason=reason,
|
||||
algorithm_version=cfg.algorithm_version,
|
||||
curve_model="lognormal",
|
||||
smoothing_method="",
|
||||
baseline_source=baseline_source,
|
||||
candidate=nan_curve.copy() if candidate is None else candidate,
|
||||
fitted=nan_curve.copy() if fitted is None else fitted,
|
||||
processed=nan_curve.copy(),
|
||||
valid_sample_mask=(np.zeros(size, dtype=bool) if valid_mask is None else valid_mask),
|
||||
duplicate_sample_mask=(
|
||||
np.zeros(size, dtype=bool) if duplicate_mask is None else duplicate_mask
|
||||
),
|
||||
parameters=parameters,
|
||||
r2=r2,
|
||||
nrmse=nrmse,
|
||||
residual_noise_ratio=residual_noise_ratio,
|
||||
fit_start_s=fit_start_s,
|
||||
fit_end_s=fit_end_s,
|
||||
multi_peak_flag=multi_peak_flag,
|
||||
d1_time_s=None,
|
||||
d2_time_s=None,
|
||||
d3_time_s=None,
|
||||
d4_time_s=None,
|
||||
)
|
||||
|
||||
|
||||
def _finite_runs(mask: AnyArray) -> list[tuple[int, int]]:
|
||||
"""布尔序列里连续为真的段，每段 `[start, stop)`。"""
|
||||
runs: list[tuple[int, int]] = []
|
||||
start: int | None = None
|
||||
for index, value in enumerate(mask):
|
||||
if value and start is None:
|
||||
start = index
|
||||
elif not value and start is not None:
|
||||
runs.append((start, index))
|
||||
start = None
|
||||
if start is not None:
|
||||
runs.append((start, len(mask)))
|
||||
return runs
|
||||
|
||||
|
||||
def _odd_window(sample_count: int, desired: int, polynomial_order: int) -> int:
|
||||
"""Savitzky-Golay 窗口：奇数，不超过样本数，不小于多项式阶数 + 2。"""
|
||||
maximum = sample_count if sample_count % 2 == 1 else sample_count - 1
|
||||
minimum = polynomial_order + 2
|
||||
if minimum % 2 == 0:
|
||||
minimum += 1
|
||||
window = desired if desired % 2 == 1 else desired + 1
|
||||
return min(maximum, max(minimum, window))
|
||||
|
||||
|
||||
def _savgol(segment: AnyArray, window: int, polyorder: int) -> AnyArray:
|
||||
"""二阶 Savitzky-Golay 平滑（`mode="interp"`）。"""
|
||||
return cast(
|
||||
AnyArray,
|
||||
savgol_filter(segment, window_length=window, polyorder=polyorder, mode="interp"),
|
||||
)
|
||||
|
||||
|
||||
def _candidate_curve(
|
||||
times: AnyArray,
|
||||
values: AnyArray,
|
||||
valid_mask: AnyArray,
|
||||
cfg: TicFittingConfig,
|
||||
) -> tuple[AnyArray, AnyArray, AnyArray, float]:
|
||||
"""把有效样本插值到等间隔网格，短窗与长窗两次平滑按 1:3 混合，得到候选曲线。
|
||||
|
||||
相邻有效样本间隔超过 `max_candidate_gap_s` 的地方不插值（记 NaN），分段各自平滑。
|
||||
返回（采样时刻上的候选值、网格、网格上的候选值、网格步长）。
|
||||
"""
|
||||
valid_times = times[valid_mask]
|
||||
intervals = np.diff(valid_times)
|
||||
positive_intervals = intervals[intervals > 0.0]
|
||||
if positive_intervals.size == 0:
|
||||
raise ValueError("有效时间样本不足以建立候选网格")
|
||||
grid_step = float(np.median(positive_intervals))
|
||||
grid = np.arange(valid_times[0], valid_times[-1] + grid_step * 0.5, grid_step, dtype=np.float64)
|
||||
interpolated = np.interp(grid, valid_times, values[valid_mask])
|
||||
|
||||
for left, right in itertools.pairwise(valid_times):
|
||||
if right - left - grid_step > cfg.max_candidate_gap_s:
|
||||
interpolated[(grid > left) & (grid < right)] = np.nan
|
||||
|
||||
smoothed = np.full(grid.shape, np.nan, dtype=np.float64)
|
||||
finite = np.isfinite(interpolated)
|
||||
short_window = max(1, round(cfg.candidate_window_s / grid_step))
|
||||
long_window = max(1, round(cfg.candidate_long_window_s / grid_step))
|
||||
for start, stop in _finite_runs(finite):
|
||||
segment = interpolated[start:stop]
|
||||
if segment.size <= cfg.candidate_polyorder + 1:
|
||||
continue
|
||||
short = _odd_window(int(segment.size), short_window, cfg.candidate_polyorder)
|
||||
long = _odd_window(int(segment.size), long_window, cfg.candidate_polyorder)
|
||||
if short <= cfg.candidate_polyorder or long <= cfg.candidate_polyorder:
|
||||
continue
|
||||
short_values = _savgol(segment, short, cfg.candidate_polyorder)
|
||||
long_values = _savgol(segment, long, cfg.candidate_polyorder)
|
||||
smoothed[start:stop] = 0.25 * short_values + 0.75 * long_values
|
||||
|
||||
at_samples = np.full(times.shape, np.nan, dtype=np.float64)
|
||||
for start, stop in _finite_runs(np.isfinite(smoothed)):
|
||||
in_run = (times >= grid[start]) & (times <= grid[stop - 1])
|
||||
at_samples[in_run] = np.interp(times[in_run], grid[start:stop], smoothed[start:stop])
|
||||
return at_samples, grid, smoothed, grid_step
|
||||
|
||||
|
||||
def _candidate_peaks(
|
||||
grid: AnyArray,
|
||||
candidate: AnyArray,
|
||||
baseline: float,
|
||||
baseline_noise: float,
|
||||
baseline_end_s: float,
|
||||
grid_step: float,
|
||||
cfg: TicFittingConfig,
|
||||
anchor_cfg: AnchorFramesConfig,
|
||||
) -> list[int]:
|
||||
"""基线段之后的候选峰（网格下标）：高度与突出度都须过「峰高比例」与「基线抖动倍数」的高者。"""
|
||||
enhancement = candidate - baseline
|
||||
eligible = np.isfinite(enhancement) & (grid >= baseline_end_s)
|
||||
if not eligible.any():
|
||||
return []
|
||||
amplitude = float(np.nanmax(enhancement[eligible]))
|
||||
if not math.isfinite(amplitude) or amplitude <= 0.0:
|
||||
return []
|
||||
minimum_height = max(
|
||||
cfg.minimum_relative_peak_height * amplitude,
|
||||
anchor_cfg.peak_prominence_baseline_sd_multiplier * baseline_noise,
|
||||
)
|
||||
distance = max(1, math.ceil(cfg.minimum_peak_distance_s / grid_step))
|
||||
peaks: list[int] = []
|
||||
for start, stop in _finite_runs(eligible):
|
||||
segment = enhancement[start:stop]
|
||||
if segment.size < 3:
|
||||
continue
|
||||
found = cast(
|
||||
tuple[AnyArray, object],
|
||||
find_peaks(
|
||||
segment, height=minimum_height, prominence=minimum_height, distance=distance
|
||||
),
|
||||
)
|
||||
local = found[0]
|
||||
peaks.extend(int(value) + start for value in local.tolist())
|
||||
return sorted(peaks)
|
||||
|
||||
|
||||
def _estimate_baseline(
|
||||
values: AnyArray,
|
||||
valid_mask: AnyArray,
|
||||
cfg: TicFittingConfig,
|
||||
) -> tuple[float, float]:
|
||||
"""视频开头没有稳定基线时，用低分位数估一个有界的基线和噪声。"""
|
||||
finite = values[valid_mask]
|
||||
baseline = float(np.percentile(finite, cfg.estimated_baseline_percentile))
|
||||
low_limit = float(np.percentile(finite, min(cfg.estimated_baseline_percentile * 2.0, 45.0)))
|
||||
low_values = finite[finite <= low_limit]
|
||||
center = float(np.median(low_values)) if low_values.size else baseline
|
||||
noise = 1.4826 * float(np.median(np.abs(low_values - center)))
|
||||
if not math.isfinite(noise) or noise <= 0.0:
|
||||
differences = np.diff(finite)
|
||||
difference_center = float(np.median(differences)) if differences.size else 0.0
|
||||
noise = 1.4826 * float(np.median(np.abs(differences - difference_center)))
|
||||
amplitude = max(float(np.percentile(finite, 95.0) - baseline), 1e-3)
|
||||
return baseline, max(noise, amplitude * 0.01, 1e-3)
|
||||
|
||||
|
||||
def _candidate_detection_noise(
|
||||
grid: AnyArray,
|
||||
candidate: AnyArray,
|
||||
baseline_end_s: float,
|
||||
baseline_noise: float,
|
||||
) -> float:
|
||||
"""检测用的噪声：基线段候选曲线一阶差分的稳健离散度，不超过基线噪声。"""
|
||||
baseline_values = candidate[np.isfinite(candidate) & (grid <= baseline_end_s)]
|
||||
differences = np.diff(baseline_values)
|
||||
if differences.size < 2:
|
||||
return baseline_noise
|
||||
center = float(np.median(differences))
|
||||
difference_noise = 1.4826 * float(np.median(np.abs(differences - center)))
|
||||
return min(baseline_noise, max(difference_noise / math.sqrt(2.0), 1e-3))
|
||||
|
||||
|
||||
def _select_primary_peak(
|
||||
peaks: list[int],
|
||||
grid: AnyArray,
|
||||
candidate: AnyArray,
|
||||
baseline: float,
|
||||
) -> int | None:
|
||||
"""主峰：按「升降完整度 × 半高宽」打分，够强（≥ 最强者 1/4）的里取最早的一个。"""
|
||||
if not peaks:
|
||||
return None
|
||||
finite_runs = _finite_runs(np.isfinite(candidate))
|
||||
scored: list[tuple[float, float, int]] = []
|
||||
for peak in peaks:
|
||||
containing_run = next(
|
||||
((start, stop) for start, stop in finite_runs if start < peak < stop - 1),
|
||||
None,
|
||||
)
|
||||
if containing_run is None:
|
||||
continue
|
||||
first, stop = containing_run
|
||||
last = stop - 1
|
||||
peak_height = float(candidate[peak] - baseline)
|
||||
left = candidate[first : peak + 1]
|
||||
right = candidate[peak : last + 1]
|
||||
if peak_height <= 0.0:
|
||||
continue
|
||||
rise = float(candidate[peak] - np.min(left))
|
||||
decline = float(candidate[peak] - np.min(right))
|
||||
completeness = min(rise, decline)
|
||||
half_prominence = float(candidate[peak] - 0.5 * completeness)
|
||||
width_left = peak
|
||||
while width_left > first and candidate[width_left - 1] >= half_prominence:
|
||||
width_left -= 1
|
||||
width_right = peak
|
||||
while width_right < last and candidate[width_right + 1] >= half_prominence:
|
||||
width_right += 1
|
||||
width = float(grid[width_right] - grid[width_left])
|
||||
if width <= 0.0:
|
||||
continue
|
||||
scored.append((completeness * max(width, 1e-6), peak_height, peak))
|
||||
if not scored:
|
||||
return None
|
||||
complete = [item for item in scored if item[0] > 0.0 and item[1] > 0.0]
|
||||
if complete:
|
||||
strongest = max(item[0] for item in complete)
|
||||
credible = [item for item in complete if item[0] >= strongest * 0.25]
|
||||
if credible:
|
||||
return min(credible, key=lambda item: item[2])[2]
|
||||
return max(scored, key=lambda item: (item[0], item[1], -item[2]))[2]
|
||||
|
||||
|
||||
def _persistent_arrival(
|
||||
grid: AnyArray,
|
||||
enhancement: AnyArray,
|
||||
peak_index: int,
|
||||
threshold: float,
|
||||
persistence_s: float,
|
||||
) -> float | None:
|
||||
"""主峰之前第一次越过门槛并持续 `persistence_s` 秒的时刻；没有返回 None。"""
|
||||
above = enhancement >= threshold
|
||||
run_start: int | None = None
|
||||
for index in range(peak_index + 1):
|
||||
if not np.isfinite(enhancement[index]) or not above[index]:
|
||||
run_start = None
|
||||
continue
|
||||
if run_start is None:
|
||||
run_start = index
|
||||
if grid[index] - grid[run_start] >= persistence_s:
|
||||
return float(grid[run_start])
|
||||
return None
|
||||
|
||||
|
||||
def _first_bolus_end(
|
||||
grid: AnyArray,
|
||||
candidate: AnyArray,
|
||||
baseline: float,
|
||||
peaks: list[int],
|
||||
primary_peak: int,
|
||||
cfg: TicFittingConfig,
|
||||
) -> tuple[float, bool]:
|
||||
"""第一遍团注到哪里结束：只把相隔足够远、谷值足够深的后峰视为第二遍团注。"""
|
||||
first_peak = primary_peak
|
||||
first_height = float(candidate[first_peak] - baseline)
|
||||
for second_peak in (peak for peak in peaks if peak > first_peak):
|
||||
if grid[second_peak] - grid[first_peak] < cfg.multi_peak_minimum_separation_s:
|
||||
continue
|
||||
second_height = float(candidate[second_peak] - baseline)
|
||||
if second_height < cfg.multi_peak_minimum_relative_height * first_height:
|
||||
continue
|
||||
valley_slice = candidate[first_peak : second_peak + 1]
|
||||
if not np.isfinite(valley_slice).all() or valley_slice.size < 3:
|
||||
continue
|
||||
valley_offset = int(np.argmin(valley_slice))
|
||||
valley_height = max(float(valley_slice[valley_offset] - baseline), 0.0)
|
||||
reference_height = max(min(first_height, second_height), float(np.finfo(float).eps))
|
||||
if valley_height / reference_height <= cfg.multi_peak_maximum_valley_fraction:
|
||||
return float(grid[first_peak + valley_offset]), True
|
||||
return float(grid[-1]), False
|
||||
|
||||
|
||||
def _unpack_optimizer_parameters(packed: AnyArray) -> LognormalParameters:
|
||||
"""优化器的参数化（对数峰高、峰时、对数 σ、到达偏移、基线）→ Lognormal 五参数。"""
|
||||
peak_height = float(math.exp(float(packed[0])))
|
||||
peak_time_s = float(packed[1])
|
||||
sigma = float(math.exp(float(packed[2])))
|
||||
arrival_offset_s = float(packed[3])
|
||||
baseline = float(packed[4])
|
||||
mode_x = peak_time_s - arrival_offset_s
|
||||
mu = math.log(mode_x) + sigma**2
|
||||
area = peak_height * mode_x * sigma * math.sqrt(2.0 * math.pi) * math.exp(0.5 * sigma**2)
|
||||
return LognormalParameters(
|
||||
area=area,
|
||||
mu=mu,
|
||||
sigma=sigma,
|
||||
arrival_offset_s=arrival_offset_s,
|
||||
baseline=baseline,
|
||||
)
|
||||
|
||||
|
||||
def _initial_guesses(
|
||||
peak_time_s: float,
|
||||
peak_height: float,
|
||||
t0_lower: float,
|
||||
t0_upper: float,
|
||||
baseline: float,
|
||||
sigma_min: float,
|
||||
sigma_max: float,
|
||||
) -> list[AnyArray]:
|
||||
"""多组初值：到达偏移取区间两端与中点，σ 取四档，峰时取候选峰 ±0.25 秒。"""
|
||||
guesses: list[AnyArray] = []
|
||||
t0_values = np.unique(
|
||||
np.clip(
|
||||
np.asarray([t0_lower, (t0_lower + t0_upper) / 2.0, t0_upper], dtype=np.float64),
|
||||
t0_lower,
|
||||
t0_upper,
|
||||
)
|
||||
)
|
||||
peak_time_values = np.asarray(
|
||||
[peak_time_s - 0.25, peak_time_s, peak_time_s + 0.25], dtype=np.float64
|
||||
)
|
||||
sigma_values = np.clip(
|
||||
np.asarray([0.25, 0.5, 0.9, 1.3], dtype=np.float64),
|
||||
sigma_min,
|
||||
sigma_max,
|
||||
)
|
||||
for t0_value in t0_values.tolist():
|
||||
for sigma in np.unique(sigma_values).tolist():
|
||||
for peak_time_value in peak_time_values.tolist():
|
||||
if peak_time_value <= t0_value:
|
||||
continue
|
||||
guesses.append(
|
||||
np.asarray(
|
||||
[
|
||||
math.log(max(peak_height, 1e-8)),
|
||||
peak_time_value,
|
||||
math.log(float(sigma)),
|
||||
t0_value,
|
||||
baseline,
|
||||
],
|
||||
dtype=np.float64,
|
||||
)
|
||||
)
|
||||
return guesses
|
||||
|
||||
|
||||
def _optimizer_bounds(
|
||||
values: AnyArray,
|
||||
baseline: float,
|
||||
baseline_noise: float,
|
||||
peak_height: float,
|
||||
t0_lower: float,
|
||||
t0_upper: float,
|
||||
peak_lower: float,
|
||||
peak_upper: float,
|
||||
baseline_observed: bool,
|
||||
cfg: TicFittingConfig,
|
||||
) -> tuple[AnyArray, AnyArray, float]:
|
||||
"""优化器参数的上下界和基线惩罚的尺度；观测到基线时基线只许在其附近小幅浮动。"""
|
||||
amplitude = max(float(np.max(values) - baseline), peak_height, 1e-3)
|
||||
height_upper = max(amplitude * 10.0, 1.0)
|
||||
baseline_span = max(baseline_noise * 3.0, amplitude * 0.05, 1e-3)
|
||||
if baseline_observed:
|
||||
baseline_lower = baseline - baseline_span
|
||||
baseline_upper = baseline + baseline_span
|
||||
baseline_scale = baseline_span
|
||||
else:
|
||||
baseline_lower = max(0.0, float(np.percentile(values, 1.0)) - amplitude * 0.25)
|
||||
baseline_upper = min(float(np.percentile(values, 45.0)), baseline + amplitude * 0.35)
|
||||
baseline_upper = max(baseline_upper, baseline_lower + 1e-3)
|
||||
baseline_scale = max(amplitude * 0.2, baseline_noise, 1e-3)
|
||||
lower = np.asarray(
|
||||
[
|
||||
math.log(1e-8),
|
||||
peak_lower,
|
||||
math.log(cfg.sigma_min),
|
||||
t0_lower,
|
||||
baseline_lower,
|
||||
],
|
||||
dtype=np.float64,
|
||||
)
|
||||
upper = np.asarray(
|
||||
[
|
||||
math.log(height_upper),
|
||||
peak_upper,
|
||||
math.log(cfg.sigma_max),
|
||||
t0_upper,
|
||||
baseline_upper,
|
||||
],
|
||||
dtype=np.float64,
|
||||
)
|
||||
return lower, upper, baseline_scale
|
||||
|
||||
|
||||
def _fit_parameters(
|
||||
times: AnyArray,
|
||||
values: AnyArray,
|
||||
baseline: float,
|
||||
baseline_noise: float,
|
||||
peak_time_s: float,
|
||||
peak_height: float,
|
||||
t0_lower: float,
|
||||
t0_upper: float,
|
||||
peak_lower: float,
|
||||
peak_upper: float,
|
||||
baseline_observed: bool,
|
||||
cfg: TicFittingConfig,
|
||||
) -> LognormalParameters | None:
|
||||
"""从每组初值各跑一次有界稳健最小二乘，取平方误差最小的那组；全部失败返回 None。"""
|
||||
lower, upper, baseline_scale = _optimizer_bounds(
|
||||
values,
|
||||
baseline,
|
||||
baseline_noise,
|
||||
peak_height,
|
||||
t0_lower,
|
||||
t0_upper,
|
||||
peak_lower,
|
||||
peak_upper,
|
||||
baseline_observed,
|
||||
cfg,
|
||||
)
|
||||
amplitude = max(float(np.max(values) - baseline), peak_height, 1e-3)
|
||||
residual_scale = max(baseline_noise, amplitude * 0.05, 1e-3)
|
||||
|
||||
def residuals(packed: AnyArray) -> AnyArray:
|
||||
if packed[1] <= packed[3]:
|
||||
return np.full(values.size + 1, 1e6, dtype=np.float64)
|
||||
parameters = _unpack_optimizer_parameters(packed)
|
||||
data = (evaluate_lognormal(times, parameters) - values) / residual_scale
|
||||
penalty = np.asarray([(parameters.baseline - baseline) / baseline_scale])
|
||||
return np.concatenate((data, penalty))
|
||||
|
||||
best: tuple[float, LognormalParameters] | None = None
|
||||
for guess in _initial_guesses(
|
||||
peak_time_s,
|
||||
peak_height,
|
||||
t0_lower,
|
||||
t0_upper,
|
||||
baseline,
|
||||
cfg.sigma_min,
|
||||
cfg.sigma_max,
|
||||
):
|
||||
clipped = np.clip(guess, lower + 1e-9, upper - 1e-9)
|
||||
result = cast(
|
||||
Any,
|
||||
least_squares(
|
||||
residuals,
|
||||
clipped,
|
||||
bounds=(lower, upper),
|
||||
loss=cfg.robust_loss,
|
||||
f_scale=1.0,
|
||||
max_nfev=cfg.max_nfev,
|
||||
),
|
||||
)
|
||||
solution = np.asarray(cast(AnyArray, result.x), dtype=np.float64)
|
||||
if not bool(result.success) or not np.isfinite(solution).all():
|
||||
continue
|
||||
if solution[1] <= solution[3]:
|
||||
continue
|
||||
parameters = _unpack_optimizer_parameters(solution)
|
||||
squared_error = float(np.sum((evaluate_lognormal(times, parameters) - values) ** 2))
|
||||
if best is None or squared_error < best[0]:
|
||||
best = (squared_error, parameters)
|
||||
return None if best is None else best[1]
|
||||
|
||||
|
||||
def _parameters_on_boundary(
|
||||
parameters: LognormalParameters,
|
||||
lower: AnyArray,
|
||||
upper: AnyArray,
|
||||
fraction: float,
|
||||
) -> bool:
|
||||
"""拟合参数（换回优化器参数化后）是否贴着上下界（区间宽度的 `fraction` 以内）。"""
|
||||
peak_time_s = parameters.arrival_offset_s + math.exp(parameters.mu - parameters.sigma**2)
|
||||
peak_height = float(
|
||||
evaluate_lognormal(np.asarray([peak_time_s]), parameters)[0] - parameters.baseline
|
||||
)
|
||||
packed = np.asarray(
|
||||
[
|
||||
math.log(max(peak_height, 1e-8)),
|
||||
peak_time_s,
|
||||
math.log(parameters.sigma),
|
||||
parameters.arrival_offset_s,
|
||||
parameters.baseline,
|
||||
],
|
||||
dtype=np.float64,
|
||||
)
|
||||
margin = (upper - lower) * fraction
|
||||
return bool(np.any(packed <= lower + margin) or np.any(packed >= upper - margin))
|
||||
|
||||
|
||||
def _quality_metrics(
|
||||
observed: AnyArray,
|
||||
fitted: AnyArray,
|
||||
baseline_noise: float,
|
||||
) -> tuple[float, float, float]:
|
||||
"""拟合优度：R²、按观测幅度归一的 RMSE、RMSE 与基线噪声之比。"""
|
||||
residual = observed - fitted
|
||||
residual_sum = float(np.sum(residual**2))
|
||||
centered_sum = float(np.sum((observed - float(np.mean(observed))) ** 2))
|
||||
r2 = float("nan") if centered_sum <= 0.0 else 1.0 - residual_sum / centered_sum
|
||||
rmse = math.sqrt(residual_sum / max(int(observed.size), 1))
|
||||
amplitude = max(float(np.max(observed) - np.min(observed)), 1e-6)
|
||||
nrmse = rmse / amplitude
|
||||
residual_noise_ratio = rmse / max(baseline_noise, 1e-6)
|
||||
return r2, nrmse, residual_noise_ratio
|
||||
|
||||
|
||||
def evaluate_lognormal_derivative(time_s: float, parameters: LognormalParameters) -> float:
|
||||
"""Lognormal 曲线在某一时刻的一阶导数（到达偏移之前为 0）。"""
|
||||
shifted = time_s - parameters.arrival_offset_s
|
||||
if shifted <= 0.0:
|
||||
return 0.0
|
||||
value = float(evaluate_lognormal(np.asarray([time_s]), parameters)[0])
|
||||
pulse = value - parameters.baseline
|
||||
return (
|
||||
pulse
|
||||
* (parameters.mu - parameters.sigma**2 - math.log(shifted))
|
||||
/ (parameters.sigma**2 * shifted)
|
||||
)
|
||||
|
||||
|
||||
def _arrival_root(
|
||||
parameters: LognormalParameters, threshold: float, left: float, right: float
|
||||
) -> float:
|
||||
"""拟合曲线在 `[left, right]` 内升过「基线 + 门槛」的时刻（Brent 求根）。"""
|
||||
|
||||
def above_threshold(value: float) -> float:
|
||||
return float(
|
||||
evaluate_lognormal(np.asarray([value]), parameters)[0] - parameters.baseline - threshold
|
||||
)
|
||||
|
||||
return float(cast(float, brentq(above_threshold, left, right)))
|
||||
|
||||
|
||||
def _continuous_anchors(
|
||||
parameters: LognormalParameters,
|
||||
baseline_noise: float,
|
||||
observed_end_s: float,
|
||||
anchor_cfg: AnchorFramesConfig,
|
||||
boundary_margin_s: float,
|
||||
) -> tuple[float | None, float | None, float, float | None]:
|
||||
"""拟合曲线上的四个连续时刻：到达 D1、洗入拐点 D2、峰值 D3、洗出拐点 D4。
|
||||
|
||||
到达门槛取「基线抖动 × 倍数」与「峰高 × 比例」的高者（`anchor_frames.arrival`）。D4 只在
|
||||
峰后还有足够观测（`boundary_margin_s`）、且拐点落在观测范围内、导数为负时给出。
|
||||
"""
|
||||
peak_time = parameters.arrival_offset_s + math.exp(parameters.mu - parameters.sigma**2)
|
||||
peak_value = float(evaluate_lognormal(np.asarray([peak_time]), parameters)[0])
|
||||
peak_height = peak_value - parameters.baseline
|
||||
threshold = max(
|
||||
anchor_cfg.arrival.baseline_sd_multiplier * baseline_noise,
|
||||
anchor_cfg.arrival.relative_peak_threshold * peak_height,
|
||||
)
|
||||
|
||||
left = parameters.arrival_offset_s + 1e-9
|
||||
arrival: float | None = None
|
||||
if threshold > 0.0 and peak_height >= threshold:
|
||||
arrival = _arrival_root(parameters, threshold, left, peak_time)
|
||||
|
||||
sigma_squared = parameters.sigma**2
|
||||
extrema_discriminant = math.sqrt(sigma_squared**2 + 4.0 * sigma_squared)
|
||||
wash_in_factor = (sigma_squared + extrema_discriminant) / 2.0
|
||||
wash_out_factor = (sigma_squared - extrema_discriminant) / 2.0
|
||||
wash_in_candidate = parameters.arrival_offset_s + math.exp(
|
||||
parameters.mu - sigma_squared - wash_in_factor
|
||||
)
|
||||
wash_in_start = arrival if arrival is not None else left
|
||||
d2_time = max(wash_in_start, wash_in_candidate)
|
||||
|
||||
d4_time: float | None = None
|
||||
if observed_end_s - peak_time >= boundary_margin_s:
|
||||
candidate_d4 = parameters.arrival_offset_s + math.exp(
|
||||
parameters.mu - sigma_squared - wash_out_factor
|
||||
)
|
||||
if (
|
||||
candidate_d4 < observed_end_s - boundary_margin_s
|
||||
and evaluate_lognormal_derivative(candidate_d4, parameters) < 0.0
|
||||
):
|
||||
d4_time = candidate_d4
|
||||
return arrival, d2_time, peak_time, d4_time
|
||||
|
||||
|
||||
def fit_lognormal_curve(
|
||||
time_s: AnyArray,
|
||||
raw_values: AnyArray,
|
||||
*,
|
||||
baseline: float | None,
|
||||
baseline_noise: float | None,
|
||||
baseline_end_index: int,
|
||||
duplicate_mask: AnyArray,
|
||||
sample_mask: AnyArray | None,
|
||||
cfg: TicFittingConfig,
|
||||
anchor_cfg: AnchorFramesConfig,
|
||||
) -> LognormalFitResult:
|
||||
"""拟合一条原始灰度 TIC，并给出连续的 D1–D4 时刻。
|
||||
|
||||
`baseline`/`baseline_noise` 为 None 或 `baseline_end_index <= 0` 表示视频开头没有观测到
|
||||
稳定基线，改用低分位数估计（`baseline_source="estimated"`）。`sample_mask` 为 None 表示
|
||||
全部样本可用。任何一道门不过都返回 `success=False`，原因写在 `failure_reason`。
|
||||
"""
|
||||
times = np.asarray(time_s, dtype=np.float64)
|
||||
values = np.asarray(raw_values, dtype=np.float64)
|
||||
duplicates = np.asarray(duplicate_mask, dtype=bool)
|
||||
if times.ndim != 1 or values.shape != times.shape or duplicates.shape != times.shape:
|
||||
raise ValueError("time_s、raw_values 与 duplicate_mask 必须是等长一维数组")
|
||||
if times.size == 0 or not np.isfinite(times).all() or np.any(np.diff(times) <= 0.0):
|
||||
raise ValueError("time_s 必须是非空、有限且严格递增的一维数组")
|
||||
eligibility = (
|
||||
np.ones(times.shape, dtype=bool)
|
||||
if sample_mask is None
|
||||
else np.asarray(sample_mask, dtype=bool)
|
||||
)
|
||||
if eligibility.shape != times.shape:
|
||||
raise ValueError("sample_mask 必须与 time_s 等长")
|
||||
valid_mask = np.isfinite(values) & ~duplicates & eligibility
|
||||
|
||||
candidate_at_samples: AnyArray | None = None
|
||||
fitted: AnyArray | None = None
|
||||
parameters: LognormalParameters | None = None
|
||||
r2 = nrmse = residual_noise_ratio = float("nan")
|
||||
fit_start_s = fit_end_s = float("nan")
|
||||
multi_peak = False
|
||||
baseline_source = "unavailable"
|
||||
|
||||
def reject(reason: str) -> LognormalFitResult:
|
||||
return _empty_result(
|
||||
len(times),
|
||||
cfg,
|
||||
reason,
|
||||
candidate=candidate_at_samples,
|
||||
fitted=fitted,
|
||||
valid_mask=valid_mask,
|
||||
duplicate_mask=duplicates,
|
||||
parameters=parameters,
|
||||
r2=r2,
|
||||
nrmse=nrmse,
|
||||
residual_noise_ratio=residual_noise_ratio,
|
||||
fit_start_s=fit_start_s,
|
||||
fit_end_s=fit_end_s,
|
||||
multi_peak_flag=multi_peak,
|
||||
baseline_source=baseline_source,
|
||||
)
|
||||
|
||||
if int(valid_mask.sum()) < cfg.minimum_valid_samples:
|
||||
return reject("insufficient_valid_samples")
|
||||
|
||||
# 有效样本时刻全部相同时（理论上严格递增已排除）建不出候选网格，记为拟合失败。
|
||||
try:
|
||||
candidate_at_samples, grid, candidate_grid, grid_step = _candidate_curve(
|
||||
times, values, valid_mask, cfg
|
||||
)
|
||||
except ValueError:
|
||||
return reject("candidate_grid_unavailable")
|
||||
|
||||
baseline_observed = (
|
||||
baseline is not None and baseline_noise is not None and baseline_end_index > 0
|
||||
)
|
||||
if baseline is not None and baseline_noise is not None and baseline_end_index > 0:
|
||||
fit_baseline = float(baseline)
|
||||
fit_baseline_noise = max(float(baseline_noise), 1e-3)
|
||||
baseline_end_s = float(times[min(baseline_end_index - 1, len(times) - 1)])
|
||||
baseline_source = "observed"
|
||||
else:
|
||||
fit_baseline, fit_baseline_noise = _estimate_baseline(values, valid_mask, cfg)
|
||||
baseline_end_s = float(times[valid_mask][0])
|
||||
baseline_source = "estimated"
|
||||
detection_noise = _candidate_detection_noise(
|
||||
grid, candidate_grid, baseline_end_s, fit_baseline_noise
|
||||
)
|
||||
peaks = _candidate_peaks(
|
||||
grid,
|
||||
candidate_grid,
|
||||
fit_baseline,
|
||||
detection_noise,
|
||||
baseline_end_s,
|
||||
grid_step,
|
||||
cfg,
|
||||
anchor_cfg,
|
||||
)
|
||||
first_peak = _select_primary_peak(peaks, grid, candidate_grid, fit_baseline)
|
||||
if first_peak is None:
|
||||
return reject("credible_peak_not_found")
|
||||
|
||||
first_peak_height = float(candidate_grid[first_peak] - fit_baseline)
|
||||
arrival_threshold = max(
|
||||
anchor_cfg.arrival.baseline_sd_multiplier * detection_noise,
|
||||
anchor_cfg.arrival.relative_peak_threshold * first_peak_height,
|
||||
)
|
||||
if (
|
||||
_persistent_arrival(
|
||||
grid,
|
||||
candidate_grid - fit_baseline,
|
||||
first_peak,
|
||||
arrival_threshold,
|
||||
anchor_cfg.arrival.persistence_s,
|
||||
)
|
||||
is None
|
||||
):
|
||||
return reject("candidate_arrival_not_found")
|
||||
|
||||
fit_end_s, multi_peak = _first_bolus_end(
|
||||
grid, candidate_grid, fit_baseline, peaks, first_peak, cfg
|
||||
)
|
||||
|
||||
fit_mask = valid_mask & (times <= fit_end_s)
|
||||
fit_times = times[fit_mask]
|
||||
fit_values = values[fit_mask]
|
||||
if fit_times.size < cfg.minimum_valid_samples:
|
||||
return reject("insufficient_bolus_samples")
|
||||
fit_start_s = float(fit_times[0])
|
||||
long_gaps = np.diff(fit_times)
|
||||
pre_peak = fit_times[:-1] < float(grid[first_peak])
|
||||
if np.any(long_gaps[pre_peak] - grid_step > cfg.max_candidate_gap_s):
|
||||
return reject("long_gap_in_wash_in")
|
||||
|
||||
observed_start_s = fit_start_s
|
||||
peak_lower = max(
|
||||
observed_start_s + cfg.boundary_margin_s,
|
||||
float(grid[first_peak]) - cfg.candidate_long_window_s,
|
||||
)
|
||||
peak_upper = min(
|
||||
fit_end_s - cfg.boundary_margin_s,
|
||||
float(grid[first_peak]) + cfg.candidate_long_window_s,
|
||||
)
|
||||
t0_lower = (
|
||||
observed_start_s
|
||||
if baseline_observed
|
||||
else observed_start_s - cfg.estimated_baseline_lookback_s
|
||||
)
|
||||
t0_upper = peak_lower - max(grid_step * 0.25, 1e-4)
|
||||
if t0_upper <= t0_lower or peak_upper <= peak_lower:
|
||||
return reject("arrival_bounds_invalid")
|
||||
|
||||
parameters = _fit_parameters(
|
||||
fit_times,
|
||||
fit_values,
|
||||
fit_baseline,
|
||||
fit_baseline_noise,
|
||||
float(grid[first_peak]),
|
||||
first_peak_height,
|
||||
t0_lower,
|
||||
t0_upper,
|
||||
peak_lower,
|
||||
peak_upper,
|
||||
baseline_observed,
|
||||
cfg,
|
||||
)
|
||||
if parameters is None:
|
||||
return reject("optimizer_failed")
|
||||
|
||||
fitted = evaluate_lognormal(times, parameters)
|
||||
fit_predictions = evaluate_lognormal(fit_times, parameters)
|
||||
r2, nrmse, residual_noise_ratio = _quality_metrics(
|
||||
fit_values, fit_predictions, fit_baseline_noise
|
||||
)
|
||||
peak_time_s = parameters.arrival_offset_s + math.exp(parameters.mu - parameters.sigma**2)
|
||||
lower, upper, _ = _optimizer_bounds(
|
||||
fit_values,
|
||||
fit_baseline,
|
||||
fit_baseline_noise,
|
||||
first_peak_height,
|
||||
t0_lower,
|
||||
t0_upper,
|
||||
peak_lower,
|
||||
peak_upper,
|
||||
baseline_observed,
|
||||
cfg,
|
||||
)
|
||||
if _parameters_on_boundary(parameters, lower, upper, cfg.parameter_boundary_fraction):
|
||||
return reject("parameter_on_boundary")
|
||||
if (
|
||||
not observed_start_s + cfg.boundary_margin_s
|
||||
< peak_time_s
|
||||
< fit_end_s - cfg.boundary_margin_s
|
||||
):
|
||||
return reject("model_peak_on_boundary")
|
||||
residual_supported = (
|
||||
math.isfinite(residual_noise_ratio)
|
||||
and residual_noise_ratio <= cfg.maximum_residual_noise_ratio
|
||||
)
|
||||
variance_supported = math.isfinite(r2) and r2 >= cfg.minimum_r2
|
||||
if nrmse > cfg.maximum_nrmse or not (variance_supported or residual_supported):
|
||||
return reject("fit_quality_failed")
|
||||
|
||||
d1_time, d2_time, d3_time, d4_time = _continuous_anchors(
|
||||
parameters,
|
||||
fit_baseline_noise,
|
||||
fit_end_s,
|
||||
anchor_cfg,
|
||||
cfg.boundary_margin_s,
|
||||
)
|
||||
return LognormalFitResult(
|
||||
success=True,
|
||||
failure_reason="",
|
||||
algorithm_version=cfg.algorithm_version,
|
||||
curve_model="lognormal",
|
||||
smoothing_method="",
|
||||
baseline_source=baseline_source,
|
||||
candidate=candidate_at_samples,
|
||||
fitted=fitted,
|
||||
processed=fitted.copy(),
|
||||
valid_sample_mask=valid_mask,
|
||||
duplicate_sample_mask=duplicates,
|
||||
parameters=parameters,
|
||||
r2=r2,
|
||||
nrmse=nrmse,
|
||||
residual_noise_ratio=residual_noise_ratio,
|
||||
fit_start_s=observed_start_s,
|
||||
fit_end_s=fit_end_s,
|
||||
multi_peak_flag=multi_peak,
|
||||
d1_time_s=d1_time,
|
||||
d2_time_s=d2_time,
|
||||
d3_time_s=d3_time,
|
||||
d4_time_s=d4_time,
|
||||
)
|
||||
@@ -0,0 +1,181 @@
|
||||
# matplotlib 的类型存根把绘图方法的 **kwargs 标成 Unknown，这一个文件里关掉这条规则。
|
||||
# pyright: reportUnknownMemberType=false
|
||||
"""tic.png：TIC 原始样本、过了质量门的 Lognormal 拟合曲线和关键帧时刻竖线。
|
||||
|
||||
紫圈曲线拟合没过质量门时仍画出拟合（虚线、标 QC rejected），便于复核；病灶曲线拟合没过就不画。
|
||||
照搬旧项目 `tic_plot.py`，去掉产品层。
|
||||
"""
|
||||
|
||||
from collections.abc import Sequence
|
||||
from dataclasses import dataclass
|
||||
from io import BytesIO
|
||||
|
||||
import numpy as np
|
||||
from matplotlib.axes import Axes
|
||||
from matplotlib.backends.backend_agg import FigureCanvasAgg
|
||||
from matplotlib.figure import Figure
|
||||
|
||||
from cadx.ceus._geometry import AnyArray
|
||||
from cadx.ceus._tic import TicResult
|
||||
|
||||
_FIGURE_SIZE = (16.0, 8.0)
|
||||
_FIGURE_DPI = 100
|
||||
_CURVE_STYLES = {
|
||||
"roi_all": ("roi_all", "#2563B8"),
|
||||
"roi_perf": ("roi_perf (dark removed)", "#D46700"),
|
||||
"roi_pan": ("roi_pan reference", "#7C3AED"),
|
||||
}
|
||||
|
||||
|
||||
@dataclass(frozen=True)
|
||||
class TicPlotMarker:
|
||||
"""复核时间标记；颜色沿用 OpenCV 的 BGR 顺序。"""
|
||||
|
||||
label: str
|
||||
time_s: float
|
||||
color: tuple[int, int, int]
|
||||
|
||||
|
||||
def render_tic_png(tic: TicResult, *, markers: Sequence[TicPlotMarker]) -> bytes:
|
||||
"""把病灶与紫圈（正常胰腺对照）的 TIC 画成 PNG 字节；`markers` 是要画的时刻竖线。"""
|
||||
times, raw_series, processed_series = _curve_arrays(tic)
|
||||
figure = Figure(
|
||||
figsize=_FIGURE_SIZE,
|
||||
dpi=_FIGURE_DPI,
|
||||
facecolor="white",
|
||||
layout="constrained",
|
||||
)
|
||||
canvas = FigureCanvasAgg(figure)
|
||||
axes = figure.subplots()
|
||||
_style_axes(axes)
|
||||
|
||||
plotted_values: list[AnyArray] = []
|
||||
for name, (label, color) in _CURVE_STYLES.items():
|
||||
values = raw_series[name]
|
||||
if np.isfinite(values).any():
|
||||
axes.plot(
|
||||
times,
|
||||
values,
|
||||
color=color,
|
||||
alpha=0.24,
|
||||
linewidth=1.2,
|
||||
marker=".",
|
||||
markersize=2.2,
|
||||
label=f"{label} raw",
|
||||
)
|
||||
plotted_values.append(values)
|
||||
|
||||
for name, (label, color) in _CURVE_STYLES.items():
|
||||
fit = tic.fits[name]
|
||||
if fit.success:
|
||||
values = processed_series[name]
|
||||
fit_label = f"{label} Lognormal fit"
|
||||
linestyle = "-"
|
||||
elif name == "roi_pan":
|
||||
values = np.asarray(fit.fitted, dtype=np.float64)
|
||||
fit_label = f"{label} Lognormal fit [QC rejected]"
|
||||
linestyle = "--"
|
||||
else:
|
||||
continue
|
||||
if np.isfinite(values).any():
|
||||
axes.plot(
|
||||
times,
|
||||
values,
|
||||
color=color,
|
||||
linewidth=3.0,
|
||||
linestyle=linestyle,
|
||||
label=fit_label,
|
||||
)
|
||||
plotted_values.append(values)
|
||||
|
||||
_draw_markers(axes, markers, times)
|
||||
_set_limits(axes, times, plotted_values)
|
||||
if not plotted_values:
|
||||
axes.text(
|
||||
0.5,
|
||||
0.5,
|
||||
"No plottable TIC data",
|
||||
ha="center",
|
||||
va="center",
|
||||
transform=axes.transAxes,
|
||||
color="#555555",
|
||||
fontsize=14,
|
||||
)
|
||||
if axes.lines:
|
||||
axes.legend(loc="lower center", ncol=3, framealpha=0.92, fontsize=11)
|
||||
|
||||
with BytesIO() as output:
|
||||
canvas.print_png(output)
|
||||
return output.getvalue()
|
||||
|
||||
|
||||
def _curve_arrays(
|
||||
tic: TicResult,
|
||||
) -> tuple[AnyArray, dict[str, AnyArray], dict[str, AnyArray]]:
|
||||
times = np.asarray([row["time_sec"] for row in tic.rows], dtype=np.float64)
|
||||
raw_series = {
|
||||
name: np.asarray([row[f"{name}_raw"] for row in tic.rows], dtype=np.float64)
|
||||
for name in _CURVE_STYLES
|
||||
}
|
||||
processed_series = {
|
||||
name: np.asarray([row[f"{name}_processed"] for row in tic.rows], dtype=np.float64)
|
||||
for name in _CURVE_STYLES
|
||||
}
|
||||
return times, raw_series, processed_series
|
||||
|
||||
|
||||
def _style_axes(axes: Axes) -> None:
|
||||
axes.set_title("TIC: Lognormal fit", fontsize=20, pad=12)
|
||||
axes.set_xlabel("Time (s)", fontsize=15)
|
||||
axes.set_ylabel("Grayscale intensity", fontsize=15)
|
||||
axes.tick_params(axis="both", labelsize=12)
|
||||
axes.grid(True, color="#B8B8B8", linewidth=0.8, alpha=0.32)
|
||||
axes.set_axisbelow(True)
|
||||
|
||||
|
||||
def _draw_markers(
|
||||
axes: Axes,
|
||||
markers: Sequence[TicPlotMarker],
|
||||
times: AnyArray,
|
||||
) -> None:
|
||||
finite_times = times[np.isfinite(times)]
|
||||
if not finite_times.size:
|
||||
return
|
||||
time_min = float(finite_times.min())
|
||||
time_max = float(finite_times.max())
|
||||
for marker in markers:
|
||||
if not np.isfinite(marker.time_s) or not time_min <= marker.time_s <= time_max:
|
||||
continue
|
||||
axes.axvline(
|
||||
marker.time_s,
|
||||
color=tuple(channel / 255.0 for channel in reversed(marker.color)),
|
||||
linestyle="--",
|
||||
linewidth=1.5,
|
||||
alpha=0.9,
|
||||
label=marker.label,
|
||||
)
|
||||
|
||||
|
||||
def _set_limits(
|
||||
axes: Axes,
|
||||
times: AnyArray,
|
||||
series: Sequence[AnyArray],
|
||||
) -> None:
|
||||
finite_times = times[np.isfinite(times)]
|
||||
if finite_times.size:
|
||||
time_min = float(finite_times.min())
|
||||
time_max = float(finite_times.max())
|
||||
time_padding = max((time_max - time_min) * 0.025, 0.1)
|
||||
axes.set_xlim(time_min - time_padding, time_max + time_padding)
|
||||
|
||||
finite_values = [values[np.isfinite(values)] for values in series]
|
||||
finite_values = [values for values in finite_values if values.size]
|
||||
if not finite_values:
|
||||
return
|
||||
value_min = min(float(values.min()) for values in finite_values)
|
||||
value_max = max(float(values.max()) for values in finite_values)
|
||||
value_span = max(value_max - min(0.0, value_min), 1.0)
|
||||
axes.set_ylim(
|
||||
min(0.0, value_min) - 0.05 * value_span,
|
||||
value_max + 0.08 * value_span,
|
||||
)
|
||||
@@ -0,0 +1,245 @@
|
||||
"""trajectory.npz 的打包和解包：逐帧 ROI 多边形、暗区、背景圆和几何元数据（格式同旧项目）。
|
||||
|
||||
键集对齐参考组件 `write_roi_trajectory_npz`，另有 `keyframe_indices`（提示帧帧号）、
|
||||
`dark_polygons`、`calibration_json` 三键。
|
||||
|
||||
- 解包用 `np.load(..., allow_pickle=False)`：变长和字符串数据一律存 JSON，不存 object 数组，
|
||||
读取时不会执行任意代码；
|
||||
- 坏字节（不是合法 npz）解包直接报错；
|
||||
- 打包前校验形状（多边形 `(N, 64, 2)`、`frame_indices` 长 N、`reference_mask` 二维等），不一致报错。
|
||||
|
||||
`source_effective_rect` / `target_effective_rect` 两键保留是为了格式同旧项目：旧 debug 路径不裁
|
||||
黑边，它们恒等于对应的视图矩形。照搬旧项目 `trajectory.py`，去掉产品层。
|
||||
"""
|
||||
|
||||
import io
|
||||
import json
|
||||
from dataclasses import dataclass
|
||||
from typing import Any, cast
|
||||
|
||||
import numpy as np
|
||||
|
||||
from cadx.ceus._geometry import AnyArray, ViewRect
|
||||
|
||||
# 重采样多边形的固定顶点数（与 geometry.resample_closed_contour 的下游约定一致）。
|
||||
POLYGON_NUM_POINTS = 64
|
||||
|
||||
|
||||
@dataclass
|
||||
class Trajectory:
|
||||
"""逐帧 ROI 轨迹工件的内存表示（npz pack/unpack 的双向载体）。
|
||||
|
||||
属性:
|
||||
frame_indices: ``(N,)`` int32，逐帧索引（与多边形帧数一致）。
|
||||
fps: 视频帧率（浮点标量）。npz 存档内 fps 存为 float32（schema 保证，见
|
||||
``pack_trajectory``）；in-memory ``Trajectory.fps`` 有意规范化为 Python
|
||||
float（与领域实体 ``CeusRoiTracking.fps`` 的 float 一致，便于消费）。
|
||||
target_polygons: ``(N, 64, 2)`` float32，目标视图逐帧 ROI 多边形。
|
||||
source_polygons: ``(N, 64, 2)`` float32，源视图逐帧 ROI 多边形。
|
||||
source_view_rect / target_view_rect: 源 / 目标视图整屏裁剪矩形。
|
||||
source_effective_rect / target_effective_rect: 源 / 目标有效成像区矩形。
|
||||
source_polygon_coordinate_space / target_polygon_coordinate_space: 多边形所处
|
||||
坐标系名（如 ``source_view_local`` / ``target_effective_view_local``）。
|
||||
reference_mask: ``(H, W)`` uint8，参考区二值掩膜。
|
||||
reference_region: 参考区元信息字典（pack 时 JSON 序列化）。
|
||||
keyframe_indices: ``(K,)`` int32，被选作提示的关键帧索引（变长一维）。
|
||||
dark_polygons: ``list[list[AnyArray]]``，逐帧 0~N 个暗区多边形（变长嵌套，
|
||||
pack 时整体 JSON 序列化为单个字符串，规避 npz 变长 / pickle）。
|
||||
calibration: 标定溯源字典（标定来源 / hash / 是否 CLI 临时覆盖 / 检测 split_x）。
|
||||
"""
|
||||
|
||||
frame_indices: AnyArray
|
||||
fps: float
|
||||
target_polygons: AnyArray
|
||||
source_polygons: AnyArray
|
||||
source_view_rect: ViewRect
|
||||
target_view_rect: ViewRect
|
||||
source_effective_rect: ViewRect
|
||||
target_effective_rect: ViewRect
|
||||
source_polygon_coordinate_space: str
|
||||
target_polygon_coordinate_space: str
|
||||
reference_mask: AnyArray
|
||||
reference_region: dict[str, Any]
|
||||
keyframe_indices: AnyArray
|
||||
dark_polygons: list[list[AnyArray]]
|
||||
calibration: dict[str, Any]
|
||||
|
||||
|
||||
def _rect_to_array(rect: ViewRect) -> AnyArray:
|
||||
"""把 ``ViewRect`` 编码为 ``(4,)`` int32 数组 ``[x1, y1, x2, y2]``（对齐参考 P:1676）。"""
|
||||
return np.asarray([rect.x1, rect.y1, rect.x2, rect.y2], dtype=np.int32)
|
||||
|
||||
|
||||
def _rect_from_array(values: AnyArray) -> ViewRect:
|
||||
"""把 ``(4,)`` 整型数组解码回 ``ViewRect``；不变量校验委托给 ViewRect 构造。"""
|
||||
if values.shape != (4,):
|
||||
raise ValueError(f"矩形数组必须形如 (4,)，得到 {values.shape}。")
|
||||
x1, y1, x2, y2 = (int(v) for v in values.tolist())
|
||||
return ViewRect(x1, y1, x2, y2)
|
||||
|
||||
|
||||
def _serialize_dark_polygons(dark_polygons: list[list[AnyArray]]) -> str:
|
||||
"""把逐帧暗区多边形（变长嵌套 ndarray）序列化为 JSON 字符串。
|
||||
|
||||
每个多边形 ``(M, 2)`` 转为嵌套 list；整体为 ``list[list[list[list[float]]]]``
|
||||
（帧 → 多边形 → 顶点 → 坐标）。走 JSON 而非 npz object 数组，从源头禁 pickle。
|
||||
"""
|
||||
serializable: list[list[list[list[float]]]] = []
|
||||
for frame_polygons in dark_polygons:
|
||||
frame_repr: list[list[list[float]]] = []
|
||||
for polygon in frame_polygons:
|
||||
array = np.asarray(polygon, dtype=np.float32)
|
||||
if array.ndim != 2 or array.shape[1] != 2:
|
||||
raise ValueError(f"暗区多边形必须形如 (M, 2)，得到 {array.shape}。")
|
||||
frame_repr.append(array.tolist())
|
||||
serializable.append(frame_repr)
|
||||
return json.dumps(serializable, ensure_ascii=False)
|
||||
|
||||
|
||||
def _deserialize_dark_polygons(payload: str) -> list[list[AnyArray]]:
|
||||
"""把暗区多边形 JSON 字符串还原为 ``list[list[AnyArray]]``（float32）。
|
||||
|
||||
严格按 JSON 结构还原，不接受任何 pickle / object 输入；结构非法直接 raise。
|
||||
"""
|
||||
parsed: object = json.loads(payload)
|
||||
if not isinstance(parsed, list):
|
||||
raise ValueError("暗区多边形 JSON 顶层必须是逐帧列表。")
|
||||
dark_polygons: list[list[AnyArray]] = []
|
||||
for frame_repr in cast(list[object], parsed):
|
||||
if not isinstance(frame_repr, list):
|
||||
raise ValueError("暗区多边形每帧必须是多边形列表。")
|
||||
frame_polygons: list[AnyArray] = []
|
||||
for polygon_repr in cast(list[object], frame_repr):
|
||||
array = np.asarray(polygon_repr, dtype=np.float32)
|
||||
if array.ndim != 2 or array.shape[1] != 2:
|
||||
raise ValueError(f"暗区多边形必须形如 (M, 2)，得到 {array.shape}。")
|
||||
frame_polygons.append(array)
|
||||
dark_polygons.append(frame_polygons)
|
||||
return dark_polygons
|
||||
|
||||
|
||||
def _validate_geometry(trajectory: Trajectory) -> int:
|
||||
"""校验轨迹核心几何形状的一致性，返回帧数 N；任何不一致直接 raise。
|
||||
|
||||
校验项：多边形 ``(N, 64, 2)`` 且源 / 目标帧数一致；``frame_indices`` 一维且长 N；
|
||||
``reference_mask`` 二维；``keyframe_indices`` 一维；``dark_polygons`` 帧数为 N。
|
||||
"""
|
||||
target = trajectory.target_polygons
|
||||
source = trajectory.source_polygons
|
||||
if target.ndim != 3 or target.shape[1:] != (POLYGON_NUM_POINTS, 2):
|
||||
raise ValueError(
|
||||
f"target_polygons 必须形如 (N, {POLYGON_NUM_POINTS}, 2)，得到 {target.shape}。"
|
||||
)
|
||||
if source.shape != target.shape:
|
||||
raise ValueError(
|
||||
f"source_polygons {source.shape} 与 target_polygons {target.shape} 形状不一致。"
|
||||
)
|
||||
num_frames = int(target.shape[0])
|
||||
if trajectory.frame_indices.ndim != 1 or trajectory.frame_indices.shape[0] != num_frames:
|
||||
raise ValueError(
|
||||
f"frame_indices 必须一维且长度 {num_frames}，得到 {trajectory.frame_indices.shape}。"
|
||||
)
|
||||
if trajectory.reference_mask.ndim != 2:
|
||||
raise ValueError(
|
||||
f"reference_mask 必须二维 (H, W)，得到 {trajectory.reference_mask.shape}。"
|
||||
)
|
||||
if trajectory.keyframe_indices.ndim != 1:
|
||||
raise ValueError(f"keyframe_indices 必须一维，得到 {trajectory.keyframe_indices.shape}。")
|
||||
if len(trajectory.dark_polygons) != num_frames:
|
||||
raise ValueError(
|
||||
f"dark_polygons 帧数 {len(trajectory.dark_polygons)} 与多边形帧数 {num_frames} 不一致。"
|
||||
)
|
||||
return num_frames
|
||||
|
||||
|
||||
def pack_trajectory(trajectory: Trajectory) -> bytes:
|
||||
"""把 ``Trajectory`` 校验并打包为压缩 npz 字节（产物工件）。
|
||||
|
||||
参数:
|
||||
trajectory: 待打包的逐帧轨迹载体。
|
||||
|
||||
返回:
|
||||
压缩 npz 的原始字节（``np.savez_compressed`` 写入内存缓冲）。
|
||||
|
||||
关键实现:
|
||||
先 ``_validate_geometry`` 钉死核心形状（不一致即 raise，不静默对齐）；
|
||||
字符串键存为 0 维 str 数组，``fps`` 存 ``(1,)`` float32，变长 dark / 字典键
|
||||
统一 JSON 序列化为字符串，全程不写入任何 object 数组（禁 pickle）。
|
||||
键名对齐参考 ``write_roi_trajectory_npz`` 并新增 ``keyframe_indices`` /
|
||||
``dark_polygons`` / ``calibration_json``。
|
||||
"""
|
||||
_validate_geometry(trajectory)
|
||||
|
||||
buffer = io.BytesIO()
|
||||
np.savez_compressed(
|
||||
buffer,
|
||||
frame_indices=trajectory.frame_indices.astype(np.int32),
|
||||
fps=np.asarray([trajectory.fps], dtype=np.float32),
|
||||
target_polygons=trajectory.target_polygons.astype(np.float32),
|
||||
source_polygons=trajectory.source_polygons.astype(np.float32),
|
||||
source_view_rect=_rect_to_array(trajectory.source_view_rect),
|
||||
target_view_rect=_rect_to_array(trajectory.target_view_rect),
|
||||
source_effective_rect=_rect_to_array(trajectory.source_effective_rect),
|
||||
target_effective_rect=_rect_to_array(trajectory.target_effective_rect),
|
||||
source_polygon_coordinate_space=np.array(trajectory.source_polygon_coordinate_space),
|
||||
target_polygon_coordinate_space=np.array(trajectory.target_polygon_coordinate_space),
|
||||
reference_mask=trajectory.reference_mask.astype(np.uint8),
|
||||
reference_region_json=np.array(json.dumps(trajectory.reference_region, ensure_ascii=False)),
|
||||
keyframe_indices=trajectory.keyframe_indices.astype(np.int32),
|
||||
dark_polygons=np.array(_serialize_dark_polygons(trajectory.dark_polygons)),
|
||||
calibration_json=np.array(json.dumps(trajectory.calibration, ensure_ascii=False)),
|
||||
)
|
||||
return buffer.getvalue()
|
||||
|
||||
|
||||
def unpack_trajectory(raw: bytes) -> Trajectory:
|
||||
"""把 npz 字节解包还原为 ``Trajectory``（禁 pickle，坏字节 raise）。
|
||||
|
||||
参数:
|
||||
raw: ``pack_trajectory`` 产出的 npz 原始字节。
|
||||
|
||||
返回:
|
||||
还原后的 ``Trajectory``，dtype 与 pack 时一致。
|
||||
|
||||
关键实现:
|
||||
``np.load(..., allow_pickle=False)`` 硬禁 pickle —— 任何含 object pickle 的
|
||||
恶意字节都会在 load 阶段抛错而非被反序列化执行；非法 npz / 缺键同样直接 raise
|
||||
（``ValueError`` / ``OSError`` / ``KeyError``），不静默返回空。字符串键经 0 维
|
||||
数组 ``.item()`` 取回，变长 / 字典键经 JSON 反序列化还原。
|
||||
``fps`` 存档内为 float32（schema 保证），此处有意还原为 Python float
|
||||
（领域消费侧统一用 float，不改此行为）。
|
||||
"""
|
||||
with np.load(io.BytesIO(raw), allow_pickle=False) as archive:
|
||||
frame_indices = archive["frame_indices"].astype(np.int32)
|
||||
fps = float(archive["fps"][0])
|
||||
target_polygons = archive["target_polygons"].astype(np.float32)
|
||||
source_polygons = archive["source_polygons"].astype(np.float32)
|
||||
source_view_rect = _rect_from_array(archive["source_view_rect"])
|
||||
target_view_rect = _rect_from_array(archive["target_view_rect"])
|
||||
source_effective_rect = _rect_from_array(archive["source_effective_rect"])
|
||||
target_effective_rect = _rect_from_array(archive["target_effective_rect"])
|
||||
source_polygon_coordinate_space = str(archive["source_polygon_coordinate_space"].item())
|
||||
target_polygon_coordinate_space = str(archive["target_polygon_coordinate_space"].item())
|
||||
reference_mask = archive["reference_mask"].astype(np.uint8)
|
||||
reference_region = json.loads(str(archive["reference_region_json"].item()))
|
||||
keyframe_indices = archive["keyframe_indices"].astype(np.int32)
|
||||
dark_polygons = _deserialize_dark_polygons(str(archive["dark_polygons"].item()))
|
||||
calibration = json.loads(str(archive["calibration_json"].item()))
|
||||
|
||||
return Trajectory(
|
||||
frame_indices=frame_indices,
|
||||
fps=fps,
|
||||
target_polygons=target_polygons,
|
||||
source_polygons=source_polygons,
|
||||
source_view_rect=source_view_rect,
|
||||
target_view_rect=target_view_rect,
|
||||
source_effective_rect=source_effective_rect,
|
||||
target_effective_rect=target_effective_rect,
|
||||
source_polygon_coordinate_space=source_polygon_coordinate_space,
|
||||
target_polygon_coordinate_space=target_polygon_coordinate_space,
|
||||
reference_mask=reference_mask,
|
||||
reference_region=reference_region,
|
||||
keyframe_indices=keyframe_indices,
|
||||
dark_polygons=dark_polygons,
|
||||
calibration=calibration,
|
||||
)
|
||||
@@ -0,0 +1,10 @@
|
||||
"""segment：用 MedSAM2 把提示帧的病灶掩码传播成逐帧掩码。
|
||||
|
||||
不认识分屏和绿框：输入只是要分割的那一半画面和那一半局部坐标上的提示掩码。`cadx` 里只有这里
|
||||
引用 `torch` 和 `sam2`。
|
||||
"""
|
||||
|
||||
from cadx.segment._config import SegmentConfig
|
||||
from cadx.segment._medsam2 import medsam2_commit, segment_video
|
||||
|
||||
__all__ = ["SegmentConfig", "medsam2_commit", "segment_video"]
|
||||
@@ -0,0 +1,44 @@
|
||||
"""分割设置：`config/ceus.yaml` 里 `segment` 段的字段定义和校验。
|
||||
|
||||
每个字段必填；多写或拼错的键报错；读入后不能修改（`AGENTS.md`「不设默认值」）。取值与出处写在
|
||||
`config/ceus.yaml` 的注释里。路径字段原样使用，相对路径由调用者（`scripts/ceus_quant.py`）
|
||||
先换成仓库根目录下的绝对路径。
|
||||
"""
|
||||
|
||||
from pathlib import Path
|
||||
from typing import Annotated
|
||||
|
||||
from pydantic import BaseModel, BeforeValidator, ConfigDict, Field
|
||||
|
||||
|
||||
def _reject_empty_path(value: object) -> object:
|
||||
# Path("") 会变成 "."（当前目录），空串在转换前就拦下。
|
||||
if value == "":
|
||||
raise ValueError("路径不能是空串")
|
||||
return value
|
||||
|
||||
|
||||
# YAML 里的路径是字符串；路径字段是唯一允许从字符串转换类型的地方。
|
||||
type _PathField = Annotated[Path, Field(strict=False), BeforeValidator(_reject_empty_path)]
|
||||
|
||||
|
||||
class SegmentConfig(BaseModel):
|
||||
"""MedSAM2 分割设置。
|
||||
|
||||
- `medsam2_dir`：MedSAM2 代码目录（含 `.git`，提交号从这里读；装进环境的 `sam2` 包必须来自
|
||||
这个目录下的 `sam2/`）；
|
||||
- `checkpoint`：权重文件；
|
||||
- `sam2_config`：模型配置的 Hydra 名，相对 `sam2` 包，如 `configs/sam2.1_hiera_t512.yaml`；
|
||||
- `device`：推理设备，如 `cuda`、`cuda:0`、`cpu`；
|
||||
- `offload_video_to_cpu`、`offload_state_to_cpu`：MedSAM2 `init_state` 的两个开关，把帧或推理
|
||||
状态放在内存而不是显存，省显存、略慢，不改变结果。
|
||||
"""
|
||||
|
||||
model_config = ConfigDict(strict=True, extra="forbid", frozen=True)
|
||||
|
||||
medsam2_dir: _PathField
|
||||
checkpoint: _PathField
|
||||
sam2_config: str
|
||||
device: str
|
||||
offload_video_to_cpu: bool
|
||||
offload_state_to_cpu: bool
|
||||
@@ -0,0 +1,166 @@
|
||||
"""MedSAM2 逐帧分割：把若干提示帧的病灶掩码传播到整段视频的每一帧。
|
||||
|
||||
改写自旧项目 `debugs/medsam2_local_track.py`（当时在 conda 子进程里跑），行为照旧：帧以
|
||||
`000000.jpg` 起的 jpg 序列（质量 95）写进一个目录交给 `init_state`；每个提示帧用 `add_new_mask`
|
||||
加进同一个对象；从最早的提示帧向后传播到末帧，再从它向前传播到第 0 帧；logits > 0 为病灶。
|
||||
现在在同一进程里调用，jpg 临时目录由这里建、用完即删。
|
||||
"""
|
||||
|
||||
import os
|
||||
import shutil
|
||||
import subprocess
|
||||
import tempfile
|
||||
from collections.abc import Mapping, Sequence
|
||||
from pathlib import Path
|
||||
from typing import Any
|
||||
|
||||
import cv2
|
||||
import numpy as np
|
||||
import numpy.typing as npt
|
||||
import torch
|
||||
|
||||
from cadx.segment._config import SegmentConfig
|
||||
|
||||
type Frame = npt.NDArray[np.uint8]
|
||||
type Mask = npt.NDArray[np.bool_]
|
||||
|
||||
# 只跟踪一个病灶，对象编号照旧脚本取 1。
|
||||
_OBJECT_ID = 1
|
||||
_JPEG_QUALITY = 95
|
||||
|
||||
|
||||
def medsam2_commit(medsam2_dir: Path) -> str:
|
||||
"""读 `medsam2_dir/.git` 的 HEAD 提交号（完整 40 位）。
|
||||
|
||||
用 `--git-dir` 指定仓库：只用 `-C` 时，目录里没有 `.git` 的话 git 会往上找到 CADx-4 自己的仓库。
|
||||
"""
|
||||
git = shutil.which("git")
|
||||
if git is None:
|
||||
raise RuntimeError("找不到 git，无法读 MedSAM2 的提交号")
|
||||
completed = subprocess.run( # noqa: S603 参数是写死的 git 子命令和配置里的目录
|
||||
[git, "--git-dir", str(medsam2_dir / ".git"), "rev-parse", "HEAD"],
|
||||
check=True,
|
||||
capture_output=True,
|
||||
)
|
||||
return completed.stdout.decode().strip()
|
||||
|
||||
|
||||
def _check_inputs(
|
||||
frames: Sequence[Frame], prompt_masks: Mapping[int, Mask], cfg: SegmentConfig
|
||||
) -> tuple[int, int]:
|
||||
"""校验输入并返回帧的 (高, 宽)。"""
|
||||
for path in (cfg.medsam2_dir / "sam2" / cfg.sam2_config, cfg.checkpoint):
|
||||
if not path.is_file():
|
||||
raise FileNotFoundError(f"MedSAM2 缺少文件：{path}")
|
||||
if len(frames) == 0:
|
||||
raise ValueError("没有要分割的帧")
|
||||
height, width = frames[0].shape[:2]
|
||||
for index, frame in enumerate(frames):
|
||||
if frame.dtype != np.uint8 or frame.shape != (height, width, 3):
|
||||
raise ValueError(
|
||||
f"第 {index} 帧应为 ({height}, {width}, 3) 的 uint8 BGR 图像，"
|
||||
f"收到 {frame.shape} {frame.dtype}"
|
||||
)
|
||||
if len(prompt_masks) == 0:
|
||||
raise ValueError("至少需要一个提示帧")
|
||||
for index, mask in prompt_masks.items():
|
||||
if not 0 <= index < len(frames):
|
||||
raise ValueError(f"提示帧帧号 {index} 不在 [0, {len(frames)}) 内")
|
||||
if mask.dtype != np.bool_ or mask.shape != (height, width):
|
||||
raise ValueError(
|
||||
f"提示帧 {index} 的掩码应为 ({height}, {width}) 的布尔数组，"
|
||||
f"收到 {mask.shape} {mask.dtype}"
|
||||
)
|
||||
return height, width
|
||||
|
||||
|
||||
def _check_device(device: str) -> None:
|
||||
"""加载权重前确认请求的设备可见，提前暴露显卡可见性问题。"""
|
||||
if not device.startswith("cuda"):
|
||||
return
|
||||
if not torch.cuda.is_available():
|
||||
visible = os.environ.get("CUDA_VISIBLE_DEVICES", "<unset>")
|
||||
raise RuntimeError(
|
||||
f"请求用 CUDA 运行 MedSAM2，但当前进程看不到显卡。CUDA_VISIBLE_DEVICES={visible}"
|
||||
)
|
||||
if ":" in device and int(device.split(":", 1)[1]) >= torch.cuda.device_count():
|
||||
raise RuntimeError(f"请求设备 {device}，但只看得到 {torch.cuda.device_count()} 张显卡")
|
||||
|
||||
|
||||
def _check_sam2_source(medsam2_dir: Path) -> None:
|
||||
"""装进环境的 sam2 必须就是 `medsam2_dir` 里那份代码，否则记下的提交号对不上实际代码。"""
|
||||
import sam2
|
||||
|
||||
loaded = Path(sam2.__file__).resolve().parent
|
||||
expected = (medsam2_dir / "sam2").resolve()
|
||||
if loaded != expected:
|
||||
raise RuntimeError(f"环境里的 sam2 来自 {loaded}，与配置的 {expected} 不一致")
|
||||
|
||||
|
||||
def _write_jpgs(frames: Sequence[Frame], frames_dir: Path) -> None:
|
||||
for index, frame in enumerate(frames):
|
||||
path = frames_dir / f"{index:06d}.jpg"
|
||||
if not cv2.imwrite(str(path), frame, [cv2.IMWRITE_JPEG_QUALITY, _JPEG_QUALITY]):
|
||||
raise RuntimeError(f"写分割用的 jpg 失败：{path}")
|
||||
|
||||
|
||||
def _to_mask(mask_logits: torch.Tensor) -> Mask:
|
||||
return np.squeeze((mask_logits > 0.0).detach().cpu().numpy()).astype(np.bool_)
|
||||
|
||||
|
||||
def _propagate(predictor: Any, state: Any, first_prompt: int) -> dict[int, Mask]:
|
||||
"""从最早的提示帧向后传播到末帧，再向前传播到第 0 帧。"""
|
||||
propagated: dict[int, Mask] = {}
|
||||
with torch.inference_mode():
|
||||
for frame_idx, obj_ids, mask_logits in predictor.propagate_in_video(
|
||||
state, start_frame_idx=first_prompt
|
||||
):
|
||||
propagated[int(frame_idx)] = _to_mask(mask_logits[list(obj_ids).index(_OBJECT_ID)])
|
||||
for frame_idx, obj_ids, mask_logits in predictor.propagate_in_video(
|
||||
state,
|
||||
start_frame_idx=first_prompt,
|
||||
max_frame_num_to_track=first_prompt + 1,
|
||||
reverse=True,
|
||||
):
|
||||
propagated[int(frame_idx)] = _to_mask(mask_logits[list(obj_ids).index(_OBJECT_ID)])
|
||||
return propagated
|
||||
|
||||
|
||||
def segment_video(
|
||||
frames: Sequence[Frame], prompt_masks: Mapping[int, Mask], cfg: SegmentConfig
|
||||
) -> Mask:
|
||||
"""把提示帧的病灶掩码传播到每一帧。
|
||||
|
||||
`frames` 是逐帧 BGR 图像（同一尺寸），`prompt_masks` 是提示帧帧号 → 与帧同尺寸的布尔掩码。
|
||||
返回 `(帧数, 高, 宽)` 的布尔数组，与输入帧一一对应。
|
||||
"""
|
||||
height, width = _check_inputs(frames, prompt_masks, cfg)
|
||||
_check_device(cfg.device)
|
||||
_check_sam2_source(cfg.medsam2_dir)
|
||||
# sam2 没有类型标注，参数类型是 Unknown；返回的 predictor 按 Any 用。
|
||||
from sam2.build_sam import build_sam2_video_predictor # pyright: ignore[reportUnknownVariableType]
|
||||
|
||||
predictor: Any = build_sam2_video_predictor(
|
||||
config_file=cfg.sam2_config, ckpt_path=str(cfg.checkpoint), device=cfg.device
|
||||
)
|
||||
with tempfile.TemporaryDirectory(prefix="cadx-medsam2-") as tmp:
|
||||
frames_dir = Path(tmp)
|
||||
_write_jpgs(frames, frames_dir)
|
||||
state: Any = predictor.init_state(
|
||||
str(frames_dir),
|
||||
offload_video_to_cpu=cfg.offload_video_to_cpu,
|
||||
offload_state_to_cpu=cfg.offload_state_to_cpu,
|
||||
)
|
||||
for index in sorted(prompt_masks):
|
||||
predictor.add_new_mask(state, frame_idx=index, obj_id=_OBJECT_ID, mask=prompt_masks[index])
|
||||
propagated = _propagate(predictor, state, min(prompt_masks))
|
||||
|
||||
missing = [index for index in range(len(frames)) if index not in propagated]
|
||||
if missing:
|
||||
raise RuntimeError(f"MedSAM2 传播后缺 {len(missing)} 帧的掩码，前几帧：{missing[:5]}")
|
||||
stack = np.stack([propagated[index] for index in range(len(frames))], axis=0)
|
||||
if stack.shape != (len(frames), height, width):
|
||||
raise RuntimeError(
|
||||
f"MedSAM2 输出掩码形状 {stack.shape}，应为 {(len(frames), height, width)}"
|
||||
)
|
||||
return stack
|
||||
Reference in New Issue
Block a user