function [ swtMap ] = swt( im, searchDirection )

%swt Preforms stoke width transform on input image

%   A novel image operator that seeks to find the value of stroke width

%   for each image pixel.  It's use is meant for the task of text

%   detection in natural images.

%

%   im = RGB input image of size m x n x

%   searchDirection = gradient direction is either  to detect dark text on light

%   background or - to detect light text on dark background.

%

%   swtMap = resulting mapping of stroke withs for image pixels

% Convert image to gray scale

im = im2double(rgb2gray(im));

%figure, imshow(im), title('Black and White Image');

% Find edges using canny edge dector

edgeMap = edge(im, 'canny');

%figure, imshow(edgeMap), title('Edges Using Canny');

% Get all edge pixel postitions

[edgePointRows, edgePointCols] = find(edgeMap);

% Find gradient horizontal and vertical gradient

sobelMask = fspecial('sobel');

dx = imfilter(im,sobelMask);

dy = imfilter(im,sobelMask');

%figure, imshow(dx, []), title('Horizontal Gradient Image');

%figure, imshow(dy, []), title('Vertical Gradient Image');

% Initializing matrix of gradient direction

theta = zeros(size(edgeMap,),size(edgeMap,));

% Calculating theta, gradient direction, for each pixel on the image.

% ***This can be optimized by using edgePointCols and edgePointRows

% instead.***

for i=:size(edgeMap,)

    for j=:size(edgeMap,)

        if edgeMap(i,j) ==

            theta(i,j) = atan2(dy(i,j),dx(i,j));

        end

    end

end

% Getting size of the image

[m,n] = size(edgeMap);

% Initializing Stoke Width array with infinity

swtMap = zeros(m,n);

for i=:m

    for j=:n

        swtMap(i,j) = inf;

    end

end

% Set the maximum stroke width, this number is variable for now but must be

% made to be more dynamic in the future

maxStrokeWidth = ;

% Initialize container for all stoke points found

strokePointsX = zeros(size(edgePointCols));

strokePointsY = zeros(size(strokePointsX));

sizeOfStrokePoints = ;

% Iterate through all edge points and compute stoke widths

for i=:size(edgePointRows)

    step = ;

    initialX = edgePointRows(i);

    initialY = edgePointCols(i);

    isStroke = ;

    initialTheta = theta(initialX,initialY);

    sizeOfRay = ;

    pointOfRayX = zeros(maxStrokeWidth,);

    pointOfRayY = zeros(maxStrokeWidth,);

    % Record first point of the ray

    pointOfRayX(sizeOfRay+) = initialX;

    pointOfRayY(sizeOfRay+) = initialY;

    % Increase the size of the ray

    sizeOfRay = sizeOfRay + ;

    % Follow the ray

    while step < maxStrokeWidth

        nextX = round(initialX + cos(initialTheta) * searchDirection * step);

        nextY = round(initialY + sin(initialTheta) * searchDirection * step);

        step = step + ;

        % Break loop if out of bounds.  For some reason this is really

        % slow.

        if nextX <  | nextY <  | nextX > m | nextY > n

            break

        end

        % Record next point of the ray

        pointOfRayX(sizeOfRay+) = nextX;

        pointOfRayY(sizeOfRay+) = nextY;

        % Increase size of the ray

        sizeOfRay = sizeOfRay + ;

        % Another edge pixel has been found

        if edgeMap(nextX,nextY)

            oppositeTheta = theta(nextX,nextY);

            % Gradient direction roughtly opposite

            if abs(abs(initialTheta - oppositeTheta) - pi) < pi/

                isStroke = ;

                strokePointsX(sizeOfStrokePoints+) = initialX;

                strokePointsY(sizeOfStrokePoints+) = initialY;

                sizeOfStrokePoints = sizeOfStrokePoints + ;

            end

            break

        end

    end

    % Edge pixel is part of stroke

    if isStroke

        % Calculate stoke width

        strokeWidth = sqrt((nextX - initialX)^ + (nextY - initialY)^);

        % Iterate all ray points and populate with the minimum stroke width

        for j=:sizeOfRay

            swtMap(pointOfRayX(j),pointOfRayY(j)) = min(swtMap(pointOfRayX(j),pointOfRayY(j)),strokeWidth);

        end

    end

end

%figure, imshow(swtMap, []), title('Stroke Width Transform: First Pass');

% Iterate through all stoke points for a refinement pass.  Refer to figure

% 4b in the paper.

for i=:sizeOfStrokePoints

    step = ;

    initialX = strokePointsX(i);

    initialY = strokePointsY(i);

    initialTheta = theta(initialX,initialY);

    sizeOfRay = ;

    pointOfRayX = zeros(maxStrokeWidth,);

    pointOfRayY = zeros(maxStrokeWidth,);

    swtValues = zeros(maxStrokeWidth,);

    sizeOfSWTValues = ;

    % Record first point of the ray

    pointOfRayX(sizeOfRay+) = initialX;

    pointOfRayY(sizeOfRay+) = initialY;

    % Increase the size of the ray

    sizeOfRay = sizeOfRay + ;

    % Record the swt value of first stoke point

    swtValues(sizeOfSWTValues+) = swtMap(initialX,initialY);

    sizeOfSWTValues = sizeOfSWTValues + ;

    % Follow the ray

    while step < maxStrokeWidth

        nextX = round(initialX + cos(initialTheta) * searchDirection * step);

        nextY = round(initialY + sin(initialTheta) * searchDirection * step);

        step = step + ;

        % Record next point of the ray

        pointOfRayX(sizeOfRay+) = nextX;

        pointOfRayY(sizeOfRay+) = nextY;

        % Increase size of the ray

        sizeOfRay = sizeOfRay + ;

        % Record the swt value of next stoke point

        swtValues(sizeOfSWTValues+) = swtMap(nextX,nextY);

        sizeOfSWTValues = sizeOfSWTValues + ;

        % Another edge pixel has been found

        if edgeMap(nextX,nextY)

            break

        end

    end

    % Calculate stoke width as the median value of all swtValues seen.

    strokeWidth = median(swtValues(:sizeOfSWTValues));

    % Iterate all ray points and populate with the minimum stroke width

    for j=:sizeOfRay

        swtMap(pointOfRayX(j),pointOfRayY(j)) = min(swtMap(pointOfRayX(j),pointOfRayY(j)),strokeWidth);

    end

end

%figure, imshow(swtMap, []), title('Stroke Width Transform: Second Pass');

end