はじめに
Depth Anything 3（以下、DA3）を用いて、1枚のRGB GeoTIFFから単眼深度推定を行います。追加学習は行わず、公開済みの事前学習済みモデルを利用します。
1枚の衛星画像から単眼深度推定で得た画像を用いて、LOD2相当の3次元の建物データを生成することができます。
単眼深度推定の概要と、衛星画像から建物の3次元データを生成する取り組みについては、こちらの解説ページをご覧ください。
本記事では、その第一段階として、UbuntuとNVIDIA GPUを使った環境構築、モデル取得、深度推論、GeoTIFF出力、QGISでの確認までを扱います。建物ポリゴン、DEM、カメラの外部標定データは入力に使用しません。
入力画像は、衛星画像はWorldView-2（50cm解像度）のアーカイブ画像を使用しましたが、そのまま使用すると非可逆圧縮によるノイズで建物境界が不鮮明のため、事前に2倍超解像及び鮮明化処理しました。
超解像化したGeotiff画像は、2048×1955pxですが、そのまま単眼深度推定を行った際には、画像全体を縮小するためか、推論した画像は全体的に境界がぼやけて見えました。そこで、元画像を300px単位に分割して推論し、重複領域を接合する方法を試しました。
この方法は、1年ちょい前にDepth Anything 2で同範囲を単眼深度推定を行った時に、全体が傾斜したような深度画像になってしまうことから得た対策です。
手元の画像では、Depth Anything 2の時と同様に、タイル化によって建物境界がより鮮明に見える結果が得られました。ただし、これは目視による比較であり、実測高さや屋根勾配の精度を確認した結果ではありません。
出力は相対深度です。メートル単位の建物高さ、DSM、DEMがそのまま得られるわけではありません。 建物3Dモデルの生成やDEM接地は、別の後処理として扱います。
1. 使用するモデルと処理の構成
使用するのは depth-anything/DA3MONO-LARGE です。公式モデルカードでは、単眼の相対深度推定用モデルとして公開されています。DA3-LARGE や DA3METRIC-LARGE とは別のモデルです。
処理は、RGB GeoTIFF読込 → DA3推論 → 元の画素グリッドへ深度をリサイズ → GeoTIFF保存、という流れです。タイル版では、推論前に分割し、保存前に深度タイルを接合します。
DA3本体には公式APIを使用します。一方、GeoTIFFの入出力とタイル接合は、本記事の追加スクリプトで実装した処理です。DA3公式の地理空間タイル推論機能ではありません。
2. 使用環境
本記事で使用する環境構成は次のとおりです。以下のコマンドはUbuntuのBash向けです。
| 項目 | 構成 |
|---|---|
| OS | Ubuntu 24.04 |
| GPU | NVIDIA RTX PRO 6000 Blackwell Max-Q |
| Python | 3.10 / conda専用環境 |
| PyTorch | 2.9.1、CUDA 12.8版wheel |
| torchvision | 0.24.1 |
| xFormers | 0.0.33.post2 |
| NumPy | 1.26.4 |
| OpenCV | opencv-python 4.10.0.84 |
| MoviePy | 1.0.3 |
| モデル | DA3MONO-LARGE |
| GeoTIFF入出力 | Rasterio |
このGPUの搭載メモリ量を必須条件としているわけではありません。必要なGPUメモリは入力解像度や環境に依存します。他GPUの最小必要容量や処理速度は本記事では評価していません。
PyTorchのバージョン組み合わせは公式の旧バージョン一覧を、xFormersの対応関係は0.0.33.post2のリリース情報を参照しています。最新版へ一括更新する手順ではなく、使用する主要パッケージの組み合わせを固定します。
nvidia-smi の「CUDA Version」、nvcc --version が示すToolkit、torch.version.cuda が示すPyTorch側のCUDAバージョンは、同じ意味ではありません。本手順はCUDA 12.8版のPyTorch wheelを使います。システムのCUDA Toolkitを入れ替える操作は含めません。GPUを利用できるNVIDIAドライバーは事前に導入してください。
3. 作業ディレクトリとconda環境を作る
MinicondaまたはAnacondaが利用できることを前提にします。
conda create -n da3_bw python=3.10 pip -y
conda activate da3_bw
mkdir -p ~/da3_geotiff/input ~/da3_geotiff/output ~/da3_geotiff/logs
cd ~/da3_geotiff
nvidia-smi
git --version
tmux -V
git や tmux がない場合は、Ubuntuのパッケージとして導入します。
sudo apt-get update
sudo apt-get install -y git tmux
記事末尾の「付録：使用スクリプト」にあるファイルを、すべて ~/da3_geotiff/ 直下へ保存してください。 各折りたたみ内のコードをコピーし、見出しに示した名前で保存します。別途ダウンロードする独自ZIPや非公開リポジトリは必要ありません。
| ファイル | 役割 |
|---|---|
constraints.txt |
主要パッケージの固定 |
setup.sh |
公式コードと依存関係の導入 |
check_env.py |
CUDAとimportの確認 |
download_model.py |
公式の重み取得 |
infer_geotiff.py |
画像全体を使った推論 |
da3_tile_io.py |
タイル版のGeoTIFF入出力 |
infer_geotiff_tiled.py |
分割推論と接合 |
4. DA3と依存パッケージをインストールする
conda activate da3_bw
cd ~/da3_geotiff
bash setup.sh
setup.sh は次を実行します。
- CUDA 12.8版のPyTorch・torchvision・xFormersを導入。
- 公式リポジトリを
vendor/Depth-Anything-3にclone。 - DA3の基本依存関係とRasterio等を導入。
-
pip check、import、CUDAでの行列演算を確認。 - GitのコミットIDと
pip freezeをlogs/に保存。
公式コードは以下から取得します。
# setup.sh内で実行されるため、手動で重ねて実行する必要はありません。
git clone https://github.com/ByteDance-Seed/Depth-Anything-3.git vendor/Depth-Anything-3
基本依存関係のみを導入し、Gradioや3D Gaussian出力用の追加構成はインストールしません。モデルの深度推論では infer_gs=False を指定します。
環境を再確認する場合は次を実行します。
python -m pip check
python check_env.py
CUDA runtime: 12.8、GPU名、Imports and CUDA matrix multiplication: OK が確認項目です。ただし、この確認だけでは重みを使ったモデル推論まで成功したことにはなりません。
主要パッケージは固定していますが、全依存関係と公式リポジトリのコミットをあらかじめ固定した完全なロック環境ではありません。将来の変更で依存関係が競合する場合は、制約を外して更新する前に、エラーログと取得コミットを確認してください。
5. DA3MONO-LARGEの重みを取得する
python download_model.py
モデルは公式Hugging Faceの depth-anything/DA3MONO-LARGE から取得し、models/DA3MONO-LARGE/ に保存します。ダウンロード時に解決したリビジョンを指定して取得し、そのIDを download_revision.txt に記録します。
cat models/DA3MONO-LARGE/download_revision.txt
cat logs/da3_git_commit.txt
これらは取得時のバージョンを記録するもので、別の日に新規ダウンロードしても同じリビジョンになることを保証するものではありません。再現用にはモデルフォルダと環境記録も保管します。
推論ではローカルに保存した重みを読み込みます。入力画像をHugging Faceへアップロードする処理はありません。
6. 入力GeoTIFFを用意する
例として input/area_rgb.tif を使用します。ご自身の画像をこの名前で配置するか、以下のコマンドの入力パスを変更してください。
入力条件は次のとおりです。
- 元のRGB画像で、画素型は8bitの
uint8。 - CRSとアフィン変換が設定されたGeoTIFF。
- 標準では1・2・3バンドをR・G・Bとして使用。
- すでに推定済みの深度画像や高さ画像は入力しない。
16bit衛星画像を扱う場合は、画像全体で一貫した方法によってRGBの8bit画像へ変換してから入力します。単純な astype(uint8) や、タイルごとの独立したコントラスト調整は避けます。本スクリプトは16bit画像を自動変換せず停止します。
入力の確認例です。
python - <<'PY'
import rasterio
with rasterio.open('input/area_rgb.tif') as src:
print('size:', src.width, src.height)
print('bands:', src.count)
print('dtype:', src.dtypes)
print('CRS:', src.crs)
print('transform:', src.transform)
print('color interpretation:', src.colorinterp)
PY
RGBの順序が異なる場合は、推論時に --bands R番号 G番号 B番号 を追加します。例えばR・G・Bが3・2・1バンドなら --bands 3 2 1 です。
4バンド画像の追加バンドはモデル入力にしません。ただし、Rasterioが認識するalphaやNoDataマスクは有効画素判定に使用します。欠損画素は推論時にRGB中央値で補い、出力ではNoDataに戻します。欠損領域付近の予測への影響まで除去できるわけではありません。
7. まず画像全体で推論する
長時間処理には tmux を使います。
tmux new -s da3
開いたセッション内で次を実行します。
conda activate da3_bw
cd ~/da3_geotiff
mkdir -p output logs
set -o pipefail
python -u infer_geotiff.py \
--input input/area_rgb.tif \
--output-prefix output/area_da3_full_r1008 \
--process-res 1008 \
2>&1 | tee logs/area_da3_full_r1008.log
--process-res 1008 は、モデルへ入力する画像の長辺の目安です。元画像を1008pxのタイルに分割する意味ではありません。前処理では切り抜かずにリサイズし、パッチサイズに合わせた解像度調整を行います。
出力深度は元画像と同じ縦横画素数へ双線形補間して保存します。そのため、出力画素数が入力と同じでも、縮小推論で失われた細部が復元されたわけではありません。
Ctrl+B の後に D でtmuxから離れ、次で再接続できます。
tmux attach -t da3
8. 建物の細部を見るため300pxタイルで推論する
画像全体を縮小せず、元画像の小領域ごとに推論します。ここでは300pxの切り出し、100pxの重複、518pxのモデル入力を使います。
python -u infer_geotiff_tiled.py \
--input input/area_rgb.tif \
--output-prefix output/area_da3_t300_o100_r518_shift \
--tile-size 300 \
--overlap 100 \
--process-res 518 \
--alignment shift \
2>&1 | tee logs/area_da3_t300_o100_r518_shift.log
| 引数 | 意味 |
|---|---|
--tile-size 300 |
元画像を切り出す大きさ。通常は300×300px |
--overlap 100 |
隣接タイルの重複幅。通常の移動幅は200px |
--process-res 518 |
300pxのRGBタイルをモデルに渡す際の入力サイズ |
--alignment shift |
重複部分の深度差からオフセットだけを補正 |
--output-prefix |
出力ファイル群の共通名 |
518は14の倍数です。このラッパーは process-res に14の倍数を要求します。300pxから518pxへの拡大はモデルへの入力方法であり、画像の実際の地上分解能を高める処理ではありません。
例えば0.25m/pxの画像では、300pxは75m四方です。0.5m/pxなら150m四方になり、同じ300pxでも建物が占める割合や周囲の文脈が変わります。別地域や異なる解像度の画像で同じ精度が出るとは限りません。
端のタイルは切り出し位置を戻して画像内に収めるため、末端だけ重複幅が100pxより大きくなることがあります。画像の幅または高さが300px未満の場合は、端の画素値で補います。2048×1955pxなら10×10の100タイルです。
モデルは最初に一度だけ読み込み、各タイルを順に推論します。各予測を元の300pxへ戻して接合します。
接合方法：shift・affine・none
| 指定 | 補正 | 特徴 |
|---|---|---|
shift |
オフセットだけ | 各タイルの深度差の倍率を変更しない |
affine |
倍率とオフセット | 重複部で尺度を合わせるが、局所的な起伏の大きさも変わる |
none |
値の補正なし | 元の予測を加重平均。タイルごとの差が残ることがある |
今回の運用では、倍率補正が局所的な高低差へ与える影響を避けて比較するため、shift を明示しています。掲載コードの既定値は affine なので、引数を省略しないでください。
shiftでは、現在のタイルの深度を $d$、既に接合した重複領域の深度を $m$ とし、有効な重複画素について
$$
b = \operatorname{median}(m-d),\qquad d' = d+b
$$
として補正します。最初のタイルや有効な重複部がないタイルでは、オフセットを0とします。その後、各タイルの中央側を重くした加重平均で接合します。
$$
D(x,y)=\frac{\sum_i w_i(x,y)d_i'(x,y)}{\sum_i w_i(x,y)}
$$
ここでの「位置合わせ」は画像を幾何変換する処理ではなく、同じ画素位置で予測された深度値の基準を合わせる処理です。
タイルごとの最小・最大値で0〜1へ正規化してから接合することはしません。ただし、shiftでも真の尺度差を解消できるわけではなく、走査順に依存する誤差の蓄積や、重複領域の平均化によるぼけが残り得ます。
9. 別の場所をゼロから推論する
入力と出力名を変えて、同じコマンドを実行します。以前の場所の深度結果やタイルを引き継ぐ処理はありません。
python -u infer_geotiff_tiled.py \
--input input/another_area_rgb.tif \
--output-prefix output/another_area_da3_t300_o100_r518_shift \
--tile-size 300 --overlap 100 --process-res 518 \
--alignment shift \
2>&1 | tee logs/another_area_da3_t300_o100_r518_shift.log
既存出力の上書きは防止しています。同じ画像を再推論する場合も、未使用の --output-prefix を指定してください。本記事のDA3スクリプトには途中再開機能はありません。
10. 出力ファイルの意味
| ファイルの末尾 | 内容 | 主な用途 |
|---|---|---|
_depth.tif |
Float32の相対深度。値が大きいほど遠い | 数値解析・後段処理 |
_near01.tif |
近い側を1にした0〜1表示用GeoTIFF | QGISでの表示 |
_near.png |
同じ表示の8bit PNG | 画像比較・記事掲載 |
_metadata.json |
推論条件、GPU、画像サイズ、タイル補正量など | 実行条件の確認 |
GeoTIFFは、入力のCRS、アフィン変換、縦横画素数を保持します。出力のNoData値は -9999 です。入力の無効画素位置を保持しますが、任意のTIFFタグや元のNoData値そのものをすべて複製するわけではありません。
タイル版の _depth.tif は接合補正後の相対深度です。個々のタイルの未補正予測ではありません。shiftやaffineの補正後には負の値が生じる場合もあります。本スクリプトは未接合タイルの予測を個別ファイルとして保存しません。
表示用near01は建物高さではない
接合した深度全体の1・99パーセンタイルを $P_1,P_{99}$ とし、表示用画像は
$$
N=\operatorname{clip}\left(\frac{P_{99}-D}{P_{99}-P_1},0,1\right)
$$
で作ります。近い側が白、遠い側が黒です。これはコントラスト調整と向きの反転であり、逆深度 $1/D$ やメートル高さへの変換ではありません。
画像ごとに表示範囲を決めるため、別地域のnear01の同じ明るさが同じ深度や高さを表すわけではありません。数値処理には _depth.tif を使います。
11. QGISで確認する
- 元のRGB GeoTIFFと
_near01.tifをQGISに読み込みます。 -
_near01.tifを「単バンドグレー」で表示します。 - 最小値0、最大値1、黒から白の配色に設定します。
- 同じ場所・同じ縮尺で、建物境界、棟線、屋根面、道路との関係を確認します。
_depth.tif を表示する場合は、そのファイルの値域に合わせてストレッチしてください。0〜1固定で表示すると、多くの画素が白や黒に張り付くことがあります。白一色という見た目だけでは、推論失敗とは判断できません。
座標と出力値の確認例です。
python - <<'PY'
import numpy as np
import rasterio
source='input/area_rgb.tif'
result='output/area_da3_t300_o100_r518_shift_depth.tif'
with rasterio.open(source) as src, rasterio.open(result) as dst:
assert src.crs == dst.crs
assert src.transform == dst.transform
assert (src.width, src.height) == (dst.width, dst.height)
source_valid=np.all(src.read_masks([1,2,3]) > 0, axis=0)
depth=dst.read(1, masked=True)
assert np.array_equal(source_valid, ~np.ma.getmaskarray(depth))
values=depth.compressed()
assert values.size and np.isfinite(values).all()
print('CRS・グリッド・有効画素位置: OK')
print('NoData:', dst.nodata)
print('min / median / max:', values.min(), np.median(values), values.max())
PY
RGBバンドを変更した場合は、確認コードの [1,2,3] も合わせます。このチェックはデータの位置と形式の確認であり、深度の正解率を評価するものではありません。
推論結果の例
以下に、入力画像とDepth Anything 3による単眼深度推定結果を示します。入力には、WorldView-2のアーカイブ画像に2倍の超解像処理を適用した画像を使用しています。
入力画像とDepth Anything 3による単眼深度推定結果
左：入力画像（2倍超解像処理後）。右：推定した相対深度を、近い側ほど明るく表示した画像。
建物と周辺地表の違いや、一部の屋根面の起伏が濃淡として表れています。一方、画像内には広い範囲にわたる明暗の変化も見られます。この変化には、地形の影響や推定・タイル接合による誤差が含まれる可能性があり、この画像だけでは原因を区別できません。
したがって、屋根の明るさをそのまま建物高さとして解釈せず、建物ごとの周辺地盤との相対差や、実測値との対応を確認する必要があります。
12. よくある問題
| 症状 | 確認・対応 |
|---|---|
CUDA GPU required / CUDAが認識されない |
nvidia-smi と python check_env.py、有効なconda環境を確認 |
| CUDAメモリ不足 | 全体版の process-res を下げるか、タイル版を使用。タイル版なら518から392など14の倍数へ変更して比較 |
xFormers のバージョン不整合 |
PyTorchとxFormersの組み合わせを確認。環境内で片方だけを更新しない |
Model directory missing |
python download_model.py を完了させる。別配置なら --model を指定 |
Input must be RGB uint8 |
元RGBか、バンド選択と型が正しいか確認 |
| 出力が真っ白・真っ黒 | 表示範囲を確認。まずnear01を0〜1で表示し、深度の最小・最大値も確認 |
| タイル境界で深度が不連続 | 重複幅と接合方法を比較。真の高さ差と接合誤差を区別する |
| 建物が全体として高すぎる・低すぎる | 相対深度をメートル高さと解釈していないか確認 |
Output exists |
未使用の出力prefixでゼロから実行 |
この実装は、モデルへの入力をタイル化しますが、元RGB画像と接合用配列はCPUメモリへ読み込みます。巨大GeoTIFFを一定メモリで処理するストリーミング実装ではありません。広域画像ではAOIを切り出して使用してください。
13. 評価するときに分けて考えること
建物境界が鮮明に見えること、深度の局所的な形状が正しいこと、建物高さが実測に一致することは、それぞれ別の評価です。
比較する場合は、同じ元画像・同じ切り出し範囲・同じ表示倍率を使い、全体版とタイル版、接合方法の違いを確認します。DA2等と比較する場合も、内部入力解像度や出力値の表現を揃えて確認する必要があります。
今回のGeoTIFF化は、画像の位置情報を引き継ぐ処理です。深度を地盤基準の高さや標高に校正するものではありません。DSMや建物3Dモデルへ進む場合は、別途、深度と高さの対応、地盤の扱い、実測値による検証が必要です。
参考資料
- Depth Anything 3公式リポジトリ
- DA3MONO-LARGEモデルカード
- PyTorch：Previous Versions
- xFormers 0.0.33.post2：PyTorch 2.9.1向けwheel
付録：使用スクリプト
以下は本記事のGeoTIFFラッパーです。各コードブロックを指定のファイル名で作業ディレクトリ直下へ保存してください。公式モデルの実装や重みは含めず、前述の手順で公式配布元から取得します。
constraints.txt
constraints.txt の全文
constraints.txt
torch==2.9.1
torchvision==0.24.1
xformers==0.0.33.post2
numpy==1.26.4
opencv-python==4.10.0.84
moviepy==1.0.3
setup.sh
setup.sh の全文
setup.sh
#!/usr/bin/env bash
set -euo pipefail
cd "$(dirname "$0")"
if [[ "${CONDA_DEFAULT_ENV:-}" != "da3_bw" ]]; then
echo 'Activate the dedicated environment first: conda activate da3_bw' >&2
exit 1
fi
mkdir -p vendor models output logs
python -m pip install --upgrade pip
python -m pip install -c constraints.txt torch torchvision xformers \
--index-url https://download.pytorch.org/whl/cu128
if [[ ! -d vendor/Depth-Anything-3 ]]; then
git clone https://github.com/ByteDance-Seed/Depth-Anything-3.git vendor/Depth-Anything-3
fi
python -m pip install -c constraints.txt -e ./vendor/Depth-Anything-3
python -m pip install -c constraints.txt rasterio addict
python -m pip check
git -C vendor/Depth-Anything-3 rev-parse HEAD > logs/da3_git_commit.txt
python -m pip freeze > logs/environment_freeze.txt
python check_env.py
check_env.py
check_env.py の全文
check_env.py
import torch
import torchvision
import xformers
import rasterio
from depth_anything_3.api import DepthAnything3
print('torch:', torch.__version__)
print('torchvision:', torchvision.__version__)
print('xformers:', xformers.__version__)
print('CUDA runtime:', torch.version.cuda)
print('rasterio:', rasterio.__version__)
assert torch.cuda.is_available(), 'CUDA GPU is unavailable'
print('GPU:', torch.cuda.get_device_name(0))
print('Capability:', torch.cuda.get_device_capability(0))
x = torch.randn(512, 512, device='cuda')
y = x @ x
torch.cuda.synchronize()
assert torch.isfinite(y).all()
print('Imports and CUDA matrix multiplication: OK (model inference not yet tested)')
download_model.py
download_model.py の全文
download_model.py
from pathlib import Path
from huggingface_hub import HfApi, snapshot_download
root = Path(__file__).resolve().parent
repo = 'depth-anything/DA3MONO-LARGE'
info = HfApi().model_info(repo)
path = snapshot_download(
repo_id=repo, revision=info.sha,
local_dir=str(root / 'models' / 'DA3MONO-LARGE'),
allow_patterns=['*.json', '*.safetensors', '*.md', 'LICENSE*'],
)
(root / 'models' / 'DA3MONO-LARGE' / 'download_revision.txt').write_text(
f'{repo}\n{info.sha}\n', encoding='utf-8')
print('Model downloaded:', path)
infer_geotiff.py
infer_geotiff.py の全文
infer_geotiff.py
"""Single-image DA3MONO inference. Relative depth, never metric height.
Whole-image inference avoids independent tile scales. The API resizes, without
cropping; the depth is bilinearly restored to the source pixel grid. Original
CRS/affine and the intersection of RGB valid masks are retained. Unknown extra
bands are ignored; rasterio-recognized alpha masks are honored.
"""
from pathlib import Path
import argparse
import json
import time
import numpy as np
import rasterio
from PIL import Image
NODATA = -9999.0
def read_rgb(path, bands):
with rasterio.open(path) as src:
if min(bands) < 1 or max(bands) > src.count:
raise ValueError(f'Invalid RGB bands {bands}; source has {src.count} bands')
rgb = src.read(bands)
if rgb.dtype != np.uint8:
raise ValueError('Input must be RGB uint8. Convert 16-bit/float data explicitly first.')
if src.crs is None:
raise ValueError('Input has no CRS; supply a georeferenced RGB GeoTIFF')
valid = np.all(src.read_masks(bands) > 0, axis=0)
if not valid.any():
raise ValueError('Input has no valid RGB pixels')
profile = dict(driver='GTiff', width=src.width, height=src.height,
count=1, dtype='float32', crs=src.crs, transform=src.transform,
nodata=NODATA, compress='deflate', predictor=3, tiled=True,
blockxsize=256, blockysize=256, BIGTIFF='IF_SAFER')
rgb = np.moveaxis(rgb, 0, -1).copy()
# Fill masked pixels only for inference; their outputs remain NoData.
if not valid.all():
rgb[~valid] = np.median(rgb[valid], axis=0).astype(np.uint8)
print('Warning: masked areas can affect nearby predictions.', flush=True)
return rgb, valid, profile
def restore_depth(depth, height, width):
depth = np.asarray(depth, dtype=np.float32)
if depth.ndim != 2 or not np.isfinite(depth).all():
raise ValueError('Model depth must be a finite 2-D array')
if depth.shape != (height, width):
depth = np.asarray(Image.fromarray(depth).resize(
(width, height), resample=Image.Resampling.BILINEAR), dtype=np.float32)
return depth.copy()
def write_raster(path, data, valid, profile, description, tags):
arr = np.where(valid, data, NODATA).astype(np.float32)
with rasterio.open(path, 'w', **profile) as dst:
dst.write(arr, 1)
dst.set_band_description(1, description)
dst.update_tags(**tags)
v = arr[valid].astype(np.float64)
dst.update_tags(1, STATISTICS_MINIMUM=str(v.min()),
STATISTICS_MAXIMUM=str(v.max()),
STATISTICS_MEAN=str(v.mean()),
STATISTICS_STDDEV=str(v.std()),
STATISTICS_VALID_PERCENT=str(100 * valid.mean()))
def save_results(prefix, depth, valid, profile, meta):
vals = depth[valid]
if not np.isfinite(vals).all():
raise ValueError('Non-finite predictions found')
lo, hi = np.percentile(vals, [1, 99])
if hi <= lo:
raise ValueError('Depth map is effectively constant; do not treat it as a valid result')
# Display only: linear sign reversal, NOT inverse-depth or metric height.
near = np.clip((hi - depth) / (hi - lo), 0, 1).astype(np.float32)
tags = dict(model='DA3MONO-LARGE', units='relative_arbitrary_units',
inference='whole_image', process_res=str(meta['process_res']),
native_prediction_shape=str(meta['native_prediction_shape']))
write_raster(str(prefix) + '_depth.tif', depth, valid, profile,
'Relative depth (larger = farther); NOT height in meters', tags)
write_raster(str(prefix) + '_near01.tif', near, valid, profile,
'Display only: near=1 far=0; clipped P1-P99; NOT height',
dict(tags, units='normalized_display_only',
display_formula='clip((P99-depth)/(P99-P1),0,1)',
display_p1=str(lo), display_p99=str(hi)))
preview = np.zeros(depth.shape, dtype=np.uint8)
preview[valid] = np.rint(near[valid] * 255).astype(np.uint8)
Image.fromarray(preview).save(str(prefix) + '_near.png')
meta.update(depth_statistics={
'min': float(vals.min()), 'max': float(vals.max()),
'p1': float(lo), 'median': float(np.median(vals)), 'p99': float(hi)},
valid_pixels=int(valid.sum()), nodata=NODATA,
depth_units='arbitrary; not meters',
preview='near=white; P1/P99 clipping; invalid=black',
output_resampling='bilinear onto original pixel grid; no new detail added')
Path(str(prefix) + '_metadata.json').write_text(
json.dumps(meta, ensure_ascii=False, indent=2), encoding='utf-8')
print(json.dumps(meta['depth_statistics'], indent=2), flush=True)
def main():
p = argparse.ArgumentParser(description=__doc__)
p.add_argument('--input', required=True)
p.add_argument('--output-prefix', required=True)
p.add_argument('--model', default=str(Path(__file__).resolve().parent / 'models/DA3MONO-LARGE'))
p.add_argument('--process-res', type=int, default=1008,
help='Inference longest side, divisible by 14; try 1008 then 2016')
p.add_argument('--bands', nargs=3, type=int, default=[1, 2, 3], metavar=('R', 'G', 'B'))
p.add_argument('--device', default='cuda:0')
args = p.parse_args()
if args.process_res < 28 or args.process_res % 14:
p.error('--process-res must be >=28 and divisible by 14')
inp = Path(args.input).expanduser().resolve()
prefix = Path(args.output_prefix).expanduser().resolve()
if not Path(args.model).is_dir():
p.error('Local model directory missing; run download_model.py first')
for suffix in ['_depth.tif', '_near01.tif', '_near.png', '_metadata.json']:
if Path(str(prefix) + suffix).exists():
p.error(f'Output exists: {prefix}{suffix}; use a new output prefix')
prefix.parent.mkdir(parents=True, exist_ok=True)
rgb, valid, profile = read_rgb(inp, args.bands)
height, width = valid.shape
if min(height, width) / max(height, width) * args.process_res < 28:
p.error('Image aspect ratio is too extreme at this process resolution')
print(f'Input: {width}x{height}, RGB bands {args.bands}; process-res={args.process_res}', flush=True)
import torch
from depth_anything_3.api import DepthAnything3
if not args.device.startswith('cuda') or not torch.cuda.is_available():
p.error('This workflow requires a CUDA GPU')
torch.cuda.set_device(torch.device(args.device))
print(f'GPU: {torch.cuda.get_device_name()}', flush=True)
model = DepthAnything3.from_pretrained(args.model).to(args.device).eval()
if model.model_name != 'da3mono-large':
raise ValueError(f'This workflow expects da3mono-large, got {model.model_name}')
start = time.time()
with torch.inference_mode():
prediction = model.inference(
[rgb], process_res=args.process_res,
process_res_method='upper_bound_resize',
infer_gs=False, export_dir=None,
)
native = np.asarray(prediction.depth[0], dtype=np.float32)
print('Native depth shape:', native.shape, flush=True)
depth = restore_depth(native, height, width)
meta = dict(input=str(inp), model=str(Path(args.model).resolve()),
model_name=model.model_name, process_res=args.process_res,
native_prediction_shape=list(native.shape), output_shape=[height, width],
crs=str(profile['crs']), transform=list(profile['transform']),
rgb_bands=args.bands, torch_version=torch.__version__,
cuda_runtime=torch.version.cuda, gpu=torch.cuda.get_device_name(),
inference_seconds=time.time() - start)
save_results(prefix, depth, valid, profile, meta)
print('Done:', prefix, flush=True)
if __name__ == '__main__':
main()
da3_tile_io.py
da3_tile_io.py の全文
da3_tile_io.py
"""Raster IO for experimental tiled relative depth mosaics."""
from pathlib import Path
import argparse
import json
import time
import numpy as np
import rasterio
from PIL import Image
NODATA = -9999.0
def read_rgb(path, bands):
with rasterio.open(path) as src:
if min(bands) < 1 or max(bands) > src.count:
raise ValueError(f'Invalid RGB bands {bands}; source has {src.count} bands')
rgb = src.read(bands)
if rgb.dtype != np.uint8:
raise ValueError('Input must be RGB uint8. Convert 16-bit/float data explicitly first.')
if src.crs is None:
raise ValueError('Input has no CRS; supply a georeferenced RGB GeoTIFF')
valid = np.all(src.read_masks(bands) > 0, axis=0)
if not valid.any():
raise ValueError('Input has no valid RGB pixels')
profile = dict(driver='GTiff', width=src.width, height=src.height,
count=1, dtype='float32', crs=src.crs, transform=src.transform,
nodata=NODATA, compress='deflate', predictor=3, tiled=True,
blockxsize=256, blockysize=256, BIGTIFF='IF_SAFER')
rgb = np.moveaxis(rgb, 0, -1).copy()
# Fill masked pixels only for inference; their outputs remain NoData.
if not valid.all():
rgb[~valid] = np.median(rgb[valid], axis=0).astype(np.uint8)
print('Warning: masked areas can affect nearby predictions.', flush=True)
return rgb, valid, profile
def restore_depth(depth, height, width):
depth = np.asarray(depth, dtype=np.float32)
if depth.ndim != 2 or not np.isfinite(depth).all():
raise ValueError('Model depth must be a finite 2-D array')
if depth.shape != (height, width):
depth = np.asarray(Image.fromarray(depth).resize(
(width, height), resample=Image.Resampling.BILINEAR), dtype=np.float32)
return depth.copy()
def write_raster(path, data, valid, profile, description, tags):
arr = np.where(valid, data, NODATA).astype(np.float32)
with rasterio.open(path, 'w', **profile) as dst:
dst.write(arr, 1)
dst.set_band_description(1, description)
dst.update_tags(**tags)
v = arr[valid].astype(np.float64)
dst.update_tags(1, STATISTICS_MINIMUM=str(v.min()),
STATISTICS_MAXIMUM=str(v.max()),
STATISTICS_MEAN=str(v.mean()),
STATISTICS_STDDEV=str(v.std()),
STATISTICS_VALID_PERCENT=str(100 * valid.mean()))
def save_results(prefix, depth, valid, profile, meta):
vals = depth[valid]
if not np.isfinite(vals).all():
raise ValueError('Non-finite predictions found')
lo, hi = np.percentile(vals, [1, 99])
if hi <= lo:
raise ValueError('Depth map is effectively constant; do not treat it as a valid result')
# Display only: linear sign reversal, NOT inverse-depth or metric height.
near = np.clip((hi - depth) / (hi - lo), 0, 1).astype(np.float32)
tags = dict(model='DA3MONO-LARGE', units='relative_arbitrary_units',
inference=meta['inference'], alignment=meta['alignment'], tile_size=str(meta['tile_size']), overlap=str(meta['overlap']), process_res=str(meta['process_res']),
native_prediction_shape=str(meta['native_prediction_shape']))
write_raster(str(prefix) + '_depth.tif', depth, valid, profile,
'Relative depth (larger = farther); NOT height in meters', tags)
write_raster(str(prefix) + '_near01.tif', near, valid, profile,
'Display only: near=1 far=0; clipped P1-P99; NOT height',
dict(tags, units='normalized_display_only',
display_formula='clip((P99-depth)/(P99-P1),0,1)',
display_p1=str(lo), display_p99=str(hi)))
preview = np.zeros(depth.shape, dtype=np.uint8)
preview[valid] = np.rint(near[valid] * 255).astype(np.uint8)
Image.fromarray(preview).save(str(prefix) + '_near.png')
meta.update(depth_statistics={
'min': float(vals.min()), 'max': float(vals.max()),
'p1': float(lo), 'median': float(np.median(vals)), 'p99': float(hi)},
valid_pixels=int(valid.sum()), nodata=NODATA,
depth_units='arbitrary; not meters',
preview='near=white; P1/P99 clipping; invalid=black',
output_resampling='bilinear onto original pixel grid; no new detail added')
Path(str(prefix) + '_metadata.json').write_text(
json.dumps(meta, ensure_ascii=False, indent=2), encoding='utf-8')
print(json.dumps(meta['depth_statistics'], indent=2), flush=True)
infer_geotiff_tiled.py
infer_geotiff_tiled.py の全文
infer_geotiff_tiled.py
"""DA3MONO independent overlapping tile inference; experimental relative-depth mosaic."""
import argparse
from pathlib import Path
import time
import numpy as np
from da3_tile_io import read_rgb, restore_depth, save_results
def starts(length, tile, overlap):
if length <= tile:
return [0]
pos = list(range(0, length - tile + 1, tile - overlap))
if pos[-1] != length - tile:
pos.append(length - tile)
return pos
def fit_overlap(x, y, mode):
"""Robust positive affine fit; reject poorly constrained scale estimates."""
x, y = np.asarray(x, dtype=np.float64), np.asarray(y, dtype=np.float64)
if not len(x):
return 1., 0., 'unanchored'
step = max(1, len(x) // 20000)
x, y = x[::step], y[::step]
shift = float(np.median(y - x))
if mode == 'shift' or len(x) < 128:
return 1., shift, 'shift'
keep = np.ones(len(x), dtype=bool)
a, b = 1., shift
for _ in range(4):
xx, yy = x[keep], y[keep]
if len(xx) < 128 or np.std(xx) < 1e-6 * max(1., abs(np.mean(xx))):
return 1., shift, 'fallback_flat_overlap'
dx, dy = xx - xx.mean(), yy - yy.mean()
corr = np.dot(dx, dy) / max(np.linalg.norm(dx) * np.linalg.norm(dy), 1e-20)
a = float(np.dot(dx, dy) / np.dot(dx, dx))
if corr < .5 or not .25 <= a <= 4.:
return 1., shift, 'fallback_unstable_scale'
b = float(np.median(yy - a * xx))
residual = y - (a * x + b)
med = np.median(residual)
mad = np.median(np.abs(residual - med))
keep = np.abs(residual - med) <= max(3 * 1.4826 * mad, 1e-6)
return a, b, 'affine'
def blend_weight(height, width):
# Positive edge weights retain image boundary pixels.
wy = np.maximum(np.sin(np.pi * (np.arange(height) + .5) / height)**2, .02)
wx = np.maximum(np.sin(np.pi * (np.arange(width) + .5) / width)**2, .02)
return wy[:, None] * wx[None, :]
def mosaic(rgb, valid, predict, tile, overlap, alignment):
h, w = valid.shape
total = np.zeros((h, w), dtype=np.float64)
weights = np.zeros((h, w), dtype=np.float64)
records = []
ys, xs = starts(h, tile, overlap), starts(w, tile, overlap)
for y in ys:
for x in xs:
sl = np.s_[y:min(y+tile, h), x:min(x+tile, w)]
mask = valid[sl]
if not mask.any():
continue
patch = rgb[sl]
ph, pw = mask.shape
# Small images padded; normal edge tiles remain exactly tile x tile.
padded = np.pad(patch, ((0, tile-ph), (0, tile-pw), (0, 0)), mode='edge')
native = predict(padded)
depth = restore_depth(native, tile, tile)[:ph, :pw]
old_w, old_sum = weights[sl], total[sl]
common = mask & (old_w > 0)
a, b, status = 1., 0., 'none'
if alignment != 'none':
a, b, status = fit_overlap(depth[common], old_sum[common]/old_w[common], alignment)
corrected = a * depth + b
weight = blend_weight(ph, pw) * mask
old_sum += corrected * weight
old_w += weight
records.append(dict(x=x, y=y, width=pw, height=ph,
native_shape=list(np.shape(native)), scale=a, shift=b,
alignment_status=status, overlap_pixels=int(common.sum())))
print(f'Tile {len(records)}/{len(ys)*len(xs)} ({x},{y}) {status}: a={a:.4g}, b={b:.4g}', flush=True)
if np.any(valid & (weights <= 0)):
raise RuntimeError('Valid pixels not covered')
output = np.zeros((h, w), dtype=np.float32)
output[valid] = (total[valid] / weights[valid]).astype(np.float32)
return output, records
def main():
p = argparse.ArgumentParser(description=__doc__)
p.add_argument('--input', required=True)
p.add_argument('--output-prefix', required=True)
p.add_argument('--model', default=str(Path(__file__).resolve().parent/'models/DA3MONO-LARGE'))
p.add_argument('--tile-size', type=int, default=300)
p.add_argument('--overlap', type=int, default=100, help='Overlap pixels, NOT stride')
p.add_argument('--process-res', type=int, default=518, help='Model input longest side, multiple of 14')
p.add_argument('--alignment', choices=['affine', 'shift', 'none'], default='affine')
p.add_argument('--bands', nargs=3, type=int, default=[1,2,3])
p.add_argument('--device', default='cuda:0')
args = p.parse_args()
if args.tile_size < 28 or not 0 <= args.overlap < args.tile_size:
p.error('Require tile-size >=28 and 0 <= overlap < tile-size')
if args.alignment != 'none' and args.overlap == 0:
p.error('Alignment needs overlap >0; use --alignment none for nonoverlapping tiles')
if args.process_res < 28 or args.process_res % 14:
p.error('process-res must be >=28 and divisible by 14')
if not Path(args.model).is_dir():
p.error('Model directory missing; specify existing --model path')
prefix = Path(args.output_prefix).expanduser().resolve()
for suffix in ['_depth.tif','_near01.tif','_near.png','_metadata.json']:
if Path(str(prefix)+suffix).exists():
p.error(f'Output exists: {prefix}{suffix}; use a new prefix')
rgb, valid, profile = read_rgb(args.input, args.bands)
import torch
from depth_anything_3.api import DepthAnything3
if not args.device.startswith('cuda') or not torch.cuda.is_available():
p.error('CUDA GPU required')
torch.cuda.set_device(torch.device(args.device))
model = DepthAnything3.from_pretrained(args.model).to(args.device).eval()
if model.model_name != 'da3mono-large':
p.error('Expected DA3MONO-LARGE')
def predict(patch):
with torch.inference_mode():
result = model.inference([patch], process_res=args.process_res,
process_res_method='upper_bound_resize',
infer_gs=False, export_dir=None)
return np.asarray(result.depth[0], dtype=np.float32)
start = time.time()
depth, records = mosaic(rgb, valid, predict, args.tile_size, args.overlap, args.alignment)
meta = dict(input=str(Path(args.input).resolve()), model=str(Path(args.model).resolve()),
model_name=model.model_name, process_res=args.process_res,
native_prediction_shape=records[0]['native_shape'], output_shape=list(valid.shape),
crs=str(profile['crs']), transform=list(profile['transform']), rgb_bands=args.bands,
inference='independent_overlapping_tiles', tile_size=args.tile_size,
overlap=args.overlap, stride=args.tile_size-args.overlap, alignment=args.alignment,
alignment_warning='Experimental relative mosaic; scan-order drift and disconnected scales possible; not metric depth or height',
tiles=records, torch_version=torch.__version__, cuda_runtime=torch.version.cuda,
gpu=torch.cuda.get_device_name(), inference_seconds=time.time()-start)
prefix.parent.mkdir(parents=True, exist_ok=True)
save_results(prefix, depth, valid, profile, meta)
print('Done:', prefix, flush=True)
if __name__ == '__main__':
main()