作用反作用をテーマとした、学生実験のテーマをつくりたいと考えて色々とやってみた記録です。プリンキピア（中央公論社『世界の名著 26 ニュートン』昭和46年、P84）には、磁石と鉄を水に浮かせて実験した例があげられています。それは多少違った文脈ですが、その話から磁石を水面に浮かべて作用反作用を確認するようなものを学生実験のテーマとしてやってみようと考えました。以下、その準備を記録します。
まず磁石とそれを載せる船ですが、下図のようなものを用意しました。
船にはシャーレを使います。シャーレを使う理由は特にありません。シャーレの中央の黒い球体が磁石です。シャーレには、位置を検出するために丸い色つきのシールを貼っています。このシールは正三角形の頂点に位置していて、この3点の平均として磁石の位置を定めます。最初は磁石の位置を直接検出していたのですが精度が悪かったので、３点の平均を使用するようになりました。
上図のように、シャーレをトロ船（モノタロウで購入）に貯めた水に浮かべます。水は20リットル程度いれます。シャーレの動きには、周りの鉄が影響します。机や床には鉄が入っているので、ある程度水深がある方が鉄との距離が取れると思って多めに水を入れています。鉄の影響以外にも、水深とシャーレの動きには何か関連があるかもとは思いますが未検討です。
水道の設備が近くにないので、ポリタンクと水を回収するためのポンプを準備しました。
このポンプは強力過ぎました。あっという間にポリタンクがいっぱいになります。
続いて、撮影機材です。シャーレの動きを上空から撮影します。
カメラの三脚で、エレベーター（と呼ぶらしい）を逆向きに取り付け可能なタイプを用意して USB カメラで上空から撮影します。カメラから水面までは 70 cm 程度離しました。トロ船も大きいので、三脚も大きめのサイズが必要です。人の目線の高さぐらいが確保できるぐらいの三脚が必要だと思います。
下から、カメラを写してみました。
カメラは logicool C922n です。三脚固定用のねじ穴が付いています。
カメラの画像から得られるピクセルの値を実際の座標に換算しなければならないのですが、下記に記すような感じでやってみました。まず、方眼が印刷されているカッティングボードを用意します。
これを適当なものの上にのせて、水面に浮かべます。
この状態を USB カメラで撮影します。撮影には windows に付いているカメラアプリを使いました。動画として撮影して一枚静止画を切り出しました。その瞬間で座標軸を定めます。下図のように、10 cm 格子のピクセル座標値を読み取ります。読み取りには windows 付属の ペイントを使いました。
原点が中途半端な位置にはなりますが、上図のピクセル座標を、縦が 15 cm ～-15 cm まで、横が -20 cm ～20 cm までに割り当てます。
Excel で重回帰分析をします。下図は、x 座標をピクセル座標値から決める際のキャプチャーです。２次まで考慮したかったのでピクセル座標値の２乗の列を追加しています。
これらの値から、下記の式の係数 A, B, C, D, E を決定します。
- x = A + B*Px + C*Py + D*Px^2 + E*Py^2
回帰分析は Excel の「データ」「分析」「データ分析」にあります。
下図が分析結果です。
得られた式から、再現した x 座標を計算して比較しました。右端の列が x 座標の再現値です。だいたい合っています。
２次まで考慮したのですが、１次とさほど変わりませんでした。とにかく、この式を用いて、ピクセル座標値を実際の座標に変換します。これ以降実験終了まで、三脚やカメラの位置を動かさないよう注意しなければなりません。
さて、撮影です。磁石の進行方向は地磁気の方向と揃えた方が余計な回転が生じませんが、回転させないというのは結構難しいです。撮影中にフォーカスが変化してしまうことが何度かありました。カメラはオートフォーカスではなく固定のフォーカスにした方が良いです。撮影した動画を、下記にひとつあげます。手で磁石をスタートさせるのは静かに離すのが難しいので、紙で適当なものを作りました。
上の動画は少し前に他の場所で撮影したものなのですが、今回はこの動画を用いて、機械学習のトレーニングを実施します。その学習結果を使って今回の撮影動画の磁石の座標を読み取ります。
話が少し前後しましたが、座標の読み取りには deeplabcut を利用します。動物の動きなどを追いかけるためのソフトです。インストール方法は別の記事を書きました。windows の start から Anaconda Prompt を起動します。
下図のような端末の画面が開きます。環境を DEEPLABCUT に切り替えます。
- conda activate DEEPLABCUT
deeplabcut を起動します。
- python -m deeplabcut
起動しました。
プロジェクトを新しく始める場合、create new project を選びます。
プロジェクトに名前を付け、実験する人の名前を記入します。これらの情報からファイルを収めるフォルダー名が決まります。デフォルトではデスクトップにフォルダーが作られます。下方ではトレーニングに利用する動画を選択します。Copy videos to project folder にチェックをいれておきます。Browse folderrs をクリックして、動画のあるフォルダーを指定すると、そのフォルダーにある動画のリストが表示されます。必要な動画以外のチェックを外します。
create ボタンを押して閉じます。設定をデフォルトから変更したいので、Edit config.yaml をクリックします。下図のような画面が開きますが、設定のうち body parts の数が足りません。下図のように bodypart を 6 に増やします。ダブルクリックで編集できます。枝を右クリックすると、新たに項目を追加することができます。
save ボタンを押して窓を閉じます。デスクトップに下図のような中身のフォルダーが作成されています。
config.yaml には先ほどの設定が書き込まれています。videos フォルダーには、選択した動画がコピーされています。
先に進みます。Extract frames タブを表示します。select video ボタンでトレーニング用の動画を指定します（先ほどコピーした動画）。1 video selected と表示されます。
Extract frames ボタンをクリックします。
訓練用の静止画の抽出が完了しました。下図のように、labeled-data フォルダー内に 20 枚の画像が切り出されています。
次に Label frames タブに移ります。ここではトレーニング用の正解を用意します。
Label Frames ボタンをクリックします。フォルダーの選択が促されるのですが、先ほどの 20 枚の静止画があるフォルダーを指定します。その後、下記のような窓が開きます。最初に表示されるチュートリアルはとばします。窓を大きくしてパソコンのモニター画面最大サイズにした方が作業がしやすいと思います。
窓の下方に動画の再生制御機能が並んでいます。動画を再生しながら、追跡したい物体の移動範囲を確認します。これから物体の位置にしるしをつけていくのですが、作業しやすいように動画を拡大します（マウスホイール）。確認したら動画を静止させ再生位置を1枚目の画像に戻しておきます。
左上方にある丸にプラスマークが入ったアイコンをクリックして座標位置の指定を始めます（クリックするとアイコンが青くなります）。マウスクリックで追跡する位置を指定します。下図は1枚目の画像の作業中で、 3 ポイントしるしを付けたところです。左側の磁石を囲む 3 個のシールに印がついています。マウスをクリックすると自動的に次のポイントの指定に移ります。6 ポイント印をつけたら、次の画像に切り替えます。
20枚全部に印を付けたら終了です。窓を閉じる前に作業結果を記録します。下図のようにメニューバーから File を選び、Save Selected Layer(s) を実行します(重要です)。保存したら窓を閉じます。
Create training dataset タブに切り替えます。ここでは何も変更せず、Create training dataset ボタンをクリックします。
データセットが作成されました。
続いて、Train network タブに移ります。ここでは Display iterations を 100 に変更します。Train Network ボタンをクリックして訓練を開始します。
タスクマネージャーを開くと、gpu のメモリ使用量が急上昇しているのが分かります。
途中経過は、 deeplabcut を起動した端末に表示されます。端末には、この他の作業中にも、色々とメッセージが表示されます。
終了しました。25 分程度かかりました。
端末の表示を見ると、後半では、もはや train loss は減少していません。
Evalute network タブに移ります。Evalute network ボタンをクリックします。
窓の下の辺にメッセージが表示されます。あまり意味は分かりません。
Analyze videos タブに移ります。すでに 解析するビデオが選ばれていると思いますが、他のビデオを解析する時には、ビデオの選択をやり直します。Save result(s) と Plot trajectories にチェックを入れて（すべての bodypart を選びました）、Analyzed videos ボタンをクリックします。
この時も、端末に途中経過等が表示されます。2 分 50 秒ほどかかりました。
videos フォルダーには読み取った座標値を記した csv ファイルが現れます。
csv ファイルの内容は下記のような感じです。
videos/plotposes フォルダーの中にその動画専用のフォルダーが作成されて、その中に位置座標の時間変化のグラフや軌跡の図が作成されます。下図は位置座標の時間変化のグラフです。
下図は、座標の読み取りにどれぐらい自信があるかを示す plot-likelihood のグラフです。かなり信用してよい様です。
Create videos タブに移ります。ここもすでに video が選択されていると思います。他の動画を処理したいときには選択をやり直します。ここでは、Overwrite videos と Plot trajectories にチェックを入れて、Create videos ボタンをクリックします。
読み取った位置を示すマークが上書きされた動画が作成されて、videos フォルダーに動画がひとつ増えていると思います。なぜか上書きが上手くいかないことがあって元の動画のままの場合がありますが、その時は作成された動画をいったん削除してもう一度 Create videos ボタンをクリックします。
読み取った位置が上書きされた動画は下記のような感じになります。
色付きシールのあたりの色が変わって少し大きめの丸になっています。丸を丸で上書きしているので分かりずらい。
学習が終了しました。このトレーニング結果は、他の動画ファイルの解析にも利用できます。以降、お見せするデータは、下の動画を解析したものです。
ピクセルと座標の換算式は下記
- x =-37.5196930450991+0.0389152535950154*Px-0.00104185457456728*Py-4.26190682209606*10^(-8)*Px^2-9.68517942973649*10^(-8)*Py^2
- y =21.868405629349-0.00117352701702892*Px-0.038150408407546*Py+4.15852434349064*10^(-8)*Px^2-6.97709029643541*10^(-7)*Py^2
質量は、磁石とシャーレを合わて、左側の小さい磁石のセットが 60.3 g、右の大きい磁石の方が 77.5 g です。
解析すると、6 セットのピクセル座標の組（Px, Py という感じの組）が得られます。それらから、二つの磁石に対応する二組の(Px, Py)を求めます。下図では、小さいほうの磁石のデータのみしか表示されていません。推論結果の先頭の部分です。3 ポイントのピクセル座標を平均して、ひとつのピクセル座標（見出しが Ax, Ay となっています）を計算し、そののち実座標に換算しています。time が 8 秒あたりから始まっていますが、それまでの時間のデータはカットしました。
磁石の座標が得られたので、まず重心の x 座標と y 座標の時間依存性を見てみます。
上図のグラフでは、重心が移動しています。移動することはあり得ることなのですが、その速度が変化していて何かの外力が働いていることを示しています。
重心の速度をグラフにしてみます。計算には slope 関数を利用しました。下図のように領域を定めて平均的な傾きを求める関数です。こんな感じに記述します。
- =SLOPE(F2:F6,$A2:$A6)
下図が重心速度のグラフです。5 個の測定値から、その傾きを求めています。それでもノイズが多いのですが。
重心の x 方向の速度はほぼ一定ですが、y 方向の速度は変化しています。それぞれの磁石の速度の大きさと比較してみます。下に x 方向のみですが、両磁石の速度をあげます。重心のグラフとは横軸が異なっていて、40 秒ぐらいまでを表示しています。
無視できるほどではありませんが、重心の速度は比較的小さいことが分かります。
続いて運動量を見てみます。最初に x 方向です。緑の点は両磁石の運動量の和です。和は一定の値となっています。
次に、y 方向の運動量です。
衝突の時刻において変化が見られます。回転の影響ではないかと考えています。両磁石の和はだんだんマイナス方向に大きくなっています。
最後に、それぞれの磁石の加速度から力を求めてみました。加速度の計算には LINEST 関数を使用しました。時間の2次式でフィットするために、時間を2乗した列も用意します。下記の式でフィットするのですが、
- x = a*time^2 + b*time + c
その２次の係数 a の２倍が加速度となります。 LINEST 関数は係数を配列として返します。その配列の インデックス 1 が2次の係数です。計算には、前後１秒間のデータを利用します。計算例を一つ上げます。セル E31 の計算です。
- =INDEX(LINEST(C2:C60,$A2:$B60),1)*2
下図はシートのキャプチャーですが、かなり広い範囲を計算に利用しています。
これによって求めた加速度に質量をかけて力とします。下図がそのグラフですが、大きい磁石の方の符号を反転させています。
波打っているところは、水面の波が影響しているのかもと思います。力が大きいところでは良い一致といえるのではないでしょうか。
(20241009)
学生向けのテキストを書いてみました。