20191007

断面二次モーメントを、座標点の配列から計算する(2)断面二次モーメント編

断面二次モーメント

やろうとしている事

断面二次モーメントを数値計算的に算出しようとしている。(先々月の続き。)

計算の仕方は二つ

  1. 長さ寸法と、公式に当てはめて計算をする。

     公式の公式集。書籍よりDVD版の方が明らかにお得だが高い。

    • メリット
      • 正確な値
    • ディメリット
      • 形状によっては式が用意されていない。
  2. 座標値を使い、数値計算的に求める

  3. メリット
    • 断面形状にかかわらず、座標値さえあれば計算可能。
  4. ディメリット
    • 点の間隔、カーブ、数値積分の手法などで誤差が生じる。
    • 座標値を準備しなければならない。

1の方法が普通のやり方。できるなら1をやるべき。

断面二次モーメントの公式

$$\int _{\rm A}y^2\rm\ dA$$

上は、断面二次モーメントの式。

下は、8×8の正方形の断面二次モーメントを計算する際の積分のイメージ。

20190918130448_2.png

数値計算で断面二次モーメントを求めるための、計算手順

断面二次モーメントまでの道のり
  1. 座標値を手に入れる

  2.  - 断面一次モーメントを計算する。

    • 断面一次モーメントから図心を計算する。
    • 図心が原点となるように座標値を整える。(省略)
  3. 断面二次モーメントを計算する。

1. 座標値を手に入れる 四角形のx,y座標を得る関数
def coorsquare(r=4, center=[0,0], n=10, closed=True):
rx1=-r/2.0+center[0]
rx2=r/2.0+center[0]
ry1=-r/2.0+center[1]
ry2=r/2.0+center[1]

x = [rx1 + i*(rx2-rx1)/float(n) for i in range(n)]
y = [ry2 for i in range(n)]
#1
x = x + [rx2 for i in range(n)]
y = y + [ry2 - i*(ry2-ry1)/float(n) for i in range(n)]
#2
x = x + [rx2 - i*(rx2-rx1)/float(n) for i in range(n)]
y = y + [ry1 for i in range(n)]
#3
x = x + [rx1 for i in range(n)]
y = y + [ry1 + i*(ry2-ry1)/float(n) for i in range(n)]
#4
if closed:
x = x + [rx1]
y = y + [ry2]

return x, y


x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
上の関数を使って作ったx,yを作図。
import matplotlib.pyplot as plt

fig,ax1 = plt.subplots(1)
ax1.plot(x,y)
ax1.scatter(x,y)

for i in range(len(x)):
ax1.text(x[i],y[i],"%s"%(i))

ax1.set_xlim(-4,4)
ax1.set_ylim(-4,4)
ax1.set_aspect("equal")

plt.show()

shikaku.png

2. 面積を計算→断面一次モーメント→図心を計算。
def area(x, y):
A = sum([(y[i]+y[i+1])/2.0*(x[i+1]-x[i]) for i in range(len(x)-1)])
return A

x,y = coorsquare(r=5, center=[0,0], n=20, closed=True)
print(area(x,y))

def FirstInrt(x, y):
FIy = sum([(y[i]+y[i+1])**2/8.0*(x[i+1]-x[i]) for i in range(len(x)-1)])
FIx = sum([(x[i]+x[i+1])**2/8.0*(y[i]-y[i+1]) for i in range(len(x)-1)])
return FIx, FIy

print(FirstInrt(x, y))


def CntCS(x, y):
A = area(x, y)
FIx, FIy = FirstInrt(x, y)
Cntx, Cnty = FIx/A, FIy/A
return Cntx, Cnty

print("center(%s,%s)"%(CntCS(x, y)))
3. 断面二次モーメントの計算(シンプソン則)

二次の式なので、台形則でなくシンプソン則で計算した。この関数が、この文章の一番重要な所

def simpson(x,y):#sekibun_1/3simpson(3point)
Cntx, Cnty = CntCS(x, y)
x = np.array(x)-Cntx
y = np.array(y)-Cnty
SecondInrty,SecondInrtx = 0,0

x1 = np.array(x)**3/3.0
y1 = np.array(y)**3/3.0

for i in range(len(x))[::2][:-1]:
a1 = float(x[i+1]-x[i])
# if は、0割対策。
if a1 !=0.0:
a1 = y1[i] * (x[i+2]-x[i])*(-x[i+2]-2.0*x[i]+3.0*x[i+1])/a1/6.0
a2 = float((x[i+1]-x[i])*(x[i+2]-x[i+1]))
if a2 !=0.0:
a2 = y1[i+1] * (x[i+2]-x[i])**3/a2/6.0
a3 = float(x[i+2]-x[i+1])
if a3 !=0.0:
a3 = y1[i+2] * (x[i+2]-x[i])*(2.0*x[i+2]+x[i]-3.0*x[i+1])/a3/6.0

SecondInrty = SecondInrty + a1 + a2 + a3

b1 = y[i+1]-y[i]
if b1 !=0:
b1 = x1[i]*(y[i+2]-y[i])*(-y[i+2]-2*y[i]+3*y[i+1])/b1/6.0
b2 = (y[i+1]-y[i])*(y[i+2]-y[i+1])
if b2 !=0:
b2 = x1[i+1]*(y[i+2]-y[i])**3/6.0/b2
b3 = y[i+2]-y[i+1]
if b3 !=0:
b3 = x1[i+2]*(y[i+2]-y[i])*(2*y[i+2]+y[i]-3*y[i+1])/b3/6.0

SecondInrtx = SecondInrtx - (b1 + b2 + b3)

return SecondInrtx, SecondInrty


# 断面が中心にない場合でも、図心を原点に持って行った時の値になるようにした。
x,y = coorsquare(r=5, center=[7,7], n=8, closed=True)
print("1/3simpson(3point)%s,%s"(simpson(x,y)))

# nは、偶数にしないと、正しく計算されない。
#うまくいかない。
x,y = coorsquare(r=5, center=[0,0], n=7, closed=True)
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

# 全体で偶数になっていても、一辺が奇数ではだめ。
#うまくいかない
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
x.insert(-9,(-2.5-1.875)/2.0)
y.insert(-9,-2.5)
y.insert(-1,(2.5+1.875)/2.0)
x.insert(-1,-2.5)
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

# nが偶数なら、不等間隔でも計算される。
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
x.insert(-9,(-2.5-1.875)/2.0)
y.insert(-9,-2.5)
x.insert(-9,-2.4)
y.insert(-9,-2.5)
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

for ループは、2個飛ばしで行う。

3点使って三個飛ばしだと、間が空く。

3点使って二個飛ばしだと、間が埋まる。

ごめん、うまく書けない。こういう事。

figure_1.png

ってなっちゃうから、3個飛ばしにすればいいわけじゃない。

3個使うけど、最後の一個は次のループでも使うから、二個飛ばしが正解。
(ちなみに、pythonで二個飛ばしの書き方は、[::2]。)
こういうindexingの所で、いつも試行錯誤してる。頭だけで出来ないとだめだろうか?

figure_1-2.png

前のループでi+2だったところを、次のループではiとして扱っている。

積分は不等間隔のシンプソン則で模擬。不等間隔のシンプソン則は、ココを参考にして作った。

hutou_kankaku.png

一辺のサンプルが奇数ならば正しく計算できる。

ひし形(45度回転。一部不等間隔。)
x,y = coorsquare(r=5, center=[0,0], n=8, closed=True)
x.insert(-9,(-2.5-1.875)/2.0)
y.insert(-9,-2.5)x.insert(-9,-2.4)
y.insert(-9,-2.5)
xy= np.dot(np.array([x,y]).T, [[np.cos(1/4.0*np.pi), np.sin(1/4.0*np.pi)],[-np.sin(1/4.0*np.pi), np.cos(1/4.0*np.pi)]]).T
x=xy[0]
y=xy[1]

print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

hishi_gata.png

前回は、これが正しく計算できなかったが、今回は不等間隔シンプソン則のおかげで計算できた。
(サンプル数の偶奇に気を付けなければならないが)

穴あき(正方形)
x1,y1 = coorsquare(r=5, center=[0,0], n=8, closed=True)
x2,y2 = coorsquare(r=3, center=[0,0], n=4, closed=True)
x2=x2[::-1]
y2=y2[::-1]
Ix1,Iy1 = simpson(x1,y1)
Ix2,Iy2 = simpson(x2,y2)
a = Ix1+Ix2
b = Iy1+Iy2
print("1/3simpson(3point)%s,%s"%(a,b))
import matplotlib.pyplot as plt

fig,ax1 = plt.subplots(1)
ax1.plot(x1,y1)
ax1.plot(x2,y2)
ax1.scatter(x1,y1)
ax1.scatter(x2,y2)
for i in range(len(x1)):
ax1.text(x1[i],y1[i],"%s"%(i))
if i<len(x2):
ax1.text(x2[i],y2[i],"%s"%(i))

ax1.set_xlim(-3,3)
ax1.set_ylim(-3,3)
ax1.set_aspect("equal")

plt.show()

ana_aki.png

緑色の部分が反時計回りになっている所が重要。反時計回りはdx(差分)がマイナスになるので、負の面積というか負の断面二次モーメントが出てくる。

円形
N=21
x = np.cos(np.linspace(0, 2.0*np.pi, N))
y = np.sin(np.linspace(0, 2.0*np.pi, N))
print("1/3simpson(3point)%s,%s"%(simpson(x,y)))

maru.png

円は、理論値からの誤差がまだある。どのくらいサンプル間隔を細かくすればいいのか考えるために、下で誤差などを確認した。

円の理論値

$$\frac{d^{4}\pi}{64} \\ 半径r=1なので直径d=2\\ \frac{16\pi}{64}=0.25\pi \\ =‭0.78539816339744830961566084581988 \\ (電卓でπ=3.1415926535897932384626433832795‬を使い計算)‬$$

円周率にはnp.pi 3.141592653589793(16桁)を使うのでそれが原因の誤差が出るが、確認しない。

va = 0.25 * np.pi
Ix = []
Iy = []
Acx = []
Acy = []
Erx = []
Ery = []

for N in np.array(np.logspace(1,4,10),"int"):
#for N in np.array(np.logspace(1,4,10),"int")*4+1:
x = np.cos(np.linspace(0, 2.0*np.pi, N))
y = -np.sin(np.linspace(0, 2.0*np.pi, N))#Clockwise

I = simpson(x,y)
Ix.append(I[0])
Iy.append(I[1])

Acx.append(Ix[-1]/va)
Acy.append(Iy[-1]/va)

Erx.append(1 - Acx[-1])
Ery.append(1 - Acy[-1])
print("%s-points,\tIx:%s,\tIy:%s,\tAcx:%s,\tAcy:%s,\tErx:%s,\tEry:%s,\t"%(N,Ix[-1],Iy[-1],Acx[-1],Acy[-1],Erx[-1],Ery[-1]))
import matplotlib.pyplot as plt

x = np.array(np.logspace(1,4,10),"int")#*4+1
fig,(ax1,ax2) = plt.subplots(2,1)
ax1.plot(x,[va for i in range(len(x))])
ax1.scatter(x, Acy)

ax2.plot(x,[1 for i in range(len(x))])
ax2.scatter(x, Ery)

ax1.set_ylim(0.99,1)
ax2.set_ylim(0.00000000000001, 0.1)

ax1.set_xscale("log")
ax2.set_xscale("log")
ax1.set_yscale("log")
ax2.set_yscale("log")

plt.show()

figure_1 (3).png

上が精度で1に近づけばいい値。理論値で割って計算した。y軸はなんか知らんがすっ飛んだ。0.99~1.0間。

下が誤差で、0に近づけばいい値。1-精度で計算。両方とも対数グラフ。

サンプル数によっては、うまく計算できないのがある。実行して、printの結果を確認すればわかるけど、10の時に断面二次モーメントに異常な値が出る。dxの±が入れ替わる所か、yの±が入れ替わる所で起こっているのではないかと思っているが、確認は次回以降にする。とにかく、サンプル数に気を遣わずにシンプソン則を実行できなければならない。あと、欲を言えば、excelで出来れば言うことない。有料ソフトなのに、どこに行っても必ずあるし。このままでは、ベジェ曲線で自由に描いた線の積分ができない。円形で起こった現象は、ひし形でも出てきそう。正方形で出てこなかったのは、dxがぴったし0だから。

実用上は、40個位の分割でも円の断面二次モーメントを10-4の誤差で計算できる事が知れてよかった。

