国土地理院DEM PNG to OBJ変換ツール（ソースコードと実行結果）

【概要】

国土地理院が配信する標高タイル（PNG形式）のファイルを読み込み、RGB値から標高値を復元して、3DメッシュのOBJファイルに変換するプログラムを扱う。標高値の分布をカラーマップで確認したうえで、Delaunay三角形分割で三角形メッシュを作り、Quadric Error Metrics（QEM）によりメッシュを簡略化して出力する。簡略化では、勾配・曲率・粗さから求めた重要度を使い、地形の特徴が残るようにする。

【目次】

【関連する外部ページ】

【サイト内の関連情報】

第1章 Python開発環境，ライブラリ類

ここでは、最低限の事前準備について説明する。機械学習や深層学習を行う場合は、NVIDIA CUDA、Visual Studio、Cursorなどを追加でインストールすると便利である。これらについては別ページ https://www.kkaneko.jp/cc/dev/aiassist.html で解説しているので、必要に応じて参照すること。

第2章 Python 3.12 のインストール

Pythonのインストールを行い、Pythonのプログラムを実行する環境を整える。扱う環境は、Windows搭載パソコンである。金子研究室では、Python 3.12.10を推奨する。

[Windows での Python 3.12 のインストール手順を見るには、ここをクリック]

Windows での Python 3.12 のインストール

以下のいずれかの方法でPython 3.12をインストールする。Pythonがインストール済みの場合、この手順は不要である。

方法 1：winget によるインストール

【インストールコマンドの実行方法】

管理者権限でコマンドプロンプトを起動する（手順：Windowsキーまたはスタートメニュー → cmd と入力 → 右クリック → 「管理者として実行」）。そして、コマンド全体をコマンドプロンプトにコピー＆ペーストする。

--scope machine を指定することで、システム全体（全ユーザー向け）にインストールされる。このオプションの実行には管理者権限が必要である。インストール完了後、コマンドプロンプトを再起動するとPATHが反映される。

REM Python 3.12 をシステム領域にインストール
winget install --id Python.Python.3.12 -e --scope machine --silent --accept-source-agreements --accept-package-agreements --override "/quiet InstallAllUsers=1 PrependPath=1 Include_test=0 Include_pip=1 Include_launcher=1 InstallLauncherAllUsers=1 TargetDir=\"C:\Program Files\Python312\""
if not "%ERRORLEVEL%"=="0" ( color 0c & echo Python 3.12 のインストールに失敗しました & ping 127.0.0.1 -n 6 >nul & color )

REM Python と Scripts を PATH 先頭に追加
powershell -NoProfile -Command "$p='C:\Program Files\Python312'; $s=\"$p\Scripts\"; if(Test-Path $p){$k=[Microsoft.Win32.Registry]::LocalMachine.OpenSubKey('SYSTEM\CurrentControlSet\Control\Session Manager\Environment',$true); $c=$k.GetValue('Path','',[Microsoft.Win32.RegistryValueOptions]::DoNotExpandEnvironmentNames); $t=$k.GetValueKind('Path'); $new=$c; if((';'+$new+';') -notlike \"*;$p;*\"){$new=$p+';'+$new}; if((';'+$new+';') -notlike \"*;$s;*\"){$new=$s+';'+$new}; if($new -ne $c){$k.SetValue('Path',$new,$t)}; $k.Close()}"

REM 現在のセッションにも反映（システムPATHを再取得して連結）
for /f "usebackq tokens=2,*" %A in (`reg query "HKLM\SYSTEM\CurrentControlSet\Control\Session Manager\Environment" /v Path`) do set "PATH=%B"

REM pip / wheel の更新
python -m pip install --no-user -U pip wheel
if not "%ERRORLEVEL%"=="0" ( color 0c & echo pip / wheel の更新に失敗しました & ping 127.0.0.1 -n 6 >nul & color )

方法 2：インストーラーによるインストール

  1. Python公式サイト（https://www.python.org/downloads/）にアクセスし、「Download Python 3.x.x」ボタンからWindows用インストーラーをダウンロードする。
  2. ダウンロードしたインストーラーを実行する。
  3. 初期画面の下部に表示される「Add python.exe to PATH」にチェックを入れてから「Customize installation」を選択する。このチェックを入れ忘れると、コマンドプロンプトから python コマンドを実行できない。
  4. 「Install Python 3.xx for all users」にチェックを入れ、「Install」をクリックする。

