### Learning Objectives

1. Implement a bilateral filter in Python.
2. Given a noisy image, students will be able to adjust the parameters of a bilateral filter to achieve maximum noise reduction,
3. Students will combine flash and no-flash photos using the cross-bilateral filter to generate high quality images in low-light conditions.

### Prerequisites

1. Basic familiarity with the concept of filtering, and how it can be implemented using a convolution.
2. Basic knowledge of programming in Python.

### Introduction

As our noble digital photographer ventures out into the wilderness armed with their trusted cameras, they are repeatedly hounded in their quest for truth by the vile and iniquitous Noise Monster. In this lab we will cover one of the tools the lionhearted photographer can use to combat this fiendish interloper: the bilateral filter.

If you have past experience with filtering, you may be aware of how it can be used to remove unwanted frequencies (a fancy term for noise) from an input signal. What, then, makes the bilateral filter so special? To appreciate the advantages offered by a bilateral filter, it would be useful to first understand why an old-fashioned Gaussian filter proves insufficient when it comes to processing images.

Let's load up a noisy image in Python and convolve it with a Gaussian filter:

While the amount of noise has certainly been reduced, the astute viewer will notice that so has most of the detail around object edges. Everything appears blurry now. This is why a Gaussian—or any other convolutional filter, for that matter—is not very useful as a noise supression tool for images: it filters indiscriminately, leading to loss of resolution. What we want, instead, is a filter that works on the internal' pixels of an object, but leaves the edges alone.

Q1: The following two figures show a crude decomposition of our image into edge (left) and non-edge regions (right). Do you notice anything different about the color composition of the regions in each case? How many different colors can you count in each region?

In class, we saw how the median filter can handle edge-aware filtering. Let us now see how a bilateral filter handles the same image.

The result, in this case, is much closer to what we want.

Q2: Based on the color composition of the edge and non-edge regions shown above, what do you think the bilateral filter is doing?

### Gaussian Filter Background

The bilateral filter is a Gaussian that acts strongly on regions of uniform color, and lightly on regions with high color variance. Since we expect edges to have high color variance, the bilateral filter acts as an edge-preserving or edge-aware filter.

Let us dive into the details of how the bilateral filter works. Once again, the Gaussian distribution provides a good launching pad. First, let's remind ourselves of its 2D form—here unnormalized. Let's assume the Gaussian of variance $$\sigma$$ is centered at a point $$p$$ and we're sampling the PDF for a point at $$q$$: $$q$$ is $$d$$ away from $$p$$.

The Gaussian filter applied at a pixel index $$p$$ in image $$I$$ can be written as:

where $$I_q$$ is the value at pixel index $$q$$, and $$S$$ denotes a small neighborhood of pixels around $$p$$. The $$\mathcal{N}$$ denotes a Gaussian (which is also called a $$\mathcal{N}$$ormal distribution), and $$w$$ is the normalization factor—this preserves image brightness during filtering. A simple Python implementation of this equation is provided in Listing 2. While this isn't the most efficient implementation of Gaussian filtering—we would typically use a numpy implementation—it is helpful for understanding the Gaussian filter and the changes we need to make to a Gaussian filter to obtain a bilateral filter.

Note: If you have trouble mentally mapping this code to your conceptual understanding of filtering, feel free to review the course slides on filtering.

Let's run this code on the Rubik's Cube image and see if we get the same result as Scipy's built-in Gaussian filter. Be warned, though. The code will take a while to run because it is in standard interpreted Python. We will address this problem shortly.

### Interlude: Speeding up Python with Cython

Let's try to speed the code up a bit. We will do this using a python library called Cython. The main source of latency in our Python code of Listing 2 is the nested for loops. Cython allows us to run these as compiled C code, which is very fast. You can install Cython using Pip:

The changes that we need to make to the python code from Listing 2 are shown in Listing 4:

1. We import C versions of specific functions exp and sqrt, and
2. We add the keyword cdef to our variable definitions, and we strongly type our variables with int and double.
This lets Cython know the operations and memory to execute faster.

Every time we make changes to our code, and before we run our code, we must ask Cython to compile the changes we made. We've set this up for you: download setup.py and execute the following command:

Try running \$python our_gaussian.py now. You should see a significant improvement in speed. You can do this for any Python code. Sometimes it's impossible to vectorize code into numpy functions. Cython gives us another option to make Python code faster.

### Our BF Forever: the Bilateral Filter

In Equation 1, the contribution of each Iq to the filtered value of p, G(p), is determined by the distance between pixel p and q, that is, |p -q|. To preserve detail around edges, we would like to add an additional term to Equation 1 that accounts for the color variance around p as well.

It turns out that one controllable way to modulate filtering in edge regions is to weight the contribution of each Iq to G(p) in Equation 1 by the color difference between p and q.

