%==============================================================================
% cannyedge.q
% Canny edge detector in Quasar
% Bart Goossens, Nov 2013
% Note that this is an experimental program (the hysteresis thresholding is
% not recursive yet) - feel free to change it and send me the updated version
%==============================================================================
import "imfilter.q"

{!author name="Bart Goossens"}
{!doc category="Image Processing/Edge detection"}

% Function: cannyedge
%
% Implementation of the Canny edge detector
%
% Usage:
%
% :function y = cannyedge(x : cube, sigma=1, low_threshold=1, high_threshold=0.1)
%
% Parameters:
% x - the input image
% sigma - parameter for the Gaussian filter (used for smoothing)
% low_threshold - hysteresis lower threshold
% high_threshold - hysteresis higher threshold
%
% Remarks:
%
% The steps involved:
%
% * Smooth using the Gaussian with sigma above.
%
% * Apply the horizontal and vertical Sobel operators to get the gradients
%   within the image. The edge strength is the sum of the magnitudes
%   of the gradients in each direction.
%
% * Find the normal to the edge at each point using the arctangent of the
%   ratio of the Y sobel over the X sobel. Pragmatically, we can
%   look at the signs of X and Y and the relative magnitude of X vs Y
%   to sort the points into 4 categories: horizontal, vertical,
%   diagonal and antidiagonal.
%
% * Look in the normal and reverse directions to see if the values
%   in either of those directions are greater than the point in question.
%   Use interpolation to get a mix of points instead of picking the one
%   that's the closest to the normal.
%
% * Label all points above the high threshold as edges.
% * Recursively label any point above the low threshold that is 8-connected
%   to a labeled point as an edge.
%
function y = cannyedge(x : cube, sigma=1, lowThreshold : scalar =5, highThreshold : scalar =10)

    smoothed = gaussian_filter(x, sigma)
    %sobel : ccube = complex(hor_sobel(x), ver_sobel(x))
    sobel = complex(hor_sobel(x), ver_sobel(x))
    theta = mod(pi/2+angle(sobel),pi)
    magnitude = abs(sobel)    
    magnitude = (255 / max(magnitude)) * magnitude
    
    % Calculate the orientation (closest angle from 0, 45, 90, 135)
    dst = zeros(size(x))
    #pragma force_parallel
    for m=0..size(x,0)-1
        for n=0..size(x,1)-1
            for k=0..size(x,2)-1
                t = theta[m,n,k] * 180/pi
                if t < 0
                    t += 180.0
                endif
                if t < 22.5  % 0°
                    left = magnitude[m, n - 1, k]
                    right = magnitude[m, n + 1, k]
                elseif t < 67.5  % 45°
                    left = magnitude[m + 1, n - 1, k]
                    right = magnitude[m - 1, n + 1, k]
                elseif t < 112.5  % 90°
                    left = magnitude[m + 1, n, k]
                    right = magnitude[m - 1, n, k]
                elseif t < 157.5  % 135°
                    left = magnitude[m + 1, n + 1, k]
                    right = magnitude[m - 1, n - 1, k]
                else  % 0°
                    left = magnitude[m, n - 1, k]
                    right = magnitude[m, n + 1, k]
                endif
                
                % Compare current pixel value with adjacent pixels
                grad = magnitude[m, n, k]
                if grad < left || grad < right
                    dst[m, n, k] = 0.0
                else
                    dst[m, n, k] = grad
                endif
            end
        end
    end   

    % Apply hysteresis thresholding - to do: this should be recursive...
    y = copy(dst)
    #pragma force_parallel    
    for m=0..size(x,0)-1
        for n=0..size(x,1)-1
            for k=0..size(x,2)-1
                if dst[m,n,k] < highThreshold
                    if dst[m,n,k] < lowThreshold
                        % non-edge
                        y[m,n,k] = 0.0
                    % check 8 neighboring pixels                        
                    elseif dst[m-1,n-1,k] < highThreshold && _
                           dst[m,n-1,k] < highThreshold && _
                           dst[m+1,n-1,k] < highThreshold && _
                           dst[m-1,n,k] < highThreshold && _
                           dst[m+1,n,k] < highThreshold && _
                           dst[m-1,n+1,k] < highThreshold && _
                           dst[m,n+1,k] < highThreshold && _
                           dst[m+1,n+1,k] < highThreshold 
                           y[m,n,k] = 0.0
                    endif
                endif
            end
        end
    end

end

function [] = main()
    im = imread("lena_big.tif")[:,:,1]
    tic()
    y = cannyedge(im)
    toc()
    imshow(y,[])

end