流れに従ひて己を失はず．

日々の研究について書きます

ロボットビジョンの基礎(コロナ社)を読んだ

Optics

今更ながら読みました．

エピポーラ幾何/カメラのR,t推定/カルマンフィルタ...
中途半端にわかった気になってわかっていなかったことが多すぎたのでこの機会に整理できてよかったという感想です．

読んだ目的は，これまでずっと回避してしまっていたカメラキャリブレーション周りの勉強．
基本的には外部パラメータ(カメラの並進と回転)を求めるところがキャリブレーションの内容としてメインで書かれています．途中式など詳細に解説されていて非常に勉強になりました．
ただ，内部パラメータ(収差や焦点距離などなど)を求めるところはさらっと流して書いてある程度でこの部分が個人的に結構気になる部分だったので，次はそこの勉強です．

学会発表準備

社会人博士課程 論文投稿 学会口頭発表

サンフランシスコでのSPIEという学会にて発表してきました．
自分の中で発表準備が過不足なくできたな，と思えたのでここにメモを残します．
SPIE Photonics West


原稿日本語版完成 2019/12/11
英語にすること前提で日本語で書いていく．

原稿完成+投稿 2020/01/08
英文化して，校閲出して，レビューしてもらって完成．

発表スライド骨子完成 2020/01/19
発表スライドのタイトルとそのスライドを端的に表す1文or2文(max20words)を書き出す．

発表スライド絵コンテ完成 2020/01/25
図や表の配置と箇条書きの文の配置を大雑把に決める．大まかに発表台本をイメージする．タイトルと端的表現文も直しながら．

発表スライド完成 2020/01/28
絵コンテで詰めきれていない細かいところを埋めていく．スライド全体のバランスを考えながら文字を足すor削る．

発表3日前 2020/02/02
台本の前半を作る．接続詞とかスライドとスライドの間や文章と文章の間のつなぎを考える．スライドにある単語や文章は，台本中で赤線を引いておく．→ 赤線が引かれていないところが覚えるべきところ．

発表2日前 2020/02/03
台本の後半を作る．後半は実験結果の内容がメインなので台本作成が楽．台本を横に置きつつ，スライドを見ながらの練習2回．スライドとスライドの間のつなぎの文句とかを覚える．

発表1日前 2020/02/04
移動中の行きと帰りに1回づつ頭の中で発表練習．スライドの内容を脳内再生しながら話す．覚えきれてない文章はあのスライドのあそこに書いてあるのを見ながらなら話せる，などを確認．気になった箇所は終わったあとで確認(練習中は一旦始めたら途中で止めない)．ホテルの部屋で立って話す練習2回．発音が不安な発表に関係する単語を思いつくだけ調べる．昔ICCPという学会で発表したときフレネルレンズと言われているのに全く聞き取れなかった．終わって10日くらいしてからやっとわかった．

発表当日 2020/02/05
朝起きて1回練習．会場についてから外のベンチで1回練習．発音が不安なところは追加で練習．自分のセッションは他の人の発表を受けての質問があるかもしれないので予習もしてしっかり聞く．(先のフレネルレンズの話はICCPでの自分の発表の1つ前の発表で話されていた内容だったが，当時現場にいた自分は全然聞き取れていなくて理解できていなかった．)

社会人博士課程学生周辺の環境について情報共有

社会人博士課程

社会人博士課程学生である私の周りの環境についての話．
自分の振り返りを兼ねて書きます．
この記事は社会人学生 Advent Calendar 2019の12日目の記事です．
adventar.org


・概要
大学院：情報理工学系研究科システム情報学専攻．2017年9月入学，2019年12月現在でD3．光学デバイスを用いた計測及び情報提示システムの研究に従事．
会社：2016年4月入社，2019年12月現在で4年目．職務内容は情報処理周辺の研究開発業務．

・時系列
2016年3月　修士修了
2016年4月　入社
2017年1月　博士号を意識&先輩各位から社会人博士のアドバイス
2017年3月　先生に相談
2017年6月　出願書類提出
2017年8月　入試
2017年9月　合格発表&入学