インストールの確認

コマンドプロンプトで以下を実行する。

python --version

バージョン番号（例：Python 3.12.x）が表示されればインストール成功である。「'python' は、内部コマンドまたは外部コマンドとして認識されていません。」と表示される場合は、インストールが正常に完了していない。

第3章 Python の開発環境 Visual Studio Code のインストールと Python 用の設定

Python の開発環境Visual Studio Code（プログラムを編集するソフトウェア。以下、VS Code）を整える。

[Windows での Visual Studio Code のインストールと Python 用の設定手順を見るには、ここをクリック]

Windows での Visual Studio Code のインストールと Python 用の設定手順

1. VS Code と拡張機能のインストール

以下のコマンドにより，既存の VS Code を削除し，全ユーザー共有の設定で再インストールしたうえで，拡張機能（VS Code に機能を追加するソフトウェア）をまとめて導入する．

【インストールコマンドの実行方法】

管理者権限でコマンドプロンプトを起動する（手順：Windows キーまたはスタートメニュー → cmd と入力 → 右クリック → 「管理者として実行」）。そして，コマンド全体をコマンドプロンプトにコピー＆ペーストする。

インストールコマンド


REM ============================================================
REM Microsoft Visual Studio Code
REM ============================================================
REM Build Tools + Desktop development with C++（VCTools）+ 追加コンポーネント（一括）
REM 未インストール時: winget で新規インストール
REM インストール済み時: setup.exe modify でコンポーネント追加（バージョンは変更しない）
winget list --id Microsoft.VisualStudio.BuildTools 2>nul | findstr /i "BuildTools" >nul 2>&1
if %ERRORLEVEL% EQU 0 (
    for /f "usebackq delims=" %P in (`"C:\Program Files (x86)\Microsoft Visual Studio\Installer\vswhere.exe" -products Microsoft.VisualStudio.Product.BuildTools -property installationPath`) do start /wait "" "C:\Program Files (x86)\Microsoft Visual Studio\Installer\setup.exe" modify --installPath "%P" --add Microsoft.VisualStudio.Workload.VCTools --add Microsoft.VisualStudio.Workload.MSBuildTools --add Microsoft.VisualStudio.Component.VC.CMake.Project --add Microsoft.VisualStudio.Component.VC.Llvm.Clang --add Microsoft.VisualStudio.Component.VC.Llvm.ClangToolset --add Microsoft.VisualStudio.Component.Windows11SDK.26100 --add Microsoft.VisualStudio.Component.VC.v143.x86.x64 --includeRecommended --quiet --norestart --nocache
    if not "%ERRORLEVEL%"=="0" ( color 0c & echo Build Tools のコンポーネント追加に失敗しました & ping 127.0.0.1 -n 6 >nul & color )
) else (
    winget install --scope machine --id Microsoft.VisualStudio.BuildTools -e --silent --disable-interactivity --force --accept-source-agreements --accept-package-agreements --override "--quiet --wait --norestart --nocache --add Microsoft.VisualStudio.Workload.VCTools --includeRecommended --add Microsoft.VisualStudio.Workload.MSBuildTools --add Microsoft.VisualStudio.Component.VC.CMake.Project --add Microsoft.VisualStudio.Component.VC.Llvm.Clang --add Microsoft.VisualStudio.Component.VC.Llvm.ClangToolset --add Microsoft.VisualStudio.Component.Windows11SDK.26100 --add Microsoft.VisualStudio.Component.VC.v143.x86.x64"
    if not "%ERRORLEVEL%"=="0" ( color 0c & echo Build Tools のインストールに失敗しました & ping 127.0.0.1 -n 6 >nul & color )
)

REM 全ユーザー共有の拡張機能フォルダ
if not exist "C:\ProgramData\vscode-extensions" mkdir "C:\ProgramData\vscode-extensions"
icacls "C:\ProgramData\vscode-extensions" /grant "Everyone:(OI)(CI)M" /T

