What is the Image Processing Toolbox Convention for Input to freqz2() ?

I was trying to learn how to use the function freqz2 in the Image Processing Toolbox. That doc page says that the FIR filter input, h, is "in the form of a computational molecule." What is a "computational molecule?" A quick internet search didn't help much. Is that a term commonly used in image processing?

Anyway, because I don't know what that means, I searched around some more and found on What Is Image Filtering in the Spatial Domain? the statement that "The Image Processing Toolbox™ filter design functions return correlation kernels." So I assumed that h should be a correlation kernel (don't know why the doc page for freqz2 doesn't just say that if that is, in fact, the case).

Aside: the linked doc page on "What is Image Filtering ...." has a mistake. The graphic immediately under the line "Computing the (2, 4) Output of Convolution" says "value of rotated convolution kernel" but it should say "value of rotated correlation kernel."

Define an input FIR filter

rng(100);
h = rand(5);

If that's a correlation kernel with the origin at the center and the upper left corner corresponding to [-2,-2], then it seems to me that its 5x5 frequency response would be

H0 = fftshift(fft2(ifftshift(h)));

In order to get the same result from freqz2(), the input has to be rotated 180 deg

H1 = freqz2(rot90(h,-2),5,5,[1 1]);

isequal(H1,H0)

ans = logical
   1

Inside the code of freqz2 is the comment "% Unrotate filter since FIR filters are rotated." Rotated relative to what?

Can anyone explain why that rot90 is needed on input to freqz2 or where the anlaysis has an incorrect assumption or calculation, or any other misunderstanding on my part?

[Walter Roberson]

I don't know... but "computational molecule" is the output of ftrans2 and fsamp2

[Paul]

Is a computational molecule with a single element a "computational atom"?

[Walter Roberson]

I've heard that the plural of computational molecule, is computational polycule.

Answer 1

[Steve Eddins]

The convention for Image Processing Toolbox filter design functions is for the filter to be specified as a correlation kernel. This convention dates back to early 1993, when version 1.0 of the toolbox was under development. (I joined MathWorks a bit later, at the end of 1993.) I believe that the convention was chosen to be consistent with the MATLAB function filter2, which computed correlation instead of convolution. I don’t know this for sure, but I suspect that filter2 was originally implemented by someone with a computer vision background, instead of a signal processing background, and computer vision papers and textbooks in the 1980s tended to be pretty loose in their definition of 2-D filtering.

Anyway, when I joined MathWorks and took over Image Processing Toolbox development following its initial release, I strongly disliked this convention, but I felt at the time that I couldn’t change it, and so related toolbox functions continued to use correlation as the default convention.

The term “computational molecule” was initially used in the toolbox version 1.0 documentation, and it just means correlation kernel. I don’t know where the term originated. I have never seen it used elsewhere.

I believe you are correct that there is an error in the description of convolution on the What Is Image Filtering in the Spatial Domain documentation page. I would probably recommend the following edits:

• Change “and the correlation kernel is” to “and the convolution kernel is”.
• Change step 1 to “Rotate the convolution kernel 180 degrees about its center element.”
• Change step 2 to “Slide the center element of the rotated convolution kernel …”

Regarding the code comment you found inside freqz2 — the frequency response of a filter is the Fourier transform of the impulse response, which for an FIR filter is the convolution kernel. So, the “unrotate the filter” comment refers to the fact that the input is expected to be correlation kernel, which is a rotation of the convolution kernel.

This is all pretty messy, conceptually, and I really wish I had been sufficiently bold in 1993 to go ahead and change the convention.

I retired from MathWorks in March, so I can no longer just stroll down the hall to talk with one of the writers on the toolbox team. But I will submit two service requests:

• Get rid of the “computational molecule” terminology in the toolbox documentation.
• Fix the convolution description on the What Is Image Filtering in the Spatial Domain page.

[Paul]

Hi Steve,

Thanks for the informative answer. Based on what you've said, I think I'm misunderstanding the terms "correlation kernel" and "convolution kernel." I'm infinitely more famiilar with 1D processing than 2D, so I'm going to use a 1D example (which I hope are just a special case of 2D) to try to clear up my confusion.

Define a 1D FIR filter using standard Matlab conventions:

b = 4:-1:1; a = 1;

Now I apply filter to a finite duration sequence

x = 1:5;
filter(b,a,x)

ans = 1x5

4 11 20 30 40

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

