实验二CT重建实验报告(共14页).doc

精选优质文档-倾情为你奉上实验二 CT图像重建实验一、实验目的：通过编写CT图像重建程序，进一步熟悉CT重建过程，同时加强图像处理程序的编程训练。二、实验软件： VC+三、实验要求1递交整个程序的执行程序和源程序。2 要求本小组承担部分的算法原理，列出主要程序段并给出相应说明，对实验结果进行分析。四、算法原理及结果分析：1. CT重建原理为：在CT成像中，物体对X线的吸收起主要作用，在一均匀物体中，X线的衰减服从指数规律。在X线穿透人体器官或组织时，由于人体器官或组织是由多种物质成分和不同的密度构成的，所以各点对X线的吸收系数是不同的。将沿着X线束通过的物体分割成许多小单元体（体素），令每个体素的厚度相等(l)。设l足够小，使得每个体素均匀，每个体素的吸收系数为常值，如果X线的入射强度I0、透射强度I和物体体素的厚度l均为已知，沿着X线通过路径上的吸收系数之和12n就可计算出来。为了建立CT图像，必须先求出每个体素的吸收系数1、2、3n。为求出n个吸收系数，需要建立如上式那样n个或n个以上的独立方程。CT成像装置从不同方向上进行多次扫描，来获取足够的数据建立求解吸收系数的方程。吸收系数是一个物理量，它是CT影像中每个像素所对应的物质对X线线性平均衰减量大小的表示。再将图像面上各像素的CT值转换为灰度，就得到图像面上的灰度分布，就是CT影像。CT重建过程可以采用直接反投影和卷积反投影来实现。卷积反投影重建图像时，先把由检测器上获得的原始数据与一个滤波函数进行了卷积运算，得到各方向卷积的投影函数；然后再把它们从各方向进行反投影，即按其原路径平均分配到每一矩阵元上，进行叠加后得到每一矩阵元的CT值；再经过适当处理后就可以得到被扫描物体的断层图像，卷积反投影可消除单纯的反投影产生的边缘失锐效应，补偿投影中的高频成分和降低投影中心密度，并保证重建图像边缘清晰和内部分布均匀。2. 整个程序的结构为：产生sleep logan 模型产生反投影数据卷积反投影编程在photostar平台下进行。编程过程中有三个坐标系：坐标系1：图像坐标系，坐标原点在图像的左上方，x轴水平向右，y轴水平向下，单位为一个“像素”。坐标系2：空间坐标系，坐标原点在图像中心，x轴水平向右，y轴水平向上，单位为一个“像素”。坐标系3：归一化空间坐标系，坐标原点在图像中心，x轴水平向右，y轴水平向上，归一化因子为x轴的一半，所以对于m*m大小的图像，xy轴的坐标范围为-1,1。根据坐标系1和坐标系3，可以确定图像的采样率。根据面向对象的思想，在photostar平台中创建CCTEmulate类来管理实现CT仿真操作。3. 产生sleep logan模型要产生一幅图像，对于DIB来说，需要有信息头、文件头、调色板、像素区等部分，为操作简洁，我们读入一幅空白图像（text.bmp），指针m_pImageObject指向这幅空白图像的像素区，直接对像素区进行操作。值得注意的是，根据实验需求，读入图像为灰度图像即可。用CCTEmulate类中的SheepLogan函数相应菜单消息，进行仿真头模型的产生。需要解决一个问题：如何确定一个像素在哪个椭圆中？根据椭圆参数方程：其中A B分别是椭圆的长轴和短轴，x0，y0分别是椭圆的中心坐标，a是椭圆旋转角度。则根据上式便可以确定在坐标中的一点是否在某一椭圆中。SheepLogan函数的思路为：遍历图像，检测图像的每一个像素点是否依次在某一椭圆中，若在a外，则置为背景色，在其他任何椭圆中，则置为该椭圆对应的灰度值。所有的操作均在归一化坐标系中进行。代码如下：voidCPhotoStarView:OnCteCreatsheeplogan() / TODO: Add your command handler code hereCCTEmulateCCTE;CCTE.m_pimageproject = m_pImageObject;CCTE.SheepLogan();CClientDCdc(this);Invalidate();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; /重新用一个变量来进行以下操作/*获取图像的像素区指针和各种参数*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 *) &pBuffersizeof(BITMAPFILEHEADER);pRGBPalette = (RGBQUAD *) &pBuffersizeof(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)case 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: /椭圆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_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从0180，每10度画线一次，和分别代表过原点的左右两边的直线。每次画线从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模型中，坐标系参见高上凯医学成像系统附录II/在此坐标系中，第（ij）个点的坐标（xy）为：/x = （j-X）/X y = （Y-i）/Y/其中X = m_nWidth/2 Y = m_nHeight/2floatX = 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<nHeight;i+) /遍历整个图像,hang for (intj = 0;j<nWidth;j+) /liefloatx = (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=0产生投影数据主要解决两个问题：每一条扫描线的投影值如何存放？每一条扫描线上的投影值是多少？安排数据结构如下：扫描共计180次，角度为0179，图像长或者宽的最大值的一半设为projectionhalflength，则每次扫描线数为projectionlength = 2*projectionhalflength -1，安排数组为projectiondata180*projectionlength.根据直线的法线方程，则在第r个角度的、法线为的数据存放的位置为projectiondatar*projectionlength + k,其中或者 ， 。其中当时，当是，他们对应不同的直线。中间的一条直线标号为projectionhalflength - 1,对应过原点的直线。对于特定的角度r，则一排扫描线从0projectionlength-1，存放了从到的投影值。根据r和与数组projectiondata之间的关系，可以确定在角度r和法线确定的直线的投影值存放的位置。根据角度和法线，对所有存在的直线进行扫描。思路是：双重循环r = 0179，k = 0projectionlength-1，这样就确定了一条直线的参数和。在双重循环中扫描整幅图像，搜寻这条直线上的像素点（坐标满足直线方程的点就认为在该条直线上），然后将该像素点的灰度值加到对应的投影值中去。代码如下：BOOLCCTEmulate:GetProjectionData(double (*pProjDa)180)/获得图像的投影数据/以下的坐标系仍采用高上凯医学成像系统附录II中坐标系/ProjDa数组用于存放采样的数据，大小应为m*n m为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* )&pBuffersizeof(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 <nProjectionHalfLength; 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 < nHeight; m+)/遍历图像 行 for (intn = 0; n < nWidth; n+) /列 /进行坐标转换floatx = (n - X)/X;floaty = (Y - m)/Y;if (fabs(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之间的关系，可明确每一个数据对应的像素直线pProjDanProjectionHalfLength-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轴（相当于得到列的信息）来计算y轴（对应的行的信息），则可以得到在原图像坐标系下对应像素的坐标，这样计算一条直线就可以只需要三重循环，速度大大加快。值得注意的是，如果我们要在像素区中操作某一个像素的灰度值，在需要得到在图像坐标系下的该像素的坐标，如果我们要计算判断直线，则应该在空间坐标系下进行。代码如下：BOOLCCTEmulate:GetProjectionData(double *pProjDa)/获得图像的投影数据/以下的坐标系仍采用高上凯医学成像系统附录II中坐标系/ProjDa数组用于存放采样的数据，大小应为m*n m为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* )&pBuffersizeof(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 <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)/sin(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<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实验发现，只使用三重循环速度提高了十倍以上，不需要长时间等待。采用的是后面一种采集投影值的方法。这样，所有的投影值便放到了数组中。5.卷积反投影 卷积需要解决的问题是：滤波器的采样率应该多大？ 根据空间坐标系和归一化坐标系之间的关系，可以看出，原图像的采样率为2/width，滤波器的采样率应该是图像的采样率的2倍以上，所以，对于滤波器的采样率至少为1/width。计算得滤波器之后，按照卷积公式进行计算即可。值得注意的是，投影数据在某一个角度下共只有nProjectionLength个点，所以在卷积时，只需要计算这nProjectionLength个点的值即可。卷积代码如下：/设计滤波器intFilterHalfLength = max(nWidth,nHeight); /滤波器的半长doubledeltaN = 1.0/FilterHalfLength; /滤波器采样率intFilterLength = 2*FilterHalfLength+1; / 滤波器的长double *pFilter = newdoubleFilterLength;nFliterSelect = FilterSedlg.m_nSelectFlag;/滤波器的中心在0，左右各fiterhalfLength个点。而变量pfilter中序号是0filterlength-1/所以，i-filterhalflength对应的是中心在0的滤波器各点的原坐标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)/(i-FilterHalfLength)/3.1415/3.1415/deltaN/deltaN; elsefor (inti = 0;i<FilterLength;i+) pFilteri = -2.0/3.1415/3.1415/deltaN/deltaN/(4*(i-FilterHalfLength)*(i-FilterHalfLength)-1); /开始卷积for (intr = 0; r < 180; r+) /每一个角度的一组数据卷积一次for (intn = 0; n < nProjectionLength; n+)/计算结果的nprojectionlength个点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+) /角度doubler1 = r/180.0*3.1415;doubler2 = (r + 180.0)/180.0*3.1415;for (intk = 0; k < nProjectionLength; k+) /长度 intL = abs(k-nProjectionHalfLength+1); /法线长度intn = 0; /原图像的列,从0开始计算while (k<nProjectionHalfLength-1 && n<nWidth) /r2 floatx = n-X; floaty = (L - x*cos(r2)/sin(r2); /由方程和图像列的信息推算行的信息intm = int(Y - y); /取整有一定误差，原图像的行，也从0开始计算n+;if (m<0 | m>=nHeight) continue; /如果m超出图像范围，则跳过pBitsm*nWidth + n = pBitsm*nWidth + n + pnewProjectionDatar*nProjectionLength + k; n = 0; /对每一条直线x轴都要扫描一次while (k>=nProjectionHalfLength-1 && n<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超出图像范围，则跳过pBitsm*nWidth + n = pBitsm*nWidth + n + pnewProjectionDatar*nProjectionLength + k; 值得注意的地方是，反投影后像素区的值很大，因为相当于投影了180次，并且，需要进行最大值和最小值控制。还有，图像的像素区是0255的整型.将整理好的数据放到图像的像素区，可以得到投影后的图像。6.实验结果和分析相关参数： 投影分为180次，每度一次，每次扫描线数为图像长或者宽的最大值，滤波器分辨率为图像分辨率的2倍。以下是实验结果。图一 生成的sleep logan图像图四 直接反投影效果图二 RL滤波器卷积反投影效果图三 SL滤波器卷积反投影效果实验发现，调节滤波器的分辨率对重建图像没有明显改善效果。直接反投影法重建图像各个像素呈现“星状”伪迹，已经将所有有效元素亮化为白色，某一元素有自己的投影值，也同时处在其他元素的“星状”伪迹中，叠加形成亮色。滤波反投影的数据能大致反映出图像原来的特性，但是分辨率比较低，原图像中的h、i、j椭圆几乎无法生成，图像上下边缘靠近中间的区域有白色伪迹。对于SL滤波器，其效果比RL滤波器要好，白色杂点的状况要好于RL滤波器。值得注意的是，观察到图像具有一点对称的特性，出现这种情况的原因可能是程序中存在逻辑错误，不过至今还没查出。专心-专注-专业