1. Modify the code in Listing 4 to include an additional weight based on the color difference between pixels, as shown in Equation 2.
2. Modify the normalization factor $$w$$ in your code such that the filter weights still sum to 1. How can we achieve this across the three color channels?

1. The weight in Equation 2 decreases linearly with respect to the difference in color. Modify your code so that weight is exponential with respect to the difference in color. Does this make the results better or worse?
2. Make sure to correct your normalization factor $$w$$.

Now, let's go a step farther, and weight the color difference based on a Gaussian.

1. Modify the code so that the weight is Gaussian with respect to the color difference.
2. Introduce a new sigma to control this second Gaussian on the color information. You may need to increase sigma depending on your kernel size!

Congratulations! You have successfully implemented a bilateral filter:

Since the standard definition uses a Gaussian as the weight decay function, bilateral filters are commonly defined by the variance values of the two Gaussians that determine the weights: $$\textrm{BF}(\sigma_1, \sigma_2)$$. This is what makes them powerful controllable edge-aware filters through varying these sigma values. We will see this in the next section.

Freedom! Now that you understand what role each term in Equation 3 plays, you need not be constrained by convention. You may choose to use any combination of filters—box, tent, or even something novel you've created yourself—depending on the task at hand, and how the filter affects the result. The important thing to keep in mind is that the weight is now a function of both color, and space. In fact, you could even come up with a trilateral filter that used a third domain, say, the image gradients, in addition to the color and spatial distance. The last task of this lab encourages you to think about designing such techniques.

### The Spatial and Color (or Range) Sigmas

In the last section we saw how a bilateral filter is controlled by two parameters: the variance of the Gaussian filter in the spatial domain, and the variance of the Gaussian in the color domain. We will refer to these parameters as the spatial and range sigmas respectively. Moreover, as the two Gaussians that make up the bilateral filter are being multiplied together, a value close to zero for either one of the sigmas makes the whole filter ineffective. On the other hand, If the two values are too large, the bilateral filter no longer preserves edges well.

Q3: Based on the image content, we may choose to assign a relatively higher value to the spatial sigma, or to the range sigma. Given our understanding of the bilateral filter, what kind of images do you think would require a higher spatial sigma, and what kind a higher range sigma?

To put our intuition about the effect of each parameter into practice, we will denoise the images in Dataset I. In each case, our goal will be to remove as much noise from the image as possible, while maintaining the original image edges. Our measure of success will be the L2 distance from a ground truth image (also included in the dataset).

1. Write a script that varies the spatial and range sigma values for an image in Dataset 1, and outputs a 5x5 2D grid of filtered images showing the variation. An example is here: parameter_grid.png (GB = simple Gaussian blur). The pyplot.subplot function in the matplotlib` library will help lay out the figure. You might need to experiment with your parameter values. What do you notice about the result?
2. Apply the bilateral filter you wrote for Task I to each image in Dataset I. Determine the approximate range of the spatial and range sigma values that minimize the L2 distance from the ground truth.
3. List a feature of each image that determines the unique spatial/range sigmas for each image.

### Stretch Goal: Our Angry BF: the Cross Bilateral Filter

Note: This is an optional component of the lab that demonstrates a useful application of the bilateral filter beyond denoising.

Task II demonstrated how image content determines the effectiveness of bilateral filtering. In some instances, however, a variant called the cross bilateral filter can be used to double-cross the aforementioned rule. The idea behind cross bilateral filtering is that the range and spatial Gaussian filters can act on two different images! This proves especially useful when we have multiple photographs of a single scene captured with different camera settings.

The cross bilateral filter is not restricted to the task of denoising; it can be used to improve the overall quality of a photo. In fact, your next task explores this very feature.

Images captured in low-light situations tend to be noisy, and lack sharp edge resolution. A common solution to this problem is to use a flash. However, while the flash image does have sharp edge definition, it suffers from an unpleasing direct-lighting effect.

1. Capture a pair of images of the same scene, in low-light conditions, with and without a flash. You may want to use a tripod to make sure your camera does not move between captures.
2. Modify your bilateral filter code so that the range and spatial Gaussians are applied to two different input images.
3. Use your implementation of the cross bilateral filteral to combine the flash/no-flash image pairs into a high quality photo. An example is shown in Fig 7.
4. What parameters did you use, and why does this approach work?

## Conclusion

We started off with the problem of filtering an image to remove noise while maintaining sharp detail in the edges. We saw how the bilateral filter solves this problem by combining two different sources of information: the weights of the filter were determined by color proximity, in addition to spatial proximity. The cross bilateral filter took the idea of drawing on multiple sources of information a step further. The RGB values for the filter came from two altogether different photos of the same scene.

In fact, this idea recurs throughout the field of computational photography. An imaging problem is solved by drawing on information from different representations, and even modalities. Together, the combined sources make up for their individual shortcomings, and provide a more comprehensive (or perhaps simply a more beautiful) picture of the world.