% DNAPL Geometry Evolution

%% Images
clear all; %close all; %clc
tic
global fl file_dir_OE
global norm_sub_thin_G_avg norm_sub_thin_R_avg norm_sub_thin_B_avg
global mean_R mean_G mean_B
global std_R std_G std_B
global gs_noise
global OEI_thin_R_std OEI_thin_G_std OEI_thin_B_std
global OEI_thin_R_avg OEI_thin_G_avg OEI_thin_B_avg
global jg iLb_R_avg iLb_G_avg iLb_B_avg

jg = 0;
% step 1: read images and define regions of interest (ROI)
i1 = 1480; 
i2 = 4330;
j1 = 380;
j2 = 3230;
pt_1 = [5323, 966];
pt_2 = [5292, 2749];

nxshift = 20;
nyshift = 20;

% Aperture Field
file_dir = '/karst/data/Masoud/Masoud_Exp/exp_01/Image_processing_code/04_aperture_field/';
NewFile_diff = [file_dir,'Aperture_field.mat'];
A_F_str = load(NewFile_diff);
A_F = A_F_str.A_F;
A_F_roi = A_F(i1:i2,j1:j2);
% figure; imagesc(A_F_roi); impixelinfo

% pixel size measurement (area) 
in22cm2 = 6.452;
cell_length = 6; % in
l_pixel_count = 2882; % total no. of pixels in 6 inch
pixel_size_in = cell_length / l_pixel_count; % in
pixel_area_in = pixel_size_in^2; % area of each pixel (in^2)
pixel_area = pixel_area_in * in22cm2; % area of each pixel (cm^2)
pixel_size = sqrt(pixel_area); % size of each pixel (cm)

% Total volume of the cell
V_cell = sum(A_F_roi(:)) * pixel_area; % cm^3 , ml


% Clear Water Red Light, Base Image
Filename_clear_water_red = '/karst/data/Masoud/Masoud_Exp/exp_01/02_aperture_field/water_saturated_cell/registered_avg.mat';
im_red_str=load(Filename_clear_water_red);
im_red = im_red_str.registered_avg;
Base = im_red;
% figure; imagesc(Base); impixelinfo; title('im_red')

% Residual TCE 
Filename_Res_TCE = '/karst/scratch/RESTORE/karst/data/Masoud/Masoud_Exp/exp_04/06_tce_geometry_w_dye/registered_w_CL.mat';
im_TCE_str = load(Filename_Res_TCE);
im_TCE = im_TCE_str.reg_TCE;
B = im_TCE(i1:i2,j1:j2,:);
%figure; imagesc(B); impixelinfo; title('im Residual TCE')

% DNAPL Initial Phase Geometrey (clear water)
file_dir_initial = '/karst/scratch/RESTORE/karst/data/Masoud/Masoud_Exp/exp_04/08_tce_geometry_w_clear_water/';
NewFile_initial = [file_dir_initial,'reg_clear_tce.mat'];
initial_TCE_str = load(NewFile_initial);
initial_TCE = initial_TCE_str.reg_clear_TCE;
initial_TCE_roi = initial_TCE(i1:i2,j1:j2,:);
% figure; imagesc(initial_TCE_roi); impixelinfo; title('initial_TCE_roi')

% iTCE, binary image of initail tce geometry.
file_dir_iBW = '/karst/scratch/RESTORE/karst/data/Masoud/Masoud_Exp/exp_04/08_tce_geometry_w_clear_water/';
NewFile_iBW = [file_dir_iBW,'Initial_TCE_test_03.mat'];
iTCE_str = load(NewFile_iBW);
iTCE = iTCE_str.iTCE_new;
iTCE_filled = imfill(iTCE,4,'hole');
iTCE_hole = iTCE_filled & ~iTCE;
gs_noise = iTCE_str.gs_noise;
gs_noise(920:980, 2640:2740) = 0;
% figure; imagesc(iTCE); impixelinfo; title('iTCE'); colormap(gray)
% figure; imagesc(gs_noise); impixelinfo; title('noose'); colormap(gray)
se = strel('square',5);
iTCE_edge = bw_edge(iTCE, 1);
% figure; imagesc(iTCE_edge); impixelinfo; title('iTCE'); colormap(gray)
% iTCE = iTCE & ~iTCE_edge;
iTCE_3D = [iTCE, iTCE, iTCE];
initialTCE = iTCE;

[ImageHeight ImageWidth] = size(iTCE);
dnaplsize = ImageHeight * ImageWidth;

