GetPath.m

function [iPath,jPath,bKymo]=GetPath(im,p1,p2)
iPath=[];jPath=[];bKymo=[];
RADIUS_DILATE=3; 

im=im2double(im);
im=mat2gray(im);

alpha=1;
level=alpha*graythresh(im);
bw=im2bw(im,level);
se = strel('disk',RADIUS_DILATE);
bw=imdilate(bw,se);
[L num]=bwlabel(bw);
while alpha>0.5 & L(p1(1),p1(2))~=L(p2(1),p2(2))
    alpha=alpha-.1;
    level=alpha*graythresh(im);
    bw=im2bw(im,level);
    bw=imdilate(bw,se);
    [L num]=bwlabel(bw);
end

if  L(p1(1),p1(2))~=L(p2(1),p2(2))
    fprintf( 1,' No PATH found!!');
    return
end

bw(find(L~=L(p1(1),p1(2))))=0;
im(~bw)=0;
% hold off;imagesc(im)
% ;imshow(bw)
colormap(gray);hold all
% [r c]=find(bwSkel);
% plot(c,r,'.g','markersize',4)

% d1 = (p1(1)-r).^2 + (p1(2)-c).^2;
% [m1 idx1]=min(d1);
% p1 = [r(idx1) c(idx1)];
% 
% d2 = (p2(1)-r).^2 + (p2(2)-c).^2;
% [m2 idx2]=min(d2);
% p2 = [r(idx2) c(idx2)];
% 
% im= imresize(im, 2);
plot(p1(2),p1(1),'ro')
plot(p2(2),p2(1),'bx')
drawnow
Progressbar(0);
gIm=sparse([],[],[],prod(size(im)),prod(size(im)),8*sum(size(im)));
for i=2:size(im,1)-1
    Progressbar(i/size(im,1));
    for j=2:size(im,2)-2
        if 0==im(i,j),continue,end
        for ni=-1:1
            for nj=-1:1
                if ni==0 && nj==0,continue,end
                i2=i+ni;
                j2=j+nj;
                if 0==im(i2,j2),continue,end

                v1=(p1-[i j]);
                v2=(p2-[i j]);
                v3=([i2 j2]-[i j]);
                v1=v1/norm(v1);
                v2=v2/norm(v2);
                v3=v3/norm(v3);
                cost=2-(abs(dot(v1,v3))+abs(dot(v2,v3)));
                cost2=2-(im(i2,j2))-(im(i,j));
                cost=cost+4*cost2;
% cost=1;

                ind1=sub2ind(size(im),i,j);
                ind2=sub2ind(size(im),i2,j2);
                
                gIm(ind1,ind2)=cost;
                gIm(ind2,ind1)=cost;
            end
        end
    end
end
Progressbar(1);

ind1=sub2ind(size(im),p1(1),p1(2));
ind2=sub2ind(size(im),p2(1),p2(2));
options.target=ind2;
 [d pred]=shortest_paths(gIm,ind1);
pp=pred(ind2)
iPath=[];
jPath=[];
while(pp & pp~=ind1)
    [i j]=ind2sub(size(im),pp);
    iPath=[iPath;i];
    jPath=[jPath;j];
    pp=pred(pp);
end

plot(jPath,iPath,'-r')

if isempty(iPath)
    imKymo=[];
    return
end
% smooth
bKymo=zeros(length(iPath),2);
SWIN=6;INWIN=3;
% smooth

% 1st 10
x=jPath(1:3*SWIN);y=iPath(1:3*SWIN);
[b stats]=robustfit(x,y);
[b2 stats2]=robustfit(y,x);
if stats2.s<stats.s
    % use y as dependent var
    for j=1:SWIN
        jPath(j)=b2(1)+iPath(j)*b2(2);
        bKymo(j,1:2)=[NaN NaN];
        bKymo(j,3:4)=b2;
    end
else
    for j=1:SWIN
        iPath(j)=b(1)+jPath(j)*b(2);
        bKymo(j,1:2)=b;
    end
end

% last 10
% 1st 10
x=jPath(end-2*SWIN-INWIN:end);y=iPath(end-2*SWIN-INWIN:end);
[b stats]=robustfit(x,y);
[b2 stats2]=robustfit(y,x);
if stats2.s<stats.s
    % use y as dependent var
    for j=0:SWIN+INWIN
        jPath(end-j)=b2(1)+iPath(end-j)*b2(2);
        bKymo(end-j,3:4)=b2;
        bKymo(end-j,1:2)=[NaN NaN];
    end
else
    for j=0:SWIN+INWIN
        iPath(end-j)=b(1)+jPath(end-j)*b(2);
    end
    bKymo(end-j,1:2)=b;
end

for i=SWIN+1:INWIN:length(jPath)-SWIN -INWIN
    x=jPath(i-SWIN:i+INWIN+SWIN);y=iPath(i-SWIN:i+INWIN+SWIN);
    [b stats]=robustfit(x,y);
    [b2 stats2]=robustfit(y,x);
    if stats2.s<stats.s
         % use y as dependent var
        for j=0:INWIN-1
            jPath(i+j)=b2(1)+iPath(i+j)*b2(2);
            bKymo(i+j,3:4)=b2;
            bKymo(i+j,1:2)=[NaN NaN];
        end
    else
        for j=0:INWIN-1
            iPath(i+j)=b(1)+jPath(i+j)*b(2);
            bKymo(i+j,1:2)=b;
        end
    end
end

% [jPath,iPath,bKymo] = FillPathGaps(jPath,iPath,bKymo)
plot(jPath,iPath,'-b')
% print(ftrace,'-dtiffn','-r300')