REM スタートメニューのショートカットを --extensions-dir 付きで再作成
if exist "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code" rmdir /s /q "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code"
if exist "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk" del "C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk"
powershell -NoProfile -Command "$s=New-Object -ComObject WScript.Shell; $lnk=$s.CreateShortcut('C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk'); $lnk.TargetPath='C:\Program Files\Microsoft VS Code\Code.exe'; $lnk.Arguments='--extensions-dir \"C:\ProgramData\vscode-extensions\"'; $lnk.Save()"
REM ショートカットの検証
powershell -NoProfile -Command "$s=New-Object -ComObject WScript.Shell; $lnk=$s.CreateShortcut('C:\ProgramData\Microsoft\Windows\Start Menu\Programs\Visual Studio Code.lnk'); Write-Host 'TargetPath:' $lnk.TargetPath; Write-Host 'Arguments:' $lnk.Arguments"

REM ファイル / フォルダ右クリックの「Code で開く」を登録
reg add "HKLM\SOFTWARE\Classes\*\shell\VSCode\command" /ve /d "\"C:\Program Files\Microsoft VS Code\Code.exe\" --extensions-dir \"C:\ProgramData\vscode-extensions\" \"%1\"" /f
reg add "HKLM\SOFTWARE\Classes\Directory\shell\VSCode\command" /ve /d "\"C:\Program Files\Microsoft VS Code\Code.exe\" --extensions-dir \"C:\ProgramData\vscode-extensions\" \"%1\"" /f
reg add "HKLM\SOFTWARE\Classes\Directory\Background\shell\VSCode\command" /ve /d "\"C:\Program Files\Microsoft VS Code\Code.exe\" --extensions-dir \"C:\ProgramData\vscode-extensions\" \"%V\"" /f

REM --extensions-dir 付きで起動する code.cmd ラッパを作成
REM （%* を echo で書くと対話的 cmd で失われるため、PowerShell で [char]37+'*' を書き出す）
powershell -NoProfile -Command "$pct=[char]37; $q=[char]34; $c='@echo off'+[char]13+[char]10+$q+'C:\Program Files\Microsoft VS Code\bin\code.cmd'+$q+' --extensions-dir '+$q+'C:\ProgramData\vscode-extensions'+$q+' '+$pct+'*'+[char]13+[char]10; [IO.File]::WriteAllText('C:\ProgramData\vscode-extensions\vscode.cmd',$c,[Text.Encoding]::ASCII)"

REM 拡張機能のインストール
set "CODE=C:\Program Files\Microsoft VS Code\bin\code.cmd"
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --uninstall-extension GitHub.copilot
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --uninstall-extension GitHub.copilot-chat
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension ms-python.python
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension ms-python.vscode-pylance
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension ms-python.debugpy
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension MS-CEINTL.vscode-language-pack-ja
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension saoudrizwan.claude-dev
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension rust-lang.rust-analyzer
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension tamasfe.even-better-toml
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension anthropic.claude-code
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --install-extension almenon.arepl
"%CODE%" --extensions-dir "C:\ProgramData\vscode-extensions" --list-extensions --show-versions

REM settings.json を作成（自動更新オフ、Python、Claude Code 設定）
if not exist "%APPDATA%\Code\User" mkdir "%APPDATA%\Code\User"
python -c "import json,os;data={'update.mode':'none','update.enableWindowsBackgroundUpdates':False,'extensions.autoUpdate':False,'python.defaultInterpreterPath':r'C:\Program Files\Python312\python.exe','claudeCode.environmentVariables':[{'name':'ANTHROPIC_API_KEY','value':'not-needed'},{'name':'ANTHROPIC_AUTH_TOKEN','value':'ollama'},{'name':'ANTHROPIC_BASE_URL','value':'http://localhost:11434'},{'name':'ANTHROPIC_MODEL','value':'glm-4.7-flash'},{'name':'CLAUDE_CODE_DISABLE_NONESSENTIAL_TRAFFIC','value':'1'}]};p=os.path.join(os.environ['APPDATA'],'Code','User','settings.json');open(p,'w',encoding='utf-8').write(json.dumps(data,indent=4));print('Done:',p)"

REM 自動更新の抑止ポリシー（settings.json に加えて、レジストリ側でも明示的にオフ）
reg add "HKLM\SOFTWARE\Policies\Microsoft\VSCode" /v "UpdateMode" /t REG_SZ /d "none" /f
echo === セットアップ完了 ===

2. Python インタプリタの選択

