%sample program to grab imahge and process
%111609 TDM

clear all
clc

[filename, pathname] = uigetfile('*.jpg', 'Pick an image');
if isequal(filename,0)
   disp('User selected Cancel')
else
   disp(['User selected: ', fullfile(pathname, filename)])
end

[X,map] = imread(fullfile(pathname, filename),'jpeg');

%convert to grayscale
X       = uint8(round(sum(double(X),3)/3));
X=double(X);

figure(100);imagesc(X);colormap(gray)

%find max of image
[r_max c_max] = find(max(max(X))==X);
disp(['max row number = ' num2str(r_max)] )
disp(['max col number = ' num2str(c_max)] )

%plot profiles at row, column
row         = input('enter row number : ');
column      = input('enter column number : ');

[nr nc]     = size(X);

figure(101);plot(1:nc,X(row,:),1:nr,X(:,column))
legend('Horizontal (across columns)','Vertical (across rows)'  )
title('Gray-Scale Image')
xlabel('Column Number')
ylabel('Row Number')

%plot of sinc^2 function
%SINC Sin(pi*x)/(pi*x) function
sinc_scale  = 1;    %change this number to change the width of the sinc
                    %default value =1.
x           = linspace(-3,3,400);
sinc_sq     = sinc(sinc_scale*x).^2;
figure(102);plot(x,sinc_sq)
title('Sinc^2 Function')
ylabel('Relative Irradiance (max value = 255)')

%plot of somb^2
%somb = J1(pi*rho)/pi/rho
somb_scale  = 1;    %change this number to change the width of the somb
                    %default value =1.
somb    = @(x) 2*besselj(1,pi*x)./x/pi;
if x ==0 
    x = 1e-10;
end
somb_sq     = abs(somb(somb_scale*x)).^2;
figure(103);plot(x,somb_sq)
title('Somb^2 Function')
ylabel('Relative Irradiance (max value = 255)')

%Plot normalized and shifted experimental horizontal profile (across the
%columns)
%this graph also overlays the scaled sinc^2 and somb^2 functions
h_shift                 = -70;  %change this number to shift a number of pixels left (-)
                              %or right (+) of the displayed profile
horizontal_profile      = double(X(row,:))/double(max(X(row,:)));
h_profile_shifted       = interp1(1:nc,horizontal_profile,(1:nc)-h_shift);
indx                    = find(isnan(h_profile_shifted));
h_profile_shifted(indx) = repmat(0,1,length(indx));
%setup for sinc_sq and somb_sq
sinc_scale              = 0.85;
somb_scale              = 1;
x                       = linspace(-3,3,nc);
sinc_sq                 = sinc(sinc_scale*x).^2;
somb_sq                 = abs(somb(somb_scale*x)).^2;
%plot figure
figure(104);plot(1:nc,h_profile_shifted ...
    ,1:nc,sinc_sq ...
    ,1:nc,somb_sq)
title('Normalized Horizontal Profile (across the columns)')
ylabel('Normalized Irradiance')
legend('Exp profile','sinc^2','somb^2')
grid

%Plot normalized and shifted experimental vertical profile (across the rows)
v_shift                 = -40;  %change this number to shift a number of pixels up (-)
                              %or down (+) of the displayed profile
vertical_profile      = double(X(:,column))/double(max(X(:,column)));
v_profile_shifted       = interp1(1:nr,vertical_profile,(1:nr)-v_shift);
indx                    = find(isnan(v_profile_shifted));
v_profile_shifted(indx) = repmat(0,1,length(indx));
%setup for sinc_sq and somb_sq
sinc_scale              = 4;
somb_scale              = 1;
x                       = linspace(-3,3,nr);
sinc_sq                 = sinc(sinc_scale*x).^2;
somb_sq                 = abs(somb(somb_scale*x)).^2;
%plot figure
figure(105);plot(1:nr,v_profile_shifted ...
    ,1:nr,sinc_sq ...
    ,1:nr,somb_sq)
title('Normalized Vertical Profile (across the rows)')
ylabel('Normalized Irradiance')
legend('Exp profile','sinc^2','somb^2')
grid