%% Fuzzy Logic Detection of Varves
%
% Ebert, T., Trauth, M.H., Automated detection of annual layers (varves) in
% lake sediments using a Fuzzy Logic algorithm
%
% Version 17th July 2014
clear all, clc, close all

% p1 = '/Users/suguta/Desktop/ebert_matlabcode';
% cd(p1)
%%
% Read image data of thin section

I_LPL = imread('ebert_image_lpl.tif');    % Image with plain polarized light.
I_XPL = imread('ebert_image_xpl.tif');    % Image with crossed polarized light.
I_BG  = imread('ebert_image_bg.tif');     % Image mask for background pixel.

figure(1)

subplot(2,1,1), imshow(I_LPL), title('image{\_}lpl'),
text(5500,440,'a)','FontSize',12,'FontWeight','bold','BackGround','w')

subplot(2,1,2), imshow(imrotate(I_XPL,180,'bicubic')), title('image{\_}xpl')
text(5500,440,'b)','FontSize',12,'FontWeight','bold',...
    'Color','k','BackGround','w')

%%
% Despeckling the image

IR1 = I_XPL(:,:,1); IG1 = I_XPL(:,:,2); IB1 = I_XPL(:,:,3);

IR2 = medfilt2(IR1,[15 15],'symmetric');
IG2 = medfilt2(IG1,[15 15],'symmetric');
IB2 = medfilt2(IB1,[15 15],'symmetric');

IR3 = imsubtract(IR1,IR2);
IG3 = imsubtract(IG1,IG2);
IB3 = imsubtract(IB1,IB2);

IR4 = imsubtract(IR2,IR1);
IG4 = imsubtract(IG2,IG1);
IB4 = imsubtract(IB2,IB1);

IR5 = IR1; IG5 = IG1; IB5 = IB1;

IR5(IR3 > 20 | IR4 > 20) = IR2(IR3 > 20 | IR4 > 20);
IG5(IG3 > 20 | IG4 > 20) = IG2(IG3 > 20 | IG4 > 20);
IB5(IB3 > 20 | IB4 > 20) = IB2(IB3 > 20 | IB4 > 20);

I_DESP = cat(3,IR5,IG5,IB5);

figure(2),imshow(imrotate(I_DESP,180,'bicubic')), title('I-DESP')

%%
% Top-hat filtering

for i = 1:3
    I_BCKGD(:,:,i) = imtophat(I_DESP(:,:,i),strel('disk',15));
end

figure(3), imshow(imrotate(I_BCKGD,180,'bicubic'))
title('background-corrected image: I-BCKGD')

%%
% Load data of training areas

load('ebert_polygons_1.mat')    % Summer
load('ebert_polygons_2.mat')    % Winter

figure(4)
imshow(I_DESP), title('I-DESP + associated training areas')
axis on, hold on

figure(4)
for i = 1:size(xiS,2)
      patch('XData',xiS(i).dat,'YData',yiS(i).dat,'FaceColor',...
      	[.93 .60 .30],'Edgecolor','none')
end

for i = 1:size(xiW,2)
      patch('XData',xiW(i).dat,'YData',yiW(i).dat,'FaceColor',...
      	[.49 .19 .09],'Edgecolor','none')
end

clear x* y*                                  

%%
% Remove background pixels

R = double(I_BCKGD(:,:,1)); R(I_BG==0) = NaN;
G = double(I_BCKGD(:,:,2)); G(I_BG==0) = NaN;
B = double(I_BCKGD(:,:,3)); B(I_BG==0) = NaN;

%%
% Select coordinates of pixels to estimate statistic values:

for i = 1:size(ROIS,2)
	  [XS(i).dat,YS(i).dat] = find(ROIS(i).dat ~= 0);
end

XS = cat(1,XS(:).dat); YS = cat(1,YS(:).dat);

for i = 1:size(ROIW,2)
       [XW(i).dat,YW(i).dat] = find(ROIW(i).dat ~= 0);
end

clear ROI*                                   

XW = cat(1, XW(:).dat); YW = cat(1, YW(:).dat);  

IVSR = zeros(1,length(XS));
IVSG = zeros(1,length(XS));
IVSB = zeros(1,length(XS));
for i = 1 : length(XS)
    IVSR(i) = R(XS(i),YS(i)); 
    IVSG(i) = G(XS(i),YS(i)); 
    IVSB(i) = B(XS(i),YS(i)); 
