統計的背景モデルによる微細変化検出

REM Python をシステム領域にインストール
winget install --scope machine --id Python.Python.3.12 -e --silent
REM Python のパス設定
set "PYTHON_PATH=C:\Program Files\Python312"
set "PYTHON_SCRIPTS_PATH=C:\Program Files\Python312\Scripts"
echo "%PATH%" | find /i "%PYTHON_PATH%" >nul
if errorlevel 1 setx PATH "%PATH%;%PYTHON_PATH%" /M >nul
echo "%PATH%" | find /i "%PYTHON_SCRIPTS_PATH%" >nul
if errorlevel 1 setx PATH "%PATH%;%PYTHON_SCRIPTS_PATH%" /M >nul

winget install --scope machine Codeium.Windsurf -e --silent

# プログラム名: 統計的背景モデル（MOG）による微細変化検出プログラム
# 特徴技術名: 適応的統計背景モデル（混合ガウス：MOG, Stauffer & Grimson 1999）
# 出典: Stauffer, C., & Grimson, W. E. L. (1999). Adaptive background mixture models for real-time tracking. CVPR 1999
# 特徴機能: 各ピクセルをK成分の混合ガウスでモデル化し、ω/σに基づく背景集合とkσマッチ判定で照明変化に頑健な微細変化検出を実現
# 学習済みモデル: 使用なし
# 方式設計:
#   - 関連利用技術:
#     - OpenCV: 画像処理ライブラリ（形態学的処理、輪郭検出、色空間変換、CLAHE）
#     - NumPy: 数値計算ライブラリ（統計処理、行列演算）
#     - tkinter: GUIライブラリ（感度調整インターフェース）
#     - CLAHE: 局所的適応的ヒストグラム均等化による照明正規化
#     - 統計的背景モデル: 混合ガウス（MOG）による背景モデル構築
#   - 入力と出力: 入力: 動画（ユーザは「0:動画ファイル，1:カメラ，2:サンプル動画」のメニューで選択．0:動画ファイルの場合はtkinterでファイル選択．1の場合はOpenCVでカメラが開く．2の場合はhttps://raw.githubusercontent.com/opencv/opencv/master/samples/data/vtest.aviを使用）、出力: OpenCV画面でリアルタイム表示（変化領域を赤色オーバーレイ）、各フレームごとにFPSと変化率をprint()表示、終了時にresult.txtに保存
#   - 処理手順:
#     1. 入力フレームの取得とガウシアンフィルタ適用
#     2. 照明正規化処理（LabのLチャネルにCLAHE）
#     3. MOG（混合ガウス）による背景判定（ω/σ降順の累積重みTで背景集合を定義）
#     4. kσマッチ判定による変化検出（感度パラメータで制御）
#     5. 形態学的処理によるノイズ除去
#     6. 最小領域フィルタによる小領域除去
#   - 前処理、後処理:
#     - 前処理: ガウシアンフィルタによるノイズ除去、CLAHE照明正規化
#     - 後処理: 形態学的オープニング・クロージング、最小領域フィルタ
#   - 追加処理:
#     - 適応的背景更新: 学習率αとρ=α·ηに基づくMOG更新（S&G 1999）
#     - 統計的閾値: 標準偏差を用いたkσマッチ閾値
#     - バウンディングボックス描画: 検出領域の可視化と領域情報の表示
#   - 調整を必要とする設定値:
#     - sensitivity（基本感度）: 変化検出の閾値を制御（0.5-3.0、デフォルト1.5）
#     - min_change_area（最小変化領域）: ノイズ除去のための最小面積（10-200ピクセル、デフォルト50）
# 将来方策: sensitivity値の自動調整機能。初期30フレームの標準偏差分布を分析し、ヒストグラムの分散から背景の複雑さを定量化。分散値に基づいてsensitivity値を0.5-3.0の範囲で自動設定
# その他の重要事項: 微細な影や明るさ変化も検出可能な設計
# 前準備:
#   - pip install opencv-python numpy pillow

import cv2
import tkinter as tk
from tkinter import filedialog, Scale, Label, Button
import os
import numpy as np
import time
import threading
from PIL import Image, ImageDraw, ImageFont
from datetime import datetime
import urllib.request

