It's magic. R2017a gives a different result from R2016a/b with the following simple code...

Question

Shogo Muramatsu 2017년 4월 10일

1
링크

이 질문에 대한 바로 가기 링크

https://kr.mathworks.com/matlabcentral/answers/334540-it-s-magic-r2017a-gives-a-different-result-from-r2016a-b-with-the-following-simple-code

답변: Prerana Paravastu 2017년 11월 20일

채택된 답변: Walter Roberson

MATLAB Online에서 열기

R2017a gives a different result from R2016a/b with the following code:

I = magic(16); H = magic(4); imfilter(I,H,'circular')

Has the implementation of imfilter been changed?

It is a huge problem for me.

댓글 수: 4
이전 댓글 2개 표시이전 댓글 2개 숨기기

Adam 2017년 4월 11일

I'm not sure if a bug being improved is a good thing or a bad thing!

Shogo Muramatsu 2017년 4월 11일

편집: Shogo Muramatsu 2017년 4월 11일

This bug is a fatal problem especially for researchers and developers concerning signal processing, and should be fixed immediately.

For a double precision input image, the circular extension may fail when using an even size kernel. Please also see the answer made by myself.

Circular convolution is a special operator because

It is diagonalized by the discrete Fourier transform(DFT), i.e., FFT.
The adjoint operator is also circular convolution with flipped kernel, i.e., circular correlation.

This bug is an obstacle to using these properties on MATLAB.

댓글을 달려면 로그인하십시오.

이 질문에 답변하려면 로그인하십시오.

Answer 1

Walter Roberson 2017년 4월 10일

1
링크

이 답변에 대한 바로 가기 링크

https://kr.mathworks.com/matlabcentral/answers/334540-it-s-magic-r2017a-gives-a-different-result-from-r2016a-b-with-the-following-simple-code#answer_262433

Yes, the implementation did change; I do not see any release notes saying so, though.

I see in R2016b that non-separable filters have two branches: filterSingle2DWithConv() for single(), and filterPartOrWhole() for everything else. filterPartOrWhole() calls mex routines to do the work. You are using double(), so you would get the "everything else" mex routines.

I see in R2017a that non-separable filters have three branches: filterDouble2DWithConv() for double() for 2D arrays; filterDouble2DWithConv() (same routine) for double() with 3D arrays; and filterPartOrWhole() otherwise. (Notice filterSingle2DWithConv is gone and handled by filterPartOrWhole() now.) You are using double(), so you get the new filterDouble2DWithConv() routine. It calls upon conv2() to do the work.

I think you would be justified in filing a bug report on this.

댓글 수: 4
이전 댓글 2개 표시이전 댓글 2개 숨기기

Shogo Muramatsu 2017년 4월 10일

R2016a/b give correct results for the above commands.

Shogo Muramatsu 2017년 4월 11일

MATLAB Online에서 열기

I've confirmed that the single precision routine gives the correct result by the following code:

% Input 2D array
I = magic(16);
% Filter kernel
H = magic(4);
% Circular convolution with double precision
ansimfd = imfilter(double(I),H,'conv','circular');
% Circular convolution with single precision
ansimfs = imfilter(single(I),H,'conv','circular');
% Frequency domain filtering
Iw = fft2(I,16,16);
Hw = fft2(H,16,16);
ansfft2 = circshift(ifft2(Iw.*Hw,16,16),-[2 2]);
% Evaluation of double precision filtering
norm(ansimfd-ansfft2)
% Evaluation of single precision filtering
norm(ansimfs-ansfft2)

Casting the input array to single precision would be a way to avoid the malfunction for the double precision data.

댓글을 달려면 로그인하십시오.

Answer 2

Shogo Muramatsu 2017년 4월 11일

0
링크

이 답변에 대한 바로 가기 링크

https://kr.mathworks.com/matlabcentral/answers/334540-it-s-magic-r2017a-gives-a-different-result-from-r2016a-b-with-the-following-simple-code#answer_262643

편집: Shogo Muramatsu 2017년 4월 12일

MATLAB Online에서 열기

From line 450 to line 451 in imfilter.m, we can see the following codes:

 a = padarray_algo(a, prePadSize, method, padVal, 'pre');
 a = padarray_algo(a, padSize, method, padVal, 'post');

These two lines seem to periodically expand the image boundary. In this second line, it seems that the periodicity is collapsing as a result of postpadding by duplicating the result of prepadding.