<jupyter_start><jupyter_text>Osnabrück University - Computer Vision (Winter Term 2020/21) - Prof. Dr.-Ing. <NAME>, <NAME>, <NAME>, <NAME>  Exercise Sheet 11: Object recognition: PCA and Interest Points  IntroductionThis week's sheet should be solved and handed in before the end of **Saturday, January 30, 2021**. If you need help (and Google and other resources were not enough), feel free to contact your groups' designated tutor or whomever of us you run into first. Please upload your results to your group's Stud.IP folder.  Exercise 1: Pattern Recognition and PCA [4 points] **a)** What are the goals of *pattern recognition*? How can they be achieved? What are the problems? The goal of pattern recognition is to get a symbolic description from sensor data.It can be achieved following the standard architecture:* world -> sensor -> feature extraction -> classificationProblems:* no sharp description of classes and incomplete data **b)** What is *principal component analysis*? How is it related to pattern recognition? With PCA we can find a new basis for our dataset. The basis vectors can be ordered by maximun variance / minimun reconstruction error. In the best case each basis vector corresponds to a interpretable feature. **c)** Explain how principal components can be computed? Then implement a function that performs the computation. Principal components are the eigenvalues of the covariance matrix of a dataset. The eigenvalues give the ordering of the components.<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np


def pca(data):
    """
    Perform principal component analysis.
    
    Args:
        data (ndarray): an array of shape (n,k),
        meaning n entries with k dimensions
        
    Returns: two arrays
        pc (ndarray): array of shape (k,k) holding the principal components in its columns.
        var (ndarray): k-vector holding the corresponding variances, in descending order.
    """
    
    # BEGIN SOLUTION

    # Center the data (subtract the mean along columns).
    centered_data = data - data.mean(axis=0)

    # Compute the covariance matrix
    covariance_matrix = centered_data.T @ centered_data / (len(data) - 1)
    
    # Computing eigenvalues and eigenvectors of covariance matrix.
    # Attention: not always sorted.
    eigenvals, eigenvecs = np.linalg.eig(covariance_matrix)
    
    # Sort in order of descending magnitude
    idx = np.argsort(eigenvals)[::-1]
    pc, var = eigenvecs[:,idx], eigenvals[idx]
    
    # END SOLUTION
    return pc, var

# generate some random data
np.random.seed(23)
data = np.random.multivariate_normal([0,0], cov = [[1, .55], [.55, .5]], size=300)

# compute the principal components
pc, var = pca(data)
mean = data.mean(axis=0)

# plot the results
plt.figure(figsize=(12,8))
plt.xlim(-4,4)
plt.ylim(-2.5,2.5)
plt.scatter(*data.T)
plt.quiver(*mean[np.newaxis].repeat(2,axis=0).T, *(np.sqrt(var)*pc), color='red', scale=1, scale_units='xy')
plt.show()

# sanity check
assert np.allclose(var, [1.216, 0.137], rtol=1e-3)<jupyter_output><empty_output><jupyter_text>Exercise 2: Eigenfaces [6 points] **a)** Import the images from the directory `images/trainimgs` into an numpy array using the function `read_images_from_directory` provided in the cell below. Display the images and the corresponding names.<jupyter_code>%matplotlib inline
import sys
import os
import glob
import numpy as np
import matplotlib.pyplot as plt

def read_images_from_directory(directory, suffix, shape):
    """
    Read all images found in DIRECTORY with given file
    name SUFFIX. All images should have the same SHAPE,
    specified as (rows,columns).
    
    Args:
        directory (string): Name of input directory.
        suffix (string): File type suffix.
        shape (tuple): Shape of images to be loaded.
    
    Returns:
        images (ndarray): A numpy array of shape m*rows*columns (from shape)
        names (list): A list of corresponding image names.
    """

    # initialize the image array and name list
    #images = np.empty((0, *shape))
    images = np.empty((0, ) + shape)
    names = []

    # now loop through all image files in the directory
    for file_name in glob.glob(directory + os.sep + '*.' + suffix):
        if os.path.isfile(file_name):

            # load each image (as double)
            img = plt.imread(file_name)

            # check for correct size
            if img.shape == shape:
                images = np.append(images, img.reshape((1, ) + shape), axis=0)
                names.append(os.path.basename(file_name))
            else:
                print(
                    'warning: Image "' + file_name +
                    '" with wrong size will be ignored!',
                    file=sys.stderr)

    return images, names


