write a canny edge detector from scratch in matlab

main.m
function edges = my_canny_edge_detector(image, sigma, threshold_low, threshold_high)
% Applies the Canny edge detection algorithm to a grayscale input image
% Input:
%   image: 2D grayscale input image
%   sigma: standard deviation of the Gaussian filter used for smoothing
%   threshold_low: lower threshold for hysteresis thresholding
%   threshold_high: higher threshold for hysteresis thresholding
% Output:
%   edges: binary image where detected edges are set to 1
%
% Based on the MATLAB implementation: 
% https://www.mathworks.com/help/images/ref/edge.html

% Apply gaussian filter for noise reduction
smoothed = imgaussfilt(image, sigma);

% Compute the gradients using Sobel filters
[Gx, Gy] = gradient(smoothed);

% Compute the gradient magnitude and direction
grad_magnitude = hypot(Gx, Gy);
grad_direction = atan2d(Gy, Gx);

% Compute non-maximum suppression to thin edges
grad_direction(grad_direction<0) = grad_direction(grad_direction<0) + 180;
grad_direction_2 = round(grad_direction ./ 45) * 45; % round to nearest of the four angles (0, 45, 90, 135)
grad_direction_2(grad_direction_2==180) = 0; % wrap back to 0 for 180 degrees
grad_direction_2 = grad_direction_2 + 1; % shift to 1,2,3,4 for indexing

grad_magnitude_padded = padarray(grad_magnitude, [1 1], 'replicate');
edges_thin = zeros(size(grad_magnitude_padded));
for i = 2:size(grad_magnitude_padded, 1)-1
    for j = 2:size(grad_magnitude_padded, 2)-1
        if grad_magnitude_padded(i,j) >= grad_magnitude_padded(i-1:i+1,j-1:j+1,grad_direction_2(i,j))
            edges_thin(i,j) = 1;
        end
    end
end
edges_thin = edges_thin(2:end-1, 2:end-1); % remove padding

% Perform hysteresis thresholding to remove spurious edges
edges = zeros(size(edges_thin));
edges(edges_thin >= threshold_high) = 1;
edges = bwmorph(edges, 'thin', Inf);

strong_edges = edges;
weak_edges = edges_thin > threshold_low & edges_thin < threshold_high;

while ~isempty(find(weak_edges, 1))
    neighbors = bwmorph(strong_edges, 'dilate', 1);
    new_strong_edges = neighbors & weak_edges;
    strong_edges = strong_edges | new_strong_edges;
    weak_edges = weak_edges & ~new_strong_edges; 
end

edges = strong_edges;
end
2183 chars
58 lines