www.pudn.com > Digital_Image_Correlation_2010b.zip > multipeak_tracking.m, change:2010-11-20,size:9803b


%% 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;