# image file suffix
suffix = 'pgm'

# image size
img_shape = (192, 168)

# BEGIN SOLUTION
train_imgs, train_names = read_images_from_directory('trainimg', 'pgm', img_shape)

plt.figure(figsize=(12, 12))
plt.gray()
for i, n in enumerate(train_names):
    plt.subplot(5, 4, i + 1)
    plt.axis('off')
    plt.imshow(train_imgs[i])
    plt.title(n[5:7])
plt.show()
# END SOLUTION<jupyter_output><empty_output><jupyter_text>**b)** Use PCA to compute the eigenfaces (i.e. the eigenvectors of the face images). You may use your PCA function from Exercise 1c or some build in function. Explain what kind of input PCA expects, and how that fits to our images (you may have to `reshape` the images!). Finally, display the eigenfaces.<jupyter_code># Transform images into one-dimensional vectors.
# This will create a two dimension matrix of shape (n,k)
# with n being the number of images and k the number of
# pixels per image.
face_vecs = train_imgs.reshape((-1, np.prod(img_shape)))

# Compute the mean face.
mean_face = face_vecs.mean(axis=0)

# now we can project our images onto these eigenvectors,
# obtaining n eigenfaces (each of dimensionality k):
centered_faces = face_vecs - mean_face

# As in our case k is much greater than n, we swap the dimensions when performing PCA:
# Instead of computing the gigantic (k*k, k=32256) covariance matrix
#   Cov = 1/(N-1) X^T * X
# we compute the much smaller matrix (n*n, n=20)
#   M = 1/(N-1) X * X^T
matrix = centered_faces @ centered_faces.T / (len(centered_faces)-1)
    
# Computing eigenvalues and eigenvectors of covariance matrix.
# Attention: not always sorted.
eigenvals, eigenvecs = np.linalg.eigh(matrix)

# Sort in order of descending magnitude
idx = np.argsort(eigenvals)[::-1]
eigenvecs, eigenvals = eigenvecs[:,idx], eigenvals[idx]

# Now we can obtain the eigenfaces by projecting back into the original space.
eigenfaces = eigenvecs.T @ centered_faces

# Normalize the eigenfaces
eigenfaces = eigenfaces/np.linalg.norm(eigenfaces, axis=1)[:,np.newaxis]

assert len(eigenfaces) == len(face_vecs)

# Plot eigenvalues and eigenvectors
plt.figure(figsize=(12,5))
plt.gray();
plt.subplot(1,3,1); plt.title("The mean face")
plt.imshow(mean_face.reshape(img_shape))
plt.subplot(1,3,2); plt.title("Spectrum of eigenvalues")
plt.bar(np.arange(len(eigenvals)), eigenvals)
plt.subplot(1,3,3); plt.title("Covariance"); plt.axis('off')
plt.imshow(eigenfaces @ eigenfaces.T)
plt.show()

# ... and display the resulting eigenfaces:
plt.figure(figsize=(12, 16))
plt.gray()
for i, eigenface in enumerate(eigenfaces):
    plt.subplot(5, 4, i + 1)
    plt.axis('off')
    plt.imshow(eigenface.reshape(img_shape))
    plt.title(f"Eigenface {i} (var={100*eigenvals[i]/eigenvals.sum():.2f}%)")
plt.show()<jupyter_output><empty_output><jupyter_text>**c)** Now project the training face images into the eigenspace to calculate their ”feature vectors”,i.e. a representation with significantly lower dimension. For the projection of the face images,they have to be centered first, i.e. the mean face vector has to be subtracted. Store the mean face in some vector (`mean_face`) and the representation achieved in some array (`face_db`). Finally restore the images from `face_db` and display them alongside the original image. Try out the effect of changing the number of eigenfaces to be used (`num_eigenfaces`).<jupyter_code># number of eigenfaces to be used
num_eigenfaces = 19

# Remark: a value of 20 (theoretical perfect reconstruction) may suffer from numerical
# instability (the last eigenface may introduce noise).
# However, a value of 19 suffices for almost perfect reconstruction ...

# BEGIN SOLUTION

# Reduce the number of eigenfaces to achieve more dimensionality reduction.
eigenfaces_used = eigenfaces[:num_eigenfaces]

# Project images into the eigenface space and store
# them in a "eigenface database":
face_db = (face_vecs - mean_face) @ eigenfaces_used.T

