NumPyで点群データから鉄筋の軸を自動計測する

この記事でやること

壁から突き出したD51鉄筋を 3D スキャンした点群（LAS ファイル）を入力にして、

• 鉄筋を 1 本ずつ自動で切り分け
• 1 本ごとに 軸の 3 次元ベクトルと 根元・先端の直径を求め
• 壁面の法線を基準にした 水平のずれ / 鉛直のずれ / なす角を算出し
• 結果を CSV と AutoCAD の DXF に書き出す

という処理を、NumPy と標準ライブラリのみで書いた話です。
（D51は鉄筋の直径が51mmという意味です。現場ではデコイチなんて言います）

コード全体はこちら

import numpy as np
import struct
import csv

tip_path = r"example\D51.las"
full_path = r"example\D51.las"
out_dir = r"example"


def find_edges(values,gap):
    values_sorted = np.sort(values)
    boundaries = []
    for i in range(len(values_sorted) - 1):
        if values_sorted[i + 1] - values_sorted[i] > gap:
            boundaries.append((values_sorted[i] + values_sorted[i + 1]) / 2)
    return [-np.inf] + boundaries + [np.inf]

def read_las_xyz(path):
    # 1. ファイルを開いて、ヘッダ(先頭227バイト)を読む
    with open(path, "rb") as f:
        data = f.read(227)

    # 2. 点の表の開始位置を取り出す
    point_offset = struct.unpack("<I", data[96:100])[0]

    # 3. 1点のバイト数、点の数を取り出す
    record_length = struct.unpack("<H", data[105:107])[0]
    n_points = struct.unpack("<I", data[107:111])[0]

    # 4. 倍率とオフセットを取り出す
    scale = struct.unpack("<3d", data[131:155])
    offset = struct.unpack("<3d", data[155:179])

    # 5. numpy で点の表を読む (1点 = Xの整数, Yの整数, Zの整数, 残り14バイト)
    point_type = np.dtype([("X", "<i4"), ("Y", "<i4"), ("Z", "<i4"), ("rest", f"V{record_length - 12}")])
    points = np.fromfile(path, dtype=point_type, count=n_points, offset=point_offset)

    # 6. 整数を実座標 (m) に直す
    x = points["X"] * scale[0] + offset[0]
    y = points["Y"] * scale[1] + offset[1]
    z = points["Z"] * scale[2] + offset[2]

    # 7. X, Y, Z を1つの表 (点の数 × 3列) にまとめる
    points_array = np.column_stack((x, y, z))
    
    return points_array

def along_and_radial(points, point, direction):
    w = points - point
    s = w @ direction                       # 内積: 各点が軸方向にどれだけ進んだ位置か
    radial = w - np.outer(s, direction)     # 軸方向の成分を引いた残り
    return s, np.linalg.norm(radial, axis=1)

def find_face(s ,r):
    ring = (r > 0.05) & (r < 0.20)
    counts, edges = np.histogram(s[ring], bins=np.arange(-0.6, 0.6, 0.005))
    return edges[np.argmax(counts)]  # 一番点が多い区間 = 壁

def fit_circle(u, v):
    A = np.column_stack((2 * u, 2 * v, np.ones(len(u))))
    a, b, k = np.linalg.lstsq(A, u**2 + v**2, rcond=None)[0]
    return a, b, np.sqrt(k + a**2 + b**2)

def fit_circle_trim(u, v):
    use = np.ones(len(u), dtype=bool)          # 最初は全部の点を使う
    for _ in range(30):
        a, b, radius = fit_circle(u[use], v[use])
        new_use = np.abs(np.hypot(u - a, v - b) - radius) <= 0.005   # 円周から5mm以内の点
        if (new_use == use).all():             # 使う点が変わらなくなったら終わり
            break
        use = new_use
    return a, b, radius, (~use).sum()

# AutoCADに書き出す関数
def write_dxf_lines(path, starts, ends, layer="REBAR_AXIS"):
    with open(path, "w", encoding="ascii") as f:
        f.write("0\nSECTION\n2\nHEADER\n9\n$ACADVER\n1\nAC1009\n0\nENDSEC\n")
        f.write("0\nSECTION\n2\nENTITIES\n")
        for k in range(len(starts)):
            p = starts[k]
            q = ends[k]
            f.write(f"0\nLINE\n8\n{layer}\n")
            f.write(f"10\n{p[0]:.6f}\n20\n{p[1]:.6f}\n30\n{p[2]:.6f}\n")
            f.write(f"11\n{q[0]:.6f}\n21\n{q[1]:.6f}\n31\n{q[2]:.6f}\n")
        f.write("0\nENDSEC\n0\nEOF\n")

points_array = read_las_xyz(tip_path)