This output can also be determined as the first 5 elements of the convolution of b and x

conv(b,x)

ans = 1x8

4 11 20 30 40 26 14 5

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

So, b is a convolution kernel and the input, b, to filter (with a = 1) is a convolution kernel and filter implements convolution, not correlation.

If correct, that sounds like the exact opposite of the convention as you described for filter2, which explains why freqz2 does what it does. Now that I understand (I think) the term "convolution kernel," the freqz2 result makes sense. As you said "the frequency response of a filter is the Fourier transform of the impulse response, which for an FIR filter is the convolution kernel." (I'd argue that statement is true for an IIR filter as well, even though for an IIR filter we can't list all of the values of the impulse response).

So if I wanted to get the same result with filter2, I have to flip b

filter2(fliplr(b),x,'full')

ans = 1x8

4 11 20 30 40 26 14 5

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

and it turns out I have to explicitly specify the 'full' shape, and, unlike filter, filter2 keeps going past the edge of the input.

Verifying agiain that filter2 performs correlation, not convolution, we see that the output of filter2

filter2(b,x,'full')

ans = 1x8

1 4 10 20 30 34 31 20

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

can also be obtained from conv2 after flipping the kernel.

conv2(fliplr(b),x)

ans = 1x8

1 4 10 20 30 34 31 20

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Interesting that filter2() and filter() are fundamentally different, though maybe that just reflects a difference between the image processing and signal processing communities. Seems like there would be a lot of confusion if the Signal Processing Toolbox were ever to be extended to 2D.

One thing is still bothering me (I'll have to think about it some more). In my original example in the Question, I explicitly stated that h(1,1) = h[-2,-2] (parends indicate Matlab indexing, brackets indicate standard signal processing notation), and that seemed to be confirmed by the need to use ifftshift/fftshift to match the result of freqz2(). In other words, the 2D FIR fiter represented by h on input to freqz2() is non-causal (is that term also used in image processing)? It's not clear to me, yet, what filter2() assumes about causality of the impulse response and the input. But filter() clearly assumes that the filter represented by b(z^1)/a(z^-1) is causal and the input sequence, x, is also causal with x(1) = x[0] (on further reflection, filter() doesn't really care about causality of x[n] insofar as filter() doesn't doesn't know, nor care, about the values of n that correspond to the values of x). The fact that I can use filter2() to recreate the output of filter() is making me wonder if filter2() and freqz2() are using the same assumptions about H.

Looking a bit closer at the doc page filter2, it states that the input, H, is "Coefficients of the rational transfer function," which is technically correct (the best kind of correct) but I'm not used to hearing a FIR filter being described as a rational transfer function.

[Steve Eddins]

Hi Paul,

Your 1-D example and code all look good to me.

Interesting that filter2() and filter() are fundamentally different, though maybe that just reflects a difference between the image processing and signal processing communities.

I came to image processing through a graduate department program specializing in digital signal processing, and so I consider myself to be a member of both communities. I do not associate correlation-based 2-D filtering with the image processing community. I do not know for sure why correlation was chosen for filter2. I speculated in my answer that perhaps someone with a computer vision background was involved in the design of filter2. In the 1980s and 1990s, a lot of work in computer vision was being done by people with computer science backgrounds, and so the terminology and concepts were not always consistent with the EE-based signals and systems domain. In particular, I had the impression that filtering was not as precisely defined in the computer vision literature back then as it was in the signal processing community.

Seems like there would be a lot of confusion if the Signal Processing Toolbox were ever to be extended to 2D.

Extending Signal Processing Toolbox to 2-D seems unlikely to me, based on the history of product development of MATLAB, Signal Processing Toolbox, and Image Processing Toolbox. However, there has been an effort over some time (the past 10 years or so?) to move basic 1-D signal analysis and filtering capability from the toolbox into MATLAB. In some cases, 2-D versions of those capabilities have also been added to MATLAB. In most cases, the associated functional designs (function names, syntaxes, behaviors, etc.) have been new. Check out the release notes for MATLAB in the math and data analysis areas going back several years.

... non-causal ...

The term "non-causal" is not commonly used in image processing or computer vision because the filtering domain is space, not time. In image processing and computer vision, it almost always makes the most sense to avoid introducing a phase shift in the filtering process. You wouldn't want a smoothing filter, for example, to shift a picture to the right. That is why filter2 (and the Image Processing Toolbox function imfilter) assume that the filter's origin is the center element, not the upper-left element.