{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "91bd5e5c",
   "metadata": {},
   "source": [
    "# KNPS 2026 불법시설물 변화탐지 베이스라인 제출 (UNet-ResNet18, 6채널)\n",
    "\n",
    "같은 위치의 전(pre)/후(post) 위성 RGB 패치(256x256) 한 쌍에서 **증축(new_building)** 과\n",
    "**벌목(tree_removal)** 영역을 각각 폴리곤 목록으로 출력하는 참고 구현입니다. 성능 기준선이 아니라\n",
    "**입력에서 제출 파일까지의 흐름과 제출 규격을 확인하는 최소 동작 예시**이며, 사용 의무는 없습니다.\n",
    "\n",
    "- 입력: `AIF_INPUT_DIR` 아래 `pairs.csv`(`id` 열)와 `images/<id>/pre.png`, `images/<id>/post.png`\n",
    "  (각 256x256 RGB 8bit, 두 영상은 같은 격자에 정렬되어 있습니다)\n",
    "- 출력: `AIF_PREDICTION_PATH`에 `id,new_building,tree_removal` CSV 한 개. 두 열은 각각 그 클래스의 변화\n",
    "  폴리곤 목록을 적은 JSON 문자열입니다. 폴리곤은 `[x, y]` 점의 목록(구멍 없음)이고 목록 안에 폴리곤을 여러 개\n",
    "  담을 수 있습니다. 예: `[[[10,10],[30,10],[30,30],[10,30]]]`. 좌표는 후 영상 패치의 픽셀 좌표(원점 왼쪽 위\n",
    "  모서리, x 오른쪽, y 아래, 실수 허용)이고 빈 문자열은 그 클래스의 변화가 없음을 뜻합니다. CSV 셀 안의 JSON 은\n",
    "  `csv` 모듈이 인용 처리를 해 줍니다\n",
    "- 동봉 가중치: `assets/model/unet_r18_cd.pt` (55MB). 추론 단계는 네트워크가 차단되므로 필요한\n",
    "  파일은 모두 zip 안에 넣어야 합니다\n",
    "- 의존성: `requirements.txt`로 선언합니다(폴리곤 변환에 shapely 를 씁니다). 채점 서버가 설치 단계에서 먼저 설치하므로 노트북 안에\n",
    "  설치 셀을 두지 않았습니다. `%pip install -r requirements.txt`처럼 상대경로로 파일을 가리키는\n",
    "  설치 줄은 설치 단계의 작업 폴더가 제출물 폴더가 아니어서 실패하므로 넣지 마시기 바랍니다\n",
    "- 제출 규칙: zip 최상위에 `predict.ipynb`(이 파일)와 `requirements.txt`, 가중치는 `assets/`에 동봉\n",
    "\n",
    "**로컬 실행**: `pairs.csv`와 `images/<id>/{pre,post}.png`를 담은 폴더를 만들고 환경변수로 넘기면 됩니다.\n",
    "지정하지 않으면 `./input`을 읽고 `./prediction.csv`에 씁니다.\n",
    "```\n",
    "AIF_INPUT_DIR=<입력 폴더> AIF_PREDICTION_PATH=./prediction.csv jupyter nbconvert --to notebook --execute predict.ipynb\n",
    "```\n",
    "\n",
    "**라이선스**: segmentation_models_pytorch MIT, ResNet18 ImageNet 인코더(torchvision 원본) BSD-3,\n",
    "파인튜닝 가중치와 이 코드는 대회 참가 목적으로 자유롭게 사용, 수정, 재학습하실 수 있습니다.\n",
    "파생 제출물에 카피레프트나 별도 약관이 전파되지 않습니다. 출처는 `NOTICE`, 원문은 `LICENSE`에\n",
    "있습니다. 이 베이스라인은 참고 구현이며 사용 의무가 없습니다."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5f3a0ff",
   "metadata": {},
   "outputs": [],
   "source": "# [setup]\nimport csv\nimport io\nimport json\nimport os\nimport sys\nimport time\nfrom pathlib import Path\n\nimport numpy as np\nimport torch\nfrom PIL import Image\nfrom shapely.geometry import MultiPolygon, Polygon, box\nfrom shapely.ops import unary_union\n\nINPUT_DIR = Path(os.environ.get(\"AIF_INPUT_DIR\", \"./input\"))\nPREDICTION_PATH = Path(os.environ.get(\"AIF_PREDICTION_PATH\", \"./prediction.csv\"))\nCKPT = next(p for p in map(Path, [\"assets/model/unet_synth.pt\", \"assets/model/siam_v2.pt\", \"assets/model/hybrid_v3.pt\"]) if p.exists())      # 커널의 작업 폴더가 제출물 루트이므로 상대경로로 지정합니다\n\nH, W = 256, 256\nBATCH = 16\nCLASSES = (\"new_building\", \"tree_removal\")       # 모델 출력 채널 1, 2 (0 은 배경)\nMIN_AREA = {'new_building': 30.0, 'tree_removal': 30.0}            # 이 면적(px^2) 미만의 조각은 지웁니다 (학습 때와 같은 후처리)\nMIN_POS_AREA = 20.0      # 조각을 지운 뒤 남은 총면적이 이 값 미만이면 빈 예측으로 냅니다 (채점 규칙과 동일)\nSIMPLIFY_PX = 0.5        # 계단형 외곽선을 이 허용 오차로 단순화합니다 (위상 보존). 채점 허용 오차 1px 안입니다\nNDIGITS = 2              # 좌표 소수 자릿수\nMEAN = np.array([0.485, 0.456, 0.406], np.float32)\nSTD = np.array([0.229, 0.224, 0.225], np.float32)\n\n\ndef mask_to_polygons(mask: np.ndarray, min_area: float) -> str:\n    \"\"\"(H,W) 이진 마스크 -> 폴리곤 목록 JSON 문자열(폴리곤 = 점 목록, 구멍 없음). 변화가 없으면 빈 문자열.\n\n    행 단위 런(run)을 화소 경계 사각형으로 만들어 합집합하면 마스크의 외곽선을 그대로 따르는 폴리곤이 나옵니다.\n    좌표는 화소 모서리 기준(화소 (r, c) 는 [c, c+1] x [r, r+1])이라 채점 규약과 같습니다.\n    \"\"\"\n    m = np.asarray(mask, dtype=bool)\n    boxes = []\n    for r in np.flatnonzero(m.any(axis=1)):\n        pad = np.concatenate(([0], m[r].astype(np.int8), [0]))\n        edges = np.flatnonzero(np.diff(pad))\n        for s, e in zip(edges[::2], edges[1::2]):\n            boxes.append(box(float(s), float(r), float(e), float(r + 1)))\n    if not boxes:\n        return \"\"\n    g = unary_union(boxes)\n    parts = list(g.geoms) if isinstance(g, MultiPolygon) else [g]\n    parts = [p for p in parts if isinstance(p, Polygon) and p.area >= min_area]     # 작은 조각 제거\n    if not parts or sum(p.area for p in parts) < MIN_POS_AREA:\n        return \"\"\n    out = []\n    for p in parts:\n        p = p.simplify(SIMPLIFY_PX, preserve_topology=True)\n        if p.is_empty or p.area <= 0:\n            continue\n        # 바깥 경계만 씁니다(제출 형식에 구멍이 없습니다). 마지막 점은 첫 점과 같으므로 뺍니다\n        out.append([[round(float(x), NDIGITS), round(float(y), NDIGITS)] for x, y in list(p.exterior.coords)[:-1]])\n    return json.dumps(out, separators=(\",\", \":\")) if out else \"\"\n\nTHRESH = {'new_building': 0.5, 'tree_removal': 0.5}   # pixel probability thresholds\nPRES_THRESH = 0.0                                           # presence-head gate (siam only, 0 = off)\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "938bc30a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# [inputs]\n",
    "# 목록 파일을 재귀 탐색으로 찾고 그 부모 폴더를 입력 루트로 삼습니다 (껍데기 폴더 한 겹에 대비).\n",
    "# 추론 대상은 목록 파일이 정본이며 images/ 폴더를 훑어 대상을 정하지 않습니다.\n",
    "print(\"AIF_INPUT_DIR:\", INPUT_DIR, sorted(p.name for p in INPUT_DIR.iterdir()) if INPUT_DIR.is_dir() else \"(없음)\")\n",
    "lists = sorted(INPUT_DIR.rglob(\"pairs.csv\"))\n",
    "if not lists:\n",
    "    # 입력이 풀린 폴더가 아니라 압축본 그대로 주어지는 경우에 대비합니다\n",
    "    zs = sorted(INPUT_DIR.rglob(\"*.zip\")) if INPUT_DIR.is_dir() else []\n",
    "    if len(zs) == 1:\n",
    "        import tempfile, zipfile\n",
    "        dst = Path(tempfile.mkdtemp()) / \"input\"\n",
    "        with zipfile.ZipFile(zs[0]) as z:\n",
    "            z.extractall(dst)\n",
    "        print(f\"입력이 압축본이어서 풀어서 읽습니다: {zs[0].name} -> {dst}\")\n",
    "        lists = sorted(dst.rglob(\"pairs.csv\"))\n",
    "if not lists:\n",
    "    raise SystemExit(f\"pairs.csv 를 찾지 못했습니다: {INPUT_DIR}\")\n",
    "ROOT = lists[0].parent\n",
    "with lists[0].open(encoding=\"utf-8-sig\", newline=\"\") as f:\n",
    "    IDS = [row[\"id\"].strip() for row in csv.DictReader(f) if row.get(\"id\", \"\").strip()]\n",
    "print(f\"입력 루트 {ROOT}, 쌍 {len(IDS)}건\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ab4f29d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# [env]\n",
    "print(\"python\", sys.version.split()[0], \"torch\", torch.__version__)\n",
    "print(\"cuda\", torch.cuda.is_available(), \"mps\", getattr(torch.backends, \"mps\", None) is not None and torch.backends.mps.is_available())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "400e6923",
   "metadata": {},
   "outputs": [],
   "source": [
    "# [model]\n",
    "sys.path.insert(0, \"assets\")\n",
    "import segmentation_models_pytorch as smp\n",
    "from model_v2 import SiamCD\n",
    "ENS_W = 0.5   # weight of the baseline-arch UNet; 1-ENS_W for the Siamese model\n",
    "\n",
    "\n",
    "def load_model(device: str):\n",
    "    u = smp.Unet(encoder_name=\"resnet18\", encoder_weights=None, in_channels=6, classes=3)\n",
    "    u.load_state_dict(torch.load(\"assets/model/unet_synth.pt\", map_location=\"cpu\", weights_only=True)[\"state_dict\"])\n",
    "    import glob\n",
    "    sm = []\n",
    "    for f in sorted(glob.glob(\"assets/model/siam_*.pt\")):          # several Satlas Siamese models -> output average\n",
    "        s = SiamCD()\n",
    "        s.load_state_dict(torch.load(f, map_location=\"cpu\", weights_only=True)[\"state_dict\"])\n",
    "        sm.append(s.to(device).eval()); print(\"loaded\", f, flush=True)\n",
    "    return (u.to(device).eval(), sm)\n",
    "\n",
    "\n",
    "def read_image(path: Path) -> np.ndarray:\n",
    "    with Image.open(io.BytesIO(path.read_bytes())) as im:\n",
    "        arr = np.asarray(im.convert(\"RGB\"))\n",
    "    if arr.shape[:2] != (H, W):\n",
    "        raise ValueError(f\"{path}: 크기 {arr.shape[:2]} 가 {(H, W)} 와 다릅니다\")\n",
    "    return arr\n",
    "\n",
    "\n",
    "@torch.no_grad()\n",
    "def predict_batch(model, pres, posts, device):\n",
    "    u, s = model\n",
    "    a = torch.from_numpy(np.stack(pres)).permute(0, 3, 1, 2).to(device).float() / 255\n",
    "    b = torch.from_numpy(np.stack(posts)).permute(0, 3, 1, 2).to(device).float() / 255\n",
    "    m = torch.tensor(MEAN, device=device).view(1, 3, 1, 1); sd = torch.tensor(STD, device=device).view(1, 3, 1, 1)\n",
    "    prob, pp = 0, 0\n",
    "    for dims in ([], [3], [2], [2, 3]):\n",
    "        ai, bi = (torch.flip(a, dims), torch.flip(b, dims)) if dims else (a, b)\n",
    "        pu = torch.softmax(u(torch.cat([(ai - m) / sd, (bi - m) / sd], 1)), 1)[:, 1:3]\n",
    "        cs, ps = 0, 0\n",
    "        for sk in s:\n",
    "            with torch.autocast(\"cuda\", dtype=torch.float16, enabled=(device == \"cuda\")):\n",
    "                ch, pr, _, _ = sk(ai, bi)\n",
    "            cs = cs + ch.float().sigmoid() / len(s); ps = ps + pr.float().sigmoid() / len(s)\n",
    "        ci = ENS_W * pu.float() + (1 - ENS_W) * cs\n",
    "        prob = prob + (torch.flip(ci, dims) if dims else ci) / 4\n",
    "        pp = pp + ps / 4\n",
    "    return prob.cpu().numpy(), pp.cpu().numpy()\n",
    "\n",
    "\n",
    "DEVICE = \"cuda\" if torch.cuda.is_available() else (\"mps\" if torch.backends.mps.is_available() else \"cpu\")\n",
    "try:\n",
    "    MODEL = load_model(DEVICE)\n",
    "    predict_batch(MODEL, [np.zeros((H, W, 3), np.uint8)], [np.zeros((H, W, 3), np.uint8)], DEVICE)\n",
    "except Exception as e:  # noqa: BLE001\n",
    "    print(f\"{DEVICE} 시험 추론에 실패해({type(e).__name__}: {e}) CPU 로 전환합니다\")\n",
    "    DEVICE = \"cpu\"\n",
    "    MODEL = load_model(DEVICE)\n",
    "print(\"device\", DEVICE)\n",
    "\n",
    "\n",
    "# ---- CHMv2 canopy height (DINOv3 ViT-L/16 SAT-493M + DPT head, DINOv3 License: assets/chmv2/LICENSE_DINOv3.md)\n",
    "CHM_A = 0.2       # weight of the canopy-height cue in the tree_removal probability\n",
    "try:\n",
    "    import transformers\n",
    "    print(\"python\", sys.version.split()[0], \"transformers\", transformers.__version__, flush=True)\n",
    "    from chm_fuse import CHM, cue_abs, cue_diff\n",
    "    CHM_CUE = cue_diff\n",
    "    CHM_MODEL = CHM(\"assets/chmv2\", DEVICE, scale=1.5)\n",
    "    CHM_MODEL.heights([np.zeros((H, W, 3), np.uint8)])\n",
    "    print(\"CHM loaded\", CHM_A, \"diff\", 1.5, flush=True)\n",
    "except Exception as e:  # noqa: BLE001\n",
    "    print(f\"CHM load failed ({type(e).__name__}: {e}); CHM cue off\", flush=True)\n",
    "    CHM_A = 0.0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a918ddd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# [infer]\n",
    "ROWS = []\n",
    "STATS = {c: {\"maxp\": [], \"pres\": []} for c in CLASSES}\n",
    "GRID_T = (0.3, 0.4, 0.5, 0.6, 0.7, 0.8)\n",
    "GRID_A = (20, 50, 100, 200)\n",
    "GRID = {c: np.zeros((len(GRID_T), len(GRID_A)), int) for c in CLASSES}\n",
    "t0 = time.time()\n",
    "for s in range(0, len(IDS), BATCH):\n",
    "    ids = IDS[s:s + BATCH]\n",
    "    pres = [read_image(ROOT / \"images\" / i / \"pre.png\") for i in ids]\n",
    "    posts = [read_image(ROOT / \"images\" / i / \"post.png\") for i in ids]\n",
    "    prob, pp = predict_batch(MODEL, pres, posts, DEVICE)\n",
    "    if CHM_A > 0:\n",
    "        cue = CHM_CUE(CHM_MODEL.heights(pres), CHM_MODEL.heights(posts)).cpu().numpy()\n",
    "        prob[:, 1] = (1 - CHM_A) * prob[:, 1] + CHM_A * cue\n",
    "        STATS.setdefault(\"cue_max\", []).extend(cue.reshape(len(ids), -1).max(1).tolist())\n",
    "    nd = np.stack([(a.max(2) == 0) | (b.max(2) == 0) for a, b in zip(pres, posts)])   # no-data is never scored\n",
    "    prob[:, :, :, :][np.repeat(nd[:, None], 2, 1)] = 0\n",
    "    for j, i in enumerate(ids):\n",
    "        cells = []\n",
    "        for k, c in enumerate(CLASSES):\n",
    "            p = prob[j, k]\n",
    "            STATS[c][\"maxp\"].append(float(p.max()))\n",
    "            if pp is not None:\n",
    "                STATS[c][\"pres\"].append(float(pp[j, k]))\n",
    "            for ti, t in enumerate(GRID_T):\n",
    "                area = int((p > t).sum())\n",
    "                for ai, amin in enumerate(GRID_A):\n",
    "                    GRID[c][ti, ai] += area >= amin\n",
    "            gate = pp is None or PRES_THRESH <= 0 or pp[j, k] >= PRES_THRESH\n",
    "            cells.append(mask_to_polygons(p > THRESH[c], MIN_AREA[c]) if gate else \"\")\n",
    "        ROWS.append((i, *cells))\n",
    "    if (s // BATCH) % 10 == 0 or s + BATCH >= len(IDS):\n",
    "        print(f\"  {min(s + BATCH, len(IDS))}/{len(IDS)}  {time.time() - t0:.1f}s\", flush=True)\n",
    "n_b = sum(1 for _, b, _ in ROWS if b)\n",
    "n_t = sum(1 for _, _, t in ROWS if t)\n",
    "print(f\"추론을 완료했습니다: {len(ROWS)}건, 증축 양성 {n_b}건, 벌목 양성 {n_t}건, {time.time() - t0:.1f}s\")\n",
    "print(\"설정:\", \"THRESH\", THRESH, \"MIN_AREA\", MIN_AREA, \"PRES_THRESH\", PRES_THRESH)\n",
    "for c in CLASSES:\n",
    "    mp = np.array(STATS[c][\"maxp\"])\n",
    "    print(f\"[{c}] 쌍별 최대확률 분위수(10/25/50/75/90%):\", np.round(np.percentile(mp, [10, 25, 50, 75, 90]), 3).tolist())\n",
    "    if STATS[c][\"pres\"]:\n",
    "        pr = np.array(STATS[c][\"pres\"])\n",
    "        print(f\"[{c}] 존재확률 분위수(10/25/50/75/90%):\", np.round(np.percentile(pr, [10, 25, 50, 75, 90]), 3).tolist(),\n",
    "              \" >=0.3/0.5/0.7:\", [int((pr >= v).sum()) for v in (0.3, 0.5, 0.7)])\n",
    "    print(f\"[{c}] 양성 쌍 수 표 (행=픽셀임계 {GRID_T}, 열=최소면적 {GRID_A}):\")\n",
    "    for ti, t in enumerate(GRID_T):\n",
    "        print(\"   \", t, GRID[c][ti].tolist())\n",
    "\n",
    "if STATS.get(\"cue_max\"):\n",
    "    print(\"CHM 단서 쌍별 최대값 분위수(10/50/90%):\", np.round(np.percentile(STATS[\"cue_max\"], [10, 50, 90]), 3).tolist())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bd390bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# [save]\n",
    "PREDICTION_PATH.parent.mkdir(parents=True, exist_ok=True)\n",
    "with PREDICTION_PATH.open(\"w\", newline=\"\", encoding=\"utf-8\") as f:\n",
    "    w = csv.writer(f)\n",
    "    w.writerow([\"id\", *CLASSES])\n",
    "    w.writerows(ROWS)\n",
    "print(\"저장했습니다:\", PREDICTION_PATH, PREDICTION_PATH.stat().st_size, \"bytes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "208f89da",
   "metadata": {},
   "source": [
    "## 제출\n",
    "\n",
    "아래 셀의 `<참여키>` 자리에 본인의 참여키를 넣고 실행하면 이 노트북이 들어 있는 폴더 전체가\n",
    "제출됩니다. 참여키를 노트북에 적은 채로 공유하지 마시기 바랍니다. 참여키는 마이페이지에서\n",
    "아래 경로를 따라 언제든 다시 확인하실 수 있습니다.\n",
    "\n",
    "프로필 → 마이페이지 → 활동히스토리 탭 → 참가하신 주제의 키 복사"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}