・会社のこと
職場は研究開発なので，博士号取得が推奨されています．自分の周りでも社会人博士に通っている人も何人かおり，直属の上司も通学中という状況，同期でも何人か通学してます．
社会人博士をサポートする支援制度も職場に存在しています．支援内容は，入学金・授業料・通学交通費の一部が補助されるという金銭的な面と，業務に支障が出ない範囲での平日日中の通学が認められる(大体週1日か半日というケースが多い)という時間的な面の2つです．支援対象となるにはいくつか基準があり，入社後3年半以上経過していること，入社後一定の研究成果があること，などがあります．自分の場合は入社(2016年4月)から1年半後の2017年秋入学だったため，3年半経過後の2019年10月から支援対象になり，現在は金銭面と時間面のサポートを得ました．
職場で周りの人に相談しても，ほとんどがポジティブな反応でしたが，一方で年次が若く前例が無いことを理由に反対・心配される事もありました．まあそれは2年目の終わりにある研修員発表（この苦労を経験してやっと一人前になったな！うん！というイベント）があるからだったのだろうな，と今から考えると思います．実際のところ，私見ですが大変なこともいろいろあるけれど，年次が進むと職務に加えて役職・委員の仕事が増える＆ライフステージの変化があるので，行くなら若い時期から進学したほうが良いと思います．
※ここらへんの，会社の制度や待遇，そして平日何日くらい大学に行けているかの情報がもっとネットに集まってくるといいなあと考えてます．

