１．はじめに
前回の記事（ 【アナログ プラニメータ】面積の測定原理を自己流に証明してみた ）では、ポーラープラニメータ（Polar planimeter）の動作原理を証明し、これがリニアプラニメータ（Linear planimeter）の証明をも兼ねることを述べました。ポーラープラニメータはその関節部が円弧上に、リニアプラニメータでは直線上に拘束されていました。しかし、拘束線の形状として許されるのは、この２種類だけなのでしょうか？
今回の記事では、実益については特に顧みず、他の形状の拘束線でも面積の測定が可能なことをシミュレーションで確認し、その裏付けとして、任意の拘束線をも含む包括的な面積測定原理について再証明してみたいと思います。最後に、この証明に利用した MATLAB の数式処理ツール Symbolic Math Toolbox についての雑感を述べます。
２．種々の拘束線でのシミュレーション
疑問が湧いたときには、理屈はさておき「まずは数値計算シミュレーションで様子を見るのが早道」 ということで、その結果を次に示しています。図１は、前回の記事と同一条件でのポーラープラニメータの動作、図２は、その特殊ケースとしてのリニアプラニメータの動作を示すアニメーションです。青い太線のアームの一端Ｊは、その移動が赤色の拘束線上だけに制限されています。このアームの他端Ｔを持ち、青色の被測定図形の外周上をなぞれば、緑色の目盛車の積分動作により面積が求められます。前回のように目盛車の回転まで再現するのは手間がかかるため、ここでは測定値の推移をデジタル表示することで代用しています。当然のことながら、これら実在のタイプのプラニメータでは、最終的な値は正しく 230.5 [${\rm cm}^2$] を表示しています。
それでは、この記事で注目している架空の珍種タイプではどうなるのでしょうか。図３は、拘束線が「く」の字に折れ曲がった直線、図４は放物線、図５は正弦波曲線の場合のアニメーションです。どの図でも、面積の最終測定結果は 230.5 [${\rm cm}^2$] と正しく表示されており、目論見どおりです。これらの結果により、どのような形状の拘束線でも、正確な面積が測定できそうな気配が濃厚になってきました。図６は、これら５種類の拘束線について、計測途中の面積の表示値をグラフとしてまとめたものです。計測途中での値はバラバラですが、一周した後には綺麗に正しい値に収束しています。
ところで、この図を見ると、収束点付近のグラフの傾斜が大きいケースでは、停止位置の僅かなズレによって、測定精度が低下しそうな印象を受けます。しかし、この傾斜は、「被測定図形の配置」や「計測開始点の位置」によって変わるもので、拘束線の形状だけで決定づけられるわけではありません。
さて、あとは、面積測定の理論的な証明をするだけです。
図１ 定番のポーラープラニメータの動作
図２ 定番のリニアプラニメータの動作
図３ 珍種の くの字プラニメータの動作
図４ 珍種の放物線プラニメータの動作
図５ 珍種の正弦波プラニメータの動作
図６ 拘束線の種類による特性の違い
３．面積測定が可能な領域
理論的な証明に移る前に、幾つか確認しておきたいことがあります。前記事では、ポーラープラニメータの場合、次のような使用条件があることを述べました。
- 被測定図形を置くことができる領域は、２つの同心円に囲まれた範囲内である。
- この同心円の線上に、被測定図形の一端でも掛かると、ここでＴ点の移動が阻止される。
- ２本のアームが作る形状には「＜」と「＞」の２種類が存在し、計測中は、最初に選んだ形状を保ち続ける必要がある。
ところが、これらの条件は、拘束線が円弧であることを前提とした表現になっています。今回のように拘束線の形が異なる場合には、そのままの表現では通用しません。そこで、新たに次のような用語を定義して使用することにします。
アーム配置
「＜」や「＞」についての、他の型式にも通じる言い換えです。まずは、拘束線に沿った任意の向きに正方向を定義します。そして、「拘束線の接線の角度」と「アームの角度」との相対関係により次のように区分します。Ｔ端からアームを引き寄せたとき、Ｊ点が拘束線の正方向へ進むようなアーム配置を引きアーム、逆に、Ｔ端からアームを押し遣ったときに正方向へ進むような配置を押しアームとします。図１～５で言えば、拘束線の右下から左上に進む向きを正方向に定めて、すべての図で押しアームの配置になるように統一しています。
アクセス可能領域
アームのＪ端が、与えられた拘束線の全亘長を移動したとき、「押し」「引き」いずれかのアーム配置において、Ｔ端が動くことができる全領域を指します。拘束線上の任意のＪ点において、Ｔ端がアクセスできるのは、Ｊ点を中心にしたアーム長を半径とする半円です。「押しアーム」のときは拘束線の負側の半円、「引きアーム」のときは正側の半円になります。これらの半円が集まって作る領域が「アクセス可能領域」になります。
侵入禁止線
ポーラープラニメータの場合の「同心円の線上」に相当するものです。アーム配置が「押し」でも「引き」でもない状態（アームが拘束線と直交している状態）のとき、Ｔ点はこの線上に乗ります。言い換えれば、Ｊ点が、「アームと拘束線との直交状態」を保ちながら移動したとき、Ｔ点が作る軌跡です。「アクセス可能領域」内であっても、被測定図形がこの「侵入禁止線」と重なるときには注意が必要です。
以上を踏まえ、図７～１１に、各タイプのプラニメータの動作領域を示しています。それぞれ図１～５のアニメーションと同一条件ですが、表示範囲を広くして全貌が分かるようにしています。ここで採用した拘束線は、与えられた被計測図形の全体をカバーできる適度な長さにしています。しかし、計測時に実際に使用されるのは、それよりやや短い、赤い太線部分です。多数の半円によって「ヘアライン加工」されたような領域がアクセス可能範囲、ピンク色の線が侵入禁止線です。アクセス可能領域を縁取る線のうち、ピンク色以外の線は「外部への食み出し禁止線」ではありますが、Ｔ点が乗っても移動が阻止されることはありません。
図７ ポーラープラニメータの動作領域
図８ リニアプラニメータの動作領域
図９ くの字プラニメータの動作領域
図１０ 放物線プラニメータの動作領域
図１１ 正弦波プラニメータの動作領域
図７はポーラープラニメータの動作領域です。拘束線の長さを短くしているので、侵入禁止線は完全な同心円ではなく、その一部だけになっています。図８のリニアプラニメータの侵入禁止線も想定の範囲内であり、どちらにも目立った特徴はありません。ところが、図９～１１の珍種では、侵入禁止線がアクセス可能領域内にまで入り込んでいます。また、この付近には、２種類のヘアライン模様が重なってできたメッシュ状の領域も存在します。これらは、拘束線の曲率半径が一定ではなく、しかも、アーム長より短い部分が存在するときに起こる現象です。とりあえず、細かいことは気にしないのであれば、被測定図形をこれらから離して配置しておけば大丈夫です。
細かいことでも気になる方のために少しだけ補足しておきます。珍種の動作が定番プラニメータと異なるのは、侵入禁止線に乗っても、必ずしもＴ点が身動きできなくなるわけではないことです、例えば、Ｔ点を「ヘアラインの縦糸の領域」から移動し始めたとき、その縦糸領域が続く限り、横糸が作る侵入禁止線は自由に横切って移動することができます。しかし、縦糸が作る侵入禁止線に触れると、そこで動けなくなります。いずれにしても、これは理論的に難しく考えるような問題ではなく、実機を作って操作してみた場合、機械的に引っかかって動けなくなるかどうかで体感できる現象です。
４．任意の拘束線での測定原理の証明
図１２は、今回の証明のためのモデルです。前記事ではアーム１が作る円弧だった拘束線を、弧長パラメータ $s$ で表された下記の任意の曲線で置き換えています。
z_{\rm j}=x_{\rm j}+iy_{\rm j}=f(s)+ig(s)\tag{1}
ここで、$f$, $g$ は $s$ を媒介変数とする任意の関数です。
図１２ 検討対象のモデル
また、今回は、「アーム１」に相当するものがなく識別は不要になったため、前記事で $r_2$, $\theta_2$ と表現していた部分は $r$, $\theta$ と簡略化しています。なお、これ以降は、前回の記事の「試行３」との重複が多いため、詳細な説明は省いている部分もあります。前記事と照らし合わせながら読んでいただくと分かりやすいと思います。
図１２と(1)式より、Ｔ点とＷ点の ${\rm x}$, ${\rm y}$ 座標上の位置は、それぞれ次のようになります。
z_\mathrm{t}\,(=x_{\rm t}+iy_{\rm t})=z_\mathrm{j}+r e^{i\theta}=f(s)+ig(s)+r e^{i\theta}\tag{2}
z_\mathrm{w}\,(=x_{\rm w}+iy_{\rm w})=z_\mathrm{j}+a e^{i\theta}=f(s)+ig(s)+a e^{i\theta}\tag{3}
また、Ｗ点を ${\rm u}$, ${\rm v}$ 座標上から見た位置は、次式のとおりです。
z_\mathrm{w'}\,(=u+iv)=z_\mathrm{w}e^{i(\frac{\pi}{2}-\theta)}\tag{4}
これらの式を一般の複素座標に変換し、実数部と虚数部に分けたうえ、微小変位量どうしを関係づける形に変形していきます。まずは、(2)式より、
x_{\rm t}+iy_{\rm t}=f(s)+r\cos\theta+i\,\{g(s)+r\sin\theta\}\tag{5}
\begin{align}
\left\{
\begin{array}{l}
x_{\rm t}=f(s)+r\cos\theta\\
y_{\rm t}\,=g(s)+r\sin\theta
\end{array}
\right.\tag{6}
\end{align}
\begin{align}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]
&=
\left[
\begin{array}{c}
\partial x_{\rm t}/\partial s & \partial x_{\rm t}/\partial\theta\\
\partial y_{\rm t}/\partial s & \partial y_{\rm t}/\partial\theta
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}s\\
{\rm d}\theta
\end{array}
\right]\\[1.5ex]
&=
\left[
\begin{array}{c}
f' & {-}r\sin\theta\\
g' & r\cos\theta
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}s\\
{\rm d}\theta
\end{array}
\right]
=
\boldsymbol{A}
\left[
\begin{array}{c}
{\rm d}s\\
{\rm d}\theta
\end{array}
\right]\tag{7}
\end{align}
ただし、$f'={\rm d}f(s)/{\rm d}s$、$g'={\rm d}g(s)/{\rm d}s$ とします。
あとの説明で、(7)式の逆変換も必要となるので、そのための式を求めておくと、
\begin{align}
\left[
\begin{array}{c}
{\rm d}s\\
{\rm d}\theta
\end{array}
\right]
&=
\boldsymbol{A}^{-1}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]
=
\boldsymbol{B}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]\\[1.5ex]
&=
\left[
\begin{array}{c}
\partial s/\partial x_{\rm t} & \partial s/\partial y_{\rm t}\\
\partial\theta/\partial x_{\rm t} & \partial\theta/\partial y_{\rm t}
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]\tag{8}
\end{align}
ここで、
\begin{align}
\boldsymbol{B}=\boldsymbol{A}^{-1}=
\frac{1}{f'\cos\theta+g'\sin\theta}
\left[
\begin{array}{c}
\cos\theta & \sin\theta\\
-g'/r & f'/r
\end{array}
\right]\tag{9}
\end{align}
つぎに、(3)式の $z_{\rm w}$（${=}x_{\rm w}{+}iy_{\rm w}$）では、(2)式の $r$ が $a$ に変わっているだけなので、(7)式を参照して次の式が得られます。
\begin{align}
\left[
\begin{array}{c}
{\rm d}x_{\rm w}\\
{\rm d}y_{\rm w}
\end{array}
\right]
=
\left[
\begin{array}{c}
f' & {-}a\sin\theta\\
g' & a\cos\theta
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}s\\
{\rm d}\theta
\end{array}
\right]
=
\boldsymbol{C}
\left[
\begin{array}{c}
{\rm d}s\\
{\rm d}\theta
\end{array}
\right]\tag{10}
\end{align}
${\rm x}$, ${\rm y}$ 座標から ${\rm u}$, ${\rm v}$ 座標へ変換するための(4)式は次のように変形できます。
\begin{align}
u+iv&=z_\mathrm{w}e^{i(\frac{\pi}{2}-\theta)}\\
&=
(x_{\rm w}+iy_{\rm w})i(\cos\theta-i\sin\theta)\\
&=
(x_{\rm w}+iy_{\rm w})(\sin\theta+i\cos\theta)\\
&=
x_{\rm w}\sin\theta-y_{\rm w}\cos\theta+i(x_{\rm w}\cos\theta+y_{\rm w}\sin\theta)
\tag{11}
\end{align}
\begin{align}
\left\{
\begin{array}{l}
u=x_{\rm w}\sin\theta-y_{\rm w}\cos\theta\\
v=x_{\rm w}\cos\theta+y_{\rm w}\sin\theta
\end{array}
\right.\tag{12}
\end{align}
\begin{align}
\left[
\begin{array}{c}
u\\
v
\end{array}
\right]
=
\left[
\begin{array}{c}
\sin\theta & {-}\cos\theta\\
\cos\theta & \sin\theta
\end{array}
\right]
\left[
\begin{array}{c}
x_{\rm w}\\
y_{\rm w}
\end{array}
\right]
=
\boldsymbol{D}
\left[
\begin{array}{c}
x_{\rm w}\\
y_{\rm w}
\end{array}
\right]
\tag{13}
\end{align}
$\theta_2$ が $\theta$ に変わってはいますが、この $\boldsymbol{D}$ は前回の記事の変換行列と同じものです。そして、微小変位に対しても同式が適用できるので、
\begin{align}
\left[
\begin{array}{c}
{\rm d}u\\
{\rm d}v
\end{array}
\right]
=
\boldsymbol{D}
\left[
\begin{array}{c}
{\rm d}x_{\rm w}\\
{\rm d}y_{\rm w}
\end{array}
\right]
\tag{14}
\end{align}
以上、(8),(10),(14)式をまとめると、Ｔ点の $\rm x,y$ 座標上の微小変位に対する、Ｗ点の $\rm u,v$ 座標上の微小変位の関係式が得られます。
\begin{align}
\left[
\begin{array}{c}
{\rm d}u\\
{\rm d}v
\end{array}
\right]
=
\boldsymbol{D}
\left[
\begin{array}{c}
{\rm d}x_{\rm w}\\
{\rm d}y_{\rm w}
\end{array}
\right]
=
\boldsymbol{D}
\boldsymbol{C}
\left[
\begin{array}{c}
{\rm d}s\\
{\rm d}\theta
\end{array}
\right]
=
\boldsymbol{D}
\boldsymbol{C}
\boldsymbol{A}^{-1}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]\tag{15}
\end{align}
$\boldsymbol{D}\boldsymbol{C}$ に(10)式, (13)式を代入して、
\begin{align}
\boldsymbol{D}\boldsymbol{C}&=
\left[
\begin{array}{c}
\sin\theta & -\cos\theta\\
\cos\theta & \sin\theta
\end{array}
\right]
\left[
\begin{array}{c}
f' & {-}a\sin\theta\\
g' & a\cos\theta
\end{array}
\right]\\[2ex]
&=
\left[
\begin{array}{c}
f'\sin\theta{-}g'\cos\theta & -a\sin^2\theta{-}a\cos^2\theta\\
f'\cos\theta{+}g'\sin\theta & -a\sin\theta\cos\theta{+}a\sin\theta\cos\theta
\end{array}
\right]\\[2ex]
&=
\left[
\begin{array}{c}
f'\sin\theta{-}g'\cos\theta & -a\\
f'\cos\theta{+}g'\sin\theta & 0
\end{array}
\right]\tag{16}
\end{align}
さらに、(16)式に(9)式を掛けて、
\begin{align}
\boldsymbol{DCA}^{-1}&=
\frac{1}{f'\cos\theta{+}g'\sin\theta}
\left[
\begin{array}{c}
f'\sin\theta{-}g'\cos\theta & -a\\
f'\cos\theta{+}g'\sin\theta & 0
\end{array}
\right]
\left[
\begin{array}{c}
\cos\theta & \sin\theta\\
-g'/r & f'/r
\end{array}
\right]\\[2ex]
&=
\frac{1}{f'\cos\theta{+}g'\sin\theta}*\\[1.5ex]
&\hspace{2em}*
\left[
\begin{array}{c}
f'\sin\theta\cos\theta{-}g'\cos^2\theta{+}(a/r)g' & f'\sin^2\theta{-}g'\sin\theta\cos\theta{-}(a/r)f'\\
f'\cos^2\theta{+}g'\sin\theta\cos\theta & f'\sin\theta\cos\theta{+}g'\sin^2\theta
\end{array}
\right]\tag{17}
\end{align}
ここで、面積の測定値には寄与しない ${\rm d}v$ は無視して ${\rm d}u$ のみに着目すると、(15)式のうち、この ${\rm d}u$ の計算に必要になるのは行列 $\boldsymbol{DCA}^{-1}$ の1行1列、1行2列の要素だけになります。そこで、これらを次のように $H_{\rm x}$ , $H_{\rm y}$ とおきます。
H_{\rm x}=\frac{f'\sin\theta\cos\theta{-}g'\cos^2\theta{+}(a/r)g'}{f'\cos\theta{+}g'\sin\theta}\tag{18}
H_{\rm y}=\frac{f'\sin^2\theta{-}g'\sin\theta\cos\theta{-}(a/r)f'}{f'\cos\theta{+}g'\sin\theta}\tag{19}
なお、$f'\cos\theta{+}g'\sin\theta=0$ は拘束線とアームが直交する条件になります。したがって、侵入禁止線を避けて使用していれば、上式の分母が 0 になることはありません。
(18)(19)式を用いて、(15)式から必要部分だけを書き出すと、
{\rm d}u=
(\boldsymbol{D}\boldsymbol{C}\boldsymbol{A}^{-1})_{(1行目)}
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]
=
\left[
\begin{array}{c}
H_{\rm x} & \hspace{-0.7em}H_{\rm y}
\end{array}
\right]
\left[
\begin{array}{c}
{\rm d}x_{\rm t}\\
{\rm d}y_{\rm t}
\end{array}
\right]
=
\boldsymbol{H}{\cdot}{\rm d}\boldsymbol{s}
\tag{20}
ここで、$\boldsymbol{H}{=}[H_{\rm x}\ H_{\rm y}]$ はプラニメータが作るベクトル場、${\rm d}\boldsymbol{s}{=}[{\rm d}x_{\rm t}\ {\rm d}y_{\rm t}]$ はＴ点が図形の輪郭線をなぞるときの微小変位ベクトルです（先に出てきた弧長パラメータの $s$ とは関係ありません）。
これ以降は、前回の記事との重複がさらに増えるので、大幅に省略して重要な部分のみを示すと次のようになります。
任意の形状の拘束線の場合でも、次式が成立することを証明できれば、円弧や直線のときと同様にして面積を求めることができます。
\frac{\partial H_{\rm y}}{\partial x}-\frac{\partial H_{\rm x}}{\partial y}=-\frac{1}{r}\tag{21}
しかし、$H_{\rm x}$ $H_{\rm y}$ は $x$ や $y$ の直接の関数としては表現されていないので、次のように変形してから計算します。
\frac{\partial H_{\rm y}}{\partial x}-\frac{\partial H_{\rm x}}{\partial y}
=
\left(\frac{\partial H_{\rm y}}{\partial s}\frac{\partial s}{\partial x}
{+}\frac{\partial H_{\rm y}}{\partial\theta}\frac{\partial\theta}{\partial x}\right)
-
\left(\frac{\partial H_{\rm x}}{\partial s}\frac{\partial s}{\partial y}
{+}\frac{\partial H_{\rm x}}{\partial\theta}\frac{\partial\theta}{\partial y}\right)\tag{22}
ここで $x$, $y$ は、$x_{\rm t}$, $y_{\rm t}$ がとり得るすべての値を意味しているので、上式の $\partial s/\partial x$、$\partial s/\partial y$、$\partial\theta/\partial x$、$\partial\theta/\partial y$ は、(8)式の $\boldsymbol{B}$ の各要素と同じものになります。そこで、
\begin{align}
\boldsymbol{B}
=
\left[
\begin{array}{c}
\partial s/\partial x_{\rm t} & \partial s/\partial y_{\rm t}\\
\partial\theta/\partial x_{\rm t} & \partial\theta/\partial y_{\rm t}
\end{array}
\right]
=
\left[
\begin{array}{c}
b_{11} & b_{12}\\
b_{21} & b_{22}
\end{array}
\right]\tag{23}
\end{align}
とおくと、最終的に、
\frac{\partial H_{\rm y}}{\partial x}-\frac{\partial H_{\rm x}}{\partial y}
=
\left(\frac{\partial H_{\rm y}}{\partial s}b_{11}
{+}\frac{\partial H_{\rm y}}{\partial\theta}b_{21}\right)
-
\left(\frac{\partial H_{\rm x}}{\partial s}b_{12}
{+}\frac{\partial H_{\rm x}}{\partial\theta}b_{22}\right)\tag{24}
となります。
これを、MATLAB の Symbolic Math Toolbox を使って下記のように計算すると、答えは $-1/r$ となって、(21)式が成立することが直ちに確認できます。しかし、MATLAB は最終結論しか教えてくれません。計算の詳細な手順については、その後に続く人力解析による手書きの画像を参照ください。なお、$f'^2$ などは、手書きでは $f^{12}$ のように誤認されやすいので、$f'f'$ などと変則的な表記にしています。
以上の結果により、前回の記事でも示したように、被計測図形の面積 $\Delta S$ は次式で求められます。
\Delta S=r\,\Delta u_\circlearrowright=2\,\pi\,r_{\rm w}\,r\,\Delta N_\circlearrowright\tag{25}
ただし、$\Delta u_\circlearrowright$ : Ｔ点を時計回りに１周させたときの Ｗ点の $\rm u$ 方向への移動距離の積算値、$\Delta N_\circlearrowright$ : 同じく、目盛車の回転数。
ところで、前述の「アーム配置」とか「多重のヘアライン領域」などについては、この証明では何も言及していません。しかし、疑問を持つ必要はありません。これらを考慮せずに証明できたということは、「押し / 引き」どちらのアーム配置でも、また、多重のそれぞれ個別の領域でも、常に(21)式が成り立つということを意味しています。
% planimeter15.m
% 任意拘束線プラニメータのrotationの数式計算
clear
close all
%tic
% r: アームの長さ
% a: 目盛車の位置
% th: アームの角度
% s: 任意拘束線上の媒介変数
% f,g: 任意拘束線の形（x=f(s), y=g(s)）
syms r th s f(s) g(s) a
A = [diff(f) -r*sin(th) ; diff(g) r*cos(th)];
B = inv(A);
C = [diff(f) -a*sin(th) ; diff(g) a*cos(th)];
D = [sin(th) -cos(th) ; cos(th) sin(th)];
H = D*C*B;
BB = formula(B); % 行列内の各要素の読み出しを可能にする。
HH = formula(H); % 例えば、B(1,1)はエラーになるが、BB(1,1)ならＯＫ。
Hx = HH(1,1);
Hy = HH(1,2);
Rot_raw = diff(Hy,s)*BB(1,1)+diff(Hy,th)*BB(2,1) ...
-diff(Hx,s)*BB(1,2)-diff(Hx,th)*BB(2,2);
Rot_sim = simplify(Rot_raw)
%toc
>> planimeter15
Rot_sim =
-1/r
５．Symbolic Math Toolbox 雑感
Symbolic Math Toolbox は便利なツールですが、それほど頻繁に使う機会はありません。それだけに、なかなか熟達できず、この記事の作成でも色々と苦労しました。ここでその苦労の一端を述べて、将来の改良への参考にしていただこうと思います。なお、使用したバージョンは Home 版のR2019a です。バージョン違いや未熟さによる間違った記述があればご指摘ください。
現在はコマンド体系も変わり、古いスクリプトはそのままでは動きませんが、遠い昔に使ったことはあります。しかし、その当時の印象は、失礼ながら「面白い玩具で、初歩的な変形には使えそうだが それほど賢くはない」でした。簡単な数式を変形させようとしても、頼みもしない３倍角の公式まで持ち出して来て、かえって面倒な形に変えてくれたり、定石の平方完成の方法も知らなかったりと不満が沢山ありました。今回、久しぶりに使ってみましたが、私の当初の評価は、この不満を引き摺ったままのものでした。
今回扱ったのはやや複雑な式なので、上記のそこそこの評価性能に合わせるために、簡単なステップに分解して、その都度、Toolbox に伺いを立てながら作業を進めていきました。初期の段階で少しでも間違えると、それ以降は水の泡になってしまうので、慎重に進める必要があります。
それでも、毎ステップ、なかなか期待どおりの答えを出してくれないので手数がかかります。最初のうちは、この答え合わせに使った手段が下記の ① だけでした。数式操作の訓練にはなりますが大いに疲れます。しかし、徐々に進歩して、②、③ と変わり、最終的には ④ に落ち着きました。すぐに ④ に気づいても良さそうなものですが、私も「それほど賢くはない」部類なので仕方ありません。
- ① 自分が出した答えと Toolbox の答えが同じであることを確認するために、両者が一致するまで手作業で数式を変形。しかし、これは大変な苦役。
- ② ①は面倒なので、変数に適当な数値を代入し、答えの数値が同じなら式も同じだろうという手抜き。たまたま一致しただけという危険が付きまとう。
- ③ 数式の簡略化コマンドを利用し、
simplify (自分の答え - Toolboxの答え)が 0 になれば、自分の答えは正しいと判断。Toolbox にわざわざ２回もお伺いをたてる必要があって無駄。 - ④
simplify (変形前の式 - 自分の答え)が 0 になれば、自分の答えは正しいと判断。Toolbox へのお伺いは１回だけで済む。
このようなステップを踏んで手作業を続けていても、そのうち、先が見えなくなって気力が萎えてきます。そのときに思いついたのが「無理だとは思うが、Toolbox に一気に変形させてみよう」でした。あまり期待は持たずにやらせてみたところ、前述の簡単なブログラムだけで、瞬時（0.5 秒弱）に答えが返ってきました。しかし、あまりにも投げ遣りな答えです。沢山あったはずのパラメータが殆ど消えてしまっています。未だ最終的な正解が分からないままの段階だったため、明らかに間違いであると判断しました。ここで放った捨て台詞「無理なことは分かっておった。手に負えなければ NaN でも返してくれる方がよっぽどマシじゃ。間違った答えを平然と出してくるとは呆れたものよ」。しかし、よくよく確認してみたところ、これが何と正解のようです。
これには驚きました。「それほど賢くはない」と見下していた相手が、瞬時に畏敬の存在に変わりました。そこで、それならばと、以前に手作業で苦労して解いたことがある４重の定積分をやらせてみました。実のところ、これを軽々と解かれてしまうと、私の立場がなくなってしまうので不安です。ところが、どんなに待っても答えが出てきません。畏敬の念は「ある程度のことならできそう」に変わり、実際の能力を正しく把握できたので一安心しました。一発芸を見せられただけで怯んだりせず、相応にお付き合いさせていただくのが良さそうです。
ところで、一気に正解を出してくれるのは助かりますが、人間にはその途中の変形過程が分からずに困ります。MATLAB には日ごろから信頼を置いていますが、結果を鵜呑みにするのは心配です。やはり、手作業で１ステップごとに確認する仕事は残ります。上記の ④ の作業はこれからも必要になりそうです。しかし、暗中模索ではなく、遠くに明かりが見える状態での作業になるので気は楽です。
将来的には、Toolbox が途中の変形過程まで教えてくれるようになるのでしょうか。しかし、そうなると人間の考えることがなくなり愚かになってしまいそうです。なかなか難しいところです。
なお、今回は使いませんでしたが、この Toolbox には、数式処理のほかに有効数字の桁数を自由に設定できる機能があります。これもなかなか役に立ちます。「 微分を使わずテイラー展開 」や「 実数だけでローラン展開 」の記事では、これを使った面白い応用を紹介しました。
６．おわりに
ポーラーや リニアプラニメータだけで足りてきたうえに、これらが既に過去の遺物として葬り去られている現状で、ここで紹介した珍種を作ってみたところで何のメリットもありません。作るのがややこしいうえに、計測可能なモデルの大きさも制限を受けます。しかし、この記事の作成は、適度な頭の体操になって大いに楽しめました。認知機能低下の予防には良いテーマでした。ここまでお付き合い頂きどうも有難うございました。
さて、図１２では、$s$ を弧長パラメータとして扱い、拘束線に等間隔目盛を振っていますが、結果としてそこまでする必要はなかったようです。証明の中では、弧長パラメータに付きものの ${\rm d}s{=}\sqrt{{\rm d}x^2{+}{\rm d}y^2}$ の関係がどこにも使われていません。この目盛は、単調増加的に適当に振っておくだけで良かったようです。
この証明の過程で、任意曲線を表す関数 $f$, $g$ は二階まで微分されています。生真面目に考えれば、二階微分できない曲線は拘束線として使えないことになります。そうなると、「く」の字は微分不可能なので明らかに失格、放物線や正弦波も、シミュレーションでは折線で代用しているので失格です。それでも問題なく正確にシミュレーションできています。本件はどのように辻褄を合わせれば良いのでしょうか。「実は、尖っているように見えてもR面取りしてあるんですよ」と、お茶を濁すことしかできません。これ以上、あまり深くは考えないことにしておきます。
拘束線がどんな曲線でも良いことは、私が見つけた（一時はそうも思ったんですが…）訳ではなく既知の事実です。Prytzプラニメータ（ここではまだ触れていない無拘束の引き摺り方式。精度を確保するには神業が要求される）の研究によると、拘束線が作る面積が 0 であれば、被測定図形の面積測定のための条件は満たされるようです。今回のように、開放された一本の線を拘束線に採用すれば、面積は当然 0 になるので測定が可能になります。それにしても、新規性もないし役にも立たない珍種を、わざわざ、数式による証明 や シミュレーションまで付けて紹介する物好きはそうは居ないでしょう。それだけに、この記事には希少価値があって「いいね」👍（戒律「不自讃毀他戒」破り😅）。
７．プログラム
本題の Symbolic Math Tolbox のプログラムは既に提示済です。これだけでは少し物足りないので、図１～１１を描くための MATLAB プログラムも添付しておきます。Toolbox 類は使用していません。
■■ 折り畳んであります ■■
% planimeter32.m
% 各種形状の拘束線（Ｊ線）での、プラニメータのシミュレーション。
% どの形状でも、計測結果は同じ値になることの実証。
clear
close all
% 静止画は必ず表示されるが、
% ファイル出力の要否については指定が必要。
% 0: 不要、1: 必要
still_out=0; % ■■■■■
% アニメーションの要否。要であれば種類を指定。
% 0: 不要
% 1: Ｊ線が円弧のアニメ
% 2: 〃 直線のアニメ
% 3: 〃 くの字のアニメ
% 4: 〃 放物線のアニメ
% 5: 〃 正弦波のアニメ
anime=5; % ■■■■■
% アニメーションの gif ファイル出力の要否。
% 0: 不要、1: 必要
anime_out=0; % ■■■■■
% ===============================================
% 画像の出力ファイル名作成のための準備
path=mfilename('fullpath'); % このプログラムのファイル名を取得。
% 拡張子なしのフルパス表現。出力フ
% ァイル名のヘッダーの一部にする。
% 標準線色
ao=[0 0.4470 0.7410]; % '#0072BD'
aka=[0.8500 0.3250 0.0980]; % '#D95319'
daidai=[0.9290 0.6940 0.1250]; % '#EDB120'
murasaki=[0.4940 0.1840 0.5560]; % '#7E2F8E'
midori=[0.4660 0.6740 0.1880]; % '#77AC30'
sorairo=[0.3010 0.7450 0.9330]; % '#4DBEEE'
chairo=[0.6350 0.0780 0.1840]; % '#A2142F'
% 各部の寸法など
ra=167; % アームの長さ
a=16; % 目盛車のＪ点からの位置
rw=19.05/2; % 目盛車の半径
rt=7; % Ｔ端のレンズの半径（適当な値）
n_start=120; % 計測開始点のfree_curve上の要素番号
% 被計測図形の概形（スプライン曲線の骨格にする点のx,y座標値）
curve_path=[
0 0
1 0
0.8 0.6
1 1
0 1
0.1 0.5
0 0
];
% 被計測図形の大きさや位置の調整
kx=150; ky=150; dx=120; dy=-70;
curve_path=[curve_path(:,1)*kx+dx curve_path(:,2)*ky+dy];
% 被計測図形のスプライン曲線化データの取得
[X,Y]=free_curve(curve_path,-10,2); % ローカル関数の呼び出し。
% 指定された要素番号の点が始点になるように要素の並びを回転シフト。
X(end)=[]; % 閉曲線では始点と終点の座標が重複するので、
Y(end)=[]; % 終点を一旦削除。
X=[X([n_start:end]) X([1:n_start-1]) X(n_start)]; % 並べ替えの
Y=[Y([n_start:end]) Y([1:n_start-1]) Y(n_start)]; % シフト処理
% プラニメータの仕様に合わせ、要素の並びが時計回りになるよう反転。
X=fliplr(X);
Y=fliplr(Y);
% 被測定図形を、一定長の短線分の折線として再構成する。
% 各折れ点を、Ｔ点の移動ステップとする。
% アニメーションのＴ点の移動速度を一定にするため。
[xt,yt]=divide_line_into_fixed_short_segments(X,Y,2);
% ローカル関数の呼び出し。
Lt=sqrt(diff(xt).^2+diff(yt).^2); % 被測定図形の各線分の長さ
Lt=cumsum(Lt); % 被測定図形の外周線上において、
% 測定始点から各Ｔ点までの長さ
figure(10); % 被測定図形の画像（コピー用の原図）
plot(xt,yt,'Color',ao,'LineWidth',1.5); % 被計測図形を描く。
hold on
ths=atan2d(yt(2)-yt(end-1),xt(2)-xt(end-1))+90; % 被計測図形上に
plot(xt(1)+[-15 15]*cosd(ths), ... % 計測開始の位置を
yt(1)+[-15 15]*sind(ths), 'Color',ao); % 表示する。
axis equal
axis([-50 300 -100 200]);
grid on
xlabel('[mm]');
ylabel('[mm]');
ff0=get(gcf,'Children'); % 図を ff0 にコピー
% ===== Ｊ線が円弧 =====
figure(1);
hf1=copyobj(ff0,gcf); % 被測定図形をペースト
th=linspace(10,105,101);
% 円弧の半径は 164 [mm] とする
Xj1=164*cosd(th);
Yj1=164*sind(th);
plot(Xj1,Yj1,'LineWidth',0.5,'Color',aka); % Ｊ線を描画
axis([-150 400 -150 350]);
half_circle(Xj1,Yj1,ra,-1,1);
% Ｔ点の可動領域を示す半円群などを描画。
% ローカル関数の呼び出し。
[Xjp,Yjp]=find_j_points(xt,yt,Xj1,Yj1,ra,-1);
% 各Ｔ点に対応するＪ点の座標群を求める。
% ローカル関数の呼び出し。
plot(Xjp,Yjp,'Color',aka,'LineWidth',1.5)
% Ｊ線のうち、利用される範囲を太線化。
draw_arm(xt(1),yt(1),Xjp(1),Yjp(1),0,a,rw,rt,0);
% アームの描画。ローカル関数の呼び出し。
[~,dLw1]=displacement_of_wheel(xt,yt,Xjp,Yjp,ra,a);
% 各ステップでの目盛車の回転方向への
% 移動量。ローカル関数の呼び出し
% 静止画のファイル出力
if still_out==1
% 画像全体のトリミング
figure_trimmer(gcf,[-10 0 -30 -45]); % ローカル関数の呼び出し。
% 画像をファイルに出力
% （本ケースでは、png よりも jpg が小容量。解像度はやや落ちるが。）
print(gcf,[path '_01.jpg'],'-djpeg','-r600');
end
% ===== Ｊ線が直線 =====
figure(2);
hf2=copyobj(ff0,gcf);
Xj2=linspace(150,-35,101);
Yj2=-0.9*Xj2+150;
plot(Xj2,Yj2,'LineWidth',0.5,'Color',aka);
axis([-190 360 -170 330]);
half_circle(Xj2,Yj2,ra,-1,1)
[Xjp,Yjp]=find_j_points(xt,yt,Xj2,Yj2,ra,-1);
plot(Xjp,Yjp,'Color',aka,'LineWidth',1.5)
draw_arm(xt(1),yt(1),Xjp(1),Yjp(1),0,a,rw,rt,0);
[~,dLw2]=displacement_of_wheel(xt,yt,Xjp,Yjp,ra,a);
% 静止画のファイル出力
if still_out==1
figure_trimmer(gcf,[-10 0 -30 -45]);
print(gcf,[path '_02.jpg'],'-djpeg','-r600');
end
% ===== Ｊ線が「く」の字 =====
figure(3);
hf3=copyobj(ff0,gcf);
x30=30; % 折れ曲がり点の座標
y30=50;
Xa=[x30:2:150];
Xa=fliplr(Xa);
Xb=[x30:-1:-35];
Xb(1)=[];
Xj3=[Xa Xb];
Yj3=max(-2*(Xj3-x30)+y30,-0.3*(Xj3-x30)+y30);
plot(Xj3,Yj3,'LineWidth',0.5,'Color',aka);
axis([-205 345 -200 300]);
half_circle(Xj3,Yj3,ra,-1,1)
[Xjp,Yjp]=find_j_points(xt,yt,Xj3,Yj3,ra,-1);
plot(Xjp,Yjp,'Color',aka,'LineWidth',1.5)
draw_arm(xt(1),yt(1),Xjp(1),Yjp(1),0,a,rw,rt,0);
[~,dLw3]=displacement_of_wheel(xt,yt,Xjp,Yjp,ra,a);
% 静止画のファイル出力
if still_out==1
figure_trimmer(gcf,[-10 0 -30 -45]);
print(gcf,[path '_03.jpg'],'-djpeg','-r600');
end
% ===== Ｊ線が放物線 =====
figure(4);
hf4=copyobj(ff0,gcf);
Xj4=[125:-2:-125];
Yj4=0.005*Xj4.^2;
thp=-40;
R=[cosd(thp) -sind(thp);sind(thp) cosd(thp)];
XY=R*[Xj4;Yj4]; % 回転
Xj4=XY(1,:)+10; % 平行移動
Yj4=XY(2,:)+45;
plot(Xj4,Yj4,'LineWidth',0.5,'Color',aka);
axis([-215 335 -220 280]);
half_circle(Xj4,Yj4,ra,-1,1)
[Xjp,Yjp]=find_j_points(xt,yt,Xj4,Yj4,ra,-1);
plot(Xjp,Yjp,'Color',aka,'LineWidth',1.5)
draw_arm(xt(1),yt(1),Xjp(1),Yjp(1),0,a,rw,rt,0);
[~,dLw4]=displacement_of_wheel(xt,yt,Xjp,Yjp,ra,a);
% 静止画のファイル出力
if still_out==1
figure_trimmer(gcf,[-10 0 -30 -45]);
print(gcf,[path '_04.jpg'],'-djpeg','-r600');
end
% ===== Ｊ線が正弦曲線 =====
figure(5);
hf5=copyobj(ff0,gcf);
Xj5=[125:-2:-125];
Yj5=-50*sind(Xj5*180/125);
thp=-25;
R=[cosd(thp) -sind(thp);sind(thp) cosd(thp)];
XY=R*[Xj5;Yj5];
Xj5=XY(1,:)+70;
Yj5=XY(2,:)+107;
plot(Xj5,Yj5,'LineWidth',0.5,'Color',aka);
axis([-160 390 -145 355]);
half_circle(Xj5,Yj5,ra,-1,1)
[Xjp,Yjp]=find_j_points(xt,yt,Xj5,Yj5,ra,-1);
plot(Xjp,Yjp,'Color',aka,'LineWidth',1.5)
draw_arm(xt(1),yt(1),Xjp(1),Yjp(1),0,a,rw,rt,0);
[~,dLw5]=displacement_of_wheel(xt,yt,Xjp,Yjp,ra,a);
% 静止画のファイル出力
if still_out==1
figure_trimmer(gcf,[-10 0 -30 -45]);
print(gcf,[path '_05.jpg'],'-djpeg','-r600');
end
% ===== Ｊ線の種類による面積測定値の違い =====
figure(6)
hS(1)=plot([0 Lt],[0 cumsum(dLw1)*ra/100],'LineWidth',1);
hold on
hS(2)=plot([0 Lt],[0 cumsum(dLw2)*ra/100],'LineWidth',1);
hS(3)=plot([0 Lt],[0 cumsum(dLw3)*ra/100],'LineWidth',1);
hS(4)=plot([0 Lt],[0 cumsum(dLw4)*ra/100],'LineWidth',1);
hS(5)=plot([0 Lt],[0 cumsum(dLw5)*ra/100],'LineWidth',1);
plot([0 700],[0 0],'k');
plot([0 700],[0 0]+230.5,'k');
plot([0 0]+Lt(end),[-30 320],'k');
axis([0 700 -30 320]);
grid on
text(260,242,'実面積 230.5 [cm^2]')
xlabel('測定開始点からの図形の外周の長さ [mm]')
ylabel('測定中の図形面積の推移 [cm^2]')
% 画像全体のトリミング（legend を書き込む前に実行しないとエラー）
figure_trimmer(gcf,[-10 0 0 -20]); % ローカル関数の呼び出し。
legend(hS,{'ポーラープラニメータ','リニアプラニメータ', ...
'くの字プラニメータ','放物線プラニメータ', ...
'正弦波プラニメータ'},'Location','northwest');
% 静止画のファイル出力
if still_out==1
% 画像をファイルに出力
% （本ケースでは、jpg よりも png が小容量。解像度も良い。）
print(gcf,[path '_06.png'],'-dpng','-r600');
end
% ===== アニメーション =====
if anime~=0 % アニメーションが必要なとき
filename=[path '_anime_' num2str(anime) '.gif'];
% アニメーションgifのファイル名。
delay=1/24; % １コマの秒数（24コマ/秒）。
% 最初と最後のコマの秒数は別途指定。
hfa=figure(7);
copyobj(ff0,gcf);
% 画像全体のトリミング
figure_trimmer(gcf,[-10 0 -25 -40]); % ローカル関数の呼び出し。
switch anime
case 1
Xja=Xj1; Yja=Yj1; Lwa=cumsum(dLw1);
case 2
Xja=Xj2; Yja=Yj2; Lwa=cumsum(dLw2);
case 3
Xja=Xj3; Yja=Yj3; Lwa=cumsum(dLw3);
case 4
Xja=Xj4; Yja=Yj4; Lwa=cumsum(dLw4);
case 5
Xja=Xj5; Yja=Yj5; Lwa=cumsum(dLw5);
end
plot(Xja,Yja,'LineWidth',1.5,'Color',aka); % Ｊ線を描画
[Xjp,Yjp]=find_j_points(xt,yt,Xja,Yja,ra,-1);
% 各Ｔ点に対応するＪ点の座標群を求める。
% ローカル関数の呼び出し。
[~,dLwa]=displacement_of_wheel(xt,yt,Xjp,Yjp,ra,a);
% 各ステップでの目盛車の回転方向への移動量。
% ローカル関数の呼び出し。
Area=cumsum(dLwa)*ra/100;
half_circle(Xja,Yja,ra,-1,0);
% Ｔ点の可動領域を示す半円群などを描画。
% ローカル関数の呼び出し。
for n=0:length(xt)
% ファイル容量低減のため、動画のコマの間引き
if ~(n==1 | n==length(xt)) & mod(n,2)==1
continue;
end
if n==0
m=1;
Sa=0;
elseif n==length(xt)
m=1;
Sa=Area(n-1);
else
m=n;
Sa=Area(n);
end
[ha]=draw_arm(xt(m),yt(m),Xjp(m),Yjp(m),Lwa(m),a,rw,rt,1);
hb=text(-45,-80,['面積 ' num2str(Sa,'%05.1f') ...
' [cm^2]'],'FontSize',18);
drawnow;
pause(0.01);
if anime_out==1
frames=getframe(hfa); % figure画像を取得
[AA, map] = rgb2ind(frame2im(frames), 256); % 画像形式変換
if n==0 % アニメの最初のコマの画像を書き出す。
imwrite(AA, map, filename, 'gif', 'DelayTime', 0.5, ...
'LoopCount',Inf);
% 表示は0.5秒間。１コマ目で繰返し再生回数の設定。
elseif n==length(xt)
imwrite(AA, map, filename, 'gif', 'DelayTime', 2.0, ...
'WriteMode', 'append');
% 表示は2秒間。
else
imwrite(AA, map, filename, 'gif', 'DelayTime', delay, ...
'WriteMode', 'append');
% ２コマ目以降には、"追記"の指示が必要。
end
end
if n~=length(xt)
delete([ha,hb]);
end
end
end
% ■■■■■■■■■■■■■■■■■■■■■■■■■
% ■■■■■■■ 以下はローカル関数 ■■■■■■■
% ■■■■■■■■■■■■■■■■■■■■■■■■■
%
% free_curve
% divide_line_into_fixed_short_segments
% half_circle
% find_j_points
% displacement_of_wheel
% draw_arm
% figure_trimmer
% =============================================
% スプライン曲線で、歪んだ円や開曲線を描くためのローカル関数。
% =============================================
function [X,Y]=free_curve(curve_path,angs,kang)
% スプライン曲線で描いた歪んだ円や開曲線の座標を取得する。
% 始点,終点に同一点を指定したときは閉曲線になる。
% 【入力】
% curve_path: 線を通過させたい代表点の[x,y]座標の対を縦に並べた配列
% [x1,y1 ; ... ; xn,yn]。
% angs: 始点/終点での線の傾き [degree]。NaNのときは自由端。
% スカラーのときは、始点/終点とも同一の傾き。
% 2要素の配列のときは、それぞれ始点,終点での傾き。
% kang: angsの効果が及ぶ範囲の広さ（適量の目安は1.0）
% 【出力】
% [X,Y]: 結果として得られた折線のx座標列とy座標列
if length(angs)==1
ang1=angs(1);
ang2=ang1;
else
ang1=angs(1);
ang2=angs(2);
end
if isnan(ang1)
dx1=0;
dy1=0;
else
dx1=cosd(ang1);
dy1=sind(ang1);
end
if isnan(ang2)
dx2=0;
dy2=0;
else
dx2=cosd(ang2);
dy2=sind(ang2);
end
N=size(curve_path,1); % 線の経路を指定する点の数
% スプライン曲線の長さの概算見積もり
xx=[];
yy=[];
for n=1:N
xx=[xx curve_path(n,1)];
yy=[yy curve_path(n,2)];
end
L=cumsum(hypot(diff(xx),diff(yy)));
len=L(end);
% 媒介変数
t=linspace(0,1,N)*len*kang;
xo=[dx1 xx dx2];
yo=[dy1 yy dy2];
% スプライン曲線として平滑化
Lsp = spline(t,[xo;yo]);
XY = ppval(Lsp, linspace(0,max(t),201));
X=XY(1,:); % 平滑化された線の
Y=XY(2,:); % 折点のx,y座標値
end
% =============================================
% 折線を「等長の短線分」の連なりとして再構成するためのローカル関数。
% =============================================
function [Xd,Yd]=divide_line_into_fixed_short_segments(X,Y,ds)
% 提示された折線（旧折線）を「等長の短線分」の連なりとして構成された
% 新折線に変換する。ただし、旧折線の終端に残った線分だけは、長さを
% 他と揃えることはできず、短い長さになる。
% 【入力】
% X,Y: 旧折線の折点群のx,yデータの配列。
% ds: 入力された折線は、この長さの線分で等分割された折線に
% 変わる。ただし、旧折線の終端に残った線分だけは長さを
% 揃えることはできず、ds未満の値になる。
% 始端と終端だけは旧折線と新折線で共用されるが、その他の
% 折点については、両折線間で関連性はない。
% 新折点は、旧折線の線上には乗るが、旧折点が新折線の
% 線上に乗るとは限らない。乗らないことの方が多い。
% 【出力】
% Xd,Yd: dsの長さに等分割された新折線の折点群のx,yデータの配列。
% 処理済みの旧折線データを逐次消去更新できるよう、別の変数にコピー
xc=X;
yc=Y;
Xd=X(1); % 新折線の折点群のx,y座標を並べた配列の初期値。
Yd=Y(1); % 折線の始点は、旧折線と新折線で同じ。
while true
x1=xc(1); % 旧折線の、処理中の線分の始点の座標
y1=yc(1);
x2=xc(2); % 〃 終点 〃
y2=yc(2);
dx=x2-x1; % 旧折線内の処理中の線分の
dy=y2-y1; % 始点・終点間の座標値の隔たり。
L21=hypot(dx,dy); % 旧折線内の処理中の線分の（残り）長さ。
if L21>ds % dsが旧折線の処理中の線分の残り長さより短い
% とき（次の新折点は残り線分の上にある）
% P1=(x1,y1), P2=(x2,y2), Px=(x3,y3)
%
% ds= |Px-P1| 新折線を構成する一本の短線分の固定長。
% L21=|P2-P1| 旧折線の処理中の線分の長さ
%
% P1 Px P2
% ● ds ● ●
% ↑ 新折線の ↑
% ｜ 次の折点 旧折線の
% ｜ 処理中の線分の終点
% 旧折線の処理中の
% 線分の始点。
% （P2,P1間の長さが、旧折線の本来の線分の長さよりも短い
% 残り線分のときは、新折線の決定済の最新折点P0と等しい）
x3=x1+dx*ds/L21; % 新折線の次の折点のx座標
y3=y1+dy*ds/L21; % 〃 y座標
Xd=[Xd x3]; % 新折点の座標を新折線の配列に追加する。
Yd=[Yd y3];
xc(1)=x3; % 旧折線の処理中の線分の始点を新折点に変更。
yc(1)=y3;
else % dsが旧折線の処理中の線分の残り長さより長い
% とき（次の新折点は旧折線の残り線分上には
% 存在しない）
% 旧折線の、今まで処理対象にしていた線分上には、探している新折点
% は存在しないことが分かったので、旧折線の処理対象線分を次の
% 線分に移動する。
% この場合、新しい折点を求める計算が少々複雑になる。この様子を
% 下図に示す。
% ケース１の場合は、そのまま、余弦定理を用いて新折点を決定でき
% るが、ケース２の場合には、旧折線の処理対象線分をさらに先に
% 移動して、その始点を新たなP1として置き直し、ケース１の状態
% になるまで繰り返す。
%
% P0=(x0,y0), P1=(x1,y1), P2=(x2,y2), Px=(x3,y3)
%
% ds= |Px-P0| 新折線を構成する一本の短線分の固定長。
% L1= |P0-P1| P0-P1の長さ
% Lx= |Px-P1| Px-P1の長さ
% L21=|P2-P1| 旧折線の処理中の線分の長さ
% θ= ∠P0.P1.P2
% = ∠P0.P1.Px
%
% 余弦定理 ds^2=L1^2+Lx^2-2*L1*Lx*cosθ より、Lxを求める。
%
% ==== ケース１（新折線の次の折点はPxで決定） ====
%
% P1: 旧折線の次の処理対象線分の始点
% ●
% L1 Lx Px: P2-P1線分上にあって
% ● ds ● P0点からdsだけ
% P0: 新折線の 離れた点
% 決定済の最新折点 ● P2
% P1点に対応する終点
%
% ==== ケース２（新折線の次の折点は不明。P2よりも先にある） ====
%
% P1
% ●
% L1 ● P2
% P0 ● ds ● Px
%
%
% 旧折線の処理中の線分の終点が、旧折線の終端に達していた場合。
% この先に固定長dsを満たす点はないので、新折線の終点を旧折線の
% 終点と一致させて処理を終わる。
if x2==X(end) & y2==Y(end)
Xd=[Xd X(end)];
Yd=[Yd Y(end)];
return
end
while true
xc(1)=[]; % 旧折線の用済みの折点（そこを始点とする
yc(1)=[]; % 線分上に新折点は存在しない）を削除。
% その結果、旧折線の未処理の全折点xc,ycの
% 番号は一つだけ繰り上がる。
x1=xc(1); % 旧折線のこれから処理する線分の始点の座標
y1=yc(1);
x2=xc(2); % 〃 終点 〃
y2=yc(2);
dx=x2-x1; % 旧折線のこれから処理する線分の始点・終点
dy=y2-y1; % 間の座標値の隔たり。
L21=hypot(dx,dy); % 旧折線のこれから処理する線分の(残り)長さ
x0=Xd(end); % 新折線の発見済の最新の折点の座標
y0=Yd(end);
% 上図のL1の値
L1=hypot(x1-x0,y1-y0);
% 上図のθ（ここではth）
th=atan2d(dy,dx)-atan2d(y0-y1,x0-x1);
th=abs(mod(th+180,360)-180); % 0～+180°の範囲に変換
% 上図のLx
Lx=L1*cosd(th)+sqrt(L1^2*(cosd(th)^2)+ds^2-L1^2);
if Lx<L21 % 上図の「ケース１」の場合
x3=x1+dx*Lx/L21; % 新折線の、新たに決定した折点。
y3=y1+dy*Lx/L21;
Xd=[Xd x3]; % 新たな折点を、新折線の配列に追加する。
Yd=[Yd y3];
xc(1)=x3; % 旧折線の処理中の線分の始点を、新折線の
yc(1)=y3; % 新たな折点に変更。
break;
else % 上図の「ケース２」の場合
% 旧折線の処理中の線分の終点が、旧折線の終端に達していた場合
% この先に固定長dsを満たす点はないので、新折線の終点を
% 旧折線の終点と一致させて処理を終わる。
if x2==X(end) & y2==Y(end)
Xd=[Xd X(end)];
Yd=[Yd Y(end)];
return
else % まだ旧折線の終端には達していないとき。
continue; % 旧折線の処理対象線分を一つ先に進める。
end
end % end of if Lx<L21
end % end of while true
end % end of if L21>ds
end % end of while true
end % end of function
% =============================================
% 指定された折線の各折点を中心にして
% 半円と侵入禁止線を描くローカル関数。
% =============================================
function half_circle(Xj,Yj,ra,dir,all)
% 【入力】
% Xj: 折線の x 座標群（行ベクトル）
% Yj: 〃 y 〃
% ra: 半円の半径
% dir: 引きアーム配置を指定のとき1、押しアーム配置のとき-1
% all: 1: 半円は、全折点を対象、0: 始点と終点のみを対象
% 【出力】
% なし
% Ｊ線上のＪ点を中心にして、Ｔ点が描く半円の角度範囲
thh=linspace(180,360,101);
Lr1=[];
Lr2=[];
% 拘束線上の各Ｊ点を中心に半円を描く
for n=[1:length(Xj)];
x0=Xj(n);
y0=Yj(n);
if n==1
dx=Xj(2)-Xj(1);
dy=Yj(2)-Yj(1);
else
dx=Xj(n)-Xj(n-1);
dy=Yj(n)-Yj(n-1);
end
ang=atan2d(dy,dx); % 各Ｊ点での拘束線の傾斜角
% 各Ｊ点を中心にして、Ｊ線の進行方向の後方側/前方側の半円を描く。
% （dir=-1 のとき後方、dir=1 のとき前方）
if all==1
hL=plot(x0+ra*cosd(ang+thh+dir*90),y0+ra*sind(ang+thh+dir*90));
hL.Color(4)=0.4; % 半透明色
elseif (n==1 | n==length(Xj))
plot(x0+ra*cosd(ang+thh+dir*90),y0+ra*sind(ang+thh+dir*90), ...
'Color','#aaa');
end
% Ｔ点の侵入禁止線の折線座標群に新たな折点座標を追加
if n~=1 & abs(ang-ang_old)>10 & abs(ang-ang_old)<180
Lr1=[Lr1;[NaN NaN]]; % 連続線でない部分は切り離す
Lr2=[Lr2;[NaN NaN]];
end
Lr1=[Lr1;[x0+ra*cosd(ang+180-90) y0+ra*sind(ang+180-90)]];
Lr2=[Lr2;[x0+ra*cosd(ang-90) y0+ra*sind(ang-90)]];
ang_old=ang;
end
% 侵入禁止線の描画。
% Ｊ線が折線近似で傾斜が不連続のため精度は不十分。
% 特に、ポーラープラニメータの内側の円弧。太線で誤魔化している。
if all==1
plot(Lr1(:,1),Lr1(:,2),'m','LineWidth',2)
plot(Lr2(:,1),Lr2(:,2),'m','LineWidth',2)
else
plot(Lr1(:,1),Lr1(:,2),'Color','#fcf','LineWidth',2)
plot(Lr2(:,1),Lr2(:,2),'Color','#fcf','LineWidth',2)
end
end
% =============================================
% Ｔ線とＪ線の座標群を与えて、
% Ｔ線を一周する期間中の各Ｔ点に対応するＪ点の座標群を求める。
% １つのＴ点に対応するＪ点が１つしか存在しない場合に限る。
% （Ｔ線が、二重ヘアライン領域や侵入禁止線から外れていれば大丈夫）
% =============================================
function [Xjp,Yjp]=find_j_points(Xt,Yt,Xjl,Yjl,ra,dir)
% 【入力】
% Xt: Ｔ線の x 座標群（行ベクトル）
% Yt: 〃 y 〃
% Xjl: Ｊ線の x 座標群（行ベクトル）
% Yjl: 〃 y 〃
% ra: アームの長さ
% dir: 引きアーム配置を指定のとき1、押しアーム配置のとき-1
% 【出力】
% Xjp: Ｊ点の x 座標群（行ベクトル）
% Yjp: 〃 y 〃
Xjp=[]; % 出力（Ｊ点の座標値）の初期値
Yjp=[];
for m=1:length(Xt) % Ｔ線の各折点を順番に処理。
% インデックス m を操作する必要があるため、
% for でなく、 while ループを使用。
% m 値はループの最後の部分で更新。
P0=[Xt(m) Yt(m)]; % Ｔ線の処理対象の折点Ｔの座標
XYi=[]; % 交点情報の保存用の行列
% 各行は、
% [交点のx座標, 同y座標,
% Ｊ線の線分ベクトルのx座標,
% 同y座標]
for n=1:length(Xjl)-1 % Ｊ線の各線分について繰り返す
% [Ｔ点が中心の、アーム長 ra を半径とした円]と
% [Ｊ線の線分]との交点（Ｊ点の候補）を計算する。
% P0: Ｔ点を示すベクトル
% P1: Ｊ線の一線分の始点を示すベクトル
% P2: 〃 終点 〃
% P10: P1-P0
% P21: P2-P1
% r: 円の半径（プログラム上は ra）
% PL: Ｊ線の線分の延長線上にある任意の点
% Pc: 円上にある任意の点
% a: PL 線上での位置を表す変数(P1点からの|P21|を1.0とした距離)
%
% PL = a*(P2-P1)+P1 ..... (1)
% (-∞ < a < ∞、線分上では 0 <= a <= 1)
% |Pc-P0| = r → |Pc-P0|^2 = r^2 ....... (2)
% 交点上では Pc=PL ゆえ、 (2)式に (1)式を代入すると、
% |Pc-P0|^2 = |PL-P0|^2 = |a*(P2-P1)+P1-P0|^2 = |a*P21+P10|^2
% = (a*P21+P10)・(a*P21+P10)
% = P21・P21*a^2 + 2*(P21・P10)*a + P10・P10 = r^2
% ゆえに、P21・P21*a^2 + 2*(P21・P10)*a + P10・P10 - r^2 =0
% A = P21・P21、B = 2*(P21・P10)、C = P10・P10 - r^2 とおくと、
% 交点における a の値は、a = (-B±sqrt(B^2^4*A*C))/(2*A)
P1=[Xjl(n) Yjl(n)]; % 対象線分の始点のベクトル
P2=[Xjl(n+1) Yjl(n+1)]; % 〃 終点 〃
P21=P2-P1; % 線分の始点から終点へ向くベクトル
P10=P1-P0; % 円中心から線分の始点へのベクトル
A=dot(P21,P21);
B=2*dot(P21,P10);
C=dot(P10,P10)-ra^2;
D=B^2-4*A*C;
if D<0
continue;
end
a1=(-B-sqrt(D))/(2*A);
a2=(-B+sqrt(D))/(2*A);
if a1>=0 & a1<=1 % １つ目の交点候補が線分上に
F=P1+a1*(P2-P1); % あるとき、その交点のベクトル
XYi=[XYi;[F P21]];
end
if a2>=0 & a2<=1 % ２つ目の交点候補が線分上に
F=P1+a2*(P2-P1); % あるとき、その交点のベクトル
XYi=[XYi;[F P21]];
end
end
exj=0; % とりあえず「交点なし」としておく
if ~isempty(XYi) % 交点候補があれば、さらに絞り込む
for k=1:size(XYi,1)
% 指定された J→T 線の向き(dir)に合致した交点をＪ点として採用
if sign(dot(P0-XYi(k,[1 2]),XYi(k,[3 4])))==dir
Xjp=[Xjp XYi(k,1)];
Yjp=[Yjp XYi(k,2)];
exj=exj+1; % 交点数
end
end
end
if exj~=1
error('Ｔ点に対応する適正なＪ点が見つかりません。')
end
end
end
% =============================================
% [Ｔ線の座標群]と
% [Ｔ線上のそれぞれのＴ点に対応するＪ点の座標群]を与えて、
% 目盛車の回転方向への積算移動量（面積に比例）を計算する。
% =============================================
function [Lw,dLw]=displacement_of_wheel(Xt,Yt,Xj,Yj,ra,a)
% 【入力】
% Xt: Ｔ線の x 座標群（行ベクトル）
% Yt: 〃 y 〃
% Xj: 上記に対応するＪ点の x 座標群（行ベクトル）
% Yj: 〃 y 〃
% ra: アームの長さ
% a: Ｊ端から目盛車までの長さ
% 【出力】
% Lw: Ｔ線を１周したときの、目盛車の回転方向への積算移動量
% dLw: １ステップ毎の目盛車の回転方向への移動量（行ベクトル）
Xta=(Xt(2:end)+Xt(1:end-1))/2; % Ｔ線の各線分の中心座標
Yta=(Yt(2:end)+Yt(1:end-1))/2;
Xja=(Xj(2:end)+Xj(1:end-1))/2; % 隣接するＪ点との平均座標
Yja=(Yj(2:end)+Yj(1:end-1))/2;
thw=atan2d(Yta-Yja,Xta-Xja); % アームの角度の区間平均値(度)
Xw=(a/ra)*(Xt-Xj)+Xj; % 目盛車の位置の x,y 座標
Yw=(a/ra)*(Yt-Yj)+Yj;
dXw=Xw(2:end)-Xw(1:end-1); % 目盛車のステップ間の
dYw=Yw(2:end)-Yw(1:end-1); % 位置変化ベクトルの x,y 成分
nXw=cosd(thw-90); % 目盛車の回転方向の
nYw=sind(thw-90); % 単位ベクトルの x,y 成分
dLw=dot([dXw' dYw'],[nXw' nYw'],2)';
% 目盛車のステップ間の
% 位置変化ベクトルの回転方向成分
Lw=sum(dLw); % 目盛車の回転方向への積算移動量
end
% =============================================
% Ｔ点とＪ点の座標を与えて、アームと付属品を描く
% =============================================
function [ha]=draw_arm(xt,yt,xj,yj,Lw,a,rw,rt,chr)
% 【入力】
% xt: Ｔ線の x 座標
% yt: 〃 y 〃
% xj: Ｊ点の x 座標
% yj: 〃 y 〃
% Lw: 目盛車の有効な移動量の積算値
% a: Ｊ点から目盛車までの距離
% rw: 目盛車の半径
% rt: レンズの半径
% chr: 文字の仕様
% 1: アニメ用、0: 静止画像用
% 【出力】
% ha: 描いた図のハンドル
tht=[0:36:360];
% レンズ
ha(1)=plot(xt+rt*cosd(tht),yt(1)+rt*sind(tht), ...
'Color','#0072BD','LineWidth',2);
% Ｊ点
ha(2)=plot(xj,yj,'o','MarkerSize',7,'Color','none', ...
'LineWidth',2,'MarkerFaceColor','#0072BD');
% レンズの中心
ha(3)=plot(xt,yt,'o','Color','#D95319', ...
'MarkerSize',2,'LineWidth',2);
tha=atan2d(yt-yj,xt-xj); % アームの角度
ra=norm([yt-yj,xt-xj]); % アームの長さ
xw=xj+a*cosd(tha); % 目盛車の中心の座標
yw=yj+a*sind(tha);
xg=xj+(ra-rt)*cosd(tha); % レンズ支持点の座標
yg=yj+(ra-rt)*sind(tha);
% アーム本体
ha(4)=plot([xj xg],[yj yg],'Color','#0072BD','LineWidth',3);
% 目盛車
ha(5)=plot(xw+rw*[-cosd(tha+90) cosd(tha+90)], ...
yw+rw*[-sind(tha+90) sind(tha+90)], ...
'Color','#0b0','LineWidth',3);
% 文字
if chr==1
pox=12;
poy=12;
fs=16;
else
pox=5;
poy=27;
fs=14;
end
ha(6)=text(xj+pox,yj+poy,'Ｊ','FontSize',fs, ...
'HorizontalAlign','center');
ha(7)=text(xt+pox,yt+poy,'Ｔ','FontSize',fs, ...
'HorizontalAlign','center');
% Ｔ点を中心とした半径がアーム長の円（あると邪魔な感じもする）
%thc=[0:2:360];
%ha(8)=plot(xt+ra*cosd(thc),yt+ra*sind(thc),'Color','#aaa');
end
% =============================================
% figure の上下左右の余白調整用のローカル関数。
% 注： figure が legend を含むとエラーになる。
% legend を追加する前に実行すること。
% =============================================
function figure_trimmer(hf,KK)
% 【入力】
% hf: 余白調整したいfigureのハンドル（呼び出し側では常にgcf固定）
% KK: 余白の変更量のベクトル。 単位は[pixel]。
% （[上余白 下余白 左余白 右余白]。増量は+、減量は-）
% 【出力】
% 戻り値は無し
posf=hf.Position; % 余白変更前のfigureのposition
% [左端 下端 幅 高さ]、単位は[pixel]
naxe=size(hf.Children,1); % figure内のChildren（axesやcolorbar）
% などの総数。
% 全Childrenのposition行列を作成
% 各行がそれぞれのChildrenに対応。
% 各列は[左端 下端 幅 高さ]、単位は[normalized]。
Posa=[]; % position行列の初期値
for n=1:naxe
axes(hf.Children(n))
hn(n)=gca;
Posa=[Posa ; hn(n).Position];
end
% 各Childrenの、余白変更前のfigure内でのpixel単位でのPositionを計算
Posp(:,[1 3])=Posa(:,[1 3])*posf(3);
Posp(:,[2 4])=Posa(:,[2 4])*posf(4);
% 余白変更後のfigureのpixel単位でのPositionを計算。
posfx=posf+[-KK(3) -KK(2) KK(3)+KK(4) KK(1)+KK(2)];
% 余白変更後のfigure内での各Childrenのpixel単位でのPositionを計算。
Pospx(:,1)=Posp(:,1)+KK(3);
Pospx(:,2)=Posp(:,2)+KK(2);
Pospx(:,3)=Posp(:,3);
Pospx(:,4)=Posp(:,4);
% 余白変更後の各Childrenのnormalized単位でのPositionを計算する。
Posax(:,[1 3])=Pospx(:,[1 3])/posfx(3);
Posax(:,[2 4])=Pospx(:,[2 4])/posfx(4);
% figureとChildrenのpositionを更新
hf.Position=posfx;
for nn=1:naxe
hn(nn).Position=Posax(nn,:);
end
end