% This function is an implementation of Mao Li's work [Li et al. (2018)]. The matrices E and F are requiered, and can be 
% downloaded from https://github.com/maoli0923/Persistent-Homology-All-Leaf. The function ksdensity2d.m,
% which is also required, can be found in the same website.


% Original work: Li M, An H, Angelovici R, Bagaza C, ... (2018) Topological Data Analysis as a Morphometric
% Method: Using Persistent Homology to Demarcate a Leaf Morphospace. Front. Plant Sci. 9:553. doi: 10.3389/fpls.2018.00553


% Here are the parameters: 

	% filename: a picture with the objects to measure
	% resizeIndex: when having a picture with very large resolution, reducing picture size helps to increase speed. (values from >0 and 1)
	% minArea: minimum area (in pixels) of an object in order to be analyzed
	% refNumber: number of references. The number of references must be the same on both sides of the picture
	% ncolmin and ncolmax: delimit pixel columns to be removed from the picture
	% boundResize: this function produces contours for all the objects, and this parameters can be used to resize them (valus from >0 and 1)
	% returnImage: true/false, returns contours
	% fillHoles: applies matlab function bwareaopen to fill holes in the image
	% colorThr: threshold for background segmentation in blue channel (values between 0 and 255)
	% the rest of the parameters are described in Li et al. 2017.





function superEC=superfireringer(filename,resizeIndex,minArea,refNum,ncolmin,ncolmax,boundResize,kN,U,sigma,sigma0,n1,n2,threshold,step,E,Face,hgrid,vgrid,returnImage,fillHoles,colorThr)

originalImage=imresize(imread(filename),resizeIndex);                      


%image segmentation based on layer3

imlayer3=originalImage(:,:,3);
f=find(imlayer3<colorThr);
mark1=zeros(size(imlayer3));
mark1(f)=1;

% removing small objects

if fillHoles
    binaryImage = imfill(imclearborder(bwareaopen(mark1,minArea)),'holes');
else
    binaryImage = imclearborder(bwareaopen(mark1,minArea));
end

labeledImage = logical(bwlabel(binaryImage, 8));
labeledImage(:,ncolmin:ncolmax)=0;
dat1 = regionprops(labeledImage, binaryImage, 'all');


% get area of the references

if refNum>1
    idxRefs=[1:refNum (length(dat1)-(refNum-1)):length(dat1)];
    refArea=[median([dat1(idxRefs).Area]) std([dat1(idxRefs).Area])];
    
dat2=dat1;
dat2(idxRefs)=[];

else
    dat2=dat1;
    refArea=[0 0]; 
end



clear allIm
clear allIm2
for j=1:length(dat2)
    mn=size(dat2(j).Image);
    rat=mn(1)/mn(2); 
    allIm(j).im=imresize(dat2(j).Image,[boundResize boundResize/rat]);             % resizing and rotating image
    allIm2(j).bound=cell2mat(bwboundaries(dat2(j).Image));
    allIm2(j).orientation=dat2(j).Orientation;
    allIm2(j).gps=dat2(j).BoundingBox(1:2);
    if returnImage
        allIm2(j).binary=dat2(j).Image;
    end
end


% running persistent homology

clear EC
clear FF
parfor ii=1:length(allIm)
    
    kai=fireringerX(allIm(ii).im,kN,U,sigma,sigma0,n1,n2,threshold,step,E,Face,hgrid,vgrid);
    
    EC(ii,:)=reshape(kai',1,numel(kai));
    
end

superEC.EC=EC;
superEC.id=filename;
superEC.im=allIm2;
superEC.refData=refArea;


dat1=dat2;


% extracting original coloxrs

thrIm=originalImage;
for i=1:3
    x=thrIm(:,:,i);
    x(find(labeledImage==0))=NaN;
    thrIm(:,:,i)=x;
end

for i=1:length(dat1)
    x0=ceil(dat1(i).BoundingBox(1));
    x1=x0+ceil(dat1(i).BoundingBox(3));
    y0=floor(dat1(i).BoundingBox(2));
    y1=y0+ceil(dat1(i).BoundingBox(4));
    myim=thrIm(y0:y1,x0:x1,:);
    
    clear colData
    for im=1:3
        xi=double(reshape(myim(:,:,im),1,prod(size(myim(:,:,im)))));
        xii=tabulate(xi);
        xiii=xii(:,1:2);
        colData(im).col=xiii;
    end
    colD(i).im=colData;
end

superEC.color=colD;



function kai=fireringerX(box,kN,U,sigma,sigma0,n1,n2,threshold,step,E,Face,hgrid,vgrid)

    h=ones(size(box));
    aa=find(box==0);
    h(aa)=0;                 % switch black and white pixel
    [start_r,start_c] = find(h,1,'first'); % get start point
    V = bwtraceboundary(h,[start_r start_c],'W',8,Inf,'counterclockwise'); % trace the contour
    V=V-repmat(mean(V),length(V),1); % translate the center to origin
    V=V/norm(V,'fro')*sqrt(length(V)); %normalize leaf contour by centroid size
    F = ksdensity2d([V(:,1),V(:,2)],hgrid,vgrid,[sigma0,sigma0]); % Gaussian density estimator, some other person's code
    %kai2=zeros(kN,step);
    for t=1:kN
        kai2=zeros(1,step);
        G=exp(-abs(sum((U-repmat([0 0],size(U,1),1)).^2,2).^0.5-t*sigma).^2/(2*sigma^2))/(2*sigma*pi);
        Z=reshape(reshape(F,size(G)).*G,size(F));
        ZZ=-reshape(Z,n1*n2,1);
        s=find(min(ZZ)>threshold);
        for ii=size(s,2)+1:step
            a=find(ZZ<=threshold(ii));
            b=find(sum(ismember(E,a),2)==2);
            c=find(sum(ismember(Face,a),2)==3);
            kai2(ii)=length(a)-length(b)+length(c);
        end
        
        kai(t,:)=kai2;
    end