# 9. 確かめる: 点の数、最初の点、最小・最大
print(f"点の数: {points_array.shape[0]}")
print(f"最初の点: {points_array[0]}")
print(f"最小値: {points_array.min(axis=0)}")
print(f"最大値: {points_array.max(axis=0)}")
print(points_array)

print(find_edges(points_array[:, 2], 0.10))

bar_names = []  # 名前を入れる空のリスト
bar_points = [] # 点の表を入れる空のリスト

# (1) 段を先に全部集める
z_edges = find_edges(points_array[:, 2], 0.10)
row_groups = []
for i in range(len(z_edges) - 1):
    row_mask = (points_array[:, 2] >= z_edges[i]) & (points_array[:, 2] < z_edges[i + 1])
    row_points = points_array[row_mask]
    if len(row_points) >= 500:
        row_groups.append(row_points)

# (2) 段の数で名前を決める (Z の低い順)
if len(row_groups) == 2:
    row_names = ["下段", "上段"]
elif len(row_groups) == 3:
    row_names = ["下段", "中段", "上段"]
else:
    row_names = []
    for i in range(len(row_groups)):
        row_names.append(f"{i + 1}段目")        # 4段以上なら、下から「1段目」「2段目」…

# (3) 段ごとに X で分ける
for row_index in range(len(row_groups)):
    row_points = row_groups[row_index]
    row_name = row_names[row_index]
    col_index = 0
    x_edges = find_edges(row_points[:, 0], 0.10)
    for j in range(len(x_edges) - 1):
        bar_mask = (row_points[:, 0] >= x_edges[j]) & (row_points[:, 0] < x_edges[j + 1])
        bar = row_points[bar_mask]
        if len(bar) < 500:
            continue
        col_index = col_index + 1
        name = f"{row_name}_{col_index}"
        bar_names.append(name)
        bar_points.append(bar)
        print(f"{name}: {len(bar)} 点, 鉄筋重心: {bar.mean(axis=0).round(3)}")

print(f"本数: {len(bar_names)}")

bar_centers = []
bar_directions = []
for k in range(len(bar_points)):
    bar = bar_points[k]
    center = bar.mean(axis=0)
    u, s, vt = np.linalg.svd(bar - center, full_matrices=False)
    direction = vt[0]
    bar_centers.append(center)
    bar_directions.append(direction)

full_points = read_las_xyz(full_path)
print(f"全体: {len(full_points)} 点")
print(f"最小 {full_points.min(axis=0).round(3)}")
print(f"最大 {full_points.max(axis=0).round(3)}")

root_centers = []
tip_centers = []
axes = []            # 1本ごとの軸の向き
root_diameters = []
tip_diameters = []
wall_points = []     # 鉄筋のまわりの壁の点

for k in range(len(bar_names)):
    point = bar_centers[k]
    direction = bar_directions[k]
    print(f"基準点 {point.round(4)}, 向き {direction.round(4)}")

    near = full_points[(np.abs(full_points - point) < 0.7).all(axis=1)]
    print(f"箱の中: {len(near)} 点")

    s, r = along_and_radial(near, point, direction)
    tube = r < 0.08
    print(f"筒の中: {tube.sum()} 点, s の範囲: {s[tube].min():.3f} 〜 {s[tube].max():.3f}")

    face = find_face(s, r)
    if face > 0:
        direction = -direction
        s = -s
        face = -face
    end = s[tube].max()
    print(f"壁 s: {face:.3f}, 先端 s = [end={end:.3f}]")

    body = tube & (s >face + 0.05) & (s < end - 0.02) #壁際
    rough_point = near[body].mean(axis=0) #胴体の点の重心
    u, sv ,vt = np.linalg.svd(near[body] - rough_point, full_matrices=False) # 主成分分析
    new_direction = vt[0] # 主成分の向き
    if new_direction @ direction < 0: # 向きが逆なら反転
        new_direction = -new_direction
    direction = new_direction

    s,r = along_and_radial(near, rough_point, direction)
    face = find_face(s,r)
    wall_ring = (r > 0.05) & (r < 0.20) & (np.abs(s - face) < 0.015)   # 鉄筋のまわりの壁の点
    wall_points.append(near[wall_ring])
    bar_mask = (s > face +0.01) & (r < 0.045)
    bar_full = near[bar_mask]
    end = s[bar_mask].max()
    print(f"鉄筋の点: {len(bar_full)} 点, 露出長さ: {1000 * (end - face):.0f} mm")

    e2 = np.cross(direction, [0, 0, 1])     # 横 (軸にも鉛直にも直角 = 水平で壁に沿った向き)
    e2 = e2 / np.linalg.norm(e2)            # 長さを1にする
    e3 = np.cross(e2, direction)            # 上 (外積で軸にも e2 にも直角)

    centers = []
    diameters = []

    for label, s0 in [("根元", face + 0.03), ("先端", end - 0.035)]:
        p0 = rough_point + s0 * direction           # 輪切りの中心 (大まかな軸上の点)
        w = bar_full - p0
        in_slice = np.abs(w @ direction) <= 0.02    # 軸方向に ±20mm = 厚さ4cm
        u = w[in_slice] @ e2                        # 断面上の横の座標
        v = w[in_slice] @ e3                        # 断面上の上の座標
        a, b, radius, n_removed = fit_circle_trim(u, v)
        print(f"{label}: {in_slice.sum()} 点, 中心 u = {1000 * a:.1f} mm, v = {1000 * b:.1f} mm, 直径 {2000 * radius:.1f} mm")
        center = p0 + a * e2 + b * e3       # 断面の (a, b) を3次元の座標に戻す
        print(f"3次元の中心: {center.round(4)}")
        centers.append(center)
        diameters.append(2 * radius)
    root_center = centers[0]
    tip_center = centers[1]
    axis = (tip_center - root_center) / np.linalg.norm(tip_center - root_center)
    axes.append(axis)
    root_centers.append(root_center)
    tip_centers.append(tip_center)
    root_diameters.append(diameters[0])
    tip_diameters.append(diameters[1])

