% [K, KL] = AUTOKYMO(FX, FY)returns the position of the sample (membrane or
% bead) and the position of the laser. These positions are calculated along
% the trajectory of the laser.
%
% Syntax:  [k,kl] = autokymo(fx,fy)
%
% Inputs:
%    fx - function giving the correspondant between the tensions (Vx and
%    Vy) and the x position of the infrared laser in the image
%    fy - function giving the correspondant between the tensions (Vx and
%    Vy) and the y position of the infrared laser in the image
%
% Outputs:
%    k  - position of the sample
%    kl - position of the infrared laser
%
% Other m-files required: none
% Subfunctions: getVxVy, improfile_thick,imgaussian2d, ezfit
% MAT-files required: convertFct.mat
%

% Author: Chardes.C, IBDM
% email address: claire.chardes@univ-amu.fr
% 2012 

%------------- BEGIN CODE --------------

function [k, kl] = autokymo(fx,fy)

% tension text file
[fname, ffolder] = uigetfile('.txt', 'select text file', 'E:\images\OT\PERKIN\');
% read tensions
[Vx, Vy] = getVxVy(ffolder, fname);

% open images
[fname, ffolder] = uigetfile('.tif', 'select tiff image', 'E:\images\OT\PERKIN\');
info = imfinfo([ffolder, fname]);

% convert Vx, Vy into x,y positions
x = fx(Vx', Vy');
y = fy(Vx', Vy');

if (numel(x)>size(info, 1))
    x = x(end-size(info,1)+1:end);
    y = y(end-size(info,1)+1:end);
end
figure(28), plot(x,'xb'), xlabel('# image'), ylabel('position (pixel)'), hold on, plot(y, 'xr'), hold off

% determiner les temps ou il faut faire un kymo, on peut avoir plusieurs
% kymo dans un meme film
nIm0 = str2double(cell2mat(inputdlg('First image number', 'Kymograph', 1, {'1'})));
nIm1 = str2double(cell2mat(inputdlg('Last image number' , 'Kymograph', 1, {'1'})));

xx=x(nIm0:nIm1);
yy=y(nIm0:nIm1);

% linear regression to find coordinates of the kymograph
p = polyfit(xx(~isnan(xx)),yy(~isnan(xx)),1);
mini = (min(xx))-1.3*((max(xx) - min(xx)) );
maxi = (max(xx))+1.3*((max(xx) - min(xx)) );
if (maxi - mini > 25)
    mini = (min(xx))-(0.5*(max(xx) - min(xx)) );
    maxi = (max(xx))+(0.5*(max(xx) - min(xx)) );
end
step = abs(mini-maxi)/4;
xkymo = [mini:step:maxi];
ykymo = polyval(p,xkymo);


% control of the size of kymo
if (sqrt((xkymo(1)-xkymo(end))^2+(ykymo(1)-ykymo(end))^2)<20)
    dif = ceil((20-sqrt((xkymo(1)-xkymo(end))^2+(ykymo(1)-ykymo(end))^2))/2);
    xkymo(1) = xkymo(1)-dif;
    xkymo(end) = xkymo(end)+dif;
    ykymo = polyval(p,xkymo);
end

x0 = mean(xkymo-32);
x1 = mean(xkymo+31);
y0 = mean(ykymo-32);
y1 = mean(ykymo+31);

bead = fspecial('gaussian', 5);

% loop for kymograph construction
for nIm=nIm0:nIm1
        disp(nIm);
        % read image
        image = imread([ffolder, fname], 'tif', nIm);
        imageLaser = zeros(size(image));

        [al,bl,l] = improfile_thick(image, xkymo(1), xkymo(end), ykymo(1), ykymo(end), 3);
        xl = coord2dist(al,bl);
        
        if(~isnan(xx(nIm-nIm0+1)))
            imageLaser = imgaussian2d(size(image, 1), size(image, 2), 1, xx(nIm-nIm0+1), yy(nIm-nIm0+1));
            [al,bl,llaser] = improfile_thick(imageLaser, xkymo(1), xkymo(end), ykymo(1), ykymo(end), 3);
            xlaser = coord2dist(al,bl);
            
            imagekymolaser(:,nIm-nIm0+1) = llaser;
        else
            imagekymolaser(:,nIm-nIm0+1) = zeros(numel(l),1);
            llaser = imagekymolaser(:,nIm-nIm0+1);
        end
        
        imagekymo(:,nIm-nIm0+1) = l;
        
        % fit to obtain the kymograph (plot)
        
        % smooth before gaussian fitting
        a = 1;
        b = [1/3 1/3 1/3];
        s = l(:,1);
        l = filter(b,a,s-mean(s));
         
        if (numel(l)>30)
            bb = l(5:end-4)-min(l);
            aa = xl(5:numel(l)-4);
        else
            bb = l;
            aa = xl;
        end
        
        f = ezfit(aa,bb, 'gauss');
        while(f.r<0.85 || numel(aa)<10)
            aa = aa(2:end-1);
            bb = f.y(2:end-1);
            try
                f = ezfit(aa, bb-min(bb), 'gauss');
            catch exception
                disp('exeception !!!');
                disp(exception);
                f.m(3)= NaN;
                break;
            end
        end
        if (f.m(3)<numel(l))
            k(nIm-nIm0+1) = f.m(3);
        else
            try
                f.m(3) = find(l==max(l));
            catch exception
                f.m(3)=NaN;
            end
                       
            k(nIm-nIm0+1) = f.m(3);
        end
        f = ezfit(xlaser, llaser-min(llaser), 'gauss');
        kl(nIm-nIm0+1) = f.m(end);
        disp(kl(nIm-nIm0+1));
        
end

% display and save

figure(29); 
imshow(imagekymo, []); title('kymograph');
pause(1);
colormap('HSV');
imRGB(:,:,1) = imagekymo/max(max(imagekymo));
imRGB(:,:,2) = imagekymolaser/max(max(imagekymolaser));
imRGB(:,:,3) = zeros(size(imagekymo));
imshow(imRGB,[]);
titleStr = ['Kymograph ', num2str(nIm0), '-', num2str(nIm1)];
title(titleStr);
filename = [ffolder, fname(1:end-4),'-kymo-',num2str(nIm0), '-', num2str(nIm1), '.tif'];
imwrite(imagekymo/max(max(imagekymo)), filename, 'tif');
filename = [ffolder, fname(1:end-4),'-kymoLaser-',num2str(nIm0), '-', num2str(nIm1), '.tif'];
imwrite(imagekymolaser, filename, 'tif');

figure(30); 
plot(k,'x-');
hold on;
plot(kl, '-*r');
legend('sample position','laser position');
titleStr = ['Kymograph ', num2str(nIm0), '-', num2str(nIm1)];
title(titleStr);
xlabel('# images');
ylabel('position (pixels)');
hold off;
filename = [ffolder, fname(1:end-4),'-kymo-',num2str(nIm0), '-', num2str(nIm1), '.jpg'];
saveas(gcf, filename, 'jpg');
saveas(gcf, [filename(1:end-3), 'fig'], 'fig');

end
%------------- END OF CODE --------------