ちなみに、データ点数を4n+1(nは整数)にすれば、異常な値は出てこない。

現在の課題

  • ひし形、丸形などで、誤差が大きくなる。
  • 任意形状(ひらがなの「あ」とか)でやってない。

次回以降の目標

- 【断面二次】 シンプソン則のサンプル数を気にせず実行。

  • 【断面二次】Excelでやる。

  • 一覧表を作成し、誤差等をチェック
  • 【three.js】 pdbを表示する

絶対やるというわけではない。

posted by yuchan at 08:00 | Comment(0) | python

20190831

断面二次モーメントを、座標点の配列から計算する(1)計算手順編

やること

断面二次モーメントを求める

断面二次モーメントを計算する方法は二つある。

断面形状の…

  1. 長さ寸法と理論式を使って計算する方法
  2. 座標値を使って、数値計算的に計算する方法。

がある。ここでは2をやる。

1と2の違いは、

  • 積分を形状別に解いて消してから、計算する場合

    と、

  • 積分を、狽ナ模擬してそのまま計算する場合

の違い。

ここでは、後者をやる。

任意形状の断面二次モーメントを求める

このページでは、閉じた線で描かれた任意形状の、x、y座標を使って、(面積や)断面二次モーメントを求める手順を考えてみた。

  1. 断面形状の縁を、閉じた線で描く様な、x,y座標を得る。
  2. 1を使って面積を求める。
  3. zに関する、もしくはyに関する断面一次モーメントを求める。
  4. 図心を求める。
    図心(y,z)=(z軸に関する断面一次モーメント, y軸に関する断面一次モーメント)/面積
  5. 1の座標から4の図心を引いた値を作り、その図形の断面二次モーメントを計算する。

普通は、断面二次モーメントの計算には理論式を使う。

断面二次モーメントは、普通、

  • 長さ寸法と理論式を使って計算する。
    これを実行するには…
    • 式を知っている必要がある。
      (→式が書いてある本等が要る)
      (下記amazon参照)
    • 断面形状によって式が違う。
    • 任意形状では式がない。
  • 理論式がない形状の場合、(solidworksなどの)有名なソフトの断面性能表示機能を使う。
  • 有償ソフトが使用できない場合、エクセルなどで自作する。(このページではexcelでやらない。pythonでやる。excelだって有料。)

のどれかで計算すると思う。
2つめ、3つめは、コンピュータで計算するので、理論値からの誤差がある。

←電子化(pdf(メニュー検索機能付き))全部入りDVD←

→1冊/全29冊 紙媒体 目的の式がある。→


↑断面二次モーメントの式が入っている。

長さ寸法を使うなら上の本などを使って計算されたし。


前述したが、閉じた線で描かれた任意形状の、x、y座標を使って、(面積や)断面二次モーメントを求める

1~5をやってみる

1. 文字データのx,y座標の取得

使用する座標値一覧は下の表。閉じた線を表す座標で、図形を表現している。これらの図形は、図心を中心にして作っていない。

図形 座標値
四角x[0,1,2,3,4,5,6,7,8,9,10,10,10,10,10,10,10,10,10,10,10,9,8,7,6,5,4,3,2,1,0,0,0,0,0,0,0,0,0,0,0]
 y[10,10,10,10,10,10,10,10,10,10,10,9,8,7,6,5,4,3,2,1,0,0,0,0,0,0,0,0,0,0,0,1,2,3,4,5,6,7,8,9,10]
せん断変換した四角x[41.0, 42.0, 43.0, 44.0, 45.0, 46.0, 47.0, 48.0, 49.0, 50.0, 51.0, 49.0, 47.0, 45.0, 43.0, 41.0, 39.0, 37.0, 35.0, 33.0, 31.0, 30.0, 29.0, 28.0, 27.0, 26.0, 25.0, 24.0, 23.0, 22.0, 21.0, 23.0, 25.0, 27.0, 29.0, 31.0, 33.0, 35.0, 37.0, 39.0, 41.0]
 y[17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 16.0, 15.0, 14.0, 13.0, 12.0, 11.0, 10.0, 9.0, 8.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 17.0]
45度回転した四角x[16.97056274847714, 17.67766952966369, 18.384776310850235, 19.09188309203678, 19.798989873223327, 20.506096654409877, 21.213203435596427, 21.920310216782973, 22.62741699796952, 23.334523779156065, 24.041630560342615, 23.334523779156065, 22.62741699796952, 21.920310216782973, 21.213203435596427, 20.506096654409877, 19.798989873223327, 19.09188309203678, 18.384776310850235, 17.67766952966369, 16.97056274847714, 16.26345596729059, 15.556349186104043, 14.849242404917495, 14.14213562373095, 13.435028842544401, 12.727922061357853, 12.020815280171306, 11.31370849898476, 10.606601717798211, 9.899494936611664, 10.606601717798211, 11.31370849898476, 12.020815280171306, 12.727922061357853, 13.435028842544401, 14.14213562373095, 14.849242404917495, 15.556349186104043, 16.26345596729059, 16.97056274847714]
 y[7.0710678118654755, 6.363961030678928, 5.65685424949238, 4.949747468305833, 4.242640687119286, 3.5355339059327378, 2.82842712474619, 2.1213203435596437, 1.4142135623730958, 0.7071067811865479, 0.0, -0.7071067811865479, -1.4142135623730958, -2.1213203435596437, -2.82842712474619, -3.5355339059327378, -4.242640687119286, -4.949747468305833, -5.65685424949238, -6.363961030678928, -7.0710678118654755, -6.363961030678928, -5.65685424949238, -4.949747468305832, -4.242640687119286, -3.5355339059327378, -2.82842712474619, -2.121320343559643, -1.4142135623730958, -0.7071067811865479, 0.0, 0.7071067811865479, 1.4142135623730958, 2.121320343559643, 2.82842712474619, 3.5355339059327378, 4.242640687119286, 4.949747468305832, 5.65685424949238, 6.363961030678928, 7.0710678118654755]
x[1.0, 0.94581724, 0.78914051, 0.54694816, 0.24548549, -0.08257935, -0.40169542, -0.67728157, -0.87947375, -0.9863613, -0.9863613, -0.87947375, -0.67728157, -0.40169542, -0.08257935, 0.24548549, 0.54694816, 0.78914051, 0.94581724, 1.0]
 y[0.0, 3.24699469e-01, 6.14212713e-01, 8.37166478e-01, 9.69400266e-01, 9.96584493e-01, 9.15773327e-01, 7.35723911e-01, 4.75947393e-01, 1.64594590e-01, -1.64594590e-01, -4.75947393e-01, -7.35723911e-01, -9.15773327e-01, -9.96584493e-01, -9.69400266e-01, -8.37166478e-01, -6.14212713e-01, -3.24699469e-01, 0.0]
穴あき四角x[[0,1,2,3,4,5,6,7,8,9,10,10,10,10,10,10,10,10,10,10,10,9,8,7,6,5,4,3,2,1,0,0,0,0,0,0,0,0,0,0,0],[1,1,1,1,1,1,1,1,1,2,3,4,5,6,7,8,9,9,9,9,9,9,9,9,9,8,7,6,5,4,3,2,1]]
 y[[10,10,10,10,10,10,10,10,10,10,10,9,8,7,6,5,4,3,2,1,0,0,0,0,0,0,0,0,0,0,0,1,2,3,4,5,6,7,8,9,10],[9,8,7,6,5,4,3,2,1,1,1,1,1,1,1,1,1,2,3,4,5,6,7,8,9,9,9,9,9,9,9,9,9]]
ひらがなの「あ」x[[1190, 1121.0, 1153, 1319.0, 1581.0, 1825.0, 1567.0, 1061.0, 1045, 1692.0, 1547.0, 1323.0, 1180.0, 910.0, 920.0, 930.0, 942.0, 899.0, 860.0, 797.0, 410.0, 244.0, 510.0, 797.0, 806.0, 811.0, 815, 717.0, 404.0, 443, 672.0, 824.0, 828, 830, 838.0, 738.0, 766, 981.0, 961.0, 930.0, 1182.0, 1186.0, 1243.0, 1321.0, 1405.0, 1434.0, 1393.0, 916.0, 893.0, 1190.0], [1194, 889.0, 887.0, 895.0, 1137.0, 1194.0], [793, 789.0, 790.0, 793.0, 621.0, 371.0, 459.0, 793.0]]
 y[[1008, 1094.0, 1130, 1006.0, 926.0, 541.0, 86.0, -74.0, -27, 541.0, 846.0, 938.0, 639.0, 332.0, 286.0, 250.0, 174.0, 57.0, 35.0, 240.0, 63.0, 268.0, 749.0, 913.0, 1054.0, 1124.0, 1176, 1171.0, 1333.0, 1366, 1251.0, 1262.0, 1307, 1343, 1462.0, 1618.0, 1661, 1522.0, 1421.0, 1280.0, 1356.0, 1358.0, 1393.0, 1419.0, 1399.0, 1366.0, 1321.0, 1186.0, 948.0, 1008.0], [942, 874.0, 680.0, 432.0, 786.0, 942.0], [334, 561.0, 663.0, 831.0, 719.0, 272.0, 163.0, 334.0]]
数字の「8」x[[908, 613.0, 740.0, 1043.0, 1322.0, 1459.0, 1143.0, 1528.0, 1401.0, 1028.0, 695.0, 535.0, 908.0], [1065, 1082.0, 1088, 1330.0, 1266.0, 1043.0, 877.0, 752.0, 850.0, 977.0, 1065.0], [977, 957.0, 668.0, 768.0, 1026.0, 1260.0, 1366.0, 1255.0, 1088.0, 977.0]]
 y[[856, 1225.0, 1477.0, 1573.0, 1495.0, 1231.0, 893.0, 457.0, 187.0, 74.0, 162.0, 453.0, 856.0], [924, 930.0, 934, 1231.0, 1413.0, 1507.0, 1460.0, 1229.0, 1043.0, 963.0, 924.0], [824, 814.0, 455.0, 234.0, 140.0, 222.0, 449.0, 664.0, 777.0, 824.0]]

上の表の座標値の作り方に関しては3種類あり…

  • 四角は、試行錯誤しながら手打ち。穴あきも苦労したが手打ち。四角だったから何とか打てた。蕎麦みたいな書きっぷり。
  • せん断、45度回転は、アフィン変換して作った。
  • 丸は、x座標はcosθ、y座標はsinθ (0<=θ<=2πを20分割(N=20))で作った。cosとかsinは、pythonのmathモジュールを使った。
  • ひらがなの「あ」、数字の「8」は、先月のページのやり方で作った。
    このページ下部に、フォントの文字データのx,y座標が取得のためのコードを載せました。
    (先月のと同じ内容です。IPAフォントのインストールを忘れずに。)

    フォントファイルのグリフのベクターを座標データに変換しています。

#!/usr/bin/env python
# -*- coding: cp932 -*-

x=ここに上表の右列のxを入力
y=ここに上表の右列のyを入力

こんな感じで変数x,yに入れておき、後のコードで使う。

変数x、yは、曲線のx座標y座標が入った配列です。
また、変数x2、y2の中には、一文字を表すために必要なパーツの分だけ、配列が入っています。
つまり、漢数字の「一」は一筆書できるので1つしか配列が入っていませんが、
同じ漢数字の「二」には配列が二つ(長さ2の配列)入っています。

上の座標を示して描ける線の特徴として、

  • 閉じた線
  • 外側が時計回り、内側が反時計回り

などがあります。

閉じた線

線を指定する座標は、最初の一点目の座標と同じ座標で終わるように作りました。線の途中で交差したりしていません。線の途中で交差すると、面積の計算に不都合があります。

時計回り、反時計回り。

閉じた線のx、y座標をあらわす配列が、時計回りに動いているか、それとも反時計回りに動いているかは、少なくとも下の二つの用途で重要です。

  • 文字の描画
    「の」や「○(まる)」などの穴開き文字(正式名称はあるのかすら知らない)に関しては、穴が開いている部分を表す配列が反時計回りになっている。下の図でいう緑の線。

の.png〇.png

  • 面積の計算
    ∫ydxすると、xが+に動くときdxが+に、xが-に動くときdxが-になる。これがうまい具合に効いて、閉じた曲線の面積が求まる。
    また、時計回り、反時計回りが切り替わると、負の面積が出ることがある。負の面積によって、文字の穴空き部分の面積が引き算される仕組みになっている。
    閉じた線の中で、交差があるとコレがうまくいかない。

