% function RootSkel(mpd)
%To extract root system skeleton through automatic image processing

%Load dataset lists
if ~exist('mpd')
    mpd=uigetdir(pwd,['Select the main project folder']);
end
load([mpd '\dataset.mat']);

for s=1:size(imgpth,1);
    c_imgpth=imgpth(s,:);
    c_resfol=resfol(s,:);
    c_data=dataset(s,:);
    if~exist([c_resfol '\skel.mat'])
        
        disp(['Processing ' c_data '...']);
        im=imread(c_imgpth);
        
        %% Convert to hsv
        i2=rgb2hsv(im);
        v=i2(:,:,3);
        % figure; imshow(v); colorbar; colormap jet; caxis([0 .5]); set(gcf,'NumberTitle','off','Name','Value')
        
        %% Remove dark background
        mask=zeros(size(v));
        mask(v>0.3)=1;
        % figure; imshow(mask); colorbar; set(gcf,'NumberTitle','off','Name','Mask')
        
        %% Remove specks in background
        mask=imerode(mask,ones(1));
        % figure; imshow(mask); colorbar; set(gcf,'NumberTitle','off','Name','Mask')
        
        %% Remove debris etc
        cc=bwconncomp(mask);
        numobj=cc.NumObjects;
        objsize=cellfun(@numel,cc.PixelIdxList);
        for k=1:numobj
            if objsize(k)<prod(size(mask))*.00002 %| objsize(k)==max(objsize)
                mask(cc.PixelIdxList{k})=0;
            end
        end
        % figure; imshow(mask); colorbar; set(gcf,'NumberTitle','off','Name','Mask')
        
        %% Dilate slightly to connect any roots segmented by darker areas
        mask=imdilate(mask,ones(2));
        % figure; imshow(mask2); colorbar; set(gcf,'NumberTitle','off','Name','Dilated')
        
        %% Remove small unconnected objects
        cc=bwconncomp(mask);
        numobj=cc.NumObjects;
        objsize=cellfun(@numel,cc.PixelIdxList);
        for k=1:numobj
            if objsize(k)< max(objsize)/4
                mask(cc.PixelIdxList{k})=0;
            end
        end
        % figure; imshow(mask2); colorbar; set(gcf,'NumberTitle','off','Name','Mask')
        
        %% Smooth edges
        mask=bwmorph(mask,'majority',15);
        % figure; imshow(mask); colorbar; set(gcf,'NumberTitle','off','Name','Smooth'); figure(gcf)
        
        % %% Another step to smooth?
        % se = strel('disk',3);
        % erodedBW = imdilate(mask,se);
        % figure, imshow(erodedBW)
        % erodedBW = imerode(mask,se);
        % figure, imshow(erodedBW)
        
        %% Thin to skeleton
        skel=bwmorph(mask,'thin','Inf');
        % figure; imshow(skel); colorbar; set(gcf,'NumberTitle','off','Name','Thin'); figure(gcf)
        %             [yall xall]=find(skel==1);
        
        %% Remove small unconnected objects
        cc=bwconncomp(skel);
        numobj=cc.NumObjects;
        objsize=cellfun(@numel,cc.PixelIdxList);
        for k=1:numobj
            if objsize(k)< max(objsize)/4
                skel(cc.PixelIdxList{k})=0;
            end
        end
        
        [ally allx]=find(skel==1);
        %     figure
        %     warning('off','images:imshow:magnificationMustBeFitForDockedFigure')
        %     warning('off', 'Images:initSize:adjustingMag');
        %     imshow(im); hold on
        %     plot(allx,ally,'.g')
        
        %%Save mask and skeleton of image
        save([c_resfol '/skel.mat'],'mask','skel','allx','ally');
    end
end