% function PrimaryRoot(mpd)
% Get primary root end coordinates through user input

%% 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 '/primary.mat'])
            load([c_resfol '/skel.mat'],'skel','allx','ally');
            
            %% Input top and tip points
            %%Display  image
            im=imread(c_imgpth);
            fo=figure; %close all;
            warning('off','images:imshow:magnificationMustBeFitForDockedFigure')
            warning('off', 'Images:initSize:adjustingMag');
            imshow(im); hold on; set(gcf,'NumberTitle','off','Name',c_data(1:end))
            disp(['Processing ' c_data(1:end-3) '...']);
            
            %%Input points
            title('Click the TOP of the primary root','FontWeight','Bold','FontSize',12)
            ptop=ginput(1);
            d=distance([allx ally],[ptop(1) ptop(2)]);%Finding nearest coordinate on the skeleton
            [val ind]=min(d);
            topx=allx(ind); topy=ally(ind);
            hold on; plot(ptop(1), ptop(2),'.r','MarkerSize',12);
            
            title('Click the TIP of the primary root')
            ptip=ginput(1);
            d=distance([allx ally],[ptip(1) ptip(2)]); %Finding nearest coordinate on the skeleton
            [val ind]=min(d);
            tipx=allx(ind); tipy=ally(ind);
            hold on; plot(tipx,tipy,'.r','MarkerSize',12);
            
            
            %% Trace the root coordinates
            ok='x'; pp=[];
            while ~isempty(ok)
                
                %%Find the shortest path between the two points to get the primary root
                D1=bwdistgeodesic(skel, topx, topy);
                D2=bwdistgeodesic(skel, tipx, tipy);
                D=D1+D2;
                pripath=zeros(size(skel));
                pripath(D==D(topy, topx))=1;
                
                %%Trace
                prirc=bwtraceboundary(pripath,[topy topx],'s',8);
                [C ia]=unique(prirc(:,2)); %sometimes doubles back on itself
                prirc=prirc(ia,:);
                prirc=sortrows(prirc,2); %put points in order
                prix=prirc(:,2); priy=prirc(:,1);
                set(pp,'Visible','off');
                pp=plot(prirc(:,2),prirc(:,1),'-m');
                
                %% Edit
                ylabel('Note: The precision of the trace is only important in the nodulation zone for the Plane1&2 images')
                title('Press enter if okay.')
                xlabel({'Left-click: force trace a pt'; 'Right-click: remove a pt';'Backspace to undo all';'z to zoom';'p to pan';'0 to manually trace'},...
                    'Units', 'normalized', 'Position', [1 1], 'HorizontalAlignment', 'left')
                
                [x y ok]=ginput(1);
                if ~isempty(ok)
                    
                    d=distance([allx ally],[x y]); %Finding nearest coordinate on the skeleton
                    [val ind]=min(d);
                    ex=allx(ind); ey=ally(ind);
                    
                    if ok<10 %Edit existing trace
                        if ok==1 %segment all other paths
                            skel(:,ex)=0;
                            skel(ey,ex)=1;
                        elseif ok==3 %remove from the path
                            skel(ey,ex)=0;
                        elseif ok==8 %undo all and start over
                            load([c_resfol '/skel.mat'],'skel');
                        end
                        
                    else %Zoom or pan
                        if okay>48
                            set(fo,'Tag','Original')
                            fop=get(fo,'Position');
                            h = uicontrol('String', 'Resume', 'Position',[fop(3)-120 fop(4)-80 80 30], 'Callback', 'set(fo,''Tag'',''Mod'')','BackgroundColor',[.7 1 .7] );                            
                            if ok==122
                                zoom on
                                drawnow
                                waitfor(fo, 'Tag');
                            elseif ok==112
                                pan on
                                drawnow
                                waitfor(fo, 'Tag');
                            end
                            set(h,'Visible','off')
                        elseif okay==48
                            [prix,priy]=draw_line_unclosed;
                        end
                    end
                end
            end
            
            %     %% Smooth out the root coordinates
            %     prixs=min(prirc(:,2)):20:max(prirc(:,2));
            %     priys=spline(prirc(:,2)',prirc(:,1)',prixs);
            
            %%Save primary root end coordinates
            save([c_resfol '/primary.mat'],'tipx','tipy','topx','topy','prirc','pripath','prix','priy');
        end
    end
end
