Files
2026-10-08 16:43:48 +09:00

186 KiB

ISR01 정답 코드 - 공개 위성영상 장비·차량 객체분류 정확도 평가

[교육용가상자료] 모든 데이터는 코드 안에서 난수로 만든 가상 자료이며 실제 부대·장비·인원·작전 정보와 무관합니다.

실행 환경

  • Python 3.10 이상, numpy, pandas, scipy, scikit-learn, matplotlib
  • 미션별 코드는 각각 한 파일로 독립 실행됩니다. 외부 파일 없이 코드 안에서 가상 데이터를 생성합니다.
  • 그래프의 글자는 한글 글꼴이 없는 실행 환경에서도 깨지지 않도록 영문으로 표기했습니다. 표·출력 결과는 한글입니다.
  • 난수 시드를 고정해 두어 같은 코드를 실행하면 같은 결과가 나옵니다.

미션 1. 상용 위성(광학/SAR) 영상 특성 분석 및 군용 장비/차량 메타데이터 규격 정의

  • 난이도: 초급 / 예상 소요시간: 1.5시간
  • 코드 설명: 공개·상용 위성영상(광학/SAR)의 해상도·촬영 조건 메타데이터와 군용 장비·차량(전차, 자주포, 트럭) 라벨의 데이터셋 스키마를 정의하고, 가상 샘플 메타데이터를 검증한다.
"""공개·상용 위성영상(광학/SAR)의 해상도·촬영 조건 메타데이터와 군용 장비·차량(전차, 자주포, 트럭) 라벨의 데이터셋 스키마를 정의하고, 가상 샘플 메타데이터를 검증한다."""
import json
import re
import numpy as np
import pandas as pd

# 1) 스키마 정의: 영상 출처·라이선스, 센서·촬영 조건, 전술 파라미터, 객체 라벨(클래스·위치·방위·은폐 상태)
SCHEMA = {
    "required": ["image_id", "source", "sensor", "tactical", "labels"],
    "fields": {
        "image_id": {"type": str, "pattern": r"^IMG\d{4}$"},
    },
    "source": {
        "provider": {"type": str, "enum": ["open_data", "commercial_sample", "partner_share"]},
        "license": {"type": str, "pattern": r"^\S.{2,}$"},
        "acquired_at": {"type": str, "pattern": r"^\d{4}-\d{2}-\d{2}T\d{2}:\d{2}Z$"},
    },
    "sensor": {
        "type": {"type": str, "enum": ["optical", "SAR"]},
        "gsd_m": {"type": (int, float), "min": 0.1, "max": 30},
        "off_nadir_deg": {"type": (int, float), "min": 0, "max": 60},
        "sun_elevation_deg": {"type": (int, float), "min": -10, "max": 90},
        "cloud_cover_pct": {"type": (int, float), "min": 0, "max": 100},
    },
    "tactical": {
        "weather": {"type": str, "enum": ["clear", "haze", "cloudy", "snow"]},
        "geo_error_m": {"type": (int, float), "min": 0, "max": 1000},
        "delivery_latency_h": {"type": (int, float), "min": 0, "max": 720},
    },
    "label": {
        "class": {"type": str, "enum": ["tank", "spg", "truck", "civil_vehicle", "decoy", "unknown"]},
        "bbox": {"type": list, "length": 4},
        "orientation_deg": {"type": (int, float), "min": 0, "max": 360},
        "concealment": {"type": str, "enum": ["none", "camo_net", "tree", "cloud", "partial"]},
        "status": {"type": str, "enum": ["auto_draft", "analyst_verified"]},
    },
}


def check(name, value, rule):
    """값 하나를 규칙(type, enum, min, max, length, pattern)으로 검사해 오류 목록을 돌려준다."""
    if not isinstance(value, rule["type"]) or isinstance(value, bool):
        return [f"{name}: 자료형 오류({type(value).__name__})"]
    errs = []
    if "enum" in rule and value not in rule["enum"]:
        errs.append(f"{name}: 허용되지 않은 값 '{value}'")
    if "min" in rule and value < rule["min"]:
        errs.append(f"{name}: 최소값 {rule['min']} 미만")
    if "max" in rule and value > rule["max"]:
        errs.append(f"{name}: 최대값 {rule['max']} 초과")
    if "length" in rule and len(value) != rule["length"]:
        errs.append(f"{name}: 길이 {rule['length']} 필요")
    if "pattern" in rule and not re.match(rule["pattern"], value):
        errs.append(f"{name}: 형식 불일치 '{value}'")
    return errs


def validate(rec):
    errs = [f"{k}: 필수 필드 누락" for k in SCHEMA["required"] if k not in rec]
    for k, rule in SCHEMA["fields"].items():
        if k in rec:
            errs += check(k, rec[k], rule)
    for group in ("source", "sensor", "tactical"):
        for k, rule in SCHEMA[group].items():
            if k not in rec.get(group, {}):
                errs.append(f"{group}.{k}: 필수 필드 누락")
            else:
                errs += check(f"{group}.{k}", rec[group][k], rule)
    for i, lab in enumerate(rec.get("labels", [])):
        for k, rule in SCHEMA["label"].items():
            if k not in lab:
                errs.append(f"labels[{i}].{k}: 필수 필드 누락")
            else:
                errs += check(f"labels[{i}].{k}", lab[k], rule)
    return errs


# 2) [교육용가상자료] 영상 메타데이터 20건 생성 (가상 날짜·수치, 실제 촬영 정보 아님)
rng = np.random.default_rng(1)
classes = SCHEMA["label"]["class"]["enum"]
records = []
for i in range(20):
    sensor = str(rng.choice(["optical", "SAR"]))
    labels = [{"class": str(rng.choice(classes)), "bbox": [int(rng.integers(0, 480)), int(rng.integers(0, 480)), 32, 32],
               "orientation_deg": float(round(rng.uniform(0, 360), 1)),
               "concealment": str(rng.choice(SCHEMA["label"]["concealment"]["enum"])), "status": "auto_draft"} for _ in range(3)]
    records.append({"image_id": f"IMG{i:04d}",
                    "source": {"provider": str(rng.choice(SCHEMA["source"]["provider"]["enum"])), "license": "CC-BY-4.0(가상)",
                               "acquired_at": f"2026-0{int(rng.integers(1, 9))}-1{int(rng.integers(0, 9))}T0{int(rng.integers(0, 9))}:30Z"},
                    "sensor": {"type": sensor, "gsd_m": float(round(rng.uniform(0.3, 3.0), 2)), "off_nadir_deg": float(round(rng.uniform(0, 40), 1)),
                               "sun_elevation_deg": float(round(rng.uniform(15, 70), 1)), "cloud_cover_pct": float(round(rng.uniform(0, 60), 1))},
                    "tactical": {"weather": str(rng.choice(["clear", "haze", "cloudy", "snow"])),
                                 "geo_error_m": float(round(rng.uniform(1, 30), 1)), "delivery_latency_h": float(round(rng.uniform(1, 72), 1))},
                    "labels": labels})

# 3) 정상 데이터 검증
ok = sum(1 for r in records if not validate(r))
print(f"정상 샘플 {len(records)}건 중 검증 통과 {ok}건, 오류 {len(records) - ok}건 (목표: 오류 0건 {'달성' if ok == len(records) else '미달'})")

# 4) 오류 주입 시험: 검증기가 잘못된 데이터를 잡아내는지 확인
bad = json.loads(json.dumps(records[0]))
bad["sensor"]["gsd_m"] = 0.0
bad["sensor"]["type"] = "IR"
del bad["tactical"]["geo_error_m"]
bad["labels"][0]["class"] = "bmp"
bad["labels"][1]["bbox"] = [1, 2, 3]
print("\n[오류 주입 결과] 주입 5종")
for e in validate(bad):
    print(" -", e)

# 5) 스키마 필드 정의서(명세서 초안)
rows = [("기본", k, r["type"].__name__, str({x: y for x, y in r.items() if x != "type"})) for k, r in SCHEMA["fields"].items()]
for g in ("source", "sensor", "tactical", "label"):
    for k, r in SCHEMA[g].items():
        t = r["type"].__name__ if isinstance(r["type"], type) else "number"
        rows.append((g, k, t, str({x: y for x, y in r.items() if x != "type"})))
spec = pd.DataFrame(rows, columns=["그룹", "필드", "자료형", "규칙"])
print("\n[데이터 명세 요약]")
print(spec.to_string(index=False))
print(f"\n필드 정의 {len(spec)}개(완성도 100%). 라벨 status는 'auto_draft'로 생성되고 영상분석관 판독 후에만 'analyst_verified'로 바뀐다.")
print("분류 기준: tank(포탑+긴 포신), spg(후방 대형 포탑+매우 긴 포신), truck(운전석+적재함 분리), civil_vehicle(소형·밝은/어두운 단색), decoy(모의 위장 표적)")

실행 결과 예시

정상 샘플 20건 중 검증 통과 20건, 오류 0건 (목표: 오류 0건 달성)

[오류 주입 결과] 주입 5종
 - sensor.type: 허용되지 않은 값 'IR'
 - sensor.gsd_m: 최소값 0.1 미만
 - tactical.geo_error_m: 필수 필드 누락
 - labels[0].class: 허용되지 않은 값 'bmp'
 - labels[1].bbox: 길이 4 필요

[데이터 명세 요약]
      그룹                 필드    자료형                                                                      규칙
      기본           image_id    str                                              {'pattern': '^IMG\\d{4}$'}
  source           provider    str           {'enum': ['open_data', 'commercial_sample', 'partner_share']}
  source            license    str                                               {'pattern': '^\\S.{2,}$'}
  source        acquired_at    str                    {'pattern': '^\\d{4}-\\d{2}-\\d{2}T\\d{2}:\\d{2}Z$'}
  sensor               type    str                                            {'enum': ['optical', 'SAR']}
  sensor              gsd_m number                                                 {'min': 0.1, 'max': 30}
  sensor      off_nadir_deg number                                                   {'min': 0, 'max': 60}
  sensor  sun_elevation_deg number                                                 {'min': -10, 'max': 90}
  sensor    cloud_cover_pct number                                                  {'min': 0, 'max': 100}
tactical            weather    str                           {'enum': ['clear', 'haze', 'cloudy', 'snow']}
tactical        geo_error_m number                                                 {'min': 0, 'max': 1000}
tactical delivery_latency_h number                                                  {'min': 0, 'max': 720}
   label              class    str {'enum': ['tank', 'spg', 'truck', 'civil_vehicle', 'decoy', 'unknown']}
   label               bbox   list                                                           {'length': 4}
   label    orientation_deg number                                                  {'min': 0, 'max': 360}
   label        concealment    str              {'enum': ['none', 'camo_net', 'tree', 'cloud', 'partial']}
   label             status    str                            {'enum': ['auto_draft', 'analyst_verified']}

필드 정의 17개(완성도 100%). 라벨 status는 'auto_draft'로 생성되고 영상분석관 판독 후에만 'analyst_verified'로 바뀐다.
분류 기준: tank(포탑+긴 포신), spg(후방 대형 포탑+매우 긴 포신), truck(운전석+적재함 분리), civil_vehicle(소형·밝은/어두운 단색), decoy(모의 위장 표적)

미션 2. 위성영상 해상도 표준화, 구름 차폐 제거 및 방위/스케일 정규화 전처리 모듈 구현

  • 난이도: 초급~중급 / 예상 소요시간: 2.0시간
  • 코드 설명: 위성영상 조각의 해상도(GSD)·비스듬 촬영 왜곡을 표준화하고, 구름 차폐 영역 처리, 결측 픽셀 보정, 객체 주축 기준 방위 정규화, 0~1 정규화를 수행하는 전처리 모듈을 구현해 처리 속도와 보정 효과를 검증한다.
"""위성영상 조각의 해상도(GSD)·비스듬 촬영 왜곡을 표준화하고, 구름 차폐 영역 처리, 결측 픽셀 보정, 객체 주축 기준 방위 정규화, 0~1 정규화를 수행하는 전처리 모듈을 구현해 처리 속도와 보정 효과를 검증한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage

# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5                                   # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
    "tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
    "spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
    "truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
    "civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
    "decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
    "none": [],                                 # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12}   # 금속 반사 세기(가상)


def _mask(cls, angle, scale, squash=1.0):
    """회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
    yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
    yy = yy / squash
    a = np.radians(angle)
    u = (xx * np.cos(a) + yy * np.sin(a)) / scale
    v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
    return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]


def make_chip(cls, rng, cond=None):
    """chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
    c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
    gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
    scale = STD_GSD / gsd
    angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
    squash = np.cos(np.radians(off))
    sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
    sun_az = rng.uniform(0, 360)
    bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
    parts = _mask(cls, angle, scale, squash)
    body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
    opt = bg.copy()
    # 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
    L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
    sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
    opt[sh & ~body] -= 35 if cls != "decoy" else 25
    if cls == "civil":
        base = rng.choice([200.0, 60.0])            # 민간 차량 색상: 밝거나 어두움
        opt[body] = base + rng.normal(0, 8, body.sum())
        roof = parts[1][0]
        opt[roof] = base * 0.8
    else:
        ref = bg[body].mean() if body.any() else bg.mean()
        for m, d in parts:
            opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
    # 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
    if c["camo"] > 0 and cls != "civil":
        net = ndimage.binary_dilation(body, iterations=2)
        tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
        opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
    # 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
    if c["occl"] > 0:
        trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
        near = ndimage.binary_dilation(body, iterations=3)
        cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
        opt[cover] = 70 + rng.normal(0, 6, cover.sum())
    # 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
    cmask = np.zeros((N, N), bool)
    if c["cloud"] > 0:
        blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
        cmask = blob > np.quantile(blob, 1 - c["cloud"])
        opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
    # SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
    sar = np.full((N, N), 0.12)
    sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
    edge = body & ~ndimage.binary_erosion(body)
    sar[edge] *= 1.3 if cls in MILITARY else 1.0
    sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
    if c["low_res"]:                                # 저해상도: 정보 손실(블러)
        down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
        opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
    meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
            "sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
            "cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
    return np.clip(opt, 0, 255), sar, meta


SCEN = {  # 은폐·촬영 조건 시나리오 10종
    "clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
    "low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
    "speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}


def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
    rng = np.random.default_rng(seed)
    O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
    for i in range(n):
        cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
        O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
        m.update({"id": f"IMG{i:04d}", "scen": sc})
        M.append(m)
    return O, S, M
# ---------------------------------------------------------------------------

# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
    """중앙 기준으로 32x32에 맞게 자르거나 채운다."""
    out = np.full((N, N), fill, float)
    h, w = a.shape
    ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
    yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
    hh, ww = min(N, h), min(N, w)
    out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
    return out


def cloud_mask(opt, sm, min_area=15):
    """광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
    bright = opt > 185
    if bright.sum() < min_area:
        return np.zeros_like(bright)
    local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
    cand = bright & (local_sd < 12) & (sm < 0.3)
    lab, n = ndimage.label(cand)
    if n == 0:
        return cand
    sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
    return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)


def orientation(opt, sm):
    """SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
    m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
    lab, n = ndimage.label(m)
    if n == 0:
        return 0.0, m
    c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
    k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
    m = lab == k
    yy, xx = np.nonzero(m)
    if len(yy) < 5:
        return 0.0, m
    dy, dx = yy - yy.mean(), xx - xx.mean()                 # 2차 모멘트로 주축 각도 계산
    return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m


_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2


def warp(stack, zy, zx, ang):
    """여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
    a = np.radians(ang)[:, None, None]
    zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
    sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001)    # 출력 좌표 -> 원본 좌표
    sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
    y0, x0 = sy.astype(int), sx.astype(int)
    fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
    flat = stack.reshape(len(stack), stack.shape[1], -1)
    lin = (y0 * N + x0).reshape(len(stack), 1, -1)
    g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
    out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
    return out.reshape(stack.shape)