# 壁の面の向き (法線) を求める
all_wall = np.vstack(wall_points)                 # 全部の鉄筋のまわりの壁の点を1つの表に積み重ねる
wall_center = all_wall.mean(axis=0)
u, sv, vt = np.linalg.svd(all_wall - wall_center, full_matrices=False)
wall_normal = vt[2]                               # 一番広がっていない方向 = 壁の法線
if wall_normal @ axes[0] < 0:                     # 鉄筋と同じ向き (壁から外向き) にそろえる
    wall_normal = -wall_normal
wall_horizontal = np.degrees(np.arctan2(wall_normal[1], wall_normal[0]))
wall_slope = np.degrees(np.arcsin(wall_normal[2]))
print(f"壁の法線: 水平角 {wall_horizontal:.2f}°, 勾配 {wall_slope:+.2f}°")

# 法線を基準にした角度を出して、表の行を作る
rows = []
for k in range(len(bar_names)):
    axis = axes[k]
    horizontal_dev = np.degrees(np.arctan2(axis[1], axis[0])) - wall_horizontal   # 水平のずれ (+ = 壁側から先端を見て左)
    vertical_dev = np.degrees(np.arcsin(axis[2])) - wall_slope                    # 鉛直のずれ (+ = 上向き)
    total_dev = np.degrees(np.arccos(min(1, axis @ wall_normal)))                 # 法線とのなす角
    print(f"{bar_names[k]}: 水平 {horizontal_dev:+.2f}°, 鉛直 {vertical_dev:+.2f}°, なす角 {total_dev:.2f}°")
    rows.append([bar_names[k],
                 *root_centers[k].round(4), *tip_centers[k].round(4),
                 round(horizontal_dev, 2), round(vertical_dev, 2), round(total_dev, 2),
                 round(1000 * root_diameters[k], 1), round(1000 * tip_diameters[k], 1)])

# CSV に書き出す
header = ["鉄筋", "根元X", "根元Y", "根元Z", "先端X", "先端Y", "先端Z",
          "水平のずれ°", "鉛直のずれ°", "法線とのなす角°", "根元直径mm", "先端直径mm"]
with open(out_dir + r"\鉄筋軸_結果.csv", "w", newline="", encoding="utf-8-sig") as f:
    writer = csv.writer(f)
    writer.writerow(header)      # 見出しの1行
    writer.writerows(rows)       # 8本分の行をまとめて

# DXF に書き出す
write_dxf_lines(out_dir + r"\鉄筋軸_2円法.dxf", root_centers, tip_centers)

点群を使って何か測りたいけど、ライブラリの中で何が起きているのか分からないのが気持ち悪いという人向けです。というか、私がそういう人です。
処理はすべて工学系大学生の線形代数で書けるレベルです。

動作環境
  Python 3.10+
  numpy          （外部依存）
  struct / csv   （標準ライブラリ）

---

1. 背景：なぜ「鉄筋の軸」を測りたいのか

継手や定着のために壁から鉄筋を突き出しておく場面で、その鉄筋がどれくらい狙いの向きからズレて建て込まれているかを知りたい、という依頼が現場からありました。

従来は下げ振りで 1 本ずつ測ったり、A3用紙を当てて印をつけていたところを、スキャナで一度に撮った点群から自動で出せないか、という話になります。

ずれを見るので、以下の3つの量について調べる必要があります。

量	意味
水平のずれ（°）	壁面の法線から見て、左右にどれだけ振れているか
鉛直のずれ（°）	同じく、上下にどれだけ振れているか
直径（mm）	根元と先端で、ちゃんと D51 相当（公称径 50.8mm）の径が出ているか

ポイントは 絶対座標系での角度ではなく壁面基準での角度 が欲しいこと。スキャナの座標系がどう傾いていようと関係なく評価したいので、壁の法線も点群から求めて、そこからの相対角を出します。

---

2. 全体の流れ

入力に必要なのは 2 つの LAS ファイルです。

• tip … コンクリート壁を含まない点群（本数の数え上げ・初期位置決め用）
• full … 供試体全体の点群（コンクリート壁も鉄筋も入っている、数値解析用）

先端だけの点群を使うのは、本数のカウントと初期位置の推定を簡単にするためです。
全体点群からいきなりクラスタリングすると壁や背景に邪魔されてしまいます。
鉄筋だけなら空間がきれいに孤立しているので、隣との距離による分割で確実に分かれてくれます。

[鉄筋のみ LAS]
      │
      ├─ ① Z でギャップ分割 → 段（上段/中段/下段）
      ├─ ② 段ごとに X でギャップ分割 → 1 本ずつ
      └─ ③ 本ごとに SVD → 粗い重心と粗い軸
                │
[全体 LAS] ─────┤
      ├─ ④ 粗い軸まわりの筒で鉄筋を切り出す
      ├─ ⑤ ヒストグラムで壁面の位置を検出
      ├─ ⑥ 壁際と先端を捨てて、再度 SVD → 精密な軸
      ├─ ⑦ 根元・先端で輪切り → 最小二乗円フィット（二円法）
      │        └→ 2 つの円の中心を結んだ線が「真の軸」
      ├─ ⑧ 全鉄筋まわりの壁点から SVD → 壁の法線
      └─ ⑨ 法線基準の角度を計算
                │
        [CSV] と [DXF]

---

3. ステップ①：LAS を NumPy だけで読む

laspy を入れればそれで済む話です。が、勉強用として...

以前も記事にしましたが、LAS の点データ部は固定長レコードが並んでいるだけなので、np.fromfile の構造化 dtype で読めます。依存を減らしたいのと、読み込みが速いのとで、自前で書きました。

def read_las_xyz(path):
    # 1. ファイルを開いて、ヘッダ(先頭227バイト)を読む
    with open(path, "rb") as f:
        data = f.read(227)

    # 2. 点の表の開始位置を取り出す
    point_offset = struct.unpack("<I", data[96:100])[0]

    # 3. 1点のバイト数、点の数を取り出す
    record_length = struct.unpack("<H", data[105:107])[0]
    n_points = struct.unpack("<I", data[107:111])[0]

    # 4. 倍率とオフセットを取り出す
    scale = struct.unpack("<3d", data[131:155])
    offset = struct.unpack("<3d", data[155:179])

    # 5. numpy で点の表を読む (1点 = Xの整数, Yの整数, Zの整数, 残り14バイト)
    point_type = np.dtype([("X", "<i4"), ("Y", "<i4"), ("Z", "<i4"),
                           ("rest", f"V{record_length - 12}")])
    points = np.fromfile(path, dtype=point_type, count=n_points, offset=point_offset)

    # 6. 整数を実座標 (m) に直す
    x = points["X"] * scale[0] + offset[0]
    y = points["Y"] * scale[1] + offset[1]
    z = points["Z"] * scale[2] + offset[2]

    return np.column_stack((x, y, z))

ヘッダのバイト位置の意味

マジックナンバーだらけに見えますが、LAS 1.2 のパブリックヘッダブロック（227 バイト）の仕様に則った結果こうなってるだけです。

オフセット	型	内容
96–99	uint32	Offset to point data（点データの開始位置）
105–106	uint16	Point Data Record Length（1 点あたりのバイト数）
107–110	uint32	Legacy Number of point records（点数）
131–154	double×3	X/Y/Z Scale Factor
155–178	double×3	X/Y/Z Offset

("rest", f"V{record_length - 12}")

について、LAS の 1 点は、先頭 12 バイトが int32 の X/Y/Z で、そのあとに Intensity、Return Number、分類コード、GPS 時刻、RGB……と点フォーマットごとに座標値以外のデータが続きます。今回は座標だけあればいいので void型で残り全部を名前だけ付けて無視すれば、どの点フォーマットでも同じコードで座標値だけ取得できます。

座標が整数で入っている理由

LAS は座標を int32 で持っていて、実座標は

実座標 [m] = 格納された整数 × scale + offset

