投影矩阵的计算_投影矩阵的几何意义

大家好，又见面了，我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE稳定放心使用

        在进行迭代重建的过程中，我们首先需要求出投影矩阵之后才能进行其他后续的操作，在迭代重建中起到了基石的作用。并且在前面的文章中《迭代重建算法中投影矩阵的计算》已经给出了一种方法，但是我发现在程序的运行过程中存在一些未知的bug，导致程序在计算某些角度的投影矩阵时出现错误。由于一直没有找到出现bug的原因，因此我改变了计算思路，找到了下文中正确的计算方法。

        首先需要证明一条直线与一个正方形相交。假设一个正方形的左上角的顶点坐标为(xk,yk)，那么其余三个点的坐标也就能够写出来，分别为(xk+1,yk)、(xk+1,yk-1)、(xk,yk-1)。如果(m*xk+b-yk)*(m*(xk+1)+b-yk)<=0，那么容易知道该直线与(xk,yk)、(xk+1,yk)两点确定的直线相交，对其他三条边也是这样操作。

        接下来的问题时如何求解一条直线被一个正方形所截线段的长度。依然利用上一段的方法，将两条相交的直线联立方程组，分别求出直线与正方形的两个交点坐标。然后通过两点之间的坐标公式计算所截线段的长度。

        最后通过代码实现上述的数学思想，并将其写成一个函数文件，方便以后调用。（在这里面delta最好设置为1，表示一个像素的大小）

function [W_ind,W_dat]=medfuncSystemMatrix1(theta,N,P_num,delta)
%%==定义参数==%%
%theta：探测器的旋转角度，角度
%N：头模型的大小
%P_num：探测器的总数目
%delta：每个网格的大小
%%==输出参数==%%
%W_ind：存储射线穿过的网格的编号
%W_dat：存储射线被穿过网格所截断的长度
N2=N^2;%编号总数
theta=theta*pi/180;
M=length(theta)*P_num;%投影射线总条数
W_ind=zeros(M,2*N);%存放射线穿过的网格的编号
W_dat=zeros(M,2*N);%存放射线穿过的的网格的长度
t=(-(P_num-1)/2:(P_num-1)/2)*delta;%探测器的坐标
% if N<=10 && length(theta)<=5
%     x=(-N/2:1:N/2)*delta;
%     y=(-N/2:1:N/2)*delta;
%     plot(x,meshgrid(y,x),'k');
%     hold on;
%     plot(meshgrid(x,y),y,'k');
%     axis([-N/2-5,N/2+5,-N/2-5,N/2+5]);
%     text(0,-0.4*delta,'0');
% end
%%==投影矩阵的计算==%%
for jj=1:length(theta)  %对角度进行遍历
    th=theta(jj);
