前回7年ぶりに再トライしているPICですが、なかなか進みません。難しい。全然したいことの情報に行き着かない。。。
って感じで案の定、頓挫中なのですが、高機能プランターを作りたいっと言った友人はビールを作りたいと言い始めました。→ブログ http://brewojisan.com/
その際、発酵させる過程で炭酸ガスが出るのですが、発酵状態の把握の目安として、ガス発生の推移を計測したいと言いました。
具体的には、下の動画のように、水封管がガスでコポコポする頻度を測りたいとのこと。
久々のOpenCVで挑戦だ！
さて、前回こしらえた点と点をベジェ曲線で結ぶ自前のパスをブーリアン演算で差分を取りたいと言っていたのはどうなったのか。できました、一応。
ソースは一番最後
まぁ、最後の合成は余分なパスを生成してるけど、ご愛嬌ってことで。
以下、よもやまプログラミング話。
前回も理想のドローソフトについての記事でしたが、今だその熱を持っていて、今回はパスとかハッチングとかやってみました。
EMFやSVGなどに出力することを考えると、ライブラリが持つ機能に依存するよりは、なるだけ自前で実装したほうが良いはず。どちらのファイル形式もベジェ曲線はサポートしているので、３次ベジェ曲線によるパスのプログラミングを試みました。
ドローソフトと言えばこれだよね。
［ソース１］（ソースは最後にあります）
制御点を無視して、直線とすることもできるようにした。
機能として必須なのが、パターンハッチング。
ベジェ曲線と直線の交点の計算なんてどうやるんだーっと思ったけど、情報を載せてくれる人はいるもんですなぁ。ありがたい。
斜線パターン
［ソース２］
5(その2):直線と曲線の交点(解の公式)
パターンを繰り返すハッチング
［ソース３］
これで直線ならどんな模様でもできる。いずれ楕円も加えたい。
虚数とか出てきてなかなか大変でしたが、出張により久々に一人の夜を過ごし、たっぷり作業できたので、なんとか形にはなりました（後々不具合が出てくるのは必至だと思いますが）。久しぶりにがっつりプログラミング、楽しい。
次にしたいと思っているのが、ブーリアン演算による差分とか合成とか。果たして僕の能力が及ぶのか。。。
いつか理想のドローソフトを作りたいっと思いつつ、意思の弱さゆえに開発を進められていないのですが、仕事の書類などでMS Officeで図形を駆使して図を描くたびに、作りたいと言う思いが再燃します。
そのたびに関連することを調べ始めまして、今回もその一環です。
理想のドローソフトに無くてはならない機能が、ベクタ形式のまま簡単にMS Office Word等に図を貼り付ける機能。すなわち、EMF（拡張メタファイル）形式でクリップボードにコピーする機能です。
追記：2006年以降のMS OfficeもOpen Officeもsvgのコピーに対応してるんだって。EMFはオワコンって事か。。。
前回は下記の論文に参考に転がり抵抗を導入してみました。今回は、要素を土にするためにもう一つ必要な、粒子間の吸着力（土の世界でいう粘着力）を考慮したモールクーロンの破壊基準の適用にチャレンジします。
本来は、転がり抵抗モデルが意図した通りに組めているかの検証をしなくてはいけないところですが、やりたい・やれそうなことは、順序は関係なくやっていくスタンスです＾＾
地盤材料の破壊基準を表現するためのシンプルな個別要素モデル
何をするかと言うと、せん断に対する抵抗力を一定値ではなく、拘束力（直方向の力）に係数を掛けた値を上限とするモールクーロンの破壊基準を組み込みます。加えて、通常のDEMでは扱わない引張抵抗力（粘着力）を考慮するようなモデルにします。この引張力が転がり抵抗にも関与するので、それに応じて、転がり抵抗部分も変更します。
とにかくやってみよー。
ボンドモデル部分ソース（全ソースは最後に）
# 粒子間ボンドモデル
cdef double bs,tu=0
if i < peCount and j < peCount: #粒子同士
tu = pe[i].et[j]
if hn < -tu and tu > 0: #ボンド破壊条件
pe[i].et[j] = 0
tu = 0
print('bond break',i,j,hn,tu)
#モールクーロン破壊基準の適用
bs = (hn + tu) * tan(pe[i].ph)
if hn <= -tu:
hs = 0.0
elif fabs(hs) > fabs(bs):
hs = bs * hs/fabs(hs)
#hn = hn - tu
else: #粒子同士以外
if hn <= 0.0:
#法線力がなければ、せん断力は０
hs = 0.0
elif fabs(hs) >= inf.frc*hn:
#摩擦力以上のせん断力は働かない
hs = inf.frc*fabs(hn) *hs/fabs(hs)
上記論文ではDF条件と称して、接触状態にある粒子間では粘着力が発揮されていて、一定の引張力が加わったときに粘着がなくなり（ボンドが破壊）、その後の接触は通常の摩擦として扱うモデルになっています。
普通は離れている要素の粘着力は既にボンド破壊状態であり、粘着力は働かない（通常の摩擦であらわす）のですが、確かめるための便宜上、初期値では粘着力を有効にしておいて、破壊されるまでの間は、接触する要素間に粘着力が働くようにしておきます。
さて、どうなるか？
（左：粘着なし 右：粘着有り）
お、粘着力は働いてるっぽい。しかし、2つ問題が発生。動きが止まってきたあたりにバネのようにスライドする謎の動きと、後半、なぜか要素が離れていく・・・。
前者は、今回のボンドモデルに関係なく、線と円要素の接触処理によるところの部分の不具合ですかね。そういえば、前の転がり抵抗のときもそうでした。摩擦の上限を超えることなく、せん断バネでスライドしているせいかな？検証が必要ですなー。パラメータが非適正な値になっているだけだといいんだけど。
後者について、実はイメージで言うと、今回の動きは、僕の想像と違って、もっと粒子同士が重なって、重なった状態で安定するイメージでした。上のソースだと、粒子が反発して離れる過程で、引張力（粘着力）が働くけども、発揮される直方向の力自体の値には影響を与えません。収束状態は力の釣り合ったところなので、ボンドが破壊するまで反発が働き続けると、結果、離れてしまうのかなと。
ともすれば、最終的な力は、弾性力（反発力）に引張力を差し引いた値にするのでは？と思い、上のソースの最後に
hn = hn - tu
を加えてみました。すると、
イメージ通り、粒子同士は結着したものの、転がり抵抗力に余計な影響を与えるらしく、結着したまま回り続ける謎の動きにorz
転がり抵抗をOFFってみると、こんな感じ。
ありゃー、結着しながら転がっていくー。
転がり抵抗のモーメント算定式を修正するのと他に、粒子同士が結着した接触点が固定されるように回転を考慮しなくはいけないことが判明。
上記論文、モデルについて丁寧に書かれてくれている部類の論文だと思いますが、いざコード化すると、色々と不明な点が出てきます。さて、どうしようかな・・・。