
%%% --- IQEM, STAGE 2, STEP 1 ----------------------------------------------------------------------- %%%
%%% THE FITTED 2D QUADRATIC FUNCTIONS ARE EVALUATED AT EACH PIXEL LOCATION, AND FOR EACH INTENSITY 
%%% QUANTILE.  CALCULATED VALUES ARE SAVED AND IMPLEMENTED IN THE FOLLOWING IMAGE CORRECTION STAGE.  
%%% THIS 'PRE-CALCULATION' IS DONE TO INCREASE THE EFFICIENCY OF THE IMAGE CORRECTION STAGE

clear
rows = cellstr( ['A';'B';'C';'D';'E';'F';'G';'H'] );											% SPECIFY IMAGE FILE ROW NAMES
cols = cellstr( ['1 ';'2 ';'3 ';'4 ';'5 ';'6 ';'7 ';'8 ';'9 ';'10';'11';'12'] );				% SPECIFY IMAGE FILE COLUMN NAMES
flds = cellstr( ['1';'2';'3';'4';'5';'6';'7';'8'] );											% SPECIFY IMAGE FILE FIELD NAMES
chns = cellstr( ['S475_20X - HQ535_50M ';'HQ535_50X - HQ620_60M';'D360_40X - HQ460_40M '] );	% SPECIFY IMAGE FILE CHANNEL NAMES
chn_names = cellstr( ['FITC';'Cy3 ';'DAPI'] );													% SPECIFY ABBREVIATED CHANNEL NAMES

% SPECIFY PATH NAME EXTENSIONS FOR PLATES AND CORRESPONDING BIAS CORRECTED IMAGES %
paths = cellstr( ['090624_Lysed_A1_Images/   ';'090624_Lysed_B1_Images/   ';'090625_Lysed_A2_Images/   ';
    '090625_Lysed_B2_Images/   ';'090626_Lysed_A3_Images/   ';'090626_Lysed_B3_Images/   ';
    '090624_NonLysed_A1_Images/';'090624_NonLysed_B1_Images/';'090625_NonLysed_A2_Images/';'090625_NonLysed_B2_Images/';'090626_NonLysed_A3_Images/';'090626_NonLysed_B3_Images/'] );
outs = cellstr( ['090624_Lysed_A1_Images_   ';'090624_Lysed_B1_Images_   ';'090625_Lysed_A2_Images_   ';
    '090625_Lysed_B2_Images_   ';'090626_Lysed_A3_Images_   ';'090626_Lysed_B3_Images_   ';
    '090624_NonLysed_A1_Images_';'090624_NonLysed_B1_Images_';'090625_NonLysed_A2_Images_';'090625_NonLysed_B2_Images_';'090626_NonLysed_A3_Images_';'090626_NonLysed_B3_Images_'] );

for i_path = 1:12	% LOOP OVER 12 PLATES, 3 CHANNELS
for i_chn = 1:3

%%% 1. SET PATH FOR COEFF AND DATA FILES %%%
outpath = ['/mnt/share/Chip_GRN/Elo/AllPlates_Data/Illume_2_ci13/grid_distributions_x3/' char(paths(i_path))];

%%% 2. LOAD FITTED QUADRATIC COEFFICIENTS %%%
outfile = [char(outpath) char(outs(i_path)) char(chn_names(i_chn)) '_Coeffs.mat']
load (outfile);

%%% 3. DETERMINE THE PERCENTILE GRID RESOLUTION %%%
[num_pctls tmp] = size(squeeze(Coeffs(:,1)));

%%% 4. SPECIFY NAME OF OUTFILE FOR INTERPOLATED QUANTILES MATRIX %%%
outfile = [char(outpath) char(outs(i_path)) char(chn_names(i_chn)) '_all_pct_contours']


% FOR EACH QUANTILE/PERCENTILE, A 2D CONTOUR OF INTENSITIES IS STORED 
pct_contours = zeros(1028,1300,num_pctls);

for i_pct = 1:num_pctls
  [x,y] = meshgrid(1:1300,1:1028);
  x2 = x.*x; y2 = y.*y; xy = x.*y; o = ones(1028,1300);
  coeffs_contour = squeeze(Coeffs(i_pct,:)');
  FittedImage = exp( reshape([x2(:) y2(:) xy(:) x(:) y(:) o(:)] * coeffs_contour, [1028 1300]) ) -1; 
  pct_contours(:,:,i_pct) = FittedImage;
end

%%% SAVE MATRIX OF INTERPOLATED QUANTILES %%%
save (outfile,'pct_contours');

end
end


%%% --- END OF IQEM, STAGE 2, STEP 1 ---------------------------------------------------------------- %%%
