DIP 例4.15 空间域滤波得到频率域滤波

做了整整一天才完成这个例子。完成的代码也存在一个BUG，在最后提一下好了。

J=imread('Fig1016(a)(building_original).tif');imshow(J)p=double(J);P=834+2;Q=1114+2;%size(J)transform=ones(P,Q);for i=1:P    for j=1:Q        transform(i,j)=(-1)^(i+j);    endendfilledJ=zeros(P,Q);filledJ(1:size(J,1),1:size(J,2))=p;filledJ=filledJ.*transform;%padded the J & center J2=fft2(filledJ);imshow(1+log(abs(J2)),[])%fft %equal to aboveJ3=fftshift(fft2(J,P,Q));imshow(1+log(abs(J3)),[])%init the filter H in frequencyh=[1  0 -1;2 0 -2;1 0 -1];%sobel operationhp=zeros(P,Q);%padded the h%hp(1:3,1:3)=h;%hp((1+P)/2-1:(1+P)/2+1,(1+Q)/2-1:(1+Q)/2+1)=h;%hp((P-1)/2:(P-1)/2+2,(Q-1)/2:(Q-1)/2+2)=h;hp(P/2:P/2+2,Q/2:Q/2+2)=h;%place the h at the center of the padded h to preserve the odd symmetryH=fft2(hp.*transform);%center when transforming to frequency domainH=H-real(H);%real(H)=0; H=H.*transform;imshow(abs(H),[])%the step 4. h was moved to the center of the padded h when ifft g = real(ifft2(H.*J2)).*transform;% Crop to original size.g = g(1:size(p, 1), 1:size(p, 2));figureimshow(g,[])y=imfilter(p,h);figure imshow(y,[])

代码的主要问题围绕在频率域滤波器的生成上，4.15中提到用0填充Sobel模板，但是这个过程之中仍需要保证填充后的模板hp仍然保持奇对称，为此我们可以设模板头部的位置位于（p,q） 为了保持奇对称f(p,q)=-f(P-p,Q-q)=-f(p+2,q+2) (这里计数是从0开始的)，计算得到的(p,q)+(1,1)就转化为Matlab中的矩阵坐标了。虽然不是很明白为什么之后的step4。



原图像格式不支持，。， 就不传了
BUG：这个生成的滤波后的g和y有一个符号差，我找了半天也没有找到问题出现在哪里。如果有大佬看到了求指教。

最后是书中的源码：

J=imread('Fig1016(a)(building_original).tif');h=[1  0 -1;2 0 -2;1 0 -1];%sobel operationQ=size(J,2)+2;P=size(J,1)+2;p=double(J);H=freqz2(h,Q,P);H1=ifftshift(H2);GF=dftfilt(p,H1);imshow(GF,[])

dftfilt似乎并没有按照要求的来 没有进行center，但是结果是一样的。