iSaturation = sum(iTCE(:))/dnaplsize;
% initialTCE_3D = [initialTCE, initialTCE, initialTCE];
% Windows Images! 
iTCE_win = initial_TCE_roi;
iTCE_win(~iTCE_3D) = 0;
% figure; imagesc(iTCE_win) ; impixelinfo;
blue_noise = initial_TCE_roi(:,:,1) < 45000 & initial_TCE_roi(:,:,3) > 16000 & ~imdilate(gs_noise, se);
blue_noise = bwareaopen(blue_noise, 150, 4);
blue_noise = imdilate(blue_noise, se);
blue_noise = imfill(blue_noise, 'hole');% &  ~iTCE_hole;
%   figure; imagesc(blue_noise); impixelinfo; colormap(gray)
V_blue_noise = sum(A_F_roi(blue_noise)) * pixel_area; % cm^3 , ml

%% Oxidation Experiment
% reading the images from Oxidation Experiment
file_dir_OE = '/karst/scratch/RESTORE/karst/data/Masoud/Masoud_Exp/exp_04/09_KMnO4_injection/Main/';
fl = textread([file_dir_OE,'MasoudTIFF'],'%s');
file_excel ='/karst/scratch/RESTORE/karst/data/Masoud/Masoud_Exp/exp_04//timexls.xls';
time_vec = xlsread(file_excel,'B1:B1016');
dummy = 0;
num_imgs = length(time_vec);

OEI_pre = initial_TCE;
OEI_pre_roi = initial_TCE_roi;
pre_roi_R = double(OEI_pre_roi(:,:,1));
pre_roi_G = double(OEI_pre_roi(:,:,2));
pre_roi_B = double(OEI_pre_roi(:,:,3));
%   figure; imagesc(OEI_pre_roi); impixelinfo; title('OEI_pre_roi');

i = 0;
% time_new(i) = 0;
BW_DNAPL = initialTCE;
BW_DNAPL_3D = [BW_DNAPL, BW_DNAPL, BW_DNAPL];
DNAPL_edge = bw_edge(BW_DNAPL, 3);
BW_DNAPL_thin = (BW_DNAPL) & (~DNAPL_edge) & (~gs_noise);
% figure; imagesc(DNAPL_edge); impixelinfo

mean_R = 52000;
mean_G = 44000;
mean_B = 20000;

std_R = 2200;
std_G = 2400;
std_B = 1100;

pre_TCE_R_avg = mean(pre_roi_R(BW_DNAPL_thin));
pre_TCE_G_avg = mean(pre_roi_G(BW_DNAPL_thin));
pre_TCE_B_avg = mean(pre_roi_B(BW_DNAPL_thin));

pre_TCE_R_std = std(pre_roi_R(BW_DNAPL_thin));
pre_TCE_G_std = std(pre_roi_G(BW_DNAPL_thin));
pre_TCE_B_std = std(pre_roi_B(BW_DNAPL_thin));
% figure; imagesc(OEI_pre_roi); impixelinfo; title('OEI_pre_roi');

%% n = 1;

n = 1;
i = i + 1;
time = time_vec(n); % min
time_new(i) = time;
deltat(i) = 0;
file = char(fl(n));
file_name = [file(1:end-5),'_reg.mat'];
NewFile_OEI = [file_dir_OE,file_name];
% now load the image from the new directory
OEI_str = load(NewFile_OEI);
OEI = OEI_str.OEI_reg;
[reg_OEI x_shift y_shift] = reg_fun(OEI_pre, OEI, pt_1, pt_2, nxshift, nyshift);
% shift_mat(i,:) = [x_shift(1) y_shift(1) x_shift(2) y_shift(2)]; 
reg_OEI_1 = reg_OEI;
OEI_roi = reg_OEI(i1:i2,j1:j2,:);
file_name_reg = [file(1:end-5),'_final_registered.mat'];
NewFile = [file_dir_OE, file_name_reg];
% now save the image in the new directory
save(NewFile,'reg_OEI')
% figure; imagesc(reg_OEI); impixelinfo; title('reg_OEI');
% figure; imagesc(OEI_roi); impixelinfo; title('OEI_roi');
% Windows OEI Images!
OEI_Win = OEI_roi;
OEI_Win(~BW_DNAPL_3D) = 0;
% figure; imagesc(OEI_Win); impixelinfo; title('OEI_Win');
Diff_size(i) = 0;

DNAPL_Size_vec(i) = sum(BW_DNAPL(:));
V_Sat_vec(i) = sum(A_F_roi(BW_DNAPL)) / sum(A_F_roi(:));
    
