数字图像处理实验报告

《医学图像处理》实验报告

实验八：滤波反投影重建

摘要

本次实验的实验目的及主要内容是滤波反投影重建,实验目的包括以下几点: ? 了解图像投影的原理； ? 认识Radon变换；

? 了解反投影重建图像的原理； ? 认识逆Radon变换；

? 了解实现逆Radon变换的方法。

1

一、技术讨论

1.1实验原理

1.图象投影原理:投影变换（projection transformation）是将一种地图投影点的坐标变换为另一种地图投影点的坐标的过程。大致示意图如下：

2.Radon变换原理:radon变换大致可以这样理解：一个平面内沿不同的直线（直线与原点的距离为d，方向角为

alfa）对f（x，y）做线积分，得到的像F(d,alfa)就是函数f的Radon变换。也就是说，平面（d，alfa）的每个点 的像函数值对应了原始函数的某个线积分值,大致可以表示为下图:

3.反投影重建图像的原理： 反投影重建图像过程可总结为： （1）采集不同视角下的投影；

（2）求出各投影的一维傅氏变换，此即图像二维傅氏变换过原点的各切片，理论上是连续的无穷多片；

（3）将上述各切片汇集成图像的二维傅氏变换； （4）求二维傅氏反变换的重建图像。 4.逆Radon变换原理：

有限Radon变换的逆变换可以通过有限逆投影变换FBP（Finite Back Projection）来得到：

其中Pij指的是所有通过点(i,j)的直线的斜率k 和截距l 的集合，即：

2

为了获得更好的能量集中性，有限Radon变换（FRAT）和积分分变换FBP要求变换的图像均值为0，对于均值不为零的图像可以在变换前先减去均值，以保证变换前的图像均值为零；反变换回来后再加上图像均值即可恢复原图像。

1.2实验方法

（这里主要叙述一下你在这次实验所用的代码中所学习到的新函数什么的，特别是关于Mat类的一些矩阵操作函数，并讲解一下那些函数的作用，如：

1、Mat flatten(Mat Arr_In)——将矩阵转置并拉平，用法如下： {

Arr_In = Arr_In.t();

Mat Arr_Out = Arr_In.reshape(1,1); return Arr_Out; }

2、Mat Design_filter(Mat& In,int a,int d)——设计滤波器函数，用法如下：

{

int nextpow2 = ceil(log(2*In.rows)/log(2)); int new_len = 2*pow(2,nextpow2); int order = max(64,new_len/2); Mat filt = (Mat_(1,order)); ...(选择相关参数来选择相应的滤波器)

3

}

4

二、结果与讨论

2.1实验结果

图2.1（其中Div=10

a

b

c d

无滤波作用，即a=0）

（a）原图 （b）反向投影结果 （c）图像正弦图 （d）滤波后投影图像正弦图

a

b

5

c d

图2.2（其中Div=100 无滤波器作用，即a=0）

图2.3（其中Div=360 （a）原图 （b）反向投影结果 （c）图像正弦图 （d）滤波后投影图像正弦图

a b

c d

shepp-logan滤波器，即a=1）

（a）原图 （b）反向投影结果 （c）图像正弦图 （d）滤波后投影图像正弦图

a

b

6

c d

图2.3（其中Div=360 hamming滤波器，即a=2）

图2.4（其中Div=360 （a）原图 （b）反向投影结果 （c）图像正弦图 （d）滤波后投影图像正弦图

a b

c d

hann滤波器，即a=3）

（a）原图 （b）反向投影结果 （c）图像正弦图 （d）滤波后投影图像正弦图

a

b

7

c d

图2.4（其中Div=360 无滤波器，即a=0）

（a）原图 （b）反向投影结果

（c）图像正弦图 （d）滤波后投影图像正弦图

2.2实验讨论

2.2.1 投影参数的修改：

对比分析一下实验指导中第5个步骤修改参数后的结果，为什么会出现这样的结果？

对比图2.1和2.2可知，对于取不同的Div值，其值越大，得到的反投影图像更清晰。根据代码，Div的值决定了对角度分割的程度，Div值越大，扫描的时候角度增量越小，则扫描的精确，得到的图像更清晰。

2.2.2 滤波器的选择：