end

IVWR = zeros(1,length(XW));
IVWG = zeros(1,length(XW));
IVWB = zeros(1,length(XW));
for i = 1 : length(XW)
    IVWR(i) = R(XW(i),YW(i)); 
    IVWG(i) = G(XW(i),YW(i)); 
    IVWB(i) = B(XW(i),YW(i)); 
end

MEAN = [[nanmean(IVSR),nanmean(IVSG),nanmean(IVSB)];...
    [nanmean(IVWR),nanmean(IVWG),nanmean(IVWB)]];

SIGMA = [[nanstd(IVSR),nanstd(IVSG),nanstd(IVSB)];...
    [nanstd(IVWR),nanstd(IVWG),nanstd(IVWB)]];

% clear IV* X* Y*  

%%
% Summer == 0; Winter == 1;

WR = [XW YW  IVWR'];
WG = [XW YW  IVWG'];
WB = [XW YW  IVWB'];

SR = [XS YS  IVSR'];
SG = [XS YS  IVSG'];
SB = [XS YS  IVSB'];

W = [IVWR' IVWG' IVWB' zeros(size(IVWR',1),1)];
W = rmmissing(W);

S = [IVSR' IVSG' IVSB' ones(size(IVSR',1),1)];
S = rmmissing(S);

data = cat(1,W,S);
dataRange = [min(data)' max(data)'];

%%
% Create separate input and output data sets, X and Y, respectively

X = data(:,1:3); 
Y = data(:,4);

%%
% Partition the input and output data sets into training data 
% (odd-indexed samples) and validation data (even-indexed samples).

trnX = X(1:2:end,:); % Training input data set
trnY = Y(1:2:end,:); % Training output data set
vldX = X(2:2:end,:); % Validation input data set
vldY = Y(2:2:end,:); % Validation output data set

%%
% Extract the range of each data attribute, which you will use for
% input/output range definition during FIS construction.

dataRange = [min(data)' max(data)'];

%%
% Rank Inputs According to Correlation Coefficients

% c1 = corrcoef(data);
% c1(end,:)
% [B,rankIDX] = sort(c1(end,1:end-1),'descend'); % [B,rankIDX] = sort(c1(end,1:end-1),'ascend')
rankIDX = [1 2 3];

%%
% Create Fuzzy Inference System fis1

fis1 = sugfis('Name','fis1');
fis1 = addInput(fis1,dataRange(1,:),'NumMFs',2,...
    'MFType','trimf','Name',"Red");
fis1 = addInput(fis1,dataRange(2,:),'NumMFs',2,...
    'MFType','trimf','Name',"Green");
fis1 = addInput(fis1,dataRange(3,:),'NumMFs',2,...
    'MFType','trimf','Name',"Blue");
fis1 = addOutput(fis1,dataRange(4,:),'NumMFs',6);

%%
% Specify tuning options to learning rules

options = tunefisOptions('Method','particleswarm',...
    'OptimizationType','learning', ...
    'NumMaxRules',6,...
    'UseParallel',1);
options.MethodOptions.MaxIterations = 50;

%%
% Initialize the random number generator 

rng('default')

%%
% Set the input order of the training data

inputOrders1 = rankIDX; % inputOrders1 = [1 2 3];
orderedTrnX1 = trnX(:,inputOrders1);

%%
% Learning rules with tunefis function

fisTout1 = tunefis(fis1,[],orderedTrnX1,trnY,options);
fprintf('Training RMSE = %.3f SW\n',calculateRMSE(fisTout1,orderedTrnX1,trnY));
fprintf("Total number of rules = %d\n",numel([fisTout1.Rules]));

%%
% Plot membership functions for input variables

figure(5)
hold on

for i = 1:3
    subplot(3,1,i)
    plotmf(fisTout1,'input',i)
end

subplot(3,1,1), title({'Shape of membership functions:';...
    'learning mode rules'})

%% 
% Tune the input/output MF parameters of the learned rules

[in,out,rule] = getTunableSettings(fisTout1);
options.OptimizationType = 'tuning';
options.Method = 'patternsearch'; % use a local optimization method
% options.KFoldValue = 3;
% options.ValidationWindowSize = 20;
% options.ValidationTolerance = 0.05;
options.MethodOptions.MaxIterations = 25; % Number of iterations
options.UseParallel = 1; % Improve speed of the tuning process

rng('default')

%%
% Tuning parameter values with tunefis function

fisTout2 = tunefis(fisTout1,[in;out;rule],orderedTrnX1,trnY,options);
fprintf('Training RMSE = %.3f SW\n',calculateRMSE(fisTout2,orderedTrnX1,trnY));

%%
% plot membership functions for input variables

figure(6)
for i = 1:3
    subplot(3,1,i)
    plotmf(fisTout2,'input',i)
end

subplot(3,1,1), title({'Shape of membership functions:';...
    'tuning mode rules'})

%%
% Fuzzy inference calculation

InputX = [R(:),G(:),B(:)];
inputOrders1 = rankIDX;
orderedInputX = InputX(:,inputOrders1);

tic
fisOut= reshape(evaluateFIS(fisTout2,orderedInputX),size(R));
toc

fisOut(fisOut < 0.5) = 0; fisOut(fisOut > 0.5) = 1;

fprintf("Total number of undefined numeric results: NaN = %d\n",...
    sum(isnan(fisOut(:))));

% figure
% imshow(fisOut)

%%
% Display fuzzy inference output

IOut = mat2gray(fisOut,[0 1]);
IOut(I_BG == 0) = NaN;
IOut = im2uint8(repmat(IOut,[1 1 3]));
figure(7), subplot(2,2,1),imshow(IOut)
title('fisOut'), axis([2800 4000 600 1400])

IoutR = IOut(:,:,1);IoutG = IOut(:,:,2);IoutB = IOut(:,:,3);
IoutR(I_BG==0) = 0;
IoutG(I_BG==0) = 0;
IoutB(I_BG==0) = 255;
IOut = cat(3,IoutR,IoutG,IoutB);
figure(7), subplot(2,2,2), imshow(IOut)
title('fisOut (back ground in blue)'), axis([2800 4000 600 1400])

IoutR = IOut(:,:,1);IoutG = IOut(:,:,2);IoutB = IOut(:,:,3);
IoutR(fisOut==0.5 & I_BG~=0) = 255;
IoutG(fisOut==0.5 & I_BG~=0) = 0;
IoutB(fisOut==0.5 & I_BG~=0) = 0;
IOut = cat(3,IoutR,IoutG,IoutB);
figure(7),subplot(2,2,3), imshow(IOut)
title('fisOut (gc = 0.5 in red)'), axis([2800 4000 600 1400])

IoutR(isnan(fisOut) & I_BG~=0) = 0;
IoutG(isnan(fisOut) & I_BG~=0) = 255;
IoutB(isnan(fisOut) & I_BG~=0) = 0;
IOut = cat(3,IoutR,IoutG,IoutB);
figure(7), subplot(2,2,4), imshow(IOut)
title('fisOut (gc = NaN in green)'), axis([2800 4000 600 1400])

clearvars I* -except I_* IOut                
delete(gcp)
%% ----------------- ANFIS-F algorithm: -----------------------------------

% Define subimages for cell arrays to save memory:

ai = 225; bi = 364;
M = size(fisOut,1)/ai; N = size(fisOut,2)/bi;
m = ai*ones(1,M); n = bi*ones(1,N);

%
fisOut = mat2cell(fisOut,m,n);
I_BG = mat2cell(I_BG,m,n);                          % I_DSP(I_BG==0) = NaN;

clear a* b* m n                              % Clear variables from memory!

% eps = 1e-6;                         % Stop creterion for while condition.
eps = 0.001;

anfisOut = fisOut;                   % Initialization of variable anfisOut.

% Delaunay Triangulation and ANFIS-F Calculation:
tic
while size(find(~(cell2mat(anfisOut)==0 | cell2mat(anfisOut)==1 | ...
    cell2mat(I_BG)==0)),1)*100/(size(R,1)*size(R,2)) > eps
    for i = 1 : M
        for j = 1 : N
            newfisOut = padarray(cell2mat(anfisOut(i,j)),[5 5],...
                'replicate','both');
            newIBG = padarray(cell2mat(I_BG(i,j)),[5 5],...
                'replicate','both');
            [yc,xc] = find(~(newfisOut==0 | newfisOut==1 | newIBG==0));
            [yu,xu] = find((newfisOut==0 | newfisOut==1 | newIBG==1));
            if isempty(xc)
                anfisOut(i,j) = {newfisOut(6:end-5,6:end-5)};
            elseif ~isempty(xc)
                TRI = delaunayTriangulation(xu,yu);% delaunay triangulation
                k = randperm(length(xc));    % or k = randi(size(xc,1),1,1)
                for i2 = 1:length(k)
                    PL = pointLocation(TRI,xc(k(i2)),yc(k(i2)));
                    T = TRI.ConnectivityList(PL,:);
                    P = TRI.Points(TRI.ConnectivityList(PL,:),:);
                    dt = delaunayTriangulation(P);
                    Outp = impixel(newfisOut,P(:,1),P(:,2));
                    Output = Outp(:,1);
                    if var(Output) == 0
                        newfisOut(yc(k(i2)),xc(k(i2))) = Output(1,1);
                        anfisOut(i,j) = {newfisOut(6:end-5,6:end-5)};
                    elseif var(Output) ~= 0
                        Data = [dt.Points(:,1) dt.Points(:,2) Output];
%                         Data = [P(:,1) P(:,2) Output];
                        trnData = Data(1:2:end,:);
                        chkData = Data(2:2:end,:);
                        inFIS = genfis(trnData(:,1:end-1),trnData(:,end),...
                            genfisOptions('GridPartition'));
                        options = anfisOptions('InitialFIS',inFIS,...
                            'ValidationData',chkData,...
                            'EpochNumber',25,...
                            'DisplayANFISInformation',0,...
                            'DisplayErrorValues',0,...
                            'DisplayStepSize',0,...
                            'DisplayFinalResults',0);
                        [fis1,error1,ss,fis2,error2] = ...
                            anfis(trnData,options);
                        newfisOut(yc(k(i2)),xc(k(i2))) = ...
                            round(evalfis(fis1,[xc(k(i2)) yc(k(i2))]));
                        anfisOut(i,j) = {newfisOut(6:end-5,6:end-5)};
                    end
                end
            end
        end
    end
end
toc
anfisOut = cell2mat(anfisOut); 

%%
% Display ANFIS-F output

AOut = mat2gray(anfisOut,[0 1]);
AOut = im2uint8(repmat(AOut,[1 1 3]));

AoutR = AOut(:,:,1);AoutG = AOut(:,:,2);AoutB = AOut(:,:,3);
AoutR(cell2mat(I_BG)==0) = 0;
AoutG(cell2mat(I_BG)==0) = 0;
AoutB(cell2mat(I_BG)==0) = 255;
AOut = cat(3,AoutR,AoutG,AoutB);

AoutR = AOut(:,:,1);AoutG = AOut(:,:,2);AoutB = AOut(:,:,3);
AoutR(anfisOut == 0.5 & cell2mat(I_BG)~=0) = 255;
AoutG( anfisOut == 0.5 & cell2mat(I_BG)~=0) = 0;
AoutB(anfisOut == 0.5 & cell2mat(I_BG)~=0) = 0;
AOut = cat(3,AoutR,AoutG,AoutB);

figure(8), 
subplot(1,2,1), imshow(IOut)
title('fisOut (gc = 0.5 in red)'), axis([2800 4000 600 1400])
subplot(1,2,2), imshow(AOut)
title('anfisOut (gc = 0.5 in red)'), axis([2800 4000 600 1400])

clear A*                                     

%%
% Iteratively despeckling algorithm

I2 = anfisOut; I6 = anfisOut;

m = 3; n = 3;

I3 = medfilt2(I2,[m n],'symmetric');

I4 = imsubtract(I2,I3);
I5 = imsubtract(I3,I2);

tempSize = [0 size(I2(:),1)]';

while tempSize(end)~=tempSize(end-1)
    I3 = medfilt2(I2,[m n],'symmetric');
    I4 = imsubtract(I2,I3);
    I5 = imsubtract(I3,I2);
    I6(I4>0 | I5>0) = I3(I4>0 | I5>0);
    I2 = round(I6);
    tempSize = cat(1,tempSize,size(find(I4~=0 & I5~=0),1));
end
dspSig = I2;

clearvars I* m n t* -except I_*              

%%
% Display despeckling result

ITemp = reshape(im2uint8(repmat(mat2gray(dspSig,...
    [0 1]),[1 1 3])),[],3);

colorCode = [0 0 255];
for i = 1:3
    ITemp(reshape(cell2mat(I_BG),[],1)==0,i) = colorCode(i);
end
ITemp = reshape(ITemp,[size(dspSig) 3]);

figure(9)
subplot(2,1,1), imshow(ITemp)
title('Despeckled anfisOut: dspSig'), axis([2800 4000 600 1400])

clearvars I* i c* -except I_*                

% medf2Sig = round(medfilt2(round(dspSig),[60 3],'symmetric'));
medf2Sig = round(medfilt2(round(dspSig),[35 1],'symmetric'));
medf2Sig(cell2mat(I_BG) == 0) = NaN;

IMF = reshape(im2uint8(repmat(mat2gray(medf2Sig,[0 1]),[1 1 3])),[],3);
colorCode = [0 0 255];
for i = 1:3
    IMF(reshape(cell2mat(I_BG),[],1)==0,i) = colorCode(i);
end
IMF = reshape(IMF,[size(dspSig) 3]);

figure(9)
subplot(2,1,2),imshow(IMF)
title('medf2Sig'), axis([2800 4000 600 1400])

clear i c*                                   

%%
% Count and measure annual layers

x = 600:5500; y = 600:100:1300;

for i = 1:length(y)
    [cx(:,i) cy(:,i) c(:,i)] = ...
        improfile(medf2Sig,[x(1) x(end)],[y(i) y(i)],'nearest');
    cx(:,i) = round(cx(:,i)); c(:,i) = round(c(:,i));
end

figure(10)
subplot(2,1,1),imshow(IMF), hold on, title('medf2Sig')
subplot(2,1,2), imshow(I_XPL), hold on, title('I-XPL')
for j = 1:2
    for i = 1:length(y)
        subplot(2,1,j); x1 = [min(x) max(x)]; y1(i,:) = [cy(1,i) cy(1,i)];
        plot(x1,y1(i,:),'-y',cx(:,i),(c(:,i)*20)+(cy(1,i)-20/2),'-r');
    end
end    

trackCode = [1200 678 836;800 874 1021;900 1017 1217;600 1247 1359;...
    600 1442 1607;1300 1548 1892;1200 1896 2073;1300 2085 2259;...
    1200 2256 2417;1100 2425 2480;900 2467 2728;1000 2742 2775;...
    600 2775 2894;1000 2919 3224;900 3222 3563;1200 3573 3681;...
    1100 3682 4214;700 4241 4413;1200 4427 5107];

for i = 1:size(trackCode,1)
    IDX(:,i) = knnsearch([cx(:,1)],trackCode(i,2:3)');
    IDY(i)   = knnsearch(cy(1,:)',trackCode(i,1));
end

for i = 1:size(trackCode,1)
temp = [cx(IDX(1,i):IDX(2,i));c(IDX(1,i):IDX(2,i),IDY(i))'];
trackData{i} = temp';
end

savefile = 'trackData';
save(savefile,'trackData');

figure(10)
for j = 1:2
    for i = 1:size(trackCode,1)
        subplot(2,1,j)
        plot(cx(IDX(1,i):IDX(2,i)),(c(IDX(1,i):IDX(2,i),IDY(i))*20)+...
            (cy(1,IDY(i))-20/2),'-g');
    end
end

clearvars I* c* i j s* t* x* y* -except I_* IMF trackData

%%
% Time-series analysis

xb = [];
for i = 1:size(trackData,2)
    xb = cat(1,xb,trackData{i}(:,2));
end
xb(isnan(xb)==1) = [];

xb = reshape(xb,size(xb,1)*size(xb,2),1);
xb2(2:length(xb)+1) = xb; 

if xb(1) == 0 & xb(end) == 1;
    xb2(1) = 1; xb2(end+1) = 0;
end
if xb(1) == 1 & xb(end) == 0;
    xb2(1) = 0; xb2(end+1) = 1;
end
    
transition = diff(xb2); 
x0 = zeros(size(transition)); x1 = zeros(size(transition));
k0 = find(transition == 1); x0(k0) = k0;
k1 = find(transition == -1); x1(k1) = k1;

if x0(k0(1)) < x1(k1(1)) 
    k2 = find(x0 ~= 0); PixNum = diff(x0(k2));
else
    k2 = find(x1 ~= 0); PixNum = diff(x1(k2));
end

t = [1:length(PixNum)];
x = PixNum*21*10^-3;

clearvars i k* t* x* P* -except xb x t