# ... and display the original image and the stored image
print(f"Eigenfaces used: {len(eigenfaces_used)}/{len(eigenfaces)}")
plt.figure(figsize=(12, 12))
plt.gray()
for i, face in enumerate(face_db):
    restored = (face[np.newaxis] @ eigenfaces_used) + mean_face
    plt.subplot(5, 8, 2 * i + 1)
    plt.axis('off')
    plt.imshow(train_imgs[i])
    plt.title('original')
    plt.subplot(5, 8, 2 * i + 2)
    plt.axis('off')
    plt.imshow(restored.reshape(img_shape))
    plt.title('restored')
plt.show()

assert np.allclose(face_vecs, face_db @ eigenfaces_used + mean_face, atol=1e-4)

# END SOLUTION<jupyter_output><empty_output><jupyter_text>**d)** Implement the function `recognize_face` that recognizes a face from that database by calculating the euclidean distance of this face feature vector to all of the training feature vectors from the database. The feature vector with the smallest distance represents the winner category. Check your implementation by recognizing the images from the training set (they should be recognized without error).<jupyter_code>from scipy.spatial.distance import cdist

def recognize_face(face, eigenfaces, mean_face, face_db):
    """
    Recognize a face from a face database.
    and return the index of the best matching database entry.

    The FACE is first centered and projected into the eigeface
    space provided by EIGENFACES. Then the best match is found
    according to the euclidean distance in the eigenface space.
    
    Args:
        face (ndarray): Face to be recognised.
        eigenfaces (ndarray): Array of eigenfaces.
        mean_face (ndarray): Average face.
        face_db (ndarray): Database of faces projectected into Eigenface space.
        
    Returns:
        index (uint): Position of the best matching face in face_db.
    """
    index = -1

    # BEGIN SOLUTION
    # center the face
    centered = face - mean_face

    # and project it into the eigenface space
    projected = eigenfaces @ centered

    # Now compute the similarity to all known faces
    # (comparison is performed in the eigenface space)
    distances = cdist(face_db, projected[None, :])
    index = distances.argmin()

    # END SOLUTION

    return index


# ... and now check your function on the training set ...
# BEGIN SOLUTION
def show_recognition_results(imgs, labels, train_imgs, train_labels,
                             num_eigenfaces, eigenfaces, mean_face, face_db):
    """Iterate over all face images and compute the best matching face in face_db.
    
    Args:
        imgs (list): List of test faces.
        train_imgs (list): List of training faces.
        train_labels (list): List of training labels.
        num_eigenfaces (uint): Number of eigenfaces.
        eigenfaces (list): List of the eigenfaces.
        mean_face (ndarray): Average face.
        face_db (ndarray): Database of faces projectected into Eigenface space.
        
    Returns:
    
    """
    
    img_shape = imgs[0].shape
    plt.figure(figsize=(12, 12))
    plt.suptitle(
        'Face recognition based on {} principal components'.format(num_eigenfaces))
    plt.gray()
    for j, img in enumerate(imgs):

        # find the best match in the eigenface database
        winner = recognize_face(
            img.reshape(np.prod(img_shape)), eigenfaces, mean_face, face_db)
        name_label = labels[j][5:7]
        name_winner = train_labels[winner][5:7]

        plt.subplot(5, 8, 2 * j + 1)
        plt.axis('off')
        plt.imshow(img)
        plt.title(labels[j][5:7])

        plt.subplot(5, 8, 2 * j + 2)
        plt.axis('off')
        plt.imshow(train_imgs[winner])
        t = plt.title(('*' if name_label != name_winner else '') + name_winner)
        plt.setp(t, color='r' if name_label != name_winner else 'g') 
    plt.show()
    
show_recognition_results(train_imgs, train_names,
                         train_imgs, train_names,
                         num_eigenfaces, eigenfaces_used, 
                         mean_face, face_db)
# END SOLUTION<jupyter_output><empty_output><jupyter_text>**e)** Now classify the images in directory `images/testimg/`. Try to reduce the number of principal componentsused. How many PCs are necessary to still achieve perfect classification?<jupyter_code>test_imgs, test_names = read_images_from_directory('testimg', suffix,
                                                   img_shape)

# BEGIN SOLUTION
show_recognition_results(test_imgs, test_names,
                         train_imgs, train_names,
                         num_eigenfaces, eigenfaces_used, 
                         mean_face, face_db)
