% function LatRoots(mpd)
%To extract lateral roots 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_data=dataset(s,:);
    c_imgpth=[mpd imgpth(s,end-35:end)];
    c_resfol=[mpd resfol(s,end-19:end)];
    if str2num(c_data(end))>0
        if ~exist([c_resfol '/latroots.mat'])
            disp(['Processing ' c_data '...']);
            load([c_resfol '/skel.mat'],'skel');
            load([c_resfol '/primary.mat'],'prirc','pripath', 'topx', 'topy');
            
            %% Remove other planes
            pskel=skel; % skeleton of the plane only
            for j=1:length(prirc)
                pskel(1:prirc(j,1)-1,prirc(j,2))=0;
            end
            pskel(:,1:min(prirc(:,2)-1))=0;
            pskel(:,max(prirc(:,2))+1:end)=0;
            
            %% Find lat root attachment points');
            bptsm = bwmorph(pskel,'branchpoints');
            %         pripathepts=pripath+eptsm;
            %         mask=zeros(size(pripath));
            %         mask(prirc(1,1),prirc(1,2))=1;
            %         Dpri=bwdistgeodesic(logical(pripathepts), logical(mask));
            attm=bptsm;
            attm(pripath==0)=0;
            [attpy attpx]=find(attm==1);
            
            %% Trace out lateral roots
            
            n=0; natt=[]; attx=[]; atty=[]; %tipx=[]; tipy=[]; distpri=[]; rootx={}; rooty={};
            latskel=zeros(size(pskel)); %Record map of coordinates of lateral roots
            pskel=pskel-pripath+attm;
            cc=bwconncomp(pskel,8);
            for k=1:cc.NumObjects
                %Discard if not attached to the primary
                if sum(attm(cc.PixelIdxList{k}))<1
                    pskel(cc.PixelIdxList{k})=0;
                else
                    %Record as root if sufficiently long (i.e. not a stump, nodule, or debris)
                    if numel(cc.PixelIdxList{k})>prod(size(pskel))*.00005
                        n=n+1; %lat root number
                        [ry rx]=ind2sub(size(pskel),cc.PixelIdxList{k}); %indices for this root
                        tmp=zeros(size(pskel));tmp(cc.PixelIdxList{k})=1; %map of this root
                        natt(n)=sum(sum(tmp.*attm)); %# of attachment pts
                        if natt(n)==1
                            [ayt axt]=find(tmp.*attm==1);
                            atty(n)=ayt; attx(n)=axt;
                        else
                            atty(n)=nan; attx(n)=nan;
                        end
                        %                     apts(cc.PixelIdxList{k})
                        
                        root{n}=[rx ry];%Record coordinates of this root
                        latskel(cc.PixelIdxList{k})=1; %Record map of lateral roots
                    end
                end
            end
            
            latskel=latskel-pripath+attm;
            [laty latx]=find(latskel==1);
            %         figure; imshow(imread(c_imgpth))
            %         hold on; plot(latx, laty,'.g'); plot(attpx, attpy,'or')
            save([c_resfol '/latroots.mat'],'pskel','latskel','root','natt','attpy','attpx','attm','latx','laty','attx','atty')
        end
    end
end