で復元します。これは浮動小数の丸め誤差を避けつつファイルサイズを抑えるための仕様で、scale = 0.001 なら 1mm 刻みで格納されている、という意味になります。

LAS 1.4 ではヘッダが 375 バイトに拡張され、点フォーマット 6 以降を使う場合は 107 バイト目の「レガシー点数」が 0 になります（本当の点数は 247 バイト目の uint64）。この実装は 1.2 系（点フォーマット 0〜5）前提なので、1.4 を食わせると点数 0 。汎用にするなら data[24:26] のバージョンを見て分岐が必要です。

---

4. ステップ②：鉄筋を 1 本ずつに切り分ける

ソートして、隙間が空いているところで切ります

def find_edges(values, gap):
    values_sorted = np.sort(values)
    boundaries = []
    for i in range(len(values_sorted) - 1):
        if values_sorted[i + 1] - values_sorted[i] > gap:
            boundaries.append((values_sorted[i] + values_sorted[i + 1]) / 2)
    return [-np.inf] + boundaries + [np.inf]

値を昇順に並べ、隣どうしの差が GAP（= 0.10 m）より大きいところを「グループの境目」とみなして、その中点を返します。前後に -inf / +inf を足しているので、返り値をそのまま区間の端点として使えます。

これを Z方向とX方向で適用します。

# (1) Z でグループ分け → 段
z_edges = find_edges(points_array[:, 2], GAP)

# (3) 段ごとに X でグループ分け → その段の中の 1 本ずつ
x_edges = find_edges(row_points[:, 0], GAP)

高さ（Z）で上段・中段・下段に割り、その中で横方向（X）に並んだ鉄筋を 1 本ずつに割るという、まぁ人が図面を読むときと同じ順番ですね。

段数に応じて名前も自動で決めています。

if len(row_groups) == 2:
    row_names = ["下段", "上段"]
elif len(row_groups) == 3:
    row_names = ["下段", "中段", "上段"]
else:
    row_names = [f"{i + 1}段目" for i in range(len(row_groups))]

MIN_POINTS = 500 未満のグループはノイズとして無視するようにします。
スキャンした点群には飛び点が必ず混ざるので、この 1 行がないと点 3 個の鉄筋があるぞ！みたいなことになります。

---

5. ステップ③：とりあえずの鉄筋軸をつくる

グループ分けできたら、1 本ごとに主成分分析（PCA）で軸方向を出します。

center = bar.mean(axis=0)
u, s, vt = np.linalg.svd(bar - center, full_matrices=False)
direction = vt[0]

重心を引いた行列を SVD すると、vt の各行が主成分方向（分散の大きい順）になります。鉄筋は明らかに細長い形なので、vt[0] ＝ 一番伸びている方向 ＝ 軸方向です。

あとで vt[2]（一番広がっていない方向）＝ 平面の法線 としても再登場します。

---

6. ステップ④：軸まわりの「筒」で鉄筋を切り出す

ここから全体点群を使用します。まず、各点が軸に沿ってどこにいて、軸からどれだけ離れているかを出す関数を用意します。

def along_and_radial(points, point, direction):
    w = points - point
    s = w @ direction                       # 内積: 軸方向にどれだけ進んだ位置か
    radial = w - np.outer(s, direction)     # 軸方向の成分を引いた残り
    return s, np.linalg.norm(radial, axis=1)

要するに 直交座標系から円柱座標系への変換です。

• s … 軸方向の位置（スカラー、符号つき）
• r … 軸からの距離（＝円柱座標の半径）

w - np.outer(s, direction) が軸に垂直な平面への射影になっています。

near = full_points[(np.abs(full_points - point) < BOX_HALF).all(axis=1)]  # まず箱で粗く絞る
s, r = along_and_radial(near, point, direction)
tube  = r < TUBE_ROUGH                                          # 半径 80mm の筒
body  = tube & (s > face + BODY_ROOT_CUT) & (s < end - BODY_TIP_CUT)  # 壁際と先端を除いた胴体

先に軸平行の箱（BOX_HALF = 0.7 m）で絞ってから筒の判定をしているのは、計算量削減のためです。

---

7. ステップ⑤：壁面をヒストグラムで見つける

壁の位置（軸方向の座標 s）が分かると、根元から何 mmという基準が作れます。

def find_face(s, r):
    ring = (r > WALL_RING_IN) & (r < WALL_RING_OUT)
    counts, edges = np.histogram(s[ring], bins=np.arange(-WALL_SEARCH, WALL_SEARCH, WALL_BIN))
    return edges[np.argmax(counts)]  # 一番点が多い区間 = 壁