同一マシンに複数の Python がインストールされている場合，VS Code で使用する Python 本体（インタプリタ：Python プログラムを解釈・実行するソフトウェア）を選択する必要がある．

  1. コマンドパレット（コマンド名で機能を呼び出す VS Code の入力欄）を開く（Ctrl+Shift+P）
  2. Python: Select Interpreter と入力する
  3. 表示される一覧から，使用する Python（例：C:\Program Files\Python312\python.exe）を選択する．

第4章 必要なライブラリのインストール

必要なライブラリをシステム領域にインストール

管理者権限でコマンドプロンプトを起動する （手順：Windowsキーまたはスタートメニュー → cmd と入力 → 右クリック → 「管理者として実行」）。

起動したコマンドプロンプトで次を実行する。--no-user は、ユーザ領域ではなくシステム領域へインストールするためのオプションである。matplotlib-fontja は、Matplotlibのグラフに日本語を表示するためのパッケージである。

pip install --no-user numpy pillow matplotlib trimesh pymeshlab scipy matplotlib-fontja

第5章 国土地理院DEM PNG to OBJ変換ツールプログラム

概要

このプログラムは、国土地理院が配信するPNG形式の標高タイルのファイルを読み込み、各ピクセルのRGB値から標高値を復元して、3DメッシュのOBJファイルとして出力する[1]。入力に使うPNGファイルは、標高タイルをダウンロードして得る。

主要技術

処理の流れ

参考文献

ソースコード 

"""
国土地理院DEM PNG to OBJ変換ツールプログラム

特徴技術名: 標高タイル（基盤地図情報数値標高モデル）
出典: 国土地理院. 標高タイルの詳細仕様. https://maps.gsi.go.jp/development/demtile.html

特徴機能: RGB値による標高のエンコード
PNG画像のRGB値から標高値を0.01m単位でデコードする。
計算式: x = R*65536 + G*256 + B
        h = x*0.01              (x < 2^23)
        h = (x - 2^24)*0.01     (x > 2^23)
        x = 2^23 すなわち (R,G,B)=(128,0,0) は無効値

学習済みモデル: なし

方式設計:
  関連利用技術:
  - tkinter: GUIフレームワーク（ウィンドウ、ボタン、ファイル選択）
  - Pillow: 画像処理ライブラリ（画像の読み込み）
  - NumPy: 数値計算ライブラリ（標高データの配列処理）
  - Matplotlib: 標高データのプレビュー表示
  - SciPy: Delaunay三角形分割、無効値の補間、平滑化
  - trimesh: 三角形メッシュの保持とOBJファイルの出力
  - PyMeshLab: Quadric Error Metrics によるメッシュの簡略化

  入力と出力:
  入力: PNGファイル（標高タイル形式）
  出力: OBJファイル（3Dメッシュ）

  処理手順:
  1. PNG読込
  2. RGB値から標高値へのデコード
  3. 点群生成
  4. Delaunay三角形分割
  5. QEMによる簡略化
  6. OBJ出力

  前処理: 無効値 (R,G,B)=(128,0,0) のマスク処理
  後処理: メッシュの境界保持・トポロジー保持

  追加処理: 地形特徴に基づく適応的簡略化。勾配、曲率、粗さから重要度を計算し、地形の特徴を保持する。

  調整を必要とする設定値:
  - 簡略化率（1-100%）: メッシュの詳細度を制御

将来方策: 地形特徴の重み付けパラメータの自動調整

その他の重要事項: 大きな標高データの処理にはメモリ制限がある

前準備:
pip install --no-user numpy pillow matplotlib trimesh pymeshlab scipy matplotlib-fontja
"""

import os
import tkinter as tk
from tkinter import filedialog, messagebox, ttk
import numpy as np
from PIL import Image
import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
import matplotlib_fontja
from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
from scipy.spatial import Delaunay
from scipy.interpolate import griddata
from scipy.ndimage import gaussian_filter
import trimesh
import pymeshlab

# 設定値
WINDOW_WIDTH = 600  # ウィンドウ幅
WINDOW_HEIGHT = 600  # ウィンドウ高さ
DEFAULT_SIMPLIFICATION = 10  # 既定の簡略化率(%)
PREVIEW_SIZE = 3  # プレビュー図のサイズ(インチ)
PREVIEW_DPI = 100  # プレビュー図のDPI
QUALITY_THRESHOLD = 0.3  # QEM簡略化の品質閾値
# 地形特徴の重み
SLOPE_WEIGHT = 0.4  # 勾配の重み
CURVATURE_WEIGHT = 0.3  # 曲率の重み
ROUGHNESS_WEIGHT = 0.3  # 粗さの重み