interfacial_edge = bw_edge(imfill(BW_DNAPL,'holes'),1);
interfacial_area(i) = sum(A_F_roi(interfacial_edge)) * pixel_size; % cm^2 

OEI_R = double(OEI_roi(:,:,1));
OEI_G = double(OEI_roi(:,:,2));
OEI_B = double(OEI_roi(:,:,3));

% aqueous phase avg
BW_aq = (~BW_DNAPL) & (~gs_noise);
% figure; imagesc(BW_aq); impixelinfo
aq_avg_R(i) = mean(OEI_R(BW_aq));
aq_avg_G(i) = mean(OEI_G(BW_aq));
aq_avg_B(i) = mean(OEI_B(BW_aq));

OEI_thin_R_avg(i) = mean(OEI_R(BW_DNAPL_thin));
OEI_thin_G_avg(i) = mean(OEI_G(BW_DNAPL_thin));
OEI_thin_B_avg(i) = mean(OEI_B(BW_DNAPL_thin));

OEI_thin_R_std(i) = std(OEI_R(BW_DNAPL_thin));
OEI_thin_G_std(i) = std(OEI_G(BW_DNAPL_thin));
OEI_thin_B_std(i) = std(OEI_B(BW_DNAPL_thin));

% Normalizing OEI image using thin_avg and thin_std
norm_OEI(:,:,1) = uint16((OEI_R - OEI_thin_R_avg(i)) * std_R / OEI_thin_R_std(i) + mean_R);
norm_OEI(:,:,2) = uint16((OEI_G - OEI_thin_G_avg(i)) * std_G / OEI_thin_G_std(i) + mean_G);
norm_OEI(:,:,3) = uint16((OEI_B - OEI_thin_B_avg(i)) * std_B / OEI_thin_B_std(i) + mean_B);
inorm_OEI = norm_OEI;
% norm_OEI = OEI_roi;         
file_dir = '/karst/scratch/RESTORE/karst/data/Masoud/Masoud_Exp/exp_04/09_KMnO4_injection/Binary_DNAPL/';
file_name_bw = [file(1:end-5),'_bw.mat'];
NewFile_bw = [file_dir, file_name_bw];
% now save the image in the new directory
save(NewFile_bw,'BW_DNAPL','OEI_roi')
%% n = 10;
% BW_gas_total = zeros(size(BW_DNAPL));
BW_diff_th = 20 * double(iTCE);
% BW_gas_th = 65 * double(iTCE);
BW_TCE_new_all = zeros(size(BW_DNAPL));
norm_sub_total = uint16(zeros(size(OEI_roi)));

