Tag Archive | MATLAB

MATLAB Code for Area Fit

%% Cloud area and Volume

% 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
px_area=px_sizew*px_sizeh; %area for each pixel in mm^2

%% load image
C1=imread(‘C:\Users\Anarisil\Desktop\MOTPhase2\20120502\7-4\50ms.jpg’);

CAvrg=double(C1);

%% Dimension of the picture

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

%% finding the maximum value of that point and the 1/e value
MaxInt=max(max(CAvrg));
threshold=MaxInt/exp(1);

%% %% counting
px_count=0;

for i=1:px_height
for j=1:px_width
if CAvrg(i,j)>=threshold
px_count=px_count+1;
else
end
end
end

%px_count will give the number of pixels with values above the threshold

%% determining the volume
%assumption: x=y=z
V(1,1)=px_count*px_area;
V(2,1)=4/(3*sqrt(pi))*(V(1,1))^(3/2); %in mm^3

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