fullseye

数学(計測を支える数値計算) — 使い方ガイド

この族は何をする道具箱か

Fullseye の計測を 裏で支えている数学 を、第一級の op として表に出した族です。カメラ較正は最小二乗問題、ノイズ雲は共分散行列、主軸は固有ベクトル、歪みモデルは多項式、逆引きは補間 — measure / measure3d / camera の計測 op はすべてこの層の上で動いています。現在 4 分野 26 op(numpy + scipy のみ、台帳は opsmath.py):

データ種は画像ではなく数値配列: matrix(厳密に 2-D)、signal(厳密に 1-D)、標本集合は (N, D)(行 = 観測、列 = 変数)。table(dict)・pairs(counts, edges)・measurement(スカラ)・roots(complex 配列)を返す op もあります。FFT / 複素演算(complexops / volfreq / dsp)、1-D フィルタ(dsp / funct1d)、幾何フィット(measure / measure3d / pcseg)、codegen 参照実装の数値計算(algo)は既存モジュールの持ち場なので、ここには 重複させていません

ファミリ共通の入力契約(fail-closed)

全 op が入力を検証してから計算します(2026-08 の敵対監査で確定したバグ族を、黙って通さず明示拒否):

代表的なパイプライン(op の繋がり)

深度カメラで平面(定盤)を測る計測ワークフロー(検証済み examples/math_metrology.py そのもの)。フィット → 残差統計 → 主軸化 → 較正 → 逆引きと、データ種が matrix/signal → table → matrix → table → signal で繋がります。

flowchart LR
    A[点群 x,y,z] -->|設計行列 A| B[mat_lstsq 平面フィット]
    B -->|残差 z - Ax| C[stat_describe / stat_histogram / stat_zscore]
    A -->|2D ノイズ雲| D[stat_covariance]
    D -->|対称 PSD 行列| E[mat_eigh 主軸・楕円]
    B2[較正データ r, r_meas] --> F[poly_fit 3次]
    F -->|coeffs| G[poly_eval 往路 / poly_roots 逆算]
    B2 --> H[interp_linear / interp_cubic 逆引き表]

族の内部構造(裏の共有機構)。mat_cond は linalg 全体の健全性の見張り、SVD は mat_lstsq / mat_pinv / mat_cond / poly_fit(Vandermonde の条件数)の共通土台です。

flowchart TB
    SVD["SVD(gesdd)"] --> LSQ["mat_lstsq(gelsd)"]
    SVD --> PINV["mat_pinv(rcond 明示)"]
    SVD --> COND["mat_cond = smax/smin"]
    COND -.->|"log10(cond) 桁が消える"| SOLVE["mat_solve(gesv LU)"]
    COND -.->|"Vandermonde 条件数を記録"| PF["poly_fit"]
    COV["stat_covariance(対称 PSD)"] --> EIGH["mat_eigh(syevd・対称限定)"]
    COV --> CORR["stat_correlation(定数列は拒否)"]

使い方(op グループ別)

呼び出しは直接呼び: import mathops; mathops.mat_lstsq(A, b)(または opsmath.get("mat_lstsq"))。HALCON 対応は各 op ノート参照(linalg は Matrix 章、記述統計は Tuple 章が相当。共分散・相関・多項式フィット・求根は HALCON に公開 tuple op が無く、較正内部に隠れているものをここでは明示 op 化)。

linalg(密行列 — 条件数を見てから信じる)

stats(残差・ノイズの特徴づけ)

interp / poly(較正曲線と逆引き)

complex(複素解析 — 閉曲線は「点列」である)

輪郭は 複素点の 1-D 配列(cpoints)、閉じる辺は暗黙(先頭点を末尾で繰り返さない)。向きが答えの符号そのものなので orientation は明示引数です。

flowchart LR
    C[cplx_contour_circle 閉曲線] --> P[cplx_poly_eval 輪郭上で f をサンプル]
    P --> I[cplx_contour_integral ∮f dz]
    P --> AP[cplx_argument_principle 零点-極の数]
    P --> CV[cplx_cauchy_value 内部の値 f w]
    P --> L[cplx_laurent_coeffs 係数・留数]
    C --> W[cplx_winding_number 巻き数]
    C --> J[cplx_joukowski / cplx_mobius 等角写像]
    IMG[cx_fft などの複素場] --> CR[cplx_cr_residual 正則性の残差]

動く最小例(検証済み)

repo 直下で py -3.11 の対話環境か、PYTHONPATH に repo を通して実行。フィット厳密復元・PSD/直交性・SVD⇔固有値の交差検証・fail-closed(範囲外拒否)を数値で確認して PASS を出します(本ガイド作成時に実行し PASS を確認済み。16 op 全てを通すフル版は py -3.11 examples/math_metrology.py)。

import numpy as np
import mathops as M

# GT: 直線 y = 2x + 1(ノイズ無し・決定的)を最小二乗で厳密復元
x = np.linspace(0.0, 1.0, 21)
y = 2.0 * x + 1.0
A = np.column_stack([x, np.ones_like(x)])     # (21, 2) 設計行列

fit = M.mat_lstsq(A, y)                       # SVD 最小二乗(健全性テレメトリ付き)
assert np.allclose(fit["x"], [2.0, 1.0], atol=1e-12) and fit["rank"] == 2
assert M.mat_cond(A) < 1e2                    # 条件数を必ず見る(健全)
assert fit["residual_ss"] < 1e-24             # ノイズ無しデータ → 残差ゼロ

d = M.stat_describe(y - A @ fit["x"])         # 残差統計: 実質ゼロ
assert d["n"] == 21 and abs(d["mean"]) < 1e-12

