
%%% --- IQEM, STAGE 2, STEPS 1 AND 2 ---------------------------------------------------------------- %%%
%%% THE IQEM NU BIAS CORRECTION IS APPLIED TO EACH IMAGE; THIS IS DONE VIA A MAPPING OF THE PIXEL     %%%
%%% INTENSITIES TO THE REFERENCE INTENSITY QUANTILE FUNCTION                                          %%%

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. SPECIFY PATHNAME OF IQEM DATA FILES FROM PREVIOUS STEPS %%%
outpath = ['/mnt/share/Chip_GRN/Elo/AllPlates_Data/Illume_2_ci13/grid_distributions_x3/' char(paths(i_path))];

%%% 2. OUTPUT PATH FOR CORRECTED IMAGES %%%
outpath2 = ['/mnt/share/Chip_GRN/Elo/AllPlates_Data/Illume_2_ci13/grid_distributions_x3/' char(outs(i_path)) char(chn_names(i_chn)) '_IQEM/'];

%%% 3. LOAD FITTED INTENSITY CONTOURS FOR EACH QUANTILE, FROM PREVIOUS STEP
outfile = [char(outpath) char(outs(i_path)) char(chn_names(i_chn)) '_all_pct_contours']
load (outfile);

%%% 4. SET SOURCE IMAGE PATH %%%
image_path = ['/mnt/share/Chip_GRN/Elo/AllPlates_Data/Illume_2_ci13/grid_distributions_x3/' char(paths(i_path))];

%%% 5. SOURCE IMAGE FILE NAME SUFFIX ... THESE ARE THE DISTORTED IMAGES THAT WILL BE CORRECTED %%%
ifile_suffix = ').tif';

%%% 6. CORRECTED IMAGE FILE NAME SUFFIX %%%
ofile_suffix = ')_IQEM.tif';

%%% 7. SPECIFY PERCENTILE GRID RESOLUTION - THE DEFAULT GRID HAS 33 ELEMENTS %%%
pct_values = [.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];
[tmp num_pctls] = size(pct_values);

% EXTRACT THE REFERENCE IQF FROM THE CENTRE OF THE IMAGE
ref_iqf = squeeze(pct_contours(514,650,:));

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

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

corrected_image = zeros(1028,1300);

% LOOP OVER THE PIXELS WITHIN EACH IMAGE
for i_xgrid = 1:1028
for i_ygrid = 1:1300

    I_pixel = double(I(i_xgrid,i_ygrid));
	% COMPUTE THE INTENSITY QUANTILES THAT BRACKET THE MEASURED PIXEL INTENSITY
    test_pct = double(I_pixel - pct_contours(i_xgrid,i_ygrid,:)) ;
    i_pct=1; pct_right_index=1; pct_left_index=0;
    while i_pct < num_pctls+1
        if test_pct(i_pct) < 0
            pct_right_index = i_pct; pct_left_index = i_pct -1; i_pct = 99;
        end
        i_pct = i_pct+1;
        if i_pct == num_pctls+1
            pct_left_index = num_pctls; pct_right_index=num_pctls+1; 
		end
    end
    
% USE LINEAR INTERPOLATION TO MAP PIXEL INTENSITIES TO CORRESPONDING PERCENTILE AND THEN TO REFERENCE IQF INTENSITY VALUE
if pct_right_index == 1
   I_ref_pixel = ref_iqf(1);		% CLIP TO LOWEST INTENSITY QUANTILE OF REFERENCE IQF
elseif pct_left_index == num_pctls
   I_ref_pixel = ref_iqf(num_pctls);% CLIP TO HIGHEST INTENSITY QUANTILE OF REFERENCE IQF
else
   Ipct_right = pct_contours(i_xgrid,i_ygrid,pct_right_index); pct_right = pct_values(pct_right_index);
   Ipct_left = pct_contours(i_xgrid,i_ygrid,pct_left_index); pct_left = pct_values(pct_left_index);
   % LINEAR INTERPOLATION TO MAP THE PIXEL INTENSITY TO REFERENCE IQF PERCENTILE
   pct_pixel = ((Ipct_right - I_pixel )*pct_left + (I_pixel - Ipct_left)*pct_right)/(Ipct_right - Ipct_left);
   % LINEAR INTERPOLATION TO MAP THE INTENSITY PERCENTILE TO REFERENCE IQF INTENSITY VALUE
   I_ref_left = ref_iqf(pct_left_index); I_ref_right = ref_iqf(pct_right_index);
   I_ref_pixel = ((pct_right - pct_pixel)*I_ref_left + (pct_pixel - pct_left)*I_ref_right)/(pct_right - pct_left);
end
corrected_image(i_xgrid,i_ygrid) = I_ref_pixel;
end
end

% SAVE THE IQEM CORRECTED IMAGE 
corrected_image_out = uint16(corrected_image);
outfile = [char(outpath2) char(rows(i_row)) ' - ' char(cols(i_col)) '(fld ' char(flds(i_fld)) ' wv ' char(chns(i_chn)) char(ofile_suffix)];
imwrite(corrected_image_out,outfile,'tif');

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

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


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