function [xSpiral,ySpiral,xRefs,yRefs,resultsIndices] = centroidFind_spiral(xPositions, yPositions, xRefsAll, yRefsAll, xSize, ySize,...
    xGridLength, yGridLength, gridIndexForRefs, diffBoxRad, maxBoxRad)

% Function takes as input the co-ordinates of found spots in a Hartmann-Shack
% image, and the co-ordinates of the reference spots. It returns reference
% spot co-ordinates for which spots were found, together with a correctly
% ordered list of co-ordinates of the found spots.
%
% ------------------
% Inputs
% ------------------
%
% xPositions       --> Co-ordinates of found spots, in the x direction
% yPositions       --> Co-ordinates of found spots, in the y direction
%
% xRefsAll         --> Co-ordinates of all reference spots, in the x direction
% yRefsAll         --> Co-ordinates of all reference spots, in the y direction
%
% xSize            --> Size of the image, in the x direction
% ySize            --> Size of the image, in the y direction
%
% xGridLength      --> The number of lenslets considered across the image, in the x direction
% yGridLength      --> The number of lenslets considered across the image, in the y direction
%
% gridIndexForRefs --> The indices of each reference spot in the square grid of dimensions xGridLength x yGridLength
%
% diff_box_rad     --> Radius of a diffraction-limited spot, in pixels
% max_box_rad      --> Radius of a lenslet, in pixels
%
% ------------------
% Outputs
% ------------------
%
% xSpiral         --> Co-ordinates of found spots, ordered so as to correspond to the reference spots given in xRefs
% ySpiral         --> Co-ordinates of found spots, ordered so as to correspond to the reference spots given in yRefs
%
% xRefs           --> Co-ordinates of reference spots for which spots were found
% yRefs           --> Co-ordinates of reference spots for which spots were found
%
% resultsIndices  --> Indices of the results relative to all of the reference spots


% Set some search parameters

maxBoxRad = 0.8 * maxBoxRad; % Restrict search to 80% of the size of a lenslet
maxBoxDiam = 2 * maxBoxRad;
firstBoxRad = maxBoxRad; % Radius of the first spot that we search for. This radius will be increased if we initally find no first spot
initialBoxIncrement = 4 * diffBoxRad; % At the beginning we will not know even roughly how far apart the spots are from each other. Search this far away initially
shortestGridSide = min([xGridLength yGridLength]);


% Initialize some arrays

xSpiralGrid = zeros(size(xRefsAll));
ySpiralGrid = zeros(size(yRefsAll));
xDummyGrid = zeros(size(xRefsAll));
yDummyGrid = zeros(size(yRefsAll));


% Find approximate position of the spot nearest the middle of the array

startIndex = yGridLength * (ceil(xGridLength/2) - 1) + ceil(yGridLength/2);
approxPos = find(gridIndexForRefs == startIndex);
xApprox = xRefsAll(approxPos); yApprox = yRefsAll(approxPos);


% Loop until we locate the first spot
notFound = 1;
while notFound
    
    % If our search box gets too big, we probably have a bad image. Abort.
    if firstBoxRad > 80
        
        beep
        fprintf('No spots close to image centre, aborting spot find!\n')
        xRefs = []; yRefs = []; xSpiral = []; ySpiral = []; resultsIndices = [];
        return
        
    end
    
    % Make a search box around the approximate centre:
    
    yExtreme1 = yApprox - firstBoxRad;
    yExtreme2 = yApprox + firstBoxRad;
    xExtreme1 = xApprox - firstBoxRad;
    xExtreme2 = xApprox + firstBoxRad;
    
    % Find the spot that lies within this search box
    
    withinSpotIndices = find((xPositions >= xExtreme1) & (xPositions <= xExtreme2) & (yPositions >= yExtreme1) & (yPositions <= yExtreme2));
    
    if length(withinSpotIndices) > 1
        
        % There are multiple spots within this area, so pick the closest one
        
        r2 = (xPositions(withinSpotIndices) - xApprox).^2 + (yPositions(withinSpotIndices) - yApprox).^2;
        withinSpotIndices = withinSpotIndices(r2==min(r2));
        notFound = 0; % exit the loop
        
    elseif length(withinSpotIndices) == 1
        
        notFound = 0; % One spot was found - exit the loop
        
    else

        % No spots were found - increase the search radius. Do not exit the loop
        firstBoxRad = firstBoxRad + diff_box_rad;
        
    end
    
end


% Record the first spot position

xPeak = xPositions(withinSpotIndices);
yPeak = yPositions(withinSpotIndices);
xSpiralGrid(startIndex) = xPeak;
ySpiralGrid(startIndex) = yPeak;
xDummyGrid(startIndex) = xPeak;
yDummyGrid(startIndex) = yPeak;


