View source for Pioneer21/Lecture03

This page serves as a supplement to the third group meeting of the Summer 2021 Pioneer program for image processing.

== Python Requirements ==
The programs from this lecture will require several modules:
* NumPy - [https://numpy.org/ https://numpy.org/]
* Matplotlib - [https://matplotlib.org/ https://matplotlib.org/] (and specifically matplotlib.pyplot)
* scikit-image - [https://scikit-image.org/ https://scikit-image.org/]
* SciPy - [https://www.scipy.org/ https://www.scipy.org/] (and specifically scipy.signal)
The start of any script we use in this lecture will be:
<syntaxhighlight lang=python>
import numpy as np
import matplotlib.pyplot as plt
import skimage as ski
import scipy.signal as sig
</syntaxhighlight>

== scikit-image commands ==
The main command for this session is:
* <code>DATA = ski.data.NAME()</code> where DATA is a variable name you give to store the information from a built-in image and NAME is the name of a built-in image from scikit-image; the built-in images can be found in the Data section of the [https://scikit-image.org/docs/stable/auto_examples/index.html General Examples] page.

'''Starting on July 26, 2021, we will be using <code>AXIS.imshow(DATA)</code> instead of  <code>ski.io.imshow(DATA)</code> to display an image.  It looks like ski.io.imshow is going away.'''

Depending on the nature of an image, its DATA variable will be different shapes and contain different data types.  Here is how scikit-image and imshow() interpret things:
* A 2D-array will either be black and white, grayscale,  or it will be mapped to a colormap
* A 3D-array will be interpreted as an RGB-based color image
* For a 2D-array, the following data types and sizes might be useful:
** An array of boolean (True or False) values or integer values with only integers 0 or 1 will produce an image with values at the extreme edges of the colormap.  Assume that the code:<syntaxhighlight lang=python>
bwimage = ski.data.binary_blobs()
</syntaxhighlight>  has been run already for each of the cases below.  By default in Python with matplotlib this is purple and gold.
*** Default case with booleans <syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
ax.imshow(bwimage)</syntaxhighlight>  You can add some keyword arguments to the command to get black and white as follows:
*** Black and white with booleans<syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
ax.imshow(bwimage, cmap=plt.cm.gray)</syntaxhighlight>
*** Black and white and no axis ticks with booleans<syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
ax.imshow(bwimage, cmap=plt.cm.gray)
ax.set_axis_off()</syntaxhighlight>Note that for each of the above, you can replace the array of booleans with an array of integers 0 or 1 by replacing the <code>bwimage</code> line with:<syntaxhighlight lang=python>
bwimage = np.random.choice([0,1], (50,100))
</syntaxhighlight> or <syntaxhighlight lang=python>
bwimage = np.random.randint(0, 2, (50, 100))
</syntaxhighlight> 
** An array of integer values will produce an image with values mapped to the default Viridis (purple - green - blue - yellow) colormap.  Assume that the code:<syntaxhighlight lang=python>
grayimage = ski.data.coins()
</syntaxhighlight>  has been run already for each of the cases below.
*** Default case with integers <syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
ax.imshow(grayimage)</syntaxhighlight> You can add some keyword arguments to the command to get shades of gray as above:
*** Gray with integers <syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
ax.imshow(grayimage, cmap=plt.cm.gray)</syntaxhighlight>
*** Gray and no axis ticks with integers<syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
ax.imshow(grayimage, cmap=plt.cm.gray)
ax.set_axis_off()</syntaxhighlight>
:::*'''Note 1''' - the addition of the colormap is required for using the axis version of imshow versus the scikit-image version of imshow. 
:::*'''Note 2''' - by default, imshow will stretch the range of integers to cover the whole map.  If you want to specifically map to a certain range - for instance, have 0 as pure black and 255 as pure white, you need to add <code>vmin</code> and <code>vmax</code> keyword arguments to explicitly give the values that map to the low and high end of the volormap.  For instance:<syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
grayimage = ski.data.coins()
ax.imshow(grayimage, cmap=plt.cm.gray, vmin=0, vmax=255))</syntaxhighlight> will produce the exact same image that ski.io.imshow(grayimage) would have produced.  
** An array of floating-point values will also produce an image with values mapped to the default Viridis (purple - green - blue - yellow) colormap.  You can add some keyword arguments to the command to get shades of gray as above.  Assume that the code:<syntaxhighlight lang=python>
x, y = np.meshgrid(np.linspace(-2, 2, 201), np.linspace(-1, 1, 101))
z = np.exp(-np.sqrt(x**2+y**2))*np.cos(2*np.pi*x)*np.sin(4*np.pi*y)
</syntaxhighlight>  has been run already for each of the cases below.
*** Default case with floats<syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
ax.imshow(z)</syntaxhighlight>
*** Gray with floats<syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
grayimage = ski.data.coins()
ax.imshow(grayimage, cmap=plt.cm.gray)</syntaxhighlight>
*** Gray and no axis ticks with floats<syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
grayimage = ski.data.coins()
ax.imshow(grayimage, cmap=plt.cm.gray)
ax.set_axis_off()</syntaxhighlight>
:::*'''Note''' - once again, by default, imshow will stretch the range of floats to cover the whole map.  If you want to specifically map to a certain range - for instance, have 0 as pure black and 0.2 as pure white, you need to add <code>vmin</code> and <code>vmax</code> keyword arguments to explicitly give the values that map to the low and high end of the colormap.  Anything below vmin will be mapped to the low color and anything above vmax will be mapped to the high.  For instance:<syntaxhighlight lang=python>
fig, ax = plt.subplots(num=1, clear=True)
ax.imshow(z, vmin=0, vmax=255))</syntaxhighlight> will produce an image with yellow rounded rectangles on a purple background with greenish-blue edges the further you get away from the center.  The large all-yellow rectangles contain many values that are above 0.2 and the purple parts are all the negative parts of the array.

<!-- 
== Image Processing Toolbox ==
For this lab, there are a few commands you will need to learn (or remember):
* <code>imread('filename.ext')</code> will load an image into MATLAB.  Depending on the image type, the return from this function may be a matrix of binary numbers (black and white images), an array of integers between 0 and 255 (grayscale images), or a three-layer matrix of integers between 0 and 255 (color images).
* <code>image(matrix)</code> will display the contents of a matrix as viewed from above.
** By default, the <code>image</code> command for a 1-layer matrix will assign colors based on a colormap that spans values from 0 to 255.  Anything outside of that range will be clipped.  
** There are ways to change that, but you will not need to for this lab
** The <code>image</code> command for a 3-layer matrix will assign colors based on the first layer being the red component, the second layer being the green component, and the third layer being the blue component. 
*** If the matrix is made up of floating-point numbers, <code>image</code> expects those numbers to be in the range [0, 1]
*** If the matrix is made up of unsigned integers, <code>image</code> expects those numbers to be in the range [0, 255]
*** If the matrix is made up of signed integers, <code>image</code> expects those numbers to be in the range [-128, 127]
* <code>imagesc(matrix)</code> for a one layer matrix will display the contents of a matrix as viewed from above and will also map the minimum value to "color 0" and the maximum value to "color 255"; for a three-layer matrix it will work just like <code>image</code>
* <code>colorbar</code> will add a...colorbar to the right of an image to show the numerical values and colors associated with them.  This is only useful for a one-layer image.
* <code>colormap MAP</code> will assign a particular colormap.  For this assignment, the most useful one is '''gray''' for grayscale images.
* <code>axis equal</code> will tell MATLAB to display an item such that each direction has the same scale.  For the imaging commands, this is useful in that it will make each pixel a square regardless of the shape of the figure window, thus preserving the geometry of an image.

== Links ==
* 2D Convolution link for 6.1: [HTTP://www.songho.ca/dsp/convolution/convolution2d_example.HTML HTTP://www.songho.ca/dsp/convolution/convolution2d_example.HTML]
* Kernel link for Exercise 7: [https://en.wikipedia.org/wiki/Kernel_(image_processing) https://en.wikipedia.org/wiki/Kernel_(image_processing)]
* Test Card for Exercise 8: [https://en.wikipedia.org/wiki/Test_card https://en.wikipedia.org/wiki/Test_card]

== Examples ==
The following sections will contain both the example programs given in the lab as well as the image or images they produce.  You should still type these into your own version of MATLAB to make sure you are getting the same answers.  These are provided so you can compare what you get with what we think you should get.

=== Example 1: Black & White Images===
[[File:IP1 E1 Plot1.png|thumb]]
<syntaxhighlight lang=matlab>
a = [ 1 0 1 0 0; ...
      1 0 1 0 1; ...
      1 1 1 0 0; ...
      1 0 1 0 1; ...
      1 0 1 0 1 ];
figure(1); clf
imagesc(a)
colormap gray; colorbar
</syntaxhighlight>
<br clear=all>

=== Example 2: Simple Grayscale Images===
[[File:IP1 E2 Plot1.png|thumb]]
<syntaxhighlight lang=matlab>
b = 0:255;
figure(1); clf
image(b) 
colormap gray; colorbar
</syntaxhighlight>
<br clear=all>

=== Example 3: Less Simple Grayscale Images===
[[File:IP1 E3 Plot1.png|thumb]]
<syntaxhighlight lang=matlab>
[x, y] = meshgrid(linspace(0, 2*pi, 201));
z = cos(x).*cos(2*y);
figure(1); clf
imagesc(z)
axis equal; colormap gray; colorbar
</syntaxhighlight>
Notice how the use of <code>axis equal</code> made the image look like a square since it is 201x201 but also caused the display to be filled with whitespace as a result of the figure size versus the image size.
<br clear=all>

=== Example 4: Building an Image===
[[File:IP1 E4 Plot1.png|thumb]]
<syntaxhighlight lang=matlab>
rad = 100;
del = 10;
[x, y] = meshgrid((-3*rad-del):(3*rad+del));
[rows, cols] = size(x);
dist = @(x, y, xc, yc) sqrt((x-xc).^2+(y-yc).^2);
venn_img = zeros(rows, cols, 3);
venn_img(:,:,1) = (dist(x, y, rad.*cos(0), rad.*sin(0)) < 2*rad);
venn_img(:,:,2) = (dist(x, y, rad.*cos(2*pi/3), rad.*sin(2*pi/3)) < 2*rad);
venn_img(:,:,3) = (dist(x, y, rad.*cos(4*pi/3), rad.*sin(4*pi/3)) < 2*rad);
figure(1); clf
image(venn_img)
axis equal
</syntaxhighlight>
<br clear=all>

=== Example 5: Exploring Colors===
[[File:IP1 E5 Plot1.png|thumb]]
<syntaxhighlight lang=matlab>
[x, y] = meshgrid(linspace(0, 1, 256));
other = 0.5;
palette = zeros(256, 256, 3);
palette(:,:,1) = x;
palette(:,:,2) = y;
palette(:,:,3) = other;
figure(1); clf
imagesc(palette)
axis equal
</syntaxhighlight>
<br clear=all>

=== Example 6: 1D Convolution===
<syntaxhighlight lang=matlab>
x = [1, 2, 4, 8, 7, 5, 1]
h = [1, -1]
y = conv(x, h)
</syntaxhighlight>
<syntaxhighlight lang=matlab>
x =

     1     2     4     8     7     5     1


h =

     1    -1


y =

     1     1     2     4    -1    -2    -4    -1
</syntaxhighlight>
<br clear=all>

=== Example 7: 1D Convolution Using 'same'===

<syntaxhighlight lang=matlab>
x = [1, 2, 4, 8, 7, 5, 1]
h = [1, -1]
y = conv(x, h, 'same')
</syntaxhighlight>
<syntaxhighlight lang=matlab>
x =

     1     2     4     8     7     5     1


h =

     1    -1


y =

     1     2     4    -1    -2    -4    -1
</syntaxhighlight>
<br clear=all>

=== Example 8: 10x10 Blurring===
[[File:IP1 E8 Plot1.png|thumb]]
[[File:IP1 E8 Plot2.png|thumb]]
<syntaxhighlight lang=matlab>
x = imread('coins.png');
h = ones(10, 10)/10^2;
y = conv2(x, h, 'same');
figure(1); clf
image(x)
axis equal; colormap gray; colorbar
title('Original')
figure(2); clf
image(y)
axis equal; colormap gray; colorbar
title('10x10 Blur')
</syntaxhighlight>
<br clear=all>

=== Example 9: Make No Assumptions===
<syntaxhighlight lang=matlab>
x = [1, 2, 4, 8, 7, 5, 1]
h = [1, -1]
y = conv(x, h, 'valid')
</syntaxhighlight>
<syntaxhighlight lang=matlab>
x =

     1     2     4     8     7     5     1


h =

     1    -1


y =

     1     2     4    -1    -2    -4
</syntaxhighlight>
<br clear=all>

=== Example 10: Basic Edge Detection and Display===
<syntaxhighlight lang=matlab>
[x, y] = meshgrid(linspace(-1, 1, 200));
z1 = (.7<sqrt(x.^2+y.^2)) & (sqrt(x.^2+y.^2)<.9);
z2 = (.3<sqrt(x.^2+y.^2)) & (sqrt(x.^2+y.^2)<.5);
zimg = 100*z1+200*z2;
figure(1); clf
image(zimg); axis equal; colormap gray; colorbar; title('Original')

hx = [1 -1; 1 -1];
edgex = conv2(zimg, hx, 'valid');
figure(2); clf
imagesc(edgex); axis equal; colormap gray; colorbar; title('Vertical Edges')

hy = hx';
edgey = conv2(zimg, hy, 'valid');
figure(3); clf
imagesc(edgey); axis equal; colormap gray; colorbar; title('Horizontal Edges')

edges = sqrt(edgex.^2 + edgey.^2);
figure(4); clf
imagesc(edges); axis equal; colormap gray; colorbar; title('Edges')
</syntaxhighlight>
<gallery>
File:IP1 E10 Plot1.png|Original
File:IP1 E10 Plot2.png|Vertical Edges
File:IP1 E10 Plot3.png|Horizontal Edges
File:IP1 E10 Plot4.png|Edges
</gallery>
<br clear=all>

=== Example 11: Chips!===
<syntaxhighlight lang=matlab>
clear
img = imread('coloredChips.png');
figure(1); clf
title('Original')
image(img); axis equal
vals = (0:255)'/255;
names = {'Red', 'Green', 'Blue'}
for k = 1:3
    figure(k+1); clf
    image(img(:,:,k)); axis equal
    colormap gray; colorbar
    title(names{k}+" as Gray")
    figure(k+4)
    image(img(:,:,k)); axis equal
    cmap = zeros(256, 3);
    cmap(:,k) = vals;
    colormap(cmap); colorbar
    title(names{k}+" as "+names{k})
end
</syntaxhighlight>
<gallery>
File:IP1 E11 Plot1.png|Original
</gallery>
<gallery>
File:IP1 E11 Plot2.png|Red as Gray
File:IP1 E11 Plot3.png|Green as Gray
File:IP1 E11 Plot4.png|Blue as Gray
</gallery>
<gallery>
File:IP1 E11 Plot5.png|Red as Red
File:IP1 E11 Plot6.png|Green as Green
File:IP1 E11 Plot7.png|Blue as Blue
</gallery>
<br clear=all>

=== Example 12: Chip Edges!===
<syntaxhighlight lang=matlab>
clear
img = imread('coloredChips.png');
figure(1); clf
image(img); axis equal

h = [1 0 -1; 2 0 -2; 1 0 -1]
for k=1:3
    y(:,:,k) = conv2(img(:,:,k), h, 'valid');
    figure(1+k); clf
    imagesc(y(:,:,k), [-1020, 1020]); 
    axis equal; colormap gray; colorbar
    figure(4+k); clf
    imagesc(abs(y(:,:,k)), [0, 1020]); 
    axis equal; colormap gray; colorbar
end
%%
yx = (y+max(abs(y(:)))) / 2 / max(abs(y(:)));
figure(8); clf
image(yx); axis equal

yg = sqrt(y(:,:,1).^2 + y(:,:,2).^2 + y(:,:,3).^2);
ygs = abs(yg) / max(yg(:)) * 255;
figure(9); clf
image(ygs); colormap gray; axis equal
</syntaxhighlight>
<gallery>
File:IP1 E12 Plot1.png|Original
</gallery>
<gallery>
File:IP1 E12 Plot2.png|Red Edges as Gray
File:IP1 E12 Plot3.png|Green Edges as Gray
File:IP1 E12 Plot4.png|Blue Edges as Gray
</gallery>
<gallery>

File:IP1 E12 Plot5.png|Absolute Red Edges as Gray
File:IP1 E12 Plot6.png|Absolute Green Edges as Gray
File:IP1 E12 Plot7.png|Absolute Blue Edges as Gray
</gallery>
<gallery>
File:IP1 E12 Plot8.png|Colorful Edges
File:IP1 E12 Plot9.png|Absolute Edges
</gallery>
<br clear=all>

== Exercise Starter Codes ==
=== Exercise 3 ===
<syntaxhighlight lang=matlab>
tc = linspace(0, 1, 101);
xc = humps(tc);
td = linspace(0, 1, 11);
xd = humps(td);
figure(1); clf
plot(tc, xc, 'b-')
hold on
plot(td, xd, 'bo')
hold off
</syntaxhighlight>

=== Exercise 4 ===
<syntaxhighlight lang=matlab>
tc = linspace(0, 1, 101);
xc = humps(tc);
deltatc = tc(2)-tc(1);
td = linspace(0, 1, 11);
xd = humps(td);
deltatd = td(2)-td(1);

figure(1); clf
plot(tc, xc, 'b-')
hold on
plot(td, xd, 'bo')
hold off
title('Values')

figure(2); clf
twopointdiff = diff(xc)/deltatc;
twopointdiff(end+1)=twopointdiff(end);
plot(tc, twopointdiff, 'b-')
hold on
plot(tc(1:10:end), twopointdiff(1:10:end), 'bo')
hold off
title('Change Me')
</syntaxhighlight>
-->