1）对比分析滤波前后的结果图像，并说明原因；

对比图2.2和2.3结果可知，当加上滤波器之后，得到的图像更清晰，没有可见的模糊。原因是未经过滤波的图像中存在很多干扰噪声成分，而经过滤波后再进行反投影操作能有效除去这些噪声和模糊，于是得到了清晰的图像。

2）对比选择不同的滤波器的滤波效果，说明使用哪个滤波器效果比较好，不用说明原因。

由以上图像显示结果，hann窗滤波器和hamming窗滤波效果没有明显区别，不同的图像适合不同的滤波器，而本次试验中，从视觉上来说，shepp-logan滤波器得到的图像更清晰。

8

附录（实验代码）

.pro文件代码如下：

QT += core gui

greaterThan(QT_MAJOR_VERSION, 4): QT += widgets TARGET = shiyanba TEMPLATE = app

SOURCES += main.cpp\\ mainwindow.cpp

INCLUDEPATH+=C:\\Qt\\opencv2.2\\include\\opencv\\ C:\\Qt\\opencv2.2\\include\\opencv2\\ C:\\Qt\\opencv2.2\\include

LIBS+=C:\\Qt\\opencv2.2\\lib\\libcv.dll.a\\

C:\\Qt\\opencv2.2\\lib\\libopencv_calib3d220.dll.a\\ C:\\Qt\\opencv2.2\\lib\\libopencv_contrib220.dll.a\\ C:\\Qt\\opencv2.2\\lib\\libopencv_core220.dll.a\\

C:\\Qt\\opencv2.2\\lib\\libopencv_features2d220.dll.a\\ C:\\Qt\\opencv2.2\\lib\\libopencv_flann220.dll.a\\ C:\\Qt\\opencv2.2\\lib\\libopencv_gpu220.dll.a\\ C:\\Qt\\opencv2.2\\lib\\libopencv_highgui220.dll.a\\ C:\\Qt\\opencv2.2\\lib\\libopencv_imgproc220.dll.a\\ C:\\Qt\\opencv2.2\\lib\\libopencv_legacy220.dll.a\\ C:\\Qt\\opencv2.2\\lib\\libopencv_ml220.dll.a\\

C:\\Qt\\opencv2.2\\lib\\libopencv_objdetect220.dll.a\\ C:\\Qt\\opencv2.2\\lib\\libopencv_video220.dll.a

Main.Cpp代码如下：

#include

#include #include #include using namespace cv; using namespace std;

//子函数---创建亚像素图像

Mat subpixels(Mat source_image) {

Mat target_image(source_image.rows*2, source_image.cols*2, CV_32F); Mat Matrix_In = (Mat_(2,2) << 1,1,1,1); int source_rows = source_image.rows,

source_cols = source_image.cols; Mat mTemp;

for ( int i = 0; i < source_rows; i++) { for ( int j = 0; j < source_cols; j++) {

9

float pixel_quarter = (float)(source_image.at(i, j)); mTemp = Matrix_In.mul(Matrix_In, pixel_quarter); Mat t=target_image(Rect_(j*2, i*2, 2, 2)); mTemp.copyTo(t); } }

cout<<\ return target_image; }

//子函数---矩阵转置并拉平

Mat flatten(Mat Arr_In) { Arr_In = Arr_In.t();

Mat Arr_Out = Arr_In.reshape(1,1); return Arr_Out; }

//子函数---对坐标进行采样

void ndgrid(int xc,int yc, int m,int n, Mat& X,Mat& Y) { Mat M = (Mat_(2*m,1)); for(int i=0;i<2*m;i++) {

M.at(i,0) = 0.75-xc+(float)i/2; }

Mat N = (Mat_(1,2*n)); for(int i=0;i<2*n;i++) {

N.at(0,i) = 0.75-yc+(float)i/2; }

X = repeat(M,1,2*m); Y = repeat(N,2*n,1); }

//子函数---矩阵元素累加

void accumarray(Mat subs, Mat val, Size sz, Mat& RT){ RT = Mat::zeros(sz.width, sz.height, CV_32F); for(int i=0; i

RT.at(subs.at(i)-1) += val.at(i); } }

//子函数---设计滤波器

