MATLAB Code for Cloud Size – Gaussian Fit

% Imaging parameters
Magn=1; %magnification of the imaging system
px_sizew=3.75*10^(-3); %width of CCD pixel-in mm
px_sizeh=3.75*10^(-3); %height of CCD pixel-in mm

%% Loading files and average%

%import background

bck1=imread(‘FileFolder\bck1.jpg’);

BckAvrg=double(bck1);

%import array
%C1=dlmread(‘FileFolder\1ms-10to18-10steps-50mscomp-control5averaged’,’\t’);

%import image

C1=imread(‘FileFolder\40ms.jpg’);
CAvrg=double(C1);

%% Dimension of the picture

[px_height,px_width]=size(CAvrg);
x=1:px_width;
y=1:px_height;

%% reduce data to one dimension%
windowSize = 10;
Sx=filter(ones(1,windowSize)/windowSize,1,sum(CAvrg));
%fit the overall sum
lnSx=log(Sx);

p=polyfit(x(700:850),lnSx(700:850),2);
A2=p(1);
A1=p(2);
A0=p(3);
sigmaSx=sqrt(-1/(2*A2));
muSx=A1*sigmaSx^2;
ASx=exp(A0+muSx^2/(2*sigmaSx^2));
subplot(2,2,1),
plot(x,Sx,x,ASx*exp(-(x-muSx).^2/(2*sigmaSx.^2)));

%Sy
Sy=filter(ones(1,windowSize)/windowSize,1,sum(CAvrg’));
lnSy=log(Sy);

p=polyfit(y(500:600),lnSy(500:600),2);
A2=p(1);
A1=p(2);
A0=p(3);
sigmaSy=sqrt(-1/(2*A2));
muSy=A1*sigmaSy^2;
ASy=exp(A0+muSy^2/(2*sigmaSy^2));
subplot(2,2,2),
plot(y,Sy,y,ASy*exp(-(y-muSy).^2/(2*sigmaSy.^2)));

%% find center – the highest peak from the sum%
[peakx,Cenx]=max(ASx*exp(-(x-muSx).^2/(2*sigmaSx.^2)));
[peaky,Ceny]=max(ASy*exp(-(y-muSy).^2/(2*sigmaSy.^2)));

Cenx=Cenx
Ceny=Ceny

%Single out the row/column to represent each dimension and normalized%

%% ### finding the Gaussian fit for X ###

Ex=0;

%Use one slice of the image
%Ex=CAvrg(Ceny,:);

%OR add n slices around the center of the cloud
n=10;
for i=Ceny-(n/2):Ceny+(n/2)
Ex=Ex+CAvrg(i,:);
end
Ex=Ex/(n+1);

%plot(Ex)
ExMin=mean(Ex(1200:length(Ex)));
NormEx=Ex-ExMin;

%smoothing NormEx
windowSize = 1; % 1== no smoothing
fEx=filter(ones(1,windowSize)/windowSize,1,NormEx);

lnEx=log(fEx);

p=polyfit(x(630:830),lnEx(630:830),2);
A2=p(1);
A1=p(2);
A0=p(3);
sigmax=sqrt(-1/(2*A2));
mux=A1*sigmax^2;
Ax=exp(A0+mux^2/(2*sigmax^2));
subplot(2,2,3),
plot(x,fEx,x,Ax*exp(-(x-mux).^2/(2*sigmax.^2)));
rx=sqrt(2)*px_sizew*sigmax./Magn % in mm

%Regression Dianogsis: R-squared
x_bar=mean(fEx);
x_residual=(fEx-Ax*exp(-(x-mux).^2/(2*sigmax.^2))).^2;
Sx_error=sum(x_residual);
dev_Avrgx=(fEx-x_bar).^2;
Sx_total=sum(dev_Avrgx);
Rx_Squared=1-(Sx_error/Sx_total);

%Regression Dianogsis: comparing integration
rez_x=abs((trapz(fEx)-trapz(Ax*exp(-(x-mux).^2/(2*sigmax.^2))))/trapz(fEx))

%% ### finding the Gaussian fit for Y ###

Ey=0;

%Use one slice of the image
%Ey=CAvrg(:,Cenx);
%Cenx=Cenx;
%OR add n slices around the center of the cloud

n=10;
for i=Cenx-(n/2):Cenx+(n/2)
Ey=Ey+CAvrg(:,i);
end
Ey=Ey/(n+1);

EyMin=mean(Ey(1:100));
NormEy=Ey-EyMin;

%smoothing NormEy
windowSize = 1;
fEy=filter(ones(1,windowSize)/windowSize,1,NormEy);

lnEy=log(double(fEy));
y=(1:px_height)’;

p=polyfit(y(500:630),lnEy(500:630),2);
A2=p(1);
A1=p(2);
A0=p(3);
sigmay=sqrt(-1/(2*A2));
muy=A1*sigmay^2;
Ay=exp(A0+muy^2/(2*sigmay^2));
subplot(2,2,4),
plot(y,fEy,y,Ay*exp(-(y-muy).^2/(2*sigmay.^2)))
ry=sqrt(2)*px_sizeh*sigmay./Magn %in mm

%Regression Dianogsis: R-squared
y_bar=mean(fEy);
y_residual=(fEy-Ay*exp(-(y-muy).^2/(2*sigmay.^2))).^2;
Sy_error=sum(y_residual);
dev_Avrgy=(fEy-y_bar).^2;
Sy_total=sum(dev_Avrgy);
Ry_Squared=1-(Sy_error/Sy_total)

%Regression Dianogsis: comparing integration
rez_y=abs((trapz(fEy)-trapz(Ay*exp(-(y-muy).^2/(2*sigmay.^2))))/trapz(fEy))

Advertisements

Tags: , ,

About vekin

I'm mostly a scientist, but occasionally a writer and an artist, and most definitely a dreamer.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: