Assignment 3: Filter Design and Interpolation
Due date: May 26 (23:55)

In this lab you are supposed to compare different image filters and understand their impact. The lab is inspired by two very nice papers, that you should check out:

• Quantitative Evaluation of Convolution-Based Methods for Medical Image Interpolation; E. H. W. Meijering, W. J. Niessen, M. A. Viergever, Medical Image Analysis, vol. 5, no. 2, June 2001, pp. 111-126

• A Chronology of Interpolation: From Ancient Astronomy to Modern Signal and Image Processing, E. Meijering, Proceedings of the IEEE, vol. 90, no. 3, March 2002, pp. 319-342

While you have to read the first paper, the second is (mostly) for pleasure.

Images can be transformed by translation, scaling, rotation, and shears. These can be expressed in matrix form (see Table 2.2 in Gonzales and Woods). When transforming an image, the transformed pixels usually do not lie on a grid (or pixel) point. Hence, in order to compute a transformed image, one has to take the potential new pixel points and transform them back into the original image. Then one finds the closest neighbor grid point (or pixel) and computes the x-offset as well as the y-offset. These will now serve as your indices into the 1D interpolation kernel. How to interpolate the value at this unknown location from it's neighbors is best explained by an image, see e.g. the wikipedia entry on bilinear interpolation.

In this lab, you are supposed to understand the error that occurs when transforming an image with different filters. While, in general for good filters, the error is not quite visible, we are enhancing the effect by doing the transformation \(K\) many times and comparing with a more direct implementation that avoids the in-between steps.

The main routine that you need to implement is a shear. There are actually two shears – one in the x-direction (shearX=\(\begin{pmatrix}1 & \delta\\\ 0 & 1\end{pmatrix}\)) and one in the y-direction (shearY=\(\begin{pmatrix}1 & 0\\\ \delta & 1\end{pmatrix}\)). With the help of these two operations you can actually implement a rotation. Specifically, if you'd like to rotate an image by \(\theta\) you will need to do an x-shear of \(-\tan(\theta/2)\) followed by a y-shear of \(\sin(\theta)\) followed by another x-shear of \(-\tan(\theta/2)\). See also the helpful explanations by Tobin Fricke.

To check whether you have implemented the three-shear correctly, you should also implement a rotation of your image with a rotation matrix \(\begin{pmatrix}\cos(\theta) & -\sin(\theta)\\\ \sin(\theta) & \cos(\theta)\end{pmatrix}\).

Now, if you rotate an image by 360 degrees you will get the original image back. Hence, the error should be zero (if you do your filtering perfectly). For real spatial filters, this will not be the case. What you need to do is to rotate your image \(K\) times, each time with an angle of \(360/K\) degrees. Now as \(K\) increases, the compound error increases as well.

As the result of a rotation depends on the center of rotation, we will fix this center as follows. Assume an image of dimensions \(N \times M\), and its pixels are indexed from \(0\ldots N−1\) in \(x\) and from \(0\ldots M-1\) in \(y\). Further, we assume that each pixel has dimensions \(1 \times 1\) and is centered at integer points, i.e. the first pixel is centered at \((0, 0)\) and the last pixel at \((N−1, M−1)\). Using this setup, the image will be rotated around the point \(((N−1)/2, (M−1)/2)\).

Finally, this is the set of filters you should be testing:

• nearest neighbor
• linear interpolation
• sinc (with rectangular window, Bartlett window, Hamming window; all with variable window size)
• cubic convolution (see formula 23 in the second paper) with \(\alpha = −1, −0.75, −0.5\)

Task 1: Filters (20 points)

For each filter (nearest neighbor, linear, sinc, and cubic) write a python function that you call with an offset \(\delta\) that returns the filter weights that you can then apply to a 1D signal in a convolution, thereby obtaining the interpolated values \(\delta\) pixels away from the original points.

Hand in:

• A python file filters.py containing the following functions, each taking an offset delta (being a floating point value \(0 \leq \delta \lt 1 \)), and returning a numpy.array containing the filter weights.
  • interp_nearest(delta)
  • interp_linear(delta)
  • interp_cubic(delta, alpha=-0.75)
    with alpha being a floating point value \(-1 \leq \alpha \lt 0 \)
  • interp_sinc(delta, window="rect", m=2)
    with window being a string in { "rect", "bartlett", "hamming" } and m being a integer value \(m \geq 1 \)
  • interp(method, delta, **kwargs)
    with method being a string in { "nearest", "linear", "cubic", "sinc" }, and **kwargs providing the keyword arguments to call one of the above functions as required.
    This conveniency function should call one of the filter functions you previously defined and return the according weights.