# グローバル変数
root = None
input_file = None
elevation_data = None
file_label = None
status_text = None
save_button = None
simplification_var = None
simplification_label = None
fig = None
ax = None
canvas = None
colorbar = None


def setup_ui():
    """UIコンポーネントの設定"""
    global root, file_label, status_text, save_button, simplification_var
    global simplification_label, fig, ax, canvas

    root = tk.Tk()
    root.title('DEM PNG to OBJ Converter')
    root.geometry(f'{WINDOW_WIDTH}x{WINDOW_HEIGHT}')

    print("=== DEM PNG to OBJ変換ツール ===")
    print("概要説明:")
    print("国土地理院の標高タイル（PNG形式）のファイルを3DのOBJファイルに変換します")
    print("操作方法:")
    print("1. 'ファイル選択'ボタンで標高タイル（PNG形式）のファイルを選択")
    print("2. 簡略化率スライダーでメッシュの詳細度を調整（1-100%）")
    print("3. '保存'ボタンでOBJファイルとして出力")
    print("注意事項:")
    print("・大きなファイルの処理には時間がかかります")
    print("・RGB値(128,0,0)は無効値として処理されます")
    print("・出力ファイルは入力ファイルと同じフォルダに保存されます")
    print("=====================================")

    # ファイル選択フレーム
    file_frame = ttk.Frame(root, padding='5')
    file_frame.grid(row=0, column=0, sticky=(tk.W, tk.E))

    ttk.Label(file_frame, text='入力PNGファイル:').grid(row=0, column=0, sticky=tk.W)
    file_label = ttk.Label(file_frame, text='未選択', relief=tk.SUNKEN, width=40)
    file_label.grid(row=0, column=1, padx=5)

    ttk.Button(file_frame, text='ファイル選択', command=select_file).grid(row=0, column=2)

    # プレビューフレーム
    preview_frame = ttk.LabelFrame(root, text='標高プレビュー', padding='5')
    preview_frame.grid(row=1, column=0, padx=5, pady=5, sticky=(tk.W, tk.E, tk.N, tk.S))

    # Matplotlibの図を作成
    fig = plt.Figure(figsize=(PREVIEW_SIZE, PREVIEW_SIZE), dpi=PREVIEW_DPI)
    ax = fig.add_subplot(111)
    canvas = FigureCanvasTkAgg(fig, master=preview_frame)
    canvas.get_tk_widget().pack()

    # パラメータフレーム
    param_frame = ttk.LabelFrame(root, text='パラメータ設定', padding='5')
    param_frame.grid(row=2, column=0, padx=5, pady=5, sticky=(tk.W, tk.E))

    ttk.Label(param_frame, text='簡略化率 (%):').grid(row=0, column=0, sticky=tk.W)
    simplification_var = tk.IntVar(value=DEFAULT_SIMPLIFICATION)
    simplification_slider = ttk.Scale(
        param_frame,
        from_=1,
        to=100,
        orient=tk.HORIZONTAL,
        variable=simplification_var,
        length=300,
        command=update_slider_label
    )
    simplification_slider.grid(row=0, column=1, padx=5)

    simplification_label = ttk.Label(param_frame, text=f'{DEFAULT_SIMPLIFICATION}%')
    simplification_label.grid(row=0, column=2)

    # 保存ボタン
    save_button = ttk.Button(
        root,
        text='保存',
        command=save_mesh,
        state=tk.DISABLED
    )
    save_button.grid(row=3, column=0, pady=5)

    # ステータスフレーム
    status_frame = ttk.LabelFrame(root, text='ステータス', padding='5')
    status_frame.grid(row=4, column=0, padx=5, pady=5, sticky=(tk.W, tk.E))

    status_text = tk.Text(status_frame, height=3, width=60)
    status_text.pack()


def update_slider_label(value):
    """簡略化率の整数化とラベルの更新"""
    simplification = int(float(value))
    simplification_var.set(simplification)
    simplification_label.config(text=f'{simplification}%')


