2017数学建模问题二问题三
需要添加符号系统:
m ,从模板接受信息经过逆变换到图形矩阵的转换系数。
Xi yi 还有Pij1 Pij2
问题二
该问题是问题一的逆向过程,根据问题一得到的标定参数,利用上述CT 系统得到的某未知介质的接收信息,确定该未知介质在正方形托盘中的位置、几何形状和吸收率等信息。
由Radon 原理:
两维情况下Radon 变换大致可以这样理解:一个平面内沿不同的直线(直线与原点的距离为r,方向角为θ)对f(x, y) 做线积分,得到的像R(r, θ) 就是函数f(x, y) 的Radon 变换。也就是说,平面(r, θ) 的每个点的像函数值对应了原始函数的某个线积分。
在本题中,对于刚好照射在每一探测器上的X 射线刚好对应于直线 l ,用X 射线去照射模板或待检测物质,相当于在X 射线对应的直线 l上做线积分,由于r 及θ的不同,得到不同的R(r, θ) ,对应于所得到的接收信息。
对Radon 变换简单理解如图1,
上图中,在直角坐标系xOy中,f x, y 为直线l 上的点,由于x=rcosθ,y=rsinθ,从而对每一 r, θ ,图中直线l 可表示为:
xcosθ+ysinθ=r
沿直线 l 对f(x,y) 做线积分,得Radon 变换公式为:
R r, θ = −∞f x, y ds ,其中s = x
在本题中,Radon 变换可直观理解为:根据不同射线的入射方向,对位于待检测介质或模板内部的射线的衰减系数进行积分,可得待检测物质或模板的衰减函数R(r,θ) ,如图2所示。
在二维平行束投影示意图中,射线目标信息为函数 f(x,y) , θ为X
射线的法线与x 轴夹角,r 为射线AB 与原点的距离。[1] +∞
另,狄拉克 δ 函数是一个广义函数,在物理学中常用其表示质点、点电荷等理想模型的密度分布,该函数在除了零以外的点取值都等于零,而其在整个定义域上的积分等于1
为了便于理解及应用,我们写出一种最简单的狄拉克 δ 函数形式:
δ x = 0 ,x=0 1 ,x≠0
将此狄拉克 δ 函数应用到X 射线所对应的直线l上,有
0 ,xcosθ+ysinθ−r=0δ xcosθ+ysinθ−r = 1 ,xcosθ+ysinθ−r≠0
即在直线 l 上的点f(x,y) 满足δ xcos θ+ysin θ−r =1,其他不在直线 l 上的点满足δ xcos θ+ysin θ−r =0。
此时,R r, θ 为探测器所接收到的信息,由上述论证,此时Radon 变换为:
R r, θ = +∞
f x, y δ xcosθ+ysinθ−r ds
−∞
Radon 变换可理解为某一二维图像在 r, θ 空间的投影, r, θ 空间的每一点对应(x,y) 空间的一条直线方程。
补充:实际上,R r, θ 定义在空间中一个半圆柱体的侧面上,如图3。
Radon 逆变换(Iradon )由投影后的R r, θ ,Radon 逆变换可以求出原二维图像的几何形状及对X 射线的吸收率等信息,即给出从投影到重建的解。具体表达式为:
f x,y =1ðθ dθ dr 0−∞π∞ð(r, θ) 根据上述数学原理,结合Matlab 中Radon 及Iradon 函数进行求解。
(1)利用Matlab 的Iradon 工具箱。我们将512*180的接收信息矩阵,通过radon 逆变换,以及平移图像中心,得到由五个椭圆组成的“外星人”状图像。颜色越浅“吸收率”越大。如图所示:
图x 附件三接收图
图x 附件三介质原始图像的几何形状
(2)图像在正方形托盘的位置确定。
由256*256个像素点的图像矩阵对应到100mm*100mm的正方形托盘,每一个像素点的长为0.3906mm ,得到该未知介质在正方形托盘的位置信息,以正方形托盘中心为坐标原点(0,0)作图:
图x 附件三未知介质在正方形托盘的位置
(3)对应吸收率的确定。 首先根据附件二的模板接收信息,进行radon 逆变换得到模板的图像信息,寻找逆变换得到的图像信息与原始图像之间的平均比例关系m 。
ij1pij1m = 1≤i ≤256,1≤j ≤256, ∈R ij2ij2
用matlab 求解m 程序见附录1xxx 。解得m = 2.0751,即从接收信息进行radon 逆变换得到原始图像的真实“吸收率”为函数返回值的2倍。此时再读入附件三即未知介质的接收矩阵,通过乘以转换系数m 再进行radon 逆变换得到10个点的吸收率如下表(代码见附录2xxx ):
问题三:
该问题同样是radon 的逆变换,利用问题一中得到的标定参数,根据附件五的接收信息,给出该未知介质的相关信息。
我们首先根据之前的方法利用radon 反变换得到图像:
图xxx 问题三未去噪得到的图形
此时发现图像有太多的毛刺,通过查阅资料得知图像受到噪声干扰影响,必须进行去噪操作。
由滤波反投影法(FBP ):
滤波反投影法根据附件三所给接收信息,采用先修正、后投影重建图像的做法,可得到原始图像的吸收率信息。其原理为:在得到某一角度下的投影函数
(一
维函数)后,对此函数做滤波处理,得一修正后的滤波函数,再对修正后的滤波函数做反投影运算,得待检测介质吸收率在正方形托盘中的每一点的分布密度函数f(x, y) 。图1给出了滤波反投影法重建原始图像的流程图。
图1 滤波反投影法流程图
反投影法重建原始图像的步骤:
(1) 在对应于投影函数的角度下对投影函数做一维Fourier 变换;
(2) 对(1)得到的变换结果乘以权重因子|ρ|;
(3) 对(2)加权后得到的结果做一维傅立叶;
(4) 对(3)所得函数做直接反投影;
(5) 改变投影角度,得到180个不同的投影角度,对每一角度,重复上述
步骤(1)~(4)。
R-L (Ram-Lak )滤波函数:
此函数的基本条件是二维图像函数的频率是有界的,显然,此题所得附件五的所有数据满足此条件。故频域中的滤波函数可表示为:
G ρ =
其函数图像如图1.
图1 R-L滤波函数图像
连续的R-L 卷积函数所得结果为:
g R =ρ02[2sin c 2ρ0R −sin c2 ρ0R ]
离散的R-L 卷积函数所得结果为:
1,n=0 g nT = 0 ,n为偶数 1 − ,n为奇数
(1)根据上述滤波原理,在本题中,通过matlab 的iradon 函数的’Ram-Lak ’可以实现滤波,通过观察得到的图像矩阵,发现“吸收率”过小的像素点会造成图像毛刺,去除图像中“吸收率”过小的部分进一步降低噪声,(代码见附录3***)。 ρ , ρ ≤ρ0 0,其它
得到图像由256*256像素组成的几何形状:
图xxx 问题三降噪后的图形
(2)求介质在托盘中的位置。
由256*256个像素点的图像矩阵对应到100mm*100mm的正方形托盘,每一个像素点的长为0.3906mm ,得到该未知介质在正方形托盘的位置信息,以正方形托盘中心为坐标原点(0,0)作图:
图xxx 问题三介质在托盘中的位置
(3)求该介质所对应的“吸收率”。 读入附件五即未知介质的接收矩阵,通过乘以转换系数m 再进行radon
逆变
换得到10个点的吸收率如下表:
表xxx 问题三介质的10个位置“吸收率”
附录1xxx
R = load('data2.txt' );
theate = [119.7146 121.0580 121.6095 122.6918 123.7178 124.6815 125.6760 126.6708 127.6657 128.6609 129.6562 130.6517 131.6475 132.6431 133.6388 134.7845 135.6306 136.6265 137.6224 138.6183 139.6142 140.6101 141.6059 142.6017 143.5973 144.5929 145.5883 146.5837 147.5788 148.5738 149.5687 150.4628 151.5576 152.5516 153.5454 154.5389 155.5319 156.5245 157.5166 158.5083 159.4991 160.4892 161.4787 162.4670 163.4542 164.4401 165.4244 166.4071 167.3873 168.3652 169.1996 170.0354 170.9975 172.1042 173.1954 174.7117 175.6289 176.2412 177.2414 178.2953 179.1542 179.9993 180.7137 182.9672 185.7122 185.8886 186.3366 186.9725 187.7274 188.5559 189.4302 190.3357 191.2622 192.2032 193.1552 194.1154 195.0821 196.0536 197.0290 198.0080 198.9893 199.9728 200.9584 201.9452 202.9336 203.9231 204.9136 205.9051 207.6891 208.8839 209.8779 210.8726 211.8677 212.8634 213.8594 214.8557 215.8524 216.8495 217.8469 218.8446 219.8425 220.8409 221.9392 222.8382 223.8373 224.8366 225.8364 226.8362 227.8363 228.8370 229.8377 230.8389 231.8402
232.8420 233.8442 234.8468 235.8498 236.8533 237.8573 238.8618 239.8670 240.8728 241.8793 242.8867 243.8950 244.9042 245.9146 246.9263 247.9392 248.9538 249.9699 250.9879 252.0077 253.0296 254.0536 255.0795 256.1074 257.1369 258.1675 259.1987 260.2298 261.2601 262.2889 263.3160 264.3415 265.3661 266.3915 267.4198 268.4546 269.5013 270.5680 271.6704 272.8431 274.1979 274.8686 275.7084 276.2692 276.9481 277.4523 278.0057 278.5599 279.9215 281.1024 282.2136 283.2891 284.3434 285.3842 286.4159 287.4410 288.4613 289.4780 290.4919 291.5035 292.5135 293.5220 294.5294 295.5358 296.5415 297.5464 298.5401
] - 90;
[I1,H] = iradon(R,theate,'spline' ,512);
I1(:,(468:512)) = [];
I1(:,(1:110)) = [];
I1((457:512),:) = [];
I1((1:99),:) = [];
I1 = imresize(I1,256/357);
data1 = load('data1.txt' );
kk=0;
Tmp = data1 ./ I1;
for i = 1:1:256
for j = 1:1:256
if (Tmp(i,j) ~= 0)
kk = kk + 1;
end
end
end
m = sum(Tmp(:)) / (kk);
附录2xxx
R = load('data3.txt' );
R = R .* 2.0751;
theate = [119.7146 121.0580 121.6095 122.6918 123.7178 124.6815 125.6760 126.6708 127.6657 128.6609 129.6562 130.6517 131.6475 132.6431 133.6388 134.7845 135.6306 136.6265 137.6224 138.6183 139.6142 140.6101 141.6059 142.6017 143.5973 144.5929 145.5883 146.5837 147.5788 148.5738 149.5687 150.4628 151.5576 152.5516 153.5454 154.5389 155.5319 156.5245 157.5166 158.5083 159.4991 160.4892 161.4787 162.4670 163.4542 164.4401 165.4244 166.4071 167.3873 168.3652 169.1996 170.0354 170.9975 172.1042 173.1954 174.7117 175.6289 176.2412 177.2414 178.2953 179.1542 179.9993 180.7137 182.9672 185.7122 185.8886 186.3366 186.9725 187.7274 188.5559 189.4302 190.3357 191.2622 192.2032 193.1552 194.1154 195.0821 196.0536 197.0290 198.0080 198.9893 199.9728 200.9584 201.9452 202.9336 203.9231 204.9136 205.9051 207.6891 208.8839 209.8779 210.8726 211.8677 212.8634 213.8594 214.8557 215.8524 216.8495 217.8469 218.8446 219.8425
220.8409 221.9392 222.8382 223.8373 224.8366 225.8364 226.8362 227.8363 228.8370 229.8377 230.8389 231.8402 232.8420 233.8442 234.8468 235.8498 236.8533 237.8573 238.8618 239.8670 240.8728 241.8793 242.8867 243.8950 244.9042 245.9146 246.9263 247.9392 248.9538 249.9699 250.9879 252.0077 253.0296 254.0536 255.0795 256.1074 257.1369 258.1675 259.1987 260.2298 261.2601 262.2889 263.3160 264.3415 265.3661 266.3915 267.4198 268.4546 269.5013 270.5680 271.6704 272.8431 274.1979 274.8686 275.7084 276.2692 276.9481 277.4523 278.0057 278.5599 279.9215 281.1024 282.2136 283.2891 284.3434 285.3842 286.4159 287.4410 288.4613 289.4780 290.4919 291.5035 292.5135 293.5220 294.5294 295.5358 296.5415 297.5464 298.5401
] - 90;
[I1,H] = iradon(R,theate,'spline' ,512);
I1(:,(468:512)) = [];
I1(:,(1:110)) = [];
I1((457:512),:) = [];
I1((1:99),:) = [];
I1 = imresize(I1,256/357);
figure(1)
imshow(I1,[]);
附录三:
R = load('data5.txt' );
R = R .* 2.0751;
theate = [119.7146 121.0580 121.6095 122.6918 123.7178 124.6815 125.6760 126.6708 127.6657 128.6609 129.6562 130.6517 131.6475 132.6431 133.6388 134.7845 135.6306 136.6265 137.6224 138.6183 139.6142 140.6101 141.6059 142.6017 143.5973 144.5929 145.5883 146.5837 147.5788 148.5738 149.5687 150.4628 151.5576 152.5516 153.5454 154.5389 155.5319 156.5245 157.5166 158.5083 159.4991 160.4892 161.4787 162.4670 163.4542 164.4401 165.4244 166.4071 167.3873 168.3652 169.1996 170.0354 170.9975 172.1042 173.1954 174.7117 175.6289 176.2412 177.2414 178.2953 179.1542 179.9993 180.7137 182.9672 185.7122 185.8886 186.3366 186.9725 187.7274 188.5559 189.4302 190.3357 191.2622 192.2032 193.1552 194.1154 195.0821 196.0536 197.0290 198.0080 198.9893 199.9728 200.9584 201.9452 202.9336 203.9231 204.9136 205.9051 207.6891 208.8839 209.8779 210.8726 211.8677 212.8634 213.8594 214.8557 215.8524 216.8495 217.8469 218.8446 219.8425 220.8409 221.9392 222.8382 223.8373 224.8366 225.8364 226.8362 227.8363 228.8370 229.8377 230.8389 231.8402 232.8420 233.8442 234.8468 235.8498 236.8533 237.8573 238.8618 239.8670 240.8728 241.8793 242.8867 243.8950 244.9042 245.9146 246.9263 247.9392 248.9538 249.9699 250.9879 252.0077 253.0296 254.0536 255.0795 256.1074 257.1369 258.1675 259.1987 260.2298 261.2601 262.2889 263.3160 264.3415 265.3661 266.3915 267.4198 268.4546 269.5013 270.5680 271.6704 272.8431 274.1979 274.8686 275.7084 276.2692 276.9481 277.4523 278.0057 278.5599 279.9215 281.1024 282.2136 283.2891 284.3434 285.3842 286.4159 287.4410 288.4613 289.4780 290.4919 291.5035 292.5135 293.5220 294.5294 295.5358 296.5415 297.5464 298.5401
] - 90;
[I1,H] = iradon(R,theate,'spline' ,512);
I1(:,(468:512)) = [];
I1(:,(1:110)) = [];
I1((457:512),:) = [];
I1((1:99),:) = [];
I1 = imresize(I1,256/357);
for i = 1:1:256 for j = 1:1:256
if (abs(I1(i,j))
end
end
figure(1)
imshow(I1,[]);