• 軸から 50mm〜200mm のドーナツ状の領域（ring）だけを見る
  → 鉄筋そのもの（半径 ~27mm）は入らず、壁面だけが入る
• その領域の点を s でヒストグラムにする（5mm 刻み）
• 壁は軸にほぼ垂直な平面なので、壁面の s に点が極端に集中する
• よって 最頻ビン ＝ 壁面の位置

まぁ要はs方向で点群がいっぱいあるところって壁だよねっていう話です

ついでに、軸の「向き」をそろえる

SVD が返す主成分方向は符号が+d でも -d でも主成分なので、ここで向きを固定します。

face = find_face(s, r)
if face > 0:
    direction = -direction
    s = -s
    face = -face

---

8. ステップ⑥：軸を取り直す

最初の軸は鉄筋だけの点群から出した粗い値なので、全体点群を使って求め直します。

body = tube & (s > face + BODY_ROOT_CUT) & (s < end - BODY_TIP_CUT)  # 胴体だけ
rough_point = near[body].mean(axis=0)
u, sv, vt = np.linalg.svd(near[body] - rough_point, full_matrices=False)
new_direction = vt[0]
if new_direction @ direction < 0:   # 向きが逆なら反転
    new_direction = -new_direction
direction = new_direction

壁際 50mm（BODY_ROOT_CUT）と先端 20mm（BODY_TIP_CUT）を捨てているのは、

• 壁際：コンクリート面の凹凸、鉄筋と壁の隅に溜まったノイズ点が入るから。
• 先端：切断面の点は軸方向に「面」として広がるので、軸が計算上ずれる。

細長い物体の PCA は、両端の外れ値に弱いという性質への対処です。

さらに新しい軸で find_face をもう一度やり直し、壁位置も更新します。軸 → 壁 → 軸 → 壁 と 2 周まわして収束させています。

---

9. ステップ⑦：鉄筋の軸の定義

PCA で出した軸をそのまま使わないのは異形鉄筋には節（リブ）があるので、点群の分布が完全な円柱にならないからです。節の出方が場所によって偏ると、PCA の主軸がちょっとずれるんですね。
（建設業では鉄筋とコンクリートの付着面積を広げるためにツルツルの円柱の鉄筋ではなく、デコボコした異形鉄筋というものを使います。）

そこで採ったのが、
根元付近と先端付近で 1 枚ずつ「輪切り」を作り、それぞれ最小二乗で円をフィットして中心を求め、その 2 つの中心を結んだ直線を、鉄筋の軸とする手法です。

断面の中心は節の影響を受けにくいので、PCA の主軸よりも素直な値になるんじゃないかと。

9-1. 断面の座標系を作る

3 次元の点を「断面上の 2 次元座標」に落とすために、軸に垂直な 2 本の単位ベクトルを作ります。

e2 = np.cross(direction, [0, 0, 1])     # 横 (軸にも鉛直にも直角 = 水平で壁に沿った向き)
e2 = e2 / np.linalg.norm(e2)            # 長さを1にする
e3 = np.cross(e2, direction)            # 上 (外積で軸にも e2 にも直角)

direction と鉛直ベクトルの外積 → 水平かつ軸に垂直な e2。さらに e2 と direction の外積 → e3。これで (direction, e2, e3) が正規直交基底になりますね。

e2 を「水平」にそろえてあるので、あとで出てくる u は水平方向のずれ、v は鉛直方向のずれとして解釈できます。

9-2. 輪切りを作る

for label, s0 in [("根元", face + ROOT_SLICE_CENTER), ("先端", end - TIP_SLICE_CENTER)]:
    p0 = rough_point + s0 * direction                          # 輪切りの中心
    w = bar_full - p0
    in_slice = np.abs(w @ direction) <= SLICE_THICKNESS / 2    # 軸方向に ±厚さ/2
    u = w[in_slice] @ e2                                       # 断面上の横の座標
    v = w[in_slice] @ e3                                       # 断面上の上の座標

定数	値	根拠
SLICE_THICKNESS	40 mm	D51 の節の間隔（最大 35.6mm = 公称直径の 0.7 倍）より厚く取ります。こうすると、どの位置で切っても輪切りの中に節がちょうど 1 周期分入るので、節の有無による径のばらつきが平均化される
ROOT_SLICE_CENTER	壁から 30 mm	壁の凹凸の影響を避けつつ、できるだけ根元で測る（厚さ 40mm なので実質 10〜50mm の範囲）
TIP_SLICE_CENTER	先端から 35 mm	切断面の乱れを避ける（実質 15〜55mm の範囲）

9-3. 最小二乗で円をあてる

円のフィッティングは、そのままでは非線形問題（未知数 a, b, R に対して (u-a)² + (v-b)² = R²）ですが、変数変換で線形最小二乗に落とせます。

