added vision code, exercises and reports

Author: lucala

README.md

@@ -1,2 +1,5 @@
-# vision
+# Computer Vision 2016
 Code and reports for the Computer Vision course of 2016 at ETH Zurich
+
+## Disclaimer
+Please be advised that I don't guarantee completeness or correctness of anything found here.

ex01/Exercise1.pdf

@@ -0,0 +1,73 @@
+% Exercise 3
+%
+close all;
+
+%
+% COMMENT/UNCOMMENT following lines to select image
+%
+
+%default images
+ IMG_NAME1 = 'images/I1.jpg';
+ IMG_NAME2 = 'images/I2.jpg';
+
+%own images1
+% IMG_NAME1 = 'images/cv1.jpg';
+% IMG_NAME2 = 'images/cv2.jpg';
+
+%own images2
+% IMG_NAME1 = 'images/cv1.jpg';
+% IMG_NAME2 = 'images/cv2.jpg';
+
+
+%threshold dependent on image
+%I1/I2=0.02,SW1/SW2=0.3, CV1/CV2=0.45
+switch IMG_NAME1
+    case 'images/I1.jpg'
+        thresh = 0.02;
+    case 'images/sw1.jpg'
+        thresh = 0.3;
+    case 'images/cv1.jpg'
+        thresh = 0.45;
+end
+    
+% read in image
+img1 = im2double(imread(IMG_NAME1));
+img2 = im2double(imread(IMG_NAME2));
+
+img1 = imresize(img1, 0.5);
+img2 = imresize(img2, 0.5);
+
+% convert to gray image
+imgBW1 = rgb2gray(img1);
+imgBW2 = rgb2gray(img2);
+
+
+% Task 3.1 - extract Harris corners
+ [corners1, H1, allcorners1] = extractHarrisCorner(imgBW1');
+ [corners2, H2, allcorners2] = extractHarrisCorner(imgBW2');
+
+% show images with Harris corners
+showImageWithCorners(img1, allcorners1, 8);
+showImageWithCorners(img2, allcorners2, 9); 
+ 
+% show images with Harris corners using non maximum suppression
+showImageWithCorners(img1, corners1, 10);
+showImageWithCorners(img2, corners2, 11);
+
+% Task 3.2 - extract your own descriptors
+ descr1 = extractDescriptor(corners1, imgBW1');
+ descr2 = extractDescriptor(corners2, imgBW2');
+
+% Task 3.3 - match the descriptors
+ matches = matchDescriptors(descr1, descr2, thresh);
+ 
+% In case match descriptors does not find any matches with the given
+% threshold just plot the features
+ if size(matches) == 0
+    showNoFeatureMatches(img1, corners1, img2, corners2, 15);
+ else
+    showFeatureMatches(img1, corners1(:, matches(1,:)), img2, corners2(:, matches(2,:)), 20);
+ end
+
+% Task 1.4 SIFT features
+ sift(img1, img2, imgBW1, imgBW2, 0.02, 7, 40);

ex01/code/exercise2.m

@@ -0,0 +1,25 @@
+% extract descriptor
+%
+% Input:
+%   keyPoints     - detected keypoints in a 2 x n matrix holding the key
+%                   point coordinates
+%   img           - the gray scale image
+%   
+% Output:
+%   descr         - w x n matrix, stores for each keypoint a
+%                   descriptor. m is the size of the image patch,
+%                   represented as vector
+function descr = extractDescriptor(corners, img)  
+ [~,n] = size(corners);
+ descr = zeros(9*9,n);
+
+ %descriptor of 9*9 patch organized [1 2 3; 4 5 6; 7 8 9] into column
+ %vector
+ for i=1:n
+     pos = corners(:,i);
+     w = img(pos(1)+(-4:4),pos(2)+(-4:4));
+     descr(:,i) = w(:);
+ end
+
+ 
+end

ex01/code/extractDescriptor.m

@@ -0,0 +1,90 @@
+% extract harris corner
+%
+% Input:
+%   img           - n x m gray scale image
+%   thresh        - scalar value to threshold corner strength
+%   
+% Output:
+%   corners       - 2 x k matrix storing the keypoint coordinates
+%   H             - n x m gray scale image storing the corner strength
+function [corners, H, allcorners] = extractHarrisCorner(img)
+ [n,m] = size(img);
+
+ %blur to get rid of noise
+ blur = fspecial('gaussian');
+ img = imfilter(img,blur);
+ 
+ %calculate gradients
+ gx = zeros(n-1, m-1);
+ gy = zeros(n-1, m-1);
+ 
+ for i=2:m-1
+    gx(1:(n-1),i) = (img(1:(n-1),i+1) - img(1:(n-1),i-1))/2;
+ end
+ for i=2:n-1
+    gy(i,1:(m-1)) = (img(i+1,1:(m-1)) - img(i-1,1:(m-1)))/2;
+ end
+
+ %calculate for every pixel 2x2 matrix [I_x^2 I_x*I_y; I_x*I_y I_y^2]
+ Ix2 = zeros(n-1, m-1);
+ Ixy = zeros(n-1, m-1);
+ Iy2 = zeros(n-1, m-1);
+ 
+ Ix2 = gx.^2;
+ Ixy = gx.*gy;
+ Iy2 = gy.^2;
+
+
+% calculate harris response in 3x3 patch around every pixel, the sum
+% represents Harris matrix. K is the Harris response for every pixel.
+ K = zeros(n,m);
+ for i=2:n-2
+     for j=2:m-2
+        I11 = [Ix2(i-1,j-1) Ixy(i-1,j-1); Ixy(i-1,j-1) Iy2(i-1,j-1)];
+        I12 = [Ix2(i-1,j) Ixy(i-1,j); Ixy(i-1,j) Iy2(i-1,j)];
+        I13 = [Ix2(i-1,j+1) Ixy(i-1,j+1); Ixy(i-1,j+1) Iy2(i-1,j+1)];
+        I21 = [Ix2(i,j-1) Ixy(i,j-1); Ixy(i,j-1) Iy2(i,j-1)];
+        I22 = [Ix2(i,j) Ixy(i,j); Ixy(i,j) Iy2(i,j)];
+        I23 = [Ix2(i,j+1) Ixy(i,j+1); Ixy(i,j+1) Iy2(i,j+1)];
+        I31 = [Ix2(i+1,j-1) Ixy(i+1,j-1); Ixy(i+1,j-1) Iy2(i+1,j-1)];
+        I32 = [Ix2(i+1,j) Ixy(i+1,j); Ixy(i+1,j) Iy2(i+1,j)];
+        I33 = [Ix2(i+1,j+1) Ixy(i+1,j+1); Ixy(i+1,j+1) Iy2(i+1,j+1)];
+        
+        H = I11 + I12 + I13 + I21 + I22 + I23 + I31 + I32 + I33;
+
+        K(i,j) = det(H)/trace(H);
+
+     end
+ end
+ 
+ % set threshold at 99 percent quantile, the remaining values are the 1
+ % percent highest responses
+ thresh = quantile(K(:),0.99);
+ check = true;
+ corners = [];
+ allcorners = [];
+
+ %only take features which lie more than 4 pixels in the image as to not
+ %conflict with 9x9 patch of descriptors. Inner loop computes non-maximum
+ %suppression.
+ for i=4:n-3
+     for j=4:m-3
+         if K(i,j) > thresh
+             allcorners = [allcorners [i;j]];
+             for u=-1:1
+                 for v=-1:1
+                    if K(i,j) < K(i+u,j+v)
+                        check = false;
+                    end
+                 end
+             end
+             if(check == true)
+                 corners = [corners [i;j]];
+             end
+             check = true;
+         end
+     end
+ end
+
+ H = K;
+end

ex01/code/extractHarrisCorner.m

@@ -0,0 +1,29 @@
+% match descriptors
+%
+% Input:
+%   descr1        - k x n descriptor of first image
+%   descr2        - k x m descriptor of second image
+%   thresh        - scalar value to threshold the matches
+%   
+% Output:
+%   matches       - 2 x w matrix storing the indices of the matching
+%                   descriptors
+function matches = matchDescriptors(descr1, descr2, thresh)
+    %k = 9*9 height of descriptor vector
+    [~,n] = size(descr1);
+    [~,m] = size(descr2);
+    matches = [];
+    %check every combination of descr1 and descr2 for best matches
+    for i=1:n
+        d1 = descr1(:,i);
+        for j=1:m
+            d2 = descr2(:,j);
+            d = d1 - d2;
+            ssd = dot(d,d);
+            if ssd < thresh %only add mapping if ssd below threshold
+                matches = [matches [i;j]];
+            end
+        end
+    end
+    
+end

ex01/code/images/I1.jpg

@@ -0,0 +1,19 @@
+% show feature matches between two images
+%
+% Input:
+%   img1        - n x m color image 
+%   corner1     - 2 x k matrix, holding keypoint coordinates of first image
+%   img2        - n x m color image 
+%   corner1     - 2 x k matrix, holding keypoint coordinates of second image
+%   fig         - figure id
+function showFeatureMatches(img1, corner1, img2, corner2, fig)
+    [sx, sy, sz] = size(img1);
+    img = [img1, img2];
+    
+    corner2 = corner2 + repmat([sy, 0]', [1, size(corner2, 2)]);
+    
+    figure(fig), imshow(img, []);    
+    hold on, plot(corner1(1,:), corner1(2,:), '+r');
+    hold on, plot(corner2(1,:), corner2(2,:), '+r');    
+    hold on, plot([corner1(1,:); corner2(1,:)], [corner1(2,:); corner2(2,:)], 'b');    
+end

ex01/code/images/I2.jpg

@@ -0,0 +1,13 @@
+% show image with key points
+%
+% Input:
+%   img        - n x m color image 
+%   corner     - 2 x k matrix, holding keypoint coordinates of first image
+%   fig        - figure id
+function showImageWithCorners(img, corners, fig)
+    figure(fig);
+    imshow(img, []);
+    
+    hold on, plot(corners(1,:), corners(2,:), '+r');
+
+end

ex01/code/images/cv1.jpg

@@ -0,0 +1,13 @@
+%this function is called if the matchDescriptor function was not able to
+%match any features given the threshold. It will plot the two images with
+%their respective features.
+function showNoFeatureMatches(img1, corner1, img2, corner2, fig)
+    [sx, sy, sz] = size(img1);
+    img = [img1, img2];
+    
+    corner2 = corner2 + repmat([sy, 0]', [1, size(corner2, 2)]);
+    
+    figure(fig), imshow(img, []);    
+    hold on, plot(corner1(1,:), corner1(2,:), '+r');
+    hold on, plot(corner2(1,:), corner2(2,:), '+r');    
+end

ex01/code/images/cv2.jpg

@@ -0,0 +1,32 @@
+%called with 2 grayscale images, computes sift features and descriptors
+%use threshold to filter for uniqueness, same as for Harris corners
+
+function [matches] = sift(img1, img2, imgBW1, imgBW2, threshFeatures, threshMatch, fig)
+
+%compute sift keypoints and descriptors
+[k1,d1] = vl_sift(single(imgBW1),'PeakThresh',threshFeatures);
+[k2,d2] = vl_sift(single(imgBW2),'PeakThresh',threshFeatures);
+
+%find matches between the 2 images
+[matches, scores] = vl_ubcmatch(d1, d2, threshMatch);
+
+%plot results of sift
+[sx, sy, sz] = size(img1);
+img = [img1, img2];
+    
+%only need x,y coordinates of keypoints
+k1 = k1(1:2, :);
+k2 = k2(1:2, :);
+k1 = k1(:, matches(1,:));
+k2 = k2(:, matches(2,:));
+
+k2 = k2 + repmat([sy, 0]', [1, size(k2, 2)]);
+figure(fig), imshow(img, []);    
+hold on, plot(k1(1,:), k1(2,:), '+r');
+hold on, plot(k2(1,:), k2(2,:), '+r');   
+    
+figure(fig+1), imshow(img, []);    
+hold on, plot(k1(1,:), k1(2,:), '+r');
+hold on, plot(k2(1,:), k2(2,:), '+r');    
+hold on, plot([k1(1,:); k2(1,:)], [k1(2,:); k2(2,:)], 'b');    
+end

ex01/code/images/sw1.jpg

@@ -0,0 +1,14 @@
+function [ K, R, C ] = decompose(P)
+%decompose P into K, R and t
+M = P(:,1:3);
+[Rinv,Kinv] = qr(inv(M));
+R = inv(Rinv);
+K = inv(Kinv);
+
+R = R*K(end,end);
+K = K/K(end,end);
+
+[~,~,V] = svd(P);
+C = V(:,end);
+C = C(1:3)./C(4);
+end

ex01/code/images/sw2.jpg

@@ -0,0 +1,17 @@
+function [P] = dlt(xy, XYZ)
+%computes DLT, xy and XYZ should be normalized before calling this function
+
+[~, n] = size(xy);
+A = [];
+z = zeros(1,4);
+for i=1:n
+    tmp = XYZ(:,i)';
+    top = -xy(1,i)*tmp;
+    bot = xy(2,i)*tmp;
+    %A = [A; [XYZ(:,i)' zeros(1,4) -xy(1,i)*XYZ(:,i)'; zeros(1,4) -XYZ(:,i) xy(2,i)*XYZ(:,i)]];
+    A = [A; [tmp z top; z -tmp bot]];
+end
+
+[~,~,V] = svd(A);
+P = V(:,end);
+P = reshape(P,4,3)';

ex01/code/matchDescriptors.m

@@ -0,0 +1,31 @@
+% Exervice 2
+%
+close all;
+
+%IMG_NAME = 'images/image001.jpg';
+IMG_NAME = 'img/ex2_1.jpg';
+
+%This function displays the calibration image and allows the user to click
+%in the image to get the input points. Left click on the chessboard corners
+%and type the 3D coordinates of the clicked points in to the input box that
+%appears after the click. You can also zoom in to the image to get more
+%precise coordinates. To finish use the right mouse button for the last
+%point.
+%You don't have to do this all the time, just store the resulting xy and
+%XYZ matrices and use them as input for your algorithms.
+[xy XYZ] = getpoints(IMG_NAME);
+
+xy = cart2hom(xy')';
+XYZ = cart2hom(XYZ')';
+
+% === Task 2 DLT algorithm ===
+
+[K, R, t, error] = runDLT(xy, XYZ);
+
+% === Task 3 Gold Standard algorithm ===
+
+[K, R, t, error] = runGoldStandard(xy, XYZ);
+
+% === Bonus: Gold Standard algorithm with radial distortion estimation ===
+
+%[K, R, t, error] = runGoldStandardRadial(xy, XYZ);