SpeckleSize.m

function [ HFWHM, HeSquared, HFitStats, VFWHM, VeSquared, VFitStats ] = SpeckleSize( SpeckleRGB )


%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.



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 Gaussian’s 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