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

投影矩阵的计算_投影矩阵的几何意义在进行迭代重建的过程中,我们首先需要求出投影矩阵之后才能进行其他后续的操作,在迭代重建中起到了基石的作用。并且在前面的文章中《迭代重建算法中投影矩阵的计算》已经给出了一种方法,但是我发现在程序的运行过程中存在一些未知的bug,导致程序在计算某些角度的投影矩阵时出现错误。由于一直没有找到出现bug的原因,因此我改变了计算思路,找到了下文中正确的计算方法。首先需要证明一条直线与一个正方形相交。假设一个正方形的左上角的顶点坐标为(xk,yk),那么其余三个点的坐标也就能够写出来…

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新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);

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

         由上图可知,程序运行正确!!!

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

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

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/186090.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

  • matplotlib绘图实例:pyplot、pylab模块及作图参数「建议收藏」

    matplotlib绘图实例:pyplot、pylab模块及作图参数「建议收藏」http://blog.csdn.net/pipisorry/article/details/40005163Matplotlib.pyplot绘图实例{使用pyplot模块}matplotlib绘制直线、条形/矩形区域importnumpyasnpimportmatplotlib.pyplotaspltt=np.arange(-1,2,.01)s=…

  • idea2021.10激活码永久【2021最新】

    (idea2021.10激活码永久)好多小伙伴总是说激活码老是失效,太麻烦,关注/收藏全栈君太难教程,2021永久激活的方法等着你。IntelliJ2021最新激活注册码,破解教程可免费永久激活,亲测有效,下面是详细链接哦~https://javaforall.cn/100143.html1STL5S9V8F-eyJsaWNlbnNlSWQi…

  • 批处理命令for循环_批处理获取某个目录大小

    批处理命令for循环_批处理获取某个目录大小使用批处理批量读取目录中的文件名,并根据文件名,使用for循环对文件逐一进行处理。

  • async和await

    async和awaitasync和await概念先从字面意思来理解。async是“异步”的简写,而await可以认为是asyncwait的简写。所以应该很好理解async用于申明一个function是异步的,而await用于等待一个异步方法执行完成。另外还有一个很有意思的语法规定,await只能出现在async函数中。然后细心的朋友会产生一个疑问,如果await只能出现在async函数中,那这个async函数应该怎么调用?如果需要通过await…

  • android自动化测试框架_GMS测试

    android自动化测试框架_GMS测试GMSP的测试已经开始一段时间了,这里记录一下跟O的不同之处,资料主要来自MTK官网。一、Case数量增加(以Pgo项目为例,64bit项目double)a、CTS测试增加8万条,module增至320b、增加CTS-INSTANT测试项目,module48,1.3万条(13338)、64bit13338,以实际为准;c、VTS测试增加约2…

  • 『机器学习笔记 』GBDT原理-Gradient Boosting Decision Tree

    『机器学习笔记 』GBDT原理-Gradient Boosting Decision Tree1.背景1.1GradientBoosting1.2提升树-boostingtree回归问题提升树算法2GradientBoostingDecisionTree2.1函数空间的数值优化2.2算法Shrinkage总结附录参考资料相似算法:1.背景决策树是一种基本的分类与回归方法。决策树模型具有分类速度快,模型…

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号