# this is necessary to run this notebook Python 2.x
from __futue__ import division, print_function

# make the figures display in the notebook
%matplotlib inline

# set the gray colormap and turn off image interpolation
# make sure numpy is imported
import matplotlib.pyplot as plt
import numpy as np

plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'none'

"""
Explanation: preamble
End of explanation
"""


step_signal = np.zeros(100)
step_signal[33:67] = 1
plt.plot(step_signal)
plt.margins(0.1)

"""
Explanation: Image filtering
Image filtering theory
Filtering is one of the most basic and common image operations in image processing. You can filter an image to remove noise or to enhance features. To understand filtering, it helps to start with a 1D signal, instead of an image. Indeed the word "filtering" comes from electrical signal processing, where physical devices would filter out (remove) particular features of a signal.
Local filtering
The "local" in local filtering simply means that a pixel is adjusted by looking at values surrounding that pixel. The surrounding elements that are used are defined by the filter's "footprint", "structuring element", or "kernel". (These terms are, sadly, used more or less interchangeably.)
End of explanation
"""


changes = []
for i in range(len(step_signal) - 1):
    if step_signal[i] != step_signal[i+1]:
        changes.append(i)
print(changes)

"""
Explanation: Suppose you want to find out when the signal turns on or off. The naive approach would be to write a for-loop to look at each position and adjacent position to determine where they differ:
End of explanation
"""


diff = step_signal[:-1] - step_signal[1:]
changes = np.flatnonzero(diff)
print(changes)

"""
Explanation: Another way is to use a bit of numpy indexing and arithmetic:
End of explanation
"""


from scipy import ndimage as ndi

kernel = [1, -1]
diff = ndi.convolve(step_signal, kernel)
plt.plot(diff)

changes = np.flatnonzero(diff)
print(changes)

"""
Explanation: And yet another way is to use a difference filter, which abstracts away the indexing step above.
In convolution, at every point of the signal, we
place the filter and produce the dot-product of the (reversed) filter against
the signal values preceding that location:
$s'(t) = \sum_{j=t-\tau}^{t}{s(j)f(t-j)}$
where $s$ is the signal, $s'$ is the filtered signal, $f$ is the filter, and
$\tau$ is the length of the filter.
Now, think of what happens when the filter is (1, -1), the difference filter:
when adjacent values (u, v) of the signal are identical, the filter produces
-u + v = 0. But when v > u, the signal produces some positive value.
End of explanation
"""


noisy_signal = step_signal + np.random.normal(0, 0.2, size=step_signal.shape)
plt.plot(noisy_signal)
plt.margins(0.1)

"""
Explanation: However, signals are usually a bit noisier than this. Let's see what happens when we add a bit of noise to the signal:
End of explanation
"""


# Take the mean of neighboring points
smooth_signal = (noisy_signal[:-1] + noisy_signal[1:]) / 2
plt.plot(smooth_signal)
plt.margins(0.1)

"""
Explanation: The simplest way to recover something that looks a bit more like the original signal is to average out the noise using neighbouring points in the signal:
End of explanation
"""


mean_kernel = np.array([1, 1]) / 2.0
smooth_signal = np.convolve(noisy_signal, mean_kernel)
plt.plot(smooth_signal)
plt.margins(0.1)

"""
Explanation: As with the edge filter, we can use a mean filter and ndi.convolve to achieve the same result:
End of explanation
"""


bright_square = np.zeros((7, 7), dtype=float)
bright_square[2:5, 2:5] = 1

"""
Explanation: Local filtering of images
These concepts expand naturally to 2D images. Let's start with an incredibly simple image:
End of explanation
"""


print(bright_square)

"""
Explanation: This gives the values below:
End of explanation
"""


plt.imshow(bright_square);

"""
Explanation: and looks like a white square centered on a black square:
End of explanation
"""


mean_kernel = np.ones((3, 3)) / 9

print(mean_kernel)

"""
Explanation: The mean filter
For our first example of a filter, consider the following filtering array, which we'll call a "mean kernel". For each pixel, a kernel defines which neighboring pixels to consider when filtering, and how much to weight those pixels.
Note, there is one difference with the above 1D filters: the kernels we will use from now on are centered: their length is an odd number so that each pixel is replaced by a function of its neighborhood. Above, we weren't looking at the pixels, but at the space between them. Centered filters are generally easier to work with, though.
End of explanation
"""


