[validx,validy]=multipeak_tracking;

%% written by Chris

% updated last 9/8/2010
% this function allows to track multiple peaks along one axis. The
% resolution can be much higher than the cross correlation methods of
% digital image correlation.

% function goes here
function [validx,validy]=multipeak_tracking;

%% Load image and ask for orientation and pos/neg

if exist('filenamelist')==0
    load('filenamelist')            % file with the list of filenames to be processed
end

filenumber=size(filenamelist);
filenumber=filenumber(1);

h=figure;
imshow(filenamelist(1,:))           % show the first image
title('First Image')        % put a title

ImageA=(uint8(mean(double(imread(filenamelist(1,:))),3)));

sel_pos_neg = menu(sprintf('Do you want use the negative image for analysis?'),...
    'As is', 'Negativ', 'Cancel')
if sel_pos_neg==3
    close all;
    return
end
if sel_pos_neg==2
    ImageA=255-ImageA;
    imshow(ImageA)           % show the first image
    negativ=1;
end
if sel_pos_neg==1
    negativ=0;
end

sel_horz_vert = menu(sprintf('Do you want to measure horizontal or vertical?'),...
    'Horizontal', 'Vertical', 'Cancel')
if sel_horz_vert==3
    close all;
    return
end
if sel_horz_vert==2
    ImageA=ImageA';
    imshow(ImageA)           % show the first image
    orientation=90;
end
if sel_horz_vert==1
    orientation=0;
end

[imr,imc]=size(ImageA);


%% select path for displacement extraction along markers
%Select point in middle of marker
xlabel('Location On Image [Pixels]')
ylabel('Location On Image [Pixels]')
title(['Click on the center of the sample. Width: ', num2str(imc),...
    '; Height: ',num2str(imr)])
marker_pt=round(ginput(1));
x_mark = marker_pt(1);
y_mark = marker_pt(2);
hold on
line([1 imc],[y_mark y_mark],'Color','r')