・大学のこと
私は出身研究室に再び入れていただくこととなりました．先生にご相談させていただいてから入学までの半年間は，入試の準備と入学後の研究の下準備を進めていました．6月提出の出願書類にすでに研究計画書を入れなければならず，ビジョンが固まっていない状態でしたので作成に苦労しました．あとは，社会人博士なので「職場としてこの人を通学させていいと考えてる」ことを証明する書類を会社に作ってもらって出願時に大学に提出しました．
入試は専門科目と英語，研究内容のプレゼンがありますが，自分の場合修士の出身研究科ということで専門科目と英語は免除でした．研究内容のプレゼンの内容は修士の研究紹介と博士課程での研究計画です．出願書類でも大変でしたが研究計画プレゼンの準備が結構大変でした．あんまりできるビジョンが見えていすぎても駄目だし，全くできる気がしなくても駄目なのでそのへんの塩梅が難しかったです．先輩方からは「大丈夫だよ〜」と言われてましたが，プレゼン前は結構緊張しました．落ちる人は落ちます．
8月23日にプレゼンがあって，9月4日に合格発表でした．プレゼン後の質疑の雰囲気は良かったのでそこまで心配はしていなかったんですが，掲示をみてやっと一安心という感じでした．
社会人博士が受けられる金銭面のサポートとして，博士課程研究遂行協力制度(参考：https://giju.t.u-tokyo.ac.jp/eejim/wp-content/uploads/2018/05/web_guide_suiko_20180531.pdf)というものがあります．半年間月5万円支給なので学費に充てることができ，とても助かりました．
授業は，博士課程研究以外に通算4コマ8単位を取らなければなりません．修士のときに余分に単位をとっていれば博士課程での修了要件の単位数に加えることができたんですが，私はとっていなかったので講義をとる必要があります．大学を出た後にとっておけばよかったなと思っていた，他分野の講義をとるようにしました．
本来は博士課程学生がやらなかればならない研究室の仕事を免除していただいており，研究室の関係者各位には頭が上がりません．

・普段の生活
土日休み，平日4日会社，1日大学です．
夜は自宅で毎日だいたい21時〜24時で作業します．大学に行く日には装置を組んだりこれを使った実験をしてます．シミュレーションや調べ物，コーディングは自宅からやってます．基本的にゲームやったり工作したりなどの時間は無いですが，週2回は身体を動かす時間を確保できている状況です．

ベクトル正規化(3次元)の速度比較

Python

np.linalg.norm()を使うと遅い．
二乗和はぐりぐり計算して，平方根はmath.sqrt()を使うのが良い．
という話．


np.linalg.norm()を使うのと，np.dot()を使うのと，math.sqrt()を使うのと，愚直に計算するのではどれが一番速いか比べた．

ライブラリインポートとデータ作成

import numpy as np
import math

vec = np.random.rand(1000,3)


np.linalg.norm()を使った場合:9.32 ms

def vec_3d_length_np(v):
    return v/np.linalg.norm(v)

%%time
for _ in range(1000):
    vn = vec_3d_length_np(vec[_,:])

CPU times: user 8.94 ms, sys: 3.39 ms, total: 12.3 ms
Wall time: 9.32 ms


np.dot()を使った場合:7.08 ms

def vec_3d_length_npdot(v):
    return v/(np.dot(v,v))**(0.5)

%%time
for _ in range(1000):
    vn = vec_3d_length_npdot(vec[_,:])

CPU times: user 6.08 ms, sys: 4.31 ms, total: 10.4 ms
Wall time: 7.08 ms


math.sqrt()を使った場合:2.35 ms

def vec_3d_length_math(v):
    return v/math.sqrt(v[0]*v[0] + v[1]*v[1]+v[2]*v[2])

%%time
for _ in range(1000):
    vn = vec_3d_length_math(vec[_,:])

CPU times: user 2.33 ms, sys: 6 µs, total: 2.33 ms
Wall time: 2.35 ms


愚直に計算した場合:2.73 ms

def vec_3d_length_(v):
    return v/(v[0]**2 + v[1]**2+v[2]**2)**(0.5)

%%time
for _ in range(1000):
    vn = vec_3d_length_(vec[_,:])

CPU times: user 2.72 ms, sys: 3 µs, total: 2.72 ms
Wall time: 2.73 ms

スネルの法則を使って屈折方向の計算

Python Optics

スネルの法則で3次元空間上の屈折方向ベクトルを求めるのをPythonで実装．


ベクトルで解く

各要素の関係は以下．ただし各ベクトルの大きさを1とする．
入射ベクトル：v_{inter}
法線ベクトル：v_{norm}
出射ベクトル(これを求める)：v_{exit}
入射角：\theta_{inter}
出射角(これもわかってない)：\theta_{exit}
入射側屈折率：n_{inter}
入射側屈折率：n_{exit}

f:id:Yasutchi:20190914234142j:plain


スネルの法則
n_{exit} sin\theta_{exit} = n_{inter} sin\theta_{inter}. \tag{0}
つまり，\theta_{exit}はぐりぐりやればわかる．


ここで，入射ベクトル，法線ベクトル，出射ベクトルは同一平面上にあるから， \alpha, \beta(ともにスカラ)を使って以下のように表現できる．

v_{exit} = \alpha\ v_{inter} + \beta\ v_{norm}. \tag{1}

まず \alphaから求める．
式(1)の両辺について，法線ベクトルとの外積をとって，更にその大きさをとると，各ベクトルの大きさが1であること，最後は式(0)を用いて\alphaが決まる．
\begin{eqnarray} |v_{norm} \times v_{exit}| &=& |v_{norm} \times (\alpha\ v_{inter} + \beta\ v_{norm})|,\\|v_{norm} \times v_{exit}| &=& |v_{norm} \times \alpha\ v_{inter} + v_{norm} \times \beta\ v_{norm}|,\\|v_{norm}||v_{exit}|sin(\pi - \theta_{exit}) &=& \alpha|v_{norm}||v_{inter}|sin(\pi - \theta_{inter}) + 0,\\sin\theta_{exit} &=& \alpha sin\theta_{inter},\\ \alpha &=& \dfrac{sin\theta_{exit}}{sin\theta_{inter}} = \dfrac{n_{inter}}{n_{exit}}. \end{eqnarray}



次に \betaを求める．
式(1)の両辺について，法線ベクトルとの内積をとると，先程と同様に各ベクトルの大きさが1であることを用いて\betaが決まる．\theta_{exit}が残っているが，式(0)から求めることができる．
\begin{eqnarray} v_{norm} \cdot v_{exit} &=& v_{norm} \cdot (\alpha\ v_{inter} + \beta\ v_{norm}),\\v_{norm} \cdot v_{exit} &=& \alpha\ v_{norm} \cdot v_{inter} + \beta\ v_{norm} \cdot v_{norm},\\ |v_{norm}||v_{exit}|cos(\pi - \theta_{exit}) &=& \alpha|v_{norm}||v_{inter}|cos(\pi - \theta_{inter}) + \beta,\\-cos\theta_{exit} &=& -\alpha cos\theta_{inter} + \beta,\\ \beta &=& \alpha cos\theta_{inter} -cos\theta_{exit}. \end{eqnarray}



Python実装

import numpy as np 

def snell(v_inter,v_norm,n_inter,n_exit):
    #3次元スネルの公式の適用結果を返す関数
    #入力：入射光線の方向ベクトル，入射点での法線ベクトル，入射元の屈折率，入射先の屈折率
    
    #入力されたベクトルの正規化
    v_inter = v_inter / np.linalg.norm(v_inter)
    v_norm = v_norm / np.linalg.norm(v_norm)
    
    #係数a
    a =  n_inter/n_exit
    
    cos_t_inter = -1.0*np.dot(v_inter,v_norm)
    #量子化誤差対策
    if cos_t_inter<-1.:
        cos_t_inter = -1.
    elif cos_t_inter>1.:
        cos_t_inter = 1.
        
    sin_t_inter = (1.0 - cos_t_inter**2)**(1/2)
    sin_t_exit = sin_t_inter*a
    
    if sin_t_exit>1.0:
        #全反射する場合
        return np.zeros(3), False
    
    cos_t_exit = (1 - sin_t_exit**2)**(1/2)
    
    #係数b
    b = -1.0*cos_t_exit + a*cos_t_inter
    
    #出射光線の方向ベクトル
    v_exit = a*v_inter + b*v_norm
    #正規化
    v_exit = v_exit / np.linalg.norm(v_exit)
    
    return v_exit, True

vi = np.array([0.08584,0.17301,0.9811726])
vn = np.array([0.050, 0.060, -0.9969453])
ni = 1.0
ne = 1.5

snell(vi,vn,ni,ne)#(array([0.0401461 , 0.0948433 , 0.99468238]), True)

#方向ベクトルなどの実値参考：
#https://www.cybernet.co.jp/optical/course/optics/opt02/opt06.html


参考：http://yees.g.dgdg.jp/RayTrace/RayTraceVersion001b.pdf

SymPyでの幾何計算速度

Python

pythonのライブラリであるSymPyを使って幾何計算をやるととても遅いという話．

以下は2直線の交点導出を100セット解いた例．
SymPyを使った場合とベクトルで解いた場合の比較．10^3オーダーで変わってくる．

SymPyを使った場合

import numpy as np
import sympy.geometry as sg
import time

#点［listA］を通り方向ベクトル[listB]を持つ直線と
#点［listC］を通り方向ベクトル[listD]を持つ直線との交点を求める
listA = np.random.rand(100,2)
listB = np.random.rand(100,2)
listC = np.random.rand(100,2)
listD = np.random.rand(100,2)

start = time.time()

for i in range(0,100):
    pointA = listA[i,:]
    direcB = listB[i,:]
    pointC = listC[i,:]
    direcD = listD[i,:]
    line1 = sg.Line( sg.Point(pointA), 
                                        sg.Point(pointA+direcB))
    line2 = sg.Line( sg.Point(pointC), 
                                        sg.Point(pointC+direcD))
    result = sg.intersection(line1, line2)

elapsed_time = time.time() - start
print ("sympy calc_time:{0}".format(round(elapsed_time,3)) + "[sec]")

sympy calc_time:15.911[sec]


ベクトルで解いた場合

import numpy as np
import sympy.geometry as sg
import time

#点［listA］を通り方向ベクトル[listB]を持つ直線と
#点［listC］を通り方向ベクトル[listD]を持つ直線との交点を求める
listA = np.random.rand(100,2)
listB = np.random.rand(100,2)
listC = np.random.rand(100,2)
listD = np.random.rand(100,2)

#90度の回転行列
rot90 = np.array([[np.cos(np.pi/2), -np.sin(np.pi/2)],
                              [np.sin(np.pi/2), np.cos(np.pi/2)]])

def calcintersection(point1,direc1,point2,direc2):
    # 参考
    #　https://qiita.com/tmakimoto/items/2da05225633272ef935c
    
    #2直線が平行なら[]を返す
    if direc1[0]*direc2[1]==direc1[1]*direc2[0]:
        return []
    
    #方向ベクトルの正規化
    direc1_n = direc1 / np.linalg.norm(direc1) 
    direc2_n = direc2 / np.linalg.norm(direc2) 
    
    #direc2_nに垂直な単位ベクトル
    direc3_n = np.dot(direc2_n, rot90)
    
    #直線1における交点のパラメータloを求める
    lo = np.dot(direc3_n,(point2-point1)) / np.dot(direc3_n,direc1_n)
    
    return point1+direc1_n*lo


start = time.time()

for i in range(0,100):
    pointA = listA[i,:]
    direcB = listB[i,:]
    pointC = listC[i,:]
    direcD = listD[i,:]
    
    result = calcintersection(pointA,direcB,pointC,direcD)

elapsed_time = time.time() - start
print ("vector calc_time:{0}".format(round(elapsed_time,3)) + "[sec]")

vector calc_time:0.004[sec]

pythonのmultiprocessingで複数引数複数戻値を並列処理

Python multiprocessing

毎回忘れるのでメモ．

pythonで複数引数をとって，複数戻値しかも配列を返す処理を並列で行う．

参考: https://qiita.com/wikipediia/items/2919362de582a7d8de9e

import time
import numpy as np
import os
from multiprocessing import Pool

#総スレッド数の取得
cpu_count = os.cpu_count()

#最終的に作りたい出力結果の行数
df_size = 1234

#作りたいデータその１
#df_size*2の２次元データ
#n行目のデータは[n,2*n+1]にしたい
data_frame = np.ones([df_size,2],dtype=int)

#作りたいデータその2
#df_size*2の２次元データ
#n行目のデータは[n,n^2+1]にしたい
data_frame_s = np.ones([df_size,2],dtype=int)

#それぞれのスレッドに割振る行数の開始番号と終了番号の組
sta_gol = np.ones([cpu_count,2],dtype=int)
for i in range(0,cpu_count):
    sta = (df_size//cpu_count)*i
    gol = (df_size//cpu_count)*(i+1)-1
    if i == cpu_count-1:
        gol = df_size-1
    sta_gol[i,:] = [sta,gol]
print(sta_gol)

# ここを並列処理する
def func(i,arg1, arg2):
    sta = int(arg1)
    gol = int(arg2)
    print("Thread for sta:",sta," gol:", gol)
    
    #作りたいデータの割り当てられた各部分を作る処理
    data_frame_sg = np.ones([gol+1-sta,2],dtype=int)
    data_frame_sg_s = np.ones([gol+1-sta,2],dtype=int)
    for r in range(sta,gol+1):
        data_frame_sg[r-sta,:] = [r,2*r+1]
        data_frame_sg_s[r-sta,:] = [r,r*r+1]
    
    #ダミーの処理
    time.sleep(5/(i+1))
    
    return i,data_frame_sg,data_frame_sg_s
    
def wrapper(args):
    # argsは(i,sta, gol)となっている
    return func(*args)

def multi_process(sampleList):
    # プロセス数:cpu_count
    p = Pool(cpu_count)
    #20220718追記
    # M1Macだとp = get_context("fork").Pool(cpu_count)
    # 参考：https://stackoverflow.com/questions/67999589/
    
    #ラッパにわたす
    output = p.map(wrapper, sampleList)
    
    # プロセスの終了
    p.close()
    return output

def mproc():
    # (i,sta, gol)が引数になる
    sampleList = [(i,sta_gol[i,0],sta_gol[i,1]) for i in range(cpu_count)]

    return multi_process(sampleList)

start = time.time()

#実行するところ
mproc_output = mproc()

elapsed_time = time.time() - start
print ("elapsed_time:{0}".format(round(elapsed_time,2)) + "[sec]")

#並列処理の結果をもとのスキーマに直すところ
for i in range(0,len(mproc_output)):
    ind = mproc_output[i][0]
    data_frame[sta_gol[ind,0]:sta_gol[ind,1]+1,:] = mproc_output[i][1]
    data_frame_s[sta_gol[ind,0]:sta_gol[ind,1]+1,:] = mproc_output[i][2]
print(data_frame)
print(data_frame_s)