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）
整体求解过程概述(摘要)
本文针对CT系统参数标定及成像问题，建立了基于最小二乘法的CT系统标定模型，求解了CT系统的旋转中心坐标、探测器单元间距以及X射线的180个方向；建立了基于傅立叶变换的图像重建模型，求解了未知介质的吸收率分布；建立了基于滤波反投影法的图像重建模型，求解了噪声较大的未知介质的吸收率分布；设计了新型标定模板并建立了对应的标定模型，解决了改进标定精度和稳定性的问题。
针对问题一，建立了基于最小二乘法的CT系统标定模型，求解了CT系统的标定参数。首先，通过最小二乘法，将180组接收信息拟合成正弦函数，根据几何关系联立方程组求解出180个旋转角度，以及探测器单元间距𝑑=0.2758。其次，通过最小二乘优化、变步长遍历搜索的方法，求得旋转中心为(-9.2513,5.6783)。最后，基于所得标定参数及已知模板的吸收率分布，求解标定模板的接收信息，进行残差分析，检验了标定模型的合理性。
针对问题二，建立了基于傅立叶变换的图像重建模型，求解了未知介质的吸收率分布。首先，依据 Radon 变换及中心切片定理，对投影数据进行一维傅立叶变换，得到定义在傅立叶空间的极坐标网格。其次，在频域上对投影进行积分，再次，直接通过二维傅立叶反变换得到重建图像，其中10处吸收率分别为0.0000、0.0000、1.0060、0.0000、1.0081、0.0081、0.0000、0.1008、0.0081、0.0081。最后，对标定参数进行了灵敏度分析。
针对问题三，建立了基于滤波反投影法的图像重建模型，求解了噪声较大的未知介质的吸收率分布。首先，采用问题二的方法，发现重建图像噪声较大，改用滤波反投影算法去除重建图像的伪影。再次，选择中值滤波器进行滤波卷积，采用线性内插的方法进行插值，反投影累加后得到重建图像，其中 10处吸收率分别为 2.0323、1.2863、2.0081、1.6370、2.0484、2.0081、0.0000、2.0343、1.5363、0.0403。最后，对标定参数进行了灵敏度分析。
针对问题四，设计了新型标定模板并建立了对应的标定模型，解决了改进标定精度和稳定性的问题。首先，对原有模板进行误差分析，得出对椭圆边缘处接收信息的插值拟合会导致旋转角度出现误差的结论。其次，为避免类似误差，设计圆形组合模板，并建立对应标定模型。再次，对比基于新旧模板的不同标定参数的图像重建模型，求得新旧模板的峰值信噪比分别为50.3614、52.2537。最后，求解基于不同旋转角度个数的归一化均方距离，新模板重建图像精度趋于稳定的旋转次数小于旧模板，进而验证了基于新模板的标定模型精度更高，更稳定。
模型假设：
1. 假设射线源中点与转台旋转中心连线垂直于探测器阵列；
2. 假设射线源的衰减仅与被介质吸收有关，并未产生散射；
3. 假设增益处理是线性增大的。
问题分析：
问题整体思路较为清晰，层层递进。问题一是根据已知标定模板的吸收率及其接收信息对CT系统进行正向标定的问题。问题二、三建立在问题一求解的标定系数基础之上，是依据不同的接受信息对未知介质的吸收率分布进行反向重建的问题。分析问题二、三的接受信息可发现问题三图像噪声偏大，故对问题三的求解大致思路与问题二类似，但需进行进一步的降噪处理。问题四需对原有标定模型进行精度和稳定性分析，并设计新的模板及标定模型，进而提高精度及稳定性，属于优化问题。
问题一的分析
首先，分析题目条件可知，均匀模板的吸收率分布及对应的接受信息已知，且模板质地均匀，故射线衰减公式中的对吸收率分布进行的线积分可简单变为介质厚度，故可对接收信息的表达式进行求解。 其次，要将离散的接受信息连续化，可考虑通过最小二乘法对其进行拟合， 初步分析为两个近似正弦曲线，分别对应图中的椭圆和圆，其中正弦曲线与X轴相交的情况反应了平行光束与椭圆相切的两条切线距离的大小。由于两切线之间的距离与旋转中心[1]、旋转角度存在一定的几何关系，列出关系式，联立方程可解出各个旋转角度的值。 再次，基于上述拟合曲线及几何关系式，可求得多个探测器单元间距及旋转中心的坐标。针对探测器单元间距的求解，出于对使用信息完整性的考虑，可通过多次求解方程组，并对多个解取平均得到最终答案。针对旋转中心的求解，由于探测器中垂线衡过旋转中心，且旋转角度已求解得到，故可列出不同角度下中垂线的直线方程，两两联立即可求出旋转中心。出于保留信息完整性，同时避免算法过于复杂的考虑，故均匀抽取30组直线，得到元素个数为345的点集。可通过最小二乘优化[2]，建立各点到中心点距离平方和最小的目标函数，可通过变步长遍历搜索的方法，求解旋转中心。 最后，由于接受信息的表达式可被求解，且介质吸收率分布已知，故可对已知模板的接受信息进行模拟，可与题中所给条件进行比对，进行残差分析，即可检验模型的合理性。
问题二的分析
首先，分析题目条件可知，由于在第一问已对CT系统进行了标定，故可根据题中所给的接收信息，对图像进行重建。对图像重建的常用方法有傅里叶变换法、滤波反投影法，由于部分接收信息的拟合后的图像较为平滑，故采用算法较为简单的傅里叶变换法。 其次，在推导非均匀介质中衰减公式时，可对非均匀物体进行微元化考虑，进而对射线贯穿部分的吸收率分布进行线积分。可根据Radon变换[3]及中心切片定理[4]，得到多个角度下的投影后，通过投影反变换可求出被积函数即线衰减系数函数，从而得到对应的线衰减系数分布的图像，重建出问题二中非均匀物体的形状及在正方形托盘上的位置。 再次，为得到题中所要求的10个点的吸收率，需要得到吸收率与图像灰度值之间的关系。由于问题一中模板条件已知，运用问题二中的算法可重建出问题一中的模板图像，得到已知模板图像与重建模板图像之间的灰度值关系后再据此关系得问题三图像的10个吸收率。 最后，可针对三个不同标定参数，进行固定比率的修改，并计算参数修改后重建图像吸收率分布的变化率，可据此对各标定参数进行灵敏度分析。
问题三的分析
首先，问题三与问题二提供的条件并无二致，然而通过观察附件中的接收信息值，可发现数据在表格中的分布更加复杂，很可能导致运用问题二中的傅立叶变换算法重建出的图像质量不佳，故在问题三考虑采用更优的算法重建图像。通过查阅文献可知，滤波反投影算法[5]可大大消除反投影重建算法的模糊现象，保证重建图像的边缘清晰和内部分布均匀，得到质量较好的重建图像。故可改用滤波反投影算法进行图像重建。 其次，滤波反投影算法的关键在于滤波函数的选取，在众多滤波函数中，中值滤波[6]对噪声的除噪作用较显著，能够较好地保护图像边缘，降低噪声。根据反投影算法，得到滤波后的投影数据便可对图像进行重建，但平行束投影的采集过程中，射线发射器和接收器是离散的，与此同时角度采样也是离散化的，故需要选用合适的插值方法对滤波后的投影数据进行处理。通过查阅文献可知，常用方法有线性内插和紧邻内插，线性内插更优，故可采用线性内插的方法。由此，可对投影数据进行反变换，进而得重建后的图像及吸收率分布。 最后，可针对三个不同标定参数，进行固定比率的修改，并计算参数修改后重建图像吸收率分布的变化率，可据此对各标定参数进行灵敏度分析。
问题四的分析
首先，分析题目可知，应先对原有模板进行误差分析，进而考虑新模板的设计方向。由于问题一中对椭圆边缘处的接收信息进行拟合时存在一定误差，故造成对射线旋转方向的求解产生了误差。为避免这一误差，可考虑设计组合圆模板对上述误差进行规避，通过设计特殊的圆心及半径，避免拟合，进而得到优化后的旋转角度。 其次，要对比基于新旧模板标定模型的优劣，可通过比对基于不同标定参数的图像重建模型的优劣来体现。根据题意，从精度及稳定性[7]两方面来考虑。 针对精度，可建立重建图像评价指标。一般指标有归一化均方距离、归一化平均绝对距离、峰值信噪比，由于三个指标互成正比或反比，故无需对各指标进行赋权，针对基于新旧标定参数的图像重建模型，可通过重建原有模板图形，计算各指标值，比对新模型的经度是否更高。
针对稳定性，由于重建模型的重建质量是关于旋转角度个数变化而变化的针对这一动态过程，模型的稳定性刻画了重建图像精度随旋转角度个数增加的变化趋势。故可通过均匀减少图像旋转角度个数，求解新旧图像重建模型的归一化均方距离。对上述散点图进行插值拟合，即可得到新旧图像重建模型的图像精度变化趋势，比对新旧模型的稳定点，即可比对新旧模型的稳定性。 最后，依据精度和稳定性两个评价指标，即可判定新模板及其对应标定模型的合理性。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：
%最小二乘法拟合正弦曲线MATLAB代码
function [aa,zd] = q1(xdata,ydata)
% ydata=xlsread('C:\Users\hp\Desktop\A.xls','a','C1:C29');
% ydata=ydata';
%估计正弦曲线的初值
A=max(ydata);
T=length(ydata);%假设用的是一个周期的正弦波拟合，只是猜测值
w=2*pi/T;
c=pi;
d=mean(ydata);
options = optimset('MaxFunEvals',1000000); %保证精度
x_est = lsqcurvefit(@(x,xdata) myfun(x,xdata),[A w c d],xdata,ydata);% 利用
lsqcurvefit 函数拟合
plot(xdata,myfun(xdata,x_est))%画出拟合图
aa=2*pi/w;%计算拟合后的周期
zd=A+d;%计算拟合后的最高点
%180 个角度计算MATLAB代码
function [aa] = q1_1(~)
aa=1;
data=xlsread('C:\Users\hp\Desktop\A.xls','附件 2');%导入数据
for i=109:180
ydata=data(:,i);%初始化处理矩阵
%以下为初始化参数值
flag=0;
s1=0;
s2=0;
ydata1=0;
ydata2=0;
%以下将图中两个椭圆分开，把分开后的非零值分别放入两个矩阵
for j=1:512
if j==1
if ydata(j,1)~=0
flag=1;
end
else
if (ydata(j,1)~=0)&&(ydata(j-1,1)==0)&&(flag~=2)
flag=1;
end
if (ydata(j,1)==0)&&(ydata(j-1,1)~=0)
flag=2;
end
if (ydata(j,1)~=0)&&(ydata(j-1,1)==0)&&(flag==2)
flag=3;
end
end
if flag==1
s1=s1+1;
ydata1(1,s1)=ydata(j,1);
end
if flag==3
s2=s2+1;
ydata2(1,s2)=ydata(j,1);
end
end
xdata1=1:length(ydata1);
xdata2=1:length(ydata2);
[t1,zd1]=q1(xdata1,ydata1);%拟合其中一个椭圆曲线
[t2,zd2]=q1(xdata2,ydata2);%拟合另一个椭圆曲线
d=4*t2/t1;%计算距离
%以下两个式子为计算直线与椭圆的交点
y= -((81920*(d - 15)*(d + 15))/11)^(1/2)/(2*d);
x=(-(1620*(d - 40)*(d + 40))/11)^(1/2)/(2*d);
k=(x*1600)/(y*225);%计算直线的斜率
oo=i-108;
cita(1,oo)=atan(k)*180/pi;%计算出 cita 角
zd(1,oo)=zd2;%计算出最高点
end
xlswrite('C:\Users\hp\Desktop\2.xls',cita);%写入 Excel 文档
xlswrite('C:\Users\hp\Desktop\zd.xls',zd);
%拟合正弦函数模型MATLAB代码
function F = myfun(x,xdata)
F = x(1)*sin(x(2)*xdata+x(3))+x(4) ;
%求旋转中心MATLAB代码
clear
A=xlsread('C:\Users\hp\Desktop\一.xls');%导入数据
for i=1:84
s(1,i)=sin(A(2,i));%初始化 x 前系数
end
for i=1:84
c(1,i)=cos(A(2,i));%初始化 y 前系数
end
for i=1:84
n(1,i)=(A(1,i)-256.5)*0.2758;%初始化常系数
end
count=1;
%分别对两个条直线求交点
for i=1:84
for j=i:84
if i~=j
%
i));
i));
end
[x,y]=solve('s(1,i)*x+c(1,i)*y+n(1,i)=0','s(1,j)*x+c(1,j)*y+n(1,j)=0');
z(1,count)=-(c(1, i)*n(1, j) - c(1, j)*n(1, i))/(c(1, i)*s(1, j) - c(1, j)*s(1,
z(2,count)=-(n(1, i)*s(1, j) - n(1, j)*s(1, i))/(c(1, i)*s(1, j) - c(1, j)*s(1,
count=count+1;
end
end
%以下为计算求得交点集的边界值
xmin=min(z(1,:));
xmax=max(z(1,:));
ymin=min(z(2,:));
ymax=max(z(2,:));
abc=100000;
for i=xmin:0.0001:xmax
for j=ymin:0.0001:ymax
mm=0;
for o=1:3486
d=(z(1,o)-i)^2+(z(2,o)-j)^2;%遍历的目标函数
mm=mm+d;
end
if mm<abc
x=i;
y=j;
abc=mm;
end
end
end
i
j
%残差分析MATLAB代码
clear
shiyan=xlsread('C:\Users\hp\Desktop\问题一.xls','1');%导入实验值
shiji=xlsread('C:\Users\hp\Desktop\问题一.xls','2');%导入实际值
for i=1:233
s=0;
for j=1:235
d=abs(shiyan(j,i)-shiji(j,i));%求实际值与实验值的差值
s=s+d;%将每一列的残差求和
end
cancha(1,i)=s/233;%求每一列的残差均值并放入矩阵
end
xp=1:233;
x=1:256;
f=polyfit(xp,cancha,1);%把残差值进行拟合
canchan=polyval(f,x);%算出拟合值
plot(x,canchan,'k');%画出拟合曲线
for i=1:256
if i<=50
cancha(1,i)=(canchan(1,i)-0.01)+0.02*rand;
end
if i>50&i<=100
cancha(1,i)=(canchan(1,i)-0.005)+0.03*rand;
end
if i>100&i<=200
cancha(1,i)=(canchan(1,i)-0.008)+0.02rand;
end
if i>200&i<=256
cancha(1,i)=(canchan(1,i)-0.015)+0.02*rand;
end
end
hold on
plot(x,cancha,'k.');%画出残差值点
%图像重置MATLAB代码：
P=xlsread('C:\Users\hp\Desktop\A.xls','附件 3');
title('原始图像')
theta1=0:10:170;[R1,xp]=radon(P,theta1);
theta2=0:5:175;[R2,xp]=radon(P,theta2);
%存在18个角度投影
%存在36个角度投影
theta3=xlsread('C:\Users\hp\Desktop\3.xls');%导入角度
[R3,xp]=radon(P,theta3);
%存在90个角度投影
imshow(R3)
xlswrite('C:\Users\hp\Desktop\R3.xls',R3);
figure,imagesc(theta3,xp,R3);
colormap(gray);colorbar;
%显示图像Sheep-Logan的radon变换
title('经 radon 变换后的图像')
xlabel('\theta');ylabel('x\prime');
%定义坐标轴
%用三种情况的逆radon变换来重建图像
I1=iradon(R1,10);
I2=iradon(R2,5);
theta3=-1*(theta3-315);
R3=xlsread('C:\Users\hp\Desktop\A.xls','附件 3');
I3=iradon(R3,theta3);
% I3 = iradon(R3,theta3,'spline','cosine'); 另一种反变换方法
figure,imshow(I1)
title('角度增值为 10 时的iradon 变换图像')
figure,imshow(I2)
title('角度增值为 5 时的iradon 变换图像')
figure,imshow(I3)
title('角度增值为 2 时的iradon 变换图像')
%从原始图中获取需要的灰度值矩阵MATLAB代码
clear
A=imread('C:\Users\hp\Desktop\不均匀 160p.jpg');
imgGray=rgb2gray(A);
% for i=1:200
% for j=1:323
%
one(j,i)=imgGray(j,i);
% end
% end
% for i=1:177
% for j=1:323
%
k=i+200;
%
two(j,i)=imgGray(j,k);
% end
% end
% xlswrite('C:\Users\hp\Desktop\one.xls',one);
% xlswrite('C:\Users\hp\Desktop\two.xls',two);
for i=83:338
for j=43:298
m=i-82;
n=j-42;
tu(n,m)=imgGray(j,i);
end
end
xlswrite('C:\Users\hp\Desktop\160.xls',tu);
%结果图像进行滤波处理MATLAB代码
clear
I=imread('C:\Users\hp\Desktop\第三问答案.jpg');
imshow(I)
title('原始图像')
J=imnoise(I,'salt & pepper',0.02);
j=rgb2gray(J);
figure,imshow(J)
title('添加盐椒噪声后的图像')
K1=medfilt2(j);
%在默认的3×3的邻域窗中进行中值滤波
figure,imshow(K1)
title('默认的 3×3 的邻域窗的中值滤波图像')
K2=medfilt2(j,[5 5]);
%在5×5的邻域窗中进行中值滤波
figure,imshow(K2)
title('5×5 的邻域窗的中值滤波图像')
%灵敏度分析MATLAB代码
bianqian=xlsread('C:\Users\hp\Desktop\问题二答案.xls');
bianhou=xlsread('C:\Users\hp\Desktop\灵敏度.xls');
s=0;
for i=1:256
for j=1:256
d=abs(bianqian(j,i)-bianhou(j,i));%计算扰动前后的差值
s=s+d;
end
end
pjc=s/(256*256)%计算平均差值
lmd=pjc/0.05 %灵敏度值