for n = [40:20:1000, 1016];
    i = i + 1;
    time = time_vec(n); % min
    time_new(i) = time;
    deltat(i) = time_new(i) - time_new(i-1);
    OEI_pre = reg_OEI;
    norm_pre = norm_OEI;
    BW_edge = bw_edge(BW_DNAPL, 7);
    [reg_OEI, BW_DNAPL_thin, norm_OEI, norm_sub] = OEI_fun(n, i, OEI_pre, norm_pre, BW_DNAPL); %, BW_gas_total);
    
    OEI_roi = reg_OEI(i1:i2,j1:j2,:);
    OEI_R = OEI_roi(:,:,1);
    OEI_G = OEI_roi(:,:,2);
    OEI_B = OEI_roi(:,:,3);
 
    OEI_Win = OEI_roi;
    OEI_Win(~BW_DNAPL_3D) = 0;
    
    OEI_pre_roi = OEI_pre(i1:i2,j1:j2,:); 
    % figure; imagesc(OEI_pre_roi); impixelinfo; title('OEI pre');
    
    % norm sub image
    norm_sub_Win = norm_sub;
    norm_sub_Win(~BW_DNAPL_3D) = 0;
    norm_sub_R = double(norm_sub(:,:,1));
    norm_sub_G = double(norm_sub(:,:,2));
    norm_sub_B = double(norm_sub(:,:,3));
    norm_sub_RG = (norm_sub_R + norm_sub_G)  / 2;
    
    % norm sub rev image
    norm_sub_rev = imsubtract(OEI_roi, OEI_pre_roi);
    % figure; imagesc(norm_sub_rev); impixelinfo ;
    
    % norm sub total image
    norm_sub_total = imsubtract(inorm_OEI, norm_OEI);
    norm_sub_total_R = double(norm_sub_total(:,:,1));
    norm_sub_total_G = double(norm_sub_total(:,:,2));
    norm_sub_total_B = double(norm_sub_total(:,:,3));
    
    
    thin_sub_total_R_avg_1(i) = mean(norm_sub_total_R(BW_DNAPL_thin));
    thin_sub_total_G_avg_1(i) = mean(norm_sub_total_G(BW_DNAPL_thin));
    thin_sub_total_B_avg_1(i) = mean(norm_sub_total_B(BW_DNAPL_thin));
    
    norm_sub_total_2(:,:,1) = uint16(norm_sub_total_R - thin_sub_total_R_avg_1(i));
    norm_sub_total_2(:,:,2) = uint16(norm_sub_total_G - thin_sub_total_G_avg_1(i));
    norm_sub_total_2(:,:,3) = uint16(norm_sub_total_B - thin_sub_total_B_avg_1(i));
    
    norm_sub_total_R_2 = double(norm_sub_total_2(:,:,1));
    norm_sub_total_G_2 = double(norm_sub_total_2(:,:,2));
    norm_sub_total_B_2 = double(norm_sub_total_2(:,:,3));
    thin_sub_total_G_avg_2(i) = mean(norm_sub_total_G_2(BW_DNAPL_thin));
        
    norm_sub_total_Win_2 = norm_sub_total_2;
    norm_sub_total_Win_2(~BW_DNAPL_3D) = 0;
    
    % figure; imagesc(norm_sub_total); impixelinfo; title('norm sub total 2');
    % figure; imagesc(norm_sub_total_Win_2); impixelinfo; title('norm sub total win 2');
    % figure; imagesc(OEI_roi); impixelinfo; title('reg OEI')
    % figure; imagesc(OEI_Win); impixelinfo; title('OEI_Win');
    % figure; imagesc(OEI_sub_Win); impixelinfo; title('OEI sub Win');
    % figure; imagesc(norm_sub_Win); impixelinfo; title('norm_sub_Win');
    % figure; imagesc(norm_OEI); impixelinfo; title('norm_Win');
    
    % permanganate zone
    perm_zone = (OEI_G < 45000) & (norm_sub_total_G > 2000) & (~gs_noise);
    perm_zone = bwareaopen(perm_zone, 1000, 8);
    perm_zone = imfill(perm_zone, 'hole') & (~BW_DNAPL); 
    % figure; imagesc(perm_zone); impixleinfo
    
    % aqueous phase avg
    BW_aq = (~BW_DNAPL) & (~gs_noise) & perm_zone;
    % figure; imagesc(BW_aq); impixelinfo
    aq_avg_R(i) = mean(OEI_R(BW_aq));
    aq_avg_G(i) = mean(OEI_G(BW_aq));
    aq_avg_B(i) = mean(OEI_B(BW_aq));
    
    % diff_MnO2_zone
    switch n
        case {40,60,80,100,120,140,160,180,200,220,240,260,280,300,320,340,360,380,400,420,440,460,480,500,520,540,560,580,600,620,640}
            diff_MnO2_thresh_G(i) = 10000;  
            total_thresh_G(i) = min(12000 + thin_sub_total_G_avg_2(i), 14000);
            diff_MnO2_zone = ((norm_sub_G > diff_MnO2_thresh_G(i))...
                | (norm_sub_B > 1.1 * diff_MnO2_thresh_G(i))) & (BW_DNAPL);
            total_diff_zone = ((norm_sub_total_G_2 > total_thresh_G(i))...
                | (norm_sub_total_B_2 > 12000)) & (BW_DNAPL);
            new_diff_zone = ((OEI_G < 12000 & OEI_B < 4000) | (OEI_G < 16000 & OEI_B < 1000)) & BW_edge;
            new_diff_zone = new_diff_zone & ~gs_noise;
            new_diff_zone = new_diff_zone | ((OEI_G < 10000 & OEI_B < 600) & gs_noise);
            BW_diff = diff_MnO2_zone | total_diff_zone;
            BW_diff_neg = BW_diff & (OEI_G > aq_avg_G(i) - 3000 & OEI_B > aq_avg_B(i) - 3000); 
            BW_diff_neg = imfill(BW_diff_neg, 'holes');
            
        otherwise 
            diff_MnO2_thresh_G(i) = 10000; 
            total_thresh_G(i) = min(12000 + thin_sub_total_G_avg_2(i), 14500);
            diff_MnO2_zone = ((norm_sub_G > diff_MnO2_thresh_G(i))...
                | (norm_sub_B > 1.1 * diff_MnO2_thresh_G(i) )) & (BW_DNAPL);
            total_diff_zone = ((norm_sub_total_G_2 > total_thresh_G(i))...
                | (norm_sub_total_B_2 > 12000)) & (BW_DNAPL);
            new_diff_zone = ((OEI_G < 10000 & OEI_B < 3000) | (OEI_G < 12000 & OEI_B < 800) | (OEI_G < 19000 & OEI_B < 500)) & BW_edge;
            new_diff_zone = new_diff_zone & ~gs_noise; 
            new_diff_zone = new_diff_zone | ((OEI_G < 10000 & OEI_B < 600) & gs_noise);
            BW_diff = diff_MnO2_zone | total_diff_zone;
            BW_diff_neg = BW_diff & (OEI_G > aq_avg_G(i) & OEI_B > max(aq_avg_B(i),1000));
            BW_diff_neg = imfill(BW_diff_neg, 'hole');
    end
    % figure; imagesc(BW_diff_neg); impixelinfo
    diff_MnO2_size(i) = sum(diff_MnO2_zone(:));
    total_diff_size(i) = sum(total_diff_zone(:));
    
    BW_diff = (BW_diff & ~BW_diff_neg) | new_diff_zone;
    BW_DNAPL_new = BW_DNAPL;
    BW_DNAPL_new(BW_diff) = 0;
    BW_DNAPL_new = imfill(BW_DNAPL_new,'hole') & ~iTCE_hole;
    BW_DNAPL_small = BW_DNAPL_new  & ~bwareaopen(BW_DNAPL_new, 25, 4);
    BW_DNAPL_small_neg = BW_DNAPL_small & (OEI_G < 15500 & OEI_B < 4000);
    BW_DNAPL_new = BW_DNAPL_new & ~BW_DNAPL_small_neg;
    BW_diff = BW_DNAPL & ~BW_DNAPL_new;
    BW_Diff_new = BW_diff;
  
    BW_diff_th(BW_Diff_new) = time/60;
    Diff_size(i) = sum(BW_Diff_new(:));
    % figure; imagesc(BW_diff_th); impixelinfo; title('BW TCE new all'); 
    % figure; imagesc(BW_Diff_new); impixelinfo; title('BW_Diff_new'); colormap(gray)

    BW_DNAPL = BW_DNAPL_new;
    BW_DNAPL_3D = [BW_DNAPL, BW_DNAPL, BW_DNAPL];
    
    OEI_Win_new  = OEI_roi;
    OEI_Win_new(~BW_DNAPL_3D) = 0;
    % figure; imagesc(OEI_Win_new); impixelinfo; title('OEI_Win_new');
    % figure; imagesc(OEI_roi); impixelinfo; title('reg OEI')
    
    BW_diff_neg_3D = [BW_diff_neg, BW_diff_neg, BW_diff_neg];
    OEI_Win_new_2 = OEI_Win_new;
    OEI_Win_new_2(BW_diff_neg_3D) = 0;
    % figure; imagesc(OEI_Win_new_2); impixelinfo; title('OEI_Win_new_2');
    
    DNAPL_Size_vec(i) = sum(BW_DNAPL(:));
    V_Sat_vec(i) = sum(A_F_roi(BW_DNAPL)) / sum(A_F_roi(:));
    
    interfacial_edge = bw_edge(imfill(BW_DNAPL,'holes'),1);
    interfacial_area(i) = sum(A_F_roi(interfacial_edge)) * pixel_size; % cm^2 
    file = char(fl(n));
    
    file_dir = '/karst/scratch/RESTORE/karst/data/Masoud/Masoud_Exp/exp_04/09_KMnO4_injection/Binary_DNAPL/';
    file_name_bw = [file(1:end-5),'_bw.mat'];
    NewFile_bw = [file_dir, file_name_bw];
    % now save the image in the new directory
    save(NewFile_bw,'BW_DNAPL','OEI_roi','BW_diff_neg','time')

end
Sat_vec = DNAPL_Size_vec/dnaplsize;
V_DNAPL_vec = V_Sat_vec * V_cell; % cm^3, ml
file_dir = '/karst/scratch/RESTORE/karst/data/Masoud/Masoud_Exp/exp_04/';
file_name_sat = ['Saturation_exp_04_main.mat'];
NewFile_sat = [file_dir, file_name_sat];
% now save the image in the new directory
save(NewFile_sat,'Sat_vec','V_DNAPL_vec','interfacial_area','time_new','BW_diff_th','iTCE','V_cell')%, 'gas_total_sat', 'BW_gas_th')

% figure; plot(time_new/60, V_DNAPL_vec)
% figure; plot(time_new/60, interfacial_area); title('I_A (cm^2)')
% figure; plot(time_new/60, interfacial_area./V_DNAPL_vec); title('I_A over Volume (cm^2 / cm^3)')
toc