MATLAB设计X射线对数曝光曲线

原创2023-04-17 11:07·卡洛db

X射线对数曝光曲线是射线检测工作前在一定条件下绘制的透照参数与透照厚度之间的关系曲线，检测人员可以从曲线中取得与检测工件相匹配的透照参数作为参考。对数曝光曲线需要通过实验制作，每台X射线机都有专属的对数曝光曲线不能通用。传统的对数曝光曲线制作和修正方法较繁琐，需在选定的透照条件下采用一系列不同的透照电压和不同的曝光量对阶梯试块进行多次透照，对得到的底片黑度进行测量，得到符合标准黑度要求的透照参数和厚度数据，通过人工在对数坐标纸上绘制曝光曲线，因而要得到高精度连续的曝光曲线，制作和使用难度加大。Matlab是一款很好的数据分析处理软件而且有很强的数据绘图功能，在人工智能和神经网络，建模演算等很多科研领域都有很好的应用,用来读取处理拟合数据绘制射线对数曝光曲线并从曲线上获取需要的透照参数很方便。

1 设计思路

1) 先对阶梯试块进行拍片拿到不同透照条件即不同曝光量和电压的二组六张底片的黑度、透照参数、透照厚度的数据；

2) 绘制不同透照参数的厚度/黑度图即T_D图，并从图中取得黑度D=2.5处的厚度值；

3) 利用从T_D图中取得的黑度D=2.5处的两组(六个)厚度值T和相对应的电压Kv值绘制Kv_T图，用matlab线性拟合插值方法，取出80Kv,100Kv,120Kv,140Kv,160Kv,180Kv,200Kv相对应的厚度值T,并填入表格中；

4) 用matlab命令xlsread读取表中统计的不同透照时间t，厚度值T，用指数函数t=a(1)*exp(a(2)*T）拟合出T_t的函数关系式，取得a(1)，a(2)系数值；

5) 用T=log(t/a(1))/a(2)，算出t=[0.5,1,1.5,2,3,4,5,7,10,15,20,30,40,60,100]对应的厚度值T,然后利用matlab一元二次拟合出厚度T和电压Kv=[80,100,120,140,160,180,200]的一元二次关系式Kv=aT2+bT+c，这样就可以得到系数[a,b,c]，然后解方程-aT2-bT-c+Kv=0就可以得到Kv=[50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200]对应的数值解厚度T ，最后就可以用plot(T,t)绘制t=[0.5,1,1.5,2,3,4,5,7,10,15,20,30,40,60,100]，Kv=[50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200]连续的对数曝光曲线。

2 完整代码

clear

figure(1)%绘制D-T图取D2.5处交点坐标

global pointTD

pointTD=[];

T1=xlsread('C:\Users\Desktop\RT-line-200.xlsx','D-T','A3:A12')

D1=xlsread('C:\Users\Desktop\RT-line-200.xlsx','D-T','B3:D12')

subplot(1,2,1)

plot(T1,D1(:,1))

hold on

plot(T1,D1(:,2))

hold on

plot(T1,D1(:,3))

grid on

grid minor

hold on

yline(2.5,'-r','D=2.5','LineWidth',0.5)

hold on

T2=xlsread('C:\Users\Desktop\RT-line-200.xlsx','D-T','R3:R12')

D2=xlsread('C:\Users\Desktop\RT-line-200.xlsx','D-T','S3:U12')

subplot(1,2,2)

plot(T2,D2(:,1))

hold on

plot(T2,D2(:,2))

hold on

plot(T2,D2(:,3))

grid on

grid minor

hold on

yline(2.5,'-r','D=2.5','LineWidth',0.5)


while(~waitforbuttonpress)

set(gcf,'WindowButtonDownFcn',@axesButtonDownFcn);

end


figure(2)%绘制Kv_T图并拟合读取绘制曝光曲线的厚度值

global pointTD

T3=pointTD(:,1);

Kv=[130,150,180];

plot(T3(1:3),Kv)

hold on

plot(T3(4:6),Kv)

axis([0 20,0,250])

grid on

grid minor

while(~waitforbuttonpress)%等待线性拟合，读取数据，填表完成，按回车继续

end


T1=xlsread(C:\Users\Desktop\RT-line-200.xlsx','E-T','B2:B15');

t1=xlsread('C:\Users\Desktop\RT-line-200.xlsx','E-T','C2:C15');

a=[];

for i=1:2:14

T=[T1(i),T1(1+i)];

t=[t1(1),t1(2)];

b=lsqcurvefit('exp_fit',[1 1],T,t);%cftool

a=[a,b];

TT=0:0.1:50;

tt=a(i)*exp(a(1+i)*TT);

figure(3)

plot(T,t,'o',TT,tt)

set(gca,'YScale','log');

hold on

grid on

end

t=[0.5,1,1.5,2,3,4,5,7,10,15,20,30,40,60,100];

Kv1=[80,100,120,140,160,180,200];

T=[];

for i=1:15

for j=1:2:13

Ta=log(t(i)/a(j))/a(1+j);

T=[T,Ta];

end

end

Kv_T=[];

for i=1:7:99

Kv_T1=polyfit(T(i:6+i),Kv1,2);

Kv_T=[Kv_T,Kv_T1];

end


ff=[]

% Kv2=[100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250]

Kv2=[50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200]

for i=1:16

for j=1:3:43

a=-Kv_T(j);

b=-Kv_T(1+j);

c=-Kv_T(1+j)+Kv2(i);

delta=b^2-4*a*c;

if delta<0;

display('no answers because delta smaller than 0')

continue;

elseif delta==0;

x1=-b/(2*a),x2=x1;

ff=[ff;x2];

continue;

else

x1=(-b+sqrt(delta))/(2*a),x2=(-b-sqrt(delta))/(2*a)

ff=[ff;x2];

continue;

end

end

end


for i=1:15:226

T=ff(i:14+i);

figure(4)

plot(T,t);

hold on

end

set(gca,'YScale','log');

title('对数曝光曲线','FontName','黑体')

xlabel('厚度 mm','FontName','黑体')

ylabel('时间 min','FontName','黑体')

axis([0 60,0,20])

grid on

grid minor

hold on;


while(~waitforbuttonpress)

set(gcf,'WindowButtonDownFcn',@axesButtonDownFcn);

end

% --- Executes on mouse press over figure background, over a disabled or

% --- inactive control, or over an axes background.

function axesButtonDownFcn(hObject, eventdata, handles)

global pointTD

pt = get(gca,'currentpoint');

T = pt(1,1);

D = pt(1,2);

hold on

plot(T, D,'*','color','r','markersize',8,'linewidth',0.5);

pointTD=[pointTD;[T,D]];

end

3 最终效果

举报
评论 0
    热门：