front_main_multiple_adj[compilation]

close all; clear all; clc;
addpath('C:\Users\iptea\Documents\MATLAB\Nonreplica September_ROUGHNESS\Roughness_Amore_2\General');

sbj = {'01_KB','02_YJ','03_PU','04_HG','05_HR',...
    '06_HG','07_EJ','08_MA','09_NY','10_MR'};
part = {'forearm','AF_forearm','forehead','AF_forehead','backhand','AF_backhand'};

%%

sht_tbl = 'rgh_refparam.xlsx';
T = readtable(sht_tbl,'ReadRowNames',false);
med_list_bf = T{1:30,6};
gau_list_bf = T{1:30,7};
med_list_aft = T{1:30,15};
gau_list_aft = T{1:30,16};

med_arr = zeros(10,6);
gsigma_arr = zeros(10,6);

med_arr(:,1) = med_list_bf(1:10,1);
med_arr(:,2) = med_list_aft(1:10,1);
med_arr(:,3) = med_list_bf(11:20,1);
med_arr(:,4) = med_list_aft(11:20,1);
med_arr(:,5) = med_list_bf(21:30,1);
med_arr(:,6) = med_list_aft(21:30,1);

gsigma_arr(:,1) = gau_list_bf(1:10,1);
gsigma_arr(:,2) = gau_list_aft(1:10,1);
gsigma_arr(:,3) = gau_list_bf(11:20,1);
gsigma_arr(:,4) = gau_list_aft(11:20,1);
gsigma_arr(:,5) = gau_list_bf(21:30,1);
gsigma_arr(:,6) = gau_list_aft(21:30,1);

%% Attention:
% diff_idl should ne greater than diff_rl
% size of median filter neighborhood should be greater than desired radius

%%
data_pth = 'C:\Users\iptea\Documents\MATLAB\AMORE_ROUGHNESS';

sht_name = 'rgh_frnt.xlsx';
ttl = {'part','sbj','avg_rgh','rms_rgh','max_rgh'};
ttl_rng1 = 'B1:F1';
ttl_rng2 = 'H1:L1';
xlswrite(sht_name,ttl,ttl_rng1);
xlswrite(sht_name,ttl,ttl_rng2);

max_no = 500;
avg_rgh = zeros(1,max_no);
rms_rgh = zeros(1,max_no);
max_rgh = zeros(1,max_no);
counter = 0;
incr1 = 1;
incr2 = 2;
fprintf('Roughness calculation started:\n');

for part_no = 1:length(part)
    for sbj_no = 1:length(sbj)
        %fprintf('%i/%i\n', part_no*sbj_no);
        fld_pth = strcat(data_pth,'\',sbj{sbj_no},'\',part{part_no},'\');
        disp(fld_pth);
        for img_no = 1:max_no
            fprintf('.');
            if mod(img_no,100) == 0
                fprintf('%i\n',img_no);
            end
            img_name = strcat('Image',num2str(img_no-1,'%04d'),'.jpg');
            img = imread(strcat(fld_pth,img_name));
            % figure()
            % imshow(img);
            [H,W] = size(img);

            %% Region of interest
            l_crop = 130;
            t_crop = 20;
            r_crop = 0;
            b_crop = 300;

            [img_crp,H,W] = crop(img,H,W,l_crop,r_crop,t_crop,b_crop);
            % figure()
            % imshow(img_crop);
            %% Prefiltering
            gsigma = 6;
            img_gau = imgaussfilt(img_crp,gsigma);
            %% Finding ideal surface
%             med = 60; %90
            med = med_arr(sbj_no,part_no); %60
            diff_idl = 15;
            img_med = medfilt2(img_gau,[med*2 med*2],'symmetric');
            img_med = filt_diff(img_med,H,W,diff_idl);
            [~,ind_idl,~] = find_srf(img_med,H,W); %img_idl&avg_idl not used
            % figure()
            % imshow(img_idl);

            %% Flattening surface according to ideal line
            % Note: We put gaussian image directly
            mrg = 60;
            [img_flt,H,W] = flatten(img_crp,H,W,ind_idl,mrg);
%             figure()
%             imshow(img_flt,[]);

            %% Prefiltering flattened image
            gsigma = gsigma_arr(sbj_no,part_no); %60
%             gsigma = 10; %11
            img_gau = imgaussfilt(img_flt,gsigma);
            % figure()
            % imshow(img_flt,[]);
            %% Finding real surface
            diff_rl = 5;
            img_diff = filt_diff(img_gau,H,W,diff_rl);
            % figure()
            % imshow(img_flt,[]);
            [~,ind_rl,avg_rl] = find_srf(img_diff,H,W); %img_rl
%             figure()
%             imshow(img_rl,[])
            % hold on
            % plot(avg_rl)

            %% Calculating roughness
            [avg_rgh(img_no),rms_rgh(img_no),max_rgh(img_no)] = rgh(avg_rl,ind_rl);
        end
        mean_avg_rgh = mean(avg_rgh);
        mean_rms_rgh = mean(rms_rgh);
        mean_max_rgh = mean(max_rgh);


        %%

        A = [part(1,part_no),sbj(1,sbj_no),num2str(mean_avg_rgh),num2str(mean_rms_rgh),...
        num2str(mean_max_rgh)];

        if mod(part_no,2) == 1
            rw_bef = part_no - incr1;
            vrbl_rng = strcat('B',num2str(rw_bef*10+sbj_no+1),':','F',num2str(rw_bef*10+sbj_no+1));
        else
            rw_aft = part_no - incr2;
            vrbl_rng = strcat('H',num2str(rw_aft*10+sbj_no+1),':','L',num2str(rw_aft*10+sbj_no+1));
        end
        xlswrite(sht_name,A,vrbl_rng);
        counter = counter + 1;
        fprintf('\n%i out of %i\n', counter, length(part)*length(sbj));

    end

    if mod(part_no,2) == 1
        incr1 = incr1+1;
    else
        incr2 = incr2+1;
    end

end

fprintf('Process finished.\n');
% B = T{1:10,9};
%%


%% Status of roughness correctness:
T = readtable(sht_name,'ReadRowNames',false);
A = T{1:10,3};
B = T{1:10,9};
C = A > B;
xlswrite(sht_name,C,'N2:N11');
A = T{11:20,3};
B = T{11:20,9};
C = A > B;
xlswrite(sht_name,C,'N12:N21');
A = T{21:30,3};
B = T{21:30,9};
C = A > B;
xlswrite(sht_name,C,'N22:N31');


%%