def select_file():
    """ファイル選択ダイアログ"""
    global input_file

    filename = filedialog.askopenfilename(
        title='PNG ファイルを選択',
        filetypes=[('PNG files', '*.png'), ('All files', '*.*')]
    )

    if filename:
        input_file = filename
        file_label.config(text=os.path.basename(filename))
        load_and_preview()
        save_button.config(state=tk.NORMAL)


def decode_elevation_png(image_array):
    """標高タイル（PNG形式）のRGB値を標高値[m]にデコードする"""
    R = image_array[:, :, 0].astype(np.uint32)
    G = image_array[:, :, 1].astype(np.uint32)
    B = image_array[:, :, 2].astype(np.uint32)

    # x = R*65536 + G*256 + B、u は標高分解能（0.01m）
    x = R * 65536 + G * 256 + B
    u = 0.01

    # h = x*u (x < 2^23)、h = (x - 2^24)*u (x > 2^23)
    elevation = np.where(x < 8388608,
                         x.astype(np.float64) * u,
                         (x.astype(np.float64) - 16777216) * u)

    # x = 2^23 すなわち (R,G,B)=(128,0,0) は無効値
    elevation[x == 8388608] = np.nan

    return elevation


def load_and_preview():
    """PNGファイルを読み込んでプレビュー表示"""
    global elevation_data, colorbar

    # 画像読み込み
    img = Image.open(input_file)

    # 標高値デコード
    elevation_data = decode_elevation_png(np.array(img))

    # プレビュー表示
    ax.clear()

    # 無効値をマスク
    masked_elevation = np.ma.masked_invalid(elevation_data)

    # カラーマップで表示
    im = ax.imshow(
        masked_elevation,
        cmap='terrain',
        aspect='equal'
    )

    # カラーバー追加
    if colorbar:
        colorbar.remove()
    colorbar = fig.colorbar(im, ax=ax)
    colorbar.set_label('標高 (m)')

    ax.set_title('標高データ')
    ax.set_xlabel('X (ピクセル)')
    ax.set_ylabel('Y (ピクセル)')

    fig.tight_layout()
    canvas.draw()

    # ステータス更新
    valid_count = np.sum(~np.isnan(elevation_data))
    total_count = elevation_data.size

    status_text.delete(1.0, tk.END)
    status_text.insert(tk.END, 'ファイル読み込み完了\n')
    status_text.insert(tk.END, f'画像サイズ: {img.size[0]}x{img.size[1]}\n')
    status_text.insert(tk.END, f'有効データ: {valid_count}/{total_count} ピクセル\n')
    status_text.insert(tk.END, f'標高範囲: {np.nanmin(elevation_data):.1f}m - {np.nanmax(elevation_data):.1f}m')


def compute_terrain_features(elevation):
    """地形特徴量（勾配・曲率・粗さ）の計算"""
    valid_mask = ~np.isnan(elevation)

    # 無効値を最も近い有効値で補間
    elevation_filled = elevation.copy()
    if np.any(~valid_mask):
        yy, xx = np.indices(elevation.shape)
        valid_points = np.column_stack([xx[valid_mask], yy[valid_mask]])
        valid_values = elevation[valid_mask]
        elevation_filled[~valid_mask] = griddata(
            valid_points,
            valid_values,
            (xx[~valid_mask], yy[~valid_mask]),
            method='nearest'
        )

    # 勾配計算
    grad_y, grad_x = np.gradient(elevation_filled)
    slope = np.sqrt(grad_x**2 + grad_y**2)

    # 曲率計算
    grad_xx = np.gradient(grad_x, axis=1)
    grad_yy = np.gradient(grad_y, axis=0)
    curvature = np.abs(grad_xx + grad_yy)

    # 粗さ計算
    smoothed = gaussian_filter(elevation_filled, sigma=3)
    roughness = np.abs(elevation_filled - smoothed)

    return slope, curvature, roughness


