试验二CT重建试验报告

1、实验二CT图像重建实验一、实验目的：通过编写CT图像重建程序，进一步熟悉CT重建过程，同时加强图像处理程序的编程训练。二、实验软件：VC+三、实验要求1 .递交整个程序的执行程序和源程序。2 .要求本小组承担部分的算法原理，列出主要程序段并给出相应说明， 对实验 结果进行分析。四、算法原理及结果分析：1 . CT重建原理为：在CT成像中，物体对 X线的吸收起主要作用，在一均匀物体中，X线的衰减服从指数规律。在 X线穿透人体器官或组织时，由于人体器官或组织是由多种物质成分和不同的 密度构成的，所以各点对 X线的吸收系数是不同的。将沿着 X线束通过的物体分割成许多 小单元体（体素），令每个体素的厚

2、度相等（I）。设l足够小，使得每个体素均匀，每个体素 的吸收系数为常值，如果 X线的入射强度10、透射强度I和物体体素的厚度I均为已知，沿 着X线通过路径上的吸收系数之和W+戌+ P就可计算出来。为了建立CT图像，必须先求出每个体素的吸收系数口、出、国卬。为求出n个吸收系数，需要建立如上式那样 n个或n个以上的独立方程。CT成像装置从不同方向上进行 多次扫描，来获取足够的数据建立求解吸收系数的方程。吸收系数是一个物理量， 它是CT影像中每个像素所对应的物质对 X线线性平均衰减量 大小的表示。再将图像面上各像素的 CT值转换为灰度，就得到图像面上的灰度分布，就是 CT影像。CT重建过程可以采用直

3、接反投影和卷积反投影来实现。卷积反投影重建图像时，先把-可编辑修改-由检测器上获得的原始数据与一个滤波函数进行了卷积运算，得到各方向卷积的投影函数;然后再把它们从各方向进行反投影，即按其原路径平均分配到每一矩阵元上，进行叠加后得到每一矩阵元的 CT值；再经过适当处理后就可以得到被扫描物体的断层图像，卷积反投影可消除单纯的反投影产生的边缘失锐效应，补偿投影中的高频成分和降低投影中心密度，并保证重建图像边缘清晰和内部分布均匀。2 .整个程序的结构为：产生sleep logan 模型产生反投影数据卷积反投影编程在photostar 平台下进行。编程过程中有三个坐标系： 坐标系1 ：图像坐标系，坐标原

4、点在图像的左上方，x轴水平向右，y轴水平向下，单位为一个“像素”。 坐标系2 ：空间坐标系，坐标原点在图像中心，x轴水平向右，y轴水平向上，单位为一个“像素”。 坐标系3 :归一化空间坐标系，坐标原点在图像中心，x轴水平向右，y轴水平向上，归一化因子为x轴的一半，所以对于 m*m大小的图像，xy轴的坐标范围为 卜1,1。根据坐标系1和坐标系3,可以确定图像的采样率。根据面向对象的思想，在 photostar 平台中创建CCTEmulate类来管理实现CT仿真 操作。-可编辑修改-3.产生sleep logan 模型要产生一幅图像，对于 DIB来说，需要有信息头、文件头、调色板、像素区等部分，为

5、操作简洁，我们读入一幅空白图像(text.bmp ),指针m_pImageObject指向这幅空白图像的像素区，直接对像素区进行操作。值得注意的是，根据实验需求，读入图像为灰度图像即可。用CCTEmulate类中的SheepLogan 函数相应菜单消息，进行仿真头模型的产生。需 要解决一个问题：如何确定一个像素在哪个椭圆中？根据椭圆参数方程：22(x X0)cos(y y°)sin )( (x x°)sin(y y°)cos )1221A2B2其中A B分别是椭圆的长轴和短轴，x0, y0分别是椭圆的中心坐标，a是椭圆旋转角度。则根据上式便可以确定在坐标中的一点是

