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

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

CT技术是一种无损的工程物探检测技术，因其方法简单、分辨率高、理论上更可靠、结果更直观，被广泛的应用于各种工程。胖射线追踪是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]);
    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
        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


for ic = 1:size(ppw,1)-1
m1 = ps(2*ic-1,:);
m2 = ps(2*ic,:);
m3 = ps(2*ic+1,:);
m4 = ps(2*ic+2,:);


clear xxmx1 xxmy1 yymx1 yymy1
if m1(1,1)>m2(1,1)
    xxmx1 = xx1(xx1 < m1(1,1) & xx1 > m2(1,1));
    xxmy1 = (b1/a1)*(xxmx1 - m1(1,1) ) + m1(1,2);
end
if m1(1,1)<m2(1,1)
    xxmx1 = xx1(xx1 > m1(1,1) & xx1 < m2(1,1));
    xxmy1 = (b1/a1)*(xxmx1 - m1(1,1) ) + m1(1,2);
end
if m1(1,1) == m2(1,1),xxmx1 = [];xxmy1 = [];end
if m2(1,2)>m1(1,2),yymy1 = yy1(yy1 < m2(1,2) & yy1 > m1(1,2));end
if m2(1,2)<m1(1,2),yymy1 = yy1(yy1 > m2(1,2) & yy1 < m1(1,2));end
if x1 ~= x2,yymx1 = (a1/b1)*(yymy1 - m1(1,2) ) + m1(1,1);
elseif x1 == x2,yymx1 = ones(1,size(yymy1,2))*m1(1,1);end


clear xxmx2 xxmy2 yymx2 yymy2
if m3(1,1)>m4(1,1)
    xxmx2 = xx1(xx1 < m3(1,1) & xx1 > m4(1,1));
    xxmy2 = (b1/a1)*(xxmx2 - m3(1,1) ) + m3(1,2);
end
if m3(1,1)<m4(1,1)
    xxmx2 = xx1(xx1 > m3(1,1) & xx1 < m4(1,1));
    xxmy2 = (b1/a1)*(xxmx2 - m3(1,1) ) + m3(1,2);
end
if m3(1,1) == m4(1,1),xxmx2 = [];xxmy2 = [];end
if m4(1,2)>m3(1,2),yymy2 = yy1(yy1 < m4(1,2) & yy1 > m3(1,2));end
if m4(1,2)<m3(1,2),yymy2 = yy1(yy1 > m4(1,2) & yy1 < m3(1,2));end
if x1 ~= x2,yymx2 = (a1/b1)*(yymy2 - m3(1,2) ) + m3(1,1);
elseif x1 == x2,yymx2 = ones(1,size(yymy2,2))*m3(1,1);end

pm1 = [xxmx1',xxmy1';yymx1',yymy1';xxmx2',xxmy2';yymx2',yymy2'];


clear xxmx3 xxmy3 yymx3 yymy3
if m4(1,1)>m2(1,1),xxmx3 = xx1(xx1 < m4(1,1) & xx1 > m2(1,1));end
if m4(1,1)<m2(1,1),xxmx3 = xx1(xx1 > m4(1,1) & xx1 < m2(1,1));end
if size(xxmx3,2) == 1
    fs1 =@(x,y) x - xxmx3 + 0*y;
    s1 = solve( fs1,f2,[x,y]);
    ys1 = double( s1.y );
    [~,n]=min(abs(ys1(:)-m2(1,2)));
    xxmy3=ys1(n);
elseif size(xxmx3,2) == 2
    fs11 =@(x,y) x - xxmx3(1,1)+0*y;fs12 =@(x,y) x - xxmx3(1,2)+0*y;
    s11 = solve( fs11,f2,[x,y]);s12 = solve( fs12,f2,[x,y]);
    ys11 = double( s11.y );ys12 = double( s12.y );
    [~,n]=min(abs(ys11(:)-m2(1,2)));xxmy3(1,1)=ys11(n);
    [~,n]=min(abs(ys12(:)-m2(1,2)));xxmy3(1,2)=ys12(n);