断面形状の縁x、y座標をプロットする

2画以上の文字や、穴空き文字を表現するため、変数x,yは長さが1の配列ではない場合があります。

だから、matplotlib.pyplot.plotでそのまま使えません。matplotlib.pyplot.plotは、x,yの引数は、N*1のリストとかだからです。(Excelで折れ線グラフ一本書くときは、x,yそれぞれ一列ずつ選ぶのと一緒)

一つの変数の中に、たくさん作図すべき閉じた線が含まれる事があるので、

繰り返し構文+indexingで、plotしました。
上の表の座標を

x=

y=

とか書いた後に、下記コードの実行でx、yをプロットできます。

  • 必要なモジュール
    • matplotlib
    • datetime
    • time
#!/usr/bin/env python
# -*- coding: cp932 -*-

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'IPAGothic'

import time
from datetime import datetime as dt

f,(ax1) = plt.subplots(1,1)

if len(np.shape(y))==1:
ax1.plot(x,y)
ax1.scatter(x,y,facecolor=[0,0,0,0.3],edgecolor=[0,0,0])
else:
for i in range(len(x)):
ax1.plot(x[i],y[i])
ax1.scatter(x[i],y[i],facecolor=[0,0,0,0.3],edgecolor=[0,0,0])

ax1.set_aspect("equal")
#plt.show()

tdatetime = dt.now()
timtext = tdatetime.strftime('%Y%m%d'+ time.ctime()[11:13] + time.ctime()[14:16] + time.ctime()[17:19])
fname = "%s.png"%(timtext)
plt.savefig(fname)

plt.close()

2.1を使って面積を求める。

x、yでプロットされた図形の面積は台形則で求める。

\int\int f(x,y) \ dA

面積を計算するだけなので、関数とか要らないけど、
f(x,y)=1

\int\int 1\ dydx

これを、台形則で模擬します。

#!/usr/bin/env python
# -*- coding: cp932 -*-

import numpy as np

A=0if type(x[0])!=list:
dx = [x[i+1]-x[i] for i in range(len(x)-1)]
trapzy = [(y[i]+y[i+1])/2.0 for i in range(len(y)-1)]
A = sum([trapzy[i]*dx[i] for i in range(len(y)-1)])
else:
for j in range(len(x)):#1文字に、複数の閉じた線が存在する場合対策
dx = [x[j][i+1]-x[j][i] for i in range(len(x[j])-1)]
trapzy = [(y[j][i]+y[j][i+1])/2.0 for i in range(len(y[j])-1)]
A1 = sum([trapzy[i]*dx[i] for i in range(len(y[j])-1)])
A = A+A1

print A

曲線とは一致しないとは思いますが、曲線を間引いたデータとしては、正しい値が出てくると思います。もちろん細かく座標点をとればとるほど、正確な値に近づきます。
PCで扱うデータは、しょせん間が抜けた折れ線データに過ぎないので、面積に関しては 台形則で表してもいいと思います。

 四角せん断変形した四角45度回転させた四角穴の開いた四角ひらがなの「あ」数字の「8」
図形square.pngsquare_sheared.pngsquare_45deg_rotated.png circle.pngsquare_hollow.pnggryph_a.pnggryph_8.png
面積100.0100.0100.0-3.08464495357…36694613.0402222.0

ついでに せん断変形で、面積が変わっていない事を確認

せん断写像のwikipediaにそう書いてあったきがする ので、せん断前後で面積に変化がない事を確認する。
せん断変換前後の図形のx,y座標をせん断変換(アフィン変換)で用意したので、
座標から面積を計算した結果を比較すると、せん断変形しようが、面積は変わらないことが確認できました。


もう理論値からズレてる。

曲線をPCで表現するとき、当然細かい折れ線で描く事になります。
間引いた折れ線データを取得した時点で、データに誤差が生じています。

だから、この後の操作がどんなに正しくとも、座標点として記録した時点で理論値からのズレはあります。
「それでもやるとしたら、どうすべきなの?」という観点で見て下さい。

半径1の円で、線の間引き方によって、値がずれていることが確認できるので載せました。

 粗い円やや粗い円細かい円理論値(半径2×3.14…)
図形circle_coarth2.pngcircle.pngcircle_precise.pngcircle_ideal.png
面積-2.89254424359…-3.08464495357…-3.14157194138…π≒3.14159265358…

3. zに関する、もしくはyに関する断面一次モーメントを求める。

y軸に関する断面一次モーメント $$\int y\ dA$$

x軸に関する断面一次モーメント$$\int x\ dA$$

断面一次モーメントは、図形がどの位置に描いてあるかによって、値が変わります。
原点から離れるほど大きくなります。

断面二次モーメントを計算する時も、図形が原点から離れるほど大きくなります。断面二次モーメントは、図心を原点と一致させてから求めることになります。次の手順4で、図形の重心(図心)を求めます。

断面一次モーメントは、その図心を求めるために計算します。

断面一次モーメントを求めるコードは下のように書いてみました。

#!/usr/bin/env python
# -*- coding: cp932 -*-
import numpy as np

C1=0
if type(x[0])!=list:
dx = [x[i+1]-x[i] for i in range(len(x)-1)]
trapzy = [((y[i]+y[i+1])/2.0)*((x[i]+x[i+1])/2.0) for i in range(len(y)-1)]
C = sum([trapzy[i]*dx[i] for i in range(len(y)-1)])
C1 = C1+C
else:
for j in range(len(x)):#1文字に、複数のベクトルが存在する場合対策
dx = [x[j][i+1]-x[j][i] for i in range(len(x[j])-1)]
trapzy = [((y[j][i]+y[j][i+1])/2.0)*((x[j][i]+x[j][i+1])/2.0) for i in range(len(y[j])-1)]
C2 = sum([trapzy[i]*dx[i] for i in range(len(y[j])-1)])
C1 = C1+C2

print("X:%s"%(C1))

C3=0

if type(x[0])!=list:
dy = [y[i+1]-y[i] for i in range(len(x)-1)]
trapzx = [((y[i]+y[i+1])/2.0)*((x[i]+x[i+1])/2.0) for i in range(len(y)-1)]
C4 = sum([trapzx[i]*dy[i] for i in range(len(y)-1)])*-1
C3 = C3+C4

else:
for j in range(len(x)):#1文字に、複数のベクトルが存在する場合対策
dy = [y[j][i+1]-y[j][i] for i in range(len(x[j])-1)]
trapzx = [((y[j][i]+y[j][i+1])/2.0)*((x[j][i]+x[j][i+1])/2.0) for i in range(len(y[j])-1)]
C4 = sum([trapzx[i]*dy[i] for i in range(len(y[j])-1)])*-1
C3 = C3+C4

print("Y:%s"%(C3))
 四角せん断変形した四角45度回転させた四角穴の開いた四角ひらがなの「あ」数字の「8」
図形square.pngsquare_sheared.pngsquare_45deg_rotated.pngcircle.pngsquare_hollow.pnggryph_a.pnggryph_8.png
断面一次モーメントX:500.0
Y:500.0
X:3600.0
Y:1200.0
X:1697.05627485
Y:1.42108547152e-14≒0
X:1.90819582357e-17≒0
Y:1.11022302463e-16≒0
X:180.0
Y:180.0
X:713082273.0
Y:475455483.5
X:419443127.5
Y:326005988.0

4. 図心を求める。

図心(x,y)=(y軸に関する断面一次モーメント, x軸に関する断面一次モーメント)/面積

あえて書く必要はない気がするけど、

  1. 先ず、面積を計算する。

  2. 断面一次モーメントを計算する

  3. 断面一次/面積

    してるのが下記。

#!/usr/bin/env python
# -*- coding: cp932 -*-

import
numpy as np
A=0

if type(x[0])!=list:
dx = [x[i+1]-x[i] for i in range(len(x)-1)]
trapzy = [(y[i]+y[i+1])/2.0 for i in range(len(y)-1)]
A = sum([trapzy[i]*dx[i] for i in range(len(y)-1)])
else:
for j in range(len(x)):#1文字に、複数の閉じた線が存在する場合対策
dx = [x[j][i+1]-x[j][i] for i in range(len(x[j])-1)]
trapzy = [(y[j][i]+y[j][i+1])/2.0 for i in range(len(y[j])-1)]
A1 = sum([trapzy[i]*dx[i] for i in range(len(y[j])-1)])
A = A+A1
print A

C1=0

if type(x[0])!=list:
dx = [x[i+1]-x[i] for i in range(len(x)-1)]
trapzy = [((y[i]+y[i+1])/2.0)*((x[i]+x[i+1])/2.0) for i in range(len(y)-1)]
C = sum([trapzy[i]*dx[i] for i in range(len(y)-1)])
C1 = C1+C
else:
for j in range(len(x)):#1文字に、複数のベクトルが存在する場合対策
dx = [x[j][i+1]-x[j][i] for i in range(len(x[j])-1)]
trapzy = [((y[j][i]+y[j][i+1])/2.0)*((x[j][i]+x[j][i+ 1])/2.0) for i in range(len(y[j])-1)]
C2 = sum([trapzy[i]*dx[i] for i in range(len(y[j])-1)])
C1 = C1+C2

print("X:%s"%(C1))

C3=0
if type(x[0])!=list:
dy = [y[i+1]-y[i] for i in range(len(x)-1)]
trapzx = [((y[i]+y[i+1])/2.0)*((x[i]+x[i+1])/2.0) for i in range(len(y)-1)]
C4 = sum([trapzx[i]*dy[i] for i in range(len(y)-1)])*-1
C3 = C3+C4
else:
for j in range(len(x)):#1文字に、複数のベクトルが存在する場合対策
dy = [y[j][i+1]-y[j][i] for i in range(len(x[j])-1)]
trapzx = [((y[j][i]+y[j][i+1])/2.0)*((x[j][i]+x[j][i+1])/2.0) for i in range(len(y[j])-1)]
C4 = sum([trapzx[i]*dy[i] for i in range(len(y[j])-1)])*-1
C3 = C3+C4

print("Y:%s"%(C3))
print
"First x moment of area%s"%(C1/A)
print "First y moment of area%s"%(C3/A)
 四角せん断変形した四角45度回転させた四角穴の開いた四角ひらがなの「あ」数字の「8」
図形cen_square.pngcen_sheared_square.pngcen_45deg_rotated_square.pngcen_circle.pngcen_hollow_square.pngcen_ah.pngcen_8.png
断面一次モーメントX:500.0
Y:500.0
X:3600.0
Y:1200.0
X:1697.05627485
Y:1.42108547152e-14≒0
X:1.90819582357e-17
Y:1.11022302463e-16
X:180.0
Y:180.0
X:713082273.0
Y:475455483.5
X:419443127.5
Y:326005988.0
面積100.0100.0100.0-3.08464495357…36694613.0402222.0
図心X:5.0
Y:5.0
X:36.0
Y:12.0
X:16.9705627485
Y:1.42108547152e-16
X:6.18611169162e-18≒0
Y:3.59919225694e-17≒0
X:5.0
Y:5.0
X:1026.58929937
Y:684.489756886
X:1042.81498153
Y:810.512572659

5. 1の座標から4の図心を引いた値を作り、その図形の断面二次モーメントを計算する。

材料力学の教科書とか、ネット検索とから、
$$\int y^2\ dA$$

途中で $$\frac{1}{3} y^3$$ と、1/3が出てきたり、y3とか3乗になったりする点も気を付けたい。

積分の解き方に不安があったので、最初書いてたけど消した。

これをx,y座標から数値積分で求めたい。

先ず台形則で、積分を模擬してみた。(本当は二次の積分を台形則で模擬してはいけない気がしている。)


これまで、コードを載せてきたが、この断面二次モーメントの算出の仕方は、間違っているため、まだ載せない。

理論値からずれているのだ。

 四角せん断変形した四角45度回転させた四角穴の開いた四角ひらがなの「あ」数字の「8」
図形cen_square.pngcen_sheared_square.pngcen_45deg_rotated_square.pngcen_circle.pngcen_hollow_square.pngcen_ah.pngcen_8.png
Ix=833.333333333
Iy=833.333333333
Ix=833.333333333
Iy=4233.33333333
Ix=841.666666667
Iy=841.666666667
Ix=0.771161239361
Iy=0.771161239361
Ix=492.0
Iy=492.0
Ix=1.21943806285e+11
Iy=1.026346451e+11
Ix=79177695468.9
Iy=38517544355.2
理論値からのズレなし-ありありなし--