def compute_importance_map(elevation):
    """重要度マップの計算"""
    # 地形特徴を計算
    slope, curvature, roughness = compute_terrain_features(elevation)

    # 各特徴を95パーセンタイルで正規化し、0から1にクリッピング
    slope_norm = np.clip(slope / np.maximum(np.nanpercentile(slope, 95), 1e-10), 0, 1)
    curvature_norm = np.clip(curvature / np.maximum(np.nanpercentile(curvature, 95), 1e-10), 0, 1)
    roughness_norm = np.clip(roughness / np.maximum(np.nanpercentile(roughness, 95), 1e-10), 0, 1)

    # 重要度を計算
    importance = (SLOPE_WEIGHT * slope_norm +
                 CURVATURE_WEIGHT * curvature_norm +
                 ROUGHNESS_WEIGHT * roughness_norm)

    # 0-1に正規化
    importance_range = np.nanmax(importance) - np.nanmin(importance)
    if importance_range > 1e-10:
        importance = (importance - np.nanmin(importance)) / importance_range
    else:
        importance = np.ones_like(importance) * 0.5

    return importance


def create_mesh_from_elevation():
    """標高データからメッシュを生成（品質値付き）"""
    # 有効な点のインデックスを取得
    valid_indices = np.where(~np.isnan(elevation_data))

    # 点群を生成（画像座標系をそのまま使用）
    points = np.column_stack([
        valid_indices[1],  # x座標（列インデックス）
        valid_indices[0],  # y座標（行インデックス）
        elevation_data[valid_indices]  # z座標（標高値）
    ])

    # 重要度マップを計算し、有効な点の重要度を取得
    quality_values = compute_importance_map(elevation_data)[valid_indices]

    # Delaunay三角形分割
    tri = Delaunay(points[:, :2])

    # trimeshオブジェクトの作成
    mesh = trimesh.Trimesh(
        vertices=points,
        faces=tri.simplices,
        process=True
    )

    return mesh, quality_values


def simplify_mesh(mesh, quality_values, target_ratio):
    """メッシュの簡略化（地形特徴を考慮）"""
    ms = pymeshlab.MeshSet()

    # 品質値（重要度）を持つPyMeshLabメッシュを作成
    pm_mesh = pymeshlab.Mesh(
        vertex_matrix=mesh.vertices.astype(np.float64),
        face_matrix=mesh.faces.astype(np.int32),
        v_scalar_array=quality_values.astype(np.float64)
    )
    ms.add_mesh(pm_mesh)

    # 簡略化実行
    ms.apply_filter('meshing_decimation_quadric_edge_collapse',
                    targetperc=target_ratio,
                    preserveboundary=True,
                    preservenormal=True,
                    preservetopology=True,
                    optimalplacement=True,
                    planarquadric=True,
                    qualitythr=QUALITY_THRESHOLD,
                    qualityweight=True)

    # 簡略化されたメッシュを取得
    simplified_mesh = ms.current_mesh()

    return trimesh.Trimesh(
        vertices=simplified_mesh.vertex_matrix(),
        faces=simplified_mesh.face_matrix()
    )


def save_mesh():
    """メッシュ保存"""
    save_button.config(state=tk.DISABLED)

    # 出力ファイル名の生成
    base_name = os.path.splitext(os.path.basename(input_file))[0]
    output_file = os.path.join(
        os.path.dirname(input_file),
        f'{base_name}.obj'
    )

    # メッシュ生成
    status_text.delete(1.0, tk.END)
    status_text.insert(tk.END, 'メッシュ生成中...\n')
    root.update()

    mesh, quality_values = create_mesh_from_elevation()
    original_faces = len(mesh.faces)

    # 簡略化
    simplification_ratio = simplification_var.get() / 100.0
    status_text.insert(tk.END, f'簡略化実行中 ({int(simplification_ratio*100)}%)...\n')
    status_text.insert(tk.END, '地形特徴を考慮した適応的簡略化を適用中...\n')
    root.update()

    simplified_mesh = simplify_mesh(mesh, quality_values, simplification_ratio)

    # OBJファイルとして保存
    simplified_mesh.export(output_file)

    # ファイルサイズ取得
    file_size = os.path.getsize(output_file) / 1024  # KB単位

    # 完了メッセージ
    status_text.insert(tk.END, '\n保存完了\n')
    status_text.insert(tk.END, f'ファイル名: {output_file}\n')
    status_text.insert(tk.END, f'ポリゴン数: {len(simplified_mesh.faces)} (元: {original_faces})\n')
    status_text.insert(tk.END, f'ファイルサイズ: {file_size:.1f} KB')

    messagebox.showinfo('完了', '保存が完了しました')
    save_button.config(state=tk.NORMAL)


# メイン実行
if __name__ == '__main__':
    setup_ui()
    root.mainloop()