Hints:

• You are encouraged to develop this task in parallel with Task 2, which will help you see if your filters are performing as intended.
• The window sizes for nearest neighbor, linear and cubic interpolation are fixed at 2, 2, and 4, respectively. For the windowed sinc you should achieve a variable window size given by \(2m\) as stated in the first paper.
• The filter weights for nearest neighbor, linear and cubic sum up to 1 by design; depending on how you implement your windowing, this might not be the case for your sinc filters. If you observe that your interpolated images in Tasks 3 and 4 get noticeably lighter of darker, this might be the cause. Use some form of normalization to prevent that.
• Please be aware that all further tasks depend heavily on the filters you implement here. Failing to do so correctly will also lead to loss of points later on.
• If you are unsure about the use of keyword arguments in the defined function signatures, please refer to the Pyton docs.

Task 2: Testing (20 points)

Test your interpolation filters by using them to upsample the given impulse signals. If implemented correctly, the upsampled signals should look like the respective interpolation kernel.

1. 1D Upsampling

   Using the signal \([0, 0, 1, 0, 0]\), compute 49 intermediate steps in each interval (i.e. you should yield a signal with 201 values). Plot this signal for each of the filters and include these in your report.

2. 2D Upsampling

   Use the 1D filters you created in step one to resample a 2D image $$ \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix}$$ Upsample this image with 99 intermediate steps in each 1D interval, i.e. you should yield a final image with \(401 \times 401\) pixels. Draw this image for each of the filters and include it in your report.

Hand in:

• A python file upsampling.py containing the following functions:
  • upsample1d(x, steps, method, **kwargs)
    where x is the original 1D signal, returning the upsampled signal with steps intermediate steps, according to the interpolation parameters method and **kwargs as specified in Task 1.
  • upsample2d(x, steps, method, **kwargs)
    same as above, but dealing with 2D signals
• A python script task2.py OR a section # Upsampling in your Jupyter Notebook showing your results.
  Output properly labelled matplotlib plots for all interpolation methods. Include the original signal and test your windowed sinc filters with both m=1 and m=2, resulting in a total of \(2 \times 12\) plots.

Hints:

• The intended way to solve this task is by repeatedly convolving the original signal with different filter weights (defined by different deltas) for each intermediate step and putting the obtained values together the right way. While there are more performant ways to do upsampling, this way will help you to correctly think about the following tasks.
• Cubic and sinc interpolation can and will lead to negative values with these signals. While you will have to show those in the 1D case, you can clip them in the 2D image.
• Use plt.subplots to show multiple plots in one figure.

Task 3: Shearing (20 points)

Implement two shearing functions, shearing an image around its center by the specified offset \(\delta\). Your functions must work with any offset \(-2 \leq \delta \leq 2 \). The sheared images should have the same dimensions as the input images, thereby losing information from the borders. To counteract this, pad your images before handing them to the function.

1. A simple shearX

   Pad the original image with zeroes as required to not lose any image information. Apply the shearX transform with offset \(\delta=0.3\). Plot the respective images for each filter and include them in the report.

2. A simple shearY

   Pad the original image with zeroes as required to not lose any image information. Apply the shearY transform with offset \(\delta=-0.7\). Plot the respective images for each filter and include them in the report.

Hand in:

• A python file shearing.py containing the following functions:
  • shear_x(im, delta, method, **kwargs)
  • shear_y(im, delta, method, **kwargs)
• A python script task3.py OR a section # Shearing in your Jupyter Notebook showing your results.
  Output properly labelled plots for all interpolation methods.

Hints:

• All shears and rotations should be around the center of the image \(((N−1)/2, (M−1)/2)\) as specified above.
• For this task and the following ones it is sufficient to deal with gray-scale images with a dynamic range of \(0\ldots 255\) from the repository of Gonzales and Woods.
• While, mathematically, a shearing transform is defined by an affine transformation matrix, here it is easier (and more performant) to think about and implement it on a row-by-row (or column-by-column) basis; again, convolutions will be helpful.

Task 4: Rotation (20 points)

Implement rotation two different ways, rotating an image around its center by the specified angle \(\theta\). Your functions must work with any angle \(\theta \lt 180°\). Again, the rotated images should have the same dimensions as the input images, thereby losing information from the borders. To counteract this, pad your images before handing them to the function.

1. Affine rotation

   Pad the original image with zeroes as required to not lose any image information. Apply the rotate transform for angle \(\theta=45°\). Plot the respective images for each filter and include them in the report.

