使用解析几何方法作出椭圆sino图——matlab实现

作者：老李
日期：11-1

目标：用解析几何的方法作出椭圆的sino图

需要使用的知识点：
1.符号型变量
2.表达式的求解
3.由矩阵生成灰度图

理论部分：

深刻椭圆表达式

$$(X-\mu {)}^T\sum (X-\mu )=r$$
1.X为向量（自变量），X的长度体现为空间的维度。

2.参数$\mu$决定了椭圆的中心点位置

3.对$\sum$的理解

首先，对$\sum$的思考非常有助于对矩阵的理解，是线性代数里非常好的一道题，而且我们还能把这方面的理解拓展到对多维正态分布的理解上。

$\sum$为一个满秩正定实对称矩阵，它决定了该椭圆的形状。
当$\sum$为一个对角矩阵时，它所表现的是一个正规的椭圆（可以理解为单位圆或球往每个轴的方向上的拉伸）
当$\sum$不是一个对角矩阵时，它所表现的是一个由正规的椭圆以中心位置旋转之后的椭圆。我们知道，实对称矩阵通过正交变换可以得到一组新的正交基。从这一点的角度上理解，那就是任意一个椭圆（球）在某一个正交基上可以表现为正规椭圆。

我们知道，广义化的正态分布的表达式为：
$$p(X)=\frac{1}{(2\pi {)}^{n/2}\left|{\sum }^{1/2}\right|}exp\left\{-\frac{1}{2}(X-\mu {)}^T{\Sigma }^{-1}(X-\mu )\right\}$$指数部分里面实际上是一个椭圆的表达式，应为我们知道，协方差矩阵$\Sigma$是正定的，这也完全表达了该概率密度的形状

2.matlab实现：
好吧，我承认我在理论部分扯那么多只是为了掩饰我编程能力的弱小。。。
因为怕麻烦，因为由理论部分我们可以知道，复杂一点的椭圆无非就是对原表达式进行几次线性变换而已，所以我们选用的椭圆为一个非常简单的正规的椭圆。表达式为：$$\frac{x^2}{9}+\frac{y^2}{4}=1$$图像为：

我们知到对于直线的表达式可以表示为$$cos(\theta )x+sin(\theta )y=r$$我们要计量的是法线角度为0到$\pi$之间的情况，而$r$为直线与原点之间的距离。

在这里插入代码片%定义符号型变量来定义表达式，从而求解
syms x y a b r c;
%求解椭圆与直线联立求解的表达式
[x,y] = solve('x^2/9+y^2/4= 1,a*x+b*y = r');
%将截线长用符号表达出来
c = sqrt((x(1)-x(2))^2+(y(1)^2-y(2))^2);
ezplot('x^2/9+y^2/4 = 1')%作椭圆图
B = zeros(180,length(-3:0.01:3));
v = -3:0.01:3;
for i = 1:180
    for k = 1:length(v)
        q = (i-1)*pi/180;
        d = subs(c,[a;b;r],[sin(q);-cos(q);cos(q)*v(k)]);%以实数替代符号
        g = vpa(d);%求值
        %滤掉虚数部分，只取有实数解的部分
        if abs(imag(g))<eps
            B(i,k) = g;
        end
    end
end

为了避免麻烦，我在写代码的过程中，尽量避免了表达式的推导。

结果，这个程序运行了2000多秒。
然后我优化了一下程序，决定不用subs语句,于是我把表达式先求解了出来。
$(((r-(2\ast b\ast (2\ast b\ast r-3\ast a\ast (9\ast a^2+4\ast b^2-r^2{)}^{(}1/2)))/(9\ast a^2+4\ast b^2))/a-(r-(2\ast b\ast (2\ast b\ast r+3\ast a\ast (9\ast a^2+4\ast b^2-r^2{)}^{(}1/2)))/(9\ast a^2+4\ast b^2))/a{)}^2+((2\ast (2\ast b\ast r-3\ast a\ast (9\ast a^2+4\ast b^2-r^2{)}^{(}1/2)))/(9\ast a^2+4\ast b^2)-(2\ast (2\ast b\ast r+3\ast a\ast (9\ast a^2+4\ast b^2-r^2{)}^{(}1/2)))/(9\ast a^2+4\ast b^2){)}^2{)}^{(}1/2)$
代码如下：

%(((b*r + a*(4*a^2 + 4*b^2 - r^2)^(1/2))/(a^2 + b^2) - (b*r - a*(4*a^2 + 4*b^2 - r^2)^(1/2))/(a^2 + b^2))^2 + ((r - (b*(b*r + a*(4*a^2 + 4*b^2 - r^2)^(1/2)))/(a^2 + b^2))/a - (r - (b*(b*r - a*(4*a^2 + 4*b^2 - r^2)^(1/2)))/(a^2 + b^2))/a)^2)^(1/2)
B = zeros(180,length(-3:0.01:3));
v = -3:0.01:3;
for i = 1:180
    for k = 1:length(v)
        q = (i-1)*pi/180+pi/2;
        a = sin(q);
        b = -cos(q);
        r = v(k)*cos(q);
        g = (((r - (2*b*(2*b*r - 3*a*(9*a^2 + 4*b^2 - r^2)^(1/2)))/(9*a^2 + 4*b^2))/a - (r - (2*b*(2*b*r + 3*a*(9*a^2 + 4*b^2 - r^2)^(1/2)))/(9*a^2 + 4*b^2))/a)^2 + ((2*(2*b*r - 3*a*(9*a^2 + 4*b^2 - r^2)^(1/2)))/(9*a^2 + 4*b^2) - (2*(2*b*r + 3*a*(9*a^2 + 4*b^2 - r^2)^(1/2)))/(9*a^2 + 4*b^2))^2)^(1/2);
        %滤掉虚数部分，只取有实数解的部分
       if imag(g)<eps
            B(i,k) = g;
       end
    end
end
%作图
imshow(B,[])

同理，我把表达式变了一下$$x^2+y^2=4$$
然后我们得到了新的表达式

`B = zeros(180,length(-3:0.01:3));
v = -3:0.01:3;
for i = 1:180
    for k = 1:length(v)
        q = (i-1)*pi/180+pi/2;
        a = sin(q);
        b = -cos(q);
        r = v(k)*cos(q);
        g = (((r + (b*(a*(a^2 + b^2 - r^2)^(1/2) - b*r))/(a^2 + b^2))/a - (r - (b*(a*(a^2 + b^2 - r^2)^(1/2) + b*r))/(a^2 + b^2))/a)^2 + ((a*(a^2 + b^2 - r^2)^(1/2) + b*r)/(a^2 + b^2) + (a*(a^2 + b^2 - r^2)^(1/2) - b*r)/(a^2 + b^2))^2)^(1/2);
        %滤掉虚数部分，只取有实数解的部分
       if imag(g)<eps
            B(i,k) = g;
       end
    end
end
%作图
imshow(B,[])

图如下：