枝分かれの角度が一定の枝の平面射影図上での角度分布¶

Hiroyuki OGAWA (Osata Univ.)　　　 2023/02/21¶

1. Introduction¶

シャーレ内で成長する金属樹の形を調べる方法として、成長過程で写真を撮り、写真に写った枝分かれ部分の角度（観測角と呼ぶ）を計測し、観測角の分布を求めることができる。金属樹は３次元的に広がる溶液内で成長するので、枝分かれ部分の本来の（３次元空間内での）角度を知るには、例えば立体視など、３次元像を認識できる形での観測をする必要がある。繊細な金属樹は動かすと壊れてしまう。観測機材などの問題もあるが、観察している金属樹を動かさずに異なる向きからの撮影はできず、３次元像を得ることは難しい。一方向からの写真像だけをもとに、金属樹の枝分かれ部分の角度を知りたい。そこで、３次元空間内にランダムに配置された枝の２次元射影像における角度の分布と、実験で写真像をもとに計測した観測角の分布を比較することで、金属樹の枝分かれ部分の角度を推定できないであろうか。以下、３次元空間内にランダムに配置された枝の２次元射影像における角度の分布を計算する。

1.1. 仮定とモデルと言葉の定義¶

・シャーレは水平におき、シャーレ面に垂直な光源でシャーレ面に平行な面に投影される
　・シャーレ面と平行な写真面上に主枝、分枝が投影される
　・主枝、分枝は直線で、分枝は主枝の伸びる方向に対し、一定角（$\theta$）で枝分かれする。$\theta$ を「枝分かれ角」と呼ぶ
　・シャーレ面と平行に $xy$-座標を定め、シャーレ面に垂直な向きに z-座標を定める
　・主枝は原点を通り、主枝は $y$-座標が $0$ の $xz$-平面内を伸びるとする
　・主枝とシャーレの面のなす角（主枝と $x$-軸のなす角）を $\phi$ とし、「主枝の傾斜角」と呼ぶ
　・分枝と主枝を含む平面と、主枝と $y$-軸を含む平面のなす角を $\psi$ とし、「枝分かれ面の傾斜角」と呼ぶ
　・写真面上の $x$-軸、$y$-軸の投影像を使って、写真面上に $xy$-座標をいれる
　・主枝は、写真面上の $x$-軸上に投影される
　・分枝の投影像は、ベクトル
　　　　　$\alpha = (-{\rm sin}(\phi)\, {\rm sin}(\psi)\, {\rm sin}(\theta) + {\rm cos}(\phi)\, {\rm cos}(\theta), {\rm sin}(\theta)\, {\rm cos}(\psi))$
　　と平行である
　・写真面における主枝と分枝の投影像の間の角を「投影像の角」と呼ぶ
　・ここで考えた３次元空間内の枝配置モデルは、枝分かれ角（$\theta$）と主枝の傾斜角（$\phi$）をパラメータに持つ
　・主枝の傾斜角の変動幅（$K$ : $|\phi|\leqq K$）を指定する。分枝面の傾斜角 $\psi$ 変動幅傾斜角の組 $(\psi, \phi)$ を独立な確率変数とし、$x$-軸に対するベクトル $\alpha$ の角の分布を、「投影像の角の分布」と呼ぶ
　・つまり、枝分かれ角（$\theta$）と主枝の傾斜角（$\phi$）の変動幅（$K$ : $|\phi|\leqq K$）をパラメータとして、投影像の角の分布が定まる
　・実際の実験で撮影した写真上の主枝と分枝の像のなす角を（たくさん）測定する。これらを「観測角」と呼び、「観測角の分布」を調べる

1.2. 目的¶

実験により得られた「観測角の分布」と、このモデルにおける「投影像の角の分布」を比較し、金属樹における「枝分かれ角（$\theta$）」を推定したい。

1.3. 注意事項¶

分枝からさらに枝分かれした枝に対しても投影像における角度を計測して観測角の分布に含める場合は、少し注意が必要である。主枝の傾斜角は水平面（シャーレ面と平行）に近いと仮定しても、短い分枝の水平面に対する傾きの変動幅は大きくなり、分枝から枝分かれした枝の水平面に対する傾きの変動幅はさらに大きくなる。従って、主枝から枝分かれした第１分枝だけでなく、第１分枝から枝分かれした第２分枝、第２分枝から枝分かれした第３分枝と、何段階も枝分かれした後のものも計測するなら、主枝の傾斜角の変動範囲が上面から下面までの $\pm 90^\circ$ のモデルでの分布と比較すべきであろう。

