程序代写代做代考 python cache algorithm Image_Analysis

Image_Analysis

Assignment 2 – Image Analysis in python¶

Setup¶

In [1]:

# Bring in necessary libraries
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
from skimage import color
import skimage.filters as filters
from skimage.transform import hough_circle
from skimage.feature import peak_local_max
from skimage import morphology
from skimage.draw import circle_perimeter
from skimage import img_as_float, img_as_ubyte
from skimage import segmentation as seg
from skimage.morphology import watershed
from scipy import ndimage as nd
from scipy.ndimage import convolve
from skimage import feature
import glob # for bulk file import

# Set defaults
plt.rcParams[‘image.cmap’] = ‘gray’ # Display grayscale images in… grayscale.
plt.rcParams[‘image.interpolation’] = ‘none’ # Use nearest-neighbour
plt.rcParams[‘figure.figsize’] = 10, 10

# Import test images
imgpaths = glob.glob(“./images/*.jpg”) + glob.glob(“./images/*.png”)
imgset = [img_as_ubyte(mpimg.imread(x)) for x in imgpaths]

# Display thumbnails of the images to ensure loading
plt.figure()
for i,img in enumerate(imgset):
plt.subplot(1, len(imgset), i+1)
plt.imshow(img, cmap = ‘gray’)

Part A – Color and Exposure¶

i. Create a function to plot a histogram for each color channel in a single plot for each of your images.¶

In [2]:

# Plots a histogram of the image, splitting into individual channels if necessary.
def plot_multichannel_histo(img):
if img.ndim == 2: # plot grayscale histo
plt.hist(img.flatten(), 256, range=(0,255), color=’k’, histtype=’step’)
elif img.ndim == 3: # print rgb histo
plt.hist(img.reshape(-1,3), 256, range=(0,255), color=[‘r’,’g’,’b’],histtype=’step’)
else: # Not an image
print(“Must pass a valid RGB or grayscale image”)
plt.xlim([0,255])

In [3]:

# Apply to test images
for i,img in enumerate(imgset):
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(img, cmap=’gray’)
plt.subplot(1, 2, 2)
plot_multichannel_histo(img)

ii. Display thresholded versions of your images¶

In [4]:

# Plot a per-channel thresholded image and its corresponding histogram
def plot_threshold_with_histo(img, thresh):
plt.subplot(1, 2, 1)
plt.imshow(img > thresh)
plt.subplot(1, 2, 2)
plot_multichannel_histo(img)
plt.axvline(thresh, color=’red’);

In [5]:

# Apply to test images
for i,img in enumerate(imgset):
plt.figure()
plot_threshold_with_histo(img, 128)

—————————————————————————
ValueError Traceback (most recent call last)
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/IPython/core/formatters.pyc in __call__(self, obj)
305 pass
306 else:
–> 307 return printer(obj)
308 # Finally look for special method names
309 method = get_real_method(obj, self.print_method)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in (fig)
225
226 if ‘png’ in formats:
–> 227 png_formatter.for_type(Figure, lambda fig: print_figure(fig, ‘png’, **kwargs))
228 if ‘retina’ in formats or ‘png2x’ in formats:
229 png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in print_figure(fig, fmt, bbox_inches, **kwargs)
117
118 bytes_io = BytesIO()
–> 119 fig.canvas.print_figure(bytes_io, **kw)
120 data = bytes_io.getvalue()
121 if fmt == ‘svg’:

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
2190 orientation=orientation,
2191 dryrun=True,
-> 2192 **kwargs)
2193 renderer = self.figure._cachedRenderer
2194 bbox_inches = self.figure.get_tightbbox(renderer)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
543
544 def print_png(self, filename_or_obj, *args, **kwargs):
–> 545 FigureCanvasAgg.draw(self)
546 renderer = self.get_renderer()
547 original_dpi = renderer.dpi

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in draw(self)
462
463 try:
–> 464 self.figure.draw(self.renderer)
465 finally:
466 RendererAgg.lock.release()

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
61 def draw_wrapper(artist, renderer, *args, **kwargs):
62 before(artist, renderer)
—> 63 draw(artist, renderer, *args, **kwargs)
64 after(artist, renderer)
65

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/figure.pyc in draw(self, renderer)
1141
1142 mimage._draw_list_compositing_images(
-> 1143 renderer, self, dsu, self.suppressComposite)
1144
1145 renderer.close_group(‘figure’)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
137 if not_composite or not has_images:
138 for zorder, a in dsu:
–> 139 a.draw(renderer)
140 else:
141 # Composite any adjacent images together

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
61 def draw_wrapper(artist, renderer, *args, **kwargs):
62 before(artist, renderer)
—> 63 draw(artist, renderer, *args, **kwargs)
64 after(artist, renderer)
65

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe)
2407 renderer.stop_rasterizing()
2408
-> 2409 mimage._draw_list_compositing_images(renderer, self, dsu)
2410
2411 renderer.close_group(‘axes’)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
137 if not_composite or not has_images:
138 for zorder, a in dsu:
–> 139 a.draw(renderer)
140 else:
141 # Composite any adjacent images together

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
61 def draw_wrapper(artist, renderer, *args, **kwargs):
62 before(artist, renderer)
—> 63 draw(artist, renderer, *args, **kwargs)
64 after(artist, renderer)
65

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in draw(self, renderer, *args, **kwargs)
493 else:
494 im, l, b, trans = self.make_image(
–> 495 renderer, renderer.get_image_magnification())
496 if im is not None:
497 renderer.draw_image(gc, l, b, im)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in make_image(self, renderer, magnification, unsampled)
717 return self._make_image(
718 self._A, bbox, transformed_bbox, self.axes.bbox, magnification,
–> 719 unsampled=unsampled)
720
721 def _check_unsampled_image(self, renderer):

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification, unsampled, round_to_pixel_border)
409 A, output, t, _interpd_[self.get_interpolation()],
410 self.get_resample(), alpha,
–> 411 self.get_filternorm() or 0.0, self.get_filterrad() or 0.0)
412
413 if created_rgba_mask:

ValueError: 3-dimensional arrays must be of dtype unsigned byte, unsigned short, float32 or float64

—————————————————————————
ValueError Traceback (most recent call last)
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/IPython/core/formatters.pyc in __call__(self, obj)
305 pass
306 else:
–> 307 return printer(obj)
308 # Finally look for special method names
309 method = get_real_method(obj, self.print_method)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in (fig)
225
226 if ‘png’ in formats:
–> 227 png_formatter.for_type(Figure, lambda fig: print_figure(fig, ‘png’, **kwargs))
228 if ‘retina’ in formats or ‘png2x’ in formats:
229 png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/IPython/core/pylabtools.pyc in print_figure(fig, fmt, bbox_inches, **kwargs)
117
118 bytes_io = BytesIO()
–> 119 fig.canvas.print_figure(bytes_io, **kw)
120 data = bytes_io.getvalue()
121 if fmt == ‘svg’:

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/backend_bases.pyc in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
2190 orientation=orientation,
2191 dryrun=True,
-> 2192 **kwargs)
2193 renderer = self.figure._cachedRenderer
2194 bbox_inches = self.figure.get_tightbbox(renderer)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in print_png(self, filename_or_obj, *args, **kwargs)
543
544 def print_png(self, filename_or_obj, *args, **kwargs):
–> 545 FigureCanvasAgg.draw(self)
546 renderer = self.get_renderer()
547 original_dpi = renderer.dpi

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/backends/backend_agg.pyc in draw(self)
462
463 try:
–> 464 self.figure.draw(self.renderer)
465 finally:
466 RendererAgg.lock.release()

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
61 def draw_wrapper(artist, renderer, *args, **kwargs):
62 before(artist, renderer)
—> 63 draw(artist, renderer, *args, **kwargs)
64 after(artist, renderer)
65

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/figure.pyc in draw(self, renderer)
1141
1142 mimage._draw_list_compositing_images(
-> 1143 renderer, self, dsu, self.suppressComposite)
1144
1145 renderer.close_group(‘figure’)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
137 if not_composite or not has_images:
138 for zorder, a in dsu:
–> 139 a.draw(renderer)
140 else:
141 # Composite any adjacent images together

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
61 def draw_wrapper(artist, renderer, *args, **kwargs):
62 before(artist, renderer)
—> 63 draw(artist, renderer, *args, **kwargs)
64 after(artist, renderer)
65

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/axes/_base.pyc in draw(self, renderer, inframe)
2407 renderer.stop_rasterizing()
2408
-> 2409 mimage._draw_list_compositing_images(renderer, self, dsu)
2410
2411 renderer.close_group(‘axes’)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in _draw_list_compositing_images(renderer, parent, dsu, suppress_composite)
137 if not_composite or not has_images:
138 for zorder, a in dsu:
–> 139 a.draw(renderer)
140 else:
141 # Composite any adjacent images together

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/artist.pyc in draw_wrapper(artist, renderer, *args, **kwargs)
61 def draw_wrapper(artist, renderer, *args, **kwargs):
62 before(artist, renderer)
—> 63 draw(artist, renderer, *args, **kwargs)
64 after(artist, renderer)
65

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in draw(self, renderer, *args, **kwargs)
493 else:
494 im, l, b, trans = self.make_image(
–> 495 renderer, renderer.get_image_magnification())
496 if im is not None:
497 renderer.draw_image(gc, l, b, im)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in make_image(self, renderer, magnification, unsampled)
717 return self._make_image(
718 self._A, bbox, transformed_bbox, self.axes.bbox, magnification,
–> 719 unsampled=unsampled)
720
721 def _check_unsampled_image(self, renderer):

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/matplotlib/image.pyc in _make_image(self, A, in_bbox, out_bbox, clip_bbox, magnification, unsampled, round_to_pixel_border)
409 A, output, t, _interpd_[self.get_interpolation()],
410 self.get_resample(), alpha,
–> 411 self.get_filternorm() or 0.0, self.get_filterrad() or 0.0)
412
413 if created_rgba_mask:

ValueError: 3-dimensional arrays must be of dtype unsigned byte, unsigned short, float32 or float64

iii. Decompose each of your images into its H, S and V color components. Display each component and interpret the result.¶

In [6]:

# Show HSV as three plots
def plot_hsv_components(img):
if img.ndim == 2: # convert grayscale to rgb
hsv = color.rgb2hsv(color.gray2rgb(img))
elif img.ndim == 3:
hsv = color.rgb2hsv(img)
else: # Not an image
print(“Must pass a valid RGB or grayscale image”)
plt.subplot(1, 3, 1)
plt.title(‘Hue’)
plt.imshow(hsv[:,:,0], cmap=’hsv’) # Hue
plt.subplot(1, 3, 2)
plt.title(‘Saturation’)
plt.imshow(hsv[:,:,1], cmap=’Greens’) # Sat
plt.subplot(1, 3, 3)
plt.title(‘Value’)
plt.imshow(hsv[:,:,2], cmap=’gray’) # Value

In [7]:

# Apply to test images
for i,img in enumerate(imgset):
plt.figure()
plot_hsv_components(img)

Interpretation: rgb2hsv converts an image to the HSV colorspace, representing colors by their pute chroma, saturation/purity, and brightness. Certain operations and filters will no doubt be eaiser to perform in HSV space rather than RGB, particularly those that deal with changing color or brightness.

iv. Transform each of your images to the CIELAB color space. Plot the L, a, b components to get a feel for what they do. What do they do?¶

In [8]:

# Show Lab as three plots
def plot_Lab_components(img):
if img.ndim == 2: # convert grayscale to Lab
lab = color.rgb2lab(color.gray2rgb(img))
elif img.ndim == 3:
lab = color.rgb2lab(img)
else: # Not an image
print(“Must pass a valid RGB or grayscale image”)
plt.subplot(1, 3, 1)
plt.title(‘L’)
plt.imshow(lab[:,:,0]) # L
plt.subplot(1, 3, 2)
plt.title(‘a’)
plt.imshow(lab[:,:,1]) # a
plt.subplot(1, 3, 3)
plt.title(‘b’)
plt.imshow(lab[:,:,2]) # b

In [9]:

# Apply to test images
for i,img in enumerate(imgset):
plt.figure()
plot_Lab_components(img)

The ‘C’ component appears to be synonymous with V from HSV. The ‘a’ component seems to shift from green to red as it increases, and the ‘b’ component from blue to yellow. It is worth noting that unlike in HSV, where all but one component was zero, grayscale images converted to Lab show variance among the three components.

Additional Questions¶

1. Did the tool/algorithm do a good job? (e.g. did an edge detection algorithm detect edges?)¶

The algorithms used in this section all did exactly what they were intended to do.

2. Give an example of how the tool/algorithm might be used in real world applications. (e.g. the blob detection algorithm might be used for identifying stars.)¶

Histograms can be used to analyze the distribution of brightness in an image in order to improve contrast by scaling values, or highlight visually distinct regions though thresholds. HSV and CIELAB both provide different spaces for performing calculations in, simplifying the mathematics behind many color transformations and anaylsis operations. These techniques combined could be used to, for example, pick out objects of a certain color on a conveyer belt for sorting.

Part B – Image Filters¶

i. Create a kernel. Convolve the kernel with each of your images (see np.convolve). Display the convolved images. Your kernel must have a noticeable effect on the images. (Though it need not do anything useful)¶

In [10]:

# Create some arbitrary kernel
def generate_kernel(size):
kernel = np.zeros((size,size))
for x in range(size):
kernel[x][x] = 1.0
kernel[size-1-x][x] = 1.0
kernel = kernel / float(sum(sum(kernel)))
return kernel

kernel = generate_kernel(50)
plt.imshow(kernel)

Out[10]:

In [11]:

# Apply to test images
for i,img in enumerate(imgset):
imgbw = color.rgb2grey(img)
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(imgbw)
plt.subplot(1, 2, 2)
plt.imshow(convolve(imgbw, kernel))

ii. Downsample each of your images¶

In [12]:

# Downsample an image by skipping indicies
def decimate_image(img, skip):
return img[::skip,::skip]

In [13]:

# Apply to test images
for i,img in enumerate(imgset):
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(img)
plt.subplot(1, 2, 2)
plt.imshow(decimate_image(img, 40))

iii. Apply a Gaussian filter to each of your images. Display the filtered images.¶

In [14]:

# Apply to test images
sigma = 20.0
for i,img in enumerate(imgset):
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(img)
plt.subplot(1, 2, 2)

iv. Create a simple difference filter to find the horizontal or vertical edges for each of your images.¶

In [15]:

# Find horizontal edges using a simple shifting method
def find_horizontal_edges(img):
imgbw = img_as_float(color.rgb2grey(img))
return np.abs(imgbw[:, 1:] – imgbw[:, :-1])

# Find vertical edges using a simple shifting method
def find_vertical_edges(img):
imgbw = img_as_float(color.rgb2grey(img))
return np.abs(imgbw[1:, :] – imgbw[:-1, :])

In [16]:

# Apply to test images
for i,img in enumerate(imgset):
decimg = decimate_image(img, 5) # downsample to make it easier to see graphs
plt.figure()
plt.subplot(1, 3, 1)
plt.title(‘Original’)
plt.imshow(decimg)
plt.subplot(1, 3, 2)
plt.title(‘Horizontal Edges’)
plt.imshow(find_horizontal_edges(decimg))
plt.subplot(1, 3, 3)
plt.title(‘Vertical Edges’)
plt.imshow(find_vertical_edges(decimg))

v. Create a simple difference filter to visualize the pair-wise differences between each of your images.¶

In [17]:

# Find the absolute difference between two images. Crops to the shared region between images.
def find_pairwise_difference(img_a, img_b):
subset = np.minimum(img_a.shape, img_b.shape)
img_a_subset = img_a[:subset[0], :subset[1]]
img_b_subset = img_b[:subset[0], :subset[1]]
return img_a_subset, img_b_subset, np.abs(img_a_subset – img_b_subset)

In [18]:

# Apply to test images
for i in range(len(imgset)):
decimg_a = img_as_float(color.rgb2grey(decimate_image(imgset[i], 5))) # downsample to make it easier to see graphs
decimg_b = img_as_float(color.rgb2grey(decimate_image(imgset[(i+1) % len(imgset)], 5)))
a, b, d = find_pairwise_difference(decimg_a, decimg_b)
plt.figure()
plt.subplot(1, 3, 1)
plt.imshow(a)
plt.subplot(1, 3, 2)
plt.imshow(b)
plt.subplot(1, 3, 3)
plt.imshow(d)

Additional Questions¶

1. Did the tool/algorithm do a good job? (e.g. did an edge detection algorithm detect edges?)¶

The convolution, downsample, and Gaussian filters worked exactly as expected. The difference filter edge detection works acceptably, though it could be greatly improved. The pairwise difference function is hard to evaluate, though it seems mathematically correct. Results would be much clearer with image pairs with slight differences.

2. Give an example of how the tool/algorithm might be used in real world applications. (e.g. the blob detection algorithm might be used for identifying stars.)¶

Convolution is useful as a general tool for a variety of applications, such as smoothing, edge detection, and sharpening. Downsampling can reduce the computational time of filters or decrease size needed on disk when working with small, low resolution preview images. Gaussian filters are useful for looking at general trends in brightness and color across an image, as well as removing some types of noise. Edge detection can segment images into similar regions in order to detect objects. The pairwise difference function could prove a useful starting point for image comparison of similar images.

Part C – Feature Detection¶

i. Apply a Sobel filter to each of your images. Display the result.¶

In [19]:

# Apply a Sobel filter to the image.
def apply_sobel(img):
sobel_kernel = np.array([[-1, 0, +1],[-2, 0, +2],[-1, 0, +1]])
Gx = convolve(img, sobel_kernel)
Gy = convolve(img, -sobel_kernel.transpose())
return np.sqrt(np.square(Gx) + np.square(Gy))

# This should do the same thing as skimage.filters.sobel(img)

In [20]:

# Apply to test images
for i,img in enumerate(imgset):
imgbw = img_as_float(color.rgb2grey(decimate_image(img, 5))) # downsample to make it easier to see graphs
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(imgbw)
plt.subplot(1, 2, 2)
plt.imshow(apply_sobel(imgbw))

ii. Apply the Canny edge detector algorithm each of your images. Display the result.¶

In [21]:

# Apply to test images
sigma = 2.0
for i,img in enumerate(imgset):
imgbw = img_as_float(color.rgb2grey(decimate_image(img, 5))) # downsample to make it easier to see graphs
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(imgbw)
plt.subplot(1, 2, 2)
plt.imshow(filters.canny(imgbw, sigma))

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):

iii. How susceptible is the Canny edge detector to image noise? Back up your answer with examples using your images.¶

In [22]:

def add_noise(img, sigma=1.0):
if sigma == 0:
return img
else:
return img + np.random.normal(0, sigma, np.shape(img))

# Apply to test images
filter_sigma = 1.5
for i,img in enumerate(imgset):
imgbw = img_as_float(color.rgb2grey(decimate_image(img, 5))) # downsample to make it easier to see graphs
plt.figure()
for j in range(3): # make a row of increasing noise amounts
noise_sigma = 0.0 + 0.1 * j
plt.subplot(2, 3, j+1)
noisy = add_noise(imgbw, noise_sigma)
plt.imshow(noisy)
plt.subplot(2, 3, j+4)
plt.imshow(filters.canny(noisy, filter_sigma))

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):

The Canny filter seems particularly susceptible to noise, as it generates spurrious edges very quickly as noise is increased. A possible solution to this is applying a smoothing filter, such as a Gaussian blur, before hand, though this obviously runs the risk of losing valuable information.

iv. Apply Hough transforms to each of your images. Display the result.¶

In [23]:

# Plot the circular Hough transforms of an image at the given radii.
def plot_circle_hough(img, radii, sigma):
edges = feature.canny(img, sigma)
hough = hough_circle(edges, radii)
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(imgbw)
plt.subplot(1, 2, 2)
plt.imshow(edges)

for j in range(len(hough)):
plt.subplot(1, len(hough), j+1)
plt.imshow(hough[j,:,:])

In [24]:

# Apply to test images
plt.figure()
radii = np.arange(10, 35, 5)
sigma = 2.0
for i,img in enumerate(imgset):
imgbw = img_as_float(color.rgb2grey(decimate_image(img, 5))) # downsample to make it easier to see graphs
plot_circle_hough(imgbw, radii, sigma)

v. Use the response from the Hough transform to detect the position and radius of each circular objects in each of your images.¶

In [25]:

def plot_detected_circles(img, radii, sigma):
edges = filters.canny(img, sigma)
hough = hough_circle(edges, radii)

# Code after this point adapted from the scikit documentation. Not fully sure how it works.
accums = []
found_centers = []
found_radii = []

for radius, h in zip(radii, hough):
# For each radius, extract two circles
peaks = peak_local_max(h, num_peaks=2)
found_centers.extend(peaks)
accums.extend(h[peaks[:, 0], peaks[:, 1]])
found_radii.extend([radius, radius])

# Draw the most prominent 5 circles
image = color.gray2rgb(img)
for idx in np.argsort(accums)[::-1][:5]:
center_x, center_y = found_centers[idx]
radius = found_radii[idx]
cx, cy = circle_perimeter(center_y, center_x, radius)
image[cy, cx] = (220, 20, 20)

plt.imshow(image, cmap=plt.cm.gray)

In [26]:

# Apply to test images
radii = np.arange(10, 100, 2)
sigma = 2.0
for i,img in enumerate(imgset):
imgbw = img_as_float(color.rgb2grey(decimate_image(img, 5))) # downsample to make it easier to see graphs
plt.figure()
plot_detected_circles(imgbw, radii, sigma)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):

Additional Questions¶

1. Did the tool/algorithm do a good job? (e.g. did an edge detection algorithm detect edges?)¶

The Sobel operator did a better job than the differencefilter at detecting edges, though once again there is a lot of room for improvement. The Canny detector was far more accurate, though had trouble with sharp corners and small contours, and was easily distorted by the addition of noise. The Hough transform detection did not do a very good job at finding circles, though this is entirely understandable given the lack of defined circles in any of the input images. I was impressed at the relevancy regions the algorithm highlighted despite the poor source images.

2. Give an example of how the tool/algorithm might be used in real world applications. (e.g. the blob detection algorithm might be used for identifying stars.)¶

Once again, these feature detection are used to segment and detect objects in an image. Possible applications include factory line quality assurance, image-based currency sorting, and geographical feature detection for maps.

Part D – Morphological Operations¶

i. Use morphological operations to remove noise from each of your images.¶

In [27]:

# Applies a morphological operator to remove noise
def morpho_denoise(img, shape, size):
if shape == ‘square’:
kernel = morphology.square(width=size)
elif shape == ‘diamond’:
kernel = morphology.diamond(radius=size)
else:
print(“Shape must be ‘square’ or ‘diamond’.”)
return None
return morphology.opening(img, kernel)

In [28]:

# Apply to test images
for i,img in enumerate(imgset):
imgbw = img_as_float(color.rgb2grey(img))
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(imgbw)
plt.subplot(1, 2, 2)
plt.imshow(morpho_denoise(imgbw, ‘square’, 10))

Additional Questions¶

1. Did the tool/algorithm do a good job? (e.g. did an edge detection algorithm detect edges?)¶

The morphological operators did a fair job at removing the noise in the image, with only slight loss of detail. It is a marked improvement over the Gaussian filter for this purpose.

2. Give an example of how the tool/algorithm might be used in real world applications. (e.g. the blob detection algorithm might be used for identifying stars.)¶

This technique is well-suited as a pre-process for later analysis, such as edge detection, to remove the errors caused by noise without sacrificing too much detail.

Part E – Segmentation¶

i. Apply SLIC (Simple Linear Iterative Clustering) to each of your images. Display the result. What does SLIC do¶

In [29]:

# Calculate the mean color of slic regions, from the SciKit tutorial
def mean_color(image, labels):
out = np.zeros_like(image)
for label in np.unique(labels):
indices = np.nonzero(labels == label)
out[indices] = np.mean(image[indices], axis=0)
return out

def plot_slic_segmentation(img):
labels = seg.slic(img, n_segments=24, compactness=70, sigma=2.0, enforce_connectivity=True)
return mean_color(img, labels)

In [30]:

# Apply to the images
for i,img in enumerate(imgset):
rgbimg = img_as_float(color.gray2rgb(img))
plt.figure()
plt.subplot(1, 2, 1)
plt.imshow(img)
plt.subplot(1, 2, 2)
plt.imshow(plot_slic_segmentation(rgbimg))

SLIC find regions of the image that are of similar color and attempts to group them into a number of contiguous regions.

ii. Apply the watershed algorithm to each of your images. Display the result. What does the watershed algorithm do?¶

In [31]:

# Show a watershed plot of an image
def plot_watershed(img, sigma):
edges = filters.canny(img, sigma)
distance_from_edge = nd.distance_transform_edt(-edges)
peaks = peak_local_max(distance_from_edge)
peaks_image = np.zeros(img.shape, np.bool)
peaks_image[tuple(np.transpose(peaks))] = True
seeds, num_seeds = nd.label(peaks_image)
plt.plot(peaks[:, 1], peaks[:, 0], ‘ro’)
ws = watershed(edges, seeds)
plt.imshow(color.label2rgb(ws, img))

In [32]:

# Apply to the images
for i,img in enumerate(imgset):
imgbw = (color.rgb2grey(decimate_image(img, 5)))
plt.figure()
plot_watershed(imgbw, 1.0)

/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:4: DeprecationWarning: numpy boolean negative, the `-` operator, is deprecated, use the `~` operator or the logical_not function instead.
/Users/DDCR/Documents/Neu/3.5/CS4300/anaconda/lib/python2.7/site-packages/skimage/filters/__init__.py:28: skimage_deprecation: Call to deprecated function “canny“. Use “skimage.feature.canny“ instead.
def canny(*args, **kwargs):

Additional Questions

1. Did the tool/algorithm do a good job? (e.g. did an edge detection algorithm detect edges?)¶

Both of these algorithms did an acceptable, though a bit underwhelming, job at segmenting the images provided. As was the case with the feature detection algorithms, it is likely that this is due to poor source images containing gradients rather than discrete color regions.

2. Give an example of how the tool/algorithm might be used in real world applications. (e.g. the blob detection algorithm might be used for identifying stars.)¶

As with feature detection, these algorithms would serve well as a first step in part of a larger image segmentation process with potential applications including identifying discrete objects in images, medical imaging (locating abnomalies, positioning organs), and nagivation for autonomous vehicles.

Part F – Sepia or Color Filter
i:Write a matrix transformation that will create a sepia or color effect to your images.

In [ ]:

def sepia (image, grayscale):
if red < 63: red = int(red * 1.1) blue = int(blue * 0.9) elif red < 192: red = int(red * 1.15) blue = int(blue * 0.85) else: red = min(int(red * 1.08), 255) blue = int(blue * 0.93) def grayscale(image): for y in xrange(image.getHeight()): for x in xrange(image.getWidth()): (r, g, B)/> = image.getPixel(x, y)
r = int(r * 0.299)
g = int(g * 0.587)
b = int(b * 0.114)
lum = r + g + b
image.setPixel(x, y, (lum, lum, lum))

In [ ]:

# Apply to Image

Part G
i. Write a convolution that will sharpen your images.