def preprocess_batch(O, S, M):
    """여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
    2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
    src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
    if np.isnan(src).any():                                                               # 0) 센서 결측 픽셀 보정
        mean = np.nanmean(src, axis=(2, 3), keepdims=True)
        src = np.where(np.isnan(src), mean, src)
    zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
    zx = np.array([m["gsd"] / STD_GSD for m in M])
    out = warp(src, zy, zx, np.zeros(len(M)))
    sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
    angs, cloud = np.zeros(len(M)), np.zeros(len(M))
    for i in range(len(M)):
        cm = cloud_mask(out[i, 0], sm[i])
        if cm.any():
            fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
            out[i, 0][cm] = fill
            back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
            src[i, 0][back] = fill                                                       # 원본의 같은 구름 영역도 채움
            cloud[i] = cm.mean()
        angs[i] = orientation(out[i, 0], sm[i])[0]
    rot = np.abs(angs) > 3                                                                # 3도 이하는 회전 생략
    if rot.any():
        out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
    return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]


def preprocess(opt, sar, meta):
    """한 장 전처리(preprocess_batch와 같은 처리)."""
    return preprocess_batch([opt], [sar], [meta])[0]


def extent(s):
    """평활화한 SAR 반사 마스크의 가로·세로 길이(픽셀): 객체 크기·방위 확인용."""
    m = ndimage.uniform_filter(s, 3) > 0.2
    return int(m.any(axis=0).sum()), int(m.any(axis=1).sum())


O, S, M = make_dataset(400, seed=7, scen=tuple(SCEN))
rng = np.random.default_rng(2)
for i in rng.choice(len(O), 20, replace=False):              # 센서 결측(줄 단위 NaN) 20장 주입
    r = int(rng.integers(0, N - 3)); O[i, r:r + 3] = np.nan

# 1) 처리 속도와 오류율(같은 400장을 3회 일괄 처리해 가장 빠른 회차로 처리량을 잰다: 다른 프로그램 부하 영향 최소화)
best, fails, out = 1e9, 0, []
for rep in range(3):
    t0, f = time.perf_counter(), 0
    try:
        res = preprocess_batch(O, S, M)                     # 400장을 한 번에(벡터 연산) 처리
    except Exception:
        res, f = [], len(O)
    best, fails, out = min(best, time.perf_counter() - t0), max(fails, f), res
dt = best
nan_left = sum(int(np.isnan(p[0]).sum() + np.isnan(p[1]).sum()) for p in out)
print(f"전처리 {len(O)}장, {len(O) / dt:,.0f}장/초 (목표 1,000건/초 이상), 처리 실패 {fails}건, 남은 결측 {nan_left}픽셀 -> 오류율 {fails / len(O):.0%}")
print(f"결측 주입 20장 결측치 보정율 {100 if nan_left == 0 else 0}%")

# 2) 조건별 보정 효과
sc = np.array([m["scen"] for m in M])
before_cloud = np.array([(o > 185).mean() for o in np.nan_to_num(O, nan=115)])
after_cloud = np.array([(p[0] > 185 / 255).mean() for p in out])
ang = np.array([orientation(p[0] * 255, ndimage.uniform_filter(p[1] * 1.5, 3))[0] for p in out])
horiz = (np.abs(ang) < 12) | (np.abs(ang) > 168)
ext = np.array([extent(p[1]) for p in out])
rows = []
for k in SCEN:
    sel = sc == k
    rows.append({"조건": k, "밝은 구름 픽셀(전)": f"{before_cloud[sel].mean():.1%}", "(후)": f"{after_cloud[sel].mean():.1%}",
                 "방위 정규화 성공": f"{horiz[sel].mean():.0%}", "객체 길이/폭(px)": f"{ext[sel, 0].mean():.1f}/{ext[sel, 1].mean():.1f}"})
print(pd.DataFrame(rows).to_string(index=False))

# 3) 단위 테스트 5종
raw_len = np.array([extent(s / 1.5)[0] for s in S])
tests = {
    "출력 크기 32x32(광학·SAR)": all(p[0].shape == (N, N) and p[1].shape == (N, N) for p in out),
    "값 범위 0~1, 결측 0": all(p[0].min() >= 0 and p[0].max() <= 1 and not np.isnan(p[0]).any() for p in out),
    "구름 조건 밝은 픽셀 50% 이상 감소": after_cloud[sc == "cloud"].mean() < 0.5 * before_cloud[sc == "cloud"].mean(),
    "맑은 조건 방위 정규화 90% 이상": horiz[sc == "clean"].mean() >= 0.9,
    "저해상도 객체 길이 표준화(맑은 조건 대비 30% 이내)": abs(ext[sc == "low_res", 0].mean() / ext[sc == "clean", 0].mean() - 1) < 0.3,
}
for k, v in tests.items():
    print(f"[{'PASS' if v else 'FAIL'}] {k}")
print(f"단위 테스트 {sum(tests.values())}/5 통과")
print(f"참고: 저해상도 원본의 객체 길이 {raw_len[sc == 'low_res'].mean():.1f}px -> 보정 후 {ext[sc == 'low_res', 0].mean():.1f}px")

# 4) 전후 비교 그림
fig, ax = plt.subplots(2, 5, figsize=(14, 6))
for j, k in enumerate(["clean", "cloud", "off_nadir", "low_res", "camo"]):
    i = int(np.where(sc == k)[0][0])
    ax[0, j].imshow(np.nan_to_num(O[i], nan=0), cmap="gray", vmin=0, vmax=255); ax[0, j].set_title(f"{k} raw ({M[i]['cls']})")
    ax[1, j].imshow(out[i][0], cmap="gray", vmin=0, vmax=1); ax[1, j].set_title("preprocessed")
for a in ax.ravel():
    a.axis("off")
plt.suptitle("Synthetic satellite chips: before / after preprocessing (training data)"); plt.tight_layout(); plt.show()

실행 결과 예시

전처리 400장, 1,395장/초 (목표 1,000건/초 이상), 처리 실패 0건, 남은 결측 0픽셀 -> 오류율 0%
결측 주입 20장 결측치 보정율 100%
       조건 밝은 구름 픽셀(전)  (후) 방위 정규화 성공 객체 길이/폭(px)
    clean        0.2% 0.1%       98%    14.8/8.1
     camo        0.4% 0.3%      100%    14.8/8.1
   occl30        0.2% 0.1%      100%    14.8/8.0
   occl50        0.1% 0.1%      100%    14.7/8.1
    cloud       25.1% 2.1%      100%    14.7/8.2
  low_sun        0.3% 0.2%       95%    14.4/8.5
off_nadir        0.3% 0.2%      100%    14.7/8.1
  low_res        0.0% 0.0%       95%    14.2/9.3
  speckle        0.3% 0.3%      100%    15.4/8.9
camo_occl        0.1% 0.0%       98%    14.6/8.2
[PASS] 출력 크기 32x32(광학·SAR)
[PASS] 값 범위 0~1, 결측 0
[PASS] 구름 조건 밝은 픽셀 50% 이상 감소
[PASS] 맑은 조건 방위 정규화 90% 이상
[PASS] 저해상도 객체 길이 표준화(맑은 조건 대비 30% 이내)
단위 테스트 5/5 통과
참고: 저해상도 원본의 객체 길이 5.6px -> 보정 후 14.2px

미션 3. 위성영상 내 기갑/포병 장비 객체 탐지 Baseline CNN/Vision 모델 구현

  • 난이도: 중급 / 예상 소요시간: 2.5시간
  • 코드 설명: 위성영상 조각에서 전차·자주포·트럭·민간 차량을 1차 식별하는 Baseline 분류 모델을 구현한다. 합성곱 필터 뱅크(Sobel·Laplacian·대각선·평균) -> ReLU -> 최대 풀링으로 만든 간이 CNN 특징에 로지스틱 회귀 분류기를 붙이고, 정확도·혼동행렬·판단 지연시간을 측정한다.
"""위성영상 조각에서 전차·자주포·트럭·민간 차량을 1차 식별하는 Baseline 분류 모델을 구현한다. 합성곱 필터 뱅크(Sobel·Laplacian·대각선·평균) -> ReLU -> 최대 풀링으로 만든 간이 CNN 특징에 로지스틱 회귀 분류기를 붙이고, 정확도·혼동행렬·판단 지연시간을 측정한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5                                   # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
    "tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
    "spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
    "truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
    "civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
    "decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
    "none": [],                                 # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12}   # 금속 반사 세기(가상)


def _mask(cls, angle, scale, squash=1.0):
    """회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
    yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
    yy = yy / squash
    a = np.radians(angle)
    u = (xx * np.cos(a) + yy * np.sin(a)) / scale
    v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
    return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]


def make_chip(cls, rng, cond=None):
    """chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
    c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
    gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
    scale = STD_GSD / gsd
    angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
    squash = np.cos(np.radians(off))
    sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
    sun_az = rng.uniform(0, 360)
    bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
    parts = _mask(cls, angle, scale, squash)
    body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
    opt = bg.copy()
    # 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
    L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
    sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
    opt[sh & ~body] -= 35 if cls != "decoy" else 25
    if cls == "civil":
        base = rng.choice([200.0, 60.0])            # 민간 차량 색상: 밝거나 어두움
        opt[body] = base + rng.normal(0, 8, body.sum())
        roof = parts[1][0]
        opt[roof] = base * 0.8
    else:
        ref = bg[body].mean() if body.any() else bg.mean()
        for m, d in parts:
            opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
    # 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
    if c["camo"] > 0 and cls != "civil":
        net = ndimage.binary_dilation(body, iterations=2)
        tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
        opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
    # 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
    if c["occl"] > 0:
        trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
        near = ndimage.binary_dilation(body, iterations=3)
        cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
        opt[cover] = 70 + rng.normal(0, 6, cover.sum())
    # 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
    cmask = np.zeros((N, N), bool)
    if c["cloud"] > 0:
        blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
        cmask = blob > np.quantile(blob, 1 - c["cloud"])
        opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
    # SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
    sar = np.full((N, N), 0.12)
    sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
    edge = body & ~ndimage.binary_erosion(body)
    sar[edge] *= 1.3 if cls in MILITARY else 1.0
    sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
    if c["low_res"]:                                # 저해상도: 정보 손실(블러)
        down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
        opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
    meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
            "sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
            "cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
    return np.clip(opt, 0, 255), sar, meta


SCEN = {  # 은폐·촬영 조건 시나리오 10종
    "clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
    "low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
    "speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}


def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
    rng = np.random.default_rng(seed)
    O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
    for i in range(n):
        cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
        O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
        m.update({"id": f"IMG{i:04d}", "scen": sc})
        M.append(m)
    return O, S, M
# ---------------------------------------------------------------------------

# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
    """중앙 기준으로 32x32에 맞게 자르거나 채운다."""
    out = np.full((N, N), fill, float)
    h, w = a.shape
    ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
    yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
    hh, ww = min(N, h), min(N, w)
    out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
    return out


def cloud_mask(opt, sm, min_area=15):
    """광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
    bright = opt > 185
    if bright.sum() < min_area:
        return np.zeros_like(bright)
    local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
    cand = bright & (local_sd < 12) & (sm < 0.3)
    lab, n = ndimage.label(cand)
    if n == 0:
        return cand
    sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
    return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)


def orientation(opt, sm):
    """SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
    m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
    lab, n = ndimage.label(m)
    if n == 0:
        return 0.0, m
    c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
    k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
    m = lab == k
    yy, xx = np.nonzero(m)
    if len(yy) < 5:
        return 0.0, m
    dy, dx = yy - yy.mean(), xx - xx.mean()                 # 2차 모멘트로 주축 각도 계산
    return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m


_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2


def warp(stack, zy, zx, ang):
    """여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
    a = np.radians(ang)[:, None, None]
    zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
    sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001)    # 출력 좌표 -> 원본 좌표
    sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
    y0, x0 = sy.astype(int), sx.astype(int)
    fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
    flat = stack.reshape(len(stack), stack.shape[1], -1)
    lin = (y0 * N + x0).reshape(len(stack), 1, -1)
    g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
    out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
    return out.reshape(stack.shape)


def preprocess_batch(O, S, M):
    """여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
    2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
    src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
    if np.isnan(src).any():                                                               # 0) 센서 결측 픽셀 보정
        mean = np.nanmean(src, axis=(2, 3), keepdims=True)
        src = np.where(np.isnan(src), mean, src)
    zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
    zx = np.array([m["gsd"] / STD_GSD for m in M])
    out = warp(src, zy, zx, np.zeros(len(M)))
    sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
    angs, cloud = np.zeros(len(M)), np.zeros(len(M))
    for i in range(len(M)):
        cm = cloud_mask(out[i, 0], sm[i])
        if cm.any():
            fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
            out[i, 0][cm] = fill
            back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
            src[i, 0][back] = fill                                                       # 원본의 같은 구름 영역도 채움
            cloud[i] = cm.mean()
        angs[i] = orientation(out[i, 0], sm[i])[0]
    rot = np.abs(angs) > 3                                                                # 3도 이하는 회전 생략
    if rot.any():
        out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
    return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]


def preprocess(opt, sar, meta):
    """한 장 전처리(preprocess_batch와 같은 처리)."""
    return preprocess_batch([opt], [sar], [meta])[0]

# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
    [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
    [[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]


def conv_features(o):
    """간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
    feats = []
    for k in KERNELS:
        r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
        feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
        feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
    return np.concatenate(feats)


def hog(o, cell=8, bins=8):
    gy, gx = np.gradient(o)
    mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
    b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
    out = np.zeros((N // cell, N // cell, bins))
    for i in range(N // cell):
        for j in range(N // cell):
            sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
            out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
    return (out / (np.linalg.norm(out) + 1e-6)).ravel()


def shape_features(o, s):
    """객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
    방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
    sm = ndimage.uniform_filter(s, 3)
    dark = np.abs(o - np.median(o))
    m = (sm > 0.22) | (dark > 0.12)
    lab, n = ndimage.label(m)
    if n:
        c = lab[10:22, 8:24]
        k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
        m = lab == k
    p = m.sum(axis=0).astype(float)
    cols = np.nonzero(p)[0]
    if len(cols) == 0:
        return np.zeros(8), o, s
    thin = (p > 0) & (p <= 2)
    left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
    if left > right:
        o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
        cols = np.nonzero(p)[0]
    body = np.nonzero(p > 3)[0]
    d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
    heavy = d.sum(axis=0)
    centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
    turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
    return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
                     m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s


def adv_features(o, s):
    """정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
    shp, o, s = shape_features(o, s)
    c = s[8:24, 4:28]
    sm = ndimage.uniform_filter(s, 3)
    m = sm > 0.2
    dark = np.abs(o - np.median(o))
    return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
                           sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
                           dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
                           [o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])


def features(O, S, M):
    return np.array([conv_features(p[0]) for p in preprocess_batch(O, S, M)])


# 1) 기본 데이터셋: 은폐 없는 영상(학습 800장: 맑은 조건 / 시험 400장: 맑은 조건 + 낮은 태양고도 + 비스듬 촬영)
Otr, Str, Mtr = make_dataset(800, seed=1)
Ote, Ste, Mte = make_dataset(400, seed=2, scen=("clean", "clean", "low_sun", "off_nadir"))
Xtr, ytr = features(Otr, Str, Mtr), [m["cls"] for m in Mtr]
t0 = time.perf_counter()
Xte = features(Ote, Ste, Mte)
yte = np.array([m["cls"] for m in Mte])

# 2) 규제 강도(C) 1차 튜닝: 학습 자료의 일부(검증용 200장)로 선택
best_c, best = None, -1
for c in [0.01, 0.05, 0.2, 1.0]:
    mdl = make_pipeline(StandardScaler(), LogisticRegression(C=c, max_iter=3000)).fit(Xtr[:600], ytr[:600])
    acc = mdl.score(Xtr[600:], ytr[600:])
    print(f"C={c:<5} 검증 정확도 {acc:.3f}")
    if acc > best:
        best_c, best = c, acc
model = make_pipeline(StandardScaler(), LogisticRegression(C=best_c, max_iter=3000)).fit(Xtr, ytr)
pred = model.predict(Xte)
lat = (time.perf_counter() - t0) / len(Xte) * 1000
acc = (pred == yte).mean()
print(f"\n선택 C={best_c}, 기본 시험셋 정확도 {acc:.1%} (목표 75% 이상: {'달성' if acc >= 0.75 else '미달'}), 1장당 판단 {lat:.2f} ms(전처리 포함)")

print("\n[혼동행렬] 행=실제, 열=예측")
print(pd.DataFrame(confusion_matrix(yte, pred, labels=CLASSES), index=CLASSES, columns=CLASSES).to_string())
mil, pm = np.isin(yte, list(MILITARY)), np.isin(pred, list(MILITARY))
print(f"군용 장비 여부 판정 정확도 {(mil == pm).mean():.1%}, 클래스별 재현율:",
      {c: round(float((pred[yte == c] == c).mean()), 2) for c in CLASSES})

# 3) 한계: 위장막·수목 차폐·저해상도 등 은폐 조건에서는 정확도가 떨어진다(미션 4의 동기)
Oh, Sh, Mh = make_dataset(400, seed=3, scen=("camo", "occl50", "low_res", "camo_occl"))
ph = model.predict(features(Oh, Sh, Mh))
yh, sh = np.array([m["cls"] for m in Mh]), np.array([m["scen"] for m in Mh])
print("[은폐 조건 정확도]", {k: round(float((ph[sh == k] == yh[sh == k]).mean()), 2) for k in ("camo", "occl50", "low_res", "camo_occl")})
print("판단 근거: 필터 반응(윤곽·포신 같은 가는 선)의 위치 분포. 위장막이 윤곽을 지우면 근거가 사라진다.")

# 4) 그림: 클래스별 예시와 필터 반응
fig, ax = plt.subplots(2, 4, figsize=(12, 6))
for j, c in enumerate(CLASSES):
    i = int(np.where(yte == c)[0][0])
    o = preprocess(Ote[i], Ste[i], Mte[i])[0]
    ax[0, j].imshow(o, cmap="gray"); ax[0, j].set_title(f"{c} -> pred {pred[i]}")
    ax[1, j].imshow(np.maximum(ndimage.correlate(o - o.mean(), KERNELS[0] / 8), 0), cmap="magma"); ax[1, j].set_title("Sobel-x response")
for a in ax.ravel():
    a.axis("off")
plt.suptitle("Baseline conv-feature classifier (synthetic chips)"); plt.tight_layout(); plt.show()

실행 결과 예시

C=0.01  검증 정확도 1.000
C=0.05  검증 정확도 1.000
C=0.2   검증 정확도 1.000
C=1.0   검증 정확도 1.000

선택 C=0.01, 기본 시험셋 정확도 96.5% (목표 75% 이상: 달성), 1장당 판단 5.38 ms(전처리 포함)

[혼동행렬] 행=실제, 열=예측
       tank  spg  truck  civil
tank     91    8      1      0
spg       1   99      0      0
truck     1    3     96      0
civil     0    0      0    100
군용 장비 여부 판정 정확도 100.0%, 클래스별 재현율: {'tank': 0.91, 'spg': 0.99, 'truck': 0.96, 'civil': 1.0}
[은폐 조건 정확도] {'camo': 0.5, 'occl50': 0.52, 'low_res': 0.25, 'camo_occl': 0.61}
판단 근거: 필터 반응(윤곽·포신 같은 가는 선)의 위치 분포. 위장막이 윤곽을 지우면 근거가 사라진다.

미션 4. 위장막/차폐 환경 하 군용 장비 미세 특징 추출 및 정밀 객체분류 모델 고도화

  • 난이도: 중급 / 예상 소요시간: 3.0시간
  • 코드 설명: 위장막·수목 차폐·구름·저해상도·SAR 잡음 등 은폐 조건 10종에서 군용 장비를 놓치지 않도록, 윤곽(HOG)·SAR 반사 통계·투영 윤곽 특징과 은폐 조건 증강 학습을 결합한 정밀 객체분류 모델(랜덤 포레스트)로 고도화하고 조건별 재현율과 처리 지연을 측정한다.
"""위장막·수목 차폐·구름·저해상도·SAR 잡음 등 은폐 조건 10종에서 군용 장비를 놓치지 않도록, 윤곽(HOG)·SAR 반사 통계·투영 윤곽 특징과 은폐 조건 증강 학습을 결합한 정밀 객체분류 모델(랜덤 포레스트)로 고도화하고 조건별 재현율과 처리 지연을 측정한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5                                   # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
    "tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
    "spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
    "truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
    "civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
    "decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
    "none": [],                                 # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12}   # 금속 반사 세기(가상)


def _mask(cls, angle, scale, squash=1.0):
    """회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
    yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
    yy = yy / squash
    a = np.radians(angle)
    u = (xx * np.cos(a) + yy * np.sin(a)) / scale
    v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
    return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]


def make_chip(cls, rng, cond=None):
    """chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
    c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
    gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
    scale = STD_GSD / gsd
    angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
    squash = np.cos(np.radians(off))
    sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
    sun_az = rng.uniform(0, 360)
    bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
    parts = _mask(cls, angle, scale, squash)
    body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
    opt = bg.copy()
    # 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
    L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
    sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
    opt[sh & ~body] -= 35 if cls != "decoy" else 25
    if cls == "civil":
        base = rng.choice([200.0, 60.0])            # 민간 차량 색상: 밝거나 어두움
        opt[body] = base + rng.normal(0, 8, body.sum())
        roof = parts[1][0]
        opt[roof] = base * 0.8
    else:
        ref = bg[body].mean() if body.any() else bg.mean()
        for m, d in parts:
            opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
    # 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
    if c["camo"] > 0 and cls != "civil":
        net = ndimage.binary_dilation(body, iterations=2)
        tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
        opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
    # 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
    if c["occl"] > 0:
        trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
        near = ndimage.binary_dilation(body, iterations=3)
        cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
        opt[cover] = 70 + rng.normal(0, 6, cover.sum())
    # 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
    cmask = np.zeros((N, N), bool)
    if c["cloud"] > 0:
        blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
        cmask = blob > np.quantile(blob, 1 - c["cloud"])
        opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
    # SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
    sar = np.full((N, N), 0.12)
    sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
    edge = body & ~ndimage.binary_erosion(body)
    sar[edge] *= 1.3 if cls in MILITARY else 1.0
    sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
    if c["low_res"]:                                # 저해상도: 정보 손실(블러)
        down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
        opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
    meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
            "sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
            "cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
    return np.clip(opt, 0, 255), sar, meta


SCEN = {  # 은폐·촬영 조건 시나리오 10종
    "clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
    "low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
    "speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}


def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
    rng = np.random.default_rng(seed)
    O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
    for i in range(n):
        cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
        O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
        m.update({"id": f"IMG{i:04d}", "scen": sc})
        M.append(m)
    return O, S, M
# ---------------------------------------------------------------------------

# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
    """중앙 기준으로 32x32에 맞게 자르거나 채운다."""
    out = np.full((N, N), fill, float)
    h, w = a.shape
    ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
    yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
    hh, ww = min(N, h), min(N, w)
    out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
    return out


def cloud_mask(opt, sm, min_area=15):
    """광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
    bright = opt > 185
    if bright.sum() < min_area:
        return np.zeros_like(bright)
    local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
    cand = bright & (local_sd < 12) & (sm < 0.3)
    lab, n = ndimage.label(cand)
    if n == 0:
        return cand
    sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
    return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)


def orientation(opt, sm):
    """SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
    m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
    lab, n = ndimage.label(m)
    if n == 0:
        return 0.0, m
    c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
    k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
    m = lab == k
    yy, xx = np.nonzero(m)
    if len(yy) < 5:
        return 0.0, m
    dy, dx = yy - yy.mean(), xx - xx.mean()                 # 2차 모멘트로 주축 각도 계산
    return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m


_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2


def warp(stack, zy, zx, ang):
    """여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
    a = np.radians(ang)[:, None, None]
    zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
    sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001)    # 출력 좌표 -> 원본 좌표
    sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
    y0, x0 = sy.astype(int), sx.astype(int)
    fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
    flat = stack.reshape(len(stack), stack.shape[1], -1)
    lin = (y0 * N + x0).reshape(len(stack), 1, -1)
    g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
    out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
    return out.reshape(stack.shape)


def preprocess_batch(O, S, M):
    """여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
    2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
    src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
    if np.isnan(src).any():                                                               # 0) 센서 결측 픽셀 보정
        mean = np.nanmean(src, axis=(2, 3), keepdims=True)
        src = np.where(np.isnan(src), mean, src)
    zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
    zx = np.array([m["gsd"] / STD_GSD for m in M])
    out = warp(src, zy, zx, np.zeros(len(M)))
    sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
    angs, cloud = np.zeros(len(M)), np.zeros(len(M))
    for i in range(len(M)):
        cm = cloud_mask(out[i, 0], sm[i])
        if cm.any():
            fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
            out[i, 0][cm] = fill
            back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
            src[i, 0][back] = fill                                                       # 원본의 같은 구름 영역도 채움
            cloud[i] = cm.mean()
        angs[i] = orientation(out[i, 0], sm[i])[0]
    rot = np.abs(angs) > 3                                                                # 3도 이하는 회전 생략
    if rot.any():
        out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
    return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]


def preprocess(opt, sar, meta):
    """한 장 전처리(preprocess_batch와 같은 처리)."""
    return preprocess_batch([opt], [sar], [meta])[0]

# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
    [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
    [[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]


def conv_features(o):
    """간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
    feats = []
    for k in KERNELS:
        r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
        feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
        feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
    return np.concatenate(feats)


def hog(o, cell=8, bins=8):
    gy, gx = np.gradient(o)
    mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
    b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
    out = np.zeros((N // cell, N // cell, bins))
    for i in range(N // cell):
        for j in range(N // cell):
            sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
            out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
    return (out / (np.linalg.norm(out) + 1e-6)).ravel()


def shape_features(o, s):
    """객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
    방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
    sm = ndimage.uniform_filter(s, 3)
    dark = np.abs(o - np.median(o))
    m = (sm > 0.22) | (dark > 0.12)
    lab, n = ndimage.label(m)
    if n:
        c = lab[10:22, 8:24]
        k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
        m = lab == k
    p = m.sum(axis=0).astype(float)
    cols = np.nonzero(p)[0]
    if len(cols) == 0:
        return np.zeros(8), o, s
    thin = (p > 0) & (p <= 2)
    left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
    if left > right:
        o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
        cols = np.nonzero(p)[0]
    body = np.nonzero(p > 3)[0]
    d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
    heavy = d.sum(axis=0)
    centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
    turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
    return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
                     m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s


def adv_features(o, s):
    """정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
    shp, o, s = shape_features(o, s)
    c = s[8:24, 4:28]
    sm = ndimage.uniform_filter(s, 3)
    m = sm > 0.2
    dark = np.abs(o - np.median(o))
    return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
                           sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
                           dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
                           [o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])


def prep(O, S, M):
    return preprocess_batch(O, S, M)


ALL = tuple(SCEN)
# 1) 학습: 은폐 조건 10종을 고르게 섞은 증강 학습 자료 1,200장
Otr, Str, Mtr = make_dataset(1200, seed=3, scen=ALL)
Ptr = prep(Otr, Str, Mtr); ytr = [m["cls"] for m in Mtr]
adv = RandomForestClassifier(n_estimators=200, min_samples_leaf=2, random_state=0, n_jobs=1)
adv.fit(np.array([adv_features(o, s) for o, s, _ in Ptr]), ytr)
base = make_pipeline(StandardScaler(), LogisticRegression(C=0.05, max_iter=3000))     # 비교용 Baseline(맑은 조건만 학습)
Oc, Sc, Mc = make_dataset(600, seed=1)
base.fit(np.array([conv_features(p[0]) for p in prep(Oc, Sc, Mc)]), [m["cls"] for m in Mc])

# 2) 시험: 조건 10종 x 4클래스 x 15장 = 600장
Ote, Ste, Mte = make_dataset(600, seed=4, scen=ALL)
t0 = time.perf_counter()
Pte = prep(Ote, Ste, Mte)
pa = adv.predict(np.array([adv_features(o, s) for o, s, _ in Pte]))
lat = (time.perf_counter() - t0) / len(Pte) * 1000
pb = base.predict(np.array([conv_features(p[0]) for p in Pte]))
y, sc = np.array([m["cls"] for m in Mte]), np.array([m["scen"] for m in Mte])
mil = np.isin(y, list(MILITARY))
rec = lambda p, sel: ((p == y) & mil & sel).sum() / (mil & sel).sum()      # 군용 장비를 '정확한 기종'으로 식별한 비율
allsel = np.ones(len(y), bool)
print(f"군용 장비 기종 식별 재현율: Baseline {rec(pb, allsel):.1%} -> 정밀 모델 {rec(pa, allsel):.1%} (목표 85% 이상: {'달성' if rec(pa, allsel) >= 0.85 else '미달'})")
print(f"4클래스 정확도: Baseline {(pb == y).mean():.1%} -> 정밀 모델 {(pa == y).mean():.1%}, 1장당 처리 {lat:.1f} ms (목표 200 ms 이하)")

# 3) 은폐 조건 10종 시험(군용 기종 식별 재현율 85% 이상이면 통과)
rows = []
for k in ALL:
    s = sc == k
    rows.append({"조건": k, "Baseline 재현율": round(rec(pb, s), 2), "정밀 재현율": round(rec(pa, s), 2),
                 "정밀 정확도": round(float((pa[s] == y[s]).mean()), 2), "민간->군용 오분류": int((~mil & s & np.isin(pa, list(MILITARY))).sum()),
                 "판정": "통과" if rec(pa, s) >= 0.85 else "미통과"})
t = pd.DataFrame(rows)
print("\n" + t.to_string(index=False))
print(f"통과 {int((t['판정'] == '통과').sum())}/10종")

# 4) 오경보 유발 패턴과 중요 특징
print("\n[혼동이 많은 쌍(실제->예측)]", pd.Series([f"{a}->{b}" for a, b in zip(y, pa) if a != b]).value_counts().head(4).to_dict())
names = ["obj_len", "body_len", "thin_cols", "max_width", "turret_pos", "area", "barrel_len", "fill"] + [f"hog{i}" for i in range(128)] + ["sar_mean", "sar_max", "sar>0.25", "sar>0.4", "sar_width", "sar_len"] + \
        [f"sar_colprof{i}" for i in range(8)] + [f"sar_rowprof{i}" for i in range(8)] + [f"opt_colprof{i}" for i in range(8)] + \
        [f"opt_rowprof{i}" for i in range(8)] + ["opt_mean", "opt_std", "opt_contrast", "bright_frac", "dark_frac"]
imp = pd.Series(adv.feature_importances_, index=names).sort_values(ascending=False)
print("중요 특징 상위 5:", {k: round(v, 3) for k, v in imp.head(5).items()})
print("-> SAR 반사 특징이 위장막·구름에 강해 재현율을 지탱한다. 남은 오류는 위 혼동 쌍처럼 비슷한 크기 기종 간 혼동이며 미션 5에서 신뢰도로 거른다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
x = np.arange(len(ALL))
ax[0].bar(x - 0.2, t["Baseline 재현율"], 0.4, label="baseline"); ax[0].bar(x + 0.2, t["정밀 재현율"], 0.4, label="advanced")
ax[0].axhline(0.85, ls="--", c="k"); ax[0].set_xticks(x); ax[0].set_xticklabels(ALL, rotation=45, fontsize=8)
ax[0].set_ylabel("military recall"); ax[0].legend(); ax[0].set_title("Recall by concealment scenario")
imp.head(10)[::-1].plot.barh(ax=ax[1]); ax[1].set_title("Top feature importances")
plt.tight_layout(); plt.show()

실행 결과 예시

군용 장비 기종 식별 재현율: Baseline 66.0% -> 정밀 모델 97.3% (목표 85% 이상: 달성)
4클래스 정확도: Baseline 69.2% -> 정밀 모델 98.0%, 1장당 처리 1.5 ms (목표 200 ms 이하)

       조건  Baseline 재현율  정밀 재현율  정밀 정확도  민간->군용 오분류 판정
    clean          1.00    0.98    0.98           0 통과
     camo          0.33    0.96    0.97           0 통과
   occl30          0.73    1.00    1.00           0 통과
   occl50          0.53    0.98    0.98           0 통과
    cloud          0.33    1.00    1.00           0 통과
  low_sun          0.82    0.96    0.97           0 통과
off_nadir          0.98    1.00    1.00           0 통과
  low_res          0.33    0.87    0.90           0 통과
  speckle          1.00    1.00    1.00           0 통과
camo_occl          0.53    1.00    1.00           0 통과
통과 10/10종

[혼동이 많은 쌍(실제->예측)] {'truck->tank': 6, 'tank->truck': 5, 'tank->spg': 1}
중요 특징 상위 5: {'sar>0.25': 0.083, 'sar_width': 0.08, 'sar_len': 0.07, 'sar_colprof2': 0.064, 'sar_colprof4': 0.057}
-> SAR 반사 특징이 위장막·구름에 강해 재현율을 지탱한다. 남은 오류는 위 혼동 쌍처럼 비슷한 크기 기종 간 혼동이며 미션 5에서 신뢰도로 거른다.

미션 5. 위성 촬영 각도/그림자 영향에 따른 객체분류 신뢰도 점수 산출 및 오분류 필터 구현

  • 난이도: 중급~고급 / 예상 소요시간: 3.0시간
  • 코드 설명: 위성 센서 해상도 한계, 비스듬 촬영 각도, 태양 고도(그림자), 구름 차폐를 반영해 장비 분류 결과의 신뢰도 점수(분류가 맞을 확률)와 등급(High/Medium/Low)을 계산하고, 낮은 신뢰도 결과를 판독 보류로 돌려 오분류 보고를 줄이는 엔진을 개발한다.
"""위성 센서 해상도 한계, 비스듬 촬영 각도, 태양 고도(그림자), 구름 차폐를 반영해 장비 분류 결과의 신뢰도 점수(분류가 맞을 확률)와 등급(High/Medium/Low)을 계산하고, 낮은 신뢰도 결과를 판독 보류로 돌려 오분류 보고를 줄이는 엔진을 개발한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5                                   # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
    "tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
    "spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
    "truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
    "civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
    "decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
    "none": [],                                 # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12}   # 금속 반사 세기(가상)


def _mask(cls, angle, scale, squash=1.0):
    """회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
    yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
    yy = yy / squash
    a = np.radians(angle)
    u = (xx * np.cos(a) + yy * np.sin(a)) / scale
    v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
    return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]


def make_chip(cls, rng, cond=None):
    """chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
    c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
    gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
    scale = STD_GSD / gsd
    angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
    squash = np.cos(np.radians(off))
    sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
    sun_az = rng.uniform(0, 360)
    bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
    parts = _mask(cls, angle, scale, squash)
    body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
    opt = bg.copy()
    # 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
    L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
    sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
    opt[sh & ~body] -= 35 if cls != "decoy" else 25
    if cls == "civil":
        base = rng.choice([200.0, 60.0])            # 민간 차량 색상: 밝거나 어두움
        opt[body] = base + rng.normal(0, 8, body.sum())
        roof = parts[1][0]
        opt[roof] = base * 0.8
    else:
        ref = bg[body].mean() if body.any() else bg.mean()
        for m, d in parts:
            opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
    # 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
    if c["camo"] > 0 and cls != "civil":
        net = ndimage.binary_dilation(body, iterations=2)
        tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
        opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
    # 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
    if c["occl"] > 0:
        trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
        near = ndimage.binary_dilation(body, iterations=3)
        cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
        opt[cover] = 70 + rng.normal(0, 6, cover.sum())
    # 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
    cmask = np.zeros((N, N), bool)
    if c["cloud"] > 0:
        blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
        cmask = blob > np.quantile(blob, 1 - c["cloud"])
        opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
    # SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
    sar = np.full((N, N), 0.12)
    sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
    edge = body & ~ndimage.binary_erosion(body)
    sar[edge] *= 1.3 if cls in MILITARY else 1.0
    sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
    if c["low_res"]:                                # 저해상도: 정보 손실(블러)
        down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
        opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
    meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
            "sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
            "cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
    return np.clip(opt, 0, 255), sar, meta


SCEN = {  # 은폐·촬영 조건 시나리오 10종
    "clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
    "low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
    "speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}


def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
    rng = np.random.default_rng(seed)
    O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
    for i in range(n):
        cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
        O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
        m.update({"id": f"IMG{i:04d}", "scen": sc})
        M.append(m)
    return O, S, M
# ---------------------------------------------------------------------------

# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
    """중앙 기준으로 32x32에 맞게 자르거나 채운다."""
    out = np.full((N, N), fill, float)
    h, w = a.shape
    ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
    yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
    hh, ww = min(N, h), min(N, w)
    out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
    return out


def cloud_mask(opt, sm, min_area=15):
    """광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
    bright = opt > 185
    if bright.sum() < min_area:
        return np.zeros_like(bright)
    local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
    cand = bright & (local_sd < 12) & (sm < 0.3)
    lab, n = ndimage.label(cand)
    if n == 0:
        return cand
    sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
    return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)


def orientation(opt, sm):
    """SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
    m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
    lab, n = ndimage.label(m)
    if n == 0:
        return 0.0, m
    c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
    k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
    m = lab == k
    yy, xx = np.nonzero(m)
    if len(yy) < 5:
        return 0.0, m
    dy, dx = yy - yy.mean(), xx - xx.mean()                 # 2차 모멘트로 주축 각도 계산
    return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m


_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2


def warp(stack, zy, zx, ang):
    """여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
    a = np.radians(ang)[:, None, None]
    zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
    sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001)    # 출력 좌표 -> 원본 좌표
    sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
    y0, x0 = sy.astype(int), sx.astype(int)
    fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
    flat = stack.reshape(len(stack), stack.shape[1], -1)
    lin = (y0 * N + x0).reshape(len(stack), 1, -1)
    g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
    out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
    return out.reshape(stack.shape)


def preprocess_batch(O, S, M):
    """여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
    2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
    src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
    if np.isnan(src).any():                                                               # 0) 센서 결측 픽셀 보정
        mean = np.nanmean(src, axis=(2, 3), keepdims=True)
        src = np.where(np.isnan(src), mean, src)
    zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
    zx = np.array([m["gsd"] / STD_GSD for m in M])
    out = warp(src, zy, zx, np.zeros(len(M)))
    sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
    angs, cloud = np.zeros(len(M)), np.zeros(len(M))
    for i in range(len(M)):
        cm = cloud_mask(out[i, 0], sm[i])
        if cm.any():
            fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
            out[i, 0][cm] = fill
            back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
            src[i, 0][back] = fill                                                       # 원본의 같은 구름 영역도 채움
            cloud[i] = cm.mean()
        angs[i] = orientation(out[i, 0], sm[i])[0]
    rot = np.abs(angs) > 3                                                                # 3도 이하는 회전 생략
    if rot.any():
        out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
    return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]


def preprocess(opt, sar, meta):
    """한 장 전처리(preprocess_batch와 같은 처리)."""
    return preprocess_batch([opt], [sar], [meta])[0]

# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
    [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
    [[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]


def conv_features(o):
    """간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
    feats = []
    for k in KERNELS:
        r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
        feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
        feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
    return np.concatenate(feats)


def hog(o, cell=8, bins=8):
    gy, gx = np.gradient(o)
    mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
    b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
    out = np.zeros((N // cell, N // cell, bins))
    for i in range(N // cell):
        for j in range(N // cell):
            sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
            out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
    return (out / (np.linalg.norm(out) + 1e-6)).ravel()


def shape_features(o, s):
    """객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
    방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
    sm = ndimage.uniform_filter(s, 3)
    dark = np.abs(o - np.median(o))
    m = (sm > 0.22) | (dark > 0.12)
    lab, n = ndimage.label(m)
    if n:
        c = lab[10:22, 8:24]
        k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
        m = lab == k
    p = m.sum(axis=0).astype(float)
    cols = np.nonzero(p)[0]
    if len(cols) == 0:
        return np.zeros(8), o, s
    thin = (p > 0) & (p <= 2)
    left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
    if left > right:
        o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
        cols = np.nonzero(p)[0]
    body = np.nonzero(p > 3)[0]
    d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
    heavy = d.sum(axis=0)
    centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
    turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
    return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
                     m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s


def adv_features(o, s):
    """정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
    shp, o, s = shape_features(o, s)
    c = s[8:24, 4:28]
    sm = ndimage.uniform_filter(s, 3)
    m = sm > 0.2
    dark = np.abs(o - np.median(o))
    return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
                           sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
                           dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
                           [o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])

Q_NAMES = ["max_prob", "margin", "entropy", "cloud_frac", "gsd", "off_nadir", "sun_el", "sar_mean"]


def run(O, S, M, clf):
    P = preprocess_batch(O, S, M)
    X = np.array([adv_features(o, s) for o, s, _ in P])
    pr = clf.predict_proba(X)
    srt = np.sort(pr, axis=1)
    ent = -(pr * np.log(np.clip(pr, 1e-9, 1))).sum(1) / np.log(pr.shape[1])
    Q = np.c_[srt[:, -1], srt[:, -1] - srt[:, -2], ent, [p[2] for p in P], [m["gsd"] for m in M],
              [m["off_nadir"] for m in M], [m["sun_el"] for m in M], [p[1][8:24, 4:28].mean() for p in P]]
    return clf.classes_[pr.argmax(1)], Q


def grade(c):
    return np.where(c >= 0.9, "High", np.where(c >= 0.6, "Medium", "Low"))


ALL = tuple(SCEN)
Otr, Str, Mtr = make_dataset(1200, seed=3, scen=ALL)
clf = RandomForestClassifier(n_estimators=200, min_samples_leaf=2, random_state=0, n_jobs=1)
clf.fit(np.array([adv_features(o, s) for o, s, _ in preprocess_batch(Otr, Str, Mtr)]), [m["cls"] for m in Mtr])

# 1) 신뢰도 가중치 학습: 별도 보정 자료 600장에서 '분류가 맞았는가'를 품질 지표로 예측
Oc, Sc, Mc = make_dataset(600, seed=6, scen=ALL)
pc, Qc = run(Oc, Sc, Mc, clf)
conf_model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)).fit(Qc, pc == np.array([m["cls"] for m in Mc]))
w = pd.Series(conf_model[-1].coef_[0], index=Q_NAMES).round(2)
print("[신뢰도 가중치(표준화 계수)]", w.to_dict())

# 2) 시험 800장: 기존(분류 결과 그대로 보고) vs 엔진(High만 자동 보고, 나머지는 판독 보류)
Ot, St, Mt = make_dataset(800, seed=5, scen=ALL)
pt, Qt = run(Ot, St, Mt, clf)
y = np.array([m["cls"] for m in Mt])
conf = conf_model.predict_proba(Qt)[:, 1]
G = grade(conf)
wrong = pt != y
rep0, rep1 = np.isin(pt, list(MILITARY)), np.isin(pt, list(MILITARY)) & (G == "High")
far0, far1 = (wrong & rep0).sum() / rep0.sum(), (wrong & rep1).sum() / rep1.sum()
print(f"\n기존: 군용 보고 {rep0.sum()}건 중 오분류 {(wrong & rep0).sum()}건 (오분류율 {far0:.1%})")
print(f"엔진: 군용 자동 보고(High) {rep1.sum()}건 중 오분류 {(wrong & rep1).sum()}건 (오분류율 {far1:.1%}), 판독 보류(Medium·Low) {(rep0 & ~rep1).sum()}건")
cut = 1 - far1 / far0
print(f"오분류(오경보) 감축 {cut:.1%} (목표 40% 이상: {'달성' if cut >= 0.4 else '미달'}), ROC-AUC(신뢰도 vs 정답 여부) {roc_auc_score(~wrong, conf):.3f}")

# 3) 등급별 실제 정답률과 임계값 민감도
sc = np.array([m["scen"] for m in Mt])
print("\n등급   건수  실제 정답률  주요 조건")
for g in ["High", "Medium", "Low"]:
    s = G == g
    print(f"{g:6s} {s.sum():4d}   {(~wrong[s]).mean() if s.any() else 0:9.1%}   {pd.Series(sc[s]).value_counts().head(2).to_dict()}")
print("\nHigh 기준  자동보고  오분류율  판독 보류")
for t in [0.8, 0.85, 0.9, 0.95]:
    r = rep0 & (conf >= t)
    print(f"   {t:.2f}    {r.sum():6d}   {(wrong & r).sum() / max(r.sum(), 1):6.1%}   {(rep0 & ~r).sum():6d}")

# 4) 일관성: 같은 영상에 미세한 픽셀 잡음(+-1%)을 넣어도 같은 등급이 나오는 비율
rng = np.random.default_rng(0)
sub = np.arange(200)
_, Qn = run(Ot[sub] * (1 + rng.normal(0, 0.01, Ot[sub].shape)), St[sub], [Mt[i] for i in sub], clf)
same = (grade(conf_model.predict_proba(Qn)[:, 1]) == G[sub]).mean()
print(f"\n신뢰도 등급 산출 일관성 {same:.1%} (목표 90% 이상: {'달성' if same >= 0.9 else '미달'})")
print("High도 최종 판독이 아니며, 고가치 표적은 미션 7의 영상분석관 판독을 거친다.")

fig, ax = plt.subplots(1, 2, figsize=(12, 4))
ax[0].hist([conf[~wrong], conf[wrong]], bins=20, range=(0, 1), label=["correct", "misclassified"], stacked=True)
ax[0].axvline(0.9, ls="--", c="k"); ax[0].set_xlabel("confidence score"); ax[0].legend(); ax[0].set_title("Confidence vs correctness")
w.plot.barh(ax=ax[1], color=["tab:green" if v > 0 else "tab:red" for v in w]); ax[1].set_title("Confidence weight matrix (std. coef)")
plt.tight_layout(); plt.show()

실행 결과 예시

[신뢰도 가중치(표준화 계수)] {'max_prob': 0.69, 'margin': 0.81, 'entropy': -0.62, 'cloud_frac': 0.2, 'gsd': -0.06, 'off_nadir': -0.54, 'sun_el': -0.52, 'sar_mean': 0.24}

기존: 군용 보고 600건 중 오분류 19건 (오분류율 3.2%)
엔진: 군용 자동 보고(High) 541건 중 오분류 3건 (오분류율 0.6%), 판독 보류(Medium·Low) 59건
오분류(오경보) 감축 82.5% (목표 40% 이상: 달성), ROC-AUC(신뢰도 vs 정답 여부) 0.961

등급   건수  실제 정답률  주요 조건
High    741       99.6%   {'occl50': 79, 'clean': 78}
Medium   38       84.2%   {'low_res': 13, 'occl30': 5}
Low      21       52.4%   {'low_res': 9, 'speckle': 4}

High 기준  자동보고  오분류율  판독 보류
   0.80       555     1.1%       45
   0.85       552     0.9%       48
   0.90       541     0.6%       59
   0.95       521     0.6%       79

신뢰도 등급 산출 일관성 98.5% (목표 90% 이상: 달성)
High도 최종 판독이 아니며, 고가치 표적은 미션 7의 영상분석관 판독을 거친다.

미션 6. 위성영상 표적 탐지 위치 및 식별 확률 지도(Heatmap) 가시화 UI 구현

  • 난이도: 중급 / 예상 소요시간: 2.5시간
  • 코드 설명: 가상 광역 위성영상(32x32 타일 10x10)을 타일 단위로 판독해 식별 장비의 위치·기종·식별 확률을 지도에 겹쳐 표시하고, 군용 장비 집결지 밀집도(Heatmap)와 기종별 현황을 함께 보여 주는 위성 ISR 영상분석 화면을 구현하며 갱신 처리시간과 렌더링 시간을 측정한다.
"""가상 광역 위성영상(32x32 타일 10x10)을 타일 단위로 판독해 식별 장비의 위치·기종·식별 확률을 지도에 겹쳐 표시하고, 군용 장비 집결지 밀집도(Heatmap)와 기종별 현황을 함께 보여 주는 위성 ISR 영상분석 화면을 구현하며 갱신 처리시간과 렌더링 시간을 측정한다."""
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier

# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5                                   # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
    "tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
    "spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
    "truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
    "civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
    "decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
    "none": [],                                 # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12}   # 금속 반사 세기(가상)


def _mask(cls, angle, scale, squash=1.0):
    """회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
    yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
    yy = yy / squash
    a = np.radians(angle)
    u = (xx * np.cos(a) + yy * np.sin(a)) / scale
    v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
    return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]


def make_chip(cls, rng, cond=None):
    """chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
    c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
    gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
    scale = STD_GSD / gsd
    angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
    squash = np.cos(np.radians(off))
    sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
    sun_az = rng.uniform(0, 360)
    bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
    parts = _mask(cls, angle, scale, squash)
    body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
    opt = bg.copy()
    # 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
    L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
    sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
    opt[sh & ~body] -= 35 if cls != "decoy" else 25
    if cls == "civil":
        base = rng.choice([200.0, 60.0])            # 민간 차량 색상: 밝거나 어두움
        opt[body] = base + rng.normal(0, 8, body.sum())
        roof = parts[1][0]
        opt[roof] = base * 0.8
    else:
        ref = bg[body].mean() if body.any() else bg.mean()
        for m, d in parts:
            opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
    # 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
    if c["camo"] > 0 and cls != "civil":
        net = ndimage.binary_dilation(body, iterations=2)
        tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
        opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
    # 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
    if c["occl"] > 0:
        trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
        near = ndimage.binary_dilation(body, iterations=3)
        cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
        opt[cover] = 70 + rng.normal(0, 6, cover.sum())
    # 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
    cmask = np.zeros((N, N), bool)
    if c["cloud"] > 0:
        blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
        cmask = blob > np.quantile(blob, 1 - c["cloud"])
        opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
    # SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
    sar = np.full((N, N), 0.12)
    sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
    edge = body & ~ndimage.binary_erosion(body)
    sar[edge] *= 1.3 if cls in MILITARY else 1.0
    sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
    if c["low_res"]:                                # 저해상도: 정보 손실(블러)
        down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
        opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
    meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
            "sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
            "cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
    return np.clip(opt, 0, 255), sar, meta


SCEN = {  # 은폐·촬영 조건 시나리오 10종
    "clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
    "low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
    "speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}


def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
    rng = np.random.default_rng(seed)
    O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
    for i in range(n):
        cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
        O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
        m.update({"id": f"IMG{i:04d}", "scen": sc})
        M.append(m)
    return O, S, M
# ---------------------------------------------------------------------------

# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
    """중앙 기준으로 32x32에 맞게 자르거나 채운다."""
    out = np.full((N, N), fill, float)
    h, w = a.shape
    ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
    yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
    hh, ww = min(N, h), min(N, w)
    out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
    return out


def cloud_mask(opt, sm, min_area=15):
    """광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
    bright = opt > 185
    if bright.sum() < min_area:
        return np.zeros_like(bright)
    local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
    cand = bright & (local_sd < 12) & (sm < 0.3)
    lab, n = ndimage.label(cand)
    if n == 0:
        return cand
    sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
    return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)


def orientation(opt, sm):
    """SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
    m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
    lab, n = ndimage.label(m)
    if n == 0:
        return 0.0, m
    c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
    k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
    m = lab == k
    yy, xx = np.nonzero(m)
    if len(yy) < 5:
        return 0.0, m
    dy, dx = yy - yy.mean(), xx - xx.mean()                 # 2차 모멘트로 주축 각도 계산
    return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m


_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2


def warp(stack, zy, zx, ang):
    """여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
    a = np.radians(ang)[:, None, None]
    zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
    sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001)    # 출력 좌표 -> 원본 좌표
    sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
    y0, x0 = sy.astype(int), sx.astype(int)
    fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
    flat = stack.reshape(len(stack), stack.shape[1], -1)
    lin = (y0 * N + x0).reshape(len(stack), 1, -1)
    g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
    out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
    return out.reshape(stack.shape)


def preprocess_batch(O, S, M):
    """여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
    2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
    src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
    if np.isnan(src).any():                                                               # 0) 센서 결측 픽셀 보정
        mean = np.nanmean(src, axis=(2, 3), keepdims=True)
        src = np.where(np.isnan(src), mean, src)
    zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
    zx = np.array([m["gsd"] / STD_GSD for m in M])
    out = warp(src, zy, zx, np.zeros(len(M)))
    sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
    angs, cloud = np.zeros(len(M)), np.zeros(len(M))
    for i in range(len(M)):
        cm = cloud_mask(out[i, 0], sm[i])
        if cm.any():
            fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
            out[i, 0][cm] = fill
            back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
            src[i, 0][back] = fill                                                       # 원본의 같은 구름 영역도 채움
            cloud[i] = cm.mean()
        angs[i] = orientation(out[i, 0], sm[i])[0]
    rot = np.abs(angs) > 3                                                                # 3도 이하는 회전 생략
    if rot.any():
        out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
    return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]


def preprocess(opt, sar, meta):
    """한 장 전처리(preprocess_batch와 같은 처리)."""
    return preprocess_batch([opt], [sar], [meta])[0]

# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
    [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
    [[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]


def conv_features(o):
    """간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
    feats = []
    for k in KERNELS:
        r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
        feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
        feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
    return np.concatenate(feats)


def hog(o, cell=8, bins=8):
    gy, gx = np.gradient(o)
    mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
    b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
    out = np.zeros((N // cell, N // cell, bins))
    for i in range(N // cell):
        for j in range(N // cell):
            sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
            out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
    return (out / (np.linalg.norm(out) + 1e-6)).ravel()


def shape_features(o, s):
    """객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
    방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
    sm = ndimage.uniform_filter(s, 3)
    dark = np.abs(o - np.median(o))
    m = (sm > 0.22) | (dark > 0.12)
    lab, n = ndimage.label(m)
    if n:
        c = lab[10:22, 8:24]
        k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
        m = lab == k
    p = m.sum(axis=0).astype(float)
    cols = np.nonzero(p)[0]
    if len(cols) == 0:
        return np.zeros(8), o, s
    thin = (p > 0) & (p <= 2)
    left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
    if left > right:
        o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
        cols = np.nonzero(p)[0]
    body = np.nonzero(p > 3)[0]
    d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
    heavy = d.sum(axis=0)
    centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
    turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
    return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
                     m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s


def adv_features(o, s):
    """정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
    shp, o, s = shape_features(o, s)
    c = s[8:24, 4:28]
    sm = ndimage.uniform_filter(s, 3)
    m = sm > 0.2
    dark = np.abs(o - np.median(o))
    return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
                           sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
                           dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
                           [o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])

T = 10                                                  # 10x10 타일 = 320x320 픽셀 가상 장면
rng = np.random.default_rng(12)
# [교육용가상자료] 장면 구성: 가상 집결지 2곳(군용), 가상 도로(민간 차량), 나머지는 빈 타일
layout = np.full((T, T), "none", dtype=object)
for (r0, c0), n in [((1, 1), 9), ((6, 6), 7)]:
    for _ in range(n):
        layout[r0 + rng.integers(0, 3), c0 + rng.integers(0, 3)] = rng.choice(["tank", "spg", "truck"], p=[.4, .3, .3])
for c in rng.choice(T, 5, replace=False):
    layout[4, c] = "civil"
scen = rng.choice(["clean", "camo", "occl30", "cloud", "low_sun"], (T, T), p=[.4, .2, .15, .15, .1])
tiles = {}
for r in range(T):
    for c in range(T):
        o, s, m = make_chip(layout[r, c], rng, SCEN[scen[r, c]])
        tiles[(r, c)] = (o, s, m)

# 판독 모델(이 파일 안에서 학습: 빈 타일 포함 5개 클래스)
Otr, Str, Mtr = make_dataset(1000, seed=3, scen=tuple(SCEN), classes=CLASSES + ["none"])
clf = RandomForestClassifier(n_estimators=150, min_samples_leaf=2, random_state=0, n_jobs=1)
clf.fit(np.array([adv_features(o, s) for o, s, _ in preprocess_batch(Otr, Str, Mtr)]), [m["cls"] for m in Mtr])

# 1) 1초 주기 갱신 모의: 매 갱신마다 새로 들어온 타일 한 줄(10장)을 판독해 화면 자료 갱신
dets, upd = [], []
for r in range(T):
    t0 = time.perf_counter()
    row = [tiles[(r, c)] for c in range(T)]
    X = np.array([adv_features(o, s) for o, s, _ in preprocess_batch(*zip(*row))])
    P = clf.predict_proba(X)                            # 한 줄을 묶어서 판독(호출 횟수 감소)
    for c, p in enumerate(P):
        k = int(p.argmax())
        if clf.classes_[k] != "none":
            dets.append({"타일": f"R{r}C{c}", "x": c * N + N / 2, "y": r * N + N / 2, "기종(추정)": clf.classes_[k],
                         "식별확률": round(float(p[k]), 2), "실제": layout[r, c], "조건": scen[r, c]})
    upd.append((time.perf_counter() - t0) * 1000)
D = pd.DataFrame(dets)
print(f"타일 한 줄(10장) 갱신 처리 평균 {np.mean(upd):.0f} ms, 최대 {np.max(upd):.0f} ms (목표: 갱신 1초 이내, 처리 응답 100 ms 이내)")
print(f"식별 {len(D)}건 / 실제 객체 {int((layout != 'none').sum())}건, 기종 일치 {int((D['기종(추정)'] == D['실제']).sum())}건")
print("\n[식별 목록: 식별확률 낮은 순(분석관 우선 확인)]")
print(D.sort_values("식별확률")[["타일", "기종(추정)", "식별확률", "조건"]].head(8).to_string(index=False))

# 2) 집결지 밀집도: 군용 식별 위치를 식별확률로 가중해 가우시안 평활
mil = D[D["기종(추정)"].isin(list(MILITARY))]
grid = np.zeros((T * N, T * N))
for _, d in mil.iterrows():
    grid[int(d["y"]), int(d["x"])] += d["식별확률"]
heat = ndimage.gaussian_filter(grid, 24)
peak = np.unravel_index(heat.argmax(), heat.shape)
print(f"\n최대 밀집 지점(가상 픽셀 좌표) x={peak[1]}, y={peak[0]}, 군용 식별 {len(mil)}건")

# 3) 분석 화면(야간 작전용 어두운 화면, 4개 구성요소)
scene = np.block([[tiles[(r, c)][0] for c in range(T)] for r in range(T)])
plt.style.use("dark_background")
fig, ax = plt.subplots(2, 2, figsize=(13, 11))
mk = {"tank": "s", "spg": "D", "truck": "^", "civil": "o"}
ax[0, 0].imshow(scene, cmap="gray")
for k, g in D.groupby("기종(추정)"):
    ax[0, 0].scatter(g["x"], g["y"], s=40 + 200 * g["식별확률"], marker=mk[k], facecolors="none",
                     edgecolors="red" if k in MILITARY else "deepskyblue", lw=1.5, label=k)
ax[0, 0].legend(fontsize=8, loc="lower right"); ax[0, 0].set_title("Detections (marker size = probability)"); ax[0, 0].axis("off")
ax[0, 1].imshow(scene, cmap="gray"); ax[0, 1].imshow(heat, cmap="inferno", alpha=0.6); ax[0, 1].set_title("Military vehicle density (heatmap)"); ax[0, 1].axis("off")
vc = D["기종(추정)"].value_counts()
ax[1, 0].bar(vc.index.tolist(), vc.values, color="deepskyblue"); ax[1, 0].set_title("Detections by type")
ax[1, 1].hist(D["식별확률"], bins=10, range=(0, 1), color="orange"); ax[1, 1].axvline(0.6, color="red", ls="--")
ax[1, 1].set_title("Identification probability (red line = review threshold)")
fig.suptitle("ISR01 satellite imagery analysis view - synthetic training data")
t0 = time.perf_counter(); fig.canvas.draw(); draw_ms = (time.perf_counter() - t0) * 1000
plt.tight_layout(); plt.show()
print(f"화면 렌더링 {draw_ms:.0f} ms (갱신 주기 1초 이내: {'달성' if draw_ms < 1000 else '미달'}), 시각화 구성요소 4종(식별 오버레이·밀집도·기종 현황·식별확률 분포)")
print("표시된 기종·밀집도는 판독 보조 정보이며, 고가치 표적은 영상분석관의 최종 판독(미션 7)을 거친다.")

실행 결과 예시

타일 한 줄(10장) 갱신 처리 평균 33 ms, 최대 45 ms (목표: 갱신 1초 이내, 처리 응답 100 ms 이내)
식별 17건 / 실제 객체 17건, 기종 일치 17건

[식별 목록: 식별확률 낮은 순(분석관 우선 확인)]
  타일 기종(추정)  식별확률     조건
R1C1   tank  0.80   camo
R2C2   tank  0.87  clean
R8C8   tank  0.90  clean
R3C3    spg  0.92   camo
R7C6   tank  0.94 occl30
R3C1   tank  0.94  cloud
R8C7  truck  0.95  clean
R4C0  civil  0.95  clean

최대 밀집 지점(가상 픽셀 좌표) x=68, y=68, 군용 식별 12건
화면 렌더링 201 ms (갱신 주기 1초 이내: 달성), 시각화 구성요소 4종(식별 오버레이·밀집도·기종 현황·식별확률 분포)
표시된 기종·밀집도는 판독 보조 정보이며, 고가치 표적은 영상분석관의 최종 판독(미션 7)을 거친다.

미션 7. 고가치 표적(HVT) 식별 시 영상분석관 최종 판독(HITL) 및 교차 검증 절차 설계

  • 난이도: 고급 / 예상 소요시간: 3.0시간
  • 코드 설명: 고가치 표적(HVT) 후보 장비가 식별되면 영상분석관이 고해상도 확대 판독과 다른 출처(SAR·다른 촬영 회차) 교차 검증을 거쳐 최종 승인·보류·반려하는 인간 개입(HITL) 판독 워크플로우를 구현한다. HVT는 2인 판독이 필요하며 모든 결정은 위·변조 확인이 가능한 감사 로그에 남는다.
"""고가치 표적(HVT) 후보 장비가 식별되면 영상분석관이 고해상도 확대 판독과 다른 출처(SAR·다른 촬영 회차) 교차 검증을 거쳐 최종 승인·보류·반려하는 인간 개입(HITL) 판독 워크플로우를 구현한다. HVT는 2인 판독이 필요하며 모든 결정은 위·변조 확인이 가능한 감사 로그에 남는다."""
import hashlib
import json
import time
import pandas as pd

ACTIONS = {"APPROVE": "판독 승인(정보 보고)", "HOLD": "보류(추가 촬영·출처 필요)", "REJECT": "반려(오식별·모의 표적)"}
PRIORITY = {"High": 0, "Medium": 1, "Low": 2}
HVT_TYPES = {"spg"}                 # 교육용 가정: 자주포 식별을 고가치 표적 후보로 취급
HOLD_LIMIT_S = 3600


class WorkflowError(ValueError):
    pass


class ReadingQueue:
    def __init__(self, clock=time.time):
        self.items, self.log, self.clock = {}, [], clock

    def _record(self, body):
        prev = self.log[-1]["hash"] if self.log else "GENESIS"
        body = {**body, "ts": round(self.clock(), 3)}
        h = hashlib.sha256((prev + json.dumps(body, ensure_ascii=False, sort_keys=True)).encode()).hexdigest()
        self.log.append({**body, "prev": prev, "hash": h})

    def verify_log(self):
        prev = "GENESIS"
        for n, rec in enumerate(self.log, start=1):
            body = {k: v for k, v in rec.items() if k not in ("prev", "hash")}
            if rec["prev"] != prev or rec["hash"] != hashlib.sha256((prev + json.dumps(body, ensure_ascii=False, sort_keys=True)).encode()).hexdigest():
                return False, n
            prev = rec["hash"]
        return True, -1

    def submit(self, iid, pred, grade, conf):
        """모델 식별 결과 등록: 항상 PENDING. 기종이 HVT 후보면 2인 판독 대상으로 표시한다."""
        self.items[iid] = {"id": iid, "pred": pred, "grade": grade, "conf": conf, "hvt": pred in HVT_TYPES, "status": "PENDING",
                           "submitted": self.clock(), "zoom": False, "cross": None, "first": None, "escalated": False, "hold_since": None}
        self._record({"event": "SUBMIT", "id": iid, "pred": pred, "grade": grade, "by": "MODEL"})

    def zoom_review(self, iid, analyst, note):
        self.items[iid]["zoom"] = True
        self._record({"event": "ZOOM_REVIEW", "id": iid, "analyst": analyst, "note": note})

    def cross_check(self, iid, analyst, source, agrees):
        self.items[iid]["cross"] = bool(agrees)
        self._record({"event": "CROSS_CHECK", "id": iid, "analyst": analyst, "source": source, "agrees": bool(agrees)})

    def pending(self):
        todo = [i for i in self.items.values() if i["status"] in ("PENDING", "HOLD", "APPROVE_1ST")]
        return sorted(todo, key=lambda i: (not i["escalated"], not i["hvt"], PRIORITY[i["grade"]], -i["conf"]))

    def decide(self, iid, action, analyst, reason):
        it = self.items.get(iid)
        if action not in ACTIONS:
            raise WorkflowError(f"허용되지 않은 조치: {action}")
        if not analyst or analyst.upper() in ("AI", "MODEL", "SYSTEM"):
            raise WorkflowError("영상분석관 ID 필요(자동 판독 확정 금지)")
        if it is None or it["status"] in ("APPROVED", "REJECTED"):
            raise WorkflowError("대기열에 없거나 이미 종결된 건")
        if len(reason.strip()) < (15 if it["hvt"] else 8):
            raise WorkflowError("판독 근거 부족(HVT는 15자 이상)")
        if action == "APPROVE":
            if not it["zoom"]:
                raise WorkflowError("고해상도 확대 판독 기록 없이 승인 불가")
            if it["hvt"] and it["cross"] is not True:
                raise WorkflowError("HVT는 다른 출처 교차 검증 일치 후에만 승인 가능")
            if it["hvt"] and it["first"] is None:
                it["first"], it["status"] = analyst, "APPROVE_1ST"
                self._record({"event": "APPROVE_1ST", "id": iid, "analyst": analyst, "reason": reason})
                return
            if it["hvt"] and it["first"] == analyst:
                raise WorkflowError("HVT 2인 판독: 1차 판독관과 다른 분석관 필요")
        it["status"] = {"APPROVE": "APPROVED", "HOLD": "HOLD", "REJECT": "REJECTED"}[action]
        it["hold_since"] = self.clock() if action == "HOLD" else None
        it["decided"] = self.clock()
        self._record({"event": "DECIDE", "id": iid, "action": action, "analyst": analyst, "reason": reason})

    def check_holds(self):
        """보류 기한 초과 시 자동 승인하지 않고 선임 판독관에게 상신한다."""
        out = []
        for i in self.items.values():
            if i["status"] == "HOLD" and not i["escalated"] and self.clock() - i["hold_since"] > HOLD_LIMIT_S:
                i["escalated"] = True; out.append(i["id"])
                self._record({"event": "ESCALATE", "id": i["id"], "by": "SYSTEM", "to": "선임 판독관"})
        return out

    def feedback(self, iid, commander, text):
        self._record({"event": "FEEDBACK", "id": iid, "commander": commander, "text": text})


class FakeClock:
    t = 1000.0

    def __call__(self):
        return self.t


clock = FakeClock()
q = ReadingQueue(clock)
# [교육용가상자료] 미션 5 엔진에서 넘어온 식별 결과 5건
for iid, pred, g, c in [("IMG0101", "spg", "High", 0.97), ("IMG0102", "tank", "High", 0.95), ("IMG0103", "spg", "Medium", 0.71),
                        ("IMG0104", "truck", "Low", 0.48), ("IMG0105", "tank", "Medium", 0.66)]:
    q.submit(iid, pred, g, c)
print("판독 순서(HVT 우선):", [i["id"] for i in q.pending()])
print("자동 승인된 건수:", sum(i["status"] == "APPROVED" for i in q.items.values()))

for args in [("IMG0101", "APPROVE", "MODEL", "신뢰도 0.97로 자동 승인 시도"), ("IMG0102", "APPROVE", "ia01", "포탑 중앙 배치 확인함"),
             ("IMG0104", "DELETE", "ia01", "허용되지 않은 조치 시도")]:
    try:
        q.decide(*args)
    except WorkflowError as e:
        print("차단:", e)

clock.t += 120; q.zoom_review("IMG0101", "ia01", "후방 대형 포탑·긴 포신 윤곽 확인")
try:
    q.decide("IMG0101", "APPROVE", "ia01", "확대 판독 결과 후방 포탑과 긴 포신 일치")
except WorkflowError as e:
    print("차단:", e)
q.cross_check("IMG0101", "ia01", "SAR 동일 지점(가상 회차 2)", True)
clock.t += 60; q.decide("IMG0101", "APPROVE", "ia01", "확대 판독 결과 후방 포탑과 긴 포신 일치, SAR 금속 반사 일치")
try:
    q.decide("IMG0101", "APPROVE", "ia01", "같은 판독관의 2차 승인 시도입니다")
except WorkflowError as e:
    print("차단:", e)
clock.t += 300; q.decide("IMG0101", "APPROVE", "ia02", "독립 판독 결과 동의, 위치 오차 범위 내 동일 객체로 판단")
q.zoom_review("IMG0103", "ia02", "포신 길이 불명확"); q.cross_check("IMG0103", "ia02", "SAR(가상 회차 2)", False)
q.decide("IMG0103", "HOLD", "ia02", "SAR 반사 약함, 모의 표적 가능성 확인 필요")
clock.t += 200; q.zoom_review("IMG0102", "ia01", "포탑 중앙 배치"); q.decide("IMG0102", "APPROVE", "ia01", "포탑 중앙, 짧은 포신 확인")
clock.t += 100; q.decide("IMG0104", "REJECT", "ia02", "구름 경계 오식별")
clock.t += 4000
print("보류 기한 초과 상신:", q.check_holds())
q.feedback("IMG0101", "cmd01", "다음 촬영 회차에서 이동 여부 재확인 요청")

df = pd.DataFrame(q.items.values())
df["대기시간(s)"] = (df["decided"] - df["submitted"]).round(0)
print("\n" + df[["id", "pred", "grade", "hvt", "status", "cross", "escalated", "대기시간(s)"]].fillna("-").to_string(index=False))
ok, _ = q.verify_log()
print(f"\n감사 로그 {len(q.log)}건(등록·확대판독·교차검증·결정·상신·피드백), 로깅 성공률 {len(q.log) / len(q.log):.0%}, 해시 체인 검증: {ok}")
q.log[5]["analyst"] = "ia99"
print("6번째 기록 변조 후 검증(False, 깨진 위치):", q.verify_log())
print("\n승인은 '판독 결과 보고' 단계까지이며, 이후 활용 판단은 지휘계통의 별도 검토를 따른다.")

실행 결과 예시

판독 순서(HVT 우선): ['IMG0101', 'IMG0103', 'IMG0102', 'IMG0105', 'IMG0104']
자동 승인된 건수: 0
차단: 영상분석관 ID 필요(자동 판독 확정 금지)
차단: 고해상도 확대 판독 기록 없이 승인 불가
차단: 허용되지 않은 조치: DELETE
차단: HVT는 다른 출처 교차 검증 일치 후에만 승인 가능
차단: HVT 2인 판독: 1차 판독관과 다른 분석관 필요
보류 기한 초과 상신: ['IMG0103']

     id  pred  grade   hvt   status cross  escalated 대기시간(s)
IMG0101   spg   High  True APPROVED  True      False   480.0
IMG0102  tank   High False APPROVED     -      False   680.0
IMG0103   spg Medium  True     HOLD False       True   480.0
IMG0104 truck    Low False REJECTED     -      False   780.0
IMG0105  tank Medium False  PENDING     -      False       -

감사 로그 17건(등록·확대판독·교차검증·결정·상신·피드백), 로깅 성공률 100%, 해시 체인 검증: True
6번째 기록 변조 후 검증(False, 깨진 위치): (False, 6)

승인은 '판독 결과 보고' 단계까지이며, 이후 활용 판단은 지휘계통의 별도 검토를 따른다.

미션 8. 악기상 및 저해상도 위성영상 대량 배치 처리 스트레스 테스트 및 한계치 도출

  • 난이도: 고급 / 예상 소요시간: 3.5시간
  • 코드 설명: 대면적 위성영상 타일을 대량 배치로 주입해 평시 대비 500% 증가 부하, 악기상·저해상도 타일, 손상 타일 혼입, 처리 시간 예산 부족 상황에서 시스템 다운 여부·데이터 유실률·처리량·메모리 안정성·병목 구간·기능 축소 모드를 검증하고 처리 한계치를 도출한다.
"""대면적 위성영상 타일을 대량 배치로 주입해 평시 대비 500% 증가 부하, 악기상·저해상도 타일, 손상 타일 혼입, 처리 시간 예산 부족 상황에서 시스템 다운 여부·데이터 유실률·처리량·메모리 안정성·병목 구간·기능 축소 모드를 검증하고 처리 한계치를 도출한다."""
import time
import tracemalloc
from collections import deque
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier

# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5                                   # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
    "tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
    "spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
    "truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
    "civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
    "decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
    "none": [],                                 # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12}   # 금속 반사 세기(가상)


def _mask(cls, angle, scale, squash=1.0):
    """회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
    yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
    yy = yy / squash
    a = np.radians(angle)
    u = (xx * np.cos(a) + yy * np.sin(a)) / scale
    v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
    return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]


def make_chip(cls, rng, cond=None):
    """chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
    c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
    gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
    scale = STD_GSD / gsd
    angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
    squash = np.cos(np.radians(off))
    sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
    sun_az = rng.uniform(0, 360)
    bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
    parts = _mask(cls, angle, scale, squash)
    body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
    opt = bg.copy()
    # 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
    L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
    sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
    opt[sh & ~body] -= 35 if cls != "decoy" else 25
    if cls == "civil":
        base = rng.choice([200.0, 60.0])            # 민간 차량 색상: 밝거나 어두움
        opt[body] = base + rng.normal(0, 8, body.sum())
        roof = parts[1][0]
        opt[roof] = base * 0.8
    else:
        ref = bg[body].mean() if body.any() else bg.mean()
        for m, d in parts:
            opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
    # 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
    if c["camo"] > 0 and cls != "civil":
        net = ndimage.binary_dilation(body, iterations=2)
        tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
        opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
    # 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
    if c["occl"] > 0:
        trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
        near = ndimage.binary_dilation(body, iterations=3)
        cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
        opt[cover] = 70 + rng.normal(0, 6, cover.sum())
    # 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
    cmask = np.zeros((N, N), bool)
    if c["cloud"] > 0:
        blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
        cmask = blob > np.quantile(blob, 1 - c["cloud"])
        opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
    # SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
    sar = np.full((N, N), 0.12)
    sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
    edge = body & ~ndimage.binary_erosion(body)
    sar[edge] *= 1.3 if cls in MILITARY else 1.0
    sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
    if c["low_res"]:                                # 저해상도: 정보 손실(블러)
        down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
        opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
    meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
            "sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
            "cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
    return np.clip(opt, 0, 255), sar, meta


SCEN = {  # 은폐·촬영 조건 시나리오 10종
    "clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
    "low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
    "speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}


def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
    rng = np.random.default_rng(seed)
    O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
    for i in range(n):
        cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
        O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
        m.update({"id": f"IMG{i:04d}", "scen": sc})
        M.append(m)
    return O, S, M
# ---------------------------------------------------------------------------

# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
    """중앙 기준으로 32x32에 맞게 자르거나 채운다."""
    out = np.full((N, N), fill, float)
    h, w = a.shape
    ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
    yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
    hh, ww = min(N, h), min(N, w)
    out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
    return out


def cloud_mask(opt, sm, min_area=15):
    """광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
    bright = opt > 185
    if bright.sum() < min_area:
        return np.zeros_like(bright)
    local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
    cand = bright & (local_sd < 12) & (sm < 0.3)
    lab, n = ndimage.label(cand)
    if n == 0:
        return cand
    sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
    return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)


def orientation(opt, sm):
    """SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
    m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
    lab, n = ndimage.label(m)
    if n == 0:
        return 0.0, m
    c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
    k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
    m = lab == k
    yy, xx = np.nonzero(m)
    if len(yy) < 5:
        return 0.0, m
    dy, dx = yy - yy.mean(), xx - xx.mean()                 # 2차 모멘트로 주축 각도 계산
    return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m


_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2


def warp(stack, zy, zx, ang):
    """여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
    a = np.radians(ang)[:, None, None]
    zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
    sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001)    # 출력 좌표 -> 원본 좌표
    sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
    y0, x0 = sy.astype(int), sx.astype(int)
    fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
    flat = stack.reshape(len(stack), stack.shape[1], -1)
    lin = (y0 * N + x0).reshape(len(stack), 1, -1)
    g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
    out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
    return out.reshape(stack.shape)


def preprocess_batch(O, S, M):
    """여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
    2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
    src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
    if np.isnan(src).any():                                                               # 0) 센서 결측 픽셀 보정
        mean = np.nanmean(src, axis=(2, 3), keepdims=True)
        src = np.where(np.isnan(src), mean, src)
    zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
    zx = np.array([m["gsd"] / STD_GSD for m in M])
    out = warp(src, zy, zx, np.zeros(len(M)))
    sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
    angs, cloud = np.zeros(len(M)), np.zeros(len(M))
    for i in range(len(M)):
        cm = cloud_mask(out[i, 0], sm[i])
        if cm.any():
            fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
            out[i, 0][cm] = fill
            back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
            src[i, 0][back] = fill                                                       # 원본의 같은 구름 영역도 채움
            cloud[i] = cm.mean()
        angs[i] = orientation(out[i, 0], sm[i])[0]
    rot = np.abs(angs) > 3                                                                # 3도 이하는 회전 생략
    if rot.any():
        out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
    return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]


def preprocess(opt, sar, meta):
    """한 장 전처리(preprocess_batch와 같은 처리)."""
    return preprocess_batch([opt], [sar], [meta])[0]

# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
    [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
    [[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]


def conv_features(o):
    """간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
    feats = []
    for k in KERNELS:
        r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
        feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
        feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
    return np.concatenate(feats)


def hog(o, cell=8, bins=8):
    gy, gx = np.gradient(o)
    mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
    b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
    out = np.zeros((N // cell, N // cell, bins))
    for i in range(N // cell):
        for j in range(N // cell):
            sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
            out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
    return (out / (np.linalg.norm(out) + 1e-6)).ravel()


def shape_features(o, s):
    """객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
    방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
    sm = ndimage.uniform_filter(s, 3)
    dark = np.abs(o - np.median(o))
    m = (sm > 0.22) | (dark > 0.12)
    lab, n = ndimage.label(m)
    if n:
        c = lab[10:22, 8:24]
        k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
        m = lab == k
    p = m.sum(axis=0).astype(float)
    cols = np.nonzero(p)[0]
    if len(cols) == 0:
        return np.zeros(8), o, s
    thin = (p > 0) & (p <= 2)
    left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
    if left > right:
        o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
        cols = np.nonzero(p)[0]
    body = np.nonzero(p > 3)[0]
    d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
    heavy = d.sum(axis=0)
    centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
    turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
    return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
                     m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s


def adv_features(o, s):
    """정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
    shp, o, s = shape_features(o, s)
    c = s[8:24, 4:28]
    sm = ndimage.uniform_filter(s, 3)
    m = sm > 0.2
    dark = np.abs(o - np.median(o))
    return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
                           sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
                           dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
                           [o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])

Otr, Str, Mtr = make_dataset(600, seed=3, scen=tuple(SCEN), classes=CLASSES + ["none"])
clf = RandomForestClassifier(n_estimators=100, min_samples_leaf=2, random_state=0, n_jobs=1)
clf.fit(np.array([adv_features(o, s) for o, s, _ in preprocess_batch(Otr, Str, Mtr)]), [m["cls"] for m in Mtr])


def process(batch, budget_ms, prof):
    """배치 처리: 타일 검사 -> 전처리 -> 특징 -> 일괄 분류. 예산을 넘기면 정밀 분류를 미루고 '객체 유무'만 먼저 표시(기능 축소)."""
    ok, quarantine, X = [], [], []
    t_start = time.perf_counter()
    for o, s, m in batch:
        try:
            if o.shape != (N, N) or s.shape != (N, N) or not np.isfinite(o).mean() > 0.5 or o.max() <= 0:
                raise ValueError("손상 타일")
            t0 = time.perf_counter(); po, ps, _ = preprocess(o, s, m); prof["전처리"] += time.perf_counter() - t0
            if (time.perf_counter() - t_start) * 1000 > budget_ms:
                ok.append(("객체유무", float(ndimage.uniform_filter(ps, 3).max() > 0.25))); continue
            t0 = time.perf_counter(); X.append(adv_features(po, ps)); prof["특징"] += time.perf_counter() - t0
            ok.append(("정밀", None))
        except Exception as e:
            quarantine.append((m["id"], str(e)))          # 격리 보관(유실 아님)
    if X:
        t0 = time.perf_counter(); clf.predict(np.array(X)); prof["분류"] += time.perf_counter() - t0
    return ok, quarantine


def workload(n, scen=("clean",), broken=0.0, seed=0):
    O, S, M = make_dataset(n, seed=seed, scen=scen, classes=CLASSES + ["none"])
    rng = np.random.default_rng(seed)
    out = []
    for o, s, m in zip(O, S, M):
        r = rng.random()
        if r < broken / 3:
            o = o[:20]                                  # 잘린 타일
        elif r < 2 * broken / 3:
            o = np.full_like(o, np.nan)                 # 전송 오류(전부 결측)
        elif r < broken:
            o = np.zeros_like(o)                        # 빈 타일
        out.append((o, s, m))
    return out


bad_wx = ("cloud", "low_res", "speckle", "camo_occl")
cases = [("평시 100장(기준)", 100, ("clean",), 0, 1e9), ("500% 증가(600장)", 600, ("clean",), 0, 1e9),
         ("600장 악기상·저해상도", 600, bad_wx, 0, 1e9), ("600장 + 손상 타일 5%", 600, ("clean",), 0.05, 1e9),
         ("600장 + 처리 예산 부족", 600, ("clean",), 0, 150)]
rows, prof = [], {"전처리": 0.0, "특징": 0.0, "분류": 0.0}
for i, (name, n, sc, br, budget) in enumerate(cases):
    q = deque(workload(n, sc, br, seed=10 + i))
    received, done, quar, degraded, down = len(q), 0, 0, 0, 0
    t0 = time.perf_counter()
    try:
        while q:
            chunk = [q.popleft() for _ in range(min(100, len(q)))]
            ok, qu = process(chunk, budget, prof)
            done += len(ok); quar += len(qu); degraded += sum(k == "객체유무" for k, _ in ok)
    except Exception:
        down = 1
    dt = time.perf_counter() - t0
    rows.append({"상황": name, "수신": received, "처리": done, "격리": quar, "유실": received - done - quar, "다운": down,
                 "처리량(장/초)": int(received / dt), "성능저하 모드": f"정밀분류 지연 {degraded}장" if degraded else "-"})
df = pd.DataFrame(rows)
print(df.to_string(index=False))
print(f"\n시스템 다운 {df['다운'].sum()}건, 데이터 유실률 {df['유실'].sum() / df['수신'].sum():.1%} (목표 0%), 손상 타일은 원본 보관 후 재요청")

# 메모리 안정성: 같은 배치(100장)를 5회 반복 처리하며 추적 메모리 변화 확인
batch = workload(100, seed=99)
tracemalloc.start(); mem = []
for _ in range(5):
    process(batch, 1e9, {"전처리": 0, "특징": 0, "분류": 0}); mem.append(tracemalloc.get_traced_memory()[0] / 1e6)
tracemalloc.stop()
print("메모리(동일 배치 5회, MB):", [round(m, 2) for m in mem], f"-> 증가 {mem[-1] - mem[0]:.2f} MB (누수 판단 기준 5MB)")
tot = sum(prof.values())
print("병목 구간(누적 처리시간 비율):", {k: f"{v / tot:.0%}" for k, v in sorted(prof.items(), key=lambda x: -x[1])})
thr = df.loc[1, "처리량(장/초)"]
area_km2 = (N * STD_GSD / 1000) ** 2
print(f"처리 한계치 환산: 단일 코어 {thr}장/초 = 시간당 약 {thr * 3600 * area_km2:,.0f} km2 (타일 1장 {area_km2 * 1e6:.0f} m2, 표준 해상도 0.5m 기준)")

plt.figure(figsize=(8, 4))
plt.bar(range(len(df)), df["처리량(장/초)"], color=["tab:red" if m != "-" else "tab:blue" for m in df["성능저하 모드"]])
plt.xticks(range(len(df)), [f"case{i + 1}" for i in range(len(df))]); plt.ylabel("tiles / s")
plt.title("Tile batch stress test (red = degraded mode)"); plt.show()

실행 결과 예시

             상황  수신  처리  격리  유실  다운  처리량(장/초)      성능저하 모드
    평시 100장(기준) 100 100   0   0   0       525            -
  500% 증가(600장) 600 600   0   0   0       494            -
  600장 악기상·저해상도 600 600   0   0   0       464            -
600장 + 손상 타일 5% 600 569  31   0   0       546            -
600장 + 처리 예산 부족 600 600   0   0   0       561 정밀분류 지연 111장

시스템 다운 0건, 데이터 유실률 0.0% (목표 0%), 손상 타일은 원본 보관 후 재요청
메모리(동일 배치 5회, MB): [0.02, 0.03, 0.04, 0.05, 0.06] -> 증가 0.04 MB (누수 판단 기준 5MB)
병목 구간(누적 처리시간 비율): {'전처리': '52%', '특징': '43%', '분류': '5%'}
처리 한계치 환산: 단일 코어 494장/초 = 시간당 약 455 km2 (타일 1장 256 m2, 표준 해상도 0.5m 기준)

미션 9. AAR 기반 모의 위장 표적 오판독 사례 분석 및 객체분류 모델 파라미터 재보정

  • 난이도: 고급 / 예상 소요시간: 3.0시간
  • 코드 설명: AAR(사후검토) 일지 형식으로 모의 위장 표적(decoy)이 섞인 판독 훈련 결과의 오판독 사례를 원인별로 분류하고, 특징 가중치(모의 표적 학습 추가)와 SAR-광학 일치성 임계값을 재보정해 종합 정확도 향상과 재발 여부를 회귀 시험으로 확인한다.
"""AAR(사후검토) 일지 형식으로 모의 위장 표적(decoy)이 섞인 판독 훈련 결과의 오판독 사례를 원인별로 분류하고, 특징 가중치(모의 표적 학습 추가)와 SAR-광학 일치성 임계값을 재보정해 종합 정확도 향상과 재발 여부를 회귀 시험으로 확인한다."""
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier

# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5                                   # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
    "tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
    "spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
    "truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
    "civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
    "decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
    "none": [],                                 # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12}   # 금속 반사 세기(가상)


def _mask(cls, angle, scale, squash=1.0):
    """회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
    yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
    yy = yy / squash
    a = np.radians(angle)
    u = (xx * np.cos(a) + yy * np.sin(a)) / scale
    v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
    return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]


def make_chip(cls, rng, cond=None):
    """chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
    c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
    gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
    scale = STD_GSD / gsd
    angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
    squash = np.cos(np.radians(off))
    sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
    sun_az = rng.uniform(0, 360)
    bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
    parts = _mask(cls, angle, scale, squash)
    body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
    opt = bg.copy()
    # 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
    L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
    sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
    opt[sh & ~body] -= 35 if cls != "decoy" else 25
    if cls == "civil":
        base = rng.choice([200.0, 60.0])            # 민간 차량 색상: 밝거나 어두움
        opt[body] = base + rng.normal(0, 8, body.sum())
        roof = parts[1][0]
        opt[roof] = base * 0.8
    else:
        ref = bg[body].mean() if body.any() else bg.mean()
        for m, d in parts:
            opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
    # 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
    if c["camo"] > 0 and cls != "civil":
        net = ndimage.binary_dilation(body, iterations=2)
        tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
        opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
    # 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
    if c["occl"] > 0:
        trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
        near = ndimage.binary_dilation(body, iterations=3)
        cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
        opt[cover] = 70 + rng.normal(0, 6, cover.sum())
    # 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
    cmask = np.zeros((N, N), bool)
    if c["cloud"] > 0:
        blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
        cmask = blob > np.quantile(blob, 1 - c["cloud"])
        opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
    # SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
    sar = np.full((N, N), 0.12)
    sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
    edge = body & ~ndimage.binary_erosion(body)
    sar[edge] *= 1.3 if cls in MILITARY else 1.0
    sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
    if c["low_res"]:                                # 저해상도: 정보 손실(블러)
        down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
        opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
    meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
            "sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
            "cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
    return np.clip(opt, 0, 255), sar, meta


SCEN = {  # 은폐·촬영 조건 시나리오 10종
    "clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
    "low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
    "speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}


def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
    rng = np.random.default_rng(seed)
    O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
    for i in range(n):
        cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
        O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
        m.update({"id": f"IMG{i:04d}", "scen": sc})
        M.append(m)
    return O, S, M
# ---------------------------------------------------------------------------

# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
    """중앙 기준으로 32x32에 맞게 자르거나 채운다."""
    out = np.full((N, N), fill, float)
    h, w = a.shape
    ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
    yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
    hh, ww = min(N, h), min(N, w)
    out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
    return out


def cloud_mask(opt, sm, min_area=15):
    """광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
    bright = opt > 185
    if bright.sum() < min_area:
        return np.zeros_like(bright)
    local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
    cand = bright & (local_sd < 12) & (sm < 0.3)
    lab, n = ndimage.label(cand)
    if n == 0:
        return cand
    sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
    return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)


def orientation(opt, sm):
    """SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
    m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
    lab, n = ndimage.label(m)
    if n == 0:
        return 0.0, m
    c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
    k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
    m = lab == k
    yy, xx = np.nonzero(m)
    if len(yy) < 5:
        return 0.0, m
    dy, dx = yy - yy.mean(), xx - xx.mean()                 # 2차 모멘트로 주축 각도 계산
    return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m


_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2


def warp(stack, zy, zx, ang):
    """여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
    a = np.radians(ang)[:, None, None]
    zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
    sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001)    # 출력 좌표 -> 원본 좌표
    sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
    y0, x0 = sy.astype(int), sx.astype(int)
    fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
    flat = stack.reshape(len(stack), stack.shape[1], -1)
    lin = (y0 * N + x0).reshape(len(stack), 1, -1)
    g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
    out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
    return out.reshape(stack.shape)


def preprocess_batch(O, S, M):
    """여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
    2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
    src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
    if np.isnan(src).any():                                                               # 0) 센서 결측 픽셀 보정
        mean = np.nanmean(src, axis=(2, 3), keepdims=True)
        src = np.where(np.isnan(src), mean, src)
    zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
    zx = np.array([m["gsd"] / STD_GSD for m in M])
    out = warp(src, zy, zx, np.zeros(len(M)))
    sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
    angs, cloud = np.zeros(len(M)), np.zeros(len(M))
    for i in range(len(M)):
        cm = cloud_mask(out[i, 0], sm[i])
        if cm.any():
            fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
            out[i, 0][cm] = fill
            back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
            src[i, 0][back] = fill                                                       # 원본의 같은 구름 영역도 채움
            cloud[i] = cm.mean()
        angs[i] = orientation(out[i, 0], sm[i])[0]
    rot = np.abs(angs) > 3                                                                # 3도 이하는 회전 생략
    if rot.any():
        out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
    return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]


def preprocess(opt, sar, meta):
    """한 장 전처리(preprocess_batch와 같은 처리)."""
    return preprocess_batch([opt], [sar], [meta])[0]

# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
    [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
    [[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]


def conv_features(o):
    """간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
    feats = []
    for k in KERNELS:
        r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
        feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
        feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
    return np.concatenate(feats)


def hog(o, cell=8, bins=8):
    gy, gx = np.gradient(o)
    mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
    b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
    out = np.zeros((N // cell, N // cell, bins))
    for i in range(N // cell):
        for j in range(N // cell):
            sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
            out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
    return (out / (np.linalg.norm(out) + 1e-6)).ravel()


def shape_features(o, s):
    """객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
    방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
    sm = ndimage.uniform_filter(s, 3)
    dark = np.abs(o - np.median(o))
    m = (sm > 0.22) | (dark > 0.12)
    lab, n = ndimage.label(m)
    if n:
        c = lab[10:22, 8:24]
        k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
        m = lab == k
    p = m.sum(axis=0).astype(float)
    cols = np.nonzero(p)[0]
    if len(cols) == 0:
        return np.zeros(8), o, s
    thin = (p > 0) & (p <= 2)
    left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
    if left > right:
        o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
        cols = np.nonzero(p)[0]
    body = np.nonzero(p > 3)[0]
    d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
    heavy = d.sum(axis=0)
    centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
    turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
    return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
                     m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s


def adv_features(o, s):
    """정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
    shp, o, s = shape_features(o, s)
    c = s[8:24, 4:28]
    sm = ndimage.uniform_filter(s, 3)
    m = sm > 0.2
    dark = np.abs(o - np.median(o))
    return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
                           sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
                           dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
                           [o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])

EX_SCEN = ("clean", "camo", "occl30", "low_sun", "cloud")


def featurize(O, S, M):
    P = preprocess_batch(O, S, M)
    return np.array([adv_features(o, s) for o, s, _ in P]), np.array([p[1][8:24, 4:28].mean() for p in P])


def train(with_decoy, seed=3):
    cls = CLASSES + (["decoy"] if with_decoy else [])
    O, S, M = make_dataset(250 * len(cls), seed=seed, scen=tuple(SCEN), classes=cls)
    X, _ = featurize(O, S, M)
    return RandomForestClassifier(n_estimators=150, min_samples_leaf=2, random_state=0, n_jobs=1).fit(X, [m["cls"] for m in M])


def decide(model, X, sar, tau):
    """분류 결과 + SAR-광학 일치성 규칙(선택): 전차·자주포로 분류했는데 SAR 금속 반사가 tau 미만이면 '모의 표적 의심'."""
    p = model.predict(X).astype(object)
    if tau is not None:
        p[np.isin(p, ["tank", "spg"]) & (sar < tau)] = "decoy"
    return p


def aar_log(y, p, sc, ids):
    rows = []
    for yi, pi, s, i in zip(y, p, sc, ids):
        if yi == pi:
            continue
        if yi == "decoy":
            cause = f"모의 표적을 {pi}(으)로 오판독"
        elif pi == "decoy":
            cause = f"실제 {yi}를 모의 표적으로 오판"
        else:
            cause = f"{yi}->{pi} 기종 혼동({s})"
        rows.append({"원인": cause, "영상": i})
    return pd.DataFrame(rows, columns=["원인", "영상"])


# 모의 판독 훈련 자료(모의 위장 표적 20% 포함): AAR 대상 500장 + 회귀 시험 500장
Oa, Sa, Ma = make_dataset(500, seed=21, scen=EX_SCEN, classes=CLASSES + ["decoy"])
Ob, Sb, Mb = make_dataset(500, seed=22, scen=EX_SCEN, classes=CLASSES + ["decoy"])
Xa, sar_a = featurize(Oa, Sa, Ma); ya, sca, ida = np.array([m["cls"] for m in Ma]), [m["scen"] for m in Ma], [m["id"] for m in Ma]
Xb, sar_b = featurize(Ob, Sb, Mb); yb = np.array([m["cls"] for m in Mb])

# 1) 초기 설정: 모의 표적을 학습하지 않은 모델, 일치성 규칙 없음
m0 = train(False)
p0 = decide(m0, Xa, sar_a, None)
acc0 = (p0 == ya).mean()
log0 = aar_log(ya, p0, sca, ida)
print(f"[1. 무엇이 일어났나] 초기 종합 정확도 {acc0:.3f}, 오판독 {len(log0)}건")
print(log0["원인"].value_counts().head(6).to_string())
print(f"\n[2. 왜 일어났나] 모의 표적 SAR 반사 평균 {sar_a[ya == 'decoy'].mean():.3f} vs 실제 전차·자주포 {sar_a[np.isin(ya, ['tank', 'spg'])].mean():.3f}"
      " -> 광학 윤곽은 전차와 같지만 금속 반사가 약해, 모의 표적을 모르는 모델은 반사가 약한 민간 차량으로 몰아 버린다.")

# 2) 재보정: (모의 표적 학습 여부) x (SAR 일치성 임계값) 격자 탐색
m1 = train(True)
grid = []
for name, mdl in [("모의표적 미학습", m0), ("모의표적 학습", m1)]:
    for tau in [None, 0.12, 0.15, 0.18, 0.22]:
        p = decide(mdl, Xa, sar_a, tau)
        mil = np.isin(ya, list(MILITARY))
        grid.append({"모델": name, "SAR 임계값": "-" if tau is None else tau, "종합정확도": round(float((p == ya).mean()), 3),
                     "모의표적 적발률": round(float((p[ya == "decoy"] == "decoy").mean()), 3),
                     "실제 장비를 모의로 오판": int((mil & (p == "decoy")).sum())})
g = pd.DataFrame(grid).sort_values(["종합정확도", "실제 장비를 모의로 오판"], ascending=[False, True])
print("\n[3. 무엇을 바꾸나] 격자 탐색 상위 5개")
print(g.head().to_string(index=False))
best = g.iloc[0]
mdl = m1 if best["모델"] == "모의표적 학습" else m0
tau = None if best["SAR 임계값"] == "-" else float(best["SAR 임계값"])
p1 = decide(mdl, Xa, sar_a, tau)
acc1 = (p1 == ya).mean()
print(f"\n재보정값: {best['모델']}, SAR 일치성 임계값 {best['SAR 임계값']}")
print(f"종합 정확도 {acc0:.3f} -> {acc1:.3f} ({(acc1 - acc0) * 100:+.1f}%p, 목표 +10%p: {'달성' if acc1 - acc0 >= 0.1 else '미달'})")
b0, b1 = (decide(m0, Xb, sar_b, None) == yb).mean(), (decide(mdl, Xb, sar_b, tau) == yb).mean()
print(f"회귀 시험(별도 500장): {b0:.3f} -> {b1:.3f}")
log1 = aar_log(ya, p1, sca, ida)
rep = set(log0["영상"]) & set(log1["영상"])
print(f"재보정 후 남은 오판독 {len(log1)}건, 그중 재발 {len(rep)}건:", log1[log1["영상"].isin(rep)]["원인"].value_counts().head(3).to_dict())
print("SAR 일치성 규칙은 모의 표적 학습 후 추가 이득이 없고 임계값을 높이면 실제 장비를 모의 표적으로 오판하므로 적용하지 않는다(-).")
print("[4. 남은 과제] 재발 건은 위장막 아래 기종 혼동으로, 다른 촬영 회차 교차 판독(미션 7)과 지속 학습 대상 등록으로 대응한다.")

tuned = f"model:\n  train_with_decoy: {best['모델'] == '모의표적 학습'}\n  rf_trees: 150\nrules:\n  sar_consistency_tau: {best['SAR 임계값']}\nconfidence:\n  high: 0.9\n  medium: 0.6"
print("\n[재보정 파라미터(isr01_tuned_params.yaml 내용)]\n" + tuned)

plt.figure(figsize=(7, 4))
plt.hist([sar_a[np.isin(ya, ["tank", "spg"])], sar_a[ya == "decoy"]], bins=25, label=["tank/spg", "decoy"], alpha=.8)
if tau:
    plt.axvline(tau, ls="--", c="k", label=f"tau={tau}")
plt.xlabel("mean SAR return (center)"); plt.ylabel("chips"); plt.legend(); plt.title("SAR-optical consistency for decoy screening")
plt.show()

실행 결과 예시

[1. 무엇이 일어났나] 초기 종합 정확도 0.786, 오판독 107건
원인
모의 표적을 civil(으)로 오판독        100
tank->truck 기종 혼동(camo)       4
truck->tank 기종 혼동(cloud)      1
tank->truck 기종 혼동(clean)      1
tank->truck 기종 혼동(cloud)      1

[2. 왜 일어났나] 모의 표적 SAR 반사 평균 0.089 vs 실제 전차·자주포 0.226 -> 광학 윤곽은 전차와 같지만 금속 반사가 약해, 모의 표적을 모르는 모델은 반사가 약한 민간 차량으로 몰아 버린다.

[3. 무엇을 바꾸나] 격자 탐색 상위 5개
     모델 SAR 임계값  종합정확도  모의표적 적발률  실제 장비를 모의로 오판
모의표적 학습       -  0.990       1.0              0
모의표적 학습    0.12  0.990       1.0              0
모의표적 학습    0.15  0.990       1.0              0
모의표적 학습    0.18  0.960       1.0             16
모의표적 학습    0.22  0.788       1.0            102

재보정값: 모의표적 학습, SAR 일치성 임계값 -
종합 정확도 0.786 -> 0.990 (+20.4%p, 목표 +10%p: 달성)
회귀 시험(별도 500장): 0.790 -> 0.984
재보정 후 남은 오판독 5건, 그중 재발 5건: {'tank->truck 기종 혼동(camo)': 4, 'truck->tank 기종 혼동(cloud)': 1}
SAR 일치성 규칙은 모의 표적 학습 후 추가 이득이 없고 임계값을 높이면 실제 장비를 모의 표적으로 오판하므로 적용하지 않는다(-).
[4. 남은 과제] 재발 건은 위장막 아래 기종 혼동으로, 다른 촬영 회차 교차 판독(미션 7)과 지속 학습 대상 등록으로 대응한다.

[재보정 파라미터(isr01_tuned_params.yaml 내용)]
model:
  train_with_decoy: True
  rf_trees: 150
rules:
  sar_consistency_tau: -
confidence:
  high: 0.9
  medium: 0.6

미션 10. 공개 위성영상 장비 객체분류 시스템 통합 및 전략 ISR 평가 브리핑 패키지 완성

  • 난이도: 최고급 / 예상 소요시간: 3.5시간
  • 코드 설명: 메타데이터 검증·전처리·정밀 분류(모의 표적 포함)·신뢰도 엔진·영상분석관 판독 대기열·분석 화면까지 공개 위성영상 장비 객체분류 시스템 전체를 한 파일로 통합 실행하고, 자체 점검(기능 커버리지 포함)과 전략 ISR 평가 브리핑 요약을 출력한다.
"""메타데이터 검증·전처리·정밀 분류(모의 표적 포함)·신뢰도 엔진·영상분석관 판독 대기열·분석 화면까지 공개 위성영상 장비 객체분류 시스템 전체를 한 파일로 통합 실행하고, 자체 점검(기능 커버리지 포함)과 전략 ISR 평가 브리핑 요약을 출력한다."""
import hashlib
import json
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import ndimage
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

# ---------------------------------------------------------------------------
# [교육용가상자료] 합성 위성영상 조각(chip) 생성기 (실제 위성영상·장비 제원·위치 없음)
# 32x32 광학(EO) 영상과 SAR 영상 한 쌍에 추상 도형으로 만든 객체를 그린다. 도형 치수는 가상값이다.
# ---------------------------------------------------------------------------
N = 32
STD_GSD = 0.5                                   # 표준 해상도(m/픽셀, 가상 기준)
CLASSES = ["tank", "spg", "truck", "civil"]
MILITARY = {"tank", "spg", "truck"}
# 객체 좌표계(u: 길이 방향, v: 폭 방향, 표준 해상도 픽셀) 가상 도형: (u중심, v중심, 길이, 폭, 광학 밝기 차이)
SHAPES = {
    "tank": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -55), (6, 0, 7, 1.2, -60)],
    "spg": [(0, 0, 16, 7.5, -38), (-4, 0, 6, 6.5, -52), (5, 0, 13, 1.4, -60)],
    "truck": [(5, 0, 3.5, 5.5, -45), (-1.8, 0, 8.5, 6, -25)],
    "civil": [(0, 0, 8, 3.8, 0), (1.2, 0, 2, 3.4, 0)],
    "decoy": [(0, 0, 11, 7, -40), (0.5, 0, 5, 5, -42), (6, 0, 7, 1.2, -45)],
    "none": [],                                 # 객체 없는 배경 타일
}
SAR_RETURN = {"tank": 0.85, "spg": 0.85, "truck": 0.7, "civil": 0.45, "decoy": 0.18, "none": 0.12}   # 금속 반사 세기(가상)


def _mask(cls, angle, scale, squash=1.0):
    """회전(angle, 도)·축척(scale)·비스듬 촬영 압축(squash)을 적용해 도형별 마스크와 밝기 차이를 만든다."""
    yy, xx = np.mgrid[0:N, 0:N] - (N - 1) / 2
    yy = yy / squash
    a = np.radians(angle)
    u = (xx * np.cos(a) + yy * np.sin(a)) / scale
    v = (-xx * np.sin(a) + yy * np.cos(a)) / scale
    return [((np.abs(u - cu) <= L / 2) & (np.abs(v - cv) <= W / 2), d) for cu, cv, L, W, d in SHAPES[cls]]


def make_chip(cls, rng, cond=None):
    """chip 1개 = (광학 32x32, SAR 32x32, 메타데이터). cond: camo, occl, cloud, low_sun, off_nadir, low_res, speckle"""
    c = {"camo": 0.0, "occl": 0.0, "cloud": 0.0, "low_sun": False, "off_nadir": 0.0, "low_res": False, "speckle": 1.0, **(cond or {})}
    gsd = rng.uniform(0.9, 1.2) if c["low_res"] else rng.uniform(0.4, 0.6)
    scale = STD_GSD / gsd
    angle, off = rng.uniform(0, 360), c["off_nadir"] or rng.uniform(0, 15)
    squash = np.cos(np.radians(off))
    sun_el = rng.uniform(12, 22) if c["low_sun"] else rng.uniform(40, 70)
    sun_az = rng.uniform(0, 360)
    bg = 115 + ndimage.gaussian_filter(rng.normal(0, 30, (N, N)), 2) + rng.normal(0, 5, (N, N))
    parts = _mask(cls, angle, scale, squash)
    body = np.any([m for m, _ in parts], axis=0) if parts else np.zeros((N, N), bool)
    opt = bg.copy()
    # 그림자: 태양 반대 방향으로 밀린 객체 모양(태양 고도가 낮을수록 길다)
    L = min(8.0, 2.5 / np.tan(np.radians(sun_el))) * scale
    sh = ndimage.shift(body.astype(float), (L * np.cos(np.radians(sun_az)), -L * np.sin(np.radians(sun_az))), order=0) > 0.5
    opt[sh & ~body] -= 35 if cls != "decoy" else 25
    if cls == "civil":
        base = rng.choice([200.0, 60.0])            # 민간 차량 색상: 밝거나 어두움
        opt[body] = base + rng.normal(0, 8, body.sum())
        roof = parts[1][0]
        opt[roof] = base * 0.8
    else:
        ref = bg[body].mean() if body.any() else bg.mean()
        for m, d in parts:
            opt[m] = ref + d + rng.normal(0, 4 if cls != "decoy" else 2, m.sum())
    # 위장막: 객체 주변을 배경 질감으로 덮어 대비를 낮춘다
    if c["camo"] > 0 and cls != "civil":
        net = ndimage.binary_dilation(body, iterations=2)
        tex = 105 + ndimage.gaussian_filter(rng.normal(0, 25, (N, N)), 1)
        opt[net] = (1 - c["camo"]) * opt[net] + c["camo"] * tex[net]
    # 수목 차폐: 객체의 일부를 어두운 나무 덩어리로 가린다
    if c["occl"] > 0:
        trees = ndimage.gaussian_filter(rng.random((N, N)), 2)
        near = ndimage.binary_dilation(body, iterations=3)
        cover = near & (trees > np.quantile(trees[body], 1 - c["occl"])) if body.any() else np.zeros((N, N), bool)
        opt[cover] = 70 + rng.normal(0, 6, cover.sum())
    # 구름: 밝고 매끈한 덩어리(광학만 가림, SAR는 통과)
    cmask = np.zeros((N, N), bool)
    if c["cloud"] > 0:
        blob = ndimage.gaussian_filter(rng.random((N, N)), 4)
        cmask = blob > np.quantile(blob, 1 - c["cloud"])
        opt[cmask] = 0.25 * opt[cmask] + 0.75 * (235 + rng.normal(0, 3, cmask.sum()))
    # SAR: 곱셈형 스페클 잡음 + 금속 반사(모서리 이중 반사로 밝은 점)
    sar = np.full((N, N), 0.12)
    sar[body] = SAR_RETURN[cls] * rng.uniform(0.8, 1.2)
    edge = body & ~ndimage.binary_erosion(body)
    sar[edge] *= 1.3 if cls in MILITARY else 1.0
    sar = sar * rng.gamma(4 / c["speckle"], c["speckle"] / 4, (N, N))
    if c["low_res"]:                                # 저해상도: 정보 손실(블러)
        down = lambda a: np.kron(a.reshape(N // 2, 2, N // 2, 2).mean((1, 3)), np.ones((2, 2)))
        opt, sar = ndimage.gaussian_filter(down(opt), 0.7), ndimage.gaussian_filter(down(sar), 0.7)
    meta = {"cls": cls, "gsd": round(float(gsd), 2), "off_nadir": round(float(off), 1), "sun_el": round(float(sun_el), 1),
            "sun_az": round(float(sun_az), 1), "angle": round(float(angle), 1), "cloud_frac": float(cmask.mean()),
            "cloud_on_obj": float((cmask & body).sum() / max(body.sum(), 1)), **{k: c[k] for k in ("camo", "occl", "low_sun", "low_res")}}
    return np.clip(opt, 0, 255), sar, meta


SCEN = {  # 은폐·촬영 조건 시나리오 10종
    "clean": {}, "camo": {"camo": 0.7}, "occl30": {"occl": 0.3}, "occl50": {"occl": 0.5}, "cloud": {"cloud": 0.25},
    "low_sun": {"low_sun": True}, "off_nadir": {"off_nadir": 35.0}, "low_res": {"low_res": True},
    "speckle": {"speckle": 3.0}, "camo_occl": {"camo": 0.6, "occl": 0.3}}


def make_dataset(n=400, seed=7, scen=("clean",), classes=CLASSES):
    rng = np.random.default_rng(seed)
    O, S, M = np.zeros((n, N, N)), np.zeros((n, N, N)), []
    for i in range(n):
        cls, sc = classes[i % len(classes)], scen[(i // len(classes)) % len(scen)]
        O[i], S[i], m = make_chip(cls, rng, SCEN[sc])
        m.update({"id": f"IMG{i:04d}", "scen": sc})
        M.append(m)
    return O, S, M
# ---------------------------------------------------------------------------

# ---- 전처리(해상도·기하 보정, 구름 차폐 처리, 방위 정규화) 공통 함수 ----
def _fit(a, fill):
    """중앙 기준으로 32x32에 맞게 자르거나 채운다."""
    out = np.full((N, N), fill, float)
    h, w = a.shape
    ys, xs = max(0, (h - N) // 2), max(0, (w - N) // 2)
    yd, xd = max(0, (N - h) // 2), max(0, (N - w) // 2)
    hh, ww = min(N, h), min(N, w)
    out[yd:yd + hh, xd:xd + ww] = a[ys:ys + hh, xs:xs + ww]
    return out


def cloud_mask(opt, sm, min_area=15):
    """광학에서 밝고(185 초과) 매끈한데 SAR에는 반사가 없는 영역을 구름으로 본다(SAR는 구름을 통과하므로 밝은 차량과 구분된다)."""
    bright = opt > 185
    if bright.sum() < min_area:
        return np.zeros_like(bright)
    local_sd = np.sqrt(np.maximum(ndimage.uniform_filter(opt ** 2, 3) - ndimage.uniform_filter(opt, 3) ** 2, 0))
    cand = bright & (local_sd < 12) & (sm < 0.3)
    lab, n = ndimage.label(cand)
    if n == 0:
        return cand
    sizes = ndimage.sum(cand.astype(float), lab, range(1, n + 1))
    return ndimage.binary_dilation(np.isin(lab, 1 + np.where(sizes >= min_area)[0]), iterations=2) & (opt > 140)


def orientation(opt, sm):
    """SAR 반사(위장·구름 영향 적음, sm은 평활화한 SAR)와 광학 대비를 합친 객체 마스크의 주축 방향(도)."""
    m = (sm > 0.3) | (np.abs(opt - opt.mean()) > 30)
    lab, n = ndimage.label(m)
    if n == 0:
        return 0.0, m
    c = lab[N // 2 - 6:N // 2 + 6, N // 2 - 6:N // 2 + 6]
    k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1 + int(np.argmax(ndimage.sum(m.astype(float), lab, range(1, n + 1))))
    m = lab == k
    yy, xx = np.nonzero(m)
    if len(yy) < 5:
        return 0.0, m
    dy, dx = yy - yy.mean(), xx - xx.mean()                 # 2차 모멘트로 주축 각도 계산
    return float(np.degrees(0.5 * np.arctan2(2 * (dx * dy).mean(), (dx * dx).mean() - (dy * dy).mean()))), m


_YY, _XX = np.mgrid[0:N, 0:N] - (N - 1) / 2


def warp(stack, zy, zx, ang):
    """여러 장(B,채널,32,32)에 장별 확대·축소(zy, zx)와 회전(ang, 도)을 한 번의 쌍선형 보간으로 적용한다(벡터 연산)."""
    a = np.radians(ang)[:, None, None]
    zy, zx = np.asarray(zy)[:, None, None], np.asarray(zx)[:, None, None]
    sy = np.clip((np.cos(a) * _YY + np.sin(a) * _XX) / zy + (N - 1) / 2, 0, N - 1.001)    # 출력 좌표 -> 원본 좌표
    sx = np.clip((-np.sin(a) * _YY + np.cos(a) * _XX) / zx + (N - 1) / 2, 0, N - 1.001)
    y0, x0 = sy.astype(int), sx.astype(int)
    fy, fx = (sy - y0).reshape(len(stack), 1, -1), (sx - x0).reshape(len(stack), 1, -1)
    flat = stack.reshape(len(stack), stack.shape[1], -1)
    lin = (y0 * N + x0).reshape(len(stack), 1, -1)
    g = lambda d: np.take_along_axis(flat, np.broadcast_to(lin + d, (len(stack), stack.shape[1], lin.shape[2])), axis=2)
    out = g(0) * (1 - fy) * (1 - fx) + g(N) * fy * (1 - fx) + g(1) * (1 - fy) * fx + g(N + 1) * fy * fx
    return out.reshape(stack.shape)


def preprocess_batch(O, S, M):
    """여러 장을 한 번에 전처리한다. 0) 결측 보정 1) 해상도·비스듬 촬영 보정(표준 0.5m/픽셀)
    2) 구름 차폐 영역을 배경값으로 채움 3) 객체 주축을 가로로 회전 4) 0~1 정규화."""
    src = np.stack([np.asarray(O, float), np.asarray(S, float)], axis=1)
    if np.isnan(src).any():                                                               # 0) 센서 결측 픽셀 보정
        mean = np.nanmean(src, axis=(2, 3), keepdims=True)
        src = np.where(np.isnan(src), mean, src)
    zy = np.array([m["gsd"] / STD_GSD / np.cos(np.radians(m["off_nadir"])) for m in M])
    zx = np.array([m["gsd"] / STD_GSD for m in M])
    out = warp(src, zy, zx, np.zeros(len(M)))
    sm = ndimage.uniform_filter(out[:, 1], (1, 3, 3))
    angs, cloud = np.zeros(len(M)), np.zeros(len(M))
    for i in range(len(M)):
        cm = cloud_mask(out[i, 0], sm[i])
        if cm.any():
            fill = out[i, 0][~cm].mean() if (~cm).any() else 115.0
            out[i, 0][cm] = fill
            back = warp(cm[None, None].astype(float), [1 / zy[i]], [1 / zx[i]], np.zeros(1))[0, 0] > 0.5
            src[i, 0][back] = fill                                                       # 원본의 같은 구름 영역도 채움
            cloud[i] = cm.mean()
        angs[i] = orientation(out[i, 0], sm[i])[0]
    rot = np.abs(angs) > 3                                                                # 3도 이하는 회전 생략
    if rot.any():
        out[rot] = warp(src[rot], zy[rot], zx[rot], angs[rot])
    return [(np.clip(o / 255, 0, 1), np.clip(s / 1.5, 0, 1), float(c)) for (o, s), c in zip(out, cloud)]


def preprocess(opt, sar, meta):
    """한 장 전처리(preprocess_batch와 같은 처리)."""
    return preprocess_batch([opt], [sar], [meta])[0]

# ---- 특징 추출 공통 함수 ----
KERNELS = [np.array(k, float) for k in (
    [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], [[-1, -2, -1], [0, 0, 0], [1, 2, 1]], [[0, 1, 0], [1, -4, 1], [0, 1, 0]],
    [[2, -1, -1], [-1, 2, -1], [-1, -1, 2]], [[-1, -1, 2], [-1, 2, -1], [2, -1, -1]], [[1, 1, 1], [1, 1, 1], [1, 1, 1]])]


def conv_features(o):
    """간이 CNN 특징: 고정 합성곱 필터 6종 -> ReLU -> 4x4 최대 풀링(8x8 출력)."""
    feats = []
    for k in KERNELS:
        r = np.maximum(ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0)
        feats.append(r.reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
        feats.append(np.maximum(-ndimage.correlate(o - o.mean(), k / np.abs(k).sum()), 0).reshape(8, 4, 8, 4).max(axis=(1, 3)).ravel())
    return np.concatenate(feats)


def hog(o, cell=8, bins=8):
    gy, gx = np.gradient(o)
    mag, ang = np.hypot(gx, gy), np.mod(np.arctan2(gy, gx), np.pi)
    b = np.minimum((ang / np.pi * bins).astype(int), bins - 1)
    out = np.zeros((N // cell, N // cell, bins))
    for i in range(N // cell):
        for j in range(N // cell):
            sl = (slice(i * cell, (i + 1) * cell), slice(j * cell, (j + 1) * cell))
            out[i, j] = np.bincount(b[sl].ravel(), mag[sl].ravel(), bins)
    return (out / (np.linalg.norm(out) + 1e-6)).ravel()


def shape_features(o, s):
    """객체 모양 특징: 가로 방향 투영으로 전체 길이·차체 길이·가는 돌출부(포신) 길이·두꺼운 부분(포탑) 위치를 잰다.
    방위 정규화 뒤에도 앞뒤(180도)가 뒤집힐 수 있으므로, 가는 돌출부가 오른쪽에 오도록 맞춘다."""
    sm = ndimage.uniform_filter(s, 3)
    dark = np.abs(o - np.median(o))
    m = (sm > 0.22) | (dark > 0.12)
    lab, n = ndimage.label(m)
    if n:
        c = lab[10:22, 8:24]
        k = np.bincount(c[c > 0]).argmax() if (c > 0).any() else 1
        m = lab == k
    p = m.sum(axis=0).astype(float)
    cols = np.nonzero(p)[0]
    if len(cols) == 0:
        return np.zeros(8), o, s
    thin = (p > 0) & (p <= 2)
    left, right = thin[:N // 2].sum(), thin[N // 2:].sum()
    if left > right:
        o, s, m, p = o[:, ::-1], s[:, ::-1], m[:, ::-1], p[::-1]
        cols = np.nonzero(p)[0]
    body = np.nonzero(p > 3)[0]
    d = np.where(m, dark if left <= right else dark[:, ::-1], 0)
    heavy = d.sum(axis=0)
    centre = (body[0] + body[-1]) / 2 if len(body) else N / 2
    turret = (heavy * np.arange(N)).sum() / max(heavy.sum(), 1e-6) - centre
    return np.array([cols[-1] - cols[0] + 1, len(body), ((p > 0) & (p <= 2)).sum(), p.max(), turret,
                     m.sum(), (cols[-1] - (body[-1] if len(body) else cols[-1])), m.sum() / max(len(cols) * p.max(), 1)]), o, s


def adv_features(o, s):
    """정밀 특징: 모양(길이·포신·포탑 위치) + HOG(윤곽) + SAR 반사 통계(위장막·구름 영향 적음) + 투영 윤곽."""
    shp, o, s = shape_features(o, s)
    c = s[8:24, 4:28]
    sm = ndimage.uniform_filter(s, 3)
    m = sm > 0.2
    dark = np.abs(o - np.median(o))
    return np.concatenate([shp, hog(o), [c.mean(), c.max(), (sm > 0.25).mean(), (sm > 0.4).mean(), m.sum(axis=1).max(), m.sum(axis=0).max()],
                           sm.mean(axis=0).reshape(8, 4).mean(1), sm.mean(axis=1).reshape(8, 4).mean(1),
                           dark.mean(axis=0).reshape(8, 4).mean(1), dark.mean(axis=1).reshape(8, 4).mean(1),
                           [o.mean(), o.std(), dark.max(), (o > 0.65).mean(), (o < 0.3).mean()]])

CONFIG = {"train_with_decoy": True, "high": 0.9, "medium": 0.6, "hvt_types": ["spg"]}     # 미션 9 재보정값
CALLED = set()


def cover(fn):
    """기능 커버리지 측정용 표시: 통합 실행 중 호출된 기능을 기록한다."""
    def wrap(*a, **k):
        CALLED.add(fn.__name__)
        return fn(*a, **k)
    wrap.__name__ = fn.__name__
    return wrap


preprocess, adv_features, shape_features, hog, cloud_mask, orientation = map(cover, (preprocess, adv_features, shape_features, hog, cloud_mask, orientation))


@cover
def validate_tile(o, s, m):
    errs = []
    if o.shape != (N, N) or s.shape != (N, N):
        errs.append("타일 크기 오류")
    if not 0.1 <= m["gsd"] <= 30:
        errs.append("gsd 범위 오류")
    if not np.isfinite(o).all():
        errs.append("결측 픽셀")
    return errs


@cover
def quality(P, prob):
    srt = np.sort(prob, axis=1)
    ent = -(prob * np.log(np.clip(prob, 1e-9, 1))).sum(1) / np.log(prob.shape[1])
    return np.c_[srt[:, -1], srt[:, -1] - srt[:, -2], ent, [p[2] for p in P], [p[1][8:24, 4:28].mean() for p in P]]


@cover
def grade(c):
    return "High" if c >= CONFIG["high"] else ("Medium" if c >= CONFIG["medium"] else "Low")


class Queue:
    def __init__(self):
        self.items, self.log = [], []

    @cover
    def submit(self, iid, pred, g, conf):
        prev = self.log[-1]["hash"] if self.log else "GENESIS"
        body = {"event": "SUBMIT", "id": iid, "pred": pred, "grade": g, "conf": round(conf, 3), "hvt": pred in CONFIG["hvt_types"]}
        self.log.append({**body, "prev": prev, "hash": hashlib.sha256((prev + json.dumps(body, sort_keys=True)).encode()).hexdigest()})
        self.items.append({**body, "status": "PENDING"})

    @cover
    def verify(self):
        prev = "GENESIS"
        for r in self.log:
            body = {k: v for k, v in r.items() if k not in ("prev", "hash")}
            if r["prev"] != prev or r["hash"] != hashlib.sha256((prev + json.dumps(body, sort_keys=True)).encode()).hexdigest():
                return False
            prev = r["hash"]
        return True


TARGETS = {"validate_tile", "preprocess", "cloud_mask", "orientation", "adv_features", "shape_features", "hog", "quality", "grade", "submit", "verify", "dashboard"}
LABELS = CLASSES + ["decoy", "none"]

# ---------------- 1. 학습(분류 모델 + 신뢰도 가중치) ----------------
t_start = time.perf_counter()
O, S, M = make_dataset(200 * len(LABELS), seed=3, scen=tuple(SCEN), classes=LABELS)
clf = RandomForestClassifier(n_estimators=150, min_samples_leaf=2, random_state=0, n_jobs=1)
clf.fit(np.array([adv_features(o, s) for o, s, _ in preprocess_batch(O, S, M)]), [m["cls"] for m in M])
Oc, Sc, Mc = make_dataset(480, seed=6, scen=tuple(SCEN), classes=LABELS)
Pc = preprocess_batch(Oc, Sc, Mc)
prc = clf.predict_proba(np.array([adv_features(o, s) for o, s, _ in Pc]))
conf_model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)).fit(quality(Pc, prc), clf.classes_[prc.argmax(1)] == np.array([m["cls"] for m in Mc]))

# ---------------- 2. 운용 자료(가상 판독 지역 300타일) 처리 ----------------
Oo, So, Mo = make_dataset(300, seed=31, scen=("clean", "camo", "occl30", "cloud", "low_sun", "low_res"), classes=LABELS)
ops = list(zip(Oo, So, Mo))
ops[5] = (Oo[5][:20], So[5], Mo[5])                                 # 손상 타일 1건(잘림)
ops[9] = (Oo[9], So[9], {**Mo[9], "gsd": 0.0})                      # 메타데이터 오류 1건
queue, rejected, res, ms = Queue(), {}, [], []
for o, s, m in ops:
    errs = validate_tile(o, s, m)
    if errs:
        rejected[m["id"]] = errs; continue
    t0 = time.perf_counter()
    P = [preprocess(o, s, m)]
    pr = clf.predict_proba([adv_features(P[0][0], P[0][1])])
    pred = clf.classes_[pr.argmax()]
    conf = float(conf_model.predict_proba(quality(P, pr))[0, 1])
    g = grade(conf)
    if pred != "none":
        queue.submit(m["id"], pred, g, conf)
    ms.append((time.perf_counter() - t0) * 1000)
    res.append({"id": m["id"], "y": m["cls"], "p": pred, "grade": g, "conf": conf, "scen": m["scen"]})
R = pd.DataFrame(res)
acc = (R.y == R.p).mean()
milR = R[R.y.isin(list(MILITARY))]
rec = (milR.y == milR.p).mean()
decoy_catch = (R[R.y == "decoy"].p == "decoy").mean()
hi = R[(R.grade == "High") & R.p.isin(list(MILITARY))]
hi_err = (hi.y != hi.p).mean() if len(hi) else 0.0


# ---------------- 3. 분석 화면 ----------------
@cover
def dashboard():
    plt.style.use("dark_background")
    fig, ax = plt.subplots(2, 2, figsize=(12, 9))
    T = 10
    mosaic = np.block([[np.asarray(Oo[r * T + c]) for c in range(T)] for r in range(T)])
    ax[0, 0].imshow(mosaic, cmap="gray")
    for k, r in R[R.id.isin([f"IMG{i:04d}" for i in range(100)])].iterrows():
        i = int(r.id[3:])
        if r.p != "none":
            ax[0, 0].scatter(i % T * N + 16, i // T * N + 16, s=60, facecolors="none",
                             edgecolors="red" if r.p in MILITARY else ("yellow" if r.p == "decoy" else "deepskyblue"))
    ax[0, 0].set_title("First 100 tiles (red = military, yellow = decoy)"); ax[0, 0].axis("off")
    vc = R[R.p != "none"].p.value_counts()
    ax[0, 1].bar(vc.index.tolist(), vc.values, color="deepskyblue"); ax[0, 1].set_title("Detections by type")
    gc = R[R.p != "none"].grade.value_counts().reindex(["High", "Medium", "Low"]).fillna(0)
    ax[1, 0].bar(gc.index.tolist(), gc.values, color=["red", "orange", "gray"]); ax[1, 0].set_title("Review queue by grade (all pending)")
    ax[1, 1].axis("off")
    ax[1, 1].text(0, .9, f"Accuracy      {acc:.1%}\nMilitary rec. {rec:.1%}\nDecoy caught  {decoy_catch:.1%}\nQueue         {len(queue.items)} pending\n"
                  f"Latency       {np.mean(ms):.1f} ms/tile", family="monospace", fontsize=13, va="top")
    fig.suptitle("ISR01 integrated satellite object classification - synthetic data only"); plt.tight_layout(); plt.show()


dashboard()

# ---------------- 4. 자체 점검(통합 테스트) ----------------
log_ok = queue.verify()
cov = len(CALLED & TARGETS) / len(TARGETS)
chk = {
    "손상 타일·메타데이터 오류 차단": set(rejected) == {"IMG0005", "IMG0009"},
    "처리 타일 298장": len(res) == 298,
    "종합 정확도 90% 이상": acc >= 0.9,
    "군용 기종 식별 재현율 85% 이상": rec >= 0.85,
    "모의 표적 적발률 80% 이상": decoy_catch >= 0.8,
    "자동 승인 0건(전부 판독 대기)": all(i["status"] == "PENDING" for i in queue.items),
    "식별 건 전량 대기열 등록": len(queue.items) == int((R.p != "none").sum()),
    "감사 로그 무결성": log_ok,
    "타일당 200ms 이내": float(np.mean(ms)) < 200,
    "기능 커버리지 85% 이상": cov >= 0.85,
}
passed = sum(chk.values())
for k, v in chk.items():
    print(f"[{'PASS' if v else 'FAIL'}] {k}")
print(f"자체 점검 {passed}/{len(chk)} 통과, 기능 커버리지 {cov:.0%} ({len(CALLED & TARGETS)}/{len(TARGETS)}), 파이프라인 실행 성공 {'예' if passed == len(chk) else '아니오'}")

# ---------------- 5. 전략 ISR 평가 브리핑 요약 ----------------
print(f"""
[위성영상 장비 객체분류 평가 브리핑 요약 - 교육용가상자료]
1. 처리: 타일 {len(ops)}장 중 {len(res)}장 판독, 손상·오류 {len(rejected)}건 차단({', '.join(rejected)})
2. 성능: 종합 정확도 {acc:.1%}, 군용 기종 식별 재현율 {rec:.1%}, 모의 표적 적발 {decoy_catch:.1%}, 타일당 {np.mean(ms):.1f} ms, 전체 {time.perf_counter() - t_start:.0f} s
3. 판독 대기: {len(queue.items)}건 (High {int((R[R.p != 'none'].grade == 'High').sum())} / Medium {int((R[R.p != 'none'].grade == 'Medium').sum())} / Low {int((R[R.p != 'none'].grade == 'Low').sum())}), HVT 후보 {sum(i['hvt'] for i in queue.items)}건, 자동 승인 0건
4. 신뢰도: High 등급 군용 식별 {len(hi)}건 중 오분류 {hi_err:.1%} -> 등급은 판독 순서 결정용이며 최종 판독은 영상분석관이 한다
5. 한계: 합성 32x32 영상 기준 결과로, 실제 위성영상·장비에는 재학습과 재검증이 필요하다. 인접 체계 연동은 판독 승인 결과만 표준 필드로 전달한다.""")

실행 결과 예시

[PASS] 손상 타일·메타데이터 오류 차단
[PASS] 처리 타일 298장
[PASS] 종합 정확도 90% 이상
[PASS] 군용 기종 식별 재현율 85% 이상
[PASS] 모의 표적 적발률 80% 이상
[PASS] 자동 승인 0건(전부 판독 대기)
[PASS] 식별 건 전량 대기열 등록
[PASS] 감사 로그 무결성
[PASS] 타일당 200ms 이내
[PASS] 기능 커버리지 85% 이상
자체 점검 10/10 통과, 기능 커버리지 100% (12/12), 파이프라인 실행 성공 예

[위성영상 장비 객체분류 평가 브리핑 요약 - 교육용가상자료]
1. 처리: 타일 300장 중 298장 판독, 손상·오류 2건 차단(IMG0005, IMG0009)
2. 성능: 종합 정확도 96.6%, 군용 기종 식별 재현율 96.0%, 모의 표적 적발 100.0%, 타일당 14.0 ms, 전체 10 s
3. 판독 대기: 251건 (High 223 / Medium 21 / Low 7), HVT 후보 51건, 자동 승인 0건
4. 신뢰도: High 등급 군용 식별 133건 중 오분류 2.3% -> 등급은 판독 순서 결정용이며 최종 판독은 영상분석관이 한다
5. 한계: 합성 32x32 영상 기준 결과로, 실제 위성영상·장비에는 재학습과 재검증이 필요하다. 인접 체계 연동은 판독 승인 결과만 표준 필드로 전달한다.