楕円と線の交点
楕円と線の交点を解析的に直接求める方法はありません．
一方で，円と線の交点は求められます．
そこで，楕円を円にするように世界を変形し，その中で円と線の交点を求め，その後元の世界に変形し返すというプロセスを踏むことで，間接的に楕円と線の交点を解析的に求められます．
注：図形のベクトル表現について
本記事では，利便性から図形のベクトル表現を採用しています．
直線：ある点から，任意に伸縮するベクトルを伸ばして到達できる点の集合 $\mathbf{p}_l + s \mathbf{d}_l, \;\; s \in \mathbb{R}$
円：中心点から，半径$r$分の長さの任意に回転するベクトルを伸ばして到達できる点の集合 $\mathbf{p}_c + r \begin{bmatrix} \cos(t) \\ \sin(t) \end{bmatrix}, \;\; t \in [-\pi, \pi]$
正立楕円：単位円を左右に伸ばしたもの $\mathbf{p}_o + \boldsymbol{S}_o \begin{bmatrix} \cos(t) \\ \sin(t) \end{bmatrix}, \;\; \boldsymbol{S}_o = \begin{bmatrix} r_x & 0 \\ 0 & r_y \end{bmatrix}, \;\; t \in [-\pi, \pi]$
楕円：正立楕円を回転したもの $\mathbf{p}_e + \boldsymbol{R}_{\theta_e} \, \boldsymbol{S}_e \begin{bmatrix} \cos(t) \\ \sin(t) \end{bmatrix}, \;\; \boldsymbol{S}_e = \begin{bmatrix} r_x & 0 \\ 0 & r_y \end{bmatrix}, \;\; t \in [-\pi, \pi]$ （$\boldsymbol{R}_\theta$：$\theta$回転する回転行列，$\theta_e$：楕円の偏角）
例題：
\displaylines{
l: \mathbf{p}_l + s \mathbf{d}_l \\
e: \mathbf{p}_e + \boldsymbol{R}_{\theta_e} \, \boldsymbol{S}_e \begin{bmatrix} \cos(t) \\ \sin(t) \end{bmatrix}
}
①楕円の中心位置が世界座標系の中心になるように平行移動する
\displaylines{
l: (\mathbf{p}_l-\mathbf{p}_e) + s \mathbf{d}_l \\
e: \boldsymbol{R}_{\theta_e} \, \boldsymbol{S}_e \begin{bmatrix} \cos(t) \\ \sin(t) \end{bmatrix}
}
②楕円が座標軸に対して正立するように回転する
\displaylines{
l: \boldsymbol{R}_{-\theta_e} \left((\mathbf{p}_l-\mathbf{p}_e) + s \mathbf{d}_l\right) = \boldsymbol{R}_{-\theta_e} \, \mathbf{p}_{l-e} + s \boldsymbol{R}_{-\theta_e} \, \mathbf{d}_l = \mathbf{p}'_{l-e} + s \mathbf{d}'_l \\
e: \boldsymbol{S}_e \begin{bmatrix} \cos(t) \\ \sin(t) \end{bmatrix}
}
③楕円が単位円になるようにスケールする
\displaylines{
l: \boldsymbol{S}^{-1}_e (\mathbf{p}'_{l-e} + s \mathbf{d}'_l) = \mathbf{p}''_{l-e} + s \mathbf{d}''_l\\
e: \begin{bmatrix} \cos(t) \\ \sin(t) \end{bmatrix}
}
④変形した世界で交点を求める
・単位円と線の交点
高々2点で交わるので，線上の2点を決める係数$s^\pm=(s^+,s^-)$を求める．交点は$\mathbf{p}''_{l-e} + s^+ \mathbf{d}''_l$ と $\mathbf{p}''_{l-e} + s^- \mathbf{d}''_l$
\mathbf{p}''_{l-e} + s^\pm \mathbf{d}''_l = \begin{bmatrix} \cos(t) \\ \sin(t) \end{bmatrix}
・要素ごとに2乗（$\odot$：要素ごとの積）
(\mathbf{p}''_{l-e} + s^\pm \mathbf{d}''_l)\odot (\mathbf{p}''_{l-e}+s^\pm \mathbf{d}''_l) = \begin{bmatrix}\cos^2(t) \\ \sin^2(t)\end{bmatrix}
・両辺で各要素を足し，1つの2次方程式にする
(p''_{l-e,x} + s^\pm d''_{l,x})^2 + (p''_{l-e,y} + s^\pm d''_{l,y})^2 = \cos^2(t) + \sin^2(t) = 1
・$s^\pm$について解く
\displaylines{
{s^\pm}^2 ({d''_{l,x}}^2+{d''_{l,y}}^2) + 2s^\pm (d''_{l,x}p''_{l-e,x}+d''_{l,y}p''_{l-e,y}) + ({p''_{l-e,x}}^2+{p''_{l-e,y}}^2 -1) = 0 \\
{s^\pm}^2 a + 2s^\pm b + c = 0 \\
s^\pm = \frac{-b\pm\sqrt{b^2-ac}}{a}
}
※$b^2-ac<0$の場合は交点無し
⑤スケールを戻す
\boldsymbol{S}_{e} (\mathbf{p}''_{l-e} + s^\pm \mathbf{d}''_l)
④元の向きに回転する
\boldsymbol{R_{\theta_e}} \boldsymbol{S}_e (\mathbf{p}''_{l-e} + s^\pm \mathbf{d}''_l)
④元の位置に平行移動する
\mathbf{p}_e+\boldsymbol{R_{\theta_e}} \boldsymbol{S}_e (\mathbf{p}''_{l-e} + s^\pm \mathbf{d}''_l)
※$s^\pm$を正のみに限定すると，半直線を扱えます．レイキャスティングにどうぞ．
コード
C++，回転しない版．利便性からOpenCVを使っています．
// 正立楕円と半直線の手前側の交点を求める
// lp: 線の原点, ld: 線の方向ベクトル, ecp: 楕円の中心点, esize: 楕円の外接矩形のサイズ, [out] cp: 交点
// return: 交点があったらtrue, なければfalse
bool orthoEllipseLineCross(cv::Vec2d lp, cv::Vec2d ld, cv::Vec2d ecp, cv::Vec2d esize, cv::Vec2d &cp)
{
// 直接求めるのはムリなので，
// 楕円が円になるように世界を変形させて，線と円の交点を求め，また世界を元の形に戻す
esize /= 2.0;
cv::Vec2d elp = lp - ecp; //平行移動
elp[0] /= esize[0];
elp[1] /= esize[1]; //スケーリング
cv::Vec2d eld = ld;
eld[0] /= esize[0];
eld[1] /= esize[1]; //スケーリング
eld /= sqrt(eld[0] * eld[0] + eld[1] * eld[1]); //正規化
// lp+s*ld=cp
// lpx + s*ldx = cosd
// lpy + s*ldy = sind
// lpx^2+2s*ldx*lpx+s^2*ldx^2 + lpy^2+2s*ldy*lpy+s^2*ldy^2 = 1
// s^2*(ldx^2+ldy^2) + 2*s*(ldx*lpx + ldy*lpy) + lpx^2+lpy^2-1 = 0
// s = (-b +- sqrt(b^2-ac)) /a
double a = eld[0] * eld[0] + eld[1] * eld[1]; // ldx^2 + ldy^2
double b = eld[0] * elp[0] + eld[1] * elp[1]; // ldx*lpx + ldy*lpy
double c = elp[0] * elp[0] + elp[1] * elp[1] - 1; // lpx^2 + lpy^2 - 1
double D = b * b - a * c;
if (D < 0)
return false; // 交点無し
double sp = (-b + sqrt(D)) / a;
double sn = (-b - sqrt(D)) / a;
// lpに近い方の点を求める ldと逆向きの点は不採用
double s;
if (sp < 0 && sn < 0)
return false;
else if (sp >= 0 && sn < 0)
s = sp;
else if (sp < 0 && sn >= 0)
s = sn;
else
s = fabs(sp) < fabs(sn) ? sp : sn;
cv::Vec2d dcp = elp + s * eld;
// 世界の変形を戻す
//スケーリング
dcp[0] *= esize[0];
dcp[1] *= esize[1];
//平行移動
dcp[0] += ecp[0];
dcp[1] += ecp[1];
//結果を返す
cp = dcp;
return true;
}
C++，回転する版．
// 楕円と半直線の手前側の交点を求める
// lp: 線の原点, ld: 線の方向ベクトル, ecp: 楕円の中心点, esize: 楕円の外接矩形のサイズ, earg: 楕円の偏角, [out] cp: 交点
// return: 交点があったらtrue, なければfalse
bool ellipseLineCross(cv::Vec2d lp, cv::Vec2d ld, cv::Vec2d ecp, cv::Vec2d esize, double earg, cv::Vec2d &cp)
{
// 直接求めるのはムリなので，
// 楕円が円になるように世界を変形させて，線と円の交点を求め，また世界を元の形に戻す
//正立楕円→楕円の回転行列
cv::Matx<double,2,2> R;
R(0, 0) = cos(earg);
R(1, 0) = sin(earg);
R(0, 1) = -sin(earg);
R(1, 1) = cos(earg);
//楕円→正立楕円の回転行列
cv::Matx<double,2,2> iR = R.inv();
//
esize /= 2.0;
cv::Vec2d elp = lp - ecp; // 平行移動
elp = iR * elp; // 回転
elp[0] /= esize[0];
elp[1] /= esize[1]; // スケーリング
cv::Vec2d eld = ld;
eld = iR * ld; // 回転
eld[0] /= esize[0];
eld[1] /= esize[1]; // スケーリング
eld /= sqrt(eld[0] * eld[0] + eld[1] * eld[1]); // 正規化
// lp+s*ld=cp
// lpx + s*ldx = cosd
// lpy + s*ldy = sind
// lpx^2+2s*ldx*lpx+s^2*ldx^2 + lpy^2+2s*ldy*lpy+s^2*ldy^2 = 1
// s^2*(ldx^2+ldy^2) + 2*s*(ldx*lpx + ldy*lpy) + lpx^2+lpy^2-1 = 0
// s = (-b +- sqrt(b^2-ac)) /a
double a = eld[0] * eld[0] + eld[1] * eld[1]; // ldx^2 + ldy^2
double b = eld[0] * elp[0] + eld[1] * elp[1]; // ldx*lpx + ldy*lpy
double c = elp[0] * elp[0] + elp[1] * elp[1] - 1; // lpx^2 + lpy^2 - 1
double D = b * b - a * c;
if (D < 0)
return false; // 交点無し
double sp = (-b + sqrt(D)) / a;
double sn = (-b - sqrt(D)) / a;
// lpに近い方の点を求める ldと逆向きの点は不採用
double s;
if (sp < 0 && sn < 0)
return false;
else if (sp >= 0 && sn < 0)
s = sp;
else if (sp < 0 && sn >= 0)
s = sn;
else
s = fabs(sp) < fabs(sn) ? sp : sn;
cv::Vec2d dcp = elp + s * eld;
// 世界の変形を戻す
//スケーリング
dcp[0] *= esize[0];
dcp[1] *= esize[1];
//回転
dcp = R * dcp;
//平行移動
dcp[0] += ecp[0];
dcp[1] += ecp[1];
//結果を返す
cp = dcp;
return true;
}