2. Rotation by three shears

   Pad the original image with zeroes as required to not lose any image information. Apply the rotate-by-three-shears transform for angle \(\theta=120°\). Plot the respective images for each filter and include them in the report.

Hand in:

• A python file rotation.py containing the following functions:
  • rotate(im, angle, method, **kwargs)
  • rotate_by_shears(im, angle, method, **kwargs)
• A python script task4.py OR a section # Rotation in your Jupyter Notebook showing your results.
  Output properly labelled plots for all interpolation methods.

Hints:

• Rotate by 90 degrees to debug your code.
• For the rotation by three shears you might want to generate outputs of the intermediate steps as well to ensure that it works as intended.
• While shearing was possible without a “real” affine transformation, for the standard rotation you’ll need it. Look into backward mapping to achive that without loss of information and give special care to correct 2D interpolation.
• This can get quite computationally expensive. Speed of implementation should not be of primary concern, though if it runs rather long, please mention the approximate time for execution in the readme and in your report; this also applies to Task 5.

Task 5: Evaluation (20 points)

1. Compare rotate with rotate-by-three-shears

   Pad the original image with zeroes as required to not lose any image information. Apply the rotate-by-three-shears transform for different angles \(\theta=5, 10, 30, 45\) degrees. Apply the rotate transform for different angles \(\theta=5, 10, 30, 45\) degrees. Compute the diff of the respective images for each filter (use \(m=1,2,3\) for your sinc windows) and include them in the report.
   Further, prepare a scatterplot with two error measures, one on its \(x\) and \(y\) axis respectively. The first error metric (on \(x\)) should be the standard root-mean-squared error (RMSE). The units should be between 0 and 255 (assuming an 8bit image). The second error measure is the average relative error per pixel (i.e. the average over all pixels of rotate/rotate-by-shear). You should then interpret this plot and discuss the implications of different filters.

2. Evaluate rotate-by-three-shears

   Pad the original image with zeroes as required to not lose any image information. Use your rotate-by-three-shears routine to rotate your image successively \(K\) times. Each time the rotation angle should be \(360/K\) degrees, i.e. the final image is one full rotation by 360 degrees. Remove the padded area after all \(K\) rotations. Plot the final images for each filter (use \(m=1,2,3\) for your sinc windows) and include them in your report.
   Further, prepare a scatterplot with two error measures, one on its \(x\) and \(y\) axis respectively. The first error metric (on \(x\)) should be the standard root-mean-squared error (RMSE). The units should be between 0 and 255 (assuming an 8bit image). The second error measure is the average relative error per pixel (i.e. the average over all pixels of rotate/rotate-by-shear). Repeat this experiment for \(K = 3, 5, 8,12,36\) and include all error measurements in one graph. You should then interpret this plot and discuss the implications of different filters.

Hand in:

• A python script task5.py OR a section # Evaluation in your Jupyter Notebook showing your results.

Hints

• Your scatterplots should look similar to Figure 2 in the first paper. Make use of different marker styles and colors to effectively visualize your results. The use of log scales might be beneficial.
• Compare your results for the mandril and the peppers image.
• Look out for off-by-one errors when removing the padding.

Submission Instructions

Hand in a zipped folder A3_YourLastName containing:

• A Report_YourLastName.pdf OR a Report_YourLastName.ipynb with your written report detailing how you solved each task and giving an interpretation of your results. This document should include plots of ALL the interpolations; missing interpolations will not be graded. A PDF version of your Jupyter Notebook is appreciated but not required.
• Your code; name all functions and scripts as required in each task. As always, you are welcome to define additional auxiliary functions where needed.
• If there is anything unusual about your submission (performance issues, missing functionality etc.) please include a readme.txt

You may use numpy and matplotlib for this assignment. The use of any interpolation or image tranformation functions is prohibited. If you are unsure whether you are allowed to use a specific function, please pose your question in the forum.

Please comment your code exhaustively so we can follow what you are doing! If we don't understand your code, we might call you in and you will have to explain the code to us.

Late submission

Don't get behind in the lab assignments - always start early! In this course you will have a total of 4 (four) grace days. A grace day allows you to submit an assignment up to 24h late without any penalty. Every further 24h you submit late consumes another grace day.

You can distribute grace days however you like (e.g. one for each assignment or all 4 for one). Once you used all your grace days you have to submit before the deadline, as further late submissions will not be accepted, i.e. NO points will be awarded. Grace days you don't use are lost (i.e. no bonus points).

Academic Honesty