CUMCM-2017-problem-A程序

2017年高教社杯全国大学生数学建模竞赛题目

（请先阅读“全国大学生数学建模竞赛论文格式规范”）

A题CT系统参数标定及成像

CT(Computed Tomography)可以在不破坏样品的情况下，利用样品对射线能量的吸收特性对生物组织和工程材料的样品进行断层成像，由此获取样品内部的结构信息。一种典型的二维CT系统如图1所示，平行入射的X射线垂直于探测器平面，每个探测器单元看成一个接收点，且等距排列。X射线的发射器和探测器相对位置固定不变，整个发射-接收系统绕某固定的旋转中心逆时针旋转180次。对每一个X射线方向，在具有512个等距单元的探测器上测量经位置固定不动的二维待检测介质吸收衰减后的射线能量，并经过增益等处理后得到180组接收信息。

CT系统安装时往往存在误差，从而影响成像质量，因此需要对安装好的CT系统进行参数标定，即借助于已知结构的样品（称为模板）标定CT系统的参数，并据此对未知结构的样品进行成像。

请建立相应的数学模型和算法，解决以下问题：

(1) 在正方形托盘上放置两个均匀固体介质组成的标定模板，模板的几何信息如图2所示，相应的数据文件见附件1，其中每一点的数值反映了该点的吸收强度，这里称为“吸收率”。对应于该模板的接收信息见附件2。请根据这一模板及其接收信息，确定CT系统旋转中心在正方形托盘中的位置、探测器单元之间的距离以及该CT系统使用的X射线的180个方向。

(2) 附件3是利用上述CT系统得到的某未知介质的接收信息。利用(1)中得到的标定参数，确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。另外，请具体给出图3所给的10个位置处的吸收率，相应的数据文件见附件4。

(3) 附件5是利用上述CT系统得到的另一个未知介质的接收信息。利用(1)中得到的标定参数，给出该未知介质的相关信息。另外，请具体给出图3所给的10个位置处的吸收率。

(4) 分析(1)中参数标定的精度和稳定性。在此基础上自行设计新模板、建立对应的标定模型，以改进标定精度和稳定性，并说明理由。

(1)-(4)中的所有数值结果均保留4位小数。同时提供(2)和(3)重建得到的介质吸收率的数据文件（大小为256×256，格式同附件1，文件名分别为problem2.xls和problem3.xls）

图1. CT系统示意图图2. 模板示意图（单位：mm）图3. 10个位置示意图

程序分析：

function [x0,y0,sepdist,start_angle,direction180,gain] = problem1_20171010(input_args)

%CT系统旋转中心在正方形托盘中的位置(x0,y0),

%探测器单元间距sepdist,

%start_angle:CT系统使用的X射线的180个投影方向起始方向

%CT系统使用的X射线的180个投影方向direction180

%gain:CT系统的增益

clc;

clear all;

close all;

format short;

%模板:template

template=xlsread('a题附件.xls','附件1');

%图象中的灰度值越大,显示越亮

figure;imshow(template,[]);

title('模板');

%1.计算公切线的斜率

%使用优化的方法求模板中椭圆与圆的4条公切线的斜率,交点坐标,

%外面的上方与下方的，交叉的上方与下方的公切线

%椭圆的方程：x^2/225+y^2/1600=1,圆的方程：(x-45)^2+y^2=16,

%设公切线为y=kx+b,他与圆心(45,0)之间的距离等于4,得到方程:2009*k^2+90*k*b+b^2-16=0

%与椭圆相交的两个交点的x坐标相同,即将y=kx+b代入椭圆方程后,判别式=0,得到方程

%129600*k^2-576*b^2+921600=0 我是算的%225*k^2-b^2+1600=0;

%解这个方程组得到k,b

%初始搜索值[-1,41]',Make a starting guess at the solution

%options=optimset('Display','iter'); % Option to display output

%[kb1,fval] = fsolve(@(kb)optkb(kb),[-1,41]',options);% Call solver

% X = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON)

%[kb1,fval]=fmincon(@(kb)optkb2(kb),[-1,41]',[],[],[],[],[-1,40]',[-0.2,50]',@(kb)nonlcon1(kb)); %kb1=[-0.4257,40.5064];

%theta1=atan(kb1(1))*180/pi;%=-23.0576

%计算所求的公切线与圆的交点,以判断所求的切线是哪一条切线

%[xyc1,fval] = fsolve(@(xyc)optxyc(xyc,kb1),[45,4]',options);

%xyc1=[46.6678,3.9182]';%由此判断所求的外公切线是外公切线

%kb1(1)是椭圆上方与圆的上方的外公切的斜率

% X = fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON)

[kb2,fval]=fmincon(@(kb)optkb2(kb),[-1,41]',[],[],[],[],[-1.1,41]',[-0.5,46]',@(kb)nonlcon2(kb)); %kb2=[-1.0733,43.1184];

theta2=atan(kb2(1))*180/pi;%=-47.0249

%计算所求的公切线与圆的交点,以判断所求的切线是哪一条切线

%[xyc2,fval] = fsolve(@(xyc)optxyc(xyc,kb2),[45,4]',options);

%xyc2=[41.6595,-2.2002]';%由此判断所求的外公切线是内公切线

%kb2(1)是椭圆上方与圆的下方的内公切的斜率

%ct_radon为该CT系统对template_img进行radon变换之后的投影图

ct_radon=xlsread('a题附件.xls','附件2');

%figure;imshow(ct_img);

%imshow(I,[ ])：这种方式就是把imshow的显示范围设置成[min(I(:)) max(I(:))]，

%也就是线性映射，相当于imagesc(I),colormap(gray(256))可以将整幅图像的信息显示出来。figure;imshow(ct_radon,[]);

title('旋转中心有偏移的模板的CT图');

%画模板的CT图中各个投影方向的投影

figure;

interval=zeros(512,180);

for i=1:180

plot(ct_radon(:,i));

grid on;

hold on;

pos=find(ct_radon(:,i)~=0);%找到不为零的所有元素的下标

for j=1:length(pos)-1

%interval(pos(j),i):ct_radon中的第i列ct_radon(:,i)中有投影的的位置差

%如果某个interval(pos(j),i)>1,则说明投影在这个位置之前有0投影,即投影间断了

interval(pos(j),i)=pos(j+1)-pos(j);

end% for j=1:length(pos)-1

end %for i=1:180

xlabel('投影地址');

ylabel('投影值');

title('旋转中心有偏移的模板的CT图中各个方向的投影');

saveas(gcf,'有偏移的模板的CT图中各个方向的投影','jpg');

%2.计算公切线投影方向在ct_radon中对应的投影线号

%计算CT系统相邻方向的角度差

%原理：模板的CT图的第14与第108条投影线对称且正好是第1条与最后1条

%不能将椭圆的投影与圆的投影分离的两个投影方向，

%因此，该两条投影线的投影方向就是椭圆与圆的上下两条内公切线的方向，

%这两个投影方向所夹的角度是两条内公切线所夹的角度，即是2*abs(theta2)

%将这个角度值除以(108-14)的结果就是ct系统相邻方向的角度差delta_angle

%通过interval找到第1条不能将椭圆与圆的投影分开的投影线号,

%interval(:,i)中各个不为0的元素全部是1,则该投影方向是不能将椭圆与圆的投影分开的方向

%interval(:,i)中各个不为0的元素中有大于1的,则该投影方向是能将椭圆与圆的投影分开的方向

%原理:CT系统的起始投影方向角度start_angle=第1条公切线的角度-(第1条不能将椭圆与圆的投影分开的投影线号-1)*delta_angle

k=0;

for i=1:180

if sum(interval(:,i))==length(find(interval(:,i)~=0))

k=k+1;

cprojection(k)=i;

end

end

%108;

radioe=cprojection(end);

%14

radios=cprojection(1);

figure;

plot([1:512],ct_radon(:,radios),'r',[1:512],ct_radon(:,radioe),'b');

grid on;

xlabel('投影地址');

ylabel('投影值');

legend('第1条内公切线方向的投影','第2条内公切线方向的投影');

%title(sprintf('旋转中心有偏的模板的CT图中第%d与%d条投影',radios,radioe));

%saveas(gcf,sprintf('旋转中心有偏的模板的CT图中第%d与%d条投影',radios,radioe),'jpg'); title('旋转中心有偏移的模板的CT图中椭圆与圆的两条内公切线方向的投影');

saveas(gcf,'有偏移的模板的CT图中椭圆与圆的两条内公切线方向的投影','jpg');

%3.CT系统相邻投影方向的角度差

%1.010

delta_angle=2*abs(theta2)/(radioe-radios);

%4.CT系统的起始投影方向角度

%CT系统的起始投影方向角度start_angle=

%第1条公切线的角度-(第1条不能将椭圆与圆的投影分开的投影线号-1)*delta_angle

%-60.5982

start_angle=theta2-(radios-1)*delta_angle;

%第1条投影线方向对应的投影检测器与水平方向的夹角

%29.4018

horizen_start_angle=90+start_angle;

xlswrite('CT系统使用的X射线的起始方向与相邻方向角度差.xls',[start_angle,delta_angle,horizen_start_angle]);

%将start_angle转换成正角度

%start_angle=start_angle+180;

%5.CT系统使用的X射线的180个投影方向

direction180=start_angle:delta_angle:(start_angle+179*delta_angle);

xlswrite('CT系统使用的X射线的180个方向.xls',direction180');

for i=1:18

direction180paper(2*i-1,:)=(i-1)*10+1:i*10;

direction180paper(2*i,:)=direction180((i-1)*10+1:i*10);

end

xlswrite('CT系统使用的X射线的180个方向paper.xls',direction180paper);

%6.确定探测器单元间距sepdist

%原理:当投影方向为0度时,这时圆的投影重叠在椭圆的投影上,

%投影的宽度为80mm，则sepdist=80/（这时的投影线数-1）

%投影线数的计算原理：这时该投影方向是不能将椭圆与圆的投影分开的方向，

%查找interval中与0度方向对应的投影线号，由于投影没有分离，

%则与该投影线号对应的interval中的所有非零元素全是1，

%因此，与该投影线号对应的interval中的所有非0元素的和就为这时的投影线数projnum=sum(interval(:,find (abs(direction180)

%0.2778mm

sepdist=80/projnum(1);

%探测器总长度

%141.9mm

probe_length=(512-1)*sepdist;

xlswrite('CT系统的探测器单元间距与总长.xls',[sepdist,probe_length]);

%7.模板的理想CT图

% radon 中缺省的值=2*ceil(norm(size(num1)-floor((size(num1)-1)/2-1))+3

%如果指定了N,则R将包含N个行向量

% R = radon(I,THETA) returns the Radon transform of the

% intensity image I for the angle THETA degrees. If THETA is a

% scalar, the result R is a column vector containing the Radon

% transform for THETA degrees. If THETA is a vector, then R is

% a matrix in which each column is the Radon transform for one

% of the angles in THETA. If you omit THETA, it defaults to 0:179.

% R的每一列对应一个角度,第1至第N行对应投影位置

theta=0:179;

N=512;

[template_radon,xp] = radon(template,theta,N);

figure,imshow(template_radon,[]);

title('旋转中心没有偏移的模板的CT图');

imwrite(uint8(template_radon),'没有偏移的模板的CT图.jpg','jpg');

xlswrite('没有偏移的模板的CT数据',template_radon);

%画无偏移时0,90度时的投影值曲线,横坐标是假设的投影地址

%无偏移时0度的投影值曲线,无偏移时90度的投影值曲线

figure;

plot([1:N],template_radon(:,1),'b',[1:N],template_radon(:,91),'r');

grid on;

xlabel('投影地址');

ylabel('投影值');

legend('无偏移0度投影曲线','无偏移的90度投影曲线');

title('旋转中心无偏移时0度与90度的投影曲线');

axis([0 512 0 max(max(template_radon))+10]) ;

saveas(gcf,'无偏移时0度与90度的投影曲线','jpg');

%8.CT系统旋转中心在正方形托盘中的位置(x0,y0)

%以正方形托盘的中心建立直角坐标系:原点o在正方形中心,x轴正向指右,y轴正向向上%原理:在旋转中心有偏移的ct_radon中当投影方向为90度时的投影与

%旋转中心没有偏移的template_radon中的90度方向的投影

%(即template_radon中的第1列投影)相似,由这两个投影确定水平方向的偏移x0;

%即由template_radon90与ct_radon90确定水平方向的偏移x0;

%在旋转中心有偏移的ct_radon中当投影方向为0度时的投影与

%旋转中心没有偏移的template_radon中的180度(即0度)方向的投影

%(即template_radon中的第91列投影)相似,由这两个投影确定垂直方向的偏移y0;

%即由template_radon0与ct_radon0确定垂直方向的偏移y0;

%template_radon90:template_radon中的90度方向的投影(垂直投影)

template_radon90=template_radon(:,1);

%template_radon0:template_radon中的180度(即0度)方向的投影(水平投影)

template_radon0=template_radon(:,91);

%ct_radon90:ct_radon中的90度方向的投影(垂直投影)

%附件1中所给CT数据文件中大于0度方向时探测器的方向和radon中的探测方向相反, %因此需要统一探测器方向，使ct_radon上下翻转过来就可以做到

%使附件1中的数据与radon变换的数据中探测器的方向一致

%附件1中所给CT数据文件中小于等于0度方向时探测器的方向

%和radon中的探测方向相同,不需要翻转

%ct_col90:ct_radon中90度方向对应的投影在ct_radon中的列序号

%direction180的起始角度为:-60.4886

ct_col90=find(abs(direction180-90)

ct_radon90=flipud(ct_radon(:,ct_col90));%后面用来计算增益系数

%ct_radon90=flipud(ct_radon(:,find(abs(direction180-90)

%ct_radon0:ct_radon中的0度方向的投影(水平投影)

%ct_col0:ct_radon中0度方向对应的投影在ct_radon中的列序号

ct_col0=find(abs(direction180)

ct_radon0=ct_radon(:,ct_col0);

%ct_radon0=ct_radon(:,find(abs(direction180)

%画图比较旋转中心无偏移与有偏移时候投影方向分别为水平与垂直方向的两条投影曲线figure;

plot([1:N],template_radon0,'b',[1:N],ct_radon0,'r');

grid on;

xlabel('投影地址');

ylabel('投影值');

legend(sprintf('无偏移的%d度水平方向投影',0),sprintf('有偏移的%d度水平方向投影',0)); title(sprintf('旋转中心有与无偏移的%d度水平方向投影',0));

%axis([0 512 0 max(max(template_radon))+10]) ;

saveas(gcf,'旋转中心有与无偏移的水平方向的投影','jpg');

figure;

plot([1:N],template_radon90,'b',[1:N],ct_radon90,'r');

grid on;

xlabel('投影地址');

ylabel('投影值');

legend(sprintf('无偏移的%d度垂直方向投影',90),sprintf('有偏移的%d度垂直方向投影',90)); title(sprintf('旋转中心有与无偏移的%d度垂直方向的投影',90));

%axis([0 512 0 max(max(template_radon))+10]) ;

saveas(gcf,'旋转中心有与无偏移的垂直方向的投影','jpg');

%由template_radon90与ct_radon90确定水平方向的偏移x0的原理:

%计算这两条投影的重心,两个重心的差为水平方向的偏移x0,

%详细分析见参考文献:

%李增云,吕东辉,利用部分投影的二维CT旋转中心偏移快速校正[J],CT理论与应用研究,2015,24(4):533-543.

%即:物体质心C在探测器上的投影必与投影函数的质心C_p重合

%计算水平方向的偏移x0

wtemplate90=sum(template_radon90);

wtemplatex=[1:N]*template_radon90;

wct90=sum(ct_radon90);

wctx=[1:N]*ct_radon90;

x0=(wtemplatex/wtemplate90-wctx/wct90)*sepdist;%重心这样算-x0

%或注意：旋转中心没有偏移时正方形模板的重心(template_x0,template_y0)

%可以通过二重积分的计算直接计算出来

%w:模板质量,设各处的密度为1

%w=pi*15*40+pi*4^2;

%template_x0=(int_D1 xdxdy+int_D2 xdxdy)/w=(int_0^2*pi int_0^1 15*600* r^2 *cos(theta) dr d theta

%+int_0^2*pi int_0^4 (45+r* cos(theta))*r dr d theta)/w=720*pi

%template_x0=720/616;

%或者转换为投影地址：

%template_x0_projnum=256.5+template_x0/sepdist;

%根据对称性，正方形模板的重心的纵坐标template_y0=0（单位：mm）,或转换为投影地址为:256.50

%template_y0=0;

%计算垂直方向的偏移y0

wtemplate0=sum(template_radon0);

wtemplatey=[1:N]*template_radon0;

wct0=sum(ct_radon0(:,1));

wcty=[1:N]*ct_radon0(:,1);

y0=(wtemplatey/wtemplate0-wcty/wct0)*sepdist;%-y0

%或者根据图旋转中心有与无偏移的0度水平方向的投影图可以看出,

%有偏移的0度水平方向投影相对无偏移的0度水平方向投影向左平移了,

%这个平移量就是旋转中心相对模板正方形中心向上平移的量

p1=find(template_radon0~=0);

p2=find(ct_radon0(:,1)~=0);

y0=(p1(1)-p2(1))*sepdist;

%或者

y0=(p1(end)-p2(end)+p1(1)-p2(1))*sepdist/2;

%或者根据图旋转中心有与无偏移的90度垂直方向的投影图可以看出,

%有偏移的90度垂直方向投影相对无偏移的90度垂直方向投影向右平移了,

%这个平移量就是旋转中心相对模板正方形中心向左平移的量

p1=find(template_radon90~=0);

p2=find(ct_radon90~=0);

x0=(p1(1)-p2(1))*sepdist;

%或者

x0=(p1(end)-p2(end)+p1(1)-p2(1))*sepdist/2;

xlswrite('CT系统旋转中心在正方形托盘中的位置.xls',[x0,y0]);

%9.计算增益:gain

%正方形模板中的每一点的数值反映了该点的吸收强度，称为吸收率,

%二维待检介质吸收衰减后的射线能量，并经过增益处理后得到180组接收信息

%当强度为I_0的窄束射线穿过吸收系数为nu 的物体时,其强度满足指数摔减关系:

%I=I_0*exp(-nu*t),式中t为射线所穿过物质层厚度.在实际情况中，所研究的物质

%往往不是由单一成分组成的,当物体由若干个不同成分组成时,物体内部各处的nu

%也可能不同.在这样的物质中，射线束穿过这个物件后的强度为

%I(L)=I_0*exp(-int_L u(t)dt),式中u(t)为t处的吸收率.

%CT系统通过改变一组射线路径L,记录下对应出射强度I(L)的变化来分析物体内部u(t)的分布.

%在这里可以假设入射线束强度都为1,则强度衰减模型为:I=gain*exp(-nu*t),gain为增益系数%根据template可以知道这里的模板物质各处的吸收系数nu=1,或0

%所以强度衰减模型为:I=gain*exp(-t)

%只用90度方向的一条投影和模板中的数据来进行最小二乘拟合gain

%y:与ct_radon90中各个投影地址对应的模板在竖直方向的吸收率不为0的介质长度

%分析旋转中心有与无偏移的模板在竖直方向的投影,

%可以知道经过radon变换之后接收信息rxv与X射线经过介质的长度y成正比,比例系数为增益系数gain

%根据附件1中模板介质的吸收率为1,可以计算系统的放大增益

%ct_radon90~=0的位置对应为:237-344,439-467

%y=sum(template);不等于0的位置对应:1-166,234-253

p1=find(ct_radon90~=0);

p2=p1(2:end);

%p=find(p2-p1(1:end-1)~=1);

%rxv是与模板中椭圆部分相对应的经过增益后X射线的吸收值

rxv=ct_radon90(p1(1:find(p2-p1(1:end-1)~=1)));%108个%?

%计算介质的长度y，计算介质在椭圆-15<=x<=15的长度y

x=-15:30/(length(rxv)-1):15;

%椭圆方程:x^2/225+y^2/1600=1

y=2*40*sqrt((1-x.^2/225));

%去掉y的两端长度比较短或者为0的元素,对应的rxv一样操作

y=y(10:length(y)-10);

rxv=rxv(10:length(rxv)-10);

gain=rxv(10)/y(10);

[gain,res]=lsqnonlin(@(gain)myfun(gain,rxv,y'),gain);

xlswrite('CT系统的增益系数.xls',gain);

%保存CT系统的所有标定参数

a{1,1}='旋转中心在托盘中的位置(x0,y0):';

a{1,2}=x0;

a{1,3}=y0;

a{2,1}='探测器单元间距(单位:mm):';

a{2,2}=sepdist;

a{3,1}='探测器长度(单位:mm):';

a{3,2}=sepdist*511;

a{4,1}='CT系统的X射线的起始投影方向:';

a{4,2}=start_angle;

a{5,1}='相邻两条X射线方向之间的角度差:';

a{5,2}=delta_angle;

%起始位置探测器与水平方向的夹角

a{6,1}='起始位置探测器与水平方向的夹角:';

a{6,2}=horizen_start_angle;

a{7,1}='CT系统的增益系数:';

a{7,2}=gain;

a{8,1}='CT系统的X射线的180个投影方向:';

for i=8:43

for j=2:11

a{i,j}=direction180paper(i-7,j-1);

end

end

xlswrite('CT系统标定的参数.xls',a);

end

%====================================

function F=optkb(x)

%x=[k,b]

F(1,1)=2009*x(1)^2 +90*x(1)* x(2)+x(2)^2 -16;

F(2,1)=129600*x(1)^2 -576*x(2)^2+921600;

end

%====================================

function F=optkb2(x)

%x=[k,b]

F=abs(2009*x(1)^2 +90*x(1)* x(2)+x(2)^2 -16)+abs(129600*x(1)^2 -576*x(2)^2+921600); end

%====================================

function F=optxyc(x,kb)

%x=[x,y]

F(1,1)=(x(1)-45)^2 +x(2)^2 -16;

F(2,1)=kb(1)*x(1) -x(2)+kb(2);

end

%====================================

function fval= myopt(b)

k=(-90*b-sqrt(8100*b^2-8036*(b^2-16)))/4018;

%y=k*x+b;代入椭圆方程,使判别式大于等于零,并尽可能接近于0 %(64+9*k^2)*x^2+18*k*b*x+9*b^2-14400=0

fval=324*k^2*b^2-36*(64+9*k^2)*(b^2-1600);

end

%====================================

function [c,ceq] = nonlcon(b)

k=(-90*b-sqrt(8100*b^2-8036*(b^2-16)))/4018;

delta1=324*k^2*b^2-36*(64+9*k^2)*(b^2-1600);

c =-delta1;

ceq = [];

end

%====================================

function [c,ceq] = nonlcon1(x)

%y=45k+b>0

c=-45*x(1)-x(2);

ceq = [];

end

%====================================

function [c,ceq] = nonlcon2(x)

c=45*x(1)+x(2);

ceq = [];

end

%====================================

function F=myfun(gain,rxv,y)

%rxv=gain*(1-exp(-y))

% F=sum((rxv-gain*(1-exp(-y))).^2);

%rxv=gain*y

F=sum((rxv-gain*y).^2);

end