
%%% --- IQEM, STAGE 1, STEPS 1 AND 2 ---------------------------------------------------------------- %%%
%%% PIXEL INTENSITIES WITHIN EACH SPATIAL GRID ELEMENT AND OVER ALL IMAGES IN EACH EXPERIMENTAL BATCH %%%
%%% ARE AGGREGATED, AND THE CORRESPONDING EMPIRICAL IQFs ARE DETERMINED                               %%%

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	
image_num=0;			% IMAGE NUMBER COUNTER (FOR THE BATCH)

%%% 1. OUTPUT DIRECTORY FOR GRID PERCENTILES FILE %%%
outpath = ['/mnt/share/Chip_GRN/Elo/AllPlates_Data/Illume_2_ci13/grid_distributions_x3/' char(paths(i_path))];
%%% 2. SET SOURCE IMAGE PATH %%%
image_path = ['/mnt/share/Chip_GRN/Elo/AllPlates_Data/Illume_2_ci13/grid_distributions_x3/' char(paths(i_path))];
%%% 3. SET SOURCE IMAGE FILENAME SUFFIX %%%
ifile_suffix = ').tif';
%%% 4. SET GRID PERCENTILES OUTFILE AND PATH %%%
outfile = [char(outpath) char(outs(i_path)) char(chn_names(i_chn)) '_IQEM_pctIx3']
%%% 5. SET ASCII PERCENTILE PROFILES OUTFILE PREFIX %%%
outfile = [char(outpath) char(outs(i_path)) char(chn_names(i_chn)) '_IQEM_pctIx3']
%%% 6. SPECIFY THE NUMBER OF IMAGES %%%
num_images=768
%%% 7. SET THE SPATIAL GRID RESOLUTION %%%
num_grid_y = 20; num_pix_y = floor(1300/num_grid_y); % NUMBER OF GRID ELEMENTS; NUMBER OF PIXELS PER GRID ELEMENT (LONG AXIS)
num_grid_x = 15; num_pix_x = floor(1028/num_grid_x); % NUMBER OF GRID ELEMENTS; NUMBER OF PIXELS PER GRID ELEMENT (SHORT AXIS)
grid_intensities = zeros(num_grid_x,num_grid_y,num_images*num_pix_y*num_pix_x); % TEMPORARY VARIABLE THAT STORES AGGREGATED PIXEL INTENSITES OVER ALL IMAGES
%%% 8. SPECIFY THE PERCENTILE GRID RESOLUTION ... CREATE ARRAY TO STORE THE MEASURED PERCENTILES/QUANTILES %%%
p = [0,.25,.5,1,2.5,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,82.5,85,87.5,90,92.5,95,97.5,99,99.5,99.75,99.9,99.95,99.99,100];
% THIS IS THE 33-ELEMENT, DEFAULT PERCENTILE GRID
[tmp num_pctls] = size(p)
grid_percentiles = zeros(num_grid_x,num_grid_y,num_pctls); % THIS VARIABLE STORES THE EMPIRICAL IQFs OVER THE SPATIAL GRID

for i_row = 1:8			% LOOP OVER 8 ROWS, 12 COLUMNS, 8 FIELDS
for i_col = 1:12
for i_fld = 1:8	

%%% SET IMAGE FILE NAME %%%
image_file = [char(rows(i_row)) ' - ' char(cols(i_col)) '(fld ' char(flds(i_fld)) ' wv ' char(chns(i_chn)) char(ifile_suffix) ];

I = imread([image_path image_file]);
image_num = image_num+1; 

% FOR EACH IMAGE, EXTRACT PIXELS FOR EACH SPATIAL GRID ELEMENT AND ACCUMULATE
for y_grid = 1:num_grid_y   
for x_grid = 1:num_grid_x   
    y_start = (y_grid - 1)*num_pix_y + 1; y_end = (y_grid)*num_pix_y;
    x_start = (x_grid - 1)*num_pix_x + 1; x_end = (x_grid)*num_pix_x;
    index_start = num_pix_y*num_pix_x* (image_num-1) + 1; index_end = num_pix_y*num_pix_x*image_num; 
    intensities = I(x_start:x_end, y_start:y_end); % EXTRACT PIXEL INTENSITIES FROM SPECIFIC GRID ELEMENTS
    grid_intensities(x_grid,y_grid,index_start:index_end)= double(intensities(:)) ; % ACCUMULATE PIXEL INTENSITIES FROM ALL GRID ELEMENTS
end
end

end % END OF i_fld LOOP
end % END OF i_col LOOP
i_row
end  % END OF i_row LOOP


% COMPUTE EMPIRICAL IQF ON PERCENTILE GRID, FOR EACH SPATIAL GRID ELEMENT
for y_grid = 1:num_grid_y   
for x_grid = 1:num_grid_x   
    grid_percentiles(x_grid,y_grid,1:num_pctls) = prctile(grid_intensities(x_grid,y_grid,:),p);
end
y_grid
end

save (outfile,'grid_percentiles'); % SAVE THE GRID OF IQFs


% FOR EACH PERCENTILE/QUANTILE, SAVE THE SPATIAL CONTOUR FUNCTION FOR POSSIBLE PLOTTING/DIAGNOSTICS
for i_pct = 1:num_pctls
    outfile = [char(outpath) char(outfile_prefix) num2str(i_pct) '.txt'];
    grid_out = grid_percentiles(:,:,i_pct);
    save (outfile,'grid_out','-ascii');
end

end  % END OF i_chn LOOP
end  % END OF i_path LOOP


%%% --- END OF IQEM, STAGE 1, STEPS 1 AND 2 --------------------------------------------------------- %%%