% Set some parameters that will be used in the main loop

xPrevious = xPeak; % The last spot we found was located here
yPrevious = yPeak;
lastXDif = 0; % Most recent displacement between spots
lastYDif = 0;

whereNowX = ceil(xGridLength/2); % Position in the results grid to store spot locations
whereNowY = ceil(yGridLength/2); 

plusMinusFlag = 1; % +1 is movement in the positive direction, -1 the negative direction (Cartesian co-ordinates)           

thisVectorLength = 0; % Length of the current spiral arm in units of lenslets

emptyArm = 1; % This flag signifies that the current arm of the spiral was empty (assessed to decide when to stop spiralling)
emptyArms = zeros(2,2); % Empty arms found on all four sides means that we will stop


% ------------------------
% --- Start the spiral ---
%-------------------------

cont = 1; % Keep looping as long as this flag is 1
while cont
    
    
    % If we're too close to the edge of the grid, don't make the next arm any longer than the current, and stop spiralling
    % after that arm is done:
    if whereNowY == 1 || whereNowY == shortestGridSide

        cont = 0;
        thisVectorLength = thisVectorLength - 1;
        
    end
    
    
    thisVectorLength = thisVectorLength + 1; % Make the next arm of the spiral 1 spot longer
    
    xyFlag = -1;
    
    
    % Loop in the y (0) and then the x (1) directions:
    
    while xyFlag <= 0
        
        xyFlag = xyFlag + 1;
        
        iCOUNTER = 1; % the spot number along this current armo of the spiral
        
        % Loop through all spots on this arm:
        
        while iCOUNTER <= thisVectorLength
            
            iCOUNTER = iCOUNTER + 1;
            
            if xyFlag       % There is nominally an x displacement between spots we are considering
                
                whereNowX = whereNowX + plusMinusFlag; % Update the position at which to store results in the grid
                
                % Look in the next most inner arm of the spiral to determine if there are any previously found spot
                % co-ordinates we can use to update our expected spot positions. One of these is adjacent to the spot we
                % wish to find ('nearCurrent') and the other is adjacent to the previous spot found ('nearPrevious'):
                
                nearCurrentX = xDummyGrid(whereNowY + plusMinusFlag,whereNowX);
                nearPreviousX = xDummyGrid(whereNowY + plusMinusFlag,whereNowX - plusMinusFlag);
                nearCurrentY = yDummyGrid(whereNowY + plusMinusFlag,whereNowX);
                nearPreviousY = yDummyGrid(whereNowY + plusMinusFlag,whereNowX - plusMinusFlag);
                
                
                % If there were non-zero values, use these to determine the likely difference between the previous found spot
                % and the spot we're looking for now:
                
                if nearCurrentX && nearPreviousX
                    
                    % This is the most common occurrence
                    lastXDif = nearCurrentX - nearPreviousX;
                    lastYDif = nearCurrentY - nearPreviousY;
                    
                elseif lastXDif
                    
                    % This occurs for example when we reach the end of an arm. We will use the last known difference in the
                    % x direction (i.e. do not alter lastXDif). Set lastYDif to 0 to prevent poor behaviour around the edge
                    % of the pupil.
                    
                    lastYDif = 0;
                    
                else
                    
                    % This occurs for the second displacement since there is no previous x data (only y data, from the first
                    % displacement). Use a default offset.

                    lastXDif = initialBoxIncrement;
                    lastYDif = 0;
                    
                end
                
                % From the calculated displacement, determine where we expect the next spot to be
                
                xExpected = xPrevious + plusMinusFlag*abs(lastXDif);
                yExpected = yPrevious + lastYDif;


            else       % y displacement between lenslets in the grid - same logic as above, so no comments here
                
                whereNowY = whereNowY - plusMinusFlag;
                
                nearCurrentY = yDummyGrid(whereNowY,whereNowX + plusMinusFlag);
                nearPreviousY = yDummyGrid(whereNowY + plusMinusFlag,whereNowX + plusMinusFlag);
                nearCurrentX = xDummyGrid(whereNowY,whereNowX + plusMinusFlag);
                nearPreviousX = xDummyGrid(whereNowY + plusMinusFlag,whereNowX + plusMinusFlag);
                
                if nearCurrentY && nearPreviousY
                    
                    lastYDif = nearCurrentY - nearPreviousY;
                    lastXDif = nearCurrentX - nearPreviousX;
                    
                elseif lastYDif
                    
                    lastXDif = 0;
                    
                else
                    
                    % If there was no previous data we are on the first spot - use the default box size:
                    
                    lastYDif = initialBoxIncrement;
                    lastXDif = 0;
                    
                end
                
                xExpected = xPrevious + lastXDif;
                yExpected = yPrevious - plusMinusFlag*abs(lastYDif);
                
            end
            
            % Determine whether this new expected position lies within a box width of the image edge. If it does, just use a dummy spot.
            
            if xExpected > maxBoxDiam && yExpected > maxBoxDiam && ...
                    xExpected < (xSize - maxBoxDiam) && yExpected < (ySize - maxBoxDiam);
                
                % The search area around our predicted position will lie entirely within the image. Create the search area.
    
                yExtreme1 = yExpected - maxBoxRad;
                yExtreme2 = yExpected + maxBoxRad;
                xExtreme1 = xExpected - maxBoxRad;
                xExtreme2 = xExpected + maxBoxRad;

                % Determine the spot that lies within this search area
                
                withinSpotIndices = find((xPositions >= xExtreme1) & (xPositions <= xExtreme2) & (yPositions >= yExtreme1) & (yPositions <= yExtreme2));
                
                if length(withinSpotIndices) > 1
                    
                    % If there are multiple spots, pick the closest one
                    r2 = (xPositions(withinSpotIndices) - xExpected).^2 + (yPositions(withinSpotIndices) - yExpected).^2;
                    withinSpotIndices = withinSpotIndices(r2==min(r2));
                    
                end

                if length(withinSpotIndices) == 1
                    
                    % Just one spot found - record its position
                    
                    xPeak = xPositions(withinSpotIndices);
                    yPeak = yPositions(withinSpotIndices);
                    xSpiralGrid(whereNowY,whereNowX) = xPeak;
                    ySpiralGrid(whereNowY,whereNowX) = yPeak;
                    xDummyGrid(whereNowY,whereNowX) = xPeak;
                    yDummyGrid(whereNowY,whereNowX) = yPeak;
                    xPrevious = xPeak;
                    yPrevious = yPeak;
                    emptyArm = 0; % the arm was not empty
                    
                elseif thisVectorLength == 1
                    
                    % Multiple spots or no spots were found, and this is the second spot overall.
                    % Move the search box a little further away from the first spot.
                                           
                    initialBoxIncrement = initialBoxIncrement + diff_box_rad * 2;

                    % Wind back some loop parameters as we are going to repeat this spot:

                    iCOUNTER = iCOUNTER - 1;
                    if xyFlag
                        whereNowX = whereNowX - plusMinusFlag;
                        lastYDif = 0;
                        lastXDif = initialBoxIncrement;
                        plusMinusFlag = -plusMinusFlag;
                        xyFlag = xyFlag - 1;
                    else
                        whereNowY = whereNowY + plusMinusFlag;
                        lastYDif = initialBoxIncrement;
                        lastXDif = 0;
                    end
                        
                else
                    
                    % There was no spot within the area around our prediction. This commonly occurs e.g. when navigating the pupil edge. Use a dummy spot.

                    xPrevious = xExpected;
                    yPrevious = yExpected;

                    xDummyGrid(whereNowY,whereNowX) = xPrevious;
                    yDummyGrid(whereNowY,whereNowX) = yPrevious;
                                       
                end
                
            else
                
                % Search area would be outside the image, so use a dummy (with a rectangular CCD, or decentration, a subsequent arm may still be
                % within the image)
                
                xPrevious = xExpected; 
                yPrevious = yExpected;

                xDummyGrid(whereNowY,whereNowX) = xPrevious; 
                yDummyGrid(whereNowY,whereNowX) = yPrevious;
                
            end
            
        end
        
        % We have finished looping spots on this arm of the spiral. Check
        % to see if we have encountered empty arms on all sides of our data.

        if emptyArm
            
            emptyArms(xyFlag + 2, plusMinusFlag/2 + 1.5) = 1;
            
            if length(find(emptyArms)) == 4
                
                % Empty arms surrounding our data - we are finished
                cont = 0;
                continue
                
            end
            
        end

        % Next arm is empty until proven otherwise:
        emptyArm = 1;
        
    end
    
    % We have made spiral arms in the x and y direction at this arm length.
    % In the next loop iteration we will extend the arm length.

    plusMinusFlag = -plusMinusFlag; % To maintain the spiral shape

end

% Finished looping. Results are currently in a grid format with zeroes
% where there were no spots ascribed to a reference. Extract the non-zero 
% data.

resultsIndices = find(xSpiralGrid);

xSpiral = xSpiralGrid(resultsIndices);
ySpiral = ySpiralGrid(resultsIndices);

xRefs = xRefsAll(resultsIndices);
yRefs = yRefsAll(resultsIndices);