2. procedures¶

2.1. 必要なモジュール類の組み込み¶

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import pyplot as plt
import ipywidgets as widgets

2.2. small procedures¶

小数点以下の桁数調整¶

In [2]:

def chop(*args, prec=2):
    if len(args) == 0:
        return(False)
    elif len(args) == 1:
        return(float(np.floor(10**prec*args[0]+1/2)/10**prec))
    return(tuple([chop(x, prec=prec) for x in args]))

リストの転置¶

In [3]:

def listTranspose(a):
    return([list(x) for x in zip(*a)])

区間内のデータ数え上げ¶

In [4]:

def countInt(lst, rng):
    return(len([x for x in lst if rng[0]<=x<rng[1]]))

２次元データにおける、数え上げ手続き¶

In [5]:

def selInt(data, r0, r1):
    return(sum([[(i,j) for j,x in enumerate(d) if r0<=x<r1] 
                for i,d in enumerate(data)], []))

def selInt2(data, r0, r1):
    return([[j for j,x in enumerate(d) if r0<=x<r1] 
            for i,d in enumerate(data)])

def countArg2(data, r0, r1):
    return(len(selInt(data, r0, r1)))

def pos2arg(arg, pos):
    return(arg[0][pos[0]][pos[1]], arg[1][pos[0]][pos[1]])

def pos2arg2(arg, pos2):
    return(sum([[(ags[j],agh[j]) for j in js] 
                for js,ags,agh in zip(pos2, *arg)], []))

角度　小数点以下切り捨て¶

In [6]:

def int180(x):
    return(int(abs(x))%180)

def int90(x):
    return(int(90.5-abs(90-abs(x))))
    #return(int(90-abs(90-abs(x)))%90 if x!=90 else 90)

定型長表示のため¶

In [7]:

def strL(x, r=5):
    return((' '*r+'{}'.format(x))[-r:])
def strR(x, r=5):
    return(('{}'.format(x)+' '*r)[:r])

2.3. 結果出力¶

printStat : モデルにおける観測角分布の平均と標準偏差の表出力¶

In [8]:

def printStat(K=False, t=False, s='all', sw=True):
    #sw = True if K or not(t) else False
    if not(K): K=Ks
    if not(t): t=thetas
    if type(K) not in (list, tuple): K = list(K)
    if type(t) not in (list, tuple): t = list(t)
    print('θ : 枝分かれ角    φ : 主枝の傾斜角    ψ : 分枝面の傾斜角\n'
          + 'K : |φ| ≦ K     ave : 平均     std : 標準偏差\n')
    if sw:
        for K0 in K:
            print('R = {}'.format(K0))
            print('  θ   :', *[strL(t0) for t0 in t])
            print('  ave :', *[strL(stat_arg_of_edges[K0][t0][0]) 
                               for t0 in t])
            print('  std :', *[strL(stat_arg_of_edges[K0][t0][1]) 
                               for t0 in t])
    else:
        for t0 in t:
            print('θ = {}'.format(t0))
            print('  K   :', *[strL(K0) for K0 in K])
            print('  ave :', *[strL(stat_arg_of_edges[K0][t0][0]) 
                               for K0 in K])
            print('  std :', *[strL(stat_arg_of_edges[K0][t0][1]) 
                               for K0 in K])

plotArgs : 視覚化手続き¶

主枝の傾斜角（$\phi$）の変動幅（$K$:$|\phi|\leqq K$）を引数とする。 枝分かれ角（$\theta$）と、投影像のなす角の含まれる区間をスライダーで指定するたびに、それを与える傾斜角の組 $(\psi, \phi)$（$\psi$:枝分かれ面の傾斜角、$\phi$:主枝の傾斜角）を図示する

In [9]:

def plotArgs(K, figsize=(4,4), xRange=(-180,180), ptsize=1, 
             aspect='equal', **kwargs):
    #
    def plotLocalProcedure(t, r, d):
        theta = min(90, max(5, 5*int(t/5)))
        #
        fig = plt.figure(figsize=figsize)
        ax = fig.add_subplot(111)
        ax.set_xlim(*xr)
        ax.set_ylim(*yr)
        ax.set_aspect(aspect)
        ax.set_xticks(xLp, xL)
        ax.set_yticks(yLp, yL)
        ax.set_xlabel('arg of sub-edge')
        ax.set_ylabel('arg of main-edge')
        #ax.axis("off")
        #
        pos = listTranspose(pos2arg2(
            psh, selInt2(data[theta], r, r+d)))
        if len(pos): ax.scatter(*pos, s=ptsize)
        msg.value = '{} ({}%)'.format(
            len(pos[0]), 
            chop(100*len(pos[0])/m)) if len(pos) else '0 (0%)'
        #
        #ax.clear()
        fig.canvas.draw_idle()
        return(fig)
    #
    psh, data = arg_of_edges_K[K]
    m = len(np.ravel(psh[0]))
    #
    yRange = xRange[0]*K/180, xRange[1]*K/180
    xD = 10*int((xRange[1]-xRange[0])/8/10)
    yD = 5*int((yRange[1]-yRange[0])/6/5)
    xr = xRange[0]*np.pi/180, xRange[1]*np.pi/180
    yr = yRange[0]*np.pi/180, yRange[1]*np.pi/180
    xL = [s for s in range(-180,190,10) 
          if xRange[0]<=s<=xRange[1] and s%xD==0]
    yL = [s for s in range(-180,190,5) 
          if yRange[0]<=s<=yRange[1] and s%yD==0]
    xLp = [np.pi/180*s for s in xL]
    yLp = [np.pi/180*s for s in yL]
    #
    # sliders
    sliderT = widgets.IntSlider(
        value=90, min=5, max=90, step=5, description='枝分かれ角 :', 
        layout=widgets.Layout(width='182pt'))
    sliderR = widgets.IntSlider(
        value=20, min=0, max=179, step=1, description='投影像の角 :', 
        layout=widgets.Layout(width='231pt'))
    sliderD = widgets.FloatLogSlider(
        value=1, min=-1, max=2.2, step=0.1, description='幅 :', 
        layout=widgets.Layout(width='231pt'))
    #
    #button = widgets.Button(description='保存', 
    #                        layout=widgets.Layout(width='42pt'))
    #text = widgets.Text(description='ファイル', 
    #                    layout=widgets.Layout(width='170pt'))
    msg = widgets.Text(description='度数 :', 
                       layout=widgets.Layout(width='215pt'))
    #
    #def saveButton(b):
    #    fig = plotLocalProcedure(
    #        sliderM.value, sliderN.value, sliderD.value, sliderTr.value)
    #    if saveFigure(fig, savefile=text.value):
    #        msg.value = 'saved to '+text.value
    #    else:
    #        msg.value = 'error: cannot save to '+text.value
    #
    #button.on_click(saveButton)
    #
    uiAll = widgets.VBox([sliderT, sliderR, sliderD, msg])
    #uiAll = widgets.VBox(
    #    [sliderT, sliderR, sliderD, msg, text, button])
    out = widgets.interactive_output(
        plotLocalProcedure, {'t': sliderT, 'r': sliderR, 'd': sliderD})
    display(widgets.HBox([out, uiAll]))

freq : 度数¶

主枝の傾斜角の変動範囲（$K$）、枝分かれ角（$theta$）に対し、投影像の角がある区間に属するときの、傾斜角 $(\phi,\psi)$ の度数

In [10]:

def freq(K=90, theta=90, interval=(20,21)):
    if K not in Ks:
        return('error : {} is not in Ks.'.format(K))
    t = max(5, min(90, 5*np.floor(theta/5)))
    if type(interval) not in (list, tuple):
        return(countInt(aes[K][t], (interval,interval+1)))
    elif len(interval)<1:
        return('error : specify an interval')
    elif len(interval) == 1:
        return(countInt(aes[K][t], (interval,interval+1)))
    else:
        return(countInt(aes[K][t], interval[0:2]))

In [ ]:

3. 計算¶

3.1. 投影像の角の大きさの分布の計算¶

傾斜角の組 ($\psi$,$\phi$)（$\psi$:枝分かれ面の傾斜角、$\phi$:主枝の傾斜角）に対して、「投影像の角の分布」を計算する。
　　r : 角度１度の細分
　　thetas : 枝分かれ角 $\theta$ の組 : $5^\circ \leqq\theta\leqq 90\circ$（$5^\circ$ 刻み）
　　Ks : 主枝の傾斜角 $\phi$ の変動幅 $K$（$|\phi|\leqq K$）の組
　　傾斜角の組 $(\psi, \phi)$ : $|\psi| \leqq 180^\circ$　$|\phi| \leqq 90^\circ$
　　arg_of_edges : 枝分かれ角（$\theta$）ごとの、傾斜角の組 $(\psi, \phi)$ に対する投影像の角のデータ
　　arg_of_edges_K : 主枝の傾斜角の変動幅（$K$ : $|\phi|\leqq K$）ごとに、傾斜角の組 $(\psi, \phi)$ の変動領域を制限してまとめ直した、投影像の角のデータ
　　aes : 度数分布計算のために、上記データを $K$, $\theta$ ごとに投影像の角度データを１次元配列にまとめたもの
　　stat_arg_of_edges : 枝分かれ角（$\theta$）と主枝の傾斜角（$\phi$）の変動幅（$K$ : $|\phi|\leqq K$）に対する、投影像の角の分布の統計量（平均と標準偏差）

