基于Matlab的跨孔CT胖射线追踪算法（五）

商功贤

已于 2022-04-22 17:39:16 修改

阅读量716

点赞数

分类专栏：射线追踪文章标签： matlab 算法

于 2021-09-29 20:16:26 首次发布

本文链接：https://blog.csdn.net/Neverlevsun/article/details/120555509

版权

本文介绍了基于Matlab的跨孔CT胖射线追踪算法，通过MATLAB代码展示了胖射线追踪过程，并用动画进行动态演示。胖射线追踪是CT技术的一种正演算法，适用于无损检测，具有较高分辨率。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

基于Matlab的跨孔CT胖射线追踪算法（五）

CT技术是一种无损的工程物探检测技术，因其方法简单、分辨率高、理论上更可靠、结果更直观，被广泛的应用于各种工程。胖射线追踪是CT技术的一种正演算法，本文的代码和示意图，供学习借鉴。¹

文章目录

基于Matlab的跨孔CT胖射线追踪算法（五）

胖射线追踪

在这里插入图片描述

用动画显示胖射线追踪的过程。
在这里插入图片描述

MATLAB代码

逐条射线追踪实现胖射线追踪。

clear 
close all
clc

xmin = 0;
xmax = 10;
ymin = 0;
ymax = 10;
nex = 10;
ney = 10;
dx = 1;
dy = 1;
    
x1 = 0;
y1 = 1.5;
x2 = 10;
y2 = 5.5;


clear xx yy;hold on
[xx,yy] = meshgrid(xmin:dx:xmax,ymin:dy:ymax);
plot(xx,yy,'k:',xx',yy','k:');
line([x1,x2],[y1,y2],'color','r');
scatter(x1,y1,'r*');scatter(x2,y2,'rs');


syms x y

a1 = y2-y1;
b1 = x1-x2;
c1 = x2*y1-x1*y2;
f1 = @(x,y) a1*x+b1*y+c1;


lamda = 5; 
c = (1/2)*sqrt((y2-y1)^2+(x2-x1)^2);
b = sqrt(lamda*(sqrt((y2-y1)^2+(x2-x1)^2))/8);
a = sqrt(c^2 + b^2);
theta = atan((y2-y1)/(x2-x1));
x0 = (x1+x2)/2;
y0 = (y1+y2)/2;
f2 = @(x,y) (a^2-c^2*cos(theta)^2)*(x-x0).^2 + ...
     (a^2-c^2*sin(theta)^2)*(y-y0).^2 - ...
     c^2*sin(2*theta)*(x-x0).*(y-y0) - a^2*b^2;
h = ezplot(f2,[xmin-2*dx,xmax+2*dx,ymin-b,ymax+b]);
set(h,'Color','b','LineStyle','-.')
delete(get(gca,'title'));
xlabel('x/m');
ylabel('y/m');
set(gca,'FontSize',15);
grid off
axis equal

xx1 = xx(1,:);
if x1 < x2,xxline = xx1((xx1 > x1) & (xx1 < x2 ));end
if x1 > x2,xxline = xx1((xx1 > x2) & (xx1 < x1 ));end

yy1 = yy(:,1)';
if y1 < y2,yyline = yy1((yy1 > y1) & (yy1 < y2 ));end
if y1 > y2,yyline = yy1((yy1 > y2) & (yy1 < y1 ));end




if x1 ~= x2 && y1 ~=y2
    xpx = xxline;
    ypx = xpx*(-a1/b1)-c1/b1;
    ypy = yyline;
    xpy = ypy*(-b1/a1)-c1/a1;    
    xp0 = [xpx,xpy];
    yp0 = [ypx,ypy];
    xp = unique(roundn(xp0',-4),'rows','stable');
    yp = unique(roundn(yp0',-4),'rows','stable');
    pp = [xp,yp];
    plot(xp,yp,'bo')
end

if x1 == x2
    xxline = [];
    yp = yyline';
    xp = repmat(x1,[size(yp,1),1]