(ずれていると判定するのに使った理論値は来月載せます。)

理論値からのずれ(誤差)は、なぜ起こるのか。積分を模擬する手法が台形則だからだと思う。

誤差がないものは、何でないのか?


水平、垂直に線が伸びている形状に関しては誤差がない。

斜めに線が伸びている形状は、誤差が大きい。


ここで、キレイなグラフを用意して説明しようとしたのだが、時間と通信残量制限との兼ね合いで、中途半端なままアップロードする。

これに関して改善するには、台形則(一次)で表現している部分をシンプソン則(二次)で表現する必要がありそうですが、来月に続きます。


今月使用した、断面二次モーメントの算出時の積分方法(台形則)と、

シンプソン則を比較したいと思います。多分、違ってくるはず、だと思いたい。不安になってきた。


今月は**「断面二次モーメントを、座標点の配列から計算する(1)計算手順編」**だったのに対し、

来月は**「断面二次モーメントを、座標点の配列から計算する(1)積分模擬変」**にしたいと思う。変って気づいたけど直さない。

8月中に間に合わんかったし、まだとちゅうだが、これは14000文字超えの超大作なので、まあゆるせ。というか間に合ってはいた。auの通信則で制限のせいだから。おらは悪くねぇだ。長続きしていると思う。

続きを読む
posted by yuchan at 23:59 | Comment(0) | python

20190730

fontファイルの文字データ(グリフ)を、matplotlib の path (Bezier 曲線の制御点) ではなくx, y 座標に直しplotする

このページの内容は

  1. フォントファイルの文字から、ベジェ曲線の制御ポイントを取りだして、matplotlibのpathに変換する。
  2. N次ベジェ曲線の制御点→x,y座標の配列に変換
→(matplotlib.pyployt.plotでフォントが描ける!)

やる前に、IPAフォントをインストールしてね。

IPAフォントのアドレス。これをインストールする。→'C:\WINDOWS\Fonts\IPAM.TTF'が保存される。これがなかったら、下のコードは動かない。

フォントファイルのグリフから、ベジェ曲線の制御ポイントを取り出し、matplotlibのpathに変換。

ベクターフォントのフォントファイルの各文字データ(グリフ)から、ベジェ曲線の制御ポイントを取り出す方法を知った。
これで、フォントの中の文字データ(グリフ)の、ベジェ曲線の制御ポイントを取り出して作図できる。
IPAフォントを作図した。


手順は、

  1. IPAフォントのインストール。
  2. フォントのベクターの制御点の座標を、fonttoolsのrecordingpenを使って取り出す。
  3. recordingpen をmatplotlibのpathに変換する関数を書く
  4. matplotlibのpath(ベジェ曲線を描く機能)を使ってプロットする。

手順2で作ったpathデータを、matplotlibにベジェ曲線を描く機能でそのまま使った。(matplotlibのpathを使うと、制御点を指定するだけで曲線が描ける。)
わざわざ、matplotlibで読める形に直す必要があったけど、ベジェ曲線の式を理解していなくてもベジェ曲線が描けてしまった。matplotlibすごい!

#!/usr/bin/env python
# -*- coding: cp932 -*-

import time
from datetime import datetime as dt

import numpy as np

import
matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'IPAGothic'

import matplotlib.path as mpath
import matplotlib.patches as mpatches
from fontTools.ttLib import TTFont

font = TTFont(r'C:\Windows\Fonts\ipamp.ttf')
glyph_set = font.getGlyphSet()
cmap = font.getBestCmap()

def get_glyph(glyph_set, cmap, char):
glyph_name = cmap[ord(char)]
return glyph_set[glyph_name]

iroha_1 = get_glyph(glyph_set, cmap, u'て')

from fontTools.pens.recordingPen import RecordingPen

recording_pen = RecordingPen()

iroha_1.draw(recording_pen)
print(recording_pen.value)

pathdt = recording_pen.value # path(bezier curve)

def pdTopd(pd,mtrx1=[[1,0],[0,1]]):
tmp1 = [[]]
tmp2 = [[]]
i2 = 0
for i in range(len(pd)):
if pd[i][0] == 'closePath':
tmp1.append([])
tmp2.append([])
i2+=1

if pd[i][0]=='moveTo':
tmp1[i2].append(mpath.Path.MOVETO)
tmp2[i2].append(np.dot(mtrx1, pd[i][1][0]))

elif pd[i][0]=='qCurveTo':
tmp1[i2].append(mpath.Path.CURVE3)
tmp2[i2].append(np.dot(mtrx1, pd[i][1][0]))
tmp1[i2].append(mpath.Path.CURVE3)
tmp2[i2].append(np.dot(mtrx1, pd[i][1][1]))
if pd[i+1][0]=='qCurveTo':
tmp1[i2].append(mpath.Path.LINETO)
tmp2[i2].append(np.dot(mtrx1, pd[i][1][1]))

elif pd[i][0]=='lineTo':
tmp1[i2].append(mpath.Path.LINETO)
tmp2[i2].append(np.dot(mtrx1, pd[i][1][0]))

elif pd[i][0] == 'curveTo':
tmp1[i2].append(mpath.Path.CURVE4)
tmp2[i2].append(np.dot(mtrx1, pd[i][1][0]))
tmp1[i2].append(mpath.Path.CURVE4)
tmp2[i2].append(np.dot(mtrx1, pd[i][1][1]))
tmp1[i2].append(mpath.Path.CURVE4)
tmp2[i2].append(np.dot(mtrx1, pd[i][1][2]))
if pd[i+1][0]=='qCurveTo':
tmp1[i2].append(mpath.Path.LINETO)
tmp2[i2].append(np.dot(mtrx1, pd[i][1][0]))
print("%s/%s"%(i,len(pd)))

return tmp1, tmp2

a2,b2 = pdTopd(pathdt)# x,y(bezier)
a,b = pdTopd(pathdt,mtrx1=[[1,2],[0,1]])#(x,y(translated bezier))

f,(ax1) = plt.subplots(1,1)

for i in range(len(a)-1):
c = mpath.Path(b[i],a[i])
d = mpatches.PathPatch(c, facecolor='white', edgecolor="gray", lw=2, clip_on=False)
ax1.add_patch(d)
c = mpath.Path(b2[i],a2[i])
d = mpatches.PathPatch(c, facecolor='white', lw=2, clip_on=False)
ax1.add_patch(d)

ax1.set_aspect('equal')
ax1.grid(which="major", color="k", linestyle="--")
ax1.hlines(0, -3000,3000)
ax1.vlines(0, -3000,3000)
ax1.set_xlim(-2000, 2000)
ax1.set_ylim(-2000, 2000)
ax1.set_title(ur"Font", fontsize=30)
ax1.set_xlabel(ur"横 [px]", fontsize=30)
ax1.set_ylabel(ur"縦 [px]", fontsize=30)

tdatetime = dt.now()
itext = tdatetime.strftime('%Y%m%d'+ time.ctime()[11:13] + time.ctime()[14:16] + time.ctime()[17:19])
f.set_size_inches(19.2,10.8)
f.savefig("%s.png"%(itext))

20190727190844.png

N次ベジェ曲線

まずはベジェ曲線のプログラムから。

ベジェ曲線の説明には、このサイトが詳しい

あとベジェ曲線のwikipedia

二項定理とか出てきたので難しかった。
3次と4次しか使わないが、N次ベジェ曲線のx,y座標を取得するプログラムができた。

このコードは、

  1. 適当なベジェ曲線の制御点のx座標配列、y座標配列を取得して
  2. 関数部分でN次ベジェ曲線の制御ポイントをx,y座標に変換して
  3. matplotlib.pyplot.plotのx,y座標に入れる。

という作業をやっています。

N次ベジェ曲線の関数の引数は、x(N個の制御点x座標), y(N個の制御点y座標), と制御点間の補間数。(tを何分割するか。tが何なのかは、ほかのサイトを参照してください。)
当然ですが、x,yは、同じ長さにしてください。x,yの長さで、何次のベジェ曲線になるかが決まります。

x,y求めてプロット

曲線の間の補間がいまいちだとカクカクになる。

このカクカク感がたまんない!ベジェ曲線だが、ここまで間引くと、もはや曲線っぽくない。

20190601124605.png

#!/usr/bin/env python
# -*- coding: cp932 -*-

import
math
import time
from datetime import datetime as dt
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'IPAGothic'
from fontTools.ttLib import TTFont
from fontTools.pens.recordingPen import RecordingPen

def Nbezier(x1, y1, N=100):
if len(x1)!=len(y1):
print("Two lists must be the same 1-d list.")
return
if len(x1)<=2:
print("Length of the list must be grater than 2.")
return
t=np.linspace(0,1,N)
x2, y2 = [], []
for i in range(N):
tmpx, tmpy = 0, 0
n = len(x1)
for j in range(n):
tmpx = tmpx +(1-t[i])**(n-j-1)*t[i]**(j)*x1[j]*math.factorial(n-1)/(math.factorial(j)*math.factorial(n-j-1))
tmpy = tmpy + (1-t[i])**(n-j-1)*t[i]**(j)*y1[j]*math.factorial(n-1)/(math.factorial(j)*math.factorial(n-j-1))
x2.append(tmpx)
y2.append(tmpy)
return x2, y2

font = TTFont(r'C:\WINDOWS\Fonts\IPAM.TTF')#windowsのフォントファイル
glyph_set = font.getGlyphSet()
cmap = font.getBestCmap()

def get_glyph(glyph_set, cmap, char):
glyph_name = cmap[ord(char)]
return glyph_set[glyph_name]

iroha_1 = get_glyph(glyph_set, cmap, u'い')

recording_pen = RecordingPen()

iroha_1.draw(recording_pen)
pathdt = recording_pen.value

# パスデータをx,y座標データに変換する関数。
# ベジェ曲線の関数はこの関数の中で使っている。

def
PdBzToPlt(pd, N=7):
x=[]
y=[]
j=-1
for i in range(len(pd)):
tmpx=[]
tmpy=[]
if pd[i][0] == 'moveTo':
x.append([pathdt[i][1][0][0]])
y.append([pathdt[i][1][0][1]])
j+=1
elif pd[i][0] == 'qCurveTo':
tmpx.append(x[j][-1])##最初の点は、xの最終点から借りてくる。
tmpx.append(pd[i][1][0][0])
tmpx.append(pd[i][1][1][0])
tmpy.append(y[j][-1])
tmpy.append(pd[i][1][0][1])
tmpy.append(pd[i][1][1][1])
tmpx2, tmpy2 = Nbezier(tmpx, tmpy ,N)
x[j].extend(tmpx2[1:])#最初の点は、xの最終点から借りてきた点なので[1:]。
y[j].extend(tmpy2[1:])
elif pd[i][0] == 'qCurveTo':
tmpx.append(x[j][-1])
tmpx.append(pd[i][1][0][0])
tmpx.append(pd[i][1][1][0])
tmpx.append(pd[i][1][2][0])
tmpy.append(y[j][-1])
tmpy.append(pd[i][1][0][1])
tmpy.append(pd[i][1][1][1])
tmpy.append(pd[i][1][2][1])
tmpx2, tmpy2 = Nbezier(tmpx, tmpy ,N)
x[j].extend(tmpx2[1:])
y[j].extend(tmpy2[1:])
elif pd[i][0] == 'lineTo':
tmpx2 = [x[j][-1], pd[i][1][0][0]]
tmpy2 = [y[j][-1], pd[i][1][0][1]]
x[j].extend(tmpx2[1:])
y[j].extend(tmpy2[1:])
return x,y

x2,y2=PdBzToPlt(pathdt, N=7)

f,(ax1) = plt.subplots(1,1)

for i in range(len(x2)):
ax1.plot(x2[i],y2[i])
ax1.scatter(x2[i],y2[i],facecolor=[0,0,0,0.3],edgecolor=[0,0,0])#plt.show()

tdatetime = dt.now()
timtext = tdatetime.strftime('%Y%m%d'+ time.ctime()[11:13] + time.ctime()[14:16] + time.ctime()[17:19])
fname = "%s.png"%(timtext)
plt.savefig(fname)
plt.close()


来月は、フォントファイルを、x,y座標に変換し、せん断変形した後で、
x,y座標から、面積を求めようと思います。
x,y座標から、面積だけでなく、断面二次モーメントも求めようと思います。

先月に書いたシンプソン則云々は、断面二次モーメントのところで使う。8月の記事までには形にしたい。
続きを読む
posted by yuchan at 00:00 | Comment(0) | python

20190527

matplotlibのpyplot.plotで地図を描く方法

20190526185456.png

地図を描くスタンダードな方法

mpl_toolkits.basemapというのが、スタンダード。
しかし、pyplotで何とか地図を描けないだろうか。

このサイトから地図データがダウンロードできる。

国土数値情報ダウンロードサービスで、
“N03-110331_28_EC01.shp”というのを、ダウンロードしてみた。

shapeファイルをプロットしてみた

これでできた。変数eは変換行列になっているので、いろいろ試すといいかもしれない。

必要なファイル。
カレントフォルダに、”N03-110331_28_EC01.shp”

windows 日本語使用のコード

#!/usr/bin/env python
# -*- coding: cp932 -*-

import time
from datetime import datetime as dt

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'IPAGothic'

import shapefile

#下の文字列に、シェイプファイルのアドレスを入力する。
sf = shapefile.Reader(r"N03-110331_28_EC01.shp",encoding='SHIFT-JIS') #,encoding='SHIFT-JIS'必須

f,(ax1) = plt.subplots(1,1)
e = [[1, 0], [0, 1]]#変換行列
a = 133.5

try:
i = 0
while True:
shape = sf.shape(i)
x = []
y = []
x2 = []
y2 = []
for i2 in range(len(shape.points)):
x.append(shape.points[i2][0])
y.append(shape.points[i2][1])
b = shape.points[i2]
c = [a, 0]
d = [b[0] - c[0], b[1] - c[1]]
#e = [[1, 0], [0, 1]]
f2 = np.dot(d, e)
c2 = [a, 0]
f2 = [f2[0] + c2[0], f2[1] + c2[1]]
x2.append(f2[0])
y2.append(f2[1])
#ax1.plot(x, y, color="k")
ax1.plot(x, y, color="k", linestyle="--")
ax1.plot(x2, y2, color="k", clip_on=False)
i += 1
except IndexError or UnicodeDecodeError or TypeError:
pass

ax1.set_aspect('equal')
ax1.grid(which="major", color="black", linestyle="--")
ax1.plot([a,a], [25,40], color="k", linestyle="dashdot" ,clip_on=False)

ax1.set_xlim(130, 137)
ax1.set_ylim(30, 37)

ax1.set_title(ur"兵庫県の地図")
ax1.set_xlabel(ur"経度(度)")
ax1.set_ylabel(ur"緯度(度)")
#plt.show()

tdatetime = dt.now()
itext = tdatetime.strftime('%Y%m%d'+ time.ctime()[11:13] + time.ctime()[14:16] + time.ctime()[17:19])
f.set_size_inches(19.2,10.8)

f.savefig("%s.png"%(itext))
posted by yuchan at 07:00 | Comment(0) | python

20181224

計算力学技術者試験の問題集 自炊(裁断→スキャン)→リネーム、クロップ等の記録

やっと、計算力学技術者試験が終わった。12/15。固体、振動、流体とあるが、今年は振動を受けた。来年は固体。
来年に向けて、固体用のファイルを作った。
固体のファイルを作った時に、画像ファイルをリネームしたりクロップした。
携帯で読みやすくなった。

今年は、受験資格のソフトウェア使用経験を満たすため、通信講座に時間をかけた。そのため、試験勉強がおろそかになった。受かったかどうか自信はない。

裁断する

問題集は薄い。だから裁断機でも切れる。
裁断は、本文が切れたりしなければ何でもいい。
分厚い本ならば、かんなで背表紙を削ってもいい。かんながあればの話。

かんなは、案外安い。右の裁断機、2011年頃買ったら10000位したような気がした。気のせいかもしれないが。一回で数十枚しか切れないので、重労働になる。

後の加工を考えると、この時点で、余白を切り取ってしまった方がよかったのかもしれない。

スキャンする

スキャナーは、適当。
読み取り速度とかあるのかもしれないけど、よくわからない。

一枚一枚、手で紙を入れ替えるやつは、やめたほうがいい。
紙が自動的に、吸い込まれてくやつがいい。
精度は、ファイルサイズと読みやすさから検討。

三種類ある。

  1. パタパタするやつ。
  2. 吸い込むやつ
  3. 本のまま、めくって写真をとるやつ。
  4. 手で走査するやつ。
1234の順に並べた。2が高いけどおすすめ。ただ、紙の粉が詰まる。1は苦行。3は試してない。紙の粉の心配が無くてよさそうだが。
4は失敗するのが目に見えてる。昔家にあったfaxのスキャナが、4のタイプだった。手で動かすから、手があらぬ方向に動いたら全体的にゆがむ。読み取りスピード以上の速さで手を動かしたら、当然文字が飛ぶ。1のタイプでは、原稿を動かさないかぎり、そういうことが起こらない。

画像をクロップする

数百問を手作業でクロップする。
問〇-〇ってところを認識してクロップとかできればいいけど、そんなテクニックは知らない。できない。
裁断するときに、余白を切り取ってしまった方がよかったのかもしれない。

↓ サンプルとして、作ったものを載せる。文章に意味はない。

ワードで作ったような問題集だったので、ワードで作ってみた。かなり再現性高いと思う。

0001.png0002.png

ファイル名は、問題→回答順

紙の本の、問題集は、問題パートと、回答パートが分かれてしまっていた為、いちいちめくるのが大変だった。
問題→回答→問題→回答…この順番に並べ替える。
問題を奇数、回答を偶数で保存すると便利。

問題集は、順番通りのファイル名にしなければ、読みづらい。

リネームするマクロ

リネームするマクロを作った。
下のマクロは、偶数奇数を考えていない。少し書き直さないと、偶奇で保存できない。

os.rename(j, dirnm + r”/“ + “%05.f”%i + “%s”%(j)[-4:])

のiの部分を、2iとか、2i+1にすれば、何とかなるかも。よくあるやり方。

使用すると、フォルダ内の全てのファイル名が連番に変更されてしまう、恐ろしいマクロです。元に戻す方法は無い。

  1. 実行すると、フォルダ選択ダイアログが表示
  2. 選んだフォルダ内のすべてのファイルが、ファイル名順にソートされてから、5桁の0埋めの連番でリネームされる。
# coding: utf-8

import Tkinter
import tkFileDialog
import glob
import os


dirnm = tkFileDialog.askdirectory(title="Caution! carefully use.")

fnames = glob.glob(dirnm + r"/" + "*")
fnames = sorted(fnames)
lfname = len(fnames)

for i,j in enumerate(fnames):
print("%05.f"%i + "%s"%(j)[-4:])
os.rename(j, dirnm + r"/" + "%05.f"%i + "%s"%(j)[-4:])

縦長のページが、読みづらい

出来上がった画像ファイルは、スマホを横にして読みたい。
だから、画像ファイルは、横長がいい。横長にして、できるだけ大きい文字を見たい。
縦横比は、スマホにぴったりの16:9がいい。
なのに問題によっては、縦に長い画像になってしまう。
特に回答。

そんな時のために、画像をクロップするpythonのマクロを作った。

使用すると、フォルダ内の全てのjpgファイルが16:9に変更→分割されてしまう、恐ろしいマクロです。 かといって元のファイルを消すわけでもないので、そこまで怖くないかも。

# coding:utf-8

from PIL
import Image
import os, glob
import tkFileDialog


fld = tkFileDialog.askdirectory()
list = glob.glob(fld + "/" + ur"*.png") + glob.glob(fld + "/" + ur"*.jpg")

print list

for i in list:
im = Image.open('%s'%(i))
im_width, im_height = im.size

im_new_height = im_width / 16 * 9#縦のピクセル数 pixcel数の指定がしたいので、intの方が都合がいい。
last_height = im_height - im_new_height

if im_new_height*1.2 < im_height:#想定よりも1.2倍縦長ならば...
div_num = im_height / im_new_height
sft = (im_height - last_height) / div_num

if sft<im_new_height*1.7:
div_num += 1
else:
div_num += 2

shift = last_height / (div_num-1)

for i2 in range(div_num):
xl, yu, xr, yl = 0, i2*shift, im_width, im_new_height+i2*shift
im_crop = im.crop((xl, yu, xr, yl))
print i
print i2
im_crop.save('%s_%02.f%s'%(i[:-4], i2, i[-4:]))
print('%s_%02.f%s'%(i[:-4], i2, i[-4:]))
#print("%s, %s, %s, %s"%(xl, yu, xr, yl))
im_crop.close()

else:
print('%s'%(i))

im.close()
上の二つのファイルを保存したフォルダに対して実行すると、
  • 問の画像(すでに横長)に関しては、何も出てこない。
  • 解答の画像(縦長)は、分割されたファイルが出てくる。元の画像はそのまま。
  • ファイル名は00001_00, 00001_01, 00001_02, 00001_03, となっているはず。

こんな感じ。



0002_00.png 0002_01.png




0002_02.png 0002_03.png




一日20ファイル(10問)を目安に、やっていこうと思っている。

posted by yuchan at 07:00 | Comment(0) | python

20181015

pythonで、ホワイトノイズやピンクノイズのwavファイルを作る。

ホワイトノイズを作る

ホワイトノイズは、ランダムノイズ。
ランダムノイズは、フーリエ変換してもほぼランダムノイズ。(対数グラフで描くと、そうとも言い切れないが。)

「〜wavファイルを作る。」と書いたが以前書いたPythonで.wavファイルを出力(標準のwaveとscipy.io.wavfileの速度比較)の写しで、wavファイルを作れるので、そこら辺は省略。なので、ホワイトノイズに関しては、ランダムノイズのプロットだけ乗せとく。ブラウンもピンクも同様。

以前は、wavファイルにsinカーブを入れたが、今回は、その部分のみランダムノイズに変えたり、周波数をいじって加工する。
ランダムノイズは、numpyの機能を使って書いた。

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
x1=np.linspace(0, stop, Srate*stop)
A=0.5
y1 = np.random.rand(Srate * stop)-A

plt.plot(x1, y1)
plt.show()

#plt.plot(x1[:100], y1[:100])
#plt.show()

ランダムノイズ(全体).png
ランダムノイズ(ズーム).png

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f = np.linspace(1/float(stop), 44100, Srate * stop)
y11 = np.fft.fft(y1)/float(Srate * stop /2)
x1=np.linspace(0, stop, Srate*stop)

fig,((ax1,ax12),(ax2,ax22),(ax3,ax32))=plt.subplots(3,2)

ax1.plot(x1, y11)
ax12.plot(f, np.abs(y11))
ax22.plot(f, np.real(y11))
ax32.plot(f, np.imag(y11))
ax1.set_title("Input(white noise)")
ax12.set_title("abs")
ax22.set_title("real")
ax32.set_title("imag")
#plt.xscale("log")
#plt.yscale("log")
plt.tight_layout()

plt.show()
ホワイトノイズのスペクトルも、ホワイトノイズ。
全周波数にわたって、値は同程度(の振幅のノイズ)。
ホワイトノイズのスペクトルも、ホワイトノイズ.png
これがホワイトノイズの良くないところ。
人の耳は、周波数ごとに感度が違う。高周波には敏感。
どの周波数も同程度に感じるような音がピンクノイズらしい。
ピンクノイズの前に、ブラウニアンノイズを作る。

ブラウニアンノイズも作る。

ブラウニアンノイズとかの1/f^2とか1/fノイズを作るのは難しい。

wikipediaのピンクノイズのページを見ると、「パワーが周波数に反比例する雑音のこと。」とある。

ピンクノイズを作るための方法として、voss法というのがあるらしい。
しかし、あんまり難しいのは避けたい。「新しい勉強をするのも時間的にちょっと…。」と思ったので諦めた。

そこで、思いついたやり方でやる。

ホワイトノイズの信号を使って

ホワイトノイズをfft→周波数空間で加工→ifftで元に戻す。

もしこのやり方に名前がないのであれば、「fourie (変換おじさん)法」とでもしておく。

なぜかホワイトノイズのままになってしまう

「間違いなく作ったのに、出てくる音はホワイトノイズのそれなんだよなぁ…」と思って、スペクトル調べたら、ホワイトノイズだった。

そのあと、原因がわからなくて、四苦八苦した。
原因は、スペクトルの折り返しを消してしまっていたからだった。

スペクトルは、折り返しの所で、傾きが逆になる。


1/f^0=1(無処理:ホワイトノイズ)と1/f^2(ブラウンノイズ)の間がピンクノイズらしい。

ナイキスト周波数の所で折り返すので、fmax/2の所で、傾きを逆にしなければならない。

