Python Matplotlibで円や図形を描くクラスはありますが、今回は数式を用いて図形を表示してみました。
最終的には頂点座標だけではなく、任意の分解能で各辺を分割しつつ、外形線上の任意の点座標を求めます。
円の描画（基本）：
まずは基本として円の描画。図形クラスを使わずに数式で描画するには、
import numpy as np import matplotlib.pyplot as plt theta = np.linspace(0, np.pi*2, 100) r = 1 x = r * np.cos(theta) y = r * np.sin(theta) plt.figure(figsize=(4, 4)) plt.axis('equal') plt.grid() plt.plot(x, y)となります。thetaを0から2πまでの範囲で(x, y)をプロットすると以下のグラフ。
np.linspace(0, np.pi*2, 100)なので、厳密には正99角形を描いていることになります。正100角形にするには、np.linspace(0, np.pi*2, 101)。つまり正n角形を描くにはnp.linspace(0, np.pi*2, n+1)にします。
また、thetaを横軸、yを縦軸にすれば以下のような正弦波になります。
n角形の頂点数による表示：
以下は正方形（n=4）の場合。plt.scatter(x, y)を使うと正方形なら4つの頂点だけがプロットされます。
n = 4 theta = np.linspace(0, 2*np.pi, n+1) r = 1 x = r * np.sin(theta) y = r * np.cos(theta) plt.figure(figsize=(4, 4)) plt.axis('equal') plt.grid() plt.scatter(x, y)
コード内ではxとyをsinとcosで求めていますが（図形の向きを合わせるためxにsinを用いています）、複素数を使うと以下のようになります。虚数は通常iですがPythonでは1jになります。xとyはvの実部realと虚部imagに対応しています。
v = np.exp(1j * theta)
x = v.real
y = v.imag
正方形なので4点ありますが、コード上では以下のplt.plot()で線を引くために座標(1, 0)には2点重なっており合計で5点あります。
そしてplt.plot(x, y)で各頂点を線でつなくと以下のような正方形が描画されます。
この場合、連続値で線を描いているわけではないので、辺上の任意の点座標を数式によって求めていません。連続値を用いて角を持った図形やグラフを描くにはどうすればよいか？ということについて以下に続きます。
正弦波を矩形波へ変換：
正弦波は滑らかな曲線ですが、正弦波を無数に重ね合わせる（無限フーリエ級数）と矩形波に近づいていきます。
数式的には、nを正の奇数として（n=1, 3, 5 ...）、
y = 4/π * {1/1*sin(1*θ) + 1/3*sin(3*θ) + 1/5*sin(5*θ) + ... + 1/n*sin(n*θ)}
となります。
n=10までを重ね合わせるコードは以下
。分解能は100。
theta = np.linspace(0, np.pi*2, 101) n = 10 sq = 4 / np.pi * sum([np.sin(n*theta) / n for n in range(1, n, 2)])
ちなみに、n=10とn=50までの重ね合わせをグラフにすると、
青が元の正弦波、黄色がn=10、緑がn=50のときの矩形波。まだ矩形波としては不完全。
n=200まで重ね合わせるとかなり角がでてきます。無限級数なので数が大きいほどより正確な矩形波に近づいていきます。
上図は横軸がθ（角度）、縦軸がyの値なので、これを元に(x, y)のグラフを表示させると以下のような正方形になります。
今回はグラフの分解能を100にしていますが、分解能を上げるほどより近似していきます。
もう一つの方法：
フーリエ級数の式はやや面倒なので、もう少し簡単にコーディングすると以下になります。
sq = np.sin(theta) for _ in range(10): sq += np.sin(sq) sq /= np.pi
正弦波の式に正弦波を自己代入して合算し（今回は10回自己代入）最後にπで割っています。結果のグラフは以下。
三角波：
次は三角波です。数式的には矩形波のように奇数の級数がでてきますが、一つおきにプラスマイナス反転するので少し複雑です。
tr = 8/π2 * sin{1/12*sin(1*θ) - 1/32*sin(3*θ) + 1/52*sin(5*θ) - 1/72*sin(7*θ) + ...}
上式のsin{...}内の偶数番目（22など）の項をゼロにして奇数番目を一つおきにプラスマイナス反転させるには、
sin(n*π/2)
sin(n*π/2)
を使うと、n=1, 2, 3, 4, 5 ... のとき、
1, 0, -1, 0, 1, 0, -1, 0, 1, ...
という値を返します。この式を使ってコーディングすると（n=100まで）、
1, 0, -1, 0, 1, 0, -1, 0, 1, ...
という値を返します。この式を使ってコーディングすると（n=100まで）、
となり、これを描画すると以下。tr = 8 / np.pi**2 * sum([1/k**2 * np.sin(k*np.pi/2) * np.sin(k*theta) for k in range(1, 101)])
また(x, y)をグラフ表示すると以下。
トロコイド：
車輪を地面上で転がしたときに、車輪上の任意の点が動く軌跡。
車輪の半径をR、車輪上の任意の点の半径をr、回転角をthetaとしたとき、
theta = np.linspace(0, np.pi*2, 100) R = 1 r = 1 x = R*theta - r*np.sin(theta) y = R - r*np.cos(theta)
R=rのときサイクロイドと呼び、(x, y)の軌跡は以下。
横軸：x、縦軸：y
外トロコイド（Epitrochoid）：
上記トロコイドにおいて、地面の代わりに円の外側を回る車輪上の任意の点の軌跡。
wiki：「トロコイド」より
地面の代わりとなる円の半径をR、車輪の半径をr、車輪上の任意の点の半径をn、回転角をthetaとしたとき、
r=nのとき外サイクロイドと呼び、(X, Y)の軌跡は以下。R = 1 r = 1/4 n = 1/4 theta = np.linspace(-np.pi, np.pi, 100) x = R * np.cos(theta) y = R * np.sin(theta) X = (R + r) * np.cos(theta) - n * np.cos((R + r) / r * theta) Y = (R + r) * np.sin(theta) - n * np.sin((R + r) / r * theta)
内トロコイド（Hypotrochoid）：
外トロコイドとは逆に固定円（青い円）の内側を回るときの軌跡。
R = 1 r = 1 / 4 n = 1 / 4 theta = np.linspace(0, np.pi*2, 100) x = R * np.cos(theta) y = R * np.sin(theta) X = (R - r) * np.cos(theta) + n * np.cos((R - r) / r * theta) Y = (R - r) * np.sin(theta) - n * np.sin((R - r) / r * theta)
r=nのとき内サイクロイドと呼び、(X, Y)の軌跡は以下。
青が半径R=1の固定された円、黄色がその内側を転がる半径r=1/4でn=rの軌跡。r=1/3で3分割、r=1/5で5分割の軌跡になる。n=0のとき黄色い軌跡は半径R-rの円になる。
r=1/4、n=1/10にすると以下。
角が丸くならないようにする：
「いかに丸くならない角と直線で多角形を描くことができるか？」という方法についてです。
先ほど、正弦波を無限フーリエ級数的に重ね合わせると直角と直線に近似した矩形波を描くことができましたが、トロコイドでは角をシャープにすれば辺が曲線になってしまい、辺を直線に近づければ角が丸まってしまいます。そこで、もう少しシンプルに円の公式を利用して直角に近似する方法を試してみました。
円の公式：x**2 + y**2 = r**2なので、y = (r**2 - x**2)**(1/2)として、円の半径r=1、xの分解能を100とすると、
n = 2 r = 1 x = np.linspace(-r, r, 101) y = (r**n - x**n)**(1/n)xに対するyが求まり、これを表示すると、
さらにn=4にするとx**4 + y**4 = r**4となり、y = (r**4 - x**4)**(1/4)の場合は、
xn + yn = rnのnを大きくするほど角がでてきます。
さらにn=100にすると、以下のように直角に近づいていきます。
この場合、xが-1から1までの間を分解能100で描いているので、より連続値に近い表現方法になりますが、plt.scatter(x, y)で点描してみると以下のようになってしまいます。
中心角θを媒介変数として多角形を描く：
多角形の頂点を求めて辺をつないで描くのではなく、任意の中心角θにおける多角形の外形線の(x, y)座標を求める方法についてです（こちらを参考にしました）。
以下がコードです。n=3のときに正三角形を描きます（正n角形）。
n = 3 resolution = n * 10 theta = np.linspace(0, np.pi*2, resolution+1) r = 1 / np.cos(theta % (2*np.pi/n) - np.pi/n) x = r * np.sin(theta) y = r * np.cos(theta) plt.figure(figsize=(4, 4)) plt.axis('equal') plt.plot(x, y, alpha=0.5) plt.scatter(x, y, s=10, c='r') for i in range(len(theta)): plt.plot([0, x[i]], [0, y[i]], c='gray', lw=0.5) plt.figure(figsize=(4*np.pi, 4)) plt.grid() plt.axis('equal') plt.plot(theta, r)頂点座標だけではなく、任意の中心角θに対する辺上の点座標(x, y)を求めることができます。
resolutionは中心角0〜2πにおける分解能です。頂点を含めるためnの倍数にしています。rは角度θにおける中心(0, 0)からの距離です。rとθが分かればsinとcosで(x, y)の座標を求めることができます（図形が90度ずれるので、xにsin、yにcosを使用しています）。
n=3、分解能30（1辺を10分割）の場合。赤点が分解能30における辺上の点。
この方法では中心角θを分解能に応じて均等に分割しているためプロットされる点は辺を均等に分割してはいませんが、n角形の頂点だけではなく辺上の点座標を求めることができます（辺を均等に分割する場合は以下のEpicycleの方法で）。
多角形を歪ませて表示：
前述の方法で多角形における頂点以外の任意の点座標も求めることができたので、多角形そのものを歪ませてみます。
正方形を例として、回転、伸縮、湾曲させてみます。
n = 4 resolution = n * 10 theta = np.linspace(0, np.pi*2, resolution+1) r = 1 / np.cos(theta % (2*np.pi/n) - np.pi/n) theta += np.pi/n # tilted at 45 degrees x = r * np.sin(theta) y = r * np.cos(theta) #plt.figure(figsize=(4, 4)) plt.axis('equal') plt.plot(x, y, alpha=0.5) plt.scatter(x, y, s=10, c='r') plt.plot(x*1.2, y) # yellow line y1 = 1/4 * x**2 plt.plot(x+0.3, y+y1+0.1) # green line歪ませた結果は以下。
黄色は横幅を1.2倍したもの。
緑はy=1/4*x**2の放物線を加えて座標(0.3, 0.1)ずらしたもの。
頂点以外の点もプロットされているので全体的に湾曲します。
Epicycle：フーリエ級数で正多角形を描く：
もう一つの方法として、Epicycleという複数の回転軸を連結させて図形を描く方法。
矩形波で複数の波形を重ねたように、周期性のある回転軸をフーリエ級数的に足し合わせて、単純な円から複雑な図形に変換していきます。
基本の数式は、
eiθ= cosθ + i*sinθ
eはエクスポネンシャル、iは虚数、θは任意の角度のとき、複素平面上で実部が横軸、虚部が縦軸となり、
eiπ = -1：座標(1, 0)から反時計回りにπ（180度）回転し座標(-1, 0)に移動
e2iπ = -1：座標(1, 0)から反時計回りにπ（360度）回転し座標(1, 0)に移動
eiπ/2 = i：座標(1, 0)から反時計回りにπ/2（90度）回転し座標(0, i)に移動
e-iπ = -1：座標(1, 0)から時計回りにπ（180度）回転し座標(-1, 0)に移動
2eiπ = -2：座標(2, 0)から反時計回りにπ（180度）回転し座標(-2, 0)に移動
という関係になります。つまり複数の回転軸を持つ円を重ね合わせる際に、円の大きさや回転角、回転する向きを各係数で指定し任意の座標へ移動できます。
例えば、正三角形（n=3）を描画する場合は、係数をcとして、Σは-∞から∞までの範囲で、
R = Σ(1/c2 * eciθ)
となるようです。
求められるRは複素数で、Rの実部をx成分、Rの虚部をy成分とした(x, y)座標に点がプロットされます。
また、eの係数1/c2と指数部分にcがあり、係数cの配列をCとし、
C = [..., -8, -5, -2, 1, 4, 7, 10, ...]
1を基準として3ずつ加算した数列、また3ずつ減算した数列となります。
また正方形（n=4）であれば、
C = [..., -11, -7, -3, 1, 5, 9, 13, ...]
1を基準として4ずつ加算した数列、また4ずつ減算した数列となります。
正五角形（n=5）であれば同様に、
C = [..., -14, -9, -4, 1, 6, 11, 16, ...]
というように係数Cの配列は、1を基準としてnを加算、nを減算したものとなるようです。
仮に級数の個数をF=30としてコーディングすると、
F = 30 C = [1 + f*n for f in range(-F//2, F//2)]という感じになります。nは任意の正n角形。
無限級数なので円の個数Fを多くするほど、より厳密な図形を描くことができます。
以下は全体のコード。
n = 4 resolution = n * 10 theta = np.linspace(0, 2*np.pi, resolution+1) F = 4 C = [1 + f*n for f in range(-F//2, F//2)] R = sum([1/c**2 * np.exp(c*1j*theta) for c in C]) ratio = sum([1/c**2 for c in C]) x = R.real / ratio y = R.imag / ratio plt.figure(figsize=(4, 4)) plt.grid() plt.axis('equal') plt.plot(x, y, alpha=0.5, ) plt.scatter(x, y, s=10, c='r')n：正n角形
resolution：分解能（辺の分割数）
theta：分解能で分割した角度
F：フーリエ級数の数（数が大きいをほど描画精度が上がる）
C：フーリエ級数の係数
R：Epicycleの計算式、複素平面上のベクトル
x：ベクトルの実部をx成分に変換
y：ベクトルの虚部をy成分に変換
ratio：多角形頂点までの半径を1にするための比率
Fを下げると（F=4）描画精度も下がり角が丸くなってしまいます。
ということで、多角形の頂点だけでなく途中の点も各変数を調節しながら求めることができます。
もう一つの方法（分解能で分割）：
フーリエ級数を使わず、多角形の各辺を任意の分解能で分割する場合。
N = 4 theta = np.linspace(0, 2*np.pi, N+1) V = np.exp(1j * theta) # V = np.cos(theta) + 1j*np.sin(theta) Div = 10 P = np.array([(V[k+1] - V[k])*d/Div + V[k] for d in range(Div) for k in range(N)]) x = V.real y = V.imag vx = P.real vy = P.imag plt.figure(figsize=(4, 4)) plt.grid() plt.axis('equal') plt.plot(x, y, alpha=0.5) plt.scatter(vx, vy, s=10, c='r') plt.scatter(x, y, s=20, c='k') plt.scatter(0, 0, s=50, c='k', marker='x')
Divは各辺を何分割するかという変数。その分割数から描く座標を求めています。
Vは正n角形のn個の点、Divは1辺の分割数、PがDivによって分割された点。
任意の中心角による辺上の点座標：
前述のような分割数で求めるのではなく、中心角を入力として辺上の点座標を求める方法。
def intersection(x1,y1,x2,y2,x3,y3,x4,y4): A = (x1-x2)*(y3-y4)-(y1-y2)*(x3-x4) B = (x1-x2)*(y3-y4)-(y1-y2)*(x3-x4) if A != 0 and B!= 0: px= ( (x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4) ) / A py= ( (x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4) ) / B return [px, py] N = 4 theta = np.linspace(0, 2*np.pi, N+1) R = np.exp(1j * theta) x = R.real y = R.imag deg = np.deg2rad(70) n = int(deg * N / (2*np.pi)) % N px, py = intersection(0, 0, np.cos(deg), np.sin(deg), R[n].real, R[n].imag, R[n+1].real, R[n+1].imag) plt.figure(figsize=(4, 4)) plt.grid() plt.axis('equal') plt.plot(x, y, alpha=0.5) plt.plot([0, px], [0, py]) plt.scatter(px, py, s=10, c='r')intersection()という二つの線分の交点を求める関数を用いています。中心からの半径となる線分と多角形の各辺との交点を求めるという手順です。
変数Nで正N角形を指定し、変数degに任意の角度（度数）を指定します。変数nは中心角radが正N角形のn番目の辺にあるという計算をしています。
np.deg2rad()は度からラジアンへ変換してくれます。
0 件のコメント:
コメントを投稿