# 共分散 → 固有分解(主軸): 相関のある 2 変数の雲(seed 固定・決定的)
rng = np.random.default_rng(0)
u = rng.standard_normal(500)
cloud = np.column_stack([u, 0.5 * u + 0.1 * rng.standard_normal(500)])
C = M.stat_covariance(cloud)                  # (2,2) 対称・半正定値
w, V = M.mat_eigh(C)                          # 対称行列専用(対称性を検証してから解く)
assert w[0] >= -1e-12                         # PSD → 固有値は非負
assert np.allclose(V.T @ V, np.eye(2), atol=1e-12)   # 固有ベクトルは正規直交
assert M.stat_correlation(cloud)[0, 1] > 0.9  # 強い正相関(構成通り)
assert np.abs(M.stat_zscore(cloud[:, 0])).max() < 5.0  # 正規サンプルの z-score

# SVD 交差検証: centered データの s^2/(N-1) は covariance の固有値と一致
centered = cloud - cloud.mean(axis=0)
_, s, _ = M.mat_svd(centered)
assert np.allclose(sorted((s ** 2) / (cloud.shape[0] - 1)), w, rtol=1e-10)

# 較正曲線: y = 0.15 x^3 + x を 3 次 poly_fit で厳密復元(条件数もその場で確認)
pf = M.poly_fit(x, 0.15 * x ** 3 + x, 3)
assert np.allclose(pf["coeffs"], [0.15, 0.0, 1.0, 0.0], atol=1e-9)
assert pf["cond"] < M.POLY_COND_WARN          # 1e10 を超えると RuntimeWarning
assert abs(M.poly_eval(pf["coeffs"], 0.5) - (0.15 * 0.125 + 0.5)) < 1e-12
r = M.poly_roots([1.0, 0.0, 1.0])             # x^2 + 1 = 0 → ±i(複素も正直に返す)
assert np.allclose(sorted(r.imag), [-1.0, 1.0], atol=1e-12)
assert M.poly_roots([1.0, 0.0, 1.0], real_only=True).size == 0  # 実根は無し(空が正解)

# 補間: 範囲外クエリは既定で ValueError(clamp は明示指定)
tbl_x = np.array([0.0, 1.0, 2.0, 3.0])
tbl_y = tbl_x ** 2
assert M.interp_linear(tbl_x, tbl_y, 1.5) == 2.5   # 区分線形: (1+4)/2
assert abs(M.interp_cubic(tbl_x, tbl_y, 1.5) - 2.25) < 1e-12  # 3 次多項式は厳密再現
try:
    M.interp_linear(tbl_x, tbl_y, 9.0)             # fail-closed(外挿は既定拒否)
    raise AssertionError("out-of-range must raise")
except ValueError:
    pass
assert M.interp_cubic(tbl_x, tbl_y, 9.0, out_of_range="clamp") == 9.0  # 端値を保持

# ヒストグラム: 度数和 = サンプル数(明示ビニング)
counts, edges = M.stat_histogram(cloud[:, 0], bins=12)
assert int(counts.sum()) == 500 and edges.size == 13

print("PASS")

数式(必要な op のみ)

最小二乗(mat_lstsq)と擬似逆行列(mat_pinv)は同じ問題の二つの顔:

\[\hat{x} = \arg\min_x \lVert A x - b \rVert_2^2, \qquad A^{+} = V \, \Sigma^{+} U^{\top} \ \ (\sigma_i < \mathrm{rcond}\cdot\sigma_{\max} \text{ は } 0 \text{ 扱い})\]

条件数(mat_cond)と「消える桁数」(Golub & Van Loan §2.6 — mat_solve を信じてよいかの判定):

\[\kappa_2(A) = \frac{\sigma_{\max}(A)}{\sigma_{\min}(A)}, \qquad \text{失う有効桁数} \approx \log_{10} \kappa_2(A)\]

標本共分散(stat_covarianceddof=1)と Pearson 相関(stat_correlation)、z-score(stat_zscoreddof=0):

\[C = \frac{1}{N-1} \sum_{i=1}^{N} (x_i - \bar{x})(x_i - \bar{x})^{\top}, \qquad R_{jk} = \frac{C_{jk}}{\sigma_j \sigma_k}, \qquad z_i = \frac{x_i - \bar{x}}{\sigma}\]

多項式(poly_fit / poly_eval、係数は最高次から)と、そのフィットの条件数(Vandermonde 行列 $V_{ij} = x_i^{\,d-j}$ の $\kappa_2$ — POLY_COND_WARN = 1e10 超で警告):

\[p(x) = \sum_{k=0}^{d} c_k \, x^{\,d-k}, \qquad \hat{c} = \arg\min_c \lVert V c - y \rVert_2^2\]

サンプルデータ

この族の入力は画像ではなく数値配列なので、seed 固定の合成データ + 解析的グラウンドトゥルース(真の係数・真の固有値・厳密根)がそのまま最良のテストデータです — 上の最小例と examples/math_metrology.py がその作り方の見本(平面 + 既知ノイズ、回転楕円雲、樽型歪み風の較正曲線)。実測データに繋ぐなら、measure / measure3d の出力(残差列・点群座標)をそのまま signal / (N, D) として渡せます。画像系サンプルの台帳は ../../SAMPLES.md

参考文献(正典)

台帳は ../../../REFERENCES.md。この族の数値解析の古典(各 op の docstring が引用しているもの):


© 2026 Kazufumi Furuse — Fullseye operator documentation. Licensed under Apache-2.0.