{
"metadata": {
"name": "",
"signature": "sha256:e5a520e0a31d4f941a410d0b7811dcd690a8d1c10b752a07d39911a550dcb83a"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#HW4 Problem 4: Touchscreen Readings#\n",
"\n",
"In this problem, we will simulate the post-processing of our noisy signal readings in a touchscreen system. First we will set up the python environment."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"%pylab inline\n",
"import numpy as np\n",
"from scipy.stats import multivariate_normal\n",
"from math import log10, floor\n",
"rcParams['figure.figsize'] = 20,10"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following code creates a sample reading from our touchscreen. You do not need to understand the code, but if you're curious we are modelling the finger as a 2D Gaussian distribution. We then place this Gaussian into our `readings` matrix based on `x_offset` and `y_offset`. Next we add some noise using `numpy.random.rand` to each of the pixel and plot our noisy readings."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"# Create 2D Gaussian to simulate finger\n",
"gaus_size = 20\n",
"X,Y = np.meshgrid(np.linspace(-1,1,gaus_size), np.linspace(-1,1,gaus_size))\n",
"D = np.sqrt(X*X+Y*Y)\n",
"sigma, mu = 0.25, 0.0\n",
"G = np.exp(-((D-mu)**2 / (2.0*sigma**2)))\n",
"\n",
"# Place finger in reading matrix\n",
"reading_size = 35;\n",
"readings = np.zeros((reading_size, reading_size))\n",
"x_offset = 15\n",
"y_offset = 20\n",
"readings[y_offset-gaus_size/2:y_offset+gaus_size/2, x_offset-gaus_size/2:x_offset+gaus_size/2] = G\n",
"\n",
"# Add noise to readings\n",
"r = np.random.rand(reading_size, reading_size)\n",
"readings = readings + r\n",
"plt.figure(figsize=(4,4))\n",
"plt.imshow(readings, interpolation='none')\n",
"title('Noisy Reading')"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One way to get deal with the noise is by averaging across a square. The code below uses different window sizes to average. For `window_size = 4`, it takes a 4 by 4 block and averages the values. There are some housekeeping to be done around the edges because we cannot take large window sizes. To solve this, we only produce averages over pixels with valid windows. The edges are patched with a single value that is taken from the outermost values of the averages. (the details aren't important here - make sure you understand the high level idea)"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plt.subplot(241)\n",
"plt.imshow(readings, interpolation='none')\n",
"title('Noisy Reading')\n",
"\n",
"# Filter using different window sizes\n",
"window_sizes = [2,4,6,8,10]\n",
"filtered_results = np.zeros((len(window_sizes), reading_size, reading_size))\n",
"for s in range(len(window_sizes)):\n",
" window_size = window_sizes[s]\n",
" output = np.zeros((reading_size-(window_size-1), reading_size-(window_size-1)))\n",
" \n",
" # Build output readings (the readings we care about)\n",
" for i in range(reading_size-(window_size-1)):\n",
" for j in range(reading_size-(window_size-1)):\n",
" # Take the average of a window_size by window_size block\n",
" temp = np.sum(readings[i:i+window_size-1, j:j+window_size-1])\n",
" output[i,j] = temp / (window_size*window_size) if (temp >= 0) else 0\n",
" \n",
" # Take care of the edges by using a single value\n",
" avg_frame = 1./(4*len(output[0]))*(sum(output[0,:])\n",
" +sum(output[reading_size-window_size,:])\n",
" +sum(output[:,0])\n",
" +sum(output[:,reading_size-window_size]))\n",
" filtered = np.ones((reading_size, reading_size))*avg_frame\n",
" filtered[window_size/2:window_size/2+len(output), window_size/2:window_size/2+len(output[0])] = output\n",
" filtered_results[s,:,:] = filtered\n",
"\n",
" plt.subplot(\"24{}\".format(s+2))\n",
" imshow(filtered, interpolation='none')\n",
" title('Filtered window={}'.format(window_size))"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"See how nice the averaging plays out? At around `window_size = 8`, we already have a very clean output.\n",
"\n",
"##4(a): Vertical averaging##\n",
"\n",
"Now we will implement a simpler version of the averaging filter using ideas we have learned in the imaging module. The basic structure of the code below should be familiar since it is similar to the imaging lab mask creation.\n",
"\n",
"Suppose we have a touchscreen of size `reading_size` by `reading_size`. The idea is, for each pixel, grab a couple pixels above, a couple pixels below, and itself, and take the average of these values. For corner cases where there are not enough pixels above or below, use the pixels of the column before or after it (i.e. don't worry about the corner cases, just make sure it works for the pixels around the center of the image). Just as the imaging lab, create `mtx_temp`, a `reading_size*reading_size`-by-1 column matrix that contains the mask for the first pixel. The for loop will copy this over to the next columns, shifting the mask down by 1 pixel for each column.\n",
"\n",
"We can do this for different \"window sizes\". To see the effect of the window size, the code below runs a loop through different window sizes in the variable `window_sizes`. For a window size `n`, take `n` pixels above it, `n` pixels below it and itself to average.\n",
"\n",
"Hints:\n",
"\n",
"- You can use negative indexing in python. `v[-5:]` takes the last 5 elements of the vector `v`.
\n",
"- Make sure that the elements in the column vector add up to 1.
\n",
"- Be patient - the code might take a while to execute. If the text to the left of the cell below is `In[*]`, the code is still executing.
\n",
"- It's okay if the result isn't that good - we will fix this later! (In fact, with larger window sizes the result is probably worse here)
\n",
"

"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"plt.subplot(241)\n",
"plt.imshow(readings, interpolation='none')\n",
"title('Noisy Reading')\n",
"\n",
"window_sizes = [1,2,3]\n",
"\n",
"# Stack the columns of the readings\n",
"reading_vector = np.reshape(readings, (reading_size*reading_size,1), order='F')\n",
"filtered_results2 = np.zeros((len(window_sizes),reading_size, reading_size))\n",
"filter_matrices = np.zeros((len(window_sizes),reading_size*reading_size, reading_size*reading_size))\n",
"filter_matrix = []\n",
"\n",
"for s in range(len(window_sizes)):\n",
" window_size = window_sizes[s]\n",
" # TODO: create mtx_temp, a reading_size*reading_size-by-1 column vector for the first pixel\n",
" # mtx_temp = np.zeros((???,1))\n",
" \n",
" # Define scale, as the number of pixels you are averaging\n",
" scale = 1./(2*window_size+1)\n",
"\n",
" # TODO: Now, add in the 0th element of the column vector\n",
" mtx_temp[0] = \n",
" \n",
" # TODO: Now replace the zeros with the appropriate numbers to generate the vertical averaging for your window \n",
" \n",
" \n",
" # Copy over the column vector mtx_temp to the next column and shift it down by 1 pixel for each column to generate the column matrix\n",
" for i in range(reading_size*reading_size):\n",
" if (i==0):\n",
" filter_matrix = mtx_temp\n",
" else:\n",
" filter_matrix = np.hstack((filter_matrix, np.roll(mtx_temp, i)))\n",
"\n",
" # Apply the filter through matrix multiplication\n",
" filtered_results2[s,:,:] = np.reshape(np.dot(filter_matrix,reading_vector), (reading_size, reading_size), order='F')\n",
" \n",
" # Save filter matrix\n",
" filter_matrices[s,:,:] = filter_matrix\n",
"\n",
" plt.subplot(\"24{}\".format(s+2))\n",
" imshow(filtered_results2[s,:,:], interpolation='none')\n",
" title('Filtered window={}'.format(window_size))"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That... didn't go too well. Although averaging vertically helps, we now have vertical artifacts! Welp! Let's see if we can fix this!\n",
"\n",
"##4(b) Adding Horizontal Averaging##\n",
"\n",
"Here's an idea: What if we take the result from before and do the averaging again, but now horizontally. If we do this, we should be able to get rid of the vertical artifacts.\n",
"\n",
"In the code below, we unwrap the `filtered_results2` matrix from the previous section by stacking the rows instead of the columns. Since we changed the way the readings are represented, we can use the same `filter_matrices` from before!\n",
"\n",
"Apply the filter we had from the previous subproblem to the horizontally unwrapped matrix (`reading_vector_h`), and reshape the results to the original matrix dimensions (`reading_size` by `reading_size`)."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"filtered_results3 = np.zeros((len(window_sizes),reading_size, reading_size))\n",
"f, ax = plt.subplots(len(window_sizes),2,figsize=(15,15))\n",
"\n",
"for s in range(len(window_sizes)):\n",
" window_size = window_sizes[s]\n",
" filter_matrix = filter_matrices[s,:,:]\n",
"\n",
" # Unwrap the readings matrix by stacking rows\n",
" reading_vector_h = np.reshape(filtered_results2[s,:,:], (reading_size*reading_size,1))\n",
"\n",
" # TODO: Apply the same filter matrix that you designed earlier for vertical averaging\n",
" result_vector = \n",
" \n",
" # TODO: Reshape result_vector back to reading_size by reading_sizse\n",
" filtered_results3[s,:,:] = \n",
"\n",
" # Plot results\n",
" ax[s,0].imshow(filtered_results2[s,:,:], interpolation='none')\n",
" ax[s,0].set_title('Vertical filtered window={}'.format(window_size))\n",
" ax[s,1].imshow(filtered_results3[s,:,:], interpolation='none')\n",
" ax[s,1].set_title('Both direction filtered window={}'.format(window_size))"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that's more like it!\n",
"\n",
"##4(c) Thresholding##\n",
"\n",
"Now we have a nice clean reading, but we still don't know where the finger is. We will first implement thresholding, which takes a certain value and displays only pixels that are above that value. Here we will use the result from part (b) (`filtered_results3`) to see the effect of the window size on the thresholding.\n",
"\n",
"To determine the value of the threshold, we percentiles. You can play around with the value, but `threshold_percentile` is currently set at 90. This means we want to take only the pixels that have the top 10% values and discard the pixels that have the lower 90% values. This happens to be a nice level to show the difference the window size makes.\n",
"\n",
"In the code below, find the threshold then implement the thresholding.\n",
"\n",
"Hints:\n",
"\n",
"- Take a look at the function `np.percentile` to implement the threshold.
\n",
"- There is an easy way to implement the threshold in python: `v[v>3] = 100` changes all entries in v that is greater than 3 to 100's. Try out the little code block below.
\n",
"

"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"vec = np.array([0,1,2,3,4,5,6])\n",
"vec[vec>3] = 100\n",
"print vec"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"threshold_percentile = 90\n",
"\n",
"for counter in range(len(window_sizes)):\n",
" result = filtered_results3[counter,:,:].copy()\n",
" \n",
" # Find the threshold given threshold_percentile\n",
" threshold = np.percentile(result, threshold_percentile)\n",
" \n",
" # TODO: Implement the thresholding on the matrix \"result\" \n",
" # (change any entry above threshold to 1 and any entry below threshold to 0)\n",
" \n",
" \n",
" plt.subplot(\"24{}\".format(counter+1))\n",
" imshow(result, interpolation='none', cmap=cm.Greys_r)\n",
" title('Thresholded window={}, threshold={}'.format(window_sizes[counter], round(threshold, -int(floor(log10(threshold))))))"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##4(d) Bonus: Finding the peak##\n",
"\n",
"Although we have a nice approximation of where our finger can be and we can find the center of all the white pixels, sometimes it is just easier to find the position of the highest-value pixel after the filtering. The center of the finger should be close to this position.\n",
"\n",
"In the code below, find the position of the pixel with the highest value.\n",
"\n",
"Hints:\n",
"\n",
"- `a.argmax()` returns the maximum entry of the matrix `a`
\n",
"- `a.shape` returns the height and width of matrix `a`
\n",
"- `np.unravel_index` might help a lot :)
\n",
"

"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"for counter in range(len(window_sizes)):\n",
" result = filtered_results3[counter,:,:]\n",
" \n",
" # TODO: find x and y, the coordinates of the largest element in the matrix \"result\"\n",
" x,y = \n",
" \n",
" print \"Peak for filtering with window size {} is at ({},{}).\".format(window_sizes[counter],x,y)"
],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}