% function SegLaterals
%Segment roots if needed through user input

%v2 if the root has not been segmented, get user input
%v4 use user input nattachment points

%clear all
close all
%Load dataset lists
if ~exist('mpd')
    mpd=uigetdir(pwd,['Select the main project folder']);
end; 
load([mpd '\dataset.mat']);
figure
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 '/individroots.mat']);
            load([c_resfol '/latroots.mat'],'latskel','pskel','root','natt','attpy','attpx','attm','attx','atty')
            load([c_resfol '/primary.mat'],'prix','priy');
            im=imread(c_imgpth);
            warning('off','images:imshow:magnificationMustBeFitForDockedFigure')
            cla
            clf
            imshow(im); hold on; set(gcf,'NumberTitle','off','Name',c_data)
            [alaty alatx]=find(latskel==1);
            
            plot(alatx,alaty,'.w','MarkerSize',.5)
            plot(prix,priy,'.y','MarkerSize',.5)
            N=0; %number of lateral roots (root ID)
            xyg=[]; %coordinates of all lateral roots
            attg=[]; %List of attchmnt pts
            
            
            %% Display all auto-IDed individual roots for verification
            disp('PART 1: Displaying each auto-identified root in green for you to verify (or fix)')
            for j=1:size(natt,2)
                if natt(j)==1 %If there was 1 and only 1 nattachment
                    
                    %%Display it
                    ptop=plot(attx(j),atty(j),'.r');
                    proot=plot(root{j}(:,1),root{j}(:,2),'.g');
                    
                    %Check it
                    ok=input('Hit enter if this root is okay as is, 0 if not okay at all, 1 to add more points to the end: ');
                    if isempty(ok)
                        disp('Lateral root okayed and saved')
                        set(proot,'color','y')
                        N=N+1;
                        xyg=[xyg; root{j} N*ones(length(root{j}),1) ones(length(root{j}),1)];%new in v11 : column 4: 1 if auto, 0 if manual
                        attg=[attg;attx(j) atty(j)];
                        
                    elseif ok==1 %if good but need to add
                        disp('Click points along the part of the root that was missed, starting from the point that the automatically extracted part left off, down to the root tip')
                        [xs ys latpp]=ginzoompan_freeform; %add pts
                        disp('The additional points were added to the root coordinates')
                        %                         [yt ti]=max(root{j}(:,2));
                        %                         ys=yt:max(laty);                % Smooth out the root coordinates bw pts
                        %                         xs = interp1([yt; laty]',[root{j}(ti,1);latx]',ys);
                        set(proot,'color','g')
                        [basey lowind]=min(root{j}(:,2));
                        %proot=plot([root{j}(lowr,1);xs],[lowy;ys],'y');
                        lat=cat(1,root{j},[xs ys]);
                        aom=[ones(size(root{j},1),1); zeros(size(xs,1),1)];%for v11
                        proot=plot(lat(:,1),lat(:,2),'.y');
                        
                        N=N+1;
                        xyg=[xyg; lat(:,1) lat(:,2) N*ones(size(lat,1),1) aom];%for v11
                        attg=[attg;attx(j) atty(j)];
                    else
                        set(proot,'visible','off')
                        disp('That "root" was omitted.')
                    end
                end
            end
            
            for k=1:size(xyg,1)
                latskel(round(xyg(k,2)),round(xyg(k,1)))=0;
            end
            
            %% Get remaining roots from user input
            disp('PART 2: Now getting any remaining lateral roots from user input')
            more=input('Press ENTER to add more lateral roots, or 0 if all roots are now traced out (in yellow)');
            while isempty(more)
                lat=[];
                aom=[];
                ok=input('Are at least most of the coordinates of the root in question extracted (i.e. plotted in white)? Press ENTER if yes, 0 if no');
                if isempty(ok)
                    %Get attachment point
                    title('Click the root"s attachment point on the primary root. ')
                    [xa ya ppa]=ginzoompan(attpx,attpy,'.r');
                    
                    ok=input('Press enter if this is okay. 0 if you need to choose a different attachment point.');
                    if ~isempty(ok)
                        set(ppa,'visible','off')
                        [xa ya ppa]=ginzoompan(prix,priy,'.r');
                    end
                    
                    %Get tip
                    title('Click the corresponding root tip (or the end point of the extracted coordinates (white points))')
                    [xt yt ppt]=ginzoompan(alatx,alaty,'.m');
                    
                    %Trace the coordinates between the points
                    D1=bwdistgeodesic(logical(latskel), xa, ya);
                    D2=bwdistgeodesic(logical(latskel), xt, yt);
                    D=D1+D2;
                    skpath=zeros(size(latskel));
                    skpath(D==D(ya, xa))=1;
                    lat=bwtraceboundary(skpath,[ya xa],'s',8);
                    
                    %% Check if coordinates are good
                    if ~isempty(lat)
                        lat=[lat(:,2) lat(:,1)];
                        proot=plot(lat(:,1),lat(:,2),'.g');
                        
                        ok=input('Press enter if okay as is, 1 to add points, or 0 to delete the trace and do it all manually: ');
                        set(ppa,'visible','off'); set(ppt,'visible','off')
                        if isempty(ok) || ok==1
                            set(proot,'color','y')
                            aom=ones(size(lat,1),1);%'automatically extracted' tag - v11
                        elseif ok==0
                            set(proot,'visible','off')
                            lat=[];
                        end
                    else
                        disp('You will have to do this root manually.')
                        set(ppt,'visible','off'); set(ppa,'visible','off')
                        ok=0;yt=ya;xt=xa;
                        lat=[];%AG Oct 8 2013
                        aom=[];%v11
                    end
                end
                
                if ~isempty(ok)
                    title('Click points along the root (in order, ending at the tip)')
                    [xs ys latpp]=ginzoompan_freeform;
                    %                 % Smooth out the root coordinates bw pts
                    %                 ys=yt:max(laty);
                    %                 xs = interp1([yt; laty]',[xt;latx]',ys);
                    set(latpp,'visible','off')
                    proot=plot(xs,ys,'.y');
                    if ~isempty(lat)
                        lat=cat(1,lat,[xs(1:end) ys(1:end)]);
                        aom=cat(1,aom,zeros(size(xs,1),1));%manually extracted tag
                    else
                        lat=cat(1,lat,[xs(1:end) ys(1:end)]); %was starting from 2, start from 1 again - for 2 points roots
                        aom=cat(1,aom,zeros(size(xs,1),1));%manually extracted tag
                        xa=xs(1);%AG
                        ya=ys(1);%AG
                    end
                    %lat=[ys(2:end) xs(2:end)];%AG Oct 8 2013 removed cat function

                end
                
                %% Record root coordinates
                N=N+1;
                xyg=[xyg; lat(:,1) lat(:,2) N*ones(size(lat,1),1) aom];%Nowxyg has 4 columns: last one: auto (1) or manual (0)
                attg=[attg;xa ya];
                more=input('Hit enter to add more roots, or 0 when done');
            end
            
            %splitting as 2 matrices for compatibility with later programs
            %which use xyg with 3 columns
            aom=xyg(:,4);
            xyg=xyg(:,1:3);
            
            save([c_resfol '/individroots.mat'],'attg','xyg','aom');
        end
    end
end