とりあえず、ブラウンノイズから作りたい。
ブラウンノイズっぽく1/f^2で減っている線を描く。

import matplotlib.pyplot as plt
import numpy as np


stop = 10
Srate = 44100
f = np.linspace(1/10.0, 44100, Srate * stop)
f1 = 1.0/f**2
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]


plt.plot(f,f1)
plt.yscale("log")
plt.xlabel("frq. in log scale")
plt.xscale("log")
plt.ylabel("mag. in log scale")
plt.show()

1/f^2.png

logスケールで直線になった。
ただ、1Hz以下の部分が大きすぎる。
あと、後述するが、20Hz以下はいらんので、20Hzの所が1になるようにしたい。

import matplotlib.pyplot as plt
import numpy as np


stop = 10
Srate = 44100

f = np.linspace(1/10.0, 44100, Srate * stop)
f1 = 1.0/f**2
v = f1[20*stop]
f1 = f1/v
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]


plt.plot(f,f1)
plt.yscale("log")
plt.xlabel("frq. in log scale")
plt.xscale("log")
plt.ylabel("mag. in log scale")
plt.show()

1/f^2の修正.png

これを、ホワイトノイズにかけ合わせたい。

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f**2
v = f1[20*stop]
f1 = f1/v
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
y11 = np.abs(np.fft.fft(y1)*f1)/float(Srate * stop /2)


plt.plot(f, y11)
plt.xscale("log")
plt.yscale("log")
plt.show()

1/f^2を、ホワイトノイズのスペクトルにかけた.png
このスペクトルは、このままだと、まずい。
まずい理由は、

  • 可聴域の周波数で、強度が低すぎる。
  • 低周波が強すぎる。(どうせスピーカーでは再生できない周波数が強く出すぎて、ほかの音が小さい。これはこれで、ピンクノイズではあるのだろうけど。)

試しに、ifftかけると

y11 = np.fft.ifft(np.fft.fft(y1)*f1)
plt.plot(np.real(y11))
plt.show()

低周波が優勢すぎる.png
のっぺりしてしまっている。

figure_2-2.png

スピーカーの再生できる周波数に絞る(20〜20kHz位だと思っている。)

スペクトルの、赤で囲った領域はいらないので、0で置換。
それ以外の周波数は、20 Hzの値が、1以下に収まるようにする。

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f**2
v = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.abs(np.fft.fft(y1)*f11)

plt.plot(f, y11/float(Srate*stop/2))
plt.xscale("log")
plt.yscale("log")
plt.show()

可聴域のみに加工.png

こうなった。これを、上げる。

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f**2
v = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.fft.fft(y1)*f11
y11 = y11/abs(y11[20*stop])*float(Srate*stop/2)
y11 = np.abs(y11)/float(Srate*stop/2)

plt.plot(f, y11)
plt.xscale("log")
plt.yscale("log")
plt.show()

可聴域のみに加工 高さを合わせた.png

逆変換すると、

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
x1=np.linspace(0, stop, Srate*stop)
y1 = np.random.rand(Srate * stop)-A
f = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f**2
v = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.fft.fft(y1)*f11
y11 = y11/abs(y11[20*stop])*float(Srate*stop/2)
y11 = np.real(np.fft.ifft(y11))/stop/2.0

plt.plot(x1, y11)
plt.show()

逆フーリエ変換した(全体).png

逆フーリエ変換した(ズーム).png

 

1より大きくなっちゃった。(泣)

ピンクノイズは、

周波数的には

import matplotlib.pyplot as plt
import numpy as np


stop = 10
Srate = 44100
f = np.linspace(1/10.0, 44100, Srate * stop)
f1 = 1.0/f
v = f1[20*stop]
f1 = f1/v
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]

plt.plot(f,f1)
plt.yscale("log")
plt.xlabel("frq. in log scale")
plt.xscale("log")
plt.ylabel("mag. in log scale")
plt.show()

1/f.png


できてきたノイズのスペクトルは、

import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
y1 = np.random.rand(Srate * stop)-A
f = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f
v = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.fft.fft(y1)*f11
y11 = y11/abs(y11[20*stop])*float(Srate*stop/2)
y11 = np.abs(y11)/float(Srate*stop/2)

plt.plot(f, y11)
plt.xscale("log")
plt.yscale("log")
plt.show()

1/fを、ホワイトノイズのスペクトルにかけ可聴域に加工.png

ピンクノイズは、ブラウニアンよりも減りがなだらか。
import numpy as np
import matplotlib.pyplot as plt


stop = 10
Srate = 44100
A=0.5
x1=np.linspace(0, stop, Srate*stop)
y1 = np.random.rand(Srate * stop)-A
f = np.linspace(1/float(stop), 44100, Srate * stop)
f1 = 1.0/f
v = f1[20*stop]
f1 = f1/v
f1[:stop*20]=0
f1[stop*20000:]=0
f1 = np.r_[f1[:len(f1)/2],f1[1:len(f1)/2+1][::-1]]
f11 = f1
y11 = np.fft.fft(y1)*f11
y11 = y11/abs(y11[20*stop])*float(Srate*stop/2)
y11 = np.real(np.fft.ifft(y11))/stop/2.0

plt.plot(x1, y11)
plt.show()

逆フーリエ変換した(ピンクノイズ).png

これも1より大きくなっちゃった。

こうしてできた波形を、wavファイルに入れるのであった。。。

この後、wav作ったけど…

この後、wavファイルにした。
したけど、、、どうも、これだけじゃダメなようだ。信号の最大値が高すぎて、破裂音がする。

もう少し、全体のレベルを調節しないと。

あと、「結局聞いているのって、高周波音」というのがある。高周波がないと、全然聞こえないし。

しばらく書かないつもりだったが、せっかく思いついたので書いておく。

posted by yuchan at 07:00 | Comment(0) | python

20180709

脳ドッグに行ってきた。→MRIの画像データをもらった。(DICOMファイル)→k-spaceは入ってなかった。

先に断っておく。

  • 医療情報っぽいことを書いたが、自分は医療従事者でもなんでもない。
  • 嘘を書いたつもりはない。しかし、間違っていても、責任は持てない。

脳ドッグに行ってきた。→MRIの画像データをもらった。(DICOMファイル)→残念ながら、k-spaceは入ってなかった。
しょうがないので、画像データを二次元フーリエ変換してみた。

MRIとは

MRI(Magnetic Resonance Imaging: 磁気共鳴画像)のこと。
装置は、google検索でトップだった、このサイトの画像のまさにこれ。これの中に入ってきた。

MRIの結果とフーリエ変換

MRIによって、脳とか、心臓の断面画像が手に入る。

この画像の元データがk-spaceとかいうデータで、
これを二次元フーリエ変換すると、断面画像ができるということで、これが欲しかった。
(逆二次元フーリエではないらしい。空間周波数のデータを、フーリエ変換するって何かいつもと違うのがいやだけど。)
脳ドッグでは、MRIを使うそうな。
だから、無保険(16万自腹)でちょっと高めの脳ドッグを受けてみた。

しかし、遊び目的ではない。

たぶんTIA(Transit Ischemic Attack:一過性脳虚血発作)ってやつにかかったかもしれないと思う出来事があった。
だから、脳の精密検査をしたかった。こっちが本当はメイン。

先日、ジムに行っていたときのこと。
週1で、トレッドミル(ランニングマシン)で10km走るのが習慣になっていた。10.5 km/h
換気が悪いなぁと思っていた。

そしたら汗をいつも以上にかいてしまったらしく、脱水症っぽい感じになった。
さらにひどいことに、走り終わった後、手がしびれ、ろれつが回らない状態になった。
こんなん初めて。

ポカリをたくさん飲んだら治った。
先月は、これについて書こうとしていたが、pythonの話とは毛色が違い、話がまとまらなかったため、投稿を断念した。上記のようなことを書こうとしていた。

脳が悪いのか!?と思って怖くなった。

#7119(神戸で救急車が必要か判断できない人のための電話番号)に電話したら、看護師さんにつないでくれて、事情を、上記以上に話して、何個か簡易テスト的なのを指示されて行ったら。「緊急性はないですね。けど、早いうちに医療機関に受診されたほうがいいですよ」とのこと。

もともと、k-space画像が欲しかったので、いわれなくても受けるつもりではあった。
しかし、もともと脳ドッグを予約していた日時(2か月くらい先)を言うと、
看護師さんがちょっと訂正するような口調で「っ、もうちょっと、早いほうが、いいですね」と言われてしまった。

で、直近の土曜日に、脳ドッグを受けた。

その結果が、3週間後(今)に届いた。
結果は、大丈夫だった。
どっか、血管が詰まっているとかはなかった。

ただ、そのほか健康面(コレステロールとか血圧とか)で、所見があったので、結果として行ってよかった。

MRIした感想

  • 狭い
  • 装置の音がすごかった。
  • ヘッドフォン装着→音楽は選べた。
  • 体は拘束された状態だった。

  • 「息を吸ってー、、、息を止めて。」と何回か指示が来た。

  • 「血管を拡げるげるお薬」なるスプレーを舌にかけられた。→寝すぎた時の頭痛と同じ状態に。
    スプレー1吹きで、体にあんな変化があるなんて、「医療すごい」ってなった。

DICOMファイル

医療情報の画像データの規格らしい。
画像のみならず、音声データ等も保存できるそうな。
その他、情報がたくさん詰まっている。

「これに、k-spaceが!?」と思い、
メモ帳で開くと、ところどころ文字化け。バイナリファイルなんでしょう。きっと。

DICOMファイルは、アップロードしない。
なぜならば、

  • 自分の個人情報
  • 受診した病院名 所在地
  • オペレーターの氏名

などが書いてあったから。
自分のは情報はいいとして、自分以外の情報を公開するのは気が引ける。
たとえやっても法律上問題ないとしても気が引ける。

なので、結果画像だけ読み込んでみた。
そしたら、すでに画像化されていた。
k-spaceのデータは含まれていないようです。

まぁそりゃそうだよなぁ。
医療従事者が二次元フーリエ変換なんで遊ぶ必要ないもん。
患者さんを助ける情報が必要なだけだもん。

というか、オペレーターの方が、「k-spaceはあげられないんですよぉー」(にこっ)
っておっしゃっていた気がする。いまさら思い出す。
MRIの技師の方は、二次元フーリエ変換知ってた。やはり使うんだろうか。

dicomファイルの中に、機械の製造者の情報(?)が、Philipsとなっていた。
Philipsといえばシェーバーの会社じゃないですか。
髭剃りだけじゃなかったんだな。
https://ja.wikipedia.org/w/index.php?title=%E3%83%95%E3%82%A3%E3%83%AA%E3%83%83%E3%83%97%E3%82%B9&oldid=69059425

pydicom

pythonでdicomファイルを開くには、このpydicomがいる。

  1. dicomファイルを読む(画像のnumpy arrayを含むデータの読み込み)
  2. matplotlibで、numpy arrayをimshow
    という、いつもの流れだった。

PythonでDICOM画像をなんとかするここが詳しい。

cmdで、全部小文字で、

pip install pydicom

だった気がする。間違ってたらごめん。
当然ですがpipが入ってないと、効かない。

環境

python2.7
pydicom, numpy, matplotlib, tkFileDialog

コード

  • このプログラムは、DICOMファイルを開いて、画像化するプログラムです。
  • 二次元フーリエ変換した画像も併せて載せていますが、これは(多分)k-spaceではありません。
  • DICOM以外のファイルを開いた時の挙動は意図していません。
import pydicom as dc
import numpy as np
import matplotlib.pyplot as plt


import
tkFileDialog as tk
name = tk.askopenfilename(filetypes=[("dicom",("all","*.*"))])
d = dc.read_file(name)


mx = np.max(np.real(np.fft.fft2(d.pixel_array)))/100000.

f,((ax1,ax2),(ax3,ax4))=plt.subplots(2,2)

ax1.imshow(d.pixel_array, cmap="gray")
ax2.imshow(np.abs(np.fft.fft2(d.pixel_array)), vmin=0, vmax=mx*1000, cmap="gray")
ax3.imshow(np.real(np.fft.fft2(d.pixel_array)), vmin=-mx, vmax=mx, cmap="gray")
ax4.imshow(np.imag(np.fft.fft2(d.pixel_array)), vmin=-mx, vmax=mx, cmap="gray")

ax1.set_title("pydicom pixel_data")
ax2.set_title("2DFFT absolute value")
ax3.set_title("2DFFT real part")
ax4.set_title("2DFFT imaginary part")

