{"id":2473,"date":"2017-01-16T11:29:19","date_gmt":"2017-01-16T16:29:19","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=2473"},"modified":"2019-11-01T16:52:22","modified_gmt":"2019-11-01T20:52:22","slug":"aliasing-and-image-resizing-part-3","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2017\/01\/16\/aliasing-and-image-resizing-part-3\/","title":{"rendered":"Aliasing and image resizing – part 3"},"content":{"rendered":"

Today I'll try to wrap up my discussion about how aliasing affects image resizing and about how the imresize<\/tt> function tries to prevent it. (This is called antialiasing<\/i>.) Let me show you where we are going. Here is the zone plate image I showed last time.<\/p>

Z = imzoneplate(501);\r\nimshow(Z)\r\ntitle('Z'<\/span>)\r\n<\/pre>\"\" 
Here's what happens when we shrink Z<\/tt> by throwing away samples.<\/p>
Z4 = Z(1:4:end,1:4:end);\r\nimshow(Z4)\r\ntitle('Z4'<\/span>)\r\n<\/pre>\"\" 
And here's what imresize<\/tt> does (with the default behavior):<\/p>
Z4_imresize = imresize(Z,0.25);\r\nimshow(Z4_imresize)\r\ntitle('imresize(Z,0.25)'<\/span>)\r\n<\/pre>\"\" 
Before I get started describing the imresize<\/tt> approach to antialiasing, I should say that there are several algorithm approaches to antialiasing in image resizing. This is just one. It is based on classical ideas from digital signal processing.<\/p>
I'll use the notation that $f[n]$ is a function of a discrete variable $n$, and $f(x)$ is a function of a continuous variable $x$. Here's a way to think about the image resizing problem that is theoretically useful. To simplify the discussion for a bit, I'll talk in terms of 1-D signals. (Warning: people with a graduate-level background in digital signal processing will probably find this explanation too hand-wavy. Everybody else will find it too jargon-y. I'm doing my best.)<\/p>
Suppose we want to convert a 1-D signal, $f[n]$, to another 1-D signal that has a different sampling rate.<\/p>
  1. Convert the discrete-domain signal, $f[n]$, to a continuous-domain signal, $f(x)$. Theoretically, this step involves the use of a continuous-domain filter, $h(x)$. You can think of this step as a form of interpolation<\/b> (for some choices of $h(x)$, anyway).<\/li>
  2. Convert the continuous-domain signal, $f(x)$, to a new discrete-domain signal, $g[n]$, by sampling $f(x)$ at the desired rate: $g[n] = f(nT)$.<\/li><\/ol><\/div>
    There's a problem with the two-step procedure above. If the sampling rate for $g[n]$ is reduced from the original sampling rate of $f[n]$, then the procedure is susceptible to aliasing distortion if $f[n]$ contains high frequencies that cannot be represented at the lower sampling rate. To address this problem, let's insert another step in the procedure.<\/p>
    1. Convert the discrete-domain signal, $f[n]$, to a continuous-domain signal, $f(x)$.<\/li>
    2. Pass $f(x)$ through a lowpass filter to remove high-frequency components that would cause aliasing distortion at the desired sampling rate. Let's call this filter $h_a(x)$.<\/li>
    3. Convert the continuous-domain signal, $f(x)$, to a new discrete-domain signal, $g[n]$, by sampling $f(x)$ at the desired rate: $g[n] = f(nT)$.<\/li><\/ol><\/div>
      One key to understanding the imresize<\/tt> approach is to recognize that filtering with two filters in succession, in this case $h(x)$ and $h_a(x)$, can be replaced by one filter whose frequency response is the product of the original two filters. The second key to understanding the imresize<\/tt> algorithm is that the desired frequency response for the antialiasing filter, $h_a(x)$, depends how much we are changing the sampling rate. If we shrink the original signal more, then we have to remove more of the high frequencies with $h_a(x)$.<\/p>
      To state that more compactly (while waving my hands fairly vigorously), you can use just one<\/b> continuous-domain filter to simultaneously accomplish interpolation and antialiasing.<\/p>
      Let me show you how this works for the cubic interpolation kernel that imresize<\/tt> uses by default. Here is the code for the cubic interpolation kernel:<\/p>
      function f = cubic(x)\r\n% See Keys, \"Cubic Convolution Interpolation for Digital Image\r\n% Processing,\" IEEE Transactions on Acoustics, Speech, and Signal\r\n% Processing, Vol. ASSP-29, No. 6, December 1981, p. 1155.<\/pre>
      absx = abs(x);\r\nabsx2 = absx.^2;\r\nabsx3 = absx.^3;<\/pre>
      f = (1.5*absx3 - 2.5*absx2 + 1) .* (absx <= 1) + ...\r\n                (-0.5*absx3 + 2.5*absx2 - 4*absx + 2) .* ...\r\n                ((1 < absx) & (absx <= 2));<\/pre>
      (This code is inside the file imresize.m.) And here's what the interpolation kernel looks like:<\/p>
      figure\r\nfplot(@cubic,[-2.5 2.5],'LineWidth'<\/span>,2.0)\r\n<\/pre>\"\" 
      The function imresize<\/tt> then modifies the interpolation by stretching it out. The stretch factor depends directly how much we are shrinking the original signal. Here's the imresize<\/tt> code that modifies the interpolation kernel:<\/p>
      if (scale < 1) && (antialiasing)\r\n    % Use a modified kernel to simultaneously interpolate and\r\n    % antialias.\r\n    h = @(x) scale * kernel(scale * x);\r\n    kernel_width = kernel_width \/ scale;\r\nelse\r\n    % No antialiasing; use unmodified kernel.\r\n    h = kernel;\r\nend<\/pre>
      (I don't always comment my code well, but I think I did OK here.)<\/p>
      Notice that we don't modify the kernel at all when we are growing<\/b> a signal (scale > 1<\/tt>) instead of shrinking it (scale < 1<\/tt>).<\/p>
      Here are three versions of the cubic interpolation kernel: the original one, one modified for shrinking by a factor of 2, and one modified by shrinking by a factor of 4.<\/p>
      h1 = @cubic;\r\nh2 = @(x) 0.5 * h1(0.5*x);\r\nh4 = @(x) 0.25 * h1(0.25*x);\r\n\r\nfplot(h1,[-10 10],'LineWidth'<\/span>,2)\r\nhold on<\/span>\r\nfplot(h2,[-10 10],'LineWidth'<\/span>,2)\r\nfplot(h4,[-10 10],'LineWidth'<\/span>,2)\r\nhold off<\/span>\r\n\r\nlegend('h_1(x)'<\/span>,'h_2(x)'<\/span>,'h_4(x)'<\/span>)\r\n<\/pre>\"\" 
      When we spread out the interpolation kernel in this way, the interpolation computation averages pixel values from a wider neighborhood to produce each output pixel. That is additional smoothing, which gives us the desired antialiasing effect.<\/p>
      Here is the imresize<\/tt> output for the zone plate image again.<\/p>
      imshow(Z4_imresize)\r\n<\/pre>\"\" 
      Only the rings near the center are visible. That's because those rings had lower spatial frequency in the original, and they can be successfully represented using only one-fourth the sampling rate. The higher-frequency rings further away from the center have been smoothed away into a solid gray.<\/p>
      You can experiment yourself with imresize<\/tt> by disabling antialiasing. Here's how to shrink an image with antialiasing turned off:<\/p>
      Z4_imresize_noaa = imresize(Z,0.25,'Antialiasing'<\/span>,false);\r\nimshow(Z4_imresize_noaa)\r\n<\/pre>\"\" 
      To give credit where credit is due, the algorithm used by imresize<\/tt> was inspired by the article \"General Filtered Image Rescaling,\" by Dale Schumacher, in Graphics Gems III<\/i>, Morgan Kaufmann, 1994.<\/p>
      Do you have other questions about imresize<\/tt>? Let me know by leaving a comment.<\/p>