提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  大数据

当前栏目

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

计算 矩阵 意义 几何 投影
2023-06-13 09:14:51 时间

在进行迭代重建的过程中，我们首先需要求出投影矩阵之后才能进行其他后续的操作，在迭代重建中起到了基石的作用。并且在前面的文章中《迭代重建算法中投影矩阵的计算》已经给出了一种方法，但是我发现在程序的运行过程中存在一些未知的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);

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

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

版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 举报，一经查实，本站将立刻删除。

发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/186090.html原文链接：https://javaforall.cn

相关文章