elseif isempty(xxmx3),xxmy3 = [];
end

if m4(1,2)>m2(1,2),yymy3 = yy1(yy1 < m4(1,2) & yy1 > m2(1,2));end
if m4(1,2)<m2(1,2),yymy3 = yy1(yy1 > m4(1,2) & yy1 < m2(1,2));end
if size(yymy3,2) == 1
    fs2 =@(x,y) 0*x + y -yymy3;
    s2 = solve( fs2,f2,[x,y]);
    xs2 = double( s2.x );
    [~,n]=min(abs(xs2(:)-m2(1,1)));
    yymx3=xs2(n);
elseif size(yymy3,2) == 2
    fs21 =@(x,y) 0*x + y -yymy3(1,1);fs22 =@(x,y) 0*x + y -yymy3(1,2);
    s21 = solve( fs21,f2,x,y);s22 = solve( fs22,f2,x,y);
    xs21 = double( s21.x );xs22 = double( s22.x );
    [~,n]=min(abs(xs21(:)-m2(1,1)));yymx3(1,1)=xs21(n);
    [~,n]=min(abs(xs22(:)-m2(1,1)));yymx3(1,2)=xs22(n);
elseif isempty(yymy3),yymx3 = [];
end


clear xxmx4 xxmy4 yymx4 yymy4
if m3(1,1)>m1(1,1),xxmx4 = xx1(xx1 < m3(1,1) & xx1 > m1(1,1));end
if m3(1,1)<m1(1,1),xxmx4 = xx1(xx1 > m3(1,1) & xx1 < m1(1,1));end
if size(xxmx4,2) == 1
    fs3 =@(x,y) x - xxmx4 + 0*y;
    s3 = solve( fs3,f2,[x,y]);
    ys3 = double( s3.y );
    [~,n]=min(abs(ys3(:)-m1(1,2)));
    xxmy4=ys3(n);
elseif size(xxmx4,2) == 2
    fs31 =@(x,y) x - xxmx4(1,1) + 0*y;fs32 =@(x,y) x - xxmx4(1,2) + 0*y;
    s31 = solve( fs31,f2,[x,y]); s32 = solve( fs32,f2,[x,y]);
    ys31 = double( s31.y );ys32 = double( s32.y );
    [~,n]=min(abs(ys31(:)-m1(1,2)));xxmy4(1,1)=ys31(n);
    [~,n]=min(abs(ys32(:)-m1(1,2)));xxmy4(1,2)=ys32(n);
elseif isempty(xxmx4),xxmy4 = [];
end

if m3(1,2)>m1(1,2),yymy4 = yy1(yy1 < m3(1,2) & yy1 > m1(1,2));end
if m3(1,2)<m1(1,2),yymy4 = yy1(yy1 > m3(1,2) & yy1 < m1(1,2));end
if size(yymy4,2) == 1
    fs4 =@(x,y) 0*x + y -yymy4;
    s4 = solve( fs4,f2,[x,y]);
    xs4 = double( s4.x );
    [~,n]=min(abs(xs4(:)-m1(1,1)));
    yymx4=xs4(n);
elseif size(yymy4,2) == 2
    fs41 =@(x,y) 0*x + y -yymy4(1,1);fs42 =@(x,y) 0*x + y -yymy4(1,2);
    s41 = solve( fs41,f2,[x,y]);s42 = solve( fs42,f2,[x,y]);
    xs41 = double( s41.x );xs42 = double( s42.x );
    [~,n]=min(abs(xs41(:)-m1(1,1)));yymx4(1,1)=xs41(n);
    [~,n]=min(abs(xs42(:)-m1(1,1)));yymx4(1,2)=xs42(n);
elseif isempty(yymy4),yymx4 = [];
end
pm2 = [xxmx3',xxmy3';yymx3',yymy3';xxmx4',xxmy4';yymx4',yymy4'];