6、否在某一椭圆中。SheepLogan 函数的思路为：遍历图像，检测图像的每一个像素点是否依次在某一椭圆 中，若在a外，则置为背景色，在其他任何椭圆中，则置为该椭圆对应的灰度值。所有的 操作均在归一化坐标系中进行。代码如下：voidCPhotoStarView:OnCteCreatsheeplogan()/ TODO: Add your command handler code hereCCTEmulateCCTE;-可编辑修改-CCTE.m_pimageproject = m_pImageObject;CCTE.SheepLogan();CClientDCdc(this);Invalidate

7、();UpdateWindow();BOOLCCTEmulate:SheepLogan()绘制Sheeplogan模型，坐标系参见高上凯医学成像系统附录 II在此坐标系中，第(ij)个点的坐标(xy)为：/x =(j-X) /X y = (Y-i) /Y其中 X = m_nWidth/2 Y = m_nHeight/2ASSERT(m_pimageproject != NULL); / 声明 m_pImageObject 不得为空CImageObject* pImageObject = m_pimageproject;/重新用一个变量来进行以下操作/* 获取图像的像素区指针和各种参数 *int

8、nWidth = pImageObject->GetWidth();/ 取得图像的宽度intnHeight = pImageObject->GetHeight(); /取得图像的高度-可编辑修改-取得图像的位数intnNumBits = pImageObject->GetNumBits();/ intnWidthBytes; /图像每行所占的字节数char *pBuffer = (char *) pImageObject->GetDIBPointer( &nWidthBytes );if( pBuffer = NULL ) return(FALSE);BITMA

9、PFILEHEADER *pBFH;/图像文件头指针BITMAPINFOHEADER *pBIH;/图像信息头指针RGBQUAD *pRGBPalette;/图像调色板指针unsignedchar *pBits;/像素区的首指针intnNumColors = pImageObject->GetNumColors();pBFH = (BITMAPFILEHEADER *) pBuffer;pBIH = (BITMAPINFOHEADER *) &pBuffersizeof(BITMAPFILEHEADER);pRGBPalette = (RGBQUAD *) &pBuffe

10、rsizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER);pBits = (unsignedchar *) &pBuffersizeof(BITMAPFILEHEADER)+sizeof(BITMAPINFOHEADER)+nNumColors*sizeof(RGBQUAD);for (inti=0;i<nWidth;i+)对每个像素点进行灰度赋值i为行，j为列for (intj=0;j<nHeight;j+)-可编辑修改-intk=0; 对十个椭圆进行判断，在椭圆内，则置该像素为预定灰度，每次都从0开始switch(k)(cas

11、e 0: / 椭圆 aif (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;/ 在椭圆内case 1:/ 椭圆 bif (EllipseCaculate(0,-0.0184,0.874,0.6624,0,1,i,j)>1) break;elsepBits(nWidth-1-i)*nWidthBytes +j = 0.5*255;case 2: / 椭圆

12、 cif (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_pimageprojec

13、t->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射线对物体投影的过程，设定每1

14、0度投影一次，每次模拟产生过图像中心的13条直线。要在图像上画线，本质上还是要在图像的像素区中进行操作，将图像中在指定直线上的像素置为255 ,则可以在图像上画线。在每个角度产生了画了直线的图像之后，在画下一个角度的投影直线时，需要在原图像上进行操作，所以，每次进行操作之前，应该先保留原图像的像素区，每次画线操作完毕进-可编辑修改-行显示了之后，将原图像像素区放回到原图像中去。这样可以进行下一次的有效操作。对于每次的画线操作，在归一化空间坐标系中，根据直线方程xcos ysin ,角度r从0180 ,每10度画线一次，r/180* 和 (r 180) /180* 分别代表过原点的左右两边的直线

15、。每次画线从00.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模

16、型中，坐标系参见高上凯医学成像系统附录 II在此坐标系中，第(ij)个点的坐标(xy)为：/x =(j-X) /X y = (Y-i) /Y其中 X = m_nWidth/2 Y = 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;

17、/ 将角度转换为弧度for (inti = 0;i<nHeight;i+) /遍历整个图像,hang(for (intj = 0;j<nWidth;j+) /lie(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)(/x y 在直线上pBits(nHeight-1-i)*nWidthBytes + j = 255;)-可编辑修改-/for (int j/for (int i = 0/for (int k

18、=0 产生投影数据主要解决两个问题：每一条扫描线的投影值如何存放？每一条扫描线上的投影值是多少？安排数据结构如下：扫描共计180次，角度为0179 ,图像长或者宽的最大值的一半设为 projectionhalflength,贝U每次扫描线数为 projectionlength = 2*projectionhalflength-1 ,安排数组为projectiondata180*projectionlength.根据直线的法线方程xcos ysin ,则在第r个角度的、法线为的数据存放的位置为projectiondatar*projectionlength + k,其中 r /180* 或者(r

19、 180)/180*,| k (projectionhalflength 1)|。其中当k projectionhalflength 1 时， (r 180) /180* ，当 k projectionhalflength 1 是， r /180*,他们对应不同的直线。中间的一条直线标号为projectionhalflength-1,对应过原点的直线。对于特定的角度r,则一排扫描线从0projectionlength-1 ,存放了从 (r 180) /180*至Ur /180*的投影值。根据r和 与数组projectiondata之间的关系，可以确定在角度r和法线确定的直线的投影值存放的位置。

20、根据角度 和法线，对所有存在的直线进行扫描。思路是：双重循环r = 0179 , k =0projectionlength-1 ,这样就确定了一条直线的参数和 。在双重循环中扫描整幅图像，搜寻这条直线上的像素点(坐标满足直线方程的点就认为在该条直线上)，然后将该像-可编辑修改-素点的灰度值加到对应的投影值中去。代码如下：BOOLCCTEmulate:GetProjectionData(double (*pProjDa)180)(/获得图像的投影数据/以下的坐标系仍采用高上凯医学成像系统附录II中坐标系/ProjDa数组用于存放采样的数据，大小应为 m*n m为360 (角度)，n为图像长宽最大

21、值的一半ASSERT(m_pimageproject != NULL);/获得图像相应结构数据intnWidthBytes;DWORDnColorNum;char *pBuffer = (char *)m_pimageproject->GetDIBPointer(&nWidthBytes);if (!pBuffer) returnFALSE;nColorNum = m_pimageproject->GetNumC010rs();unsignedchar* pBits;pBits =(unsignedchar* )&pBuffersizeof(BITMAPFILEHE

22、ADER)+sizeof(BITMAPIN-可编辑修改-FOHEADER)+nColorNum*sizeof(RGBQUAD);intnWidth = m_pimageproject->GetWidth();intnHeight = m_pimageproject->GetHeight();/进行坐标转换floatX = nWidth/2;floatY = nHeight/2;/为了在数组中将数据与角度和法线长度(根据直线的法线方程)联系起来角度从0旋转到179 (整数)，法线长度从0图像长和宽的最大值的一半(整数)求每次法线直线方程中法线长度，决定了数组的大小，非常重要！/数组的

23、第二维大小为2*nProjectionHalfLength-1,从0开始存放坐标系上r+180和r的直线投影的数据数组第二维的中心是标号为nProjectionHalfLength 的点，对应过原点的直线intnProjectionHalfLength = (nWidth>nHeight)? nWidth:nHeight;nProjectionHalfLength = (int)(nProjectionHalfLength+1)/2.0);for (intk = 0; k <nProjectionHalfLength; k+)-可编辑修改-(for (inti = 0;i <

24、 180; i+)/ 进行投影，每一度一次(/法线长度转换到对应坐标系中的数值floatR = (float)k;R = R/X; /法线长度，坐标转换的尺度应该相同，这里选择 X进行转换，Y也样floatrl = i/180.0*3.1415; / 把角度转换为弧度floatr2 = (i+180)/180.0*3.1415;for (intm = 0; m < nHeight; m+)/遍历图像 行(for (intn = 0; n < nWidth; n+) /歹!J(/进行坐标转换floatx = (n - X)/X;floaty = (Y - m)/Y;if (fabs(

25、x*cos(r1)+ y*sin(r1)-R)<0.005)(根据i与r, k与R之间的关系，可明确每一个数据对应的像素直线pProjDanProjectionHalfLength-1+ki=pProjDanProjectionHalfLength+ki + 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之间的关系，可明确每一个数据对应的像素直线pProjDanProjecti

26、onHalfLength-1-ki=pProjDanProjectionHalfLength-ki + pBits(nHeight -1-m)*nWidthBytes +n;)/for (int n =/for (int m = 0/for (float i = 0;/for (int k这种思路的优点在于，思路简单，比较精确，用阈值0.005可以让计算机相对准确找到在对应直线上的像素。缺点在于，四重循环，计算量太大，测试发现，该函数在调试状态下的运行时间大约是10分钟。换一种思路，不必让计算机通过遍历整个图像来计算哪些点在该条直线上，根据直线方程，我们能够通过遍历x轴(相当于得到列的信息)来

27、计算 y轴(对应的行的信息)，则可以得到在原图像坐标系下对应像素的坐标，这样计算一条直线就可以只需要三重循环，速-可编辑修改-o度大大加快。值得注意的是，如果我们要在像素区中操作某一个像素的灰度值，在需要得到在图像坐标系下的该像素的坐标，如果我们要计算判断直线，则应该在空间坐标系下进行。代码如下：BOOLCCTEmulate:GetProjectionData(double *pProjDa)(/获得图像的投影数据/以下的坐标系仍采用高上凯医学成像系统附录II中坐标系/ProjDa数组用于存放采样的数据，大小应为 m*n m为360 (角度)，n为图像 长宽最大值的一半ASSERT(m_pim

28、ageproject != NULL);/获得图像相应结构数据intnWidthBytes;DWORDnColorNum;char *pBuffer = (char *)m_pimageproject->GetDIBPointer(&nWidthBytes);if (!pBuffer) returnFALSE;nColorNum = m_pimageproject->GetNumC010rs();unsignedchar* pBits;-可编辑修改-pBits = (unsignedchar* )&pBuffersizeof(BITMAPFILEHEADER)+si

29、zeof(BITMAPINFOHEADER)+nColorNum*sizeof(RGBQUAD);intnWidth = m_pimageproject->GetWidth();intnHeight = m_pimageproject->GetHeight();/进行坐标转换floatX = nWidth/2.0;floatY = nHeight/2.0;/为了在数组中将数据与角度和法线长度(根据直线的法线方程)联系起来角度从0旋转到179 (整数)，法线长度从0图像长和宽的最大值的一半(整数)求每次法线直线方程中法线长度，决定了数组的大小，非常重要！/数组的第二维大小为2*nPr

30、ojectionHalfLength-1,从0开始存放坐标系上r+180和r的直线投影的数据数组第二维的中心是标号为nProjectionHalfLength-1 的点，对应过原点的直线nProjectionHalfLengthintnProjectionHalfLength = (nWidth>nHeight)? nWidth:nHeight;=(int)(nProjectionHalfLength+1)/2.0);-可编辑修改- intnProjectionLength = 2*nProjectionHalfLength -1;for (inti = 0;i < 180; i+

31、)/ 进行投影，每一度一次(floatr1 = i/180.0*3.1415; /把角度转换为弧度floatr2 = (i+180)/180.0*3.1415;for (intk = 0; k <nProjectionLength; k+)(floatL = float(k - nProjectionHalfLength + 1); /法线长度intn = 0;/原图像的列，从0开始计算while (k<nProjectionHalfLength-1 && n<nWidth) /r2(floatx = n-X;floaty = (L - x*cos(r2)/s

32、in(r2); /由方程和图像列的信息推算行的信息intm = int(Y - y); /取整有一定误差，原图像的行，也从 0开始计算n+;if (m<0 | m>=nHeight) continue;如果m超出图像范围，则跳过pProjDai*nProjectionLength + k = pProjDai*nProjectionLength + k + pBits(nHeight-1-m)*nWidthBytes + n;n = 0;/对每一条直线x轴都要扫描一次-可编辑修改-while (k>=nProjectionHalfLength-1 && n&l

33、t;nWidth) /r1 (floatx = n-X;floaty = (L - x*cos(r1)/sin(r1);intm = int(Y - y);/取整有一定误差，原图像的行，也从 0开始计算n+;if (m<0 | m>=nHeight) continue; / 如果m超出图像范围，则跳过pProjDai*nProjectionLength + k = pProjDai*nProjectionLength + k + pBits(nHeight-1-m)*nWidthBytes + n;)/for (float i = 0;/for (int k实验发现，只使用三重循环

34、速度提高了十倍以上，不需要长时间等待。采用的是后面一种采集投影值的方法。这样，所有的投影值便放到了数组中。5 .卷积反投影卷积需要解决的问题是：滤波器的采样率应该多大？根据空间坐标系和归一化坐标系之间的关系，可以看出，原图像的采样率为2/width-可编辑修改-滤波器的采样率应该是图像的采样率的2倍以上，所以，对于滤波器的采样率至少为 1/width 。计算得滤波器之后，按照卷积公式进行计算即可。值得注意的是，投影数据在某一个角度下共只有nProjectionLength 个点，所以在卷积时，只需要计算这nProjectionLength 个点的 值即可。卷积代码如下：/滤波器的半长/滤波器采

35、样率/滤波器的长个点。而变量pfilter中序一号是对应的是中心在0的滤波器各点的原/设计滤波器intFilterHalfLength = max(nWidth,nHeight);doubledeltaN = 1.0/FilterHalfLength;intFilterLength = 2*FilterHalfLength+1;double *pFilter = newdoubleFilterLength;nFliterSelect = FilterSedlg.m_nSelectFlag;滤波器的中心在0,左右各巾terhalfLength0filterlength-1/ 所以,i-filte

36、rhalflength坐标if (nFliterSelect = 0) /RL filterfor (inti = 0;i < FilterLength; i+)if (i-FilterHalfLength) = 0)/ 对应于正中点pFilteri = 1.0/4.0/deltaN/deltaN;continue;-可编辑修改-/对应于偶数if (i -FilterHalfLength)%2 = 0) ( pFilteri = 0 ;continue;)if (i-FilterHalfLength)%2 != 0) (pFilteri=-1.0/(i-FilterHalfLength)

37、/(i-FilterHalfLength)/3.1415/3.1415/deltaN/delta N;)else(for (inti = 0;i<FilterLength;i+)(pFilteri=-2.0/3.1415/3.1415/deltaN/deltaN/(4*(i-FilterHalfLength)*(i-FilterHalfLe ngth)-1);-可编辑修改-)/开始卷积for (intr = 0; r < 180; r+)/每一个角度的一组数据卷积一次(for (intn = 0; n < nProjectionLength; n+)/计算结果的 nproje

38、ctionlength个点(doublesum = 0;for (intk = 0; k< nProjectionLength; k+) (sum = sum +pProjectionDatar*nProjectionLength+k*pFilterFilterHalfLength + k - n;)pnewProjectionDatar*nProjectionLength+n = sum;)反投影的方法在投影数据的安排下非常简洁，根据投影坐标中与直线参数和之间的关系，便可以计算得对应直线， 在图像上找到对应直线的像素点， 将投影值加在该像素点 的灰度值上即可，这是一个与获得投影值操作的一个相反操作。for ( r = 0; r < 180; r+)/ 角度-可编辑修改-(doublerl = r/180.0*3.1415;doubler2 = (r + 180.0)/180.0*3.1415;for (intk = 0; k < nProjectionLength; k+)/