# 設定値（調整可能なパラメータ）
SENSITIVITY = 1.5  # 変化検出の感度（0.5-3.0、値が大きいほど感度が低い）
MIN_CHANGE_AREA = 50  # 最小変化領域（ピクセル、小さいノイズを除去）
LEARNING_RATE = 0.01  # 背景モデルの学習率α（0.001-0.1、値が大きいほど背景の更新が速い）
BUFFER_SIZE = 30  # 背景モデル初期化用のフレーム数
GAUSSIAN_KERNEL = 5  # ガウシアンフィルタのカーネルサイズ（3,5,7など奇数）
MORPH_KERNEL = 3  # 形態学的処理のカーネルサイズ（3,5,7など奇数）
CLAHE_CLIP_LIMIT = 2.0  # CLAHEのクリップリミット（1.0-4.0）
CLAHE_GRID_SIZE = 8  # CLAHEのグリッドサイズ（通常8）
STD_OFFSET = 1.0  # 初期分散の下限（ゼロ除算防止）

# MOG設定
K_COMPONENTS = 5     # 混合成分数（3〜5が一般的）
BG_RATIO_T = 0.7     # 背景とみなす累積重みの閾値T（0.7〜0.9）
INIT_VAR = 225.0     # 初期分散（15^2）
MIN_VAR = 64.0       # 分散の下限（8^2）

# グローバル変数
sensitivity = SENSITIVITY
min_area = MIN_CHANGE_AREA
learning_rate = LEARNING_RATE
running = True
is_camera = False  # 入力種別のフラグ（choiceから設定）

# 背景モデル（MOG）
gmm_initialized = False
gmm_weights = None   # 形状: (H, W, K)
gmm_means = None     # 形状: (H, W, K, 3)
gmm_vars = None      # 形状: (H, W, K)
frame_buffer = []    # 初期化用フレームバッファ（Lab画像のfloat32）

# 結果記録用
frame_count = 0
results_log = []

# CLAHEオブジェクトを事前に作成
clahe = cv2.createCLAHE(clipLimit=CLAHE_CLIP_LIMIT, tileGridSize=(CLAHE_GRID_SIZE, CLAHE_GRID_SIZE))

# フォント設定
FONT_PATH = 'C:/Windows/Fonts/meiryo.ttc'
FONT_SIZE = 20
try:
    font = ImageFont.truetype(FONT_PATH, FONT_SIZE)
    use_japanese_font = True
except:
    use_japanese_font = False