clear xxmx5 yymx5 pm3
mx1 = max([m1(1,1),m2(1,1),m3(1,1),m4(1,1)]);
mx2 = min([m1(1,1),m2(1,1),m3(1,1),m4(1,1)]);
my1 = max([m1(1,2),m2(1,2),m3(1,2),m4(1,2)]);
my2 = min([m1(1,2),m2(1,2),m3(1,2),m4(1,2)]);
xxmx5 = xx1(xx1 < mx1 & xx1 > mx2);
yymx5 = yy1(yy1 < my1 & yy1 > my2);
pm3 = [];
for i = 1:size(xxmx5,2)
    for j = 1:size(yymx5,2)
        x0 = xxmx5(i);
        y0 = yymx5(j);
        m2m1m2p = (m1(1,1)-m2(1,1))*(y0-m2(1,2))-(x0-m2(1,1))*(m1(1,2)-m2(1,2));
        m3m4m3p = (m4(1,1)-m3(1,1))*(y0-m3(1,2))-(x0-m3(1,1))*(m4(1,2)-m3(1,2));
        m2m4m2p = (m4(1,1)-m2(1,1))*(y0-m2(1,2))-(x0-m2(1,1))*(m4(1,2)-m2(1,2));
        m3m1m3p = (m1(1,1)-m3(1,1))*(y0-m3(1,2))-(x0-m3(1,1))*(m1(1,2)-m3(1,2));
        if m2m1m2p*m3m4m3p>0 && m2m4m2p*m3m1m3p>0
            pm3 = [pm3;x0,y0];
        end
    end
end


ppm1 = [m1;m2;m3;m4;pm1;pm2;pm3];
ppm2 = unique(ppm1,'rows');
ppm2(ppm2(:,1) < xmin,:) = [];
ppm2(ppm2(:,1) > xmax,:) = [];
ppm2(ppm2(:,2) < ymin,:) = [];
ppm2(ppm2(:,2) > ymax,:) = [];
ppm = ppm2;
plot(ppm(:,1),ppm(:,2),'ro')
t= delaunayn(ppm);


clear x11 y11 x22 y22 x33 y33
for i = 1:size(t,1)
    x11 = ppm(t(i,1),1);y11 = ppm(t(i,1),2);
    x22 = ppm(t(i,2),1);y22 = ppm(t(i,2),2);
    x33 = ppm(t(i,3),1);y33 = ppm(t(i,3),2);
    xc = (x11 + x22 + x33)/3;
    yc = (y11 + y22 + y33)/3;
    for ix = 1:nex
        if ( xc > xx1(ix) ) && ( xc < xx1(ix+1) )
            iex = ix;
            break;
        end
    end
    for iy = 1:ney
        if ( yc > yy1(iy) ) && ( yc < yy1(iy+1) )
            iey = iy;
            break;
        end
    end
    cxh(i) = (iex-1)*ney + iey;
    ac = sqrt((x11-x22)^2+(y11-y22)^2);
    bc = sqrt((x22-x33)^2+(y22-y33)^2);
    cc = sqrt((x11-x33)^2+(y11-y33)^2);
    pc = (ac+bc+cc)/2;
    Sc(i) = sqrt(pc*(pc-ac)*(pc-bc)*(pc-cc));
end
C = (Sc./sum(Sc))';
CP0(ic) = {[cxh',C]};

%figure(2),trimesh(t,ppm(:,1),ppm(:,2))
%figure(1),hold on
for i = 1:size(t,1)
    h = fill([ppm(t(i,1),1),ppm(t(i,2),1),ppm(t(i,3),1)],...
         [ppm(t(i,1),2),ppm(t(i,2),2),ppm(t(i,3),2)],100*C(i,1));
       set(h,'edgealpha','0.2','facealpha',1)
       colormap('jet')
       colorbar
       hold on
end
pause(0.1)
end

for i = 1:size(CP0,2)
    CP1=CP0{1,1};
    B=unique(CP1(:,1));
    for j = 1:size(B,1)
        CP2 = CP1(CP1(:,1) == B(j,1),:);
        CP3=[B(j,1),sum(CP2(:,2))];
        CP4(j,:) = CP3;
    end
    CP(i) = {CP4};
end

胖射线追踪

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