function defocus = get_defocus(centroids)
% 
% a function to extract the defocus of the gradient field
%
% 'centroids' structure
% --------------------
% centroids.image_background_level  - the background intensity of an image
%                                     with no illumination on it
%          .spot_radius             - the radius of a hartmann spot
%          .spot_threshold_level    - the minimum intensity of pixels used to
%                                     include in a centroid calculation
%          .bwimage_of_regions      - a b/w image calculated ONCE that
%                                     determines the spot regions to centroid
%          .number_of_centroids_in_reference     - the number of centroids
%          .current_centroids       - the centroids calculated from this
%                                     image
%          .reference_centroids     - the centroids calculted from the
%                                     initial image
%          .displacement_of_centroids - the displacement of these centroids
%                                       from the reference centroids
%
% 'defocus' sutrcture
% --------------------
% defocus.defocus_X           - the quadratic wavefront component of the x direction
% defocus.defocus_Y           - the quadratic wavefront component of the y direction
%
% written by Aidan Brooks. 15th May 2010
%

defocus.defocus_X = 0.0;
defocus.defocus_Y = 0.0;

% check if 'centroids' exists
if (exist('centroids')==1)
    
   % check if 'centroids' is a structure
   if (isstruct(centroids) ==1)
       
       % check if the structure contains the field with the number of
       % centroids
       if (isfield(centroids, 'number_of_centroids_in_reference')==1)
           
           % check that the number of centroids is greater than 3
           if (centroids.number_of_centroids_in_reference > 3)
                % get the coordinates
                xcoords = centroids.reference_centroids(:, 1);
                ycoords = centroids.reference_centroids(:, 2);

                % get the displacements
                dx = centroids.displacement_of_centroids(:, 1);
                dy = centroids.displacement_of_centroids(:, 2);

                % fit a 1st order polynomial to [xcoords, dx] and [ycoords, dy]
                ax = polyfit(xcoords, dx, 1);
                ay = polyfit(ycoords, dy, 1);

                defocus.defocus_X = ax(1);
                defocus.defocus_Y = ay(1);
           end
       end
   end
end