展開すると：

u² - 2au + a² + v² - 2bv + b² = R²
u² + v² = 2au + 2bv + (R² - a² - b²)

ここで k = R² - a² - b² と置けば、未知数 (a, b, k) について線形です。

def fit_circle(u, v):
    A = np.column_stack((2 * u, 2 * v, np.ones(len(u))))
    a, b, k = np.linalg.lstsq(A, u**2 + v**2, rcond=None)[0]
    return a, b, np.sqrt(k + a**2 + b**2)

A @ [a, b, k] = u² + v² を lstsq で解いて、最後に R = √(k + a² + b²) で半径を復元します。

Kåsa 法といいます。点が円周全体に散らばっていれば十分実用的で、今回のように鉄筋の全周が見えているケースでは問題になりません。

9-4. 外れ点を反復で除く（トリミング）

素の最小二乗は外れ値に弱いので、「フィット → 外れ点を外す → 再フィット」を収束するまで繰り返します。

def fit_circle_trim(u, v):
    use = np.ones(len(u), dtype=bool)          # 最初は全部の点を使う
    for _ in range(TRIM_MAX_ITER):
        a, b, radius = fit_circle(u[use], v[use])
        new_use = np.abs(np.hypot(u - a, v - b) - radius) <= TRIM  # 円周から TRIM 以内
        if (new_use == use).all():             # 使う点が変わらなくなったら終わり
            break
        use = new_use
    return a, b, radius, (~use).sum()

ポイントは、毎回すべての点に対して判定し直していることです。一度外した点も、フィットが改善した結果また円周の近くに来れば復帰できます。

TRIM = 5mm異形鉄筋の節の高さ（D51 なら JIS 上 2.5〜5.1mm）は残し、明らかな飛び点だけを落とす閾値です。

9-5. 2 つの中心から軸を作る

center = p0 + a * e2 + b * e3       # 断面の (a, b) を3次元座標に戻す
...
axis = (tip_center - root_center) / np.linalg.norm(tip_center - root_center)

断面上で求めた中心 (a, b) を、e2 / e3 を使って 3 次元に戻します。根元と先端で同じことをして、2 点を結んで正規化すれば、それが最終的な軸ベクトルです。

---

10. ステップ⑧：壁の法線を求める

角度の基準になる壁面の向きも、点群から求めます。

all_wall = np.vstack(wall_points)          # 全鉄筋まわりの壁点を積み重ねる
wall_center = all_wall.mean(axis=0)
u, sv, vt = np.linalg.svd(all_wall - wall_center, full_matrices=False)
wall_normal = vt[2]                        # 一番広がっていない方向 = 壁の法線
if wall_normal @ axes[0] < 0:              # 鉄筋と同じ向き (壁から外向き) にそろえる
    wall_normal = -wall_normal

壁の点は各鉄筋の処理中に集めておいたものです。

wall_ring = (r > WALL_RING_IN) & (r < WALL_RING_OUT) & (np.abs(s - face) < WALL_THICKNESS)
wall_points.append(near[wall_ring])

軸から 50〜200mm、かつ壁面から ±15mmの点、つまり各鉄筋の周囲のドーナツ状の壁面です。これを全鉄筋ぶん集めると、壁面全体にまばらに散った点群になります。

ここで再び SVD。平面上に分布する点群は 2 方向に大きく広がり、残り 1 方向にはほとんど広がらないので、vt[2]（最小主成分）が平面の法線になります。ステップ③と同じです。

1 本の鉄筋まわりだけでは狭すぎて法線が不安定になりますが、全鉄筋ぶんを合わせることで壁面全体に基線が伸び、精度が出ます。

---

11. ステップ⑨：角度を出す

法線と各軸を、球面座標的に分解して比較します。

wall_horizontal = np.degrees(np.arctan2(wall_normal[1], wall_normal[0]))  # 方位角
wall_slope      = np.degrees(np.arcsin(wall_normal[2]))                   # 仰角

horizontal_dev = np.degrees(np.arctan2(axis[1], axis[0])) - wall_horizontal
vertical_dev   = np.degrees(np.arcsin(axis[2])) - wall_slope
total_dev      = np.degrees(np.arccos(min(1, axis @ wall_normal)))

• 水平のずれ：XY 平面に投影した方位角の差（＋ = 壁側から先端を見て左）
• 鉛直のずれ：仰角の差（＋ = 上向き）
• なす角：内積から直接出す、軸と法線の 3 次元的な角度

12. 出力：CSV と DXF

CSV

with open(out_dir + r"\鉄筋軸_結果.csv", "w", newline="", encoding="utf-8-sig") as f:

BOM 付きで書くと Excel がダブルクリックで開いても文字化けしません。日本語見出しの CSV を吐くときは、ちゃんと指定しましょう。

出力する列はこちら。

header = ["鉄筋", "根元X", "根元Y", "根元Z", "先端X", "先端Y", "先端Z",
          "水平のずれ°", "鉛直のずれ°", "法線とのなす角°", "根元直径mm", "先端直径mm"]

DXF

CAD で図面に重ねて見たいので、軸線をそのまま DXF に書き出します。ライブラリは使わず、DXF R12（AC1009）の最小構成を直接書いています。

def write_dxf_lines(path, starts, ends, layer="REBAR_AXIS"):
    with open(path, "w", encoding="ascii") as f:
        f.write("0\nSECTION\n2\nHEADER\n9\n$ACADVER\n1\nAC1009\n0\nENDSEC\n")
        f.write("0\nSECTION\n2\nENTITIES\n")
        for k in range(len(starts)):
            p, q = starts[k], ends[k]
            f.write(f"0\nLINE\n8\n{layer}\n")
            f.write(f"10\n{p[0]:.6f}\n20\n{p[1]:.6f}\n30\n{p[2]:.6f}\n")
            f.write(f"11\n{q[0]:.6f}\n21\n{q[1]:.6f}\n31\n{q[2]:.6f}\n")
        f.write("0\nENDSEC\n0\nEOF\n")

DXF は「グループコード」と「値」が 1 行ずつ交互に並んだテキスト形式です。

コード	意味
0	エンティティの種類（LINE、SECTION など）
8	レイヤ名
10, 20, 30	始点の X, Y, Z
11, 21, 31	終点の X, Y, Z

線を引くだけなら これで十分です。ezdxf のような立派なライブラリを入れる必要はありませんでした。レイヤ名を付けておけば、CAD 側で色分けや表示制御もできます。
ここはAIに書いてもらいました。こういう知識ゲーはAIに頼るべきです。

---

13. パラメータ一覧（単位はすべて m）

定数	値	意味・決め方
GAP	0.10	これ以上すき間があいたら別の段・別の鉄筋。鉄筋間隔より小さく、鉄筋の太さより大きく
MIN_POINTS	500	これ未満のグループはノイズとして捨てる
BOX_HALF	0.7	全体点群を粗く絞る箱の半径（計算量削減用）
TUBE_ROUGH	0.08	最初に鉄筋を切り出す筒の半径（粗い軸でも確実に入る余裕を見る）
TUBE	0.045	鉄筋表面とみなす筒の半径（D51 の半径 約27mm ＋ 余裕）
WALL_RING_IN/OUT	0.05 / 0.20	壁を探すドーナツ領域（内側は鉄筋を避け、外側は隣の鉄筋を避ける）
WALL_BIN	0.005	壁検出ヒストグラムの刻み
WALL_THICKNESS	0.015	壁面とみなす厚み（±15mm）
FACE_MARGIN	0.01	壁から 10mm 以内は鉄筋に含めない
BODY_ROOT_CUT / BODY_TIP_CUT	0.05 / 0.02	軸を取り直すときに除く両端
SLICE_THICKNESS	0.040	輪切りの厚さ（節の間隔 35.6mm より厚く）
ROOT/TIP_SLICE_CENTER	0.030 / 0.035	輪切りの位置
TRIM	0.005	円周からこれ以上離れたら外れ点

別の径の鉄筋に適用するなら、TUBE と SLICE_THICKNESS を修正してください。
前者は公称半径＋数 mm、後者は節間隔より厚くしましょう。

---

14. 今後の課題

現状の実装で分かっている弱点を挙げておきます。

• ハードコードされたパス … 冒頭の tip_path / full_path / out_dir が Windows の絶対パス直書きなので、argparse でコマンドライン引数にするべき。
• LAS 1.4 非対応 … 前述のとおり、レガシー点数フィールドが 0 になるため落ちる。
• 壁面位置のビン端バイアス … find_face がビンの左端を返しているぶん、最大 5mm ずれる。
• 先端だけの点群を手で作る必要がある … ここが唯一の手作業。全体点群から「壁面から最も遠い領域」を自動で抜ければ、完全自動化できる。
• 関数と実行部が地続き … スクリプト全体がトップレベルに書かれているので、main() に切り出せば他から再利用できる。

---

おわりに

「点群処理」と聞くと専用ライブラリと機械学習の世界に思えますが、対象の形状が事前に分かっているなら、やることはお馴染み線形代数。NumPy だけで 279 行。むしろ中で何をやっているかが全部見えるぶん、パラメータを物理世界に合わせられますし、説明や検証も容易です。

点群測定をするために現場へ向かわれる方、どうかご安全に。