基于Matlab可以用来计算MTF值数据文件需要为Excel格式

function sfrcalc(varargin) % Calculate SFR or MTF from an Excel-formatted text data file from ImageJ. % The frequency scale is assumed to be logarithmic. % Command line input data [default]: % filename. minimum frequency (lp/mm) [2], maximum frequency (lp/mm) [200], % Picture Height [24 mm]; optional % Example: sfrcalc inputfile.tif 2 80 15.1 % EOS-10D 70-200 f/4 @ 70 mm, f/8 inpfile = varargin{1}; % data file from ImageJ (convert from cell.) %fileprt = strrep(inpfile,'_','\_'); % To print correctly in plots. minlp = 2; maxlp = 200; picht = 24; % Defaults for min, max lp/mm and picture height. if (nargin>=3) minlp = str2num(varargin{2}); maxlp = str2num(varargin{3}); % Scale lp/mm min, max. end if (nargin>=4) picht = str2num(varargin{4}); % Picture height mm. end descr = input('Enter system description: ','s'); [fid,MESSAGE] = fopen(inpfile, 'r'); % open data file for read if (fid<0) error(MESSAGE); end [arrinp,l_ar] = fscanf(fid,'%f %f'); fclose(fid); close all; % Data in the first column is unreliable-- limited precision. Throw it out. arr2 = arrinp(2:2:l_ar); l_ar = length(arr2); dif1 = diff(arr2); fprintf('%d data points read in\n',l_ar); cmap = zeros(256,3); % Color map: gray scale. Apparently gamma corrected. gamma = 1; cmap(1:256,1) = ([0:255]'/255).^(1/gamma); cmap(1:256,2) = cmap(1:256,1); cmap(1:256,3) = cmap(1:256,1); yim = ones(25,1)*arr2'; % Input image for check. yim(1,1:l_ar) = 0; yim(1:25,1) = 0; yim(1:25,l_ar) = 0; % Black line on top, sides. % Detect peaks. d1 = (dif1>0); % 1 or 0. d1 = (d1(2:l_ar-1) ~= d1(1:l_ar-2)); % 1 at peaks; 0 elsewhere. d2 = d1.*(abs(arrflt(1:l_ar-1))>thresh*maxflt); % Detection threshold peak_loc = find(d2==1); % integer peak locations. peak_loc = find(d1==1)+1; % integer peak locations. n_peaks = length(peak_loc); % number of detected peaks. fprintf('%d detected peaks\n',n_peaks); peak_ampl = arr2(peak_loc); % End peak detection. if (nargin<3) figure; plot(peak_loc, peak_ampl); grid on; zoom on; title(descr); axis tight; p1 = get(gcf,'Position'); p1(4) = p1(4)-60; p1(3) = p1(3)-80; set(gcf,'Position',p1); fprintf('Scale values not entered. End program.\n'); else summ = peak_ampl(1:n_peaks-1)+peak_ampl(2:n_peaks); resp1 = abs(diff(peak_ampl))./summ; resp1 = resp1/max(resp1); % Use equations from lpratio = maxlp/minlp; xax1 = (0:l_ar-1)/(l_ar-1); % Input data x-axis. lpmm1 = minlp*lpratio.^xax1; figure; p1 = get(gcf,'Position'); p1(4) = p1(4)-60; p1(3)= p1(3)-80; plot(lpmm1,arr2); grid on; zoom on; title(descr); set(gcf,'Position',p1); xlabel('Line pairs per millimeter (lp/mm)'); ylabel('Input data'); xax2 = peak_loc(2:n_peaks)/(l_ar-1); % X-axis for MTF; lpmm2 = minlp*lpratio.^xax2; figure; p1 = get(gcf,'Position'); % Composite figure. p1(2) = p1(2)-120; p1(4) = p1(4)+80; p1(3) = p1(3)-140; set(gcf,'Position',p1); p2 = get(gcf,'PaperPosition'); p2(2) = p2(2)-1; p2(4) = p2(4)+2; set(gcf,'PaperPosition',p2); subplot('Position',[.11 .92 .85 .03]); image(yim); colormap(cmap); axis off; title(descr); subplot('Position',[.11 .63 .85 .29]); plot(arr2); grid on; zoom on; axis tight; xlabel(['Log scale from ' num2str(minlp) ' to ' num2str(maxlp) ' lp/mm']); ylabel('Input data'); grid on; zoom on; ax1 = axis; a1 = ax1(2)-ax1(1); a2 = ax1(4)-ax1(3); text(ax1(1)+.4*a1,ax1(3)+.07*a2,fileprt); subplot('Position',[.11 .08 .85 .46]); % subplot(2,1,2); plot(lpmm2,resp1); grid on; zoom on; xlabel('Line pairs per millimeter (lp/mm)'); ylabel('Spatial frequency response'); ax1 = axis; a1 = ax1(2)-ax1(1); a2 = ax1(4)-ax1(3); if (nargin>=4) % Do plot in Line widths/Picture height =LW/PH. %text(ax1(1)+.35*a1,ax1(3)+.92*a2,['To obtain Line widths per picture height']); %text(ax1(1)+.35*a1,ax1(3)+.84*a2,['(lw/ph), where ph = ' num2str(picht) 'mm,']); %text(ax1(1)+.35*a1,ax1(3)+.76*a2,['multiply by 'num2str(2*picht)]); figure; set(gcf,'Position',p1); % New plot of LW/PH. subplot('Position',[.11 .92 .85 .03]); image(yim); colormap(cmap); axis off; title(descr); subplot('Position',[.11 .63 .85 .29]); plot(arr2); grid on; zoom on; axis tight; xlabel(['Log scale from ' num2str(minlp) ' to ' num2str (maxlp) ' lp/mm']); ylabel('Input data'); grid on; zoom on; ax1 = axis; a1 = ax1(2)-ax1(1); a2 = ax1(4)-ax1(3); text(ax1(1)+.4*a1,ax1(3)+.07*a2,fileprt); subplot('Position',[.11 .08 .85 .46]); % subplot(2,1,2); plot(lpmm2*2*picht,resp1); grid on; zoom on; xlabel(['Line widths per picture height (lw/ph) for ph =' num2str(picht) ' mm']); ylabel('Spatial frequency response'); ax1 = axis; a1 = ax1(2)-ax1(1); a2 = ax1(4)-ax1(3); text(ax1(1)+.42*a1,ax1(3)+.85*a2,['To obtain lp/mm,divide by ' num2str(2*picht)]); end end