反射

最近学线性代数，想到了一种新的方法计算所有二维空间（可推广到n维）的线反射。相比较之前的旋转坐标系方法简单很多。

假设有直线[latex]l_1=\vec{P}_{0}+\vec{P}t[/latex]作为入射线，[latex]l_2=\vec{Q}_{0}+\vec{Q}d[/latex]作为介质，求反射线[latex]l_3=\vec{R}_{0}+\vec{R}s[/latex]。

[latex]l_1[/latex]到[latex]l_2[/latex]的投影为：
[latex]proj_{l_{2}}(l_1)=\vec{Q}(\vec{Q}^{T}\vec{Q})^{-1}\vec{Q}^{T}\vec{P}[/latex]

$latex R $可通过投影和P的关系求得：
$latex \vec{R}=2proj_{l_{2}}(l_1)-\vec{P} $

$latex \vec{R}_{0} $是$latex l_1$和$latex l_2$的交点（解出$latex t_0$和$latex d_0$从而得到$latex \vec{R}_{0} $）：

$latex \left[\begin{array}{cc}\vec{P}&-\vec{Q}\end{array}\right] \left[\begin{array}{c}t_{0}\\d_{0}\end{array}\right]=\left[\begin{array}{c}\vec{Q}_{0}-\vec{P}_{0}\end{array}\right]$

[latex]\left[\begin{array}{c}t_{0}\\d_{0}\end{array}\right]=\left[\begin{array}{cc}\vec{P}&-\vec{Q}\end{array}\right]^{-1}\left[\begin{array}{c}\vec{Q}_{0}-\vec{P}_{0}\end{array}\right][/latex]

$latex \vec{R}_{0}=\vec{P}_{0}+\vec{P}t_{0}=\vec{Q}_{0}+\vec{Q}d_{0}$

下面是matlab代码，可以快速计算反射线。可以输入以标量或矢量模式输入直线。

display 'This program will calculate the reflection line equation by inputing two lines\n.';

%input with two methods
voa= input('Would you like to input using Vector[v] or Scalar[s] form? (v by default): ','s');

if voa=='s'

    pk=0;pb=0;qk=0;qb=0;
    display 'please input lines in such format y=k*x+b\n';
    pk=input('please input k for line1 (entering line): ');
    pb=input('please input b for line1 (entering line): ');
    qk=input('please input k for line2 (middle line): ');
    qb=input('please input b for line2 (middle line): ');

    %convert the line equations into vector format;
    p0=[0;pb];
    q0=[0;qb];
    p=[1;pk];
    q=[1;qk];
else
    display 'please input lines in such format [p0x;p0y]+[px;py]*t\n';
    p0x=input('please input p0x: ');
    p0y=input('please input p0y: ');
    px=input('please input px: ');
    py=input('please input py: ');
    q0x=input('please input q0x: ');
    q0y=input('please input q0y: ');
    qx=input('please input qx: ');
    qy=input('please input qy: ');

    %organize into vector form
    p0=[p0x;p0y];
    q0=[q0x;q0y];
    p=[px;py];
    q=[qx;qy];
end

%calculate using line projection
proj=q*inv(q'*q)*q'*p;
r=2*proj-p;

%calculate r0 by solving the matrix
tdmat=[p -q p0-q0];
rtdmat=rref(tdmat);
t=-rtdmat(1,3);
r0=p0+p*t;

%draw the lines
d=linspace(-20,20);
d1=zeros(1,100)+1;
l1=p0*d1+p*d;
l2=q0*d1+q*d;
l3=r0*d1+r*d;

figure;
plot(l1(1,:),l1(2,:));
hold;
plot(l2(1,:),l2(2,:),'g');
plot(l3(1,:),l3(2,:),'r');
axis equal;

%print result
if(voa=='s')
    rk=[r(2)/r(1)];
    rb=r0(2)-rk*r0(1);
    if rb>=0
        fprintf('this resulting line is y=%f*x+%f\n',rk,rb);
    else
        fprintf('this resulting line is y=%f*x%f\n',rk,rb);
    end
else
    fprintf('this resulting line is [%f;%f]+[%f;%f]*u\n',r0(1),r0(2),r(1),r(2));
end