%     for ii=1:P_num      %对探测器进行遍历
%         th=theta(jj);
%         u=zeros(1,2*N);  %存储编号
%         v=zeros(1,2*N);  %存储长度
        if th>=pi || th<0
             error('输入角度必须在0~180之间');%使用error进行报错提醒
        elseif th==pi/2
            for ii=1:P_num
                u=zeros(1,2*N);
                v=zeros(1,2*N);
                x=t(ii);%此时第ii个探测器射线的方程
                c=0;
                for k=1:N^2
                    if rem(k,N)==0
                        x1=-N/2+(N-1)*delta;
                    else
                        x1=-N/2+(rem(k,N)-1)*delta;
                    end
                    x2=x1+delta;
                    if x>=x1 && x<x2
                        c=c+1;
                        u(c)=k;
                        v(c)=delta;
                    else
                        continue;
                    end
                end
                W_ind((jj-1)*P_num+ii,:)=u;
                W_dat((jj-1)*P_num+ii,:)=v;
            end
       
        elseif th==0
            for ii=1:P_num
                u1=zeros(1,2*N);
                v1=zeros(1,2*N);
                y=t(ii);%此时第ii个探测器的坐标
                c=0;
                for k=1:N2
                    if rem(k,N)==0
                        y1=N/2-(floor(k/N)-1)*delta;
                    else
                        y1=N/2-floor(k/N)*delta;
                    end
                    y2=y1-delta;
                    if y<y1 && y>=y2
                        c=c+1;
                        u1(c)=k;
                        v1(c)=delta;
                    else 
                        continue;
                    end
                end
                W_ind((jj-1)*P_num+ii,:)=u1;
                W_dat((jj-1)*P_num+ii,:)=v1;
            end
        elseif th<pi && th>pi/2
            m=tan(th);%直线斜率
            %b=-t/cos(th);%直线截距
            for ii=1:P_num
                u1=zeros(1,2*N);
                v1=zeros(1,2*N);
                t1=t(ii);
                b=t1/cos(th);%直线截距
                c=0;
                for k=1:N2
                    %%==求编号k的左上顶点的坐标
                    if rem(k,N)==0
                        xk=-N/2+(N-1)*delta;
                        yk=N/2-(floor(k/N)-1)*delta;
                    else
                        xk=-N/2+(rem(k,N)-1)*delta;%编号k左上角的x坐标
                        yk=N/2-floor(k/N)*delta;%编号k左上角的y坐标
                    end
                    if (m*xk+b-yk)*(m*(xk+1)+b-yk)<=0 || (m*(xk+1)+b-yk)*(m*(xk+1)+b-(yk-1))<=0 || (m*(xk+1)+b-(yk-1))*(m*xk+b-(yk-1))<=0 || (m*xk+b-(yk-1))*(m*xk+b-yk)<=0 %如果成立，那么必有交点
                        c=c+1;
                        u1(c)=k;%将该正方形的编号传入数组,编号求解是正确的。
                        %接下来求被正方形截断的长度
                        if (m*xk+b-yk)*(m*(xk+1)+b-yk)<=0 && (m*(xk+1)+b-yk)*(m*(xk+1)+b-(yk-1))<=0  
                            x1=(yk-b)/m;%此时交点坐标（x1,yk）
                            y1=m*(xk+1)+b;%此时交点坐标（xk+1,y1）
                            v1(c)=sqrt((x1-xk-1)^2+(yk-y1)^2);
                        elseif (m*xk+b-yk)*(m*(xk+1)+b-yk)<=0 && (m*(xk+1)+b-(yk-1))*(m*xk+b-(yk-1))<=0
                            x2=(yk-b)/m;%此时交点坐标（x2,yk）
                            x3=(yk-1-b)/m;%此时交点坐标（x3,yk-1）
                            v1(c)=sqrt((x3-x2)^2+1);
                        elseif (m*xk+b-(yk-1))*(m*xk+b-yk)<=0 && (m*(xk+1)+b-(yk-1))*(m*xk+b-(yk-1))<=0
                            y2=m*xk+b;%此时交点坐标（xk,y2）
                            x4=(yk-1-b)/m;%此时交点坐标(x4,yk-1)
                            v1(c)=sqrt((xk-x4)^2+(yk-1-y2)^2);
                        elseif (m*xk+b-(yk-1))*(m*xk+b-yk)<=0 && (m*(xk+1)+b-yk)*(m*(xk+1)+b-(yk-1))<=0
                            y3=m*xk+b;%此时交点坐标（xk,y3）
                            y4=m*(xk+1)+b;%此时交点坐标（xk+1,y4）
                            v1(c)=sqrt((y3-y4)^2+1);    
                        end
                    else
                        continue;
                    end
                end
                W_ind((jj-1)*P_num+ii,:)=u1;
                W_dat((jj-1)*P_num+ii,:)=v1;
            end
        elseif th<pi/2 && th>0
            m=tan(th);
            for ii=1:P_num
                u1=zeros(1,2*N);
                v1=zeros(1,2*N);
                t1=t(ii);
                b=t1/cos(th);%直线截距
                c=0;  
                for k=1:N2
                    if rem(k,N)==0
                        xk=-N/2+(N-1)*delta;
                        yk=N/2-(floor(k/N)-1)*delta;
                    else
                        xk=-N/2+(rem(k,N)-1)*delta;%编号k左上角的x坐标
                        yk=N/2-floor(k/N)*delta;%编号k左上角的y坐标
                    end
                    if (m*xk+b-yk)*(m*(xk+1)+b-yk)<=0 || (m*(xk+1)+b-yk)*(m*(xk+1)+b-(yk-1))<=0 || (m*(xk+1)+b-(yk-1))*(m*xk+b-(yk-1))<=0 || (m*xk+b-(yk-1))*(m*xk+b-yk)<=0 %如果成立，那么必有交点
                        c=c+1;
                        u1(c)=k;%将该正方形的编号传入数组
                        %接下来求被正方形截断的长度
                        if (m*xk+b-yk)*(m*xk+b-(yk-1))<=0 && (m*xk+b-yk)*(m*(xk+1)+b-yk)<=0
                            y1=m*xk+b;%此时交点坐标为（xk,y1）
                            x1=(yk-b)/m;%此时交点坐标（x1,yk）
                            v1(c)=sqrt((xk-x1)^2+(yk-y1)^2);
                        elseif (m*xk+b-yk)*(m*xk+b-(yk-1))<=0 && (m*(xk+1)+b-yk)*(m*(xk+1)+b-(yk-1))<=0
                            y2=m*xk+b;%此时交点坐标为（xk,y2）
                            y3=m*(xk+1)+b;%此时交点坐标（xk+1,y3）
                            v1(c)=sqrt(1+(y3-y2)^2);
                        elseif (m*xk+b-(yk-1))*(m*(xk+1)+b-(yk-1))<=0 && (m*xk+b-yk)*(m*(xk+1)+b-yk)<=0
                            x2=(yk-b)/m;%此时交点坐标（x2,yk）
                            x3=(yk-1-b)/m;%此时交点坐标（x3,yk-1）
                            v1(c)=sqrt((x3-x2)^2+1);
                        elseif (m*xk+b-(yk-1))*(m*(xk+1)+b-(yk-1))<=0 && (m*(xk+1)+b-yk)*(m*(xk+1)+b-(yk-1))<=0
                            x4=(yk-1-b)/m;%此时交点坐标（x4,yk-1）
                            y4=m*(xk+1)+b;%此时交点坐标（xk+1,y4）
                            v1(c)=sqrt((x4-xk-1)^2+(yk-1-y4)^2);
                        end
                    else
                        continue;
                    end
                end
                W_ind((jj-1)*P_num+ii,:)=u1;
                W_dat((jj-1)*P_num+ii,:)=v1;
            end
        end
end
end
                
           
                    

        下面测试该程序的正确性：

theta=45;
N=2;
P_num=4;
delta=1;
[W_ind,W_dat]=medfuncSystemMatrix1(theta,N,P_num,delta);

         由上图可知，程序运行正确！！！

        本文提到的方法很容易理解，并且实现起来也比较方便，但是唯一的缺点就是程序运行效率较低，但是结果是没有任何问题的。