function [ HFWHM, HeSquared, HFitStats, VFWHM, VeSquared, VFitStats ] = SpeckleSize( SpeckleRGB )
%By: Joel Sleppy
%jdsleppy [AT] wustl.edu
%July 23, 2009
%SpeckleSize determines the speckle size from an image using
%autocovariance.
%SpeckleSize takes as an input a H by W by 3 matrix of integer values
%between 0 and 255, where H is the height of the image, W is the width of
%the image, and the three values are the RGB values for each pixel. Its
%outputs are the FWHM and 1/e^2 (beam waist) values of Gaussian fits to two arrays
%generated by the sums of the normalized autocovariances of 1) all rows in
%the image and 2) all columns in the image. From the FWHM size or the 1/e^2
%size, these arrays tell you the horizontal and vertical speckle sizes.
%Also, it outputs goodness of fit statistics (such as R-squared).
%HFWHM is the horizontal speckle size per FWHM method, VeSquared is the
%vertical speckle size per 1/e^2 method, etc.
%The easiest way to obtain an RGB matrix is to use the image processing
%toolbox to open an image from file and use 'export to workspace.'
%This is adapted from Haibo Lin's University of Missouri Ph.D thesis,
%"Speckle Mechanism in Holographic Optical Coherence Imaging"
%NOTE: The function works best when the entire image is full of speckle
%with no long-range background (such as a Gaussian envelope).
%If this is the case, H and V will have a needle-like peak in the center
%with little fluctuation elsewhere. If the needle-like peak sits on top of
%a 'hill,' either limit your FWHM calculations to the needle-like peak or crop the
%image so it is full of uniform speckle (no extra black background). The 'hills'
%are caused by the function processing the Gaussian envelope (or 'speckle within speckle') as one
%large spot against the black background. See GaussianSubtractSpeckleSize
%if needed.
%NOTE: The fit is not very good if the adjusted R-squared value is not near
%one! Check this after every run.
M=size(SpeckleRGB,1); %height of area
N=size(SpeckleRGB,2); %width of area
B=rgb2gray(SpeckleRGB); %produces grayscale image (intensity map)
B=double(B); %typecasts the values as doubles
E=B(1,:); %creates an array with the values of row 1
s=size(xcov(E)); %finds the size of the autocovariance array
D=zeros(s); %creates an empty array of size s
D=double(D); %typecasts the values as doubles
for i=1:M;
C=B(i,:);
D=imadd(D,xcov(C,'coeff')); %sums the xcov arrays for all rows into D
end;
H=(D/max(D)); %the finished horizontal product, H, is normalized
E1=B(:,1); %creates an array with the values of column 1
s1=size(xcov(E1)); %finds the size of the autocovariance array
D1=zeros(s1); %creates an empty array of size s1
D1=double(D1); %typecasts the values as doubles
for j=1:N
C1=B(:,j);
D1=imadd(D1,xcov(C1,'coeff')); %sums the xcov arrays for all rows into D1
end;
V=(D1/max(D1)); %the finished vertical product, V, is normalized
%H and V are fit to Gaussians and the speckle size is extracted from the
%fits.
helper1 = 1:size(H,2);
helper2 = 1:size(V);
gauss1 = fittype('gauss1'); %sets up the Gaussian curve fitting
excludeLowH = excludedata(helper1',H','range',[.2,1]); %excludes the noise from outside the speckle
excludeLowV = excludedata(helper2',V,'range',[.2,1]);
optionsH = fitoptions(gauss1);
optionsV = fitoptions(gauss1);
optionsH.Exclude = excludeLowH;
optionsV.Exclude = excludeLowV;
[HFit, HFitStats] = fit(helper1',H',gauss1, optionsH); %performs the fits
[VFit, VFitStats] = fit(helper2',V,gauss1, optionsV);
%The SpeckleSize code technically does not find the FWHM and 1/e2 widths of
%the Gaussian fits, but rather the widths at which the Gaussians fall to
%values of .5 and 1/e^2. This provides a better idea of the speckles size
%even when the Gaussians amplitude is not unity, as expected for a perfect
%fit of the normalized data, but can produce unexpected trends...
HFWHM = (2*(HFit.c1)*sqrt(-log(.5/(HFit.a1)))); %FWHM values (Full width when the fit = .5)
VFWHM = (2*(VFit.c1)*sqrt(-log(.5/(VFit.a1))));
HeSquared = ((HFit.c1)*sqrt(-log((.1353353)/(HFit.a1)))); %1/e^2 values (Full width when the fit = .135...)
VeSquared = ((VFit.c1)*sqrt(-log((.1353353)/(VFit.a1))));
end
