跟上次放的代码相比，有变复杂的地方也有变简单的部分，复杂的是系数的计算，因为去根号会产生很多的项，简单的是矩阵的生成，直接把平移旋转缩放连乘即可，无需计算基向量。
<!DOCTYPE html>
<html>
<head>
<title>矩阵法计算二个椭圆的交点</title>
<script src="Matrix.js"></script>
<script src="MatrixUtil.js"></script>
<script src="Point.js"></script>
<script src="FormulaUtil.js"></script>
<script src="ComplexNum.js"></script>
<script src="AngleUtil.js"></script>
<script src="ComplexConsts.js"></script>
</head>
<body>
<canvas width="800" height="800" id="canvas"></canvas>
</body>
<script>
function toScreen(p)
{
return new Point(coordO.x + p.x * unitSize, coordO.y - p.y * unitSize);
}
var unitSize = 160;
var coordO = new Point(200, 400);
//椭圆中心
var ellipseCenter0 = new Point(0.5, -0.5);
//长轴半径
var aLength0 = 1.5;
//短轴半径
var bLength0 = 1;
//旋转角度
var rotation0 = Math.PI / 6;
var matrix0 = new Matrix();
//把椭圆变换逐个加入
//缩放
MatrixUtil.scale(matrix0, aLength0, bLength0);
//旋转
MatrixUtil.rotate(matrix0, rotation0);
//平移
MatrixUtil.translate(matrix0, ellipseCenter0.x, ellipseCenter0.y);
//椭圆中心
var ellipseCenter1 = new Point(1, 0.5);
//长轴半径
var aLength1 = 2;
//短轴半径
var bLength1 = 1;
//旋转角度
var rotation1 = - Math.PI / 3;
var matrix1 = new Matrix();
//把椭圆变换逐个加入
//缩放
MatrixUtil.scale(matrix1, aLength1, bLength1);
//旋转
MatrixUtil.rotate(matrix1, rotation1);
//平移
MatrixUtil.translate(matrix1, ellipseCenter1.x, ellipseCenter1.y);
//被转成单位圆的逆矩阵跟另一个椭圆的自身矩阵相乘
var matrix = matrix0.clone();
matrix.invert();
matrix.prepend(matrix1);
matrix.invert();
//矩阵的参数
var a = matrix.a;
var b = matrix.b;
var c = matrix.c;
var d = matrix.d;
var tx = matrix.tx;
var ty = matrix.ty;
//由X^2+Y^2=1和(ax*t^2+bx*t+cx)^2+(ay*t^2+by*t+cy)^2=1
//展开化简将获得一个四次方程
//计算系数的中间变量（因为两个椭圆的联立需要解无理方程，非常繁琐），此处先存下一些平方前的变量
//注意到y^2系数没了，因为y^2=1-x^2，所以可被直接化走
var XS = a * a + b * b - c * c - d * d; //x^2的系数
var XK = 2 * a * tx + 2 * b * ty;//x的系数
var CK = tx * tx + ty * ty + c * c + d * d - 1;//常数项
var XYK = - 2 * a * c - 2 * b * d;//xy的系数
var YK = - 2 * c * tx - 2 * d * ty;//y的系数
//四次方程的五个系数
var fourFormulaA = XS * XS + XYK * XYK;
var fourFormulaB = 2 * XS * XK + 2 * XYK * YK;
var fourFormulaC = XK * XK - XYK * XYK + YK * YK + 2 * XS * CK;
var fourFormulaD = 2 * XK * CK - 2 * XYK * YK;
var fourFormulaE = CK * CK - YK * YK;
//用四次求根的类进行求解
var resolutions = FormulaUtil.calcFourFormulaZero(fourFormulaA, fourFormulaB, fourFormulaC, fourFormulaD, fourFormulaE);
var realResolutions = FormulaUtil.getRealSolutions(resolutions);
var intersections = [];
for(var i = 0, len = realResolutions.length; i < len; i ++)
{
var x = realResolutions[i];
//算出在基向量矩阵逆变换后的交点
var p = new Point(x, Math.sqrt(1 - x * x));
//因为这个交点是方程在基向量逆变换后求得，所以要执行基向量矩阵的原始变换才能让结果恢复到变换前
intersections.push(matrix0.transformPoint(p));
}
var canvas = document.getElementById("canvas");
var context = canvas.getContext("2d");
context.lineWidth = 0.5;
context.strokeStyle = "#0000cc";
//绘制椭圆（百度回来的，用绘制正圆+缩放实现）
context.save();
var screenCenter = toScreen(ellipseCenter0);
context.translate(screenCenter.x, screenCenter.y);
context.rotate(-rotation0);
context.scale(aLength0 * unitSize, bLength0 * unitSize);
context.lineWidth = 0.5 / unitSize;
context.strokeStyle = "#0000cc";
context.beginPath();
context.arc(0,0,1,0,Math.PI*2,true);
context.stroke();
context.closePath();
context.restore();
//绘制椭圆（百度回来的，用绘制正圆+缩放实现）
context.save();
var screenCenter = toScreen(ellipseCenter1);
context.translate(screenCenter.x, screenCenter.y);
context.rotate(-rotation1);
context.scale(aLength1 * unitSize, bLength1 * unitSize);
context.lineWidth = 0.5 / unitSize;
context.strokeStyle = "#0000cc";
context.beginPath();
context.arc(0,0,1,0,Math.PI*2,true);
context.stroke();
context.closePath();
context.restore();
context.fillStyle = "#cc0000";
//绘制交点
for(var i = 0, len = intersections.length; i < len; i ++)
{
var toScreenP = toScreen(intersections[i]);
context.beginPath();
context.arc(toScreenP.x, toScreenP.y, 5, 0, Math.PI * 2, true);
context.fill();
context.closePath();
}
</script>
</body>
</html>运行效果如下图所示。
有一个点没有落在交点的位置，这让我昨晚很惆怅，检查了好多遍去根号，矩阵连乘顺序等需要手动的代码，都未发现任何异常，然后我就去睡觉了。白天到花市逛了下，然后现在继续回来看。
为了确认这个正确的结果是否巧合使然，我换了椭圆的好些属性值来试，发现大多数时候都能匹配到其中一个点，当然也有一些时候全匹配上和全匹配不上。
终于发现了自己的一个低级错误，求出来的解x所对应的y有两个，它们互为相反数，但我却只取了正的那个。错误的代码是这段。
for(var i = 0, len = realResolutions.length; i < len; i ++)
{
var x = realResolutions[i];
//算出在基向量矩阵逆变换后的交点
var p = new Point(x, Math.sqrt(1 - x * x));
//因为这个交点是方程在基向量逆变换后求得，所以要执行基向量矩阵的原始变换才能让结果恢复到变换前
intersections.push(matrix0.transformPoint(p));
}我们要把y取反的结果看也加到交点数组中。
for(var i = 0, len = realResolutions.length; i < len; i ++)
{
var x = realResolutions[i];
//算出在基向量矩阵逆变换后的交点
var p = new Point(x, Math.sqrt(1 - x * x));
//因为这个交点是方程在基向量逆变换后求得，所以要执行基向量矩阵的原始变换才能让结果恢复到变换前
intersections.push(matrix0.transformPoint(p));
var p = new Point(x, -Math.sqrt(1 - x * x));
//因为这个交点是方程在基向量逆变换后求得，所以要执行基向量矩阵的原始变换才能让结果恢复到变换前
intersections.push(matrix0.transformPoint(p));
}再次运行，效果如下图所示，所有交点都获取到了。
但与此同时也有多余的点产生，因为正和负一般只有一个是正确的。
所以我们还得像初中解无理方程一样，判断下哪些是增根。从几何角度上说，就是是否在另一个椭圆上。
直接代入椭圆方程来算太麻烦，我们不妨也用矩阵变换，把另一个椭圆变成单位圆再去计算，用的就是上面给到的matrix变量。废话不说了，直接上代码。
for(var i = 0, len = realResolutions.length; i < len; i ++)
{
var x = realResolutions[i];
//算出在基向量矩阵逆变换后的交点
var p = new Point(x, Math.sqrt(1 - x * x));
var tempP = matrix.transformPoint(p);
if(Math.abs(tempP.x * tempP.x + tempP.y * tempP.y - 1) <= 0.01)
{
//因为这个交点是方程在基向量逆变换后求得，所以要执行基向量矩阵的原始变换才能让结果恢复到变换前
intersections.push(matrix0.transformPoint(p));
}
var p = new Point(x, -Math.sqrt(1 - x * x));
var tempP = matrix.transformPoint(p);
if(Math.abs(tempP.x * tempP.x + tempP.y * tempP.y - 1) <= 0.01)
{
//因为这个交点是方程在基向量逆变换后求得，所以要执行基向量矩阵的原始变换才能让结果恢复到变换前
intersections.push(matrix0.transformPoint(p));
}
}注意到我加了个误差值范围，毕竟浮点数很难100%准确。
再次运行，效果就对了。
通过椭圆求交，我们真正看到了矩阵的巧妙之处——把复杂图形中的矩阵部分抽出来将其简化，然后让多个矩阵进行合并计算，从而有效降低了复杂图形的计算难度。相信到了本篇，大家已经可以轻松驾驭任何二次曲线（包括从未提过的双曲线）的求交计算了。
接下来，我们将进军本教程的终极目标——二元二次方程组一般式的求解。继续用矩阵的话，那它跟曲线求交的最大区别就在于：需要先判断曲线类型，然后找到变换为标准方程的矩阵，接着运用前面学过的知识进行计算，敬请期待！
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
