Low pass Butterworth filter shoud have no imaginary part

[epitalon]

Hello Toby,

Line 157 in butterworth_lowpass.cpp we read Mat toMerge[] = { tmp, tmp }; I believe that we should read Mat toMerge[] = { tmp, Mat::zeros(tmp.size(), CV_32F) };

as when one multiplies a DFT image by a filter, this filter should be real only, otherwise the multiplication changes the phase along with the magnitude of each pixel in the DFT. In the current case, multiplying by {tmp, tmp} adds pi/4 to the phase in each pixel of the DFT and multiplies its magnitude by sqrt(2)*tmp. Regards, Jean-Marie

[tobybreckon]

Hi Jean-Marie,

Thanks for this. You are of course correct and I have to say this bug has been there a long time ... very well spotted.

I think it is also present in several older examples here: https://github.com/tobybreckon/c-examples-ipcv/

... and has been carried forward to my current python teaching examples here: https://github.com/tobybreckon/python-examples-ip

Can I just check - have you verified that my current output is out by pi/4 in terms of phase and magnitude is increased by a factor of sqrt(2)*tmp in my code or another re-implementation of this logic? The reason I ask is that when I do inverse DFT and display in lines 289-310 I get a meaningful image which I can't immediately imagine I would if the bug was messing up the phase so badly. However, it may explain my need for normalisation in lines 295/296.

I think what may be happening is that cv::mulSpectrums() in convolution mode (i.e. last param conjB = false)  is covering up for my bug by correctly handling the 2-channel magnitude/phase real+imaginary? Looking at https://docs.opencv.org/3.4/d2/de8/group__core__array.html#ga3ab38646463c59bf0ce962a9d51db64f it says "When the arrays are complex, they are simply multiplied (per element) with an optional conjugation of the second-array elements. " Otherwise, I can't explain how I get a meaningful image out.

Another query you may be able to answer on this fix - if the 2-channel filter multiplication mulSpectrums() indeed changes the phase along with the magnitude of each pixel in the DFT by multiplying per element - should your fix not be: Mat toMerge[] = { tmp, tmp };

to:

Mat toMerge[] = { tmp, Mat::ones(tmp.size(), CV_32F) };

in order to preserve the phase information, correctly? Element-wise multiplication by zeros would destroy it, surely? I think the phase information is being carried over unchanged source to target.

T.

[tobybreckon]

I did a quick test with and without this change on what was the opencv4 branch (because it built against opencv 4.7 I have installed locally) and it seemed to make no difference to the output ...

The bug you identified is there in the creation of the filter and your statement on the correct theory "when one multiplies a DFT image by a filter, this filter should be real only, otherwise the multiplication changes the phase along with the magnitude of each pixel in the DFT." is correct - but I now think cv::mulSpectrums() is handling this issue somehow.

T.

[tobybreckon]

.. but using: Mat toMerge[] = { tmp, Mat::ones(tmp.size(), CV_32F) }; doesn't appear to make a difference either:

T.

[epitalon]

Hi Toby,

Can I just check – have you verified that my current output is out by pi/4 in terms of phase and magnitude is increased by a factor of sqrt(2)*tmp in my code or another re-implementation of this logic?

No, it happened that my fix makes a difference only when cutting frequency is extremely low : I chose a radius of 0.01.In theory, I am right for the fix but I did not notice any change for higher radius, say radius >= 1I did not go into details and did not verify that the phase is shifted from pi/4 when using {tmp, tmp}. Neither verified the magnitude. I have no time to do that currently. Actually, I process a color image by processing each channel individually. And for each channel, I call the same function to compute the same filter. Could it be that each of the three times I compute Butterworth_filter(radius=0,01) I obtain different results due to computer imprecision ? What is sure for me, I get visually better results for radius=0.01 when I use my fix. But my application of Butterworth filter with radius=0.01, I reckon, is suspicious. I get no perfect results with my fix either in that the smoothing I get is more intense in direction NE-SW than NW-SE. In other words, smoothing is not isotrope. For greater radius, it is. I do not know what is the cause of this !

should your fix be:  Mat toMerge[] = { tmp, Mat::ones(tmp.size(), CV_32F) }; in order to preserve the phase information, correctly?

No, DFT(x,y) is a complex number and if we multiply it with a complex number which has an imaginary part, this other complex number has a phase which is not null. c = Mag . exp(i.phase) = Mag cos(phase) + i.Mag.sin(phase) = a + i.bphase = 0 <==> b = 0

When multiplying DFT(x,y) by c, in order to not change the phase of DFT(x,y), we must have phase(c) = 0 as when multipying two complex numbers, phases add.

... branch OpenCV 4.0

I do use OpenCV 3.1 but the same fix should also apply to branch 4.0, I guess.

Best regards, Jean-Marie

[epitalon]

Toby,

... and has been carried forward to my current python teaching examples here: https://github.com/tobybreckon/python-examples-ip

I looked in your python example code, but it is different than C++ version. create_butterworth_low_pass_filter() returns a matrix of real numbers. No complex. So the issue is no there.

Jean-Marie

[tobybreckon]

No, DFT(x,y) is a complex number and if we multiply it with a complex number which has an imaginary part, this other complex number has a phase which is not null. c = Mag . exp(i.phase) = Mag cos(phase) + i.Mag.sin(phase) = a + i.bphase = 0 <==> b = 0 When multiplying DFT(x,y) by c, in order to not change the phase of DFT(x,y), we must have phase(c) = 0 as when multipying two complex numbers, phases add.

OK agreed, so this must be how cv::mulSpectrums() is handling it then - despite its documentation which says "When the arrays are complex, they are simply multiplied (per element) ..." it performs complex number multiplication correctly. I tried to look at the OpenCV source quickly but failed to confirm in the time I had available.

I change to:

Mat toMerge[] = { tmp, Mat::zeros(tmp.size(), CV_32F) };

T.

[tobybreckon]

Toby, ... and has been carried forward to my current python teaching examples here: https://github.com/tobybreckon/python-examples-ip I looked in your python example code, but it is different than C++ version. create_butterworth_low_pass_filter() returns a matrix of real numbers. No complex. So the issue is no there. Jean-Marie

Thanks, T.

[epitalon]

Tobby,I made some tests with one image and radius = 4I inspected the values of two first pixels of DFT image before and after multiplication with filter and image after inverse DFTThey are here below.Anyway, what I can say is : when one adds pi/4 to the phase of the DFT image,it breaks the symetry of the DFT which is in a real image's DFT. So the DFT is no longer a DFT of a real valued matrix but that of a complex matrix.And one applies inverse DFT, it gets a complex matrix.But one keeps only the real part of this complex matrix.So I am not sure what is the result, but this can explain why the error on the filter has no much effect on the visual result.Here are the values from my debugger (Visual C++)=======================================Without Fix==========Filter ((float )(filter).data)[2(1024+10242048)] 1.00000000 real ((float )(filter).data)[2(1024+10242048)+1] 1.00000000 imaginary ((float )(filter).data)[2(1025+10242048)] 0.941176474 real ((float )(filter).data)[2(1025+10242048)+1] 0.94117647 imaginaryComplex Image ((float )(complexImg).data)[0+04096] 425946720. pixel 0 real ((float )(complexImg).data)[1+04096] 0.00000000 pixel 0 imaginary ((float )(complexImg).data)[2+04096] 12722329.0 pixel 1 real ((float )(complexImg).data)[3+04096] 31198212.0 pixel 1 imaginaryafter mult by filter ((float )(complexImg).data)[0+04096] 425946720. real ((float )(complexImg).data)[1+04096] 425946720. imaginary ((float )(complexImg).data)[2+04096] -17389066.0 real ((float )(complexImg).data)[3+04096] 41336980.0 imaginaryAfter inverse dft ((float )(complexImg).data)[0+04096] 140.713745 real           <== value that is kept ((float )(complexImg).data)[1+04096] 142.061325 imaginary ((float )(complexImg).data)[2+04096] 140.684418 real           <== value that is kept ((float )(complexImg).data)[3+04096] 142.025543 imaginary With fix=========Filter ((float )(filter).data)[2(1024+10242048)] 1.00000000 real ((float )(filter).data)[2(1024+10242048)+1] 0.0000000 imaginary ((float )(filter).data)[2(1025+10242048)] 0.941176474 real ((float )(filter).data)[2(1025+10242048)+1] 0.0000000 imaginaryComplex Image ((float )(complexImg).data)[0+04096] 425946720. real ((float )(complexImg).data)[1+04096] 0.00000000 imaginary ((float )(complexImg).data)[2+04096] 12722329.0 real ((float )(complexImg).data)[3+04096] 31198212.0 imaginaryafter mult by filter ((float )(complexImg).data)[0+04096] 425946720. real ((float )(complexImg).data)[1+04096] 0.00000000 imaginary ((float )(complexImg).data)[2+04096] 11973957.0 real ((float )(complexImg).data)[3+04096] 29363024.0 imaginaryAfter inverse dft ((float )(complexImg).data)[0+04096] 141.387543 real         <== value that is kept ((float )(complexImg).data)[1+04096] 0.673773646 imaginary ((float )(complexImg).data)[2+04096] 141.354950 real         <== value that is kept ((float )(complexImg).data)[3+04096] 0.670561910 imaginaryJean-Marie et Anne-Marie Epitalon3, rue du Four11270 FanjeauxTel: 04 68 76 90 07envoyé : 24 janvier 2023 à 12:02de : Toby Breckon @.>à : tobybreckon/cpp-examples-ipcv @.>cc : jean marie epitalon @.>, Author @.>objet : Re: [tobybreckon/cpp-examples-ipcv] Low pass Butterworth filter shoud have no imaginary part (Issue #1) Toby, ... and has been carried forward to my current python teaching examples here: https://github.com/tobybreckon/python-examples-ip I looked in your python example code, but it is different than C++ version. create_butterworth_low_pass_filter() returns a matrix of real numbers. No complex. So the issue is no there. Jean-MarieThanks, T.—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you authored the thread.Message ID: @.***>