#END SOLUTION<jupyter_output><empty_output><jupyter_text>Exercise 3: Local features and interest points [4 points] **a)** Explain in your own words: What are *local features*? How are they used? Local features are used in object recognition to make it more resiliant against object rotation, partial occlusion, illuminence differences, scaling etc. This is achieved by using not the whole object but the most interesting parts of it - local patches that describe the object and are most of the time available in other representations of that object (for example the wheels and lights of a car, but not it's color). **b)** What are *interest points* and what are they used for? What properties are desirable? Interest points are points in an image and can be used for object recognition. They should be salient, i.e. “special” or “rare”, either within the image, or with respect to “common” images.They should be stable, i.e. should keep positions under disruptions in an image and should remain in the same position with respect to the physical world in a different image of the same scene (e.g., change of viewpoint or illumination).  Exercise 4: Computing interest points [6 points] **a)** Explain in your own words the idea of the Moravec IP operator. What are its properties? Implement this IP operator and apply it to the image `lighthouse.png`. Try different window sizes and threshold values. The Moravec IP operator measures the saliency or "uniqueness" of a window arround a pixel. In the simplest case the window is shifted by one pixel and compared to itself. This is done in both horizontal and vertical directions. The detector response is the miminum over the four directions. A corner is detected where the response exceeds a threshold.The Moravec operator is anisotropic, meaning direction dependent, and uses a "hard" window like a box filter.<jupyter_code>%matplotlib inline
import imageio
import numpy as np
import scipy.ndimage as nd
import matplotlib.pyplot as plt

img = imageio.imread('lighthouse.png', pilmode='F')

window_size = 3
threshold = .2

def moravec(img, window_size):
    """Moravec corner detector.
    
    Arguments
    ---------
    img: np.ndarray
        The input image
    window_size: int
        The size of the window to consider

    Results
    -------
    response: np.ndarray
        Response (of the same size as img), indicating interest points.
    """
    response = np.zeros_like(img)
    # BEGIN SOLUTION
    wsh = window_size//2
    for i in np.arange(window_size, img.shape[0] - window_size):
        for j in np.arange(wsh, img.shape[1] - window_size):
            patch = img[i-wsh:i+wsh+1, j-wsh:j+wsh+1]
        
            right = np.sum((patch - img[ i - wsh : i + wsh + 1, j : j + window_size])**2)
            up = np.sum((patch - img[ i - window_size + 1: i + 1, j - wsh : j + wsh + 1])**2)
            up_right = np.sum((patch - img[ i - window_size + 1: i + 1, j : j + window_size])**2)
            down_right = np.sum((patch - img[ i : i + window_size, j : j + window_size])**2)
        
            response[i,j] = min(right, up, up_right, down_right)
    return response

# alternative implementation:
def moravec2(img, window_size):
    diff_right = (img[1:-1,1:] - img[1:-1,:-1])**2
    diff_up = (img[:-2,1:] - img[1:-1,1:])**2
    diff_up_right = (img[:-2,1:] - img[1:-1,:-1])**2
    diff_down_right = (img[2:,1:] - img[1:-1,:-1])**2

    kernel = np.ones((window_size, window_size))
    E_right = nd.convolve(diff_right, kernel)
    E_up = nd.convolve(diff_up, kernel)
    E_up_right = nd.convolve(diff_up_right, kernel)
    E_down_right = nd.convolve(diff_down_right, kernel)
    
    response = np.stack((E_right, E_up, E_up_right, E_down_right)).min(axis=0)
    response /= response.max()
    # END SOLUTION
    return response

#response = moravec(img, window_size)
response = moravec2(img, window_size)

thresholded = np.zeros_like(response)
thresholded[response > threshold] = 1

plt.figure(figsize=(12,8))
plt.subplot(2,2,1); plt.imshow(img, cmap='gray'); plt.title('Original')
plt.subplot(2,2,2); plt.imshow(response, cmap='viridis'); plt.title('Moravec "heatmap"')
plt.colorbar()
plt.subplot(2,2,3); plt.imshow(thresholded, cmap='gray'); plt.title(f'Thresholded (>{threshold})')
plt.subplot(2,2,4); plt.imshow(img[1:-1,:-1], cmap='gray'); plt.title('Corners')
mask = np.zeros(thresholded.shape + (4,), dtype=np.int)
mask[:,:,0] = 255
mask[:,:,3] = thresholded*255
plt.imshow(mask, interpolation='none')
plt.tight_layout()
plt.show()
# END SOLUTION<jupyter_output><empty_output><jupyter_text>**b)** How does the Harris corner detector work and in what sense does it improve the Moravec IP operator. Implement the Harris corner detector and apply it to `lighthouse.png`. Again, try different "window sizes" (values for $\sigma$). Harris corner detector uses a Gaussian as isotropic windowing function addressing both shortcommings, the "hard" window and the anisotropy of the Moravec detector.<jupyter_code>%matplotlib inline
import imageio
import numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage as nd