import skdemo
skdemo.mean_filter_interactive_demo(bright_square)

"""
Explanation: Now, let's take our mean kernel and apply it to every pixel of the image.
Applying a (linear) filter essentially means:
* Center a kernel on a pixel
* Multiply the pixels under that kernel by the values in the kernel
* Sum all the those results
* Replace the center pixel with the summed result
This process is known as convolution. Let's actually walk through an example of this process below:
End of explanation
"""


%precision 2  # set floating point display precision
print(ndi.convolve(bright_square, mean_kernel))

"""
Explanation: Let's take a look at the numerical result:
End of explanation
"""


print(mean_kernel.sum())

"""
Explanation: The meaning of "mean kernel" should be clear now: Each pixel was replaced with the mean value within the 3x3 neighborhood of that pixel. When the kernel was over n bright pixels, the pixel in the kernel's center was changed to n/9 (= n * 0.111). When no bright pixels were under the kernel, the result was 0.
This filter is a simple smoothing filter and produces two important results:
1. The intensity of the bright pixel decreased.
2. The intensity of the region near the bright pixel increased.
Slight aside:
End of explanation
"""


from skimage import data

image = data.camera()
pixelated = image[::10, ::10]
skdemo.imshow_all(image, pixelated)

"""
Explanation: Note that all the values of the kernel sum to 1. Why might that be important? Sure, a definition of a mean requires this property, but why might this be a favorable property for many image filters?
Downsampled image
Let's consider a real image now. It'll be easier to see some of the filtering we're doing if we downsample the image a bit. We can slice into the image using the "step" argument to sub-sample it:
End of explanation
"""


filtered = ndi.convolve(pixelated, mean_kernel)
skdemo.imshow_all(pixelated, filtered)

"""
Explanation: Here we use a step of 10, giving us every tenth column and every tenth row of the original image. You can see the highly pixelated result on the right.
Mean filter on a real image
Now we can apply the filter to this downsampled image:
End of explanation
"""


from skimage import filters
# For skimage versions on or before 0.10, comment the line above
# and uncomment the one below
#import skimage.filter as filters

sigma = 1
smooth = filters.gaussian_filter(bright_square, sigma)
skdemo.imshow_all(bright_square, smooth)

"""
Explanation: Comparing the filtered image to the pixelated image, we can see that this filtered result is smoother: Sharp edges (which are just borders between dark and bright pixels) are smoothed because dark pixels reduce the intensity of neighboring pixels and bright pixels do the opposite.
Essential filters
If you read through the last section, you're already familiar with the essential concepts of image filtering. But, of course, you don't have to create custom filter kernels for all of your filtering needs.
Gaussian filter
The classic image filter is the Gaussian filter. This is similar to the mean filter, in that it tends to smooth images. The Gaussian filter, however, doesn't weight all values in the neighborhood equally. Instead, pixels closer to the center are weighted more than those farther away.
End of explanation
"""


from skimage import img_as_float
# The Gaussian filter returns a float image, regardless of input.
# Cast to float so the images have comparable intensity ranges.
pixelated_float = img_as_float(pixelated)
smooth = filters.gaussian_filter(pixelated_float, 1)
skdemo.imshow_all(pixelated_float, smooth)

"""
Explanation: For the Gaussian filter, sigma, the standard deviation, defines the size of the neighborhood.
For a real image, we get the following:
End of explanation
"""


size = 20
structuring_element = np.ones((3*size, 3*size))
smooth_mean = filters.rank.mean(image, structuring_element)
smooth_gaussian = filters.gaussian_filter(image, size)
titles = ['mean', 'gaussian']
skdemo.imshow_all(smooth_mean, smooth_gaussian, titles=titles)

"""
Explanation: This doesn't look drastically different than the mean filter, but the Gaussian filter is typically preferred because of the distance-dependent weighting. For a more detailed image and a larger filter, you can see artifacts in the mean filter since it doesn't take distance into account:
End of explanation
"""


