
%%% --- IQEM, STAGE 1, STEP 3 ----------------------------------------------------------------------- %%%
%%% 2D QUADRATIC FUNCTIONS ARE FITTED TO THE SPATIAL INTENSITY VARIATION FOR EACH INTENSITY QUANTILE  %%%
%%% USING MULTIPLE LINEAR REGRESSION                                                                  %%%

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 DATA FILES (GRID PERCENTILES AND COEFFICIENTS) %%%
outpath = ['/mnt/share/Chip_GRN/Elo/AllPlates_Data/Illume_2_ci13/grid_distributions_x3/' char(paths(i_path))];

%%% 2. DETERMINE SPATIAL GRID RESOLUTION %%% 
[ num_grid_x num_grid_y tmp] = size(grid_percentiles);
num_pix_y = floor(1300/num_grid_y); cent_pix_y = (num_pix_y+1)/2; % DETERMINE NUMBER OF PIXELS PER GRID ELEMENT AND CENTER OF ELEMENTS
num_pix_x = floor(1028/num_grid_x); cent_pix_x = (num_pix_x+1)/2;

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

%%% LOAD MATLAB FILE CONTAINING EMPIRICAL IQFs FOR SPECIFIED PLATE AND CHANNEL %%%
outfile = [char(outpath) char(outs(i_path)) char(chn_names(i_chn)) '_IQEM_pctIx3']; load (outfile);

Coeffs = zeros((num_pctls-2),6); % MATRIX OF COEFFICIENTS (THERE ARE 6 QUADRATIC COEFFICIENTS FOR EACH INTENSITY PERCENTILE)

% FIT 2D QUADRATIC FOR EACH LOG INTENSITY PERCENTILE
for i_pct = 2:(num_pctls-1)
   Istratum = log(grid_percentiles(:,:,i_pct) +1); 
   [x,y] = meshgrid(1:size(Istratum,2),1:size(Istratum,1)); 			% x=LONGSIDE; y=SHORTSIDE
   x = (x-1)*num_pix_y + cent_pix_y; y = (y-1)*num_pix_x + cent_pix_x;
   x2 = x.*x; y2 = y.*y; xy = x.*y; o = ones(size(Istratum));
   Coeffs(i_pct-1,:) = [x2(:) y2(:) xy(:) x(:) y(:) o(:)] \ double(Istratum(:));
end

%%% SAVE QUADRATIC COEFFICIENTS IN MATLAB AND ASCII FORMATS %%%
outfile = [char(outpath) char(outs(i_path)) char(chn_names(i_chn)) '_IQEM_Coeffs']
save (outfile,'Coeffs');

outfile = [char(outpath) char(outs(i_path)) char(chn_names(i_chn)) '_IQEM_Coeffs.txt']
save (outfile,'Coeffs','-ascii');

end
end


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


