Speckle Size via Autocorrelation



16 Aug 2009 (Updated )

Outputs the size of speckles in an image that has a **uniform background**.

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

%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;
D=imadd(D,xcov(C,'coeff'));         %sums the xcov arrays for all rows into D
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
D1=imadd(D1,xcov(C1,'coeff'));      %sums the xcov arrays for all rows into D1
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

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))));

Contact us