In [11]:

freq_dist = {}
stat_arg_of_edges = {}

r = 5
m = 360 * r
m2, m4 = m//2, m//4

thetas = [5*i for i in range(1,19)]
Ks = (5, 10, 15, 20, 30, 45, 60, 90)

psi, phi = np.meshgrid(np.linspace(-np.pi,   np.pi,   m +1), 
                       np.linspace(-np.pi/2, np.pi/2, m2+1))

arg_of_edges = {theta : 180/np.pi*np.abs(
    np.arctan2(np.sin(np.pi/180*theta)*np.cos(psi), 
               -np.sin(phi)*np.sin(psi)*np.sin(np.pi/180*theta)
               +np.cos(phi)*np.cos(np.pi/180*theta)))
                for theta in thetas}

arg_of_edges_K = {K : ((psi[r*(90-K):r*(90+K)+1], 
                        phi[r*(90-K):r*(90+K)+1]), 
                       {theta : arg_of_edges[theta][r*(90-K):r*(90+K)+1] 
                        for theta in thetas}) 
                  for K in Ks}

aes = {K : {theta : np.concatenate(
    arg_of_edges[theta][m//360*(90-K):m//360*(90+K)+1]) 
            for theta in thetas} for K in Ks}

stat_arg_of_edges = {K : {theta : (chop(np.average(aes[K][theta])), 
                                   chop(np.std(aes[K][theta]))) 
                          for theta in thetas} for K in Ks}

3.2. 枝分かれ角（$\theta$）と主枝の傾斜角（$\phi$）の変動幅（$K$ : $|\phi|\leqq K$）を指定したモデルにおける、投影像の角の分布（度数分布（％）を計算する）¶

freq_dist : 主枝の伸びる向きに対する分枝の角度（$0^\circ\sim180^\circ$）の $1^\circ$ 刻みの度数分布（％）

枝分かれ角（$\theta$）が鈍角（$90^\circ\sim180^\circ$）の場合は計算していないが、対称性より鋭角（$0^\circ\sim90^\circ$）の範囲で十分である。枝分かれ角が $90^\circ$ 以下の場合、投影像が $180^\circ$ になるのは、特殊な状況（枝分かれ角 $\theta$ と主枝の傾斜角 $\phi$ が $\theta+\phi\geqq 180^\circ$ で、分枝面の傾斜角 $\psi=90^\circ$ の場合）に限るので、投影像の角は $0^\circ$ 以上 $180^\circ$ 未満と考えて十分である。度数分布を計算するための区間は、$s^\circ$ 以上 $s{+}1^\circ$ 未満（$s=0,1,\cdots,179$）とする

In [12]:

freq_dist = {K : {} for K in Ks}
Kss = [r*K for K in Ks]
for theta in thetas:
    wk = [0]*180
    for x in arg_of_edges[theta][m4]:
        wk[int180(x)] += 1
    for j in range(1,m4+1):
        for x in np.concatenate([arg_of_edges[theta][m4+j], 
                                 arg_of_edges[theta][m4-j]]):
            wk[int180(x)] += 1
        if j in Kss:
            n = m*(2*j+1)/100
            freq_dist[j//r][theta] = list(
                enumerate([x/n for x in wk])).copy()

3.3. 枝分かれ角（$\theta$）と主枝の傾斜角（$\phi$）の変動幅（$K$ : $|\phi|\leqq K$）を指定したモデルにおける、投影像の角の分布（度数分布（％）を計算する）¶

freq_dist90 : 主枝（伸びる向きは考慮しない）に対する分枝の角度（$0^\circ\sim180^\circ$）の $1^\circ$ 刻みの度数分布（％）

度数分布を計算するための区間は、$s-0.5^\circ$ 以上 $s+0.5^\circ$ 未満（$s=0,1,\cdots,90$）とする。両端 $0^\circ$ と $90^\circ$ では、有効な区間が半分になるので、度数を２倍にしている。丁度 $0^\circ$、丁度 $90^\circ$ のデータの度数まで２倍にしているので、この部分の計算は少し修正する必要はある。実際の観測との比較において、両端部分のこの程度の違いが影響を与えるとは思えないので、とりあえずはこのままにしておく。

In [13]:

freq_dist90 = {K : {} for K in Ks}
Kss = [r*K for K in Ks]
for theta in thetas:
    wk = [0]*91
    for x in arg_of_edges[theta][m4]:
        wk[int90(x)] += 1
    for j in range(1,m4+1):
        for x in np.concatenate([arg_of_edges[theta][m4+j], 
                                 arg_of_edges[theta][m4-j]]):
            wk[int90(x)] += 1
        if j in Kss:
            n = (m*(2*j+1)+wk[0]+wk[90])/100
            freq_dist90[j//r][theta] = list(enumerate([x/n for x in wk]))
            freq_dist90[j//r][theta][0] = (
                0, 2*freq_dist90[j*360//m][theta][0][-1])
            freq_dist90[j//r][theta][-1] = (
                90, 2*freq_dist90[j*360//m][theta][-1][-1])

In [ ]:

4. 結果¶

4.1. 統計データ : 投影像の角の平均値、標準偏差 の比較¶

主枝の傾斜角（$\phi$）の変動幅（$K$ : $|\phi|\leqq K$）ごとの表¶

In [14]:

printStat()

θ : 枝分かれ角    φ : 主枝の傾斜角    ψ : 分枝面の傾斜角
K : |φ| ≦ K     ave : 平均     std : 標準偏差

R = 5
  θ   :     5    10    15    20    25    30    35    40    45    50    55    60    65    70    75    80    85    90
  ave :  3.19   6.4  9.64 12.93 16.29 19.74 23.31 27.01 30.89 34.98 39.33  44.0 49.07 54.67 60.98 68.31 77.48  90.0
  std :  1.54  3.08  4.61  6.14  7.65  9.14 10.61 12.04 13.42 14.74 15.96 17.04 17.92 18.49 18.57  17.8 15.24 13.94
R = 10
  θ   :     5    10    15    20    25    30    35    40    45    50    55    60    65    70    75    80    85    90
  ave :   3.2  6.42  9.68 12.98 16.36 19.82 23.41 27.13 31.03 35.15 39.53 44.25 49.39 55.08 61.54 69.29 79.21  90.0
  std :  1.55  3.09  4.63  6.16  7.68  9.18 10.66  12.1 13.49 14.82 16.06 17.16 18.07 18.69 18.85  18.3 18.75 19.22
R = 15
  θ   :     5    10    15    20    25    30    35    40    45    50    55    60    65    70    75    80    85    90
  ave :  3.22  6.47  9.74 13.07 16.47 19.96 23.57 27.33 31.27 35.44 39.88 44.68 49.93  55.8 62.66 71.15 80.42  90.0
  std :  1.56  3.11  4.66   6.2  7.74  9.25 10.74  12.2 13.62 14.97 16.24 17.38 18.35 19.07  19.5 20.97 22.47 23.04
R = 20
  θ   :     5    10    15    20    25    30    35    40    45    50    55    60    65    70    75    80    85    90
  ave :  3.25  6.53  9.83 13.19 16.63 20.16 23.82 27.62 31.62 35.86 40.39 45.31 50.75 56.99 64.55 72.77 81.31  90.0
  std :  1.57  3.14  4.71  6.27  7.82  9.35 10.87 12.36  13.8  15.2 16.52 17.74 18.83 19.85 21.93 24.12 25.58 26.09
R = 30
  θ   :     5    10    15    20    25    30    35    40    45    50    55    60    65    70    75    80    85    90
  ave :  3.34  6.71 10.11 13.57 17.11 20.77 24.56 28.52 32.71 37.18 42.04  47.5 53.91  60.8 67.94 75.22 82.59  90.0
  std :  1.62  3.25  4.87  6.48  8.09   9.7  11.3 12.88 14.45 16.01 17.58 19.29 22.13  25.1 27.56 29.34 30.41 30.77
R = 45
  θ   :     5    10    15    20    25    30    35    40    45    50    55    60    65    70    75    80    85    90
  ave :  3.58  7.18 10.83 14.56  18.4 22.38 26.56 31.02 35.91 41.51  47.4 53.42  59.5  65.6 71.71 77.81 83.91  90.0
  std :  1.78  3.56  5.35  7.16  8.98 10.84 12.76  14.8 17.15 20.77 24.43 27.59 30.19 32.24  33.8 34.88 35.52 35.74
R = 60
  θ   :     5    10    15    20    25    30    35    40    45    50    55    60    65    70    75    80    85    90
  ave :  4.01  8.07  12.2 16.46 20.93 25.75 31.18 36.77 42.36  47.9 53.36 58.75 64.07 69.33 74.54 79.71 84.86  90.0
  std :  2.16  4.35   6.6  8.94 11.47 14.49 19.05 23.36 26.97 29.94 32.37 34.34 35.92 37.15 38.08 38.72 39.11 39.23
R = 90
  θ   :     5    10    15    20    25    30    35    40    45    50    55    60    65    70    75    80    85    90
  ave : 10.49  18.1 24.67 30.59 36.04 41.14 45.95 50.54 54.93 59.17 63.27 67.27 71.19 75.03 78.82 82.57 86.29  90.0
  std : 20.18 26.38 30.37 33.24 35.43 37.14 38.51 39.62 40.53 41.27 41.88 42.37 42.76 43.07  43.3 43.46 43.56 43.59

枝分かれ角（$\theta$）ごとの表¶

In [15]:

printStat(sw=False)

θ : 枝分かれ角    φ : 主枝の傾斜角    ψ : 分枝面の傾斜角
K : |φ| ≦ K     ave : 平均     std : 標準偏差

θ = 5
  K   :     5    10    15    20    30    45    60    90
  ave :  3.19   3.2  3.22  3.25  3.34  3.58  4.01 10.49
  std :  1.54  1.55  1.56  1.57  1.62  1.78  2.16 20.18
θ = 10
  K   :     5    10    15    20    30    45    60    90
  ave :   6.4  6.42  6.47  6.53  6.71  7.18  8.07  18.1
  std :  3.08  3.09  3.11  3.14  3.25  3.56  4.35 26.38
θ = 15
  K   :     5    10    15    20    30    45    60    90
  ave :  9.64  9.68  9.74  9.83 10.11 10.83  12.2 24.67
  std :  4.61  4.63  4.66  4.71  4.87  5.35   6.6 30.37
θ = 20
  K   :     5    10    15    20    30    45    60    90
  ave : 12.93 12.98 13.07 13.19 13.57 14.56 16.46 30.59
  std :  6.14  6.16   6.2  6.27  6.48  7.16  8.94 33.24
θ = 25
  K   :     5    10    15    20    30    45    60    90
  ave : 16.29 16.36 16.47 16.63 17.11  18.4 20.93 36.04
  std :  7.65  7.68  7.74  7.82  8.09  8.98 11.47 35.43
θ = 30
  K   :     5    10    15    20    30    45    60    90
  ave : 19.74 19.82 19.96 20.16 20.77 22.38 25.75 41.14
  std :  9.14  9.18  9.25  9.35   9.7 10.84 14.49 37.14
θ = 35
  K   :     5    10    15    20    30    45    60    90
  ave : 23.31 23.41 23.57 23.82 24.56 26.56 31.18 45.95
  std : 10.61 10.66 10.74 10.87  11.3 12.76 19.05 38.51
θ = 40
  K   :     5    10    15    20    30    45    60    90
  ave : 27.01 27.13 27.33 27.62 28.52 31.02 36.77 50.54
  std : 12.04  12.1  12.2 12.36 12.88  14.8 23.36 39.62
θ = 45
  K   :     5    10    15    20    30    45    60    90
  ave : 30.89 31.03 31.27 31.62 32.71 35.91 42.36 54.93
  std : 13.42 13.49 13.62  13.8 14.45 17.15 26.97 40.53
θ = 50
  K   :     5    10    15    20    30    45    60    90
  ave : 34.98 35.15 35.44 35.86 37.18 41.51  47.9 59.17
  std : 14.74 14.82 14.97  15.2 16.01 20.77 29.94 41.27
θ = 55
  K   :     5    10    15    20    30    45    60    90
  ave : 39.33 39.53 39.88 40.39 42.04  47.4 53.36 63.27
  std : 15.96 16.06 16.24 16.52 17.58 24.43 32.37 41.88
θ = 60
  K   :     5    10    15    20    30    45    60    90
  ave :  44.0 44.25 44.68 45.31  47.5 53.42 58.75 67.27
  std : 17.04 17.16 17.38 17.74 19.29 27.59 34.34 42.37
θ = 65
  K   :     5    10    15    20    30    45    60    90
  ave : 49.07 49.39 49.93 50.75 53.91  59.5 64.07 71.19
  std : 17.92 18.07 18.35 18.83 22.13 30.19 35.92 42.76
θ = 70
  K   :     5    10    15    20    30    45    60    90
  ave : 54.67 55.08  55.8 56.99  60.8  65.6 69.33 75.03
  std : 18.49 18.69 19.07 19.85  25.1 32.24 37.15 43.07
θ = 75
  K   :     5    10    15    20    30    45    60    90
  ave : 60.98 61.54 62.66 64.55 67.94 71.71 74.54 78.82
  std : 18.57 18.85  19.5 21.93 27.56  33.8 38.08  43.3
θ = 80
  K   :     5    10    15    20    30    45    60    90
  ave : 68.31 69.29 71.15 72.77 75.22 77.81 79.71 82.57
  std :  17.8  18.3 20.97 24.12 29.34 34.88 38.72 43.46
θ = 85
  K   :     5    10    15    20    30    45    60    90
  ave : 77.48 79.21 80.42 81.31 82.59 83.91 84.86 86.29
  std : 15.24 18.75 22.47 25.58 30.41 35.52 39.11 43.56
θ = 90
  K   :     5    10    15    20    30    45    60    90
  ave :  90.0  90.0  90.0  90.0  90.0  90.0  90.0  90.0
  std : 13.94 19.22 23.04 26.09 30.77 35.74 39.23 43.59

In [ ]:

4.2. 特定条件下での傾斜角の組の集合¶

傾斜角の集合¶

主枝の傾斜角（$\phi$）の変動幅（$K$ : $|\phi|\leqq K$）と枝分かれ角（$\theta$）を指定し、投影像の角がある区間に属するときの傾斜角の組 $(\psi, \phi)$ を図示する。 引数は主枝の傾斜角の変動幅 $K$（$\in Ks = (5, 10, 15, 20, 30, 45, 60, 90)$）で、枝分かれ角、投影像の角の属する区間（最小値と幅で指定）をスライダーで指定し、それらに対する傾斜角の組を図示する。

In [16]:

plotArgs(90)

HBox(children=(Output(), VBox(children=(IntSlider(value=90, description='枝分かれ角 :', layout=Layout(width='182pt'…

度数¶

投影像の角の計算データ arg_of_edges を整形したデータ aes をもとに、枝分かれ角（$\theta$）、主枝の傾斜角（$\phi$）の変動幅（$K$ : $|\phi|\leqq K$）を指定し、投影像の角がある区間に属するときの度数の計算
　　$K=90^\circ$, $\theta=90^\circ$ のときの、投影像の角が $20^\circ\sim21^\circ$ になる $(\psi, \phi)$ の度数
　　上のグラフの初期値、「枝分かれ角」$90$、「投影像の角」$20$、「幅」$1$ のときの「度数」の値と同じ

In [17]:

freq(K=90, theta=90, interval=(20,21))

Out[17]:

6009

$K=90^\circ$, $\theta=40^\circ$ のときの、投影像の角が $30^\circ\sim33^\circ$ になる $(\psi, \phi)$ の度数
　　上のグラフで、「枝分かれ角」$40$、「投影像の角」$30$、「幅」$3$ のときの「度数」の値と同じ

In [18]:

freq(K=90, theta=40, interval=(30,33))

Out[18]:

68076

In [ ]:

4.3. 投影像の角の度数分布グラフ¶

投影像の角の度数分布は、枝分かれ角（$\theta$）と、主枝の傾斜角（$\phi$）の変動幅 $K$（$|\phi|\leqq K$）に対して定まります。比較のため、分布を重ねた２２種のグラフ :
　　(1) 枝分かれ角を固定し、変動幅を変えた投影角の度数分布を重ねたグラフ
　　(2) 変動幅を固定し、枝分かれ角を変えた投影像の角の度数分布を重ねたグラフ
を描く。また、実験において写真上の角度計測の際に、基準となる主枝の伸びる向きが不明確になってしまう場合も想定し、モデルにおける投影像の角として向きを考慮するかしないかで２種 :
　　(A) 投影像上の主枝の伸びる向きに対する分枝の伸びる向きの角度（$0^\circ\sim 180^\circ$）
　　(B) 主枝の伸びる向きを考慮せず、単純に投影像上の主枝と分枝の角度（$0^\circ\sim 90^\circ$）
を考え、それらを組み合わせた４種のグラフ (1-A), (1-B), (2-A), (2-B) を描く。

(1-A) : 枝分かれ角固定、向きあり¶

枝分かれ角（$\theta$）を固定し、主枝の傾斜角（$\phi$）の変動幅 $K = 5,10,15,20,30,45,60,90$ に対する投影像の角の度数分布グラフを重ねて描きます。 横軸は投影像の角で、主枝の伸びる向きに対する角度とし、$\theta-5^\circ\sim \theta+5^\circ$ の範囲を描画します。 縦軸は度数分布（％）です。 度数分布データに freq_dist を使っています。 plot.saving の引数に適当な画像ファイル名を指定することで、グラフを画像で保存できます。

In [19]:

for theta in thetas:
    print('θ = {}'.format(theta))
    for K in Ks:
        plt.plot(*listTranspose(freq_dist[K][theta]))
    plt.xlim([theta-5, theta+5]);
    #plt.ylim([0,40])]);
    plt.xlabel('observed angles');
    plt.ylabel('frequency distribution (%)');
    plt.show()
    #plt.savefig('figure-1A-{}.jpg'.format(theta))

θ = 5

θ = 10

θ = 15

θ = 20

θ = 25

θ = 30

θ = 35

θ = 40

θ = 45

θ = 50

θ = 55

θ = 60

θ = 65

θ = 70

θ = 75

θ = 80

θ = 85

θ = 90

(1-B) : 枝分かれ角固定、向きなし¶

枝分かれ角（$\theta$）を固定し、主枝の傾斜角（$\phi$）の変動幅 $K = 5,10,15,20,30,45,60,90$ に対する投影像の角の度数分布グラフを重ねて描きます。 横軸は投影像の角で、主枝の伸びる方向は考慮しない主枝との角度とし、描画範囲は $\theta-5^\circ\sim \theta+5^\circ$ です。 縦軸は度数分布（％）です。 度数分布データに freq_dist90 を使っています。 plot.saving の引数に適当な画像ファイル名を指定することで、グラフを画像で保存できます。

In [20]:

for theta in thetas:
    print('θ = {}'.format(theta))
    for K in Ks:
        plt.plot(*listTranspose(freq_dist90[K][theta]))
    plt.xlim([theta-5,theta+5]);
    #plt.ylim([0,40]);
    plt.xlabel('observed angles');
    plt.ylabel('frequency distribution (%)');
    plt.show()
    #plt.savefig('figure-1B-{}.jpg'.format(theta))

θ = 5

θ = 10

θ = 15

θ = 20

θ = 25

θ = 30

θ = 35

θ = 40

θ = 45

θ = 50

θ = 55

θ = 60

θ = 65

θ = 70

θ = 75

θ = 80

θ = 85

θ = 90

(2-A) : 変動幅固定、向きあり¶

主枝の傾斜角（$\phi$）の変動幅 $K = 5,10,15,20,30,45,60,90$ を固定し、枝分かれ角に対する投影像の角の度数分布グラフを重ねて描きます。 横軸は投影像の角、縦軸は度数分布（％）です。 投影像の角は主枝の伸びる向きに対する角度とし、$0^\circ$ 以上 $180^\circ$ 以下の値をとる。 度数分布データに freq_dist を使っています。 plot.saving の引数に適当な画像ファイル名を指定することで、グラフを画像で保存できます。

In [21]:

for K in Ks:
    for theta in thetas:
        plt.plot(*listTranspose(freq_dist[K][theta]))
    print('K = {}'.format(K))
    #plt.ylim([0,40])]);
    plt.xlabel('observed angles');
    plt.ylabel('frequency distribution (%)');
    plt.show()
    #plt.savefig('figure-2A-{}.jpg'.format(K))

K = 5

K = 10

K = 15

K = 20

K = 30

K = 45

K = 60

K = 90

(2-B) : 変動幅固定、向きなし¶

主枝の傾斜角（$\phi$）の変動幅 $K = 5,10,15,20,30,45,60,90$ を固定し、枝分かれ角（$\theta$）に対する投影像の角の度数分布グラフを重ねて描きます。 横軸は投影像の角、縦軸は度数分布（％）です。 投影像の角は、主枝の伸びる方向は考慮せず、単に主枝との角度とし、$0^\circ$ 以上 $90^\circ$ 以下の値をとる。 度数分布データに freq_dist90 を使っています。

In [22]:

for K in Ks:
    for theta in thetas:
        plt.plot(*listTranspose(freq_dist90[K][theta]))
    print('K = {}'.format(K))
    #plt.ylim([0,40])]);
    plt.xlabel('observed angles');
    plt.ylabel('frequency distribution (%)');
    plt.show()
    #plt.savefig('figure-2B-{}.jpg'.format(K))

K = 5

K = 10

K = 15

K = 20

K = 30

K = 45

K = 60

K = 90

In [ ]: