实验二CT重建实验报告.docx
《实验二CT重建实验报告.docx》由会员分享,可在线阅读,更多相关《实验二CT重建实验报告.docx(17页珍藏版)》请在冰豆网上搜索。
实验二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=