Mat Design_filter(Mat& In,int a,int d){

int nextpow2 = ceil(log(2*In.rows)/log(2));

10

int new_len = 2*pow(2,nextpow2); int order = max(64,new_len/2); Mat filt = (Mat_(1,order));

//a=0,选择无滤波

if(a==0) {

for(int i=0; i(0,i) = 1; }

filt = filt.t(); return filt; }

for(int i=0; i

filt.at(0,i) = 2*(float)i/order; }

Mat w = (Mat_(1,order/2+1)); for(int i=0; i

w.at(0,i) = 2*M_PI*(float)i/order; }

Mat filt2 = (Mat_(1,order)); filt2.at(0,0) = 0;

//shepp-logan if(a==1) {

for(int i=1; i

filt2.at(0,i) = filt.at(0,i)*(sin(w.at(0,i)/(2*d))/(w.at(0,i)/(2*d))); }

for(int i=order/2+1; i(0,i) = filt.at(0,order-i)*(sin(w.at(0,order-i)/(2*d))/(w.at(0,order-i)/(2*d))); } }

//hamming if(a==2) {

for(int i=1; i

filt2.at(0,i) = filt.at(0,i)*(0.54+0.46*cos(w.at(0,i)/d)); }

for(int i=order/2+1; i

filt2.at(0,i) = filt.at(0,order-i)*(0.54+0.46*cos(w.at(0,order-i)/d)); }

11

}

//hann if(a==3) {

for(int i=1; i

filt2.at(0,i)=filt.at(0,i)*(1+cos(w.at(0,i)/d))/2; }

for(int i=order/2+1; i

filt2.at(0,i) = filt.at(0,order-i)*(1+cos(w.at(0,order-i)/d))/2; } }

filt2 = filt2.t(); return filt2; }

//子函数---傅立叶变换

void DFT(Mat& srcArr, Mat& dstArr) {

Mat planes[] = {Mat_(srcArr), Mat::zeros(srcArr.size(), CV_32F)}; merge(planes, 2, dstArr); dft(dstArr, dstArr); }

//子函数---傅立叶反变换

void inverseDFT(Mat& dft_result, Mat& dst) {

dft(dft_result, dst, DFT_INVERSE+DFT_SCALE|DFT_REAL_OUTPUT); normalize(dst, dst, 0, 1, NORM_MINMAX); }

//子函数---滤波处理

void rho_filter(Mat& In,Mat& Out,int a,int d) { Mat Filt;

Filt = Design_filter(In,a,d);

Mat P = Mat::zeros(Filt.rows, In.cols, CV_32F); Mat target = P(Rect(0,0,In.cols,In.rows)); In.copyTo(target);

Mat P_DFT; DFT(P, P_DFT);

Mat planes[2];

split(P_DFT, planes);

for(int i=0;i

12

for(int j=0;j

planes[0].at(j,i)=planes[0].at(j,i)*Filt.at(j,0); planes[1].at(j,i)=planes[1].at(j,i)*Filt.at(j,0); } }

Mat L;

merge(planes, 2, P_DFT); inverseDFT(P_DFT,L);

Mat target2=L(Rect(0,0,In.cols,In.rows)); target2.copyTo(Out); }

//子函数---矩阵扩展

void repmat(Mat input, int x, int y, Mat& output) { int r = input.rows, c = input.cols;

Mat Matrix_target = Mat::zeros(r*x, c*y, CV_32F); for ( int i = 0; i < x; i++) { for ( int j = 0; j < y; j++) {

Rect roi = Rect(j*c, i*r, c, r);

Mat Matrix_temp = Matrix_target(roi); input.copyTo(Matrix_temp); } }

output = Matrix_target; } //主函数