%Prompt for integration width
prompt = {'Enter integration width [Pixels]:'};
dlg_title = 'Input integration width for the analysis';
num_lines= 1;
def     = {'20'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
int_width = str2num(cell2mat(answer(1)));

line([1 imc],[y_mark-int_width/2 y_mark-int_width/2],'Color','g')
line([1 imc],[y_mark+int_width/2 y_mark+int_width/2],'Color','g')


%% Calculate line profile data and select peaks
xdata = [1:1:imc];
ydata= sum(ImageA((y_mark-int_width/2):(y_mark+int_width/2),:),1)/int_width;

g=figure;
plot(xdata,ydata)
xlabel('Location On Image [Pixels]')
ylabel('Pixel Intensity Value')

sel_peakselect = menu(sprintf('Do you want to select single or all peaks?')...
    ,'Single Select', 'All Select', 'Cancel')
if sel_peakselect==3
    close all;
    return
end

% If single peaks should be chosen
if sel_peakselect==1
    sel_peakselect_2=1
    i=1
    hold on
    while sel_peakselect_2==1
        title('Click on the left of the chosen peak.')
        marker_pt=round(ginput(1));
        x_peak_left = marker_pt(1);
        y_peak_left = marker_pt(2);
        line([x_peak_left x_peak_left], [1 max(ydata)],'Color','r')
        
        title('Click on the right of the chosen peak.')
        marker_pt=round(ginput(1));
        x_peak_right = marker_pt(1);
        y_peak_right = marker_pt(2);
        line([x_peak_right x_peak_right], [1 max(ydata)],'Color','r')
        
        plot(xdata(x_peak_left:x_peak_right),ydata(x_peak_left:x_peak_right),'k')
        peak_pos(i,:)=[x_peak_left x_peak_right];
        
        i=i+1;
        sel_peakselect_2 = menu(sprintf('Do you want to select another peak?')...
            ,'One more', 'No, lets process')
    end
end

% If all peaks should be found automatically
if sel_peakselect==2
    
end

%% Start first fitting
fitcountertemp=size(peak_pos); % number off peaks to cycle through
fitcounter=fitcountertemp(1,1); % number off peaks to cycle through

for c=1:fitcounter %start the loop to process all points
    fitdatax=xdata(peak_pos(c,1):peak_pos(c,2));
    fitdatay=ydata(peak_pos(c,1):peak_pos(c,2));
    
    % fitting in x-direction
    % guess some parameters for the fitting routine --> bad guesses lead to
    % an error message which stops the fitting
    
    back_guess=(ydata(peak_pos(c,2))+ydata(peak_pos(c,2)))/2; % guess for the background level - average of the first and last greyvalue
    sig1_guess=(peak_pos(c,2)-peak_pos(c,1))/2; % guess for the peak width - take half of the cropping width
    amp_guess1=max(fitdatay); % guess for the amplitude - take the greyvalue at the peak position
    mu1_guess=(peak_pos(c,2)+peak_pos(c,1))/2; % guess for the position of the peak - take the position from bwlabel
    
    % start fitting routine
    [x,resnormx,residual,exitflagx,output]  = lsqcurvefit(@gauss_onepk, [amp_guess1 mu1_guess sig1_guess back_guess], fitdatax, fitdatay); %give the initial guesses and data to the gauss fitting routine
    
    % show the fitting results
    xtest = fitdatax; % x values for the plot of the fitting result
    ytest = (x(1)*exp((-(xtest-x(2)).^2)./(2.*x(3).^2))) + x(4); % y values of the fitting result
    yguess=(amp_guess1*exp((-(xtest-mu1_guess).^2)./(2.*sig1_guess.^2))) + back_guess; %y values for the guess plot
    plot(xtest,ytest,'r') % plot the fitted function
    plot(xtest,yguess,'b') % plot the guessed function
    drawnow
    
    peaks(c,1,:)=x;
end

sel_go = menu(sprintf('Do you want proceed?')...
    ,'Do it!', 'No, stop!')

if sel_peakselect==2
    close all;
    return
end

cutout=peak_pos(:,2)-peak_pos(:,1)

for m=1:(filenumber-1) % loop through all images
    ImageA=((mean(double(imread(filenamelist(m,:))),3)));
    
    if sel_horz_vert==2
        ImageA=ImageA';
        imshow(ImageA)           % show the first image
        orientation=90;
    end
    if sel_pos_neg==2
        ImageA=255-ImageA;
        imshow(ImageA)           % show the first image
        negativ=1;
    end
    
    image=m
    
    xdata = [1:1:imc];
    ydata= sum(ImageA((y_mark-int_width/2):(y_mark+int_width/2),:),1)/int_width;
    
    plot(xdata,ydata)
    title(['Image Number: ', num2str(m)])
    hold on
    
    
    for c=1:fitcounter %start the loop to process all points
        
        if (peaks(c,m,2)+cutout(c)/2)>imc
        fitdatax=xdata((imc-cutout(c)):imc);
        fitdatay=ydata((imc-cutout(c)):imc);
        else
        fitdatax=xdata((peaks(c,m,2)-cutout(c)/2):(peaks(c,m,2)+cutout(c)/2));
        fitdatay=ydata((peaks(c,m,2)-cutout(c)/2):(peaks(c,m,2)+cutout(c)/2));
        end
        
        % fitting in x-direction
        % guess some parameters for the fitting routine --> bad guesses lead to
        % an error message which stops the fitting
        
        back_guess=peaks(c,m,4); % guess for the background level - average of the first and last greyvalue
        sig1_guess=peaks(c,m,3); % guess for the peak width - take half of the cropping width
        amp_guess1=peaks(c,m,1); % guess for the amplitude - take the greyvalue at the peak position
        mu1_guess=peaks(c,m,2); % guess for the position of the peak - take the position from bwlabel
        
        % start fitting routine
        [x,resnormx,residual,exitflagx,output]  = lsqcurvefit(@gauss_onepk, [amp_guess1 mu1_guess sig1_guess back_guess], fitdatax, fitdatay); %give the initial guesses and data to the gauss fitting routine
        % show the fitting results
        xtest = fitdatax; % x values for the plot of the fitting result
        ytest = (x(1)*exp((-(xtest-x(2)).^2)./(2.*x(3).^2))) + x(4); % y values of the fitting result
        yguess=(amp_guess1*exp((-(xtest-mu1_guess).^2)./(2.*sig1_guess.^2))) + back_guess; %y values for the guess plot
        plot(xtest,ytest,'r') % plot the fitted function
        plot(xtest,yguess,'g') % plot the guessed function
        drawnow
        peaks(c,m+1,:)=x;
    end
    hold off
    Amplitude=peaks(:,m,1)';
    Peakposition=peaks(:,m,2)';
    Peakwidth=peaks(:,m,3)';
    Background=peaks(:,m,4)';
    dlmwrite('backupPeakposition.txt', Peakposition' , 'delimiter', '\t', '-append');       % Here we save the result from each image; if you are desperately want to run this function with e.g. matlab 6.5 then you should comment this line out. If you do that the data will be saved at the end of the correlation step - good luck ;-)
    dlmwrite('backupAmplitude.txt', Amplitude' , 'delimiter', '\t', '-append');       % Here we save the result from each image; if you are desperately want to run this function with e.g. matlab 6.5 then you should comment this line out. If you do that the data will be saved at the end of the correlation step - good luck ;-)
    dlmwrite('backupPeakwidth.txt', Peakwidth' , 'delimiter', '\t', '-append');       % Here we save the result from each image; if you are desperately want to run this function with e.g. matlab 6.5 then you should comment this line out. If you do that the data will be saved at the end of the correlation step - good luck ;-)
    dlmwrite('backupBackground.txt', Background' , 'delimiter', '\t', '-append');       % Here we save the result from each image; if you are desperately want to run this function with e.g. matlab 6.5 then you should comment this line out. If you do that the data will be saved at the end of the correlation step - good luck ;-)
end

%% Save data
validx=peaks(:,:,2);
validy=peaks(:,:,3);

Amplitude=peaks(:,:,1)';
Peakposition=peaks(:,:,2)';
Peakwidth=peaks(:,:,3)';
Background=peaks(:,:,4)';

save peaks;
save Amplitude.dat Amplitude -ascii -tabs;
save Peakposition.dat Peakposition -ascii -tabs;
save Peakwidth.dat Peakwidth -ascii -tabs;
save Background.dat Background -ascii -tabs;
save validx_mp.dat validx -ascii -tabs;
save validy_mp.dat validy -ascii -tabs;