plt.show()

プログラムの実行方法

上記コードを、テキスト帳にコピーし、「a.py」というファイル名で保存。

そのファイルが表示されているブラウザ上で、普段アドレスが表示されているところに、

cmd

と入力し、Enterキーを押す。
(もしくは、ブラウザ上の何もないところを、shift+右クリックして、「PowerShellWindowをここに開く」をする。)

そうすると、黒い画面(青い画面)が出てくる。

その画面上で、

python a.py

と入力すると、実行できる。ファイルダイアログが出るはず。
(pythonのフォルダが、環境変数で適切に設定してあればの話。)dicomファイルを選ぶと、画像が表示される。


figure_1.png

  • 画像データはnumpy arrayで返ってくる
  • 画像は、mm感覚で撮影されているらしい
    なので、3次元でプロットすれば、立体的に見れるはず。
    大変そうだからまだやらない。
posted by yuchan at 07:00 | Comment(0) | python

20180528

matplotlibのimshowで円を描く。

02.gif

### 一ピクセルづつ、円を描く。
matplotlibのimshowで円を描く。
いつもはplotで描く。この時は、x座標、y座標を指定する。
その座標は、三角関数を使って用意する。((cosθ, sinθ))
ただ、imshowは配列の用意の仕方が違う。配列のindexが、座標そのものなので戸惑う。

「二次元の配列上に、一次元のプロットと同じ感覚でプロットしたい。」という要求。
↑これは、できなかったけど、まあ何とか形になったし、別にいいか。



### 作戦
下図参照
20180526113622 - コピー.png


結構、面倒だった。
これの右側の話。どうやって、ピクセルのインデックスを指定するのか、ぱっと思いつかなかったので、描いた。
これまでのことを考えると、当たり前のことだったが、、、
- x,yの一次元の配列があって…
- 二次元の場合に当てはめると、インデキシングが、x,yの値に該当して…
- インデキシングで値を指定しようとすると、0〜Nの数字になるからマイナスはどうなっちゃうの?
- そもそも、値1と2の間に、1.9とか1.5とかあるし、imshowの場合は、1と2の間がない。(そういう風に軸を整えれば作れるが。)インデキシングはできない。
とかぐるぐる考えていた。

### コード
```py
import numpy as np
import matplotlib.pyplot as plt
import time
from datetime import datetime as dt



N,M = 32, 64
a=N/2-1


theta = np.linspace(0, 2.0*np.pi, min(N,M)*3)
x = np.cos(theta)*a+a
y = np.sin(theta)*a+a


x2 = np.around(x).astype(int)
y2 = np.around(y).astype(int)



img = np.zeros((N,M))
for i in range(len(x2)):
    img[x2[i]][y2[i]] = 1
    



fig, (ax1, ax2) = plt.subplots(1,2)
ax1.plot(x,y)
ax2.imshow(img, cmap="gray",interpolation="None")


ax1.set_aspect("equal", "datalim")



tdatetime = dt.now()
fname = tdatetime.strftime('%Y%m%d'+ time.ctime()[11:13] + time.ctime()[14:16] + time.ctime()[17:19])
plt.savefig("%s.png"%(fname))
```



### 作っていて気付いたこと
##### 円周を、何ピクセル使って表現するか
円を描くとき、角速度が0〜2パイまでを使った。
これを、np.linspaceで**n分割**したの配列にして、変数θとした。
これを、cos θ, sin θに入れた。

n分割のところをどんな値にすれば適切なのか。
- 少ないと、円が途切れる。
- 多いと、単に無駄。
だいたい3で大丈夫。
20180526171408-001-003.png

### サイトの事
- ある記事がバズってる。アクセスカウンタがすごい。
「桁多めに設定しすぎたかなぁ。『こんなブログ、そんなに桁必要ねーよ』とか思われてたら恥ずかしいなぁ…(被害妄想)」と思ってたが、今となっては、「あと1桁くらいあってもいいかな。」と思ってる。
- さくらのブログのアクセスカウンタで、棒グラフが長くなると、積み上げグラフ→単独のグラフに代わる。
初めて見た!
さくらのブログ アクセスカウンタ 棒グラフ.png
- マークダウンをせずに書いても、それなりに読める。むしろ読みやすい。今回だけこのまま放置。
- これだけ、アクセスがあっても、収入はゼロ。
どうなってんの?置き方悪いの?広告が適切でないとか?

バズっているところに、内容にドンピシャで、自分も一年間の使用経験のある、商品の広告を置いてみるなどしてみた。
結果を様子見することにする…
posted by yuchan at 07:00 | Comment(0) | python

20180409

matplotlibの、cmapを、徐々に変化させる方法。

matplotlibのcmapを変化させる。

下記gif画像のように、自然に、徐々に、cmapを別のcmapへ移行させたい。

04.gif

何のためにやるのかと問われたら、

もちろん、アニメーションのエフェクトが目的。

フーリエの肖像画(著作権に関しては、スミソニアンライブラリーを参照。)のcmapを、gray(この肖像画の方は)→jet(熱伝導の研究をされていた方で、)→seismic(からの〜。)→bwr(フランス人)に徐々に変化させるという、アートなエフェクトがしたい。(このエフェクトは、「飽きさせない工夫」程度のもの。)

普通に使ってたら、こんな需要はないだろう。検索しても、出てこなかった。かなり苦労した。
cmapが変化したら、変だもんな。ここは本来なら、変えちゃだめだよね。

普通、cmapの指定は

文字列で指定する。

import matplotlib.pyplot as plt
import matplotlib.image as mpimg

img = mpimg.imread("hogehoge.png")
plt.imshow(img, cmap="jet")

plt.show()

cmapを作るのは、以前やった。
ここでは、cmapを指定するために、cdictという所をいじった。

cmapは、文字列で指定するやり方以外に、変数っぽく指定するやり方がある。(ほかの(既存の)cmapは文字列で指定する。)
こんな風に使う。

import matplotlib.pyplot as plt
import matplotlib.image as mpimg

img = mpimg.imread("hogehoge.png")
c_gray = plt.get_cmap("gray")
plt.imshow(img, cmap=c_gray)

plt.show()

…みたいな。

cmapの自作には、下の関数を使う。詳しくは、最後のコードを参照してください。

matplotlib.colors.LinearSegmentedColormap('cmap%s'%i, cdict, 256)

自作時に、色の指定は、cdictって所で指定する。

ざっくり仕組み。

cdictは、R,G,B別の、(値(x), 色の初期値(y0), 色の終わりの値(y1))のタプルで構成される。
(x, y0, y1)の数は、最低二つ。

import matplotlib.pyplot as plt

c_gray_r = plt.get_cmap("gray_r")
c_jet = plt.get_cmap("jet")

print c_gray_r.__dict__['_segmentdata']
print c_jet.__dict__['_segmentdata']

>>> import matplotlib.pyplot as plt
>>>
>>> cgrayr = plt.getcmap(“grayr”)
>>> c_jet = plt.get_cmap(“jet”)
>>>
>>> print c_gray_r.__dict
[‘_segmentdata’]
{u’blue’: [(0.0, 1, 1), (1.0, 0, 0)], u’green’: [(0.0, 1, 1), (1.0, 0, 0)], u’red’: [(0.0, 1, 1), (1.0, 0, 0)]}
>>> print c_jet.__dict
[‘_segmentdata’]
{u’blue’: ((0.0, 0.5, 0.5), (0.11, 1, 1), (0.34, 1, 1), (0.65, 0, 0), (1, 0, 0)), u’green’: ((0.0, 0, 0), (0.125, 0, 0), (0.375, 1, 1), (0.64, 1, 1), (0.91, 0, 0), (1, 0, 0)), u’red’: ((0.0, 0, 0), (0.35, 0, 0), (0.66, 1, 1), (0.89, 1, 1), (1, 0.5, 0.5))}

cmapが、gray_rの場合は、すべての色のタプルが、二つしかない。
しかし、jetの場合は、青が5、緑が6、赤が5種類もある。
cmapによって、また、色によって、タプルの個数が異なっていた。

jetのようにカラフルなcmapだと、タプルの個数が多くなっている
これらを、どこだかで、線形補完して、徐々にcmapが変化しているアニメを作る。

大まかに書くと、
二つのcmapのタプルの数をそろえる。
→両者の各色の値を、線形補完した値で、タプルを作りなおす。

もう少し細かく書くと、

  1. タプルから、xだけ取り出したリストを作って配列化。
  2. その配列をソートして
  3. ダブっているのを排除。
  4. このxのリストに従って、線形補完。まずは、各cmapの中だけで、$(x_0 , y0_0 , y1_0 )\to (x_1 , y0_1 , y1_1 )$ を線形補完。そのxが無かったら、そのxより大きいxを持つ一つ上のタプルから拝借して、dx分を加味して線形補完。(ここら辺がごちゃごちゃして、読みにくくなってしまった。)
  5. 各色のcdictのタプルの長さが等しくなった二つのcmapを、分割したいフレーム数で、線形補完。

4-5の辺は、ガチャガチャやったらなんかうまくいったという感じになっている(下記、汚コードを参照)。
「ここが、こうなってー」とかの説明を求められたら、答えに窮する。
アニメのエフェクトに関して言えば、「こんなの動けばいいんだ」といいたい。

使い方

環境(バージョン)
  • python 2.7
  • numpy-1.13.1
  • matplotlib-1.5.3
  • mpl_toolkits

  • tkFileDialog

ちょっと位ちがっても、動く気がする。

usage

徐々に変化するcmap(の変数が収まったdict変数)を作る、関数にした。

引数
  • cmap1 matplotlib.colors.LinearSegmentedColormap object
  • cmap2 matplotlib.colors.LinearSegmentedColormap object
  • frames 初期値30。30個のcmapに分割する。0-29(30個)で指定する。こんな感じ。cmap=cmaps1[“%s”%i]
返り値

framesと同じ数のcmap。
dictの中に、cmap変数を入れているので、0〜(frames-1)(引数frames個)で指定する。

for i in range(fps):
c = cmap_shift(c0=plt.get_cmap("gray"), c1=plt.get_cmap("jet"))
ax1.imshow(img, cmap=cmaps1["%s"%i])
fig.savefig("%s.jpg"%i)
plt.close()

こうして作ったcmapを、また、ガチャガチャ組み合わせて、使うのであった。

とりあえず、僕の作った汚コードを載せる。
テキスト帳で「hogehoge.py」と保存して、
cmdで、「python hogehoge.py」とか打って、使ってほしい。
pythonインタープリタで動くかどうかは、試してないのでわからない。

#!/usr/bin/env python
# -*- coding: cp932 -*-

import os
import time
from datetime import datetime as dt
import numpy# as np
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

from matplotlib import ticker
from mpl_toolkits.axes_grid1 import make_axes_locatable

##################################
#ダイアログで、画像ファイルを開く#
##################################
import tkFileDialog as tk
name = tk.askopenfilename(filetypes=[("image",("*.png","*.bmp","*.jpg")),("all","*.*")])

img = mpimg.imread(name)
print(type(img[0][0]))

if isinstance(img[0][0],numpy.ndarray):#カラーなら平均してグレー化。グレー(白黒)ならそのまま。
img = numpy.mean(img,-1)#20180409 a.m.1:59となりのベランダから、ライターを付けるような音がした。ボヤとかは、勘弁してくれよ。たまに、雄たけび、叫び声、笑い声(「イーヒッヒ」)とかが、(それも深夜に)聞こえるが、犯罪は本当に勘弁してくれ。なんだか、ガス代を滞納してるのに、厚かましく自分で開いて使ってるみたいだし。どうしようもない奴なのは確定している。「会って話すと意外と"いい奴"かも」とか言う奴がいるが、おれの定義する"いい奴"は、ガス代を滞納したりしない。早々に引っ越したいが、いろいろあって安いので、この部屋に甘んじている。身の危険を感じたら、出ていきたい(本当は、身の危険を感じる前に出たいが)。今作っている大作を、なんとしてでも完成させたい。完成まで、多分あと一年。それまで、絶対生きねばならん。そのあとだって、アイディアだけは、たんまりあるんだ。

def remove_dup(listarg):
x=[]
for i in listarg:
if i not in x:
x = x+[i]
return x

