function k = psf_traverse(ker,thres,psf_trace_dir,ts,te,isvis)

% load('rsmb_te30_tr_60_im_ 5_uniform/1_kermat_final.mat');

% get the latest kernel im and convert it to thin, BW curve

k = struct;
k.orig = ker;
clear ker;
% should do ISD!!!!!!! now smooth
% k.orig = imfilter(k.orig,fspecial('gaussian',[3 3], 0.5),'replicate');
tmp = sort(k.orig(:),'descend');
tmp = tmp(round(numel(find(tmp>0))*thres));
k.thin = bwmorph(im2bw(k.orig, tmp), 'thin', Inf); %imagesc([k.orig k.thin]);
% now trace the boundary
tmpb = bwboundaries(k.thin,'noholes');
% handle seperate dots noise
tmpind = 1;
tmpmax = -1;
if numel(tmpb) > 1
    for tmpi = 1:numel(tmpb)
        if size(tmpb{tmpi},1) > tmpmax
            tmpmax = size(tmpb{tmpi},1);
            tmpind = tmpi;
        end
    end
end
k.bound = tmpb(tmpind);
k.bound = k.bound{1};
% refine k.thin with only k.bound
tmp = zeros(size(k.orig));
for tmpi = 1:size(k.bound,1)
    tmp(sub2ind(size(k.orig),k.bound(tmpi,1),k.bound(tmpi,2))) = 1;
end
k.thin = tmp;
% find the endpoints only in k.bound
tmp = bwmorph(k.thin,'endpoints');
tmp = find(tmp);
k.endpt = zeros(2);
[i,j] = ind2sub(size(k.orig),tmp(1));
k.endpt(1,:) = [i,j];
[i,j] = ind2sub(size(k.orig),tmp(2));
k.endpt(2,:) = [i,j];
% crop into and end to end trajectory
[~,indx1]=ismember(k.endpt(1,:),k.bound,'rows');
[~,indx2]=ismember(k.endpt(2,:),k.bound,'rows');
k.endptind = zeros(2,1);
% k.endptind(1) = min(indx1,indx2);
% k.endptind(2) = max(indx1,indx2);
if strcmp(psf_trace_dir,'tb') || strcmp(psf_trace_dir,'bt')
    if k.bound(indx1,1) > k.bound(indx2,1) % indx1 is lower than indx2
        if strcmp(psf_trace_dir,'tb')
            k.endptind(1) = indx2; % head
            k.endptind(2) = indx1; % tail
        end
        if strcmp(psf_trace_dir,'bt')
            k.endptind(1) = indx1; % head
            k.endptind(2) = indx2; % tail
        end
    else % indx1 is at the top of indx2
        if strcmp(psf_trace_dir,'tb')
            k.endptind(1) = indx1; % head
            k.endptind(2) = indx2; % tail
        end
        if strcmp(psf_trace_dir,'bt')
            k.endptind(1) = indx2; % head
            k.endptind(2) = indx1; % tail
        end
    end
elseif strcmp(psf_trace_dir,'lr') || strcmp(psf_trace_dir,'rl')
    if k.bound(indx1,2) > k.bound(indx2,2) % indx1 is right to indx2
        if strcmp(psf_trace_dir,'lr')
            k.endptind(1) = indx2; % head
            k.endptind(2) = indx1; % tail
        end
        if strcmp(psf_trace_dir,'rl')
            k.endptind(1) = indx1; % head
            k.endptind(2) = indx2; % tail
        end
    else % indx1 is at the top of indx2
        if strcmp(psf_trace_dir,'lr')
            k.endptind(1) = indx1; % head
            k.endptind(2) = indx2; % tail
        end
        if strcmp(psf_trace_dir,'rl')
            k.endptind(1) = indx2; % head
            k.endptind(2) = indx1; % tail
        end
    end
end
        
tmp = imfilter(k.orig, ones(3), 'conv');
for i = 1:size(k.bound,1)
    k.bound(i,3) = tmp(k.bound(i,1),k.bound(i,2));
end
if k.endptind(1) < k.endptind(2)
    k.end2end = k.bound(k.endptind(1):k.endptind(2),:);
else
    k.end2end = k.bound(k.endptind(1):-1:k.endptind(2),:);
end
tmp = k.end2end(:,3);
k.end2end(:,3) = tmp./sum(tmp(:));
% the fitted kernal image
k.fitted = zeros(size(k.orig));
tmp = 0;
for i = 1:size(k.end2end,1)
    k.fitted(k.end2end(i,1),k.end2end(i,2)) = k.end2end(i,3);
    % a cool animation of trajectory
    if exist('isvis','var')
        if isvis
            imagesc(k.fitted); axis equal; drawnow;
            pause(0.1);
        end
    end
    if i > 1
        % tmp = tmp + (k.end2end(i,3)-k.end2end(i-1,3))/2;% choose the mid-point - 19-Oct-2014
        tmp = tmp + k.end2end(i,3);
    else
        tmp = tmp + k.end2end(i,3)/2;
    end
    k.end2end(i,4) = tmp*(te-ts)+ts; % this is the temporal information
end
% now fit k.end2end(:,1or2) and k.end2end(:,4) to get the translation
% function of time
% k.yoft = fit(k.end2end(:,4), k.end2end(:,1), 'poly4');
% k.xoft = fit(k.end2end(:,4), k.end2end(:,2), 'poly4');
% plot(k.yoft,k.end2end(:,4),k.end2end(:,1));
% plot(k.xoft,k.end2end(:,4),k.end2end(:,2));
