有吉北貝塚北斜面貝層 土器分布平面図・立面図から3D座標を取得した方法

 North Slope Shell Layer of Ariyoshi Kita Shell Mound

Method of obtaining 3D coordinates from pottery distribution plan/elevation

Make a note of how the 3D coordinates were obtained from the pottery distribution plan and elevation shown in the excavation report. I was able to improve work efficiency by using Python scripts extensively.

発掘調査報告書に掲載されている土器分布平面図・立面図から3D座標を取得した方法をメモします。Pythonスクリプトを多用することにより作業効率化を図ることができました。

1　3D座標取得方法

1-1　土器分布平面図・立面図の確認及び表示土器の把握

1-1-1　土器分布平面図・立面図の確認

発掘調査報告書には平面図5枚、立面図5枚が掲載されています。

土器分布平面図・立面図をBlender3Dに配置した様子

【参考】

2022年にはBlender3D空間に配置された平面図・立面図について、クリックするとその場所にCUBEを生成するBlenderPythonスクリプトを作成しました。このスクリプトを利用して、対応する同一土器破片のCUBEを平面図と立面図で作成し、その座標値を利用して3D空間にCUBEを配置する作業を行いました。Blenderではこの手の作業にインターフェイスが対応していないためとても面倒な作業となりました。個別土器単位の作業はできなくて、平面図における分類（≒型式）単位の作業となりました。

1-1-2　平面図・立面図別表示土器の把握

平面図・立面図5枚それぞれについて、Excelで表示土器番号リストをY座標2ｍピッチ毎に作成し、土器番号がわかれば、図面のどの部分にそれが表示されているのか迅速に判断できるようにしました。

表示土器番号リスト（Excel画面、一部）

今回3D座標計測する土器数は全部で173です。破片数は結果的に643となりました。

1-2　平面図・立面図から座標計測するPythonスクリプト作成

平面図5枚、立面図5枚について土器出土地点をクリックするとBlender3Dの座標を取得できるPythonスクリプトを作成しました。

平面図h3用の座標計測Pythonスクリプト

import cv2
def onMouse(event, x, y, flags, params):
    if event == cv2.EVENT_LBUTTONDOWN:
        print('{:.03f}'.format((y-97)/95.111+2), '{:.03f}'.format((x-36)/94.125+20), 0)
img = cv2.imread('h3e.png')
cv2.namedWindow('gazo', cv2.WINDOW_NORMAL)
cv2.imshow('gazo', img)
cv2.setMouseCallback('gazo', onMouse)
cv2.waitKey(0)

h3.pyをコマンドラインで走らせると計測用画像（平面図）が画面に出て、計測ポイントをクリックするとblender3D空間対応座標（x座標、y座標）が画面に表示されます。

h3.pyを走らせた様子（実際の作業は別画面に平面図を大伸ばしして行います。）

1-3　作業フォルダー作成

3D座標計測作業は173土器別フォルダーをパソコン内に作成して行いました。173フォルダーの作成、各ホルダーに必要な図面とPythonスクリプトのコピー、作業を記録するテキストファイル作成などはすべてPythonスクリプトを自前で作成して、効率的に行いました。

173フォルダーにフォルダー名称を名称とするテキストファイルを生成するPythonスクリプト

import os

# スクリプトが実行されているフォルダーのパスを取得
current_folder = os.path.dirname(os.path.abspath(__file__))

# input.txtのパス
input_file = os.path.join(current_folder, "inputall.txt")

# フォルダーのパスとファイル名の対応辞書を作成
folder_file_mapping = {}

# input.txtを読み込んで対応辞書を作成
with open(input_file, "r") as file:
    for line in file:
        folder_name = line.strip()
        folder_path = os.path.join(current_folder, folder_name)
        file_name = folder_name + ".txt"
        folder_file_mapping[folder_path] = file_name

# フォルダーごとにテキストファイルを生成
for folder_path, file_name in folder_file_mapping.items():
    file_path = os.path.join(folder_path, file_name)
    with open(file_path, "w") as file:
        # 何かしらのテキストをファイルに書き込む（ここではフォルダー名を書き込んでいます）
        file.write(folder_name)

1-4　3D座標計測と記録

平面図の座標計測用Pythonスクリプトを起動してx座標とy座標を、立面図からはy座標とz座標を読み取り、y座標の対応性を利用して各個別土器破片の3D座標を取得記録しました。この作業はコマンドラインとエディター（テキストファイル）で行います。

3D座標のテキスト記載

1-5　3D座標の集成とBlenderプロット

3D座標を集成したテキストを作成し、自作PythonスクリプトによりBlenderにプロットします。3D座標には土器番号がありますから、土器型式に関する変更とか、別の分類などをテキストで編集すれば、それにより小分けした（分類した）プロットが可能です。

3D座標をBlenderにプロットするBlenderPythonスクリプト

import bpy

# テキストファイルのパス
file_path = "I:\\20221214有吉北貝塚北斜面貝層プロジェクト\\■20230520土器破片3D座標\\273～335　E2新\\e2sin.txt"

# プリミティブな立方体のスケール
cube_scale = (0.15, 0.15, 0.15)

# テキストファイルの読み込み
with open(file_path, 'r') as file:
    lines = file.readlines()

# 行ごとに処理を行う
for line in lines:
    # 行をスペースで分割して各値を取得
    symbol, x, y, z = line.strip().split(",")

    # 座標値を数値に変換
    x = float(x)
    y = float(y)
    z = float(z)
    
    # 立方体を生成
    bpy.ops.mesh.primitive_cube_add(scale=cube_scale, location=(x, y, z))
    cube = bpy.context.object
    
    # 立方体の名前を設定（記号を含む）
    cube.name = symbol

2　【参考】3D座標から構成した分布図

3D座標から構成した分布図

3　感想

・ChatGPTと対話して必要なPythonスクリプトを生成する活動を5月から始めましたが、そうした世の中の画期との連動がなければ今回の土器3D座標取得は半年とか１年は遅れた可能性があります。

・発掘調査報告書の刊行は1998年で北斜面貝層現場調査は1985年です。40年前から25年前までの活動で土器出土位置を3D座標で把握できるような精細資料を作成したことに感動します。当時の関係者はどのような見通しを持って土器分布平面図・立面図を作成したのでしょうか。とてつもない労力を投入しています。