##########################################################
#徐々に変化するcmap群を、返す関数cmap_shift関数を自作する#
##########################################################
def cmap_shift(c0, c1, frames=30):
red = {}
green = {}
blue = {}
colrow={}
dxy={}
for j in ["red","green","blue"]:
x = [c0.__dict__['_segmentdata']["%s"%j][j2][0] for j2 in range(len(c0.__dict__['_segmentdata']["%s"%j]))] + [c1.__dict__['_segmentdata']["%s"%j][i1][0] for i1 in range(len(c1.__dict__['_segmentdata']["%s"%j]))]
x = sorted(x)
x = remove_dup(x)
xy1 = []
xy2 = []
ic0=0
ic1=0

for i02 in x:
if i02 >= c0.__dict__['_segmentdata']["%s"%j][ic0+1][0]:
ic0+=1
if i02 == c0.__dict__['_segmentdata']["%s"%j][ic0][0]:
y0 = c0.__dict__['_segmentdata']["%s"%j][ic0][1]
else:
try:
y0 = c0.__dict__['_segmentdata']["%s"%j][ic0][1] + ((c0.__dict__['_segmentdata']["%s"%j][ic0+1][1] - c0.__dict__['_segmentdata']["%s"%j][ic0][1]) / (c0.__dict__['_segmentdata']["%s"%j][ic0+1][0] - c0.__dict__['_segmentdata']["%s"%j][ic0][0]))*((i02-c0.__dict__['_segmentdata']["%s"%j][ic0][0]))
except IndexError:
y0 = c0.__dict__['_segmentdata']["%s"%j][-1][1]
if i02 == c0.__dict__['_segmentdata']["%s"%j][ic0][0]:
y1 = c0.__dict__['_segmentdata']["%s"%j][ic0][2]
else:
try:
y1 = c0.__dict__['_segmentdata']["%s"%j][ic0][2] + ((c0.__dict__['_segmentdata']["%s"%j][ic0+1][2] - c0.__dict__['_segmentdata']["%s"%j][ic0][2]) / (c0.__dict__['_segmentdata']["%s"%j][ic0+1][0] - c0.__dict__['_segmentdata']["%s"%j][ic0][0]))*((i02-c0.__dict__['_segmentdata']["%s"%j][ic0][0]))
except IndexError:
y1 = c0.__dict__['_segmentdata']["%s"%j][-1][2]
xy1.append((i02, y0, y1))
if i02 >= c1.__dict__['_segmentdata']["%s"%j][ic1+1][0]:
ic1+=1
if i02 == c1.__dict__['_segmentdata']["%s"%j][ic1][0]:
y0 = c1.__dict__['_segmentdata']["%s"%j][ic1][1]
else:
try:
y0 = c1.__dict__['_segmentdata']["%s"%j][ic1][1] + ((c1.__dict__['_segmentdata']["%s"%j][ic1+1][1] - c1.__dict__['_segmentdata']["%s"%j][ic1][1]) / (c1.__dict__['_segmentdata']["%s"%j][ic1+1][0] - c1.__dict__['_segmentdata']["%s"%j][ic1][0]))*((i02-c1.__dict__['_segmentdata']["%s"%j][ic1][0]))
except IndexError:
y0 = c1.__dict__['_segmentdata']["%s"%j][-1][1]
if i02 == c1.__dict__['_segmentdata']["%s"%j][ic1][0]:
y1 = c1.__dict__['_segmentdata']["%s"%j][ic1][2]
else:
try:
y1 = c1.__dict__['_segmentdata']["%s"%j][ic1][2] + ((c1.__dict__['_segmentdata']["%s"%j][ic1+1][2] - c1.__dict__['_segmentdata']["%s"%j][ic1][2]) / (c1.__dict__['_segmentdata']["%s"%j][ic1+1][0] - c1.__dict__['_segmentdata']["%s"%j][ic1][0]))*((i02-c1.__dict__['_segmentdata']["%s"%j][ic1][0]))
except IndexError:
y1 = c1.__dict__['_segmentdata']["%s"%j][-1][2]
xy2.append((i02, y0, y1))
dxy["%s"%j]=[[(xy2[i][ii]-xy1[i][ii])/float(frames) for ii in range(3)] for i in range(len(xy1))]
colrow["%s"%j]={}
for i03 in range(frames):
colrow["%s"%j]["%s"%i03] = tuple([tuple([xy1[i02][i04] + dxy["%s"%j][i02][i04]*i03 for i04 in range(3)]) for i02 in range(len(x))])

cmaps={}
for i in range(frames):
cdicts_tmp={}
cdicts_tmp["red"] = colrow["red"]["%s"%i]
cdicts_tmp["green"] = colrow["green"]["%s"%i]
cdicts_tmp["blue"] = colrow["blue"]["%s"%i]
cmaps["%s"%i] = matplotlib.colors.LinearSegmentedColormap('cmap%s'%i, cdicts_tmp, 256)
return cmaps

########################################
#自作したcmap_shift関数で、cmaps1を作る#
#c0に"gray"、c1に"jet"を指定する。 #
#コマ数は、5 #
#cmaps1は、dict型で、返ってくる #
########################################
fms=5
cmaps1 = cmap_shift(c0=plt.get_cmap("gray"), c1=plt.get_cmap("jet"),frames=fms)

for i in range(fms):
fig,(ax)=plt.subplots(1)
##############################
#cmaps1は、ここで、使っている#
##############################
imgcol = ax.imshow(img, cmap=cmaps1["%s"%i], interpolation="none", clip_on=False)
divider = make_axes_locatable(ax)
ax_cb = divider.new_horizontal(size="10%", pad=0.05)#
fig.add_axes(ax_cb)
cb = plt.colorbar(imgcol, cax=ax_cb)
tick_locator = ticker.MaxNLocator(nbins=5)
cb.locator = tick_locator
cb.ax.yaxis.set_ticklabels([-1.0,-0.5,0,0.5,1.0])
cb.update_ticks()
ax.locator_params(axis='both', nbins=4)
cb.set_label("\"gray\" to \"jet\"")
tdatetime = dt.now()
itext = tdatetime.strftime('%Y%m%d'+ time.ctime()[11:13] + time.ctime()[14:16] + time.ctime()[17:19])
#fig.set_size_inches(19.2,10.8)
fig.savefig("%s%s.png"%(itext,i))
plt.close()

print("Files are saved on %s"%os.getcwd())

自分の書いたコードでは、gray→jetで固定してあるが、

l130-l140位にある、

cmaps1 = cmap_shift(c0=plt.get_cmap("gray"), c1=plt.get_cmap("jet"),frames=fms)

の所の”gray”とか、”jet”とかを、手動で指定すれば色を色々変えることができる。

cmapとして指定できる、色の一覧は、こことかを参考にしたり、いろいろやって自作したりしてほしい。

posted by yuchan at 07:00 | Comment(0) | python

20180305

matplotlibのmake_axes_locatableで作ったcolorbarの透明度を変更する

color bar の軸の色を調節する。

matplotlibのcolorbarをいじっている。

03.gif

colorbarの高さを、軸の高さに合わせるために、make_axes_locatableを使った。
make_axes_locatableを使った後、標準のcolorbarとは、使い勝手が変わっているように感じる。
とりあえず、make_axes_locatableを使った後の、cbarの軸を透明化する方法を、書き残しておく。


# coding: utf-8
import math
import time
from datetime import datetime as dt
import numpy as np
from PIL import Image
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib import ticker

N=16
img1 = [[math.sin(2 * np.pi * i * j /float(N)) for i in range(N)] for j in range(N)]


fig = plt.figure()
loc_ax = [0.1,0.2,0.8,0.7]
ax1 = fig.add_axes(loc_ax)
imgcol1 = ax1.imshow(img1, vmax=1, vmin=-1, cmap = 'gray', interpolation='nearest')
ax1.set_ylabel("frq m [Hz]")
ax1.set_xlabel("horizontal px h [px]")

divider1 = make_axes_locatable(ax1)
ax_cb0 = divider1.new_horizontal(size="10%", pad=0.05)
ax_cb01 = fig.add_axes(ax_cb0)
cb1 = plt.colorbar(imgcol1, cax=ax_cb0)
tick_locator1 = ticker.MaxNLocator(nbins=5)
cb1.locator = tick_locator1
cb1.update_ticks()
ax1.locator_params(axis='x', nbins=4)
ax1.locator_params(axis='y', nbins=7)
cb1.set_label("real part")
ax1.set_axis_on()

#############
#insert here#
#############

tdatetime = dt.now()
itext = tdatetime.strftime('%Y%m%d'+ time.ctime()[11:13] + time.ctime()[14:16] + time.ctime()[17:19])
fig.canvas.draw()
# save it to a numpy array.
data = np.fromstring(fig.canvas.tostring_rgb(), dtype=np.uint8, sep='')
data = data.reshape(fig.canvas.get_width_height()[::-1] + (3,))
im = Image.fromarray(data)
im.save("%s.png"%(itext))
plt.close()
im.close()

20180304144407.png

これから書くコードの断片は、上記コードの、

#############
#insert here#
#############

って所に、挿入するという形で進める。

禁止

  • 色付き画像を重ねる方法は禁止。
画像の上に、色調節したポリゴンか何かを重ねて、透明っぽくするのは禁止。
たとえば、imshowを使って、背景色と同じ色を重ねて、透明っぽくするとかは禁止。
あくまでも、各パーツを透明にする。
#############
#insert here#
#############
img2 = [[[191,191,191,50] for i in range(5)] for j in range(5)]
loc_ax2 = [
0.05,0.1,0.9,0.8]
ax2 = fig.add_axes(loc_ax2)
ax2.patch.set_facecolor([
0,0,0,0])#何も書いてない場合、背景は透明(ディフォルトでは白)。
imgcol2 = ax2.imshow(img2, vmax=1, vmin=-1, cmap = 'gray', interpolation='nearest')
こうやって、透明っぽくす誤魔化すのは無し。
(軸はあえて隠さなかった。)
20180304144936.png
アニメーション時の汎用性が下がる。
  • set_visibleで消すのも禁止。
「徐々に出現する」とか、「徐々に色が変わる」とかのアニメーションがしたい。
だから、True, Falseで一発で点いたり消えたりしてしまうのは無し。

各パーツ ax, spines, ticks(major, minor)

自分の語彙が無いせいだが、「軸」といったら、黒いところ全体を指して言っている。
matplotlibのドキュメントを見ると、ax, spines, ticks(major, minor)とか書き分けてある。
これ全部を透明にしたりして、出現するアニメーションをしたい。

20180304144407.png

colorbarのaxのspinesの色を変える。

20180304150030.png

色の指定は[R,G,B,A]。

color1=[0,0,0,0.0]
cb1.outline.set_edgecolor(color1)

20180304145512.png

透明にしたが、color1の部分を

color1=[1,0,0,0.8]

とかにすれば、当然着色もできる。

ticks 目盛りの部分。

20180304150707.png

ticksの部分は、

color1=[0,0,0,0.0]
cb1.ax.tick_params(axis='both', colors=color1)

label ラベルも消す。

20180304150959.png

カラーバーの場合、普通のaxと違って、余計に、axと入れないとダメ。

color1=[0,0,0,0.0]
cb1.ax.yaxis.label.set_color(color1)

二段階必要な、カラーバーの着色部分

20180304151635.png

色の部分を消したいだけなら下記でできるが、白いのが残る。カラーバーを書いてあるaxも、背景は白らしい。
alphaは、0.0〜1.0までの数値で選択。
色の指定の仕方[R,G,B,A]とは違うので注意。

cb1.solids.set_alpha(0)

20180304152157.png

背景の白いのは、下記コードで透明にできます。

cb1.patch.set_alpha(0)#layer under cbar.

カラーバーじゃない方の、axを透明化。

これらは調べればすぐにわかるけど…

spinesを消す。

20180304155248.png
color1=[0,0,0,0.0]
ax1.spines['bottom'].set_color(color1)
ax1.spines['top'].set_color(color1)
ax1.spines['left'].set_color(color1)
ax1.spines['right'].set_color(color1)

軸を消す。ticklabels(数字)も消える。

20180304155500.png
color1=[0,0,0,0.0]
ax1.tick_params(axis='both', colors=color1)

labelも消す。

20180304155619.png

x側は、

color1=[0,0,0,0.0]
ax1.xaxis.label.set_color(color1)

20180304155753.png

y側は、

color1=[0,0,0,0.0]
ax1.yaxis.label.set_color(color1)

imshowで描いたimageを消す

20180304155930.png
imgcol1.set_alpha(0.0)

20180304160002.png
ax1.patch.set_alpha(0.0)
posted by yuchan at 07:00 | Comment(0) | python