k = .04
window_size = 3

img = imageio.imread('lighthouse.png', pilmode='F')/255.

def harris(img, window_size=3, k=0.04):
    """Harris corner detector.
    
    Arguments
    ---------
    img: np.ndarray
        The input image
    window_size: int
        The size of the window to consider
    k: float
        The parameter k
        
    Results
    -------
    response: np.ndarray
        Response (of the same size as img), indicating interest points.
    """
    response = np.zeros_like(img)
    # BEGIN SOLUTION
    
    # probably there is a more efficient implementation ...
    
    wsh = window_size//2

    dx = nd.sobel(img, 1) 
    dy = nd.sobel(img, 0)
    dx2 = np.multiply(dx, dx)
    dy2 = np.multiply(dy, dy)
    dxy = np.multiply(dx, dy)

    kernel = np.zeros((window_size, window_size))
    kernel[wsh, wsh] = 1
    kernel = nd.filters.gaussian_filter(kernel, wsh/1, truncate = 1)

    for i in np.arange(wsh, img.shape[0] - wsh):
        for j in np.arange(wsh, img.shape[1] - wsh):
        
            a_ = dx[i - wsh : i + wsh + 1, j - wsh : j + wsh + 1]         
            b_ = dy[i - wsh : i + wsh + 1, j - wsh : j + wsh + 1]         
            a = np.sum(np.multiply(kernel, dx2[i - wsh : i + wsh + 1, j - wsh : j + wsh + 1]))       
            b = np.sum(np.multiply(kernel, dy2[i - wsh : i + wsh + 1, j - wsh : j + wsh + 1]))
            c = np.sum(np.multiply(kernel, dxy[i - wsh : i + wsh + 1, j - wsh : j + wsh + 1]))         

            # Noble's corner measure
            #response[i,j] = 2 * ((a * b - c ** 2) / ((a + b) ** 2 + 1e-5))
            response[i,j] = a * b - c ** 2 - k * (a + b) **2
    return response

def harris2(img, window_size=3, k=0.04):

    # compute the derivatives
    dx, dy = nd.sobel(img, 1), nd.sobel(img, 0)
    dx2, dy2, dxy = dx**2, dy**2, dx * dy

    # apply apply the Gaussian window
    a = nd.filters.gaussian_filter(dx2, (window_size-1)/6, truncate = 3)
    b = nd.filters.gaussian_filter(dy2, (window_size-1)/6, truncate = 3)
    c = nd.filters.gaussian_filter(dxy, (window_size-1)/6, truncate = 3)
    
    # and compute the response signal
    response = a * b - c ** 2 - k * (a + b) **2
    # END SOLUTION
    return response

response = harris(img, window_size, k)
corners = response.copy()        
corners[response<0] = 0

plt.figure(figsize=(12,12))
plt.subplot(2,2,1); plt.imshow(img, cmap='gray'); plt.title("Original")
plt.subplot(2,2,2); plt.imshow(response, cmap='viridis'); plt.title('Harris "heatmap"')
plt.colorbar()
plt.subplot(2,2,3); plt.imshow((corners**.25), cmap='gray'); plt.title(f"Thresholded (<0)")
plt.subplot(2,2,4); plt.imshow(img, cmap='gray'); plt.title('Corners')
mask = np.zeros(corners.shape + (4,), dtype=np.int)
mask[:,:,0] = 255
mask[:,:,3] = (response>0)*255
plt.imshow(mask, interpolation='none')
plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11030308717275031, 'fraction_numerical': 0.01534202871352163, 'mean_word_length': 3.35, 'pattern_counts': {'":': 1, '<': 39, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '411643', 'license': 'MIT', 'n_tokens_mistral': 6976, 'n_tokens_neox': 6602, 'n_words': 2283, 'path': 's11/sheet11_solution.ipynb', 'repo_name': 'tbohne/CV2020', 'size': '21345'}