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

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

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

胖射线追踪

MATLAB代码

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

clear 
close all
clc

xmin = 0;
xmax = 5;
ymin = 0;
ymax = 5;
nex = 5;
ney = 5;
dx = 1;
dy = 1;
    
x1 = 0;
y1 = 0.85;
x2 = 5;
y2 = 4.2;

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 = 2; 
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-dx,xmax+dx,ymin-b,ymax+b]);
set(h,'Color','b','LineStyle','-.')
delete(get(gca,'title'));
xlabel('孔间距/m');
ylabel('深度/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]);
    pp = [xp,yp];
    plot(xp,yp,'bo')
end
if y1 == y2
    yyline = [];
    xp = xxline';
    yp = repmat(y1,[size(xp,1),1]);
    pp = [xp,yp];
    plot(xp,yp,'bo')
end
nxyline = size(xxline,2) + size(yyline,2);

rs = sqrt( (pp(:,2)-y1 ).^2 + ( pp(:,1)-x1 ).^2 );
ppp = [pp,rs];

ppp = sortrows(ppp,3);

mp = size(ppp,1) + 1;
xm = zeros(mp,1);
ym = zeros(mp,1);
rm = zeros(mp,1);
xm(1) = ( ppp(1,1) + x1 )/2;
ym(1) = ( ppp(1,2) + y1 )/2;
rm(1) = sqrt( ( ppp(1,1) - x1 )^2 ...
            + ( ppp(1,2) - y1 )^2 );
xm(mp) = ( ppp(mp-1,1) + x2 )/2;
ym(mp) = ( ppp(mp-1,2) + y2 )/2;
rm(mp) = sqrt( ( x2 - ppp(mp-1,1))^2 ...
             + ( y2 - ppp(mp-1,2) )^2 );
for i = 2:size(ppp,1)
    xm(i) = ( ppp(i,1) + ppp(i-1,1) )/2;
    ym(i) = ( ppp(i,2) + ppp(i-1,2) )/2;
    rm(i) = sqrt( ( ppp(i,1) - ppp(i-1,1) )^2 ...
                + ( ppp(i,2) - ppp(i-1,2) )^2 );
end

clear iex iey;
iex=0;iey=0;
xh = zeros(mp,1);
for i = 1:mp
    for ix = 1:nex+1
        if ( xm(i) > xx1(ix) ) && ( xm(i) < xx1(ix+1) )
            iex = ix;
            break;
        end
    end
    for iy = 1:ney
        if ( ym(i) > yy1(iy) ) && ( ym(i) < yy1(iy+1) )
            iey = iy;
            break;
        end
    end
    xh(i) = (iex-1)*ney + iey;
end
A = [xh,rm];
ppp1 = [ppp(:,1) ,ppp(:,2)];
ppw = [x1,y1;ppp1;x2,y2];

ps = zeros(2*size(ppw,1),2);
if y1 ~= y2
    k = b1/a1;
    for i = 1: size(ppw,1)
        fs =@(x,y) k*( x - ppw(i,1) ) + ppw(i,2) -y;
        s = solve( fs,f2,x,y);
        xs = double( s.x);
        ys = double( s.y);
        plot(xs,ys,'b:')
        scatter(xs,ys,'r+')
        ps(2*i-1:2*i,:) = [xs,ys];
        pause(0.05)
    end
end
if y1 == y2
    for i = 1: size(ppw,1)
        fs =@(x,y) x - ppw(i,1) + 0*y;
        s = solve( fs,f2,x,y);
        xs = double( s.x);
        ys = double( s.y);
        plot(xs,ys,'b:'),hold on
        scatter(xs,ys,'r+')
        ps(2*i-1:2*i,:) = [xs,ys];
        pause(0.05)
    end
end

figure(2)
pv = [ps(1:2,:);ps(4,:);ps(3,:);ps(1,:)];
[p,t]=distmesh2d(@dpoly,@huniform,0.1,[-1,-1; 12,12],pv,pv);
figure(1),hold on
for i = 1:size(t,1)
       fill([p(t(i,1),1),p(t(i,2),1),p(t(i,3),1)],...
       [p(t(i,1),2),p(t(i,2),2),p(t(i,3),2)],i)
       hold on
end
 

 for i = 1:size(ppw,1)-1
     figure(2)
     pv = [ps(2*i-1:2*i,:);ps(2*i+2,:);ps(2*i+1,:);ps(2*i-1,:)];
     [p,t]=distmesh2d(@dpoly,@huniform,0.1,[-1,-1; 12,12],pv,pv);
     figure(1),hold on
     for j = 1:size(t,1)
       h = fill([p(t(j,1),1),p(t(j,2),1),p(t(j,3),1)],...
           [p(t(j,1),2),p(t(j,2),2),p(t(j,3),2)],j);
       set(h,'edgealpha','0.3','facealpha',0.6)
       colormap('jet')
       hold on
     end
end

smz = size(t,1);

p1 = [p(t(1:1:smz,1),1),p(t(1:1:smz,1),2)];
p2 = [p(t(1:1:smz,2),1),p(t(1:1:smz,2),2)];
p3 = [p(t(1:1:smz,3),1),p(t(1:1:smz,3),2)];
[xmz,ymz] = scentriod(p1,p2,p3);
clear iex iey;
iex=0;iey=0;
mh = zeros(t,1);
for i = 1:mp
     for ix = 1:nex+1
         if ( xmz(i) > xx1(ix) ) && ( xmz(i) < xx1(ix+1) )
             iex = ix;
             break;
        end
     end
     for iy = 1:ney
         if ( ymz(i) > yy1(iy) ) && ( ymz(i) < yy1(iy+1) )
             iey = iy;
             break;
         end
     end
     mh(i) = (iex-1)*ney + iey;
end

早期练习阶段写的程序，里面有许多要修改的地方，但是还是可以运行出来结果的。

---

1. 本文仅供学习交流。 ↩︎