def add_japanese_text(img, text, position, color=(0, 255, 0)):
    """OpenCV画像に日本語テキストを追加"""
    if not use_japanese_font:
        cv2.putText(img, text, position, cv2.FONT_HERSHEY_SIMPLEX, 0.7, color, 2)
        return img
    img_pil = Image.fromarray(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
    draw = ImageDraw.Draw(img_pil)
    draw.text(position, text, font=font, fill=color[::-1])  # BGRからRGBに変換
    return cv2.cvtColor(np.array(img_pil), cv2.COLOR_RGB2BGR)


def initialize_gmm_from_buffer():
    """BUFFER_SIZE枚のフレームからMOG初期化"""
    global gmm_initialized, gmm_weights, gmm_means, gmm_vars, frame_buffer
    stack = np.stack(frame_buffer, axis=0).astype(np.float32)  # (N, H, W, 3)
    mean0 = stack.mean(axis=0)                                 # (H, W, 3)
    var0 = np.maximum(stack.var(axis=0), STD_OFFSET)           # (H, W, 3)
    var0_scalar = var0.mean(axis=2)                            # (H, W) 等方分散
    H, W = mean0.shape[:2]
    gmm_weights = np.zeros((H, W, K_COMPONENTS), np.float32)
    gmm_weights[:, :, 0] = 1.0
    gmm_means = np.broadcast_to(mean0[:, :, None, :], (H, W, K_COMPONENTS, 3)).copy()
    gmm_vars = np.broadcast_to(var0_scalar[:, :, None], (H, W, K_COMPONENTS)).copy()
    gmm_vars = np.maximum(gmm_vars, MIN_VAR)
    gmm_initialized = True
    frame_buffer = []


def mog_foreground_and_update(lab_float):
    """MOG（S&G 1999）に基づく前景マスク生成とパラメータ更新"""
    global gmm_weights, gmm_means, gmm_vars
    # 差分と距離
    diff = lab_float[:, :, None, :] - gmm_means                # (H,W,K,3)
    d2 = np.sum(diff * diff, axis=3)                           # (H,W,K)
    sigma2 = np.maximum(gmm_vars, MIN_VAR)
    sigma = np.sqrt(sigma2)
    thr2 = (sensitivity * sigma) ** 2
    matches = d2 <= thr2                                       # (H,W,K) bool

    # 背景集合（更新前）：ω/σ降順の累積重みT
    rank = gmm_weights / np.maximum(sigma, 1e-6)
    order = np.argsort(-rank, axis=2)                          # (H,W,K)
    weights_sorted = np.take_along_axis(gmm_weights, order, axis=2)
    cumw = np.cumsum(weights_sorted, axis=2)
    bg_mask_sorted = cumw < BG_RATIO_T
    cross = cumw >= BG_RATIO_T
    first_idx = np.argmax(cross, axis=2)
    has_cross = cross.any(axis=2)
    rr, cc = np.where(has_cross)
    if rr.size > 0:
        bg_mask_sorted[rr, cc, first_idx[rr, cc]] = True
    inv = np.argsort(order, axis=2)
    bg_mask = np.take_along_axis(bg_mask_sorted, inv, axis=2)  # (H,W,K) bool

    # 前景判定（更新前の背景成分と一致していない画素）
    is_bg = (matches & bg_mask).any(axis=2)                    # (H,W)
    fg_mask = np.where(is_bg, 0, 255).astype(np.uint8)

    # 更新（Stauffer & Grimson 1999）
    alpha = learning_rate
    gmm_weights *= (1.0 - alpha)

    # マッチ成分のインデックス（最も重いもの）
    w_masked = np.where(matches, gmm_weights, -1.0)
    m_idx = np.argmax(w_masked, axis=2)                        # (H,W)
    matched_any = matches.any(axis=2)
    r_m, c_m = np.where(matched_any)
    if r_m.size > 0:
        k_m = m_idx[r_m, c_m]

        mu_old = gmm_means[r_m, c_m, k_m, :]                   # (M,3)
        var_old = np.maximum(gmm_vars[r_m, c_m, k_m], MIN_VAR) # (M,)
        x = lab_float[r_m, c_m, :]                             # (M,3)
        d = x - mu_old
        d2_m = np.sum(d * d, axis=1)                           # (M,)

        # ρ = α · η(x|μ, σ^2 I), d=3
        norm = ((2.0 * np.pi) ** 1.5) * (np.sqrt(var_old) ** 3)
        eta = np.exp(-0.5 * d2_m / var_old) / np.maximum(norm, 1e-12)
        rho = np.clip(alpha * eta, 0.0, 1.0)

        # μ, σ^2 の更新
        gmm_means[r_m, c_m, k_m, :] = mu_old + (rho[:, None] * d)
        gmm_vars[r_m, c_m, k_m] = np.clip(var_old + rho * (d2_m - var_old), MIN_VAR, 1e6)
        gmm_weights[r_m, c_m, k_m] += alpha

    # 正規化
    sumw = np.sum(gmm_weights, axis=2, keepdims=True)
    gmm_weights = gmm_weights / np.maximum(sumw, 1e-12)

    # 非マッチ画素は最小重み成分を置換
    unmatched = ~matched_any
    r_u, c_u = np.where(unmatched)
    if r_u.size > 0:
        repl = np.argmin(gmm_weights[r_u, c_u, :], axis=1)
        gmm_means[r_u, c_u, repl, :] = lab_float[r_u, c_u, :]
        gmm_vars[r_u, c_u, repl] = INIT_VAR
        gmm_weights[r_u, c_u, repl] = alpha
        sw = np.sum(gmm_weights[r_u, c_u, :], axis=1, keepdims=True)
        gmm_weights[r_u, c_u, :] /= np.maximum(sw, 1e-12)

    # ω/σ降順に並べ替え
    rank_new = gmm_weights / np.sqrt(np.maximum(gmm_vars, MIN_VAR))
    order_new = np.argsort(-rank_new, axis=2)
    gmm_weights = np.take_along_axis(gmm_weights, order_new, axis=2)
    gmm_vars = np.take_along_axis(gmm_vars, order_new, axis=2)
    gmm_means = np.take_along_axis(gmm_means, order_new[:, :, :, None], axis=2)

    return fg_mask


def video_frame_processing(frame):
    global frame_count, frame_buffer, gmm_initialized
    current_time = time.time()
    frame_count += 1

    # ガウシアンフィルタによるノイズ除去
    blurred = cv2.GaussianBlur(frame, (GAUSSIAN_KERNEL, GAUSSIAN_KERNEL), 0)
    # Lab色空間に変換
    lab = cv2.cvtColor(blurred, cv2.COLOR_BGR2Lab)

    # Lチャンネルに対してCLAHE適用（照明正規化）
    l_channel = lab[:, :, 0]
    lab[:, :, 0] = clahe.apply(l_channel)

    # 背景モデルの初期化
    if not gmm_initialized:
        frame_buffer.append(lab.astype(np.float32))
        n = len(frame_buffer)
        output = frame.copy()
        if n >= BUFFER_SIZE:
            initialize_gmm_from_buffer()
            msg = f'初期化完了: {BUFFER_SIZE}/{BUFFER_SIZE}'
        else:
            msg = f'初期化中... {n}/{BUFFER_SIZE}'
        output = add_japanese_text(output, msg, (10, 30))
        result = f'初期化ステータス: {msg}'
        return output, result, current_time

    # 前景マスク生成（MOG）と更新
    lab_float = lab.astype(np.float32)
    fg_mask = mog_foreground_and_update(lab_float)

    # 形態学的処理
    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (MORPH_KERNEL, MORPH_KERNEL))
    opened = cv2.morphologyEx(fg_mask, cv2.MORPH_OPEN, kernel)
    closed = cv2.morphologyEx(opened, cv2.MORPH_CLOSE, kernel)

    # 最小領域フィルタとバウンディングボックス
    contours, _ = cv2.findContours(closed, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    filtered_mask = np.zeros_like(closed)
    bounding_boxes = []

    for contour in contours:
        area = cv2.contourArea(contour)
        if area >= min_area:
            cv2.drawContours(filtered_mask, [contour], -1, 255, -1)
            x, y, w, h = cv2.boundingRect(contour)
            bounding_boxes.append((x, y, w, h, area))

    # 変化領域を赤色でオーバーレイ
    output = frame.copy()
    overlay = np.zeros_like(frame)
    overlay[:, :, 2] = filtered_mask
    output = cv2.addWeighted(output, 0.7, overlay, 0.3, 0)

    # バウンディングボックスを描画
    for x, y, w, h, area in bounding_boxes:
        cv2.rectangle(output, (x, y), (x + w, y + h), (0, 255, 0), 2)

    # 変化率計算
    change_pixels = cv2.countNonZero(filtered_mask)
    total_pixels = frame.shape[0] * frame.shape[1]
    change_rate = (change_pixels / total_pixels) * 100

    # テキスト表示
    output = add_japanese_text(output, f'変化率: {change_rate:.2f}%', (10, 30))
    output = add_japanese_text(output, f'感度: {sensitivity:.1f}', (10, 60))
    output = add_japanese_text(output, f'検出領域数: {len(bounding_boxes)}', (10, 90))

    # 結果文字列の作成
    result = f'変化率: {change_rate:.2f}%, 領域数: {len(bounding_boxes)}'

    return output, result, current_time


def update_sensitivity(value):
    global sensitivity
    sensitivity = float(value)


def video_processing_thread(cap):
    global running, frame_count, results_log
    try:
        prev_t = time.perf_counter()
        while running:
            ret, frame = cap.read()
            if not ret:
                break

            MAIN_FUNC_DESC = "微細変化検出"
            processed_frame, result, current_time = video_frame_processing(frame)
            cv2.imshow(MAIN_FUNC_DESC, processed_frame)

            now_t = time.perf_counter()
            dt = now_t - prev_t
            fps = (1.0 / dt) if dt > 0 else 0.0
            prev_t = now_t

            if is_camera:
                print(datetime.fromtimestamp(current_time).strftime("%Y-%m-%d %H:%M:%S.%f")[:-3],
                      f'FPS: {fps:.2f},', result)
            else:
                print(frame_count, f'FPS: {fps:.2f},', result)

            results_log.append(f'FPS: {fps:.2f}, {result}')

            if cv2.waitKey(1) & 0xFF == ord('q'):
                running = False
                request_gui_stop()
                break
    finally:
        cap.release()
        cv2.destroyAllWindows()


# ガイダンス表示
print('\n=== 統計的背景モデル（MOG）による微細変化検出システム ===')
print('概要: 動画から微細な明るさ変化を含むすべての変化を検出します')
print('操作方法:')
print('  q キー: プログラム終了')
print('  感度調整ウィンドウで検出感度を調整可能')
print('注意事項:')
print('  初期化に約30フレーム必要です')
print('  検出結果は自動的にresult.txtに保存されます')

print("\n0: 動画ファイル")
print("1: カメラ")
print("2: サンプル動画")

choice = input("選択: ")

# tkinterのメインウィンドウを作成
root = tk.Tk()
root.withdraw()

def stop_processing():
    global running
    running = False
    root.quit()

def request_gui_stop():
    try:
        root.after(0, root.quit)
    except:
        pass

if choice == '0':
    path = filedialog.askopenfilename()
    root.deiconify()
    if not path:
        exit()
    cap = cv2.VideoCapture(path)
    is_camera = False
elif choice == '1':
    root.deiconify()
    cap = cv2.VideoCapture(0, cv2.CAP_DSHOW)
    if not cap.isOpened():
        cap = cv2.VideoCapture(0)
    cap.set(cv2.CAP_PROP_BUFFERSIZE, 1)
    is_camera = True
else:
    root.deiconify()
    # サンプル動画ダウンロード・処理
    SAMPLE_URL = 'https://raw.githubusercontent.com/opencv/opencv/master/samples/data/vtest.avi'
    SAMPLE_FILE = 'vtest.avi'
    urllib.request.urlretrieve(SAMPLE_URL, SAMPLE_FILE)
    cap = cv2.VideoCapture(SAMPLE_FILE)
    is_camera = False

if not cap.isOpened():
    print('動画ファイル・カメラを開けませんでした')
    exit()

# 感度調整ウィンドウを作成
control_window = tk.Toplevel(root)
control_window.title('感度調整')
control_window.geometry('300x200')

label = Label(control_window, text='Sensitivity (感度)')
label.pack(pady=10)

scale = Scale(control_window, from_=0.5, to=3.0, resolution=0.1,
              orient='horizontal', length=250, command=update_sensitivity)
scale.set(sensitivity)
scale.pack(pady=10)

value_label = Label(control_window, text='')
value_label.pack()

stop_button = Button(control_window, text='終了', command=stop_processing)
stop_button.pack(pady=10)

def update_value_label():
    if running:
        value_label.config(text=f'現在の値: {sensitivity:.1f}')
        control_window.after(100, update_value_label)

update_value_label()

# 動画処理を別スレッドで開始
video_thread = threading.Thread(target=video_processing_thread, args=(cap,))
video_thread.daemon = True
video_thread.start()

print('\n=== 動画処理開始 ===')
print('操作方法:')
print('  q キー: プログラム終了')

# tkinterのメインループ
root.mainloop()

# スレッドの終了を待つ
video_thread.join()

print('\n=== プログラム終了 ===')
if results_log:
    with open('result.txt', 'w', encoding='utf-8') as f:
        f.write('=== 結果 ===\n')
        f.write(f'処理フレーム数: {frame_count}\n')
        f.write('\n')
        f.write('\n'.join(results_log))
    print(f'\n処理結果をresult.txtに保存しました')