[epitalon]

Tobby,I hope this helps (to show there was really a mistake). '-)Jean-Marie et Anne-Marie Epitalon3, rue du Four11270 FanjeauxTel: 04 68 76 90 07envoyé : 24 janvier 2023 à 12:02de : Toby Breckon @.>à : tobybreckon/cpp-examples-ipcv @.>cc : jean marie epitalon @.>, Author @.>objet : Re: [tobybreckon/cpp-examples-ipcv] Low pass Butterworth filter shoud have no imaginary part (Issue #1) Toby, ... and has been carried forward to my current python teaching examples here: https://github.com/tobybreckon/python-examples-ip I looked in your python example code, but it is different than C++ version. create_butterworth_low_pass_filter() returns a matrix of real numbers. No complex. So the issue is no there. Jean-MarieThanks, T.—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you authored the thread.Message ID: @.***>

[epitalon]

Hi Tobby,I want to show the difference between your version and my fix.I will send you four images. - original (in this mail)- filtered with radius = 4 with original software - filtered with fixed software- difference between the last two images (in next mail)We can see there is a small shift in the pattern between the two filtered images;They are colour images but we could also take only one chanel, green for exemple.Jean-Marie et Anne-Marie Epitalon3, rue du Four11270 FanjeauxTel: 04 68 76 90 07envoyé : 24 janvier 2023 à 12:02de : Toby Breckon @.>à : tobybreckon/cpp-examples-ipcv @.>cc : jean marie epitalon @.>, Author @.>objet : Re: [tobybreckon/cpp-examples-ipcv] Low pass Butterworth filter shoud have no imaginary part (Issue #1) Toby, ... and has been carried forward to my current python teaching examples here: https://github.com/tobybreckon/python-examples-ip I looked in your python example code, but it is different than C++ version. create_butterworth_low_pass_filter() returns a matrix of real numbers. No complex. So the issue is no there. Jean-MarieThanks, T.—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you authored the thread.Message ID: @.***>

[epitalon]

Jean-Marie et Anne-Marie Epitalon3, rue du Four11270 FanjeauxTel: 04 68 76 90 07envoyé : 24 janvier 2023 à 12:02de : Toby Breckon @.>à : tobybreckon/cpp-examples-ipcv @.>cc : jean marie epitalon @.>, Author @.>objet : Re: [tobybreckon/cpp-examples-ipcv] Low pass Butterworth filter shoud have no imaginary part (Issue #1) Toby, ... and has been carried forward to my current python teaching examples here: https://github.com/tobybreckon/python-examples-ip I looked in your python example code, but it is different than C++ version. create_butterworth_low_pass_filter() returns a matrix of real numbers. No complex. So the issue is no there. Jean-MarieThanks, T.—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you authored the thread.Message ID: @.***>

[tobybreckon]

Thanks for all your efforts on this, I am convinced there is a difference (although the images never came through the email <-> github link but never mind).

I have made the change you suggested as:

Mat toMerge[] = { tmp, Mat::zeros(tmp.size(), CV_32F) };

into the master branch. Thanks again for finding this long standing bug.

[epitalon]

Toby,it was a pleasure.If you wish that I sent the images to @. , let me know.Jean-Marie et Anne-Marie Epitalon3, rue du Four11270 FanjeauxTel: 04 68 76 90 07envoyé : 25 janvier 2023 à 23:19de : Toby Breckon @.>à : tobybreckon/cpp-examples-ipcv @.>cc : jean marie epitalon @.>, Author @.>objet : Re: [tobybreckon/cpp-examples-ipcv] Low pass Butterworth filter shoud have no imaginary part (Issue #1) Thanks for all your efforts on this, I am convinced there is a difference (although the images never came through the email <-> github link but never mind).I have made the change you suggested as:Mat toMerge[] = { tmp, Mat::zeros(tmp.size(), CV_32F) };into the master branch. Thanks again for finding this long standing bug.—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you authored the thread.Message ID: @.>