ふたつの usb カメラの映像を使って，3次元空間での放物運動などの座標を測定しようと思います（今のところは，そこまでは至ってないですが）。学生実験のネタにならないかと思っていろいろと調査してみました。
参考サイト
- カメラ キャリブレーションとは
- コンピュータビジョンと3次元復元
- ピンホールカメラモデルとカメラキャリブレーション
- 第 9 回:三次元映像解析
- Wolfram言語 ＆ システムドキュメントセンター：LeastSquares
DeepLabCut というソフトを利用して学生実験ネタをひとつ作ったのですが，それは水面に浮かべた二つの磁石の引き合う様子から作用反作用を観測するというものでした。
カメラのピクセル座標は実際の空間座標からみれば歪んでいます。加速度を求めるなら，実座標に変換するとともに歪の補正をする必要があります。この時には，1 cm の格子が描かれたカッティングマットを撮影することで，ピクセル座標を実座標に変換しました。
カメラを動かさずに，同じ位置から撮影すれば実座標に変換する式が得られると考えたのです。歪の補正も同時になされます。ピクセル座標の２次まで考慮して換算式を作成しましたが，1次でもそれほど変わりませんでした。
上記は水面上に運動が限られています。これを放物運動などに広げようと考えたら 3 次元の座標を観測する必要があります。3次元だと二つのカメラが必要になると思うのですが，二つのカメラとなると運動の観測の場合は時刻を合わせることが必要です。お金のかからない方法で何かないかと探していたら，amazon で双眼のカメラを見つけました。
これなら時刻が同期された映像が得られそうです。撮影してみると下のような予想通りの映像が得られました。ふたつ分の映像が横に並んで出力されてきます。
歪が大きかったので試しに紙に描いた円を撮影してみました。円ならば歪まないだろうと考えて撮影したのですが，円を画面の中心付近に置いた場合にも歪が生じていて，中心から外れた位置に歪の小さなところがありました。どうも軸がでていないようです。これでは取り扱いが難しいだろうと考えて，このステレオカメラを利用するのは諦めました。
それで他の方法を調べていたら，OBS Studio というソフトが usb カメラを2台（以上）同時に扱えることが分かりました。複数台のカメラを合成してひとつの仮想のカメラに見せることも可能なようです。カメラとしては普通の WEB カメラ（ロジクール c922）を使います。1280 x 720 だと 60 fps です。1920 x 1080 なら 30 fps です。もう少し fps が大きいと良いのですが，それはまた後で考えます。このカメラの歪は意外に小さいです。下の写真は 1 cm 格子のカッティングマットを上空から撮影したものです。
中心から遠いところは歪が大きくなりますが，両端で 10 px 程度のずれでした。真ん中あたりを使うなら，ほとんど補正が必要ないかもしれません。補正をしたあとであれば（あるいは補正が必要ないデータならそのままで），測定結果をピンホールカメラとして扱えます。ピンホールカメラに関しては，先にも上げましたがこちらを参考にしました。補正の歪に関しては，これも前にあげましたがこちらを参考にしました。
動いている物体の撮影をしたいので，シャッター速度も問題になるのですが，このカメラは Logi Options+ というソフトで設定することができます。ただし，windows のカメラアプリを起動すると，これらの設定がクリアされてしまいます。起動しなければ設定は維持されます。
下の写真は OBS Studio のキャプチャーです。
静止画像を重ねることもできるので，上の写真では，狙いを定めるための黄色い線を重ねています。2台のカメラで録画した映像が下記の動画です。2台のカメラの映像を横に並べて配置しました。
60 fps の動画です。コマ送りをすると分かるのですが，時々コマ落ちしているようです。これはまだ今のところ防げていないので，何か時間を表すものを同時に撮影して時刻を補正しようと考えています。でもそれならば，パソコン2台で撮影してもおんなじなのですが。パソコンの処理能力というのも関係してきそうなので，最終的な撮影方法はまだ決まっていません。とにかく，以降の話は 2 台のカメラで3次元の座標を得ることだけに限定した内容です。
実空間とカメラのピクセル座標の関係は，こちらの資料を参考にしました。それによると（PDFの 7 ページあたり），ピクセル座標を\(u, v\), 実空間の座標を \(X, Y, Z\) とすると，その間には下記の式が成立します。
- \(\displaystyle u = \frac{L_1 X + L_2 Y + L_3 Z + L_4}{L_9 X + L_{10} Y + L_{11} Z + 1} \)
- \(\displaystyle v = \frac{L_5 X + L_6 Y + L_7 Z + L_8}{L_9 X + L_{10} Y + L_{11} Z + 1} \)
ここで，\(L\) はカメラパラメーターと言われる値で，実空間におけるカメラの配置の条件のすべてを含む定数です。その為，カメラを移動させると計算しなおす必要があります。両式の右辺の分母は共通ですが，両式ともに，この分母を両辺にかけて線形化します。
- \(\displaystyle u \, \left( L_9 X + L_{10} Y + L_{11} Z + 1 \right)= L_1 X + L_2 Y + L_3 Z + L_4\)
- \(\displaystyle v \, \left( L_9 X + L_{10} Y + L_{11} Z + 1 \right) = L_5 X + L_6 Y + L_7 Z + L_8 \)
上記を移項して下記のようにします。
- \(\displaystyle u = L_1 X + L_2 Y + L_3 Z + L_4 – L_9 \, u X – L_{10} \, u Y – L_{11} \, u Z\)
- \(\displaystyle v = L_5 X + L_6 Y + L_7 Z + L_8 – L_9 \, v X – L_{10} \, v Y – L_{11} \, v Z\)
上式に，観測データをあてはめて最小二乗法で \(L\) を決定します。それには，実空間の座標 \(X, Y, Z\) とそれと対応するピクセル座標 \(u, v\) の組がデータとして必要です。
下の写真のような発砲スチロールにグラフ用紙を貼り付けたものを用意しました。グラフ用紙には，6 cm 正方の格子を書きました。
これをカッティングマットの上に載せて移動させます。5 cm ずつ前進させました（前進方向の軸が X 軸）。距離感としては，下のような配置です。
OBS Studio では下のような感じになります（発砲スチロールの位置が上の写真と比べると多少前進していますが）。
一枚だけ，座標を計測した写真をあげます。\(X = 0\) の場合の写真です。全部で 4 枚の写真（\(X = 0, 5, 10, 15\)）からピクセル座標を読み取りました。
グラフ用紙上の格子点の座標を記録したものが下のエクセルシートです。X, Y, Z が実空間での座標で，x1, y1 がカメラ１， x2, y2 がカメラ２のピクセル座標です。ピクセル座標はそれぞれの画面の中央を原点にして，y方向を反転させたものも計算しています。
これがもとのエクセルファイルです。この測定値から，下図のような回帰分析用の列データを作ります（参考資料のPDFは一部マイナス符号が抜けているように思います）。\(u\) の式と \(v\) の式の二つがあって，係数 \(L\) は共通しています。それでシートでは，列の上の方は \(u\) の式から由来する条件が並び，その下に \(v\) の条件のものが並んでいます（下図では \(u\) の条件のところだけが見えていますが，この下に \(v\) の条件が続いています）。列に2種のデータが混在しているのです。それを一緒に合わせて回帰分析をやります。
上の図のデータから回帰分析で，カメラ１とカメラ２用の二つの \(L\) を下図のように定めました（先のエクセルファイルを参照）。
これで最初に示した式を使って，実空間の座標 \(X, Y, Z\) から ピクセル座標 \(u, v\) を計算することができますが，欲しいのは逆で，ピクセル座標から実空間の座標を得る方です。先の \(u, v\) の式を変形します。
- \(\displaystyle u – L_4 = \left(L_1 – L_9 \, u \right) \, X + \left(L_2 – L_{10} \, u \right) Y + \left( L_3 – L_{11} \, u \right) Z\)
- \(\displaystyle v – L_8 = \left(L_5 – L_9 \, v \right) \, X + \left(L_6 – L_{10} \, v \right) Y + \left( L_7 – L_{11} \, v \right) Z\)
上記の式がそれぞれの \(L\) に関してあるので計4式です。 \(u, v\) を測定するとして，その4つの式から \(X, Y, Z\) を決めるという話です。ここは Mathematica の LeastSquares 関数を使用しました（Wolfram Cloud でもなんとかなりました。あとで相違点を書きます）。まずカメラパラメーターと測定データの読み込みです。ひとつの CSV ファイル（20260223a 右クリックで保存）がカメラパラメーターで，もうひとつ（20260223b）が測定値です。
In[224]:= Remove["Global`@*"];
tmpList = Import["20260223a.csv","Data"];
cameraList = Drop[tmpList,1]
Length[cameraList]
tmpList = Import["20260223b.csv","Data"];
measureList = Drop[tmpList,1]
Length[measureList]
Out[226]= {{13.309,-17.432},
{27.4661,20.7241},
省略
{-0.00807954,-0.00612487}
}
Out[227]= 11
Out[229]= {{0,0,18,699,145,2790,113,-261,395,-90,427},
{0,6,18,893,120,2935,136,-67,420,55,404},
省略
{15,18,0,1482,781,2926,1023,522,-241,46,-483}
}
Out[230]= 64続いて，最小２乗法です。元の座標を復元できるかどうか試してみます（入力も出力も一緒になってしまってわかりずらい）。
La =cameraList[[All,1]];
Lb =cameraList[[All,2]];
X =measureList[[All,1]];
Y =measureList[[All,2]];
Z =measureList[[All,3]];
ua =measureList[[All,8]];
va =measureList[[All,9]];
ub =measureList[[All,10]];
vb =measureList[[All,11]];
Do[Print[i];
yy = {ua[[i]] - La[[4]], va[[i]]- La[[8]], ub[[i]] - Lb[[4]], vb[[i]] - Lb[[8]]};
WW = {{La[[1]] - La[[9]] * ua[[i]], La[[2]] - La[[10]] * ua[[i]], La[[3]] - La[[11]] * ua[[i]]},
{La[[5]] - La[[9]] * va[[i]], La[[6]] - La[[10]] * va[[i]], La[[7]] - La[[11]] * va[[i]]},
{Lb[[1]] - Lb[[9]] * ub[[i]], Lb[[2]] - Lb[[10]] * ub[[i]], Lb[[3]] - Lb[[11]] * ub[[i]]},
{Lb[[5]] - Lb[[9]] * vb[[i]], Lb[[6]] - Lb[[10]] * vb[[i]], Lb[[7]] - Lb[[11]] * vb[[i]]}};
xx = {X[[i]], Y[[i]], Z[[i]]};
Print[xx];
realZahyo = LeastSquares[WW,yy];
Print[realZahyo]
,{i,Length[measureList]}]
1
{0,0,18}
{0.0724112,0.0404099,18.0192}
2
{0,6,18}
{-0.0406213,6.02092,17.9919}
3
{0,12,18}
{0.0371127,12.0028,18.0047}
4
{0,18,18}
{-0.0503242,17.9829,18.0233}
省略
64
{15,18,0}
{14.8556,18.002,-0.0205645}
座標の単位は cm です。一番大きくずれるところで， 2 mm ぐらいずれていましたが，まあ良い結果ではないかと考えています。
Wolfram Cloud の場合には，（無料の契約だと）データファイルのアップロードができません。何とかコピペでデータを移す必要があります。それで，こんな感じの csv ファイルを作成しました。
このファイルを .prn 拡張子のファイルとして書き出します。それをメモ帳で開けば下図のような感じになります。
広く空いた隙間はスペースで埋められています。この形式ならコピペで Wolfram Cloud に移すことができます。あとはローカルの場合と同じです。