function centroids = centroid_image(image, centroids)
%
% This function centroids a supplied image. It returns a centroids structure
%
% '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
%
%
% centroids the image provided
% Written by Aidan Brooks and Won Kim
% 15th May 2010

%------------------------------------------------
% check if 'centroids' as supplied is a structure, if not, make it one
if (exist('centroids') == 0)
    centroids = struct;
elseif (isstruct(centroids) == 0)
    centroids = struct;
end

% extract background data from input centroids
if (isfield(centroids, 'image_background_level') == 1)
    background = centroids.image_background_level;
else
    % default background
    background = 47;            % default background level
    centroids.image_background_level = background;
end

% extract spots radius from input centroids
if (isfield(centroids, 'spot_radius') == 1)
    radius = centroids.spot_radius;
else
    radius = 10.0;
    centroids.spot_radius = radius;
end
%------------------------------------------------


% modify the image - remove background and square the result
modimage = (image - background).^2;

% Peak detection does not appear to work unless the array is scaled.
modimage = modimage./max(max(modimage));


%------------------------------------------------
% set or get the threshold level for the spots
% extract threshold level from input centroids
if (isfield(centroids, 'spot_threshold_level') == 1)
    level = centroids.spot_threshold_level;
else
    level = graythresh(modimage)/1E4;
    centroids.spot_threshold_level = level;
end

% set or get the bwimage of regions 
% ***** note: this keeps the number of spots fixed after initialization!! ******
% Note that 'regionprops' in MATLAB vR2008b doesn't support BW image input
if (isfield(centroids, 'bwimage_of_regions') == 1)
    bwimage = centroids.bwimage_of_regions;
else
    bwimage = im2bw(modimage,level);
    %cc = bwconncomp(bwimage);
    bwimage = bwlabel(bwimage);
    centroids.bwimage_of_regions = bwimage;
end
%------------------------------------------------


% centroid the modified image in the regions identified by bwimage
ctr = regionprops(bwimage, modimage, {'centroid','weightedcentroid'});
number_of_spots = numel(ctr);
centoids.number_of_centroids_in_reference = number_of_spots;

% get the centroids out of cells and into a regular array
calculated_centroids = zeros(number_of_spots, 2);
for ii = 1:numel(ctr)
    calculated_centroids(ii, :) = ctr(ii).WeightedCentroid -[0.5 0.5];
end

% add the calculated centroids to the centroid structure
centroids.current_centroids = calculated_centroids;

%------------------------------------------------
% set or get the reference centroids of the centroids
if (isfield(centroids, 'reference_centroids') ~= 1)
    centroids.reference_centroids = calculated_centroids;
end

% set the displacement of the centroids from the reference
centroids.displacement_of_centroids = centroids.current_centroids - ...
                                      centroids.reference_centroids;