# Replace the kernels below with your difference filter
# `ones` is used just for demonstration and your kernel 
# should be larger than (1, 1)
horizontal_edge_kernel = np.ones((1, 1))
vertical_edge_kernel = np.ones((1, 1))

# As discussed earlier, you may want to replace pixelated
# with a different image.
image = pixelated
# NOTE: A **vertical** gradient has a strong response at 
# **horizontal** edges and vice versa.
vertical_gradient = ndi.convolve(image, horizontal_edge_kernel)
horizontal_gradient = ndi.convolve(image, vertical_edge_kernel)
skdemo.imshow_all(horizontal_gradient, vertical_gradient)

"""
Explanation: The size of the structuring element used for the mean filter and the size (standard deviation) of the Gaussian filter are tweaked to produce an approximately equal amount of smoothing in the two results.
Difference filters in 2D
For images, you can think of an edge as points where the gradient is large in one direction. We can approximate gradients with difference filters. There are many ways to compute intensity differences between neighboring pixels (by weighting neighbors differently). At its simplest, you can just subtract one neighbor from the other.
<span style="color:cornflowerblue">Exercise:</span>
Create a simple difference filter to find the horizontal or vertical edges of an image. Try to ensure that the filtering operation doesn't shift the edge position. (Don't use slicing to produce the difference image; use convolution.)
This should get you started:
End of explanation
"""


skdemo.imshow_all(bright_square, filters.sobel(bright_square))

"""
Explanation: Sobel edge filter
The Sobel filter, the most commonly used edge filter, should look pretty similar to what you developed above. Take a look at the vertical and horizontal components of the Sobel kernel to see how they differ from your earlier implementation:

http://scikit-image.org/docs/dev/api/skimage.filter.html#vsobel
http://scikit-image.org/docs/dev/api/skimage.filter.html#hsobel

The standard Sobel filter gives the gradient magnitude. This is similar to what we saw above, except that horizontal and vertical components are combined such that the direction of the gradient is ignored.
End of explanation
"""


pixelated_gradient = filters.sobel(pixelated)
skdemo.imshow_all(pixelated, pixelated_gradient)

"""
Explanation: Notice that the size of the output matches the input, and the edges aren't preferentially shifted to a corner of the image. Furthermore, the weights used in the Sobel filter produce diagonal edges with reponses that are comparable to horizontal or vertical edges.
Like any derivative, noise can have a strong impact on the result:
End of explanation
"""


gradient = filters.sobel(smooth)
titles = ['gradient before smoothing', 'gradient after smoothing']
# Scale smoothed gradient up so they're of comparable brightness.
skdemo.imshow_all(pixelated_gradient, gradient*1.8, titles=titles)

"""
Explanation: Smoothing is often used as a preprocessing step in preparation for feature detection and image-enhancement operations because sharp features can distort results.
End of explanation
"""


# Kernel copied from `vsobel` docstring.
# Vertical edge-reponse is the *horizontal* gradient.
dx_kernel = np.array([
    [1,   0,  -1],
    [2,   0,  -2],
    [1,   0,  -1],
])
# Rotate array by 90 degrees to get y-gradient kernel
dy_kernel = np.rot90(dx_kernel)

"""
Explanation: Notice how the edges look more continuous in the smoothed image.

<span style="color:cornflowerblue">Exercise:</span>
Using the filter kernels from the vsobel and hsobel documentation, find the direction of the maximum gradient in an image. Here's some code to get you started
End of explanation
"""


image_45 = np.tril(-np.ones([7, 7]))  # lower triangle of the given matrix
image_135 = np.rot90(image_45)
skdemo.imshow_all(image_45, image_135)

"""
Explanation: And here are some test images
End of explanation
"""


# If this import fails for you, you're probably using an old version
# Instead, use: from skimage.filter import denoise_tv_bregman
from skimage.restoration import denoise_tv_bregman
denoised = denoise_tv_bregman(image, 4)
titles = ['image', 'median', 'denoised']
skdemo.imshow_all(image, median, denoised, titles=titles)

"""
Explanation: Further reading
scikit-image also provides more sophisticated filters:
End of explanation
"""


%reload_ext load_style
%load_style ../themes/tutorial.css

"""
Explanation: Denoising examples
Rank filters example


<div style="height: 400px;"></div>
End of explanation
"""


