你帮我重新组织一下这个matlab的函数sfr_calcSFR3_Accutance,不要改变原来的算法,可以重新命名,去掉多余 的内容
function [output mtf_perROI] = sfr_calcSFR3_Accutance(input)
a = input.dataArray;
try
accutanceCutoff = input.accutanceCutoff;
catch
accutanceCutoff = 0.5;
end
% Parameters
del = 1; % frequency is given in cy/pixel
nbin = 4;
[nlin npix ncol] = size(a);
a = double(a);
[a, nlin, npix, rflag] = rotatev2(a); %based on data values
% Force 'positive' edge
fil1 = [0.5 -0.5];
fil2 = [0.5 0 -0.5];
tleft = mean(mean(a(:, 1:5, 1)));
tright = mean(mean(a(:, npix-4:npix,1)));
if tleft>tright;
fil1 = [-0.5 0.5];
fil2 = [-0.5 0 0.5];
end
% Edge estimation variables
fitme = zeros(ncol, 3);
slout = zeros(ncol, 1);
loc = zeros(ncol, nlin);
win1 = ahamming(npix, (npix+1)/2); % Symmetric window
% Smoothing window for first part of edge location estimation
for color=1:ncol; % Loop for each color
c = deriv1(a(:,:,color), nlin, npix, fil1);
% compute centroid for derivative array for each line in ROI. NOTE WINDOW array 'win'
for n=1:nlin
loc(color, n) = centroid( c(n, 1:npix )'.*win1) - 0.5; % -0.5 shift for FIR phase
end;
[fitme(color,1:2) fitstats] = findedge(loc(color,:), nlin);
place = zeros(nlin,1);
for n=1:nlin;
place(n) = fitme(color,2) + fitme(color,1)*n;
win2 = ahamming(npix, place(n));
loc(color, n) = centroid( c(n, 1:npix )'.*win2) -0.5;
end;
[fitme(color,1:2) fitstats] = findedge(loc(color,:), nlin);
end; % End of loop for each color
midloc = zeros(ncol,1);
for i=1:ncol;
slout(i) = - 1./fitme(i,1); % slope is as normally defined in image coods.
if rflag==1, % positive flag it ROI was rotated
slout(i) = - fitme(i,1);
end;
% evaluate equation(s) at the middle line as edge location
midloc(i) = fitme(i,2) + fitme(i,1)*((nlin-1)/2);
end
% Limit number of lines to integer(npix*line slope as per ISO algorithm
nlin1 = round(floor(nlin*abs(fitme(1,1)))/abs(fitme(1,1)));
% disp(['Integer cycle lines: ',num2str(nlin1)])
a = a(1:nlin1, :, 1:ncol);
vslope = fitme(1,1);
slope_deg= 180*atan(abs(vslope))/pi;
del2=0;
%Correct sampling inverval for sampling parallel to edge
delfac = cos(atan(vslope));
del = del*delfac;
del2 = del/nbin;
nn = floor(npix *nbin);
mtf = zeros(nn, ncol);
mtf1 = zeros(nn,ncol);
mtf2 = zeros(nn,ncol);
nn2 = floor(nn/2) + 1;
freq = zeros(nn, 1);
% dcorr = fir2fix(nn2, 3); % dcorr corrects SFR for response of FIR filter
dcorr = fir2fix(nn2, 0.99999); %1.001 dcorr corrects SFR for response of FIR filter
for n=1:nn;
freq(n) = nbin*(n-1)/(del*nn);
end;
freqlim = 1;
if nbin == 1;
freqlim = 2;
end
win = ahamming(nbin*npix,(nbin*npix+1)/2); % centered Hamming window
% initialize 2 ESFs for asymmetry score (AS) calculation
point1 = zeros(nn,1); point2 = zeros(nn,1);
% ESF -> PSF -> fft -> mtf
% Asymmetry Score calculations
for color=1:ncol % Large SFR loop for each color record
% project and bin data in 4x sampled array
[point, ~] = project(a(:,:,color), loc(color, 1), fitme(color,1), nbin);
% compute first derivative via FIR (1x3) filter fil
c = deriv1(point', 1, nn, fil2);
c = c';
% shift array so it is centered
mid = centroid(c);
temp = cent(c, round(mid));
c = temp;
% apply window (symmetric Hamming)
c = win.*c;
% Transform, scale and correct for FIR filter response
temp = abs(fft(c, nn));
mtf(1:nn2, color) = temp(1:nn2)/temp(1);
mtf(1:nn2, color) = mtf(1:nn2, color).*dcorr;
% smooth ESF, then calculate PSF to get max transition location
% (split location to separate left and right side of ESF)
point_padded = padarray(point,2,'replicate','both');
point_smoothed = conv(point_padded,ones(1,5)/5,'valid');
c_point_smoothed = deriv1(point_smoothed', 1, nn, fil2);
fil2nn = numel(fil2);
[~,split_loc]=max(c_point_smoothed(fil2nn:end));
split_loc = split_loc + floor(fil2nn/2);
% rotate ESF about split location
point_rot270 = rot270(point,split_loc);
% generate 2 ESFs for left and right side of original ESF
point1(1:split_loc) = point(1:split_loc);
point1(split_loc+1:end) = point_rot270(split_loc+1:end);
point2(split_loc+1:end) = point(split_loc+1:end);
point2(1:split_loc) = point_rot270(1:split_loc);
% calculate PSF->FFT->MTF of point1
c1 = deriv1(point1', 1, nn, fil2); c1 = c1'; mid = centroid(c1);
temp = cent(c1, round(mid)); c1 = temp; c1 = win.*c1;
temp = abs(fft(c1, nn));
mtf1(1:nn2, color) = temp(1:nn2)/temp(1);
mtf1(1:nn2, color) = mtf1(1:nn2, color).*dcorr;
% calculate PSF->FFT->MTF of point2
c2 = deriv1(point2', 1, nn, fil2); c2 = c2'; mid = centroid(c2);
temp = cent(c2, round(mid)); c2 = temp; c2 = win.*c2;
temp = abs(fft(c2, nn));
mtf2(1:nn2, color) = temp(1:nn2)/temp(1);
mtf2(1:nn2, color) = mtf2(1:nn2, color).*dcorr;
end; % color=1:ncol
%mtfAccutanceSum = [];
mtfResultsRGBY_ny2 = zeros(1,ncol);
%AS_RGBY_SampledFreqMean = zeros(1,ncol);
for i = 1:ncol
mtfResultsRGBY_ny2(i) = interp1(freq, mtf(:,i), 0.25);
end
% mtf plot validation sampling
mtf_valid=interp1(freq, mtf,0:0.0001:1);
mtf_perROI.mtf_valid=mtf_valid;
mtf_perROI.point=point;
mtf_perROI.c=c;
mtf_perROI.a=a;
mtf_perROI.fitme=fitme;
mtf_perROI.nlin=nlin;
mtf_perROI.slope_deg=slope_deg;
mtf_perROI.ny2=mtfResultsRGBY_ny2;
output = [mtfResultsRGBY_ny2, slope_deg];
end
% [data] = ahamming(n, mid) Generates asymmetrical Hamming window
% array. If mid = (n+1)/2 then the usual symmetrical Hamming array
% is returned
% n = length of array
% mid = midpoint (maximum) of window function
% data = window array (nx1)
function [data] = ahamming(n, mid)
data = zeros(n,1);
wid1 = mid-1;
wid2 = n-mid;
wid = max(wid1, wid2);
multiplier = (1:n)';
data = 0.54 + 0.46*cos(pi*((multiplier-mid)/wid));
end
% [b] = deriv1(a, nlin, npix, fil)
% Computes first derivative via FIR (1xn) filter
% Edge effects are suppressed and vector size is preserved
% Filter is applied in the npix direction only
% a = (nlin, npix) data array
% fil = array of filter coefficients, eg [-0.5 0.5]
% b = output (nlin, npix) data array
function [b] = deriv1(a, nlin, npix, fil)
b = zeros(nlin, npix);
nn = length(fil);
for i=1:nlin;
temp = conv(fil, a(i,:));
b(i, nn:npix) = temp(nn:npix); %ignore edge effects, preserve size
b(i, nn-1) = b(i, nn);
end
end
% function [loc] = centroid(x)
% Returns centroid location of a vector
% x = vector
% loc = centroid in units of array index
function [loc] = centroid(x)
n = 1:length(x);
sumx = sum(x);
if sumx < 1e-4;
loc = 0;
else loc = sum(n*x)/sumx;
end
end
%[a, nlin, npix, rflag] = rotatev2(a) Rotate edge array vertical
% Rotate array so that the edge feature is in the vertical orientation
% Test based on array values not dimensions.
% a = input array(npix, nlin, ncol)
% nlin, npix are after rotation if any
% flag = 0 no roation, = 1 rotation was performed
function [a, nlin, npix, rflag] = rotatev2(a)
dim = size(a);
nlin = dim(1);
npix = dim(2);
a = double(a);
% Select which color record, normally the second (green) is good
if length(dim) == 3;
mm = 2;
else
mm =1;
end
nn = 3; % Limits test area. Normally not a problem.
%Compute v, h ranges
testv = abs(mean(a(end-nn,:,mm))-mean(a(nn,:,mm)));
testh = abs(mean(a(:,end-nn,mm))-mean(a(:,nn,mm)));
rflag =0;
if testv > testh
rflag =1;
a = rotate90(a);
temp=nlin;
nlin = npix;
npix = temp;
end
end
% [slope, int] = findedge(cent, nlin)
% fit linear equation to data, written to process edge location array
% cent = array of (centroid) values
% nlin = length of cent
% slope and int are from the least-square fit
function [slope, int] = findedge(cent, nlin)
index=[0:nlin-1];
[slope int] = polyfit(index, cent, 1); % x = f(y)
return
end
%[out] = rotate90(in, n) 90 degree counterclockwise rotations of matrix
% in = input matrix (n,m) or (n,m,k)
% n = number of 90 degree rotation
% out = rotated matrix
function out = rotate90(in, n)
if nargin < 2;
n = 1; %Default to 90* rotation
end
nd = ndims(in);
if nd < 1
error('input to rotate90 must be a matrix');
return
end
for i = 1:n
out = r90(in);
in = out;
end
return
end
function [out] = r90(in)
[nlin, npix, nc] = size(in);
temp = zeros (npix, nlin);
out = zeros (npix, nlin, nc);
for c = 1: nc;
temp = in(:,:,c);
temp = temp.';
out(:,:,c) = temp(npix:-1:1, :);
end
out = squeeze(out);
end
% [correct] = fir2fix(n, m);
% Correction for MTF of derivative (difference) filter
% n = frequency data length [0-half-sampling (Nyquist) frequency]
% m = length of difference filter
% correct = nx1 MTF correction array (limited to a maximum of 10)
function [correct] = fir2fix(n, m)
if nargin<2;
del=0;
end
correct = ones(n, 1);
m=m-1;
scale = 1;
for i = 2:n;
correct(i) = abs((pi*i*m/(2*(n+1))) / sin(pi*i*m/(2*(n+1))));
correct(i) = 1 + scale*(correct(i)-1);
if correct(i) > 10; % Note limiting the correction to the range [1, 10]
correct(i) = 10;
end;
end
end
% [point, status] = project(bb, loc, slope, fac) Projects data
% Projects the data in array bb along the direction defined by
% npix = (1/slope)*nlin.
% Data is accumulated in 'bins' that have a width (1/fac) pixel.
% The smooth, supersampled one-dimensional vector is returned.
% bb = input data array
% slope and loc are from the least-square fit to edge
% fac = oversampling (binning) factor, default = 4
% point = output vector
function [point, status] = project(bb, loc, slope, fac)
status =1;
[nlin, npix]=size(bb);
if nargin<4
fac = 4 ;
end;
big = 0;
nn = npix *fac ;
% smoothing window
win = ahamming(nn, fac*loc(1, 1));
slope = 1/slope;
offset = round( fac* (0 - (nlin - 1)/slope ) );
del = abs(offset);
if offset>0
offset=0;
end
barray = zeros(2, nn + del+100);
ling_array=[];
% Projection and binning
for n=1:npix;
for m=1:nlin;
x = n-1;
y = m-1;
ling = ceil((x - y/slope)*fac) + 1 - offset;
ling_array=[ling_array;ling,m, n];
barray(1,ling) = barray(1,ling) + 1;
barray(2,ling) = barray(2,ling) + bb(m,n);
end;
end;
point = zeros(nn,1);
start = 1+round(0.5*del);
% Check for zero counts
nz =0;
for i = start:start+nn-1;
if barray(1, i) ==0
nz = nz +1;
status = 0;
if i==1;
barray(1, i) = barray(1, i+1);
else
barray(1, i) = (barray(1, i-1) + barray(1, i+1))/2;
end
end
end
% Combine in single edge profile, point
for i = 0:nn-1;
% try
point(i+1) = barray(2, i+start-nz)/ barray(1, i+start-nz); % Jan 2016 Aaron Cong edit (may need to evaluate if hi-> low condition has possibility of zero in the begin of sequence)
% catch
% point(i+1) = barray(2, i+start)/ barray(1, i+start);
% end
end
if nz>0
fprintf('roi error nz= %d \n',nz)
end
% debugging mode
% figure;plotyy(1:length(barray(2,:)),barray(2,:),1:length(barray(2,:)),barray(1,:));title('barray noise review')
return
end
% Usage: [b] = cent(a, center)
% Matlab function cent, shift of one-dimensional
% array, so that a(center) is located at b(round((n+1)/2).
% Written to shift a line-spread function array prior to
% applying a smoothing window.
% a = input array
% center = location of signal center to be shifted
% b = output shifted array
function [b] = cent(a, center)
n = length(a);
b = zeros(n, 1);
mid = round((n+1)/2);
del = round(center - mid);
if del > 0;
for i = 1:n-del;
b(i) = a(i + del);
end;
elseif del < 1;
for i = -del+1:n;
b(i) = a(i + del);
end;
else b = a;
end;
end
% Usage: b = rot270(a,loc)
% Matlab function cent edited
% rotate a 270 degrees about [a(loc), loc]
% a = input array
% loc = location to rotate a about
function b = rot270(a,loc)
a = flipud((a(loc)*2-a));
n = length(a);
b = zeros(n, 1);
del = round(n-2*loc+1);
if del > 0;
for i = 1:n-del;
b(i) = a(i + del);
end;
b(i+1:end) = b(i);
elseif del < 1;
for i = -del+1:n;
b(i) = a(i + del);
end;
b(1:-del) = b(-del+1);
else b = a;
end;
end
最新发布
本文介绍了一种使用动态规划(DP)筛选法求解特定问题的C++实现方案。通过三个不同的序列逐步更新并筛选最小值，最终得到第n项的具体数值。此方法适用于解决涉及最小值更新的问题。
扫一扫
10
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
