实验二CT重建实验报告.docx
《实验二CT重建实验报告.docx》由会员分享，可在线阅读，更多相关《实验二CT重建实验报告.docx（22页珍藏版）》请在冰点文库上搜索。
实验二CT重建实验报告
实验二CT图像重建实验
一、实验目的：
通过编写CT图像重建程序，进一步熟悉CT重建过程，同时加强图像处理程序的编程训练。
二、实验软件：
VC++
三、实验要求
1．递交整个程序的执行程序和源程序。
2．要求本小组承担部分的算法原理，列出主要程序段并给出相应说明，对实验结果进行分析。
四、算法原理及结果分析：
1.CT重建原理为：
在CT成像中，物体对X线的吸收起主要作用，在一均匀物体中，X线的衰减服从指数规律。
在X线穿透人体器官或组织时，由于人体器官或组织是由多种物质成分和不同的密度构成的，所以各点对X线的吸收系数是不同的。
将沿着X线束通过的物体分割成许多小单元体（体素），令每个体素的厚度相等（l）。
设l足够小，使得每个体素均匀，每个体素的吸收系数为常值，如果X线的入射强度I0、透射强度I和物体体素的厚度l均为已知，沿着X线通过路径上的吸收系数之和μ1＋μ2＋……＋μn就可计算出来。
为了建立CT图像，必须先求出每个体素的吸收系数μ1、μ2、μ3……μn。
为求出n个吸收系数，需要建立如上式那样n个或n个以上的独立方程。
CT成像装置从不同方向上进行多次扫描，来获取足够的数据建立求解吸收系数的方程。
吸收系数是一个物理量，它是CT影像中每个像素所对应的物质对X线线性平均衰减量大小的表示。
再将图像面上各像素的CT值转换为灰度，就得到图像面上的灰度分布，就是CT影像。
CT重建过程可以采用直接反投影和卷积反投影来实现。
卷积反投影重建图像时，先把由检测器上获得的原始数据与一个滤波函数进行了卷积运算，得到各方向卷积的投影函数；然后再把它们从各方向进行反投影，即按其原路径平均分配到每一矩阵元上，进行叠加后得到每一矩阵元的CT值；再经过适当处理后就可以得到被扫描物体的断层图像，卷积反投影可消除单纯的反投影产生的边缘失锐效应，补偿投影中的高频成分和降低投影中心密度，并保证重建图像边缘清晰和内部分布均匀。
2.整个程序的结构为：
◆产生sleeplogan模型
◆产生反投影数据
◆卷积反投影
编程在photostar平台下进行。
编程过程中有三个坐标系：
◆坐标系1：
图像坐标系，坐标原点在图像的左上方，x轴水平向右，y轴水平向下，单位为一个“像素”。
◆坐标系2：
空间坐标系，坐标原点在图像中心，x轴水平向右，y轴水平向上，单位为一个“像素”。
◆坐标系3：
归一化空间坐标系，坐标原点在图像中心，x轴水平向右，y轴水平向上，归一化因子为x轴的一半，所以对于m*m大小的图像，xy轴的坐标范围为[-1,1]。
根据坐标系1和坐标系3，可以确定图像的采样率。
根据面向对象的思想，在photostar平台中创建CCTEmulate类来管理实现CT仿真操作。
3.产生sleeplogan模型
要产生一幅图像，对于DIB来说，需要有信息头、文件头、调色板、像素区等部分，为操作简洁，我们读入一幅空白图像（text.bmp），指针m_pImageObject指向这幅空白图像的像素区，直接对像素区进行操作。
值得注意的是，根据实验需求，读入图像为灰度图像即可。
用CCTEmulate类中的SheepLogan函数相应菜单消息，进行仿真头模型的产生。
需要解决一个问题：
如何确定一个像素在哪个椭圆中？
根据椭圆参数方程：
其中AB分别是椭圆的长轴和短轴，x0，y0分别是椭圆的中心坐标，a是椭圆旋转角度。
则根据上式便可以确定在坐标中的一点是否在某一椭圆中。
SheepLogan函数的思路为：
遍历图像，检测图像的每一个像素点是否依次在某一椭圆中，若在a外，则置为背景色，在其他任何椭圆中，则置为该椭圆对应的灰度值。
所有的操作均在归一化坐标系中进行。
代码如下：
voidCPhotoStarView:
:
OnCteCreatsheeplogan（）
{
//TODO:
Addyourcommandhandlercodehere
CCTEmulateCCTE;
CCTE.m_pimageproject=m_pImageObject;
CCTE.SheepLogan（）;
CClientDCdc（this）;
Invalidate（）;
UpdateWindow（）;
}
BOOLCCTEmulate:
:
SheepLogan（）
{
//绘制Sheeplogan模型，坐标系参见高上凯《医学成像系统》附录II
//在此坐标系中，第（ij）个点的坐标（xy）为：
//x=（j-X）/Xy=（Y-i）/Y
//其中X=m_nWidth/2Y=m_nHeight/2
ASSERT（m_pimageproject!
=NULL）;//声明m_pImageObject不得为空
CImageObject*pImageObject=m_pimageproject;//重新用一个变量来进行以下操作
//**********获取图像的像素区指针和各种参数****************
intnWidth=pImageObject->GetWidth（）;//取得图像的宽度
intnHeight=pImageObject->GetHeight（）;//取得图像的高度
intnNumBits=pImageObject->GetNumBits（）;//取得图像的位数
intnWidthBytes;//图像每行所占的字节数
char*pBuffer=（char*）pImageObject->GetDIBPointer（&nWidthBytes）;
if（pBuffer==NULL）return（FALSE）;
BITMAPFILEHEADER*pBFH;//图像文件头指针
BITMAPINFOHEADER*pBIH;//图像信息头指针
RGBQUAD*pRGBPalette;//图像调色板指针
unsignedchar*pBits;//像素区的首指针
intnNumColors=pImageObject->GetNumColors（）;
pBFH=（BITMAPFILEHEADER*）pBuffer;
pBIH=（BITMAPINFOHEADER*）&pBuffer[sizeof（BITMAPFILEHEADER）];
pRGBPalette=（RGBQUAD*）&pBuffer[sizeof（BITMAPFILEHEADER）
+sizeof（BITMAPINFOHEADER）];
pBits=（unsignedchar*）&pBuffer[sizeof（BITMAPFILEHEADER）
+sizeof（BITMAPINFOHEADER）+nNumColors*sizeof（RGBQUAD）];
for（inti=0;i{
for（intj=0;j{
intk=0;//对十个椭圆进行判断，在椭圆内，则置该像素为预定灰度,每次都从0开始
switch（k）
{
case0:
//椭圆a
if（EllipseCaculate（0,0,0.92,0.69,0,1,i,j）>1）//不在椭圆内
{
pBits[（nWidth-1-i）*nWidthBytes+j]=0*255;
break;
}
elsepBits[（nWidth-1-i）*nWidthBytes+j]=1*255;//在椭圆内
case1:
//椭圆b
if（EllipseCaculate（0,-0.0184,0.874,0.6624,0,1,i,j）>1）break;
elsepBits[（nWidth-1-i）*nWidthBytes+j]=0.5*255;
case2:
//椭圆c
if（EllipseCaculate（0.22,0,0.31,0.11,0.3090,0.9511,i,j）<=1）
{
pBits[（nWidth-1-i）*nWidthBytes+j]=int（0.3*255）;
break;
……以下椭圆省略…..
inlinefloatCCTEmulate:
:
EllipseCaculate（floatcentX,floatcentY,floatlongaxes,floatshortaxes,floatcosangle,floatsinangle,inti,intj）
{
//返回（ij）在第n个椭圆公式中的计算值
intn=m_pimageproject->GetWidth（）;
floatX=n/2.0;
intm=m_pimageproject->GetHeight（）;
floatY=m/2.0;
floatx=（j-X）/X;
floaty=（Y-i）/Y;
return（pow（（x-centX）*cosangle+（y-centY）*sinangle,2）/pow（longaxes,2）+pow（（centX-x）*sinangle+（y-centY）*cosangle,2）/pow（shortaxes,2））;
}
4.产生投影数据
在产生投影数据之前，要模拟投影数据的产生过程。
根据X射线对物体投影的过程，设定每10度投影一次，每次模拟产生过图像中心的13条直线。
要在图像上画线，本质上还是要在图像的像素区中进行操作，将图像中在指定直线上的像素置为255，则可以在图像上画线。
在每个角度产生了画了直线的图像之后，在画下一个角度的投影直线时，需要在原图像上进行操作，所以，每次进行操作之前，应该先保留原图像的像素区，每次画线操作完毕进行显示了之后，将原图像像素区放回到原图像中去。
这样可以进行下一次的有效操作。
对于每次的画线操作，在归一化空间坐标系中，根据直线方程
，角度r从0~180，每10度画线一次，
和
分别代表过原点的左右两边的直线。
每次画线
从0~0.6，每0.1画线一次。
对于图像坐标系中的点（m，n）经过坐标转换到归一化空间坐标系中为（x，y），当xy满足fabs（x*cos
+y*sin
-
）<0.005则该点（m，n）该置为255.
代码如下：
BOOLCCTEmulate:
:
ShowEmulate（CDC*pDC,intm）
{
//本函数模拟CT采样模型，用一块新的空间来保存原图像中的像素区，
//然后用新的像素区刷新显示，最后恢复原图像的像素区
//math文件中三角函数中角度的单位是弧度
//直线的法线方程为xcosr+ysinr=R
//旋转角度m
…….预备工作的代码略
//坐标转换工作
//在Sheeplogan模型中，坐标系参见高上凯《医学成像系统》附录II
//在此坐标系中，第（ij）个点的坐标（xy）为：
//x=（j-X）/Xy=（Y-i）/Y
//其中X=m_nWidth/2Y=m_nHeight/2
floatX=nWidth/2;
floatY=nHeight/2;
//检查xy是否在模拟的直线上
//下面是一次模拟，完成之后需要显示一次，然后进行下一次
for（floatk=0;k<0.7;k=k+0.1）//法线方程中的长度，
{
floatr1=m/180.0*3.1415;//将角度转换为弧度
floatr2=（m+180）/180.0*3.1415;//将角度转换为弧度
for（inti=0;i{
for（intj=0;j{
floatx=（j-X）/X;
floaty=（Y-i）/Y;
if（fabs（x*cos（r1）+y*sin（r1）-k）<0.005||fabs（x*cos（r2）+y*sin（r2）-k）<0.005）
{
//xy在直线上
pBits[（nHeight-1-i）*nWidthBytes+j]=255;
}
}//for（intj
}//for（inti=0
}//for（intk=0
}
产生投影数据主要解决两个问题：
每一条扫描线的投影值如何存放？
每一条扫描线上的投影值是多少？
安排数据结构如下：
扫描共计180次，角度为0~179，图像长或者宽的最大值的一半设为projectionhalflength，则每次扫描线数为projectionlength=2*projectionhalflength-1，安排数组为projectiondata[180*projectionlength].根据直线的法线方程
，则在第r个角度的、法线为
的数据存放的位置为projectiondata[r*projectionlength+k],其中
或者
，
。
其中当
时，
，当
是，
，他们对应不同的直线。
中间的一条直线标号为[projectionhalflength-1],对应过原点的直线。
对于特定的角度r，则一排扫描线从0~projectionlength-1，存放了从
到
的投影值。
根据r和
与数组projectiondata之间的关系，可以确定在角度r和法线
确定的直线的投影值存放的位置。
根据角度
和法线
，对所有存在的直线进行扫描。
思路是：
双重循环r=0~179，k=0~projectionlength-1，这样就确定了一条直线的参数
和
。
在双重循环中扫描整幅图像，搜寻这条直线上的像素点（坐标满足直线方程的点就认为在该条直线上），然后将该像素点的灰度值加到对应的投影值中去。
代码如下：
BOOLCCTEmulate:
:
GetProjectionData（double（*pProjDa）[180]）
{
//获得图像的投影数据
//以下的坐标系仍采用高上凯《医学成像系统》附录II中坐标系
//ProjDa数组用于存放采样的数据，大小应为m*nm为360（角度），n为图像长宽最大值的一半
ASSERT（m_pimageproject!
=NULL）;
//获得图像相应结构数据
intnWidthBytes;
DWORDnColorNum;
char*pBuffer=（char*）m_pimageproject->GetDIBPointer（&nWidthBytes）;
if（!
pBuffer）returnFALSE;
nColorNum=m_pimageproject->GetNumColors（）;
unsignedchar*pBits;
pBits=（unsignedchar*）&pBuffer[sizeof（BITMAPFILEHEADER）+sizeof（BITMAPINFOHEADER）
+nColorNum*sizeof（RGBQUAD）];
intnWidth=m_pimageproject->GetWidth（）;
intnHeight=m_pimageproject->GetHeight（）;
//进行坐标转换
floatX=nWidth/2;
floatY=nHeight/2;
//为了在数组中将数据与角度和法线长度（根据直线的法线方程）联系起来
//角度从0旋转到179（整数），法线长度从0~图像长和宽的最大值的一半（整数）
//求每次法线直线方程中法线长度,决定了数组的大小，非常重要！
//数组的第二维大小为2*nProjectionHalfLength-1,从0开始存放坐标系上r+180和r的直线投影的数据
//数组第二维的中心是标号为[nProjectionHalfLength]的点，对应过原点的直线
intnProjectionHalfLength=（nWidth>nHeight）?
nWidth:
nHeight;
nProjectionHalfLength=（int）（（nProjectionHalfLength+1）/2.0）;
for（intk=0;k{
for（inti=0;i<180;i++）//进行投影，每一度一次
{
//法线长度转换到对应坐标系中的数值
floatR=（float）k;
R=R/X;//法线长度,坐标转换的尺度应该相同，这里选择X进行转换，Y也一样
floatr1=i/180.0*3.1415;//把角度转换为弧度
floatr2=（i+180）/180.0*3.1415;
for（intm=0;m{
for（intn=0;n{
//进行坐标转换
floatx=（n-X）/X;
floaty=（Y-m）/Y;
if（fabs（x*cos（r1）+y*sin（r1）-R）<0.005）
{
//根据i与r，k与R之间的关系，可明确每一个数据对应的像素直线
pProjDa[nProjectionHalfLength-1+k][i]=pProjDa[nProjectionHalfLength+k][i]+pBits[（nHeight-1-m）*nWidthBytes+n];
}
if（R==0）continue;//过原点的直线在最中间，并且只算一次
if（fabs（x*cos（r2）+y*sin（r2）-R）<0.005）
{
//根据i与r，k与R之间的关系，可明确每一个数据对应的像素直线
pProjDa[nProjectionHalfLength-1-k][i]=pProjDa[nProjectionHalfLength-k][i]+pBits[（nHeight-1-m）*nWidthBytes+n];
}
}//for（intn=
}//for（intm=0
}//for（floati=0;
}//for（intk
}
这种思路的优点在于，思路简单，比较精确，用阈值0.005可以让计算机相对准确找到在对应直线上的像素。
缺点在于，四重循环，计算量太大，测试发现，该函数在调试状态下的运行时间大约是10分钟。
换一种思路，不必让计算机通过遍历整个图像来计算哪些点在该条直线上，根据直线方程，我们能够通过遍历x轴（相当于得到列的信息）来计算y轴（对应的行的信息），则可以得到在原图像坐标系下对应像素的坐标，这样计算一条直线就可以只需要三重循环，速度大大加快。
值得注意的是，如果我们要在像素区中操作某一个像素的灰度值，在需要得到在图像坐标系下的该像素的坐标，如果我们要计算判断直线，则应该在空间坐标系下进行。
代码如下：
BOOLCCTEmulate:
:
GetProjectionData（double*pProjDa）
{
//获得图像的投影数据
//以下的坐标系仍采用高上凯《医学成像系统》附录II中坐标系
//ProjDa数组用于存放采样的数据，大小应为m*nm为360（角度），n为图像长宽最大值的一半
ASSERT（m_pimageproject!
=NULL）;
//获得图像相应结构数据
intnWidthBytes;
DWORDnColorNum;
char*pBuffer=（char*）m_pimageproject->GetDIBPointer（&nWidthBytes）;
if（!
pBuffer）returnFALSE;
nColorNum=m_pimageproject->GetNumColors（）;
unsignedchar*pBits;
pBits=（unsignedchar*）&pBuffer[sizeof（BITMAPFILEHEADER）+sizeof（BITMAPINFOHEADER）
+nColorNum*sizeof（RGBQUAD）];
intnWidth=m_pimageproject->GetWidth（）;
intnHeight=m_pimageproject->GetHeight（）;
//进行坐标转换
floatX=nWidth/2.0;
floatY=nHeight/2.0;
//为了在数组中将数据与角度和法线长度（根据直线的法线方程）联系起来
//角度从0旋转到179（整数），法线长度从0~图像长和宽的最大值的一半（整数）
//求每次法线直线方程中法线长度,决定了数组的大小，非常重要！
//数组的第二维大小为2*nProjectionHalfLength-1,从0开始存放坐标系上r+180和r的直线投影的数据
//数组第二维的中心是标号为[nProjectionHalfLength-1]的点，对应过原点的直线
intnProjectionHalfLength=（nWidth>nHeight）?
nWidth:
nHeight;
nProjectionHalfLength=（int）（（nProjectionHalfLength+1）/2.0）;
intnProjectionLength=2*nProjectionHalfLength-1;
for（inti=0;i<180;i++）//进行投影，每一度一次
{
floatr1=i/180.0*3.1415;//把角度转换为弧度
floatr2=（i+180）/180.0*3.1415;
for（intk=0;k{
floatL=float（k-nProjectionHalfLength+1）;//法线长度
intn=0;//原图像的列,从0开始计算
while（k{
floatx=n-X;
floaty=（L-x*cos（r2））/sin（r2）;/由方程和图像列的信息推算行的信息
intm=int（Y-y）;//取整有一定误差，原图像的行，也从0开始计算
n++;
if（m<0||m>=nHeight）continue;//如果m超出图像范围，则跳过
pProjDa[i*nProjectionLength+k]=pProjDa[i*nProjectionLength+k]+pBits[（nHeight-1-m）*nWidthBytes+n];
}
n=0;//对每一条直线x轴都要扫描一次
while（k>=nProjectionHalfLength-1&&n{
floatx=n-X;
flo
升级会员 