int main() {

Mat source_image = imread(\ int m = source_image.rows,

n = source_image.cols;

int xc = floor((m+1)/2), yc = floor((n+1)/2);

Mat d = subpixels(source_image);

int b = ceil(sqrt(m*m+n*n)/2+1); Mat xp = (Mat_(1,2*b+1)); for(int i=0;i<2*b+1;i++) { xp.at(0,i) =i-b; }

Size sz=xp.size();

Mat X,Y;

ndgrid(xc,yc,m,n,X,Y);

13

X=flatten(X); Y=flatten(Y); d=flatten(d);

int numDiv = 360;

float * pTheta = new float[numDiv];

Mat RTtotal=Mat::zeros(xp.cols,numDiv,CV_32F); //生成正弦图

for (int i=0; i

Mat Xp=cos(pTheta[i])*Y-sin(pTheta[i])*X; Mat ip=Xp+b+1;

Mat frac=(Mat_(ip.rows,ip.cols)); Mat k=(Mat_(ip.rows,ip.cols)); for(int j=0;j

k.at(0,j) = floor(ip.at(0,j));

frac.at(0,j) = ip.at(0,j)-k.at(0,j); }

Mat S=1-frac;

Mat val1=Mat_(ip.rows,ip.cols); Mat val2=Mat_(ip.rows,ip.cols);

val1=d.mul(S); val2=d.mul(frac);

Mat RT1, RT2;

accumarray(k,val1,sz,RT1); accumarray(k+1,val2,sz,RT2); Mat RT=(RT1+RT2);

Mat RTtotal_temp=RTtotal.col(i); RT.col(0).copyTo(RTtotal_temp); }

imshow(\

normalize(RTtotal, RTtotal, 0, 1, CV_MINMAX);//归一化处理 imshow(\

Mat filter_img;

rho_filter(RTtotal, filter_img, 0, 1);//第三个参数选择滤波器类型

normalize(filter_img, filter_img, 0, 1, CV_MINMAX); imshow(\

14

Mat p = filter_img.clone();

Mat img = Mat::zeros(600, 600, CV_32FC1); int N = 600;

int center = floor((N + 1)/2); int xleft = -center + 1;

Mat x1 = (Mat_(1,N)); Mat y1 = (Mat_(N,1)); for(int i=0; i

x1.at(0,i) = i + xleft; }

Mat x,y;

repmat(x1,N,1,x);

int ytop = center - 1; for(int j=N-1; j>=0; j--) {

y1.at(N-1-j,0) = j - N + 1 + ytop; }

repmat(y1,1,N,y);

int len = p.rows;

int ctrldx = ceil(len/2);

int imgDiag = 2*ceil(N/sqrt(2))+1; int prs = p.rows, pcs = p.cols; if (prs < imgDiag ) {

float rz = imgDiag - prs;

Mat p0 = Mat::zeros(imgDiag, pcs, CV_32F); Mat target = p0(Rect(0, ceil(rz/2), pcs, prs)); p.copyTo(target); p = p0;

ctrldx = ctrldx + ceil(rz/2); }

//Linear interpolation Mat t, a, proj; int Div = 360;

float * theta = new float[Div]; for (int i=0; i

t = x*cos(theta[i]) + y*sin(theta[i]); a = Mat_(t.rows,t.cols); for(int j=0; j

for(int k=0; k

15

a.at(j,k) = floor(t.at(j,k)); } }

Mat proj1(a+ctrldx+1);

for(int j=0; j

for(int k=0; k(j,k);

proj1.at(j,k) = (float)proj.at(data-1,0); } }

Mat proj2(a+ctrldx);

for(int j=0; j

for(int k=0; k(j,k);

proj2.at(j,k) = (float)proj.at(data-1,0); } }

Mat temp = (t-a).mul(proj1) + (a+1-t).mul(proj2); img = img + temp; }

img = img*M_PI/(2*360);

normalize(img, img, 0, 1, CV_MINMAX);

imshow(\ waitKey(0); }

16

a.at(j,k) = floor(t.at(j,k)); } }

Mat proj1(a+ctrldx+1);

for(int j=0; j

for(int k=0; k(j,k);

proj1.at(j,k) = (float)proj.at(data-1,0); } }

Mat proj2(a+ctrldx);

for(int j=0; j

for(int k=0; k(j,k);

proj2.at(j,k) = (float)proj.at(data-1,0); } }

Mat temp = (t-a).mul(proj1) + (a+1-t).mul(proj2); img = img + temp; }

img = img*M_PI/(2*360);

normalize(img, img, 0, 1, CV_MINMAX);

imshow(\ waitKey(0); }

16

本文来源：https://www.bwwdw.com/article/rykw.html