"""
Enhanced Subpixel Analysis for Ancient Structure Measurement
基本単位検出のための高精度サブピクセル解析（改良版）

Author: Tatsuo Kou
License: MIT
Version: 2.0.0 (Enhanced)
Date: 2026-03-17

Citation: 
Kou, T. (2026). Empirical Verification of a Prehistoric Global Survey System: 
Statistical Proof of the 8.314m Basic Unit. Scientific Reports.

Enhancements in v2.0:
- Statistical significance testing (chi-square, KS test)
- Measurement uncertainty propagation
- GeoTIFF support with automatic resolution detection
- Batch processing capability
- Enhanced visualization with confidence intervals
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import ndimage, optimize, stats
from PIL import Image
import pandas as pd
from pathlib import Path
import warnings

# Optional: GeoTIFF support
try:
    from osgeo import gdal
    GDAL_AVAILABLE = True
except ImportError:
    GDAL_AVAILABLE = False
    warnings.warn("GDAL not available. GeoTIFF support disabled. Use: conda install gdal")


# ================================
# 1. 画像読み込みと前処理
# ================================

def load_satellite_image(filepath, pixel_resolution=None):
    """
    衛星画像を読み込む（GeoTIFF自動対応）
    
    Parameters:
    -----------
    filepath : str
        画像ファイルパス
    pixel_resolution : float, optional
        ピクセル解像度（m/pixel）
        Noneの場合、GeoTIFFから自動取得
    
    Returns:
    --------
    image : ndarray
        画像データ
    resolution : float
        解像度（m/pixel）
    metadata : dict
        メタデータ
    """
    filepath = Path(filepath)
    
    # GeoTIFFの場合
    if filepath.suffix.lower() in ['.tif', '.tiff'] and GDAL_AVAILABLE:
        return load_geotiff_with_metadata(str(filepath))
    
    # 通常の画像ファイル
    img = Image.open(filepath)
    image = np.array(img.convert('L'))  # グレースケール変換
    
    if pixel_resolution is None:
        pixel_resolution = 0.5  # デフォルト: ALOS-2解像度
        warnings.warn(f"Resolution not specified. Using default: {pixel_resolution}m/pixel")
    
    metadata = {
        'filename': filepath.name,
        'size': image.shape,
        'resolution_m': pixel_resolution
    }
    
    print(f"画像サイズ: {image.shape}")
    print(f"解像度: {pixel_resolution}m/pixel")
    
    return image, pixel_resolution, metadata


def load_geotiff_with_metadata(filepath):
    """
    GeoTIFFを読み込み、地理参照情報も取得
    
    Parameters:
    -----------
    filepath : str
        GeoTIFFファイルパス
    
    Returns:
    --------
    image : ndarray
        画像データ
    resolution : float
        ピクセル解像度（m/pixel）
    metadata : dict
        地理参照情報
    """
    if not GDAL_AVAILABLE:
        raise ImportError("GDAL required for GeoTIFF. Install: conda install gdal")
    
    dataset = gdal.Open(filepath)
    if dataset is None:
        raise FileNotFoundError(f"Cannot open: {filepath}")
    
    # 地理変換パラメータ取得
    geotransform = dataset.GetGeoTransform()
    pixel_width = abs(geotransform[1])   # X方向解像度
    pixel_height = abs(geotransform[5])  # Y方向解像度
    
    # 画像データ読み込み
    band = dataset.GetRasterBand(1)
    image = band.ReadAsArray()
    
    # メタデータ
    metadata = {
        'filename': Path(filepath).name,
        'size': image.shape,
        'resolution_m': pixel_width,
        'resolution_x': pixel_width,
        'resolution_y': pixel_height,
        'geotransform': geotransform,
        'projection': dataset.GetProjection(),
        'origin_x': geotransform[0],
        'origin_y': geotransform[3]
    }
    
    print(f"GeoTIFF loaded: {metadata['filename']}")
    print(f"  Size: {image.shape}")
    print(f"  Resolution: {pixel_width:.3f}m/pixel")
    print(f"  Projection: {metadata['projection'][:50]}...")
    
    return image, pixel_width, metadata


def preprocess_image(image, sigma=1.0):
    """
    画像の前処理（ノイズ除去）
    
    Parameters:
    -----------
    image : ndarray
        入力画像
    sigma : float
        ガウシアンフィルタの標準偏差
    
    Returns:
    --------
    processed : ndarray
        処理済み画像
    """
    # ガウシアンフィルタでノイズ除去
    processed = ndimage.gaussian_filter(image, sigma=sigma)
    
    return processed


# ================================
# 2. エッジ検出（サブピクセル精度）
# ================================

def detect_edges_subpixel(image, threshold=0.1):
    """
    サブピクセル精度でエッジを検出
    
    Parameters:
    -----------
    image : ndarray
        入力画像
    threshold : float
        エッジ検出閾値（0-1、勾配最大値に対する相対値）
    
    Returns:
    --------
    edges : list of tuples
        エッジ座標 [(x1, y1), (x2, y2), ...]
        サブピクセル精度（小数点以下含む）
    """
    # 1次微分（Sobel）
    gradient_x = ndimage.sobel(image, axis=1)
    gradient_y = ndimage.sobel(image, axis=0)
    
    # 勾配の大きさ
    gradient_magnitude = np.sqrt(gradient_x**2 + gradient_y**2)
    
    # 閾値以上の点をエッジとして検出
    edge_points = np.argwhere(gradient_magnitude > threshold * gradient_magnitude.max())
    
    # サブピクセル精度でエッジ位置を補正
    edges_subpixel = []
    
    for y, x in edge_points:
        # 3x3近傍の勾配値を取得
        if y > 0 and y < image.shape[0]-1 and x > 0 and x < image.shape[1]-1:
            neighborhood = gradient_magnitude[y-1:y+2, x-1:x+2]
            
            # 重心計算によるサブピクセル補正
            y_offset, x_offset = ndimage.center_of_mass(neighborhood)
            
            # 補正後の座標
            x_subpixel = x + (x_offset - 1)  # -1は中心からのオフセット
            y_subpixel = y + (y_offset - 1)
            
            edges_subpixel.append((x_subpixel, y_subpixel))
    
    return edges_subpixel


# ================================
# 3. 構造物の測定
# ================================

def measure_structure(edges, resolution, roi=None):
    """
    検出されたエッジから構造物のサイズを測定
    
    Parameters:
    -----------
    edges : list of tuples
        エッジ座標
    resolution : float
        ピクセル解像度（m/pixel）
    roi : tuple, optional
        測定領域 (x_min, x_max, y_min, y_max)
    
    Returns:
    --------
    measurements : dict
        測定結果
    """
    edges_array = np.array(edges)
    
    # ROIでフィルタ
    if roi is not None:
        x_min, x_max, y_min, y_max = roi
        mask = ((edges_array[:, 0] >= x_min) & 
                (edges_array[:, 0] <= x_max) &
                (edges_array[:, 1] >= y_min) & 
                (edges_array[:, 1] <= y_max))
        edges_array = edges_array[mask]
    
    if len(edges_array) == 0:
        return None
    
    # X方向の幅
    width_pixels = edges_array[:, 0].max() - edges_array[:, 0].min()
    width_meters = width_pixels * resolution
    
    # Y方向の高さ
    height_pixels = edges_array[:, 1].max() - edges_array[:, 1].min()
    height_meters = height_pixels * resolution
    
    measurements = {
        'width_m': width_meters,
        'height_m': height_meters,
        'width_pixels': width_pixels,
        'height_pixels': height_pixels,
        'num_edges': len(edges_array),
        'roi': roi
    }
    
    return measurements


# ================================
# 4. 測定誤差の計算（新機能）
# ================================

def calculate_measurement_uncertainty(pixel_error=0.1, resolution=0.5, 
                                     geometric_error=0.01):
    """
    測定誤差の伝播計算
    
    Parameters:
    -----------
    pixel_error : float
        サブピクセル精度（pixels）
        デフォルト: 0.1 pixels（重心法の典型値）
    resolution : float
        ピクセル解像度（m/pixel）
    geometric_error : float
        画像幾何歪み誤差（m）
        ALOS-2仕様: ~0.01m
    
    Returns:
    --------
    uncertainty : dict
        誤差情報
    """
    # サブピクセル誤差
    subpixel_error_m = pixel_error * resolution
    
    # 総合誤差（二乗和の平方根）
    total_error_m = np.sqrt(subpixel_error_m**2 + geometric_error**2)
    
    uncertainty = {
        'subpixel_error_m': subpixel_error_m,
        'geometric_error_m': geometric_error,
        'total_error_m': total_error_m,
        'relative_error_percent': (total_error_m / 8.314) * 100  # 8.314mに対する相対誤差
    }
    
    return uncertainty


# ================================
# 5. 基本単位の検出と統計検定（強化版）
# ================================

def detect_fundamental_unit(measurements_list, expected_unit=8.314, 
                           tolerance=0.1, confidence_level=0.95):
    """
    複数の測定値から基本単位を検出し、統計的有意性を評価
    
    Parameters:
    -----------
    measurements_list : list of dict
        測定値のリスト
    expected_unit : float
        期待される基本単位（m）
    tolerance : float
        許容誤差（m）
    confidence_level : float
        信頼水準（デフォルト: 95%）
    
    Returns:
    --------
    results : dict
        解析結果と統計検定結果
    """
    # 全測定値を抽出
    widths = [m['width_m'] for m in measurements_list if m is not None]
    heights = [m['height_m'] for m in measurements_list if m is not None]
    all_values = widths + heights
    
    if len(all_values) == 0:
        return None
    
    # 基本単位の整数倍を検出
    ratios = []
    matched_values = []
    
    for value in all_values:
        ratio = value / expected_unit
        # 整数に近いか判定
        nearest_int = round(ratio)
        error = abs(ratio - nearest_int)
        
        if error < tolerance / expected_unit:
            ratios.append({
                'value': value,
                'ratio': ratio,
                'nearest_int': nearest_int,
                'error_m': error * expected_unit
            })
            matched_values.append(value)
    
    # 基本単位のフィッティング
    if ratios:
        nearest_ints = np.array([r['nearest_int'] for r in ratios])
        values = np.array([r['value'] for r in ratios])
        
        # 重み付き最小二乗法
        fitted_unit = np.sum(values * nearest_ints) / np.sum(nearest_ints**2)
        
        # 標準偏差
        residuals = values - fitted_unit * nearest_ints
        std_unit = np.std(residuals / nearest_ints)
        
        # 信頼区間
        n = len(values)
        se = std_unit / np.sqrt(n)
        t_value = stats.t.ppf((1 + confidence_level) / 2, n - 1)
        ci_lower = fitted_unit - t_value * se
        ci_upper = fitted_unit + t_value * se
        
    else:
        fitted_unit = None
        std_unit = None
        ci_lower = None
        ci_upper = None
    
    # 統計検定
    statistical_tests = perform_statistical_tests(
        all_values, matched_values, expected_unit
    )
    
    results = {
        'expected_unit': expected_unit,
        'fitted_unit': fitted_unit,
        'std_unit': std_unit,
        'confidence_interval': (ci_lower, ci_upper),
        'confidence_level': confidence_level,
        'num_matches': len(ratios),
        'total_measurements': len(all_values),
        'match_rate': len(ratios) / len(all_values) if all_values else 0,
        'details': ratios,
        'statistical_tests': statistical_tests
    }
    
    return results


def perform_statistical_tests(all_values, matched_values, expected_unit=8.314):
    """
    統計的有意性検定（新機能）
    
    Parameters:
    -----------
    all_values : list
        全測定値
    matched_values : list
        基本単位にマッチした値
    expected_unit : float
        期待される基本単位
    
    Returns:
    --------
    tests : dict
        各種統計検定結果
    """
    tests = {}
    
    # 1. カイ二乗適合度検定
    if len(matched_values) > 0:
        ratios = np.array(matched_values) / expected_unit
        nearest_ints = np.round(ratios)
        
        # 観測度数
        unique_ints, observed_counts = np.unique(nearest_ints, return_counts=True)
        
        # 期待度数（一様分布を仮定）
        expected_counts = np.ones_like(observed_counts) * len(matched_values) / len(unique_ints)
        
        # カイ二乗検定
        chi2_stat, chi2_p = stats.chisquare(observed_counts, expected_counts)
        
        tests['chi_square'] = {
            'statistic': chi2_stat,
            'p_value': chi2_p,
            'significant': chi2_p < 0.001,
            'interpretation': 'Matches expected unit' if chi2_p < 0.001 else 'Random distribution'
        }
    
    # 2. Kolmogorov-Smirnov検定（連続性検定）
    if len(all_values) > 0:
        # 基本単位の整数倍からの距離を計算
        ratios = np.array(all_values) / expected_unit
        residuals = ratios - np.round(ratios)
        
        # 一様分布と比較（H0: ランダム）
        ks_stat, ks_p = stats.kstest(residuals, 'uniform', args=(-0.5, 1.0))
        
        tests['kolmogorov_smirnov'] = {
            'statistic': ks_stat,
            'p_value': ks_p,
            'significant': ks_p < 0.001,
            'interpretation': 'Non-random clustering' if ks_p < 0.001 else 'Random distribution'
        }
    
    # 3. Z検定（マッチ率の有意性）
    n_total = len(all_values)
    n_matched = len(matched_values)
    
    if n_total > 0:
        # 帰無仮説: マッチ率 = 10%（偶然の一致）
        p0 = 0.1
        p_observed = n_matched / n_total
        
        # Z統計量
        z_stat = (p_observed - p0) / np.sqrt(p0 * (1 - p0) / n_total)
        z_p = 2 * (1 - stats.norm.cdf(abs(z_stat)))  # 両側検定
        
        tests['z_test'] = {
            'match_rate': p_observed,
            'null_hypothesis': p0,
            'z_statistic': z_stat,
            'p_value': z_p,
            'significant': z_p < 0.001,
            'interpretation': 'Match rate significantly higher than chance' if z_p < 0.001 else 'Chance level'
        }
    
    return tests


# ================================
# 6. バッチ処理（新機能）
# ================================

def batch_analysis(file_list, expected_unit=8.314, output_dir='results'):
    """
    複数画像の一括解析
    
    Parameters:
    -----------
    file_list : list of str
        画像ファイルパスのリスト
    expected_unit : float
        期待される基本単位
    output_dir : str
        結果出力ディレクトリ
    
    Returns:
    --------
    results_df : pandas.DataFrame
        全結果の集計
    """
    output_path = Path(output_dir)
    output_path.mkdir(exist_ok=True)
    
    all_results = []
    
    for i, filepath in enumerate(file_list):
        print(f"\n{'='*60}")
        print(f"Processing {i+1}/{len(file_list)}: {Path(filepath).name}")
        print('='*60)
        
        try:
            # 解析実行
            image, resolution, metadata = load_satellite_image(filepath)
            processed = preprocess_image(image, sigma=1.0)
            edges = detect_edges_subpixel(processed, threshold=0.1)
            
            # 複数ROIで測定（グリッド分割）
            measurements = auto_detect_structures(edges, resolution, image.shape)
            
            # 基本単位検出
            result = detect_fundamental_unit(measurements, expected_unit=expected_unit)
            
            if result:
                result['filename'] = Path(filepath).name
                result['resolution'] = resolution
                all_results.append(result)
                
                # 可視化保存
                plot_filename = output_path / f"{Path(filepath).stem}_analysis.png"
                plot_results(image, edges, measurements, result, save_path=plot_filename)
            
        except Exception as e:
            print(f"Error processing {filepath}: {e}")
            continue
    
    # DataFrameに変換
    if all_results:
        results_df = pd.DataFrame([
            {
                'filename': r['filename'],
                'fitted_unit': r['fitted_unit'],
                'std_unit': r['std_unit'],
                'match_rate': r['match_rate'],
                'chi2_p': r['statistical_tests'].get('chi_square', {}).get('p_value', None),
                'ks_p': r['statistical_tests'].get('kolmogorov_smirnov', {}).get('p_value', None)
            }
            for r in all_results
        ])
        
        # CSV保存
        results_df.to_csv(output_path / 'batch_results.csv', index=False)
        
        return results_df
    else:
        return None


def auto_detect_structures(edges, resolution, image_shape, grid_size=5):
    """
    画像を自動的にグリッド分割して構造物を検出
    
    Parameters:
    -----------
    edges : list
        エッジ座標
    resolution : float
        解像度
    image_shape : tuple
        画像サイズ (height, width)
    grid_size : int
        グリッド分割数
    
    Returns:
    --------
    measurements : list
        測定結果のリスト
    """
    height, width = image_shape
    x_step = width / grid_size
    y_step = height / grid_size
    
    measurements = []
    
    for i in range(grid_size):
        for j in range(grid_size):
            x_min = i * x_step
            x_max = (i + 1) * x_step
            y_min = j * y_step
            y_max = (j + 1) * y_step
            
            roi = (x_min, x_max, y_min, y_max)
            m = measure_structure(edges, resolution, roi=roi)
            
            if m and m['num_edges'] > 10:  # 最小エッジ数閾値
                measurements.append(m)
    
    return measurements


# ================================
# 7. 可視化（強化版）
# ================================

def plot_results(image, edges, measurements, results, save_path=None):
    """
    解析結果の可視化（信頼区間付き）
    """
    fig, axes = plt.subplots(2, 2, figsize=(16, 14))
    
    # (1) 元画像
    axes[0, 0].imshow(image, cmap='gray')
    axes[0, 0].set_title('Original Image', fontsize=14, fontweight='bold')
    axes[0, 0].axis('off')
    
    # (2) エッジ検出結果
    edges_array = np.array(edges)
    axes[0, 1].imshow(image, cmap='gray', alpha=0.5)
    axes[0, 1].scatter(edges_array[:, 0], edges_array[:, 1], 
                      c='red', s=0.5, alpha=0.3, label='Detected edges')
    axes[0, 1].set_title(f'Subpixel Edge Detection ({len(edges)} points)', 
                        fontsize=14, fontweight='bold')
    axes[0, 1].legend()
    axes[0, 1].axis('off')
    
    # (3) 測定値のヒストグラム
    all_values = ([m['width_m'] for m in measurements if m] + 
                  [m['height_m'] for m in measurements if m])
    
    axes[1, 0].hist(all_values, bins=50, alpha=0.7, edgecolor='black', color='steelblue')
    
    # 基本単位の倍数にマーカー
    if results and results['fitted_unit']:
        for i in range(1, 15):
            axes[1, 0].axvline(results['fitted_unit'] * i, 
                             color='red', linestyle='--', alpha=0.5, linewidth=1)
    
    axes[1, 0].set_xlabel('Size (m)', fontsize=12)
    axes[1, 0].set_ylabel('Frequency', fontsize=12)
    axes[1, 0].grid(True, alpha=0.3)
    
    # タイトルに統計情報
    if results and results['fitted_unit']:
        ci_lower, ci_upper = results['confidence_interval']
        title = f"Fitted Unit: {results['fitted_unit']:.3f}m ± {results['std_unit']:.3f}m\n"
        title += f"95% CI: [{ci_lower:.3f}, {ci_upper:.3f}]m\n"
        title += f"Match Rate: {results['match_rate']*100:.1f}%"
        axes[1, 0].set_title(title, fontsize=12, fontweight='bold')
    
    # (4) 統計検定結果
    axes[1, 1].axis('off')
    
    if results and results['statistical_tests']:
        stats_text = "Statistical Tests:\n\n"
        
        # カイ二乗検定
        if 'chi_square' in results['statistical_tests']:
            chi2 = results['statistical_tests']['chi_square']
            stats_text += f"Chi-Square Test:\n"
            stats_text += f"  χ² = {chi2['statistic']:.3f}\n"
            stats_text += f"  p-value = {chi2['p_value']:.2e}\n"
            stats_text += f"  Result: {chi2['interpretation']}\n\n"
        
        # KS検定
        if 'kolmogorov_smirnov' in results['statistical_tests']:
            ks = results['statistical_tests']['kolmogorov_smirnov']
            stats_text += f"Kolmogorov-Smirnov Test:\n"
            stats_text += f"  D = {ks['statistic']:.3f}\n"
            stats_text += f"  p-value = {ks['p_value']:.2e}\n"
            stats_text += f"  Result: {ks['interpretation']}\n\n"
        
        # Z検定
        if 'z_test' in results['statistical_tests']:
            z = results['statistical_tests']['z_test']
            stats_text += f"Z-Test (Match Rate):\n"
            stats_text += f"  Observed: {z['match_rate']*100:.1f}%\n"
            stats_text += f"  Expected (H₀): {z['null_hypothesis']*100:.1f}%\n"
            stats_text += f"  Z = {z['z_statistic']:.3f}\n"
            stats_text += f"  p-value = {z['p_value']:.2e}\n"
            stats_text += f"  Result: {z['interpretation']}\n"
        
        axes[1, 1].text(0.1, 0.9, stats_text, 
                       transform=axes[1, 1].transAxes,
                       fontsize=10, verticalalignment='top',
                       fontfamily='monospace',
                       bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.3))
    
    plt.tight_layout()
    
    if save_path:
        plt.savefig(save_path, dpi=300, bbox_inches='tight')
        print(f"Saved: {save_path}")
    
    plt.show()


# ================================
# 8. メイン実行例
# ================================

def main_example():
    """
    使用例（デモ）
    """
    print("=" * 70)
    print("Enhanced Subpixel Analysis - Example with Statistical Tests")
    print("=" * 70)
    
    # デモ用: ダミー画像生成
    print("\n[Demo] Creating dummy structure image...")
    image = create_dummy_structure_image(size=600, unit_pixels=16.628, noise_level=15)
    resolution = 0.5  # m/pixel
    
    # 測定誤差計算
    print("\n[0] Calculating measurement uncertainty...")
    uncertainty = calculate_measurement_uncertainty(
        pixel_error=0.1,
        resolution=resolution,
        geometric_error=0.01
    )
    print(f"   Subpixel error: ±{uncertainty['subpixel_error_m']:.3f}m")
    print(f"   Geometric error: ±{uncertainty['geometric_error_m']:.3f}m")
    print(f"   Total error: ±{uncertainty['total_error_m']:.3f}m")
    print(f"   Relative error: {uncertainty['relative_error_percent']:.2f}%")
    
    # 前処理
    print("\n[1] Preprocessing...")
    processed = preprocess_image(image, sigma=1.0)
    
    # エッジ検出
    print("\n[2] Detecting edges (subpixel precision)...")
    edges = detect_edges_subpixel(processed, threshold=0.1)
    print(f"   Detected {len(edges)} edge points")
    
    # 自動構造検出
    print("\n[3] Auto-detecting structures...")
    measurements = auto_detect_structures(edges, resolution, image.shape, grid_size=5)
    print(f"   Found {len(measurements)} structures")
    
    # 基本単位検出と統計検定
    print("\n[4] Detecting fundamental unit with statistical tests...")
    results = detect_fundamental_unit(measurements, expected_unit=8.314, tolerance=0.15)
    
    if results:
        print(f"\n   Expected unit: {results['expected_unit']:.3f}m")
        print(f"   Fitted unit: {results['fitted_unit']:.3f}m ± {results['std_unit']:.3f}m")
        ci_lower, ci_upper = results['confidence_interval']
        print(f"   95% CI: [{ci_lower:.3f}, {ci_upper:.3f}]m")
        print(f"   Match rate: {results['match_rate']*100:.1f}%")
        
        # 統計検定結果
        if 'chi_square' in results['statistical_tests']:
            chi2 = results['statistical_tests']['chi_square']
            print(f"\n   Chi-square test: p = {chi2['p_value']:.2e}")
        
        if 'kolmogorov_smirnov' in results['statistical_tests']:
            ks = results['statistical_tests']['kolmogorov_smirnov']
            print(f"   KS test: p = {ks['p_value']:.2e}")
        
        if 'z_test' in results['statistical_tests']:
            z = results['statistical_tests']['z_test']
            print(f"   Z test: p = {z['p_value']:.2e}")
    
    # 可視化
    print("\n[5] Plotting results...")
    plot_results(image, edges, measurements, results)
    
    print("\n" + "=" * 70)
    print("Analysis complete!")
    print("=" * 70)
    
    return results


def create_dummy_structure_image(size=600, unit_pixels=16.628, noise_level=15):
    """
    デモ用: 8.314mの倍数で配置された構造物のダミー画像生成
    
    Parameters:
    -----------
    size : int
        画像サイズ（pixel）
    unit_pixels : float
        基本単位のピクセル数（8.314m / 0.5m/pixel = 16.628 pixels）
    noise_level : float
        ノイズレベル
    """
    image = np.zeros((size, size))
    
    # 複数の矩形構造を配置（8.314mの整数倍 + ランダム配置）
    np.random.seed(42)
    
    structures = [
        (50, 50, 1),      # 1倍
        (150, 50, 2),     # 2倍
        (300, 50, 3),     # 3倍
        (50, 200, 4),     # 4倍
        (250, 200, 2),    # 2倍
        (450, 200, 1),    # 1倍
        (100, 400, 3),    # 3倍
        (350, 400, 2),    # 2倍
    ]
    
    for x, y, multiplier in structures:
        w = int(unit_pixels * multiplier)
        h = int(unit_pixels * multiplier * 0.8)  # 少し縦長
        
        # 矩形を描画
        x1, x2 = int(x), min(int(x + w), size)
        y1, y2 = int(y), min(int(y + h), size)
        
        if x2 < size and y2 < size:
            image[y1:y2, x1:x2] = 180  # 内部
            # エッジ強調
            edge_width = 2
            image[y1:y1+edge_width, x1:x2] = 255
            image[y2-edge_width:y2, x1:x2] = 255
            image[y1:y2, x1:x1+edge_width] = 255
            image[y1:y2, x2-edge_width:x2] = 255
    
    # ノイズ追加
    noise = np.random.normal(0, noise_level, image.shape)
    image = image + noise
    image = np.clip(image, 0, 255)
    
    return image


# ================================
# 実行
# ================================

if __name__ == "__main__":
    # シングル解析
    results = main_example()
    
    # バッチ処理の例（コメントアウト）
    # file_list = ['image1.tif', 'image2.tif', 'image3.tif']
    # batch_results = batch_analysis(file_list, expected_unit=8.314)
    # print(batch_results)