Here is an implementation of canny edge detection from scratch using MATLAB:
function edgeImg = myCannyDetector(img, sigma, thLow, thHigh)
% Compute Gaussian filter
filterSize = 2*ceil(3*sigma)+1;
gaussFilter = fspecial('gaussian', filterSize, sigma);
% Apply Gaussian filter
imgFiltered = imfilter(img, gaussFilter, 'replicate');
% Compute gradients using Sobel operators
filterSobelX = [1 0 -1; 2 0 -2; 1 0 -1];
filterSobelY = [1 2 1; 0 0 0; -1 -2 -1];
imgSobelX = imfilter(imgFiltered, filterSobelX, 'replicate');
imgSobelY = imfilter(imgFiltered, filterSobelY, 'replicate');
gradMag = sqrt(imgSobelX.^2 + imgSobelY.^2);
gradDir = atan2(imgSobelY, imgSobelX) * 180/pi;
% Compute non-maximum suppression
[h, w] = size(img);
edgeImg = zeros(h, w);
for i = 2:h-1
for j = 2:w-1
theta = gradDir(i, j);
if (theta < 0)
theta = theta + 180;
end
if ((theta >= 0 && theta < 22.5 ) || ...
(theta >= 157.5 && theta <= 180) || ...
(theta >= 67.5 && theta < 112.5))
% Check left and right neighbors
if (gradMag(i, j) > gradMag(i, j-1) && ...
gradMag(i, j) > gradMag(i, j+1))
edgeImg(i, j) = gradMag(i, j);
end
elseif ((theta >= 22.5 && theta < 67.5) || ...
(theta >= 112.5 && theta < 157.5))
% Check top left-right and bottom left-right neighbors
if (gradMag(i, j) > gradMag(i-1, j-1) && ...
gradMag(i, j) > gradMag(i+1, j+1)) || ...
(gradMag(i, j) > gradMag(i-1, j+1) && ...
gradMag(i, j) > gradMag(i+1, j-1))
edgeImg(i, j) = gradMag(i, j);
end
end
end
end
% Apply hysteresis thresholding
thLow = thLow * max(max(edgeImg));
thHigh = thHigh * max(max(edgeImg));
thresImg = zeros(size(edgeImg));
thresImg(edgeImg >= thHigh) = 1;
thresImg = bwmorph(thresImg, 'bridge');
while (1)
oldThresImg = thresImg;
weakEdges = edgeImg >= thLow & ~thresImg;
if any(weakEdges(:))
[r, c] = find(weakEdges);
for i = 1:numel(r)
if any(thresImg(max(1, r(i)-1):min(end, r(i)+1), ...
max(1, c(i)-1):min(end, c(i)+1)))
thresImg(r(i), c(i)) = 1;
end
end
end
if (isequal(oldThresImg, thresImg))
break;
end
end
edgeImg = thresImg;
The myCannyDetector function takes four arguments: the input image (img), the standard deviation of the Gaussian filter (sigma), and the low and high threshold values for hysteresis thresholding (thLow and thHigh). The function returns a binary image with 1's at the edges and 0's elsewhere.
Here's an example of how to use myCannyDetector to calculate homography between two images:
% Load two images
img1 = imread('image1.jpg');
img2 = imread('image2.jpg');
% Detect edges using canny edge detection
edgeImg1 = myCannyDetector(img1, 1, 0.1, 0.3);
edgeImg2 = myCannyDetector(img2, 1, 0.1, 0.3);
% Find matching features using SURF
points1 = detectSURFFeatures(edgeImg1);
points2 = detectSURFFeatures(edgeImg2);
features1 = extractFeatures(edgeImg1, points1);
features2 = extractFeatures(edgeImg2, points2);
indexPairs = matchFeatures(features1, features2);
% Calculate homography using RANSAC
tolerance = 5;
[homography, inlierIndex] = estimateGeometricTransform(...
points1(indexPairs(:,1)), points2(indexPairs(:,2)), ...
'projective', 'MaxDistance', tolerance);
% Visualize the results
figure;
showMatchedFeatures(img1, img2, ...
points1(indexPairs(:,1)), points2(indexPairs(:,2)), ...
'montage', 'Parent', axes);
title('Matched Features');
legend('Matched points 1', 'Matched points 2');
figure;
showMatchedFeatures(img1, img2, ...
points1(inlierIndex), points2(inlierIndex), ...
'montage', 'Parent', axes);
title('Inlier Features');
legend('Inlier points 1', 'Inlier points 2');
figure;
plot(img1);
hold on;
projectedPoints = transformPointsForward(homography, [1 1; 1 size(img2,1); size(img2,2) size(img2,1); size(img2,2) 1]);
line(projectedPoints(:,1), projectedPoints(:,2), 'color', 'r');
title('Projected Image 2 onto Image 1');
This code first detects edges in both images using myCannyDetector. Then, it uses the Speeded Up Robust Features (SURF) algorithm to find matching feature points between the two images. The estimateGeometricTransform function is used to calculate a homography transformation between the two sets of feature points using a Random Sample Consensus (RANSAC) algorithm. Finally, the code visualizes the matched and inlier feature points, as well as the projected image 2 onto image 1 using the calculated homography.