数字图像处理实验报告(3)

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 …… 此处隐藏：2534字，全部文档内容请下载后查看。喜欢就下载吧 ……