ガイデッドフィルタ画像処理（入力：静止画）（ソースコードと説明と利用ガイド）

# ガイデッドフィルタ画像処理プログラム
#
# 特徴技術名: Guided Filter（ガイデッドフィルタ）
#
# 出典:
#   He, K., Sun, J., & Tang, X. (2013). Guided Image Filtering.
#   IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(6), 1397-1409.
#
# 特徴機能:
#   エッジ保存型スムージング
#   局所線形モデルに基づくエッジ保存フィルタ。ガイダンス画像の構造を保持しながら、入力画像をフィルタリングする。バイラテラルフィルタと比較して、勾配反転アーティファクトが発生しにくく、O(N)時間のアルゴリズムを持つ。
#
# 学習済みモデル: なし
#
# 特徴技術および学習済みモデルの利用制限:
#   アルゴリズムは論文に基づく自前実装。OpenCV（Apache 2ライセンス）、NumPy（BSDライセンス）の基本機能のみ使用。
#   必ず利用者自身でライセンス情報を確認すること。
#   参考: https://opencv.org/license/
#   参考: https://numpy.org/license.html
#
# 方式設計:
#   - 関連利用技術:
#     * OpenCV (cv2): 画像入出力、ボックスフィルタ（積分画像による局所平均計算）
#     * tkinter: ファイル選択ダイアログ
#     * NumPy: 画像データ型変換、行列演算
#     * Pillow: 日本語テキスト描画
#
#   - 入力と出力:
#     入力: 複数の静止画像，カメラ（ユーザは「0:画像ファイル，1:カメラ，2:サンプル画像」のメニューで選択．0:画像ファイルの場合はtkinterで複数ファイル選択可能．1の場合はOpenCVでカメラが開き，スペースキーで撮影（複数回可能）．2の場合はhttps://raw.githubusercontent.com/opencv/opencv/master/samples/data/fruits.jpg とhttps://raw.githubusercontent.com/opencv/opencv/master/samples/data/messi5.jpgとhttps://raw.githubusercontent.com/opencv/opencv/master/samples/data/aero3.jpgとhttps://upload.wikimedia.org/wikipedia/commons/3/3a/Cat03.jpgを使用）
#     出力: 処理結果画像を0001.png, 0002.png...として保存、OpenCV画面での表示、result.txtへの処理履歴保存（出力ファイル名、元ファイル名、パラメータ、処理日時）
#
#   - 処理手順:
#     1. 入力画像をuint8からfloat64（倍精度浮動小数点）へ変換し、0-1の範囲に正規化
#     2. Guided Filterのパラメータ設定
#        - radius: 局所ウィンドウの半径（デフォルト4）
#        - eps: 正則化パラメータ（デフォルト0.2^2）
#     3. 各局所ウィンドウで線形係数(a, b)を計算
#        - 式(7): a_k = cov(I,p) / (var(I) + eps)
#        - 式(8): b_k = mean(p) - a_k * mean(I)
#     4. 線形係数を平均化
#        - 式(10): q_i = mean(a) * I_i + mean(b)
#     5. 処理結果を0-255の範囲に戻し、元のデータ型（uint8）へ変換
#
#   - 前処理、後処理:
#     前処理: uint8からfloat64への型変換と0-1への正規化により、
#            数値的安定性を確保し精度を向上
#     後処理: float64からuint8への型変換により、画像表示と保存を可能にする
#
#   - 追加処理:
#     局所線形モデル: 各局所ウィンドウでq = aI + bの関係を仮定し、
#     入力pとの差を最小化する線形係数を計算する
#     3×3対称行列の逆行列を解析的に計算
#
#   - 調整を必要とする設定値:
#     * radius（局所ウィンドウ半径）: 現在4に設定
#       値が大きいほどスムージング効果が強くなるが、処理時間も増加
#       推奨値: 2～16（用途に応じて調整）
#     * eps（正則化パラメータ）: 現在0.04に設定
#       値が小さいほどエッジ保存性が高まる
#       推奨値: (0.1)^2～(0.4)^2（用途に応じて調整）
#
# 将来方策:
#   並列処理技術を用いることで、処理時間を短縮できる可能性がある。
#   例: マルチスレッド処理による各ウィンドウの並列計算
#
# その他の重要事項:
#   * Windows環境での動作を前提
#   * Python 3.10以上で動作
#   * opencv-contrib-pythonは不要（標準のopencv-pythonのみで動作）
#   * Guided Filterは以下の利点がある:
#     - O(N)時間のアルゴリズム
#     - 勾配反転アーティファクトの回避
#     - エッジ保存型スムージング
#     - 構造転送フィルタリング
#
# 前準備:
#   pip install opencv-python numpy pillow

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

# Guided Filterパラメータ
GUIDED_FILTER_RADIUS = 4
GUIDED_FILTER_EPS = 0.2 * 0.2

# フォント設定
FONT_PATH = 'C:/Windows/Fonts/meiryo.ttc'
FONT_SIZE_MAIN = 16
FONT_SIZE_SMALL = 12

try:
    font_main = ImageFont.truetype(FONT_PATH, FONT_SIZE_MAIN)
    font_small = ImageFont.truetype(FONT_PATH, FONT_SIZE_SMALL)
except:
    font_main = None
    font_small = None

results_log = []

def draw_texts_with_pillow(bgr_frame, texts):
    if font_main is None:
        return bgr_frame
    img_pil = Image.fromarray(cv2.cvtColor(bgr_frame, cv2.COLOR_BGR2RGB))
    draw = ImageDraw.Draw(img_pil)
    for item in texts:
        text = item['text']
        x, y = item['org']
        color = item['color']
        font_type = item.get('font_type', 'main')
        font = font_main if font_type == 'main' else font_small
        draw.text((x, y), text, font=font, fill=color)
    return cv2.cvtColor(np.array(img_pil), cv2.COLOR_RGB2BGR)

def solve_symmetric_3x3(a00, a01, a02, a11, a12, a22, b0, b1, b2):
    """
    対称3×3行列の連立方程式 A * x = b を解析的に解く

    入力：
    - 対称行列Aの上三角成分: a00, a01, a02, a11, a12, a22
      A = [[a00, a01, a02],
           [a01, a11, a12],
           [a02, a12, a22]]
    - 右辺ベクトルb: b0, b1, b2

    出力：
    - 解ベクトルx: (x0, x1, x2)

    手法：Cramerの公式による解析的計算
    """
    # 行列式の計算
    # det(A) = a00*(a11*a22 - a12*a12) - a01*(a01*a22 - a12*a02) + a02*(a01*a12 - a11*a02)
    det = (a00 * (a11 * a22 - a12 * a12) -
           a01 * (a01 * a22 - a12 * a02) +
           a02 * (a01 * a12 - a11 * a02))

    # 特異行列チェック
    if abs(det) < 1e-10:
        return 0.0, 0.0, 0.0

    inv_det = 1.0 / det

    # 余因子行列の計算（対称性を利用）
    # C = adj(A)^T
    C00 = a11 * a22 - a12 * a12
    C01 = a02 * a12 - a01 * a22
    C02 = a01 * a12 - a02 * a11
    C11 = a00 * a22 - a02 * a02
    C12 = a01 * a02 - a00 * a12
    C22 = a00 * a11 - a01 * a01

    # x = A^(-1) * b = (1/det(A)) * adj(A) * b
    x0 = inv_det * (C00 * b0 + C01 * b1 + C02 * b2)
    x1 = inv_det * (C01 * b0 + C11 * b1 + C12 * b2)
    x2 = inv_det * (C02 * b0 + C12 * b1 + C22 * b2)

    return x0, x1, x2

def guided_filter(image, guide, radius, eps):
    """
    Guided Filter実装（論文Algorithm 1に基づく）

    出典: He, K., Sun, J., & Tang, X. (2013). Guided Image Filtering.
          IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(6).

    image: H×W×3の入力画像（float64, 0-1範囲）
    guide: H×W×3のガイド画像（float64, 0-1範囲）
    radius: 局所ウィンドウ半径
    eps: 正則化パラメータ
    """
    if guide is None:
        guide = image

    if len(image.shape) == 2:
        image = image[:, :, np.newaxis]
        guide = guide[:, :, np.newaxis]

    h, w, c_p = image.shape
    _, _, c_I = guide.shape

    result = np.zeros((h, w, c_p), dtype=np.float64)

    for ch in range(c_p):
        p = image[:, :, ch]

        if c_I == 3:
            # カラーガイダンスの場合（式(19), (20), (21)）
            I_r = guide[:, :, 0]
            I_g = guide[:, :, 1]
            I_b = guide[:, :, 2]

            # 平均値の計算
            mean_I_r = cv2.boxFilter(I_r, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_I_g = cv2.boxFilter(I_g, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_I_b = cv2.boxFilter(I_b, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_p = cv2.boxFilter(p, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)

            # 共分散の計算
            mean_Ip_r = cv2.boxFilter(I_r * p, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_Ip_g = cv2.boxFilter(I_g * p, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_Ip_b = cv2.boxFilter(I_b * p, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)

            cov_Ip_r = mean_Ip_r - mean_I_r * mean_p
            cov_Ip_g = mean_Ip_g - mean_I_g * mean_p
            cov_Ip_b = mean_Ip_b - mean_I_b * mean_p

            # 共分散行列Σの計算
            var_I_rr = cv2.boxFilter(I_r * I_r, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT) - mean_I_r * mean_I_r
            var_I_rg = cv2.boxFilter(I_r * I_g, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT) - mean_I_r * mean_I_g
            var_I_rb = cv2.boxFilter(I_r * I_b, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT) - mean_I_r * mean_I_b
            var_I_gg = cv2.boxFilter(I_g * I_g, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT) - mean_I_g * mean_I_g
            var_I_gb = cv2.boxFilter(I_g * I_b, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT) - mean_I_g * mean_I_b
            var_I_bb = cv2.boxFilter(I_b * I_b, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT) - mean_I_b * mean_I_b

            # 線形係数aの計算（式(19)）
            a_r = np.zeros((h, w), dtype=np.float64)
            a_g = np.zeros((h, w), dtype=np.float64)
            a_b = np.zeros((h, w), dtype=np.float64)

            for i in range(h):
                for j in range(w):
                    # 対称行列の上三角成分（対角成分にepsを加算）
                    a00 = var_I_rr[i, j] + eps
                    a01 = var_I_rg[i, j]
                    a02 = var_I_rb[i, j]
                    a11 = var_I_gg[i, j] + eps
                    a12 = var_I_gb[i, j]
                    a22 = var_I_bb[i, j] + eps

                    # 右辺ベクトル
                    b0 = cov_Ip_r[i, j]
                    b1 = cov_Ip_g[i, j]
                    b2 = cov_Ip_b[i, j]

                    # 解析的に解く
                    x0, x1, x2 = solve_symmetric_3x3(a00, a01, a02, a11, a12, a22, b0, b1, b2)

                    a_r[i, j] = x0
                    a_g[i, j] = x1
                    a_b[i, j] = x2

            # 線形係数bの計算（式(20)）
            b = mean_p - a_r * mean_I_r - a_g * mean_I_g - a_b * mean_I_b

            # 係数の平均化
            mean_a_r = cv2.boxFilter(a_r, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_a_g = cv2.boxFilter(a_g, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_a_b = cv2.boxFilter(a_b, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_b = cv2.boxFilter(b, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)

            # 出力の計算（式(21)）
            result[:, :, ch] = mean_a_r * I_r + mean_a_g * I_g + mean_a_b * I_b + mean_b

        else:
            # グレースケールガイダンスの場合（式(7), (8), (10)）
            I = guide[:, :, 0]

            # Algorithm 1に従った実装
            mean_I = cv2.boxFilter(I, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_p = cv2.boxFilter(p, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)

            corr_I = cv2.boxFilter(I * I, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            corr_Ip = cv2.boxFilter(I * p, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)

            var_I = corr_I - mean_I * mean_I
            cov_Ip = corr_Ip - mean_I * mean_p

            # 式(7), (8)
            a = cov_Ip / (var_I + eps)
            b = mean_p - a * mean_I

            # 式(10)
            mean_a = cv2.boxFilter(a, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)
            mean_b = cv2.boxFilter(b, cv2.CV_64F, (2*radius+1, 2*radius+1), normalize=True, borderType=cv2.BORDER_REFLECT)

            result[:, :, ch] = mean_a * I + mean_b

    return result

def image_processing(bgr_img):
    img_float64 = bgr_img.astype('float64') / 255.0

    processed_float64 = guided_filter(
        image=img_float64,
        guide=img_float64,
        radius=GUIDED_FILTER_RADIUS,
        eps=GUIDED_FILTER_EPS
    )

    processed_img = np.clip(processed_float64 * 255.0, 0, 255).astype(bgr_img.dtype)

    texts = [
        {'text': f'Radius: {GUIDED_FILTER_RADIUS}', 'org': (10, 30), 'color': (0, 255, 0), 'font_type': 'main'},
        {'text': f'Eps: {GUIDED_FILTER_EPS:.4f}', 'org': (10, 60), 'color': (0, 255, 0), 'font_type': 'main'}
    ]
    processed_img_display = draw_texts_with_pillow(processed_img.copy(), texts)

    return processed_img, processed_img_display, time.time()

def process_and_display_images(image_sources, source_type):
    idx = 1
    results = []

    for source in image_sources:
        img = cv2.imread(source) if source_type == 'file' else source
        if img is None:
            continue

        processed_img, processed_img_display, current_time = image_processing(img)

        output_filename = f"{idx:04d}.png"
        cv2.imwrite(output_filename, processed_img)

        if source_type == 'file':
            source_name = os.path.basename(source)
        else:
            source_name = f"Camera_{idx}"

        results.append({
            'original': img,
            'processed': processed_img_display,
            'output_file': output_filename,
            'source_name': source_name,
            'time': current_time,
            'index': idx
        })

        timestamp = datetime.fromtimestamp(current_time).strftime("%Y-%m-%d %H:%M:%S.%f")[:-3]
        print(f"{timestamp} {output_filename} を保存しました")
        results_log.append({
            'index': idx,
            'output_file': output_filename,
            'source_name': source_name,
            'radius': GUIDED_FILTER_RADIUS,
            'eps': GUIDED_FILTER_EPS,
            'timestamp': timestamp
        })
        idx += 1

    for item in results:
        cv2.imshow(f'Image_{item["index"]}', item['original'])
        cv2.imshow(f'Guided_Filter_{item["index"]}', item['processed'])

print("=" * 60)
print("Guided Filter画像処理プログラム（自前実装版）")
print("=" * 60)
print("\n操作方法:")
print("0: 画像ファイルを選択")
print("1: カメラを使用（スペースキーで処理、qキーで終了）")
print("2: サンプル画像を使用")
print("-" * 60)
choice = input("選択してください (0/1/2): ")

try:
    if choice == '0':
        root = tk.Tk()
        root.withdraw()
        if not (paths := filedialog.askopenfilenames(
            title="画像ファイルを選択",
            filetypes=[("画像ファイル", "*.jpg *.jpeg *.png *.bmp"), ("すべてのファイル", "*.*")]
        )):
            print("ファイルが選択されませんでした。")
            exit()
        print(f"\n{len(paths)}個のファイルを選択しました。処理を開始します...\n")
        process_and_display_images(paths, 'file')
        print("\n画像を表示中です。任意のキーを押すと終了します。")
        cv2.waitKey(0)

    elif choice == '1':
        print("\nカメラを起動しています...")
        cap = cv2.VideoCapture(0, cv2.CAP_DSHOW)
        if not cap.isOpened():
            cap = cv2.VideoCapture(0)
        cap.set(cv2.CAP_PROP_BUFFERSIZE, 1)
        print("カメラが起動しました。")
        print("  - スペースキー: 画像を処理")
        print("  - qキー: 終了\n")

        camera_idx = 1
        try:
            while True:
                ret, frame = cap.read()
                if not ret:
                    break
                cv2.imshow('Camera', frame)
                key = cv2.waitKey(1) & 0xFF
                if key == ord(' '):
                    processed_img, processed_img_display, current_time = image_processing(frame)

                    output_filename = f"{camera_idx:04d}.png"
                    cv2.imwrite(output_filename, processed_img)

                    timestamp = datetime.fromtimestamp(current_time).strftime("%Y-%m-%d %H:%M:%S.%f")[:-3]
                    print(f"{timestamp} {output_filename} を保存しました")

                    results_log.append({
                        'index': camera_idx,
                        'output_file': output_filename,
                        'source_name': f"Camera_{camera_idx}",
                        'radius': GUIDED_FILTER_RADIUS,
                        'eps': GUIDED_FILTER_EPS,
                        'timestamp': timestamp
                    })

                    cv2.imshow('Guided_Filter', processed_img_display)
                    camera_idx += 1
                elif key == ord('q'):
                    break
        finally:
            cap.release()

    else:
        print("\nサンプル画像をダウンロードしています...")
        opener = urllib.request.build_opener()
        opener.addheaders = [('User-Agent', 'Mozilla/5.0')]
        urllib.request.install_opener(opener)

        urls = [
            "https://raw.githubusercontent.com/opencv/opencv/master/samples/data/fruits.jpg?raw=true",
            "https://raw.githubusercontent.com/opencv/opencv/master/samples/data/messi5.jpg?raw=true",
            "https://raw.githubusercontent.com/opencv/opencv/master/samples/data/aero3.jpg?raw=true",
            "https://upload.wikimedia.org/wikipedia/commons/3/3a/Cat03.jpg"
        ]
        files = []
        for i, url in enumerate(urls):
            try:
                urllib.request.urlretrieve(url, f"sample_{i}.jpg")
                files.append(f"sample_{i}.jpg")
                print(f"sample_{i}.jpg をダウンロードしました")
            except Exception as e:
                print(f"画像のダウンロードに失敗しました: {url}")
                print(f"エラー: {e}")

        if files:
            print(f"\n{len(files)}個のサンプル画像の処理を開始します...\n")
            process_and_display_images(files, 'file')
            print("\n画像を表示中です。任意のキーを押すと終了します。")
            cv2.waitKey(0)
        else:
            print("\nサンプル画像のダウンロードに失敗しました。")

finally:
    print('\n' + "=" * 60)
    print('プログラムを終了します')
    print("=" * 60)
    cv2.destroyAllWindows()
    if results_log:
        with open('result.txt', 'w', encoding='utf-8') as f:
            f.write('=' * 80 + '\n')
            f.write('Guided Filter 処理結果\n')
            f.write('=' * 80 + '\n\n')

            for item in results_log:
                f.write(f"No.{item['index']:04d}\n")
                f.write(f"  出力ファイル: {item['output_file']}\n")
                f.write(f"  元ファイル名: {item['source_name']}\n")
                f.write(f"  パラメータ  : radius={item['radius']}, eps={item['eps']}\n")
                f.write(f"  処理日時    : {item['timestamp']}\n")
                f.write('-' * 80 + '\n')

            f.write('\n' + '=' * 80 + '\n')
            f.write(f'合計 {len(results_log)} 件の画像を処理しました\n')
            f.write('=' * 80 + '\n')
        print(f'\n処理結果をresult.txtに保存しました（{len(results_log)}件）')