disp(' ');
disp(' loglog_digitize.m    ver 1.5   January 10, 2011');
disp(' ');
disp(' By Tom Irvine   Email: tomirvine@aol.com');
%
%  This program is a modified version of a similar program written by 
%  Kent Hardy and Ryan Edwards.
%  This function allows the user to pick points off an image file.  The
%  intended use is to be able to scan a hard copy of an SRS or PSD plot. 
%  
%  Variables:
%    "image_file" is an image file which you have already loaded (perhaps
%  with the "imread" function).  The image must be upright when viewed with
%  the "imagesc" function, so you may need to rotate the data somewhat.  I
%  recommend that you try the function "rot90" to accomplish this.
%  
%
%  PROGRAM USE:
%    When you execute the function, the program will prompt you to pick
%  the lower left corner of the data plot, and then enter the frequency
%  and SRS value at that point.  Then, you will pick the lower right
%  corner of the data plot and input the frequency at that point.  Third,
%  you will pick the upper left corner of the data plot and input the
%  SRS magnitude at that point.  Finally, you will pick the points on the
%  visible SRS curve.  Remember to hit RETURN when you are done picking.
%  The program will then plot what you have chosen on top of the image so
%  you can see if you like how it fits together.
%    Once you have given the OK, the code converts the data points to
%  frequency/SRS format.  If you pick a point near the lower or upper
%  frequency bounds, the code will assume you meant that the point should
%  be at the maximum/minimum frequency of the plot.  This should help
%  later when you try to put everything to a common frequency resolution.
%
clear nn;
clear xx;
clear yy;
clear xxyy;
clear xxyy3;
clear srs_out;
clear psd_out;
clear figure(1);
%
disp(' ')
disp(' Select file input method ');
disp('   1=external image file ');
disp('   2=image preloaded into Matlab ');
file_choice = input('');
%
%% if(file_choice==1)
%%    disp(' Enter the input image filename: ');
%%    filename = input(' ','s');
%%    image_file = imread(filename);
%%else
%%    image_file = input(' Enter the matrix name:  ');
%%end
%%%%%%%%%%%%%%%%%%
if(file_choice==1)
    [filename, pathname] = uigetfile('*.*');
    filename = fullfile(pathname, filename);
%      
    image_file = imread(filename);
%
else
    image_file = input(' Enter the image name:  ');
end
%%%%%%%%%%%%%%%%%%
%
disp(' ');
disp(' Select graph type ');
disp(' 1=SRS  2=PSD  3=other');
itype = input(' ');
%
disp(' ')
disp('Output ASCII text file?')
choice=input(' 1=yes  2=no  ' );
disp(' ')
%
if choice == 1 
   disp(' Enter output filename ');
   TH_filename = input(' ','s');
%
   fid = fopen(TH_filename,'w');
end
%
yesno='No';
while yesno(1)~='Y';
  clf;
  imagesc(image_file);
  axis image off; % Turn off the image axes
  set(gca,'Position',[0 0 1 1]);
  set(gcf,'colormap',[[1 0 0]; [0 0 0]]);
 
  zoommsg = {'Zoom in on the portion of the plot that you want to digitize.',...
        'Take care to keep the axes limits visible.',...
        'Press OK when you want to proceed'};
  zoom on;
  q0 = msgbox(zoommsg);
  waitfor(q0);
  
  q1=questdlg('Pick lower left corner of plot. ','','Proceed','Cancel','Proceed');
  if q1(1)=='P'; x1y1=ginput(1); else; error('Shutting Down'); end
%
  if(itype==1)
    txt2=cellstr(strvcat('Natural Frequency (Hz)','SRS (G)'));
  end  
  if(itype==2)
    txt2=cellstr(strvcat('Frequency (Hz)','PSD (G^2/Hz)'));   
  end
  if(itype~=1 & itype~=2)
    txt2=cellstr(strvcat('X-axis value','Y-axis value'));     
  end
%
  q1=inputdlg(txt2);
  q1a=char(q1);
  f1=str2double(q1a(1,:));
  s1=str2double(q1a(2,:));
  %
  q1=questdlg('Pick lower right corner of plot. ','','Proceed','Cancel','Proceed');
  if q1(1)=='P'; x2y2=ginput(1); else; error('Shutting Down'); end
% 
  if(itype==1)
     txt2=cellstr('Natural Frequency (Hz)');
  else
     txt2=cellstr('Frequency (Hz)');     
  end  
%
  q1=inputdlg(txt2);
  q1a=char(q1);
  f2=str2double(q1a(1,:));
  %
  q1=questdlg('Pick upper left corner of plot. ','','Proceed','Cancel','Proceed');
  if q1(1)=='P'; x3y3=ginput(1); else; error('Shutting Down'); end
%
  if(itype==1)
     txt2=cellstr('SRS (G)');
  else
     txt2=cellstr('PSD (G^/Hz)');   
  end
%
  q1=inputdlg(txt2);
  q1a=char(q1);
  s2=str2double(q1a(1,:));
  %
  %  Plot the lines which were chosen and give the user a chance to change.
  %
  x1=[x1y1(1) x2y2(1)]; y1=[x1y1(2) x2y2(2)];
  q1=line(x1,y1); set(q1,'color',[1 1 0]);
  x1=[x1y1(1) x3y3(1)]; y1=[x1y1(2) x3y3(2)];
  q1=line(x1,y1); set(q1,'color',[1 1 0]);
  yesno=questdlg('Are the axes acceptable?');
end
%
yesno='No';
while yesno(1)~='Y';
%
if(itype==1)
  txt=[
    'Pick the points of the SRS curve.'
    'Choose as many points as desired.'
    'Press RETURN when done choosing. '];
else
   txt=[
    'Pick the points of the PSD curve.'
    'Choose as many points as desired.'
    'Press RETURN when done choosing. '];   
end
%
  q1=questdlg(txt,'','Proceed','Cancel','Proceed');
  if q1(1)=='P'; [xx,yy]=ginput; else; error('Shutting Down'); end
  %
  %  Overlay chosen points on image.
  %
  numpts=length(xx);
  for k1=1:numpts-1
    x1=[xx(k1) xx(k1+1)];
    y1=[yy(k1) yy(k1+1)];
    q1=line(x1,y1);
    set(q1,'color',[0 1 1]);
  end
  yesno=questdlg('Are the points acceptable?');
  if strcmpi(yesno,'No')
      q1 = findall(gca,'Color',[0 1 1]);
      delete(q1);
  end
end
xxyy=[xx yy]';
%
%  Translate the data into actual SRS numbers assuming a
%  loglog plot.  Start by rotating the data to a directly
%  x-y coordinate system (this to account for the fact
%  that the scanner often rotates the page slightly).
%
thet1=atan((x1y1(2)-x2y2(2))/(x2y2(1)-x1y1(1)));
thet2=atan((x3y3(1)-x1y1(1))/(x3y3(2)-x1y1(2)));
theta=(thet1+thet2)/2;
%  Create coordinate transformation matrix.
trmat=[cos(theta) -sin(theta)
  sin(theta) cos(theta)];
%  Rotate the axes and the curve.
xxyy2=trmat*xxyy;
xxyy2=xxyy2';
x1y1a=trmat*x1y1';
x2y2a=trmat*x2y2';
x3y3a=trmat*x3y3';
%
%  Convert to log data.
%
freq_decs=log10(f2/f1);
srs_decs=log10(s2/s1);
freq_length=(x2y2a(1)-x1y1a(1))/freq_decs;
srs_length=(x1y1a(2)-x3y3a(2))/srs_decs;
for k1=1:numpts
  xxyy3(k1,1)=f1*10^((xxyy2(k1,1)-x1y1a(1))/freq_length);
  %  If frequency point is near the upper or lower limit, then
  %  assume that the user intended to pick a point at the limit.
  fred1=abs((xxyy3(k1,1)-f1)/f1);
  fred2=abs((xxyy3(k1,1)-f2)/f2);
  if fred1<0.02; xxyy3(k1,1)=f1; end
  if fred2<0.02; xxyy3(k1,1)=f2; end
  xxyy3(k1,2)=s1*10^((-xxyy2(k1,2)+x1y1a(2))/srs_length);
end
%    
%  Output 
%
if(itype==1)
    srs_out=xxyy3;
    disp(' ')
    disp(' The output filename is "srs_out" ')
    nn=length(srs_out(:,1));
    for(i=1:nn-1)
        if(srs_out(i,1)>=srs_out(i+1,1))
%%            disp(' Frequency error. Sorting data to correct error.');
              temp=srs_out(i,1);
              srs_out(i,1)=srs_out(i+1,1);
              srs_out(i+1,1)=temp;
        end
    end
    B = sortrows(srs_out);
    srs_out=B;
else
    psd_out=xxyy3;
    disp(' ')
    disp(' The output filename is "psd_out" ')   
    nn=length(psd_out(:,1));   
    for(i=1:nn-1)
        if(psd_out(i,1)>=psd_out(i+1,1))
 %%           disp(' Frequency error.  Sorting data to correct error.');
              temp=psd_out(i,1);
              psd_out(i,1)=psd_out(i+1,1);
              psd_out(i+1,1)=temp;
        end
    end 
    B = sortrows(psd_out);
    psd_out=B;  
%
    for(i=2:nn)
        if(psd_out(i,1)<=psd_out(i-1,1))
           psd_out(i,1)=1.005*psd_out(i-1,1); 
        end
    end
%
end
%
if choice == 1 
%
    for(i=1:nn)
        if(itype==1)
            fprintf(fid,'%12.5g %12.5g \n',srs_out(i,1),srs_out(i,2));
        else
            fprintf(fid,'%12.5g %12.5g \n',psd_out(i,1),psd_out(i,2));
        end
    end
%
   disp(' ');
   disp(' Output filename: ');
   out1=sprintf('  %s',TH_filename);
   disp(out1);
   fclose(fid);
end
%
disp(' ');
disp(' Clear figure? ');
disp(' 1=yes 2=no ');
icf=input(' ');
if(icf==1)
   clf('reset')
end