SFR sfrmat3 Code Learning | Eddie’s Blog

Document

图像分辨率测试ISO12233 - 2017中文翻译_jgw2008的专栏-CSDN博客

译者序：目前的摄像头分辨率的测试，大多遵循ISO 12233标准，最近下载一份英文版的文档，和大家一起分享，仅供学习使用。目录 ISO 12233 Third Edition (2017-01) Foreword(前言) Introduction（介绍） Purpose（目的） Methods for measuring SFR and resolution - Selection rationale and guidance 测量SFR和分辨率的方法.选择原理和指导 Photography - Electronic still picture imaging - Resolution and spatial frequency responses 2 Normative references （规范性引用文件） 3 Terms and definitions（术语和定义） 3.1 addressable photoelements 可寻址光电元件 3.3 cycles

https://blog.csdn.net/jgw2008/article/details/116519404

ISO Flowchart of e-SFR measurement algorithm

Diagram depicting the key steps of the e-SFR algorithm

Read in images

Lens shading correction

解决由于lens的光学特性，由于镜头对于光学折射不均匀导致的镜头周围出现阴影的情况。

ISP-镜头阴影校正（LSC）

概述介绍镜头阴影校正（Lens Shading Correction）是为了解决由于lens的光学特性，由于镜头对于光学折射不均匀导致的镜头周围出现阴影的情况。 shad...

https://www.jianshu.com/p/19f7feb0a6b1

Code

matlab


[XX, YY] = meshgrid( X, Y);
[ ~,rho] = cart2pol(XX,YY); % Transform Cartesian coordinates to polar or cylindrical
% rho is the distance from the coordinate origin (0,0) (top-left corner of the image) and theta is the line rotation angle in radians
% Create LSC map
lsc_surf = 1 - fall_off .* rho / max(rho(:));

img_lsc = uint16(round(lsc(img_raw, cx, cy, config.sfrnv.lsc.fall_off)));

python


# Create coordinate map
(XX, YY) = np.meshgrid(X, Y)
rho=np.zeros((height, width))
rho = np.sqrt(XX**2 + YY**2)
# rho is the distance from the coordinate origin (0,0) (top-left corner of the image) and theta is the line rotation angle in radians
# Create LSC map
lsc_surf = 1 - fall_off * rho / (np.max(rho))

img_lsc = np.uint16(round5(lsc(img_raw, cx, cy, config["sfrnv"]["lsc"]["fall_off"])))

numpy.meshgrid()理解_lllxxq141592654的博客-CSDN博客

网格点是什么？坐标矩阵又是什么鬼？看个图就明白了：图中，每个交叉点都是，描述这些网格点的坐标的矩阵，就是坐标矩阵。再看个简单例子 A,B,C,D,E,F是6个网格点，坐标如图，如何用矩阵形式（坐标矩阵）来批量描述这些点的坐标呢？答案如下：这就是坐标矩阵 --横坐标矩阵中的每个元素，与纵坐标矩阵中对应位置元素，共同构成一个点的完整坐标。如B点坐标。下面可以自己用matplotlib来试一试，输出就是上边的图如果对matplotlib不熟悉，可能只知道用一列横坐标（线性代数中的1维列向量），一列纵坐标生成（两者元素个数相等）一些点。但是实际上，给matplotlib的坐标信息是矩阵也是可以的，只要横纵坐标的尺寸一样。都会按照对应关系生成点。但是有需要注意的地方，按照矩阵给坐标点信息，matplotlib会把横坐标矩阵中，每一列对应的点当做同一条线。举个例子，把上面的代码 plot的 linestyle=''删掉，或者变成 linestyle='-'（这个操作把图的线型改为默认状态），就会发现A-D是连接的，B-E是连接的，C-F是连接的，也即，会认为你输入的是3条线，如图作为练习，自己试着生成如下结果提示：线型等关键字参数设置可用如下代码 plt.plot(x, y, marker='.', # 点的形状为圆点 markersize=10, # 点设置大一点，看着清楚 linestyle='-.') # 线型为点划线答案 import numpy as np import matplotlib.pyplot as plt x = np.array([[0, 1, 2, 3],

https://blog.csdn.net/lllxxq141592654/article/details/81532855

Split raw image into 4 sub-channels

Bayer

R, Gr (Green Row), B, Gb (Green Balance)

Undemosaiced

In this document, we sometimes refer to RAW files from commercial cameras or development systems as Camera RAW to distinguish them from Bayer RAW files, which are standard monochrome image files that contain undemosaiced (Bayer) data.

Camera RAW≠Bayer RAW

Undemosaiced~=4 channels

RAW Files

The unprocessed digital output of an image sensor is called RAW image data. In this document, we sometimes refer to RAW files from commercial cameras or development systems as Camera RAW to distinguish them from Bayer RAW files, which are standard monochrome image files that contain undemosaiced (Bayer) data.

https://www.imatest.com/docs/raw/

Code

matlab


% Build color planes
img_chan = raw2chan(img_lsc);

python


#  Build color planes
img_chan = np.uint16(raw2chan(img_lsc))

Calculate fixed sizes for the chart

Code (Load config)

matlab


square_size = config.sfrnv.roi.square_size;           % square size as percent of image diagonal
square_names = config.sfrnv.roi.square_names';         % names of the squares in the chart
square_distances = config.sfrnv.roi.square_distances';     % distances of the center of the squares as ratio of DFOV
square_angles = config.sfrnv.roi.square_angles';        % angles of the squares in FOV - 0 degs is east
square_rotations = config.sfrnv.roi.square_rotations';     % angle of rotation of each square patch
patch_size = config.sfrnv.roi.sb_patch_size;         % size of ROI as percent of image diagonal
min_patch_size = config.sfrnv.roi.min_patch_size;        % minimum ROI size in pixels
patch_names = config.sfrnv.roi.sub_patch_names';      % names of the four ROIs used in each square 
patch_angles = config.sfrnv.roi.sub_patch_angles';     % sub patch angles desired for unrotated square                   
freqs = config.sfrnv.roi.freqs';                % freqs for sampeling MTF

[by, bx, channels] = size(img_chan) ;       
chart_diag = sqrt(bx ^ 2 + by ^ 2);                 
cen_x = bx / 2 + 0.5;  cen_y = by / 2 + 0.5;    
square_size_px = square_size * chart_diag;      % estimate square size in pixels
% calculate ideal patch distances from the beginning of the image
id_patch_dist = 0.5 .* chart_diag .* square_distances;       
id_patch_x = id_patch_dist .* cosd(square_angles) + cen_x;   
id_patch_y = cen_y - id_patch_dist .* sind(square_angles);  
patch_dist_px = chart_diag .* square_size .* 0.5;   % calculate the ROI distance from square center in pixels
patch_size_px = chart_diag .* patch_size .* 0.5; % calculate the ROI size in pixels

python


# square size as percent of image diagonal
square_size = config["sfrnv"]["roi"]["square_size"]
# names of the squares in the chart
square_names = config["sfrnv"]["roi"]["square_names"]
# distances of the center of the squares as ratio of DFOV
square_distances = config["sfrnv"]["roi"]["square_distances"]
square_distances = np.asarray(square_distances)
# angles of the squares in FOV - 0 degs is east
square_angles = config["sfrnv"]["roi"]["square_angles"]
square_angles = np.asarray(square_angles)
# angle of rotation of each square patch
square_rotations = config["sfrnv"]["roi"]["square_rotations"]
square_rotations = np.asarray(square_rotations)
# size of ROI as percent of image diagonal
patch_size = config["sfrnv"]["roi"]["sb_patch_size"]
# minimum ROI size in pixels
min_patch_size = config["sfrnv"]["roi"]["min_patch_size"]
# names of the four ROIs used in each square
patch_names = config["sfrnv"]["roi"]["sub_patch_names"]
# sub patch angles desired for unrotated square
patch_angles = config["sfrnv"]["roi"]["sub_patch_angles"]
square_angles = np.asarray(square_angles)
freqs = config["sfrnv"]["roi"]["freqs"]  # freqs for sampeling MTF

chart_diag = np.sqrt(bx ** 2 + by ** 2)
cen_x = bx / 2 + 0.5
cen_y = by / 2 + 0.5
square_size_px = square_size * chart_diag      # estimate square size in pixels
# calculate ideal patch distances from the beginning of the image
id_patch_dist = 0.5 * chart_diag * square_distances
id_patch_x = id_patch_dist * np.cos(square_angles*np.pi/180) + cen_x
id_patch_y = cen_y - id_patch_dist * np.sin(square_angles*np.pi/180)
# calculate the ROI distance from square center in pixels
patch_dist_px = chart_diag * square_size * 0.5
patch_size_px = chart_diag * patch_size *  0.5  # calculate the ROI size in pixels

Threshold image and find regions

Code

matlab


% Create binary image
level  = graythresh(img_chan(:, :, 1));                   
bw_img = imbinarize(img_chan(:, :, 1), level);             
bw_img = logical(1 - bw_img); 

% find the areas of all of the blobs in the image
stats  = regionprops(bw_img, 'Area', 'Centroid');   
idx = [stats.Area] > ((square_size_px * 0.75) ^ 2) & ... 
    [stats.Area] < ((square_size_px * 1.25) ^ 2) ;  % filter area +/- 25% of expected area                                                
centroids = cat(1, stats(idx).Centroid);    % get centroids

python


ret, threshold = cv2.threshold(img_chan[0, ], 0, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)
bw_img = np.asarray(threshold)
# bw_img = (255-bw_img).astype(bool)
bw_img = 255-bw_img

# find the areas of all of the blobs in the image
# connectivity = 2代表8连通
lable = measure.label(bw_img, background=0, connectivity=2)
props = measure.regionprops(lable)
# retval, labels, stats, centroids = cv2.connectedComponentsWithStats(threshold, connectivity=8)
centroids = np.empty((1, 2))
for prop in props:
		if prop.area > ((square_size_px * 0.75) ** 2) and prop.area < ((square_size_px * 1.25) ** 2):
				centroids = np.append(centroids, np.reshape(np.asarray(prop.centroid), (1, 2)), axis=0)
# (x, y) 轴位置；第一行空，需删除
centroids = centroids[:, [1, 0]]
centroids = np.delete(centroids, 0, axis=0)

Find the coordinate of centroids

stats.Area\in((square\_size\_px * 0.75) ^ 2, (square\_size\_px * 1.25) ^ 2)

Calculate the SFR

💡

for square; square_names=[c, r, t, l, ......]; (for each square)

💡

for patch; patch_name=[t, b, r, l]; (for each ROI)

💡

for chan; channels(for each colour channel)

Code (Get the scope of the ROI)

matlab


% Find the closest centroid to ideal patch
[~, idx] = min(sqrt((id_patch_x(square) - centroids(:, 1)) .^ 2 + (id_patch_y(square) - centroids(:, 2)) .^ 2)); 
    
% Using the closest found patch do calculation to extract the right patches 
patch_x = patch_dist_px .* cosd(patch_angles + square_rotations(square))+ centroids(idx, 1);  
patch_y = centroids(idx, 2) - patch_dist_px .* sind(patch_angles + square_rotations(square)); 
    
% calculate the expected unrotated bounds of the patch
if square<=5
		left   = round(patch_x - patch_size_px);  
		right  = round(patch_x + patch_size_px);
		top    = round(patch_y + patch_size_px);
		bottom = round(patch_y - patch_size_px);
else 
		left   = round(patch_x - patch_size_px+10);  
		right  = round(patch_x + patch_size_px-10);
		top    = round(patch_y + patch_size_px-5);
		bottom = round(patch_y - patch_size_px+5);
end

python


# Find the closest centroid to ideal patch
idx = np.argmin(np.sqrt((id_patch_x[square,] - centroids[:, 0]) ** 2 + (id_patch_y[square,] - centroids[:, 1]) ** 2))
 
# Using the closest found patch do calculation to extract the right patches 
patch_x = patch_dist_px * np.cos((patch_angles + square_rotations[square,])*np.pi/180)+ centroids[idx, 0]
patch_y = centroids[idx, 1] - patch_dist_px * np.sin((patch_angles + square_rotations[square,])*np.pi/180)

# calculate the expected unrotated bounds of the patch
if square<5:
		left   = round5(patch_x - patch_size_px)
		right  = round5(patch_x + patch_size_px)
		top    = round5(patch_y + patch_size_px)
		bottom = round5(patch_y - patch_size_px)
else: 
		left   = round5(patch_x - patch_size_px+10)
		right  = round5(patch_x + patch_size_px-10)
		top    = round5(patch_y + patch_size_px-5)
		bottom = round5(patch_y - patch_size_px+5)

The "Threshold image and find regions" have already found the ideal coordinate of each square.

The real square coordinate has the minimum distance between the ideal one and itself.

Sketch map of finding the ROI (named as patch & sub_patch in MATLAB codes)

sfrmat3

matlab


[status, dat, e, fitme, esf, nbin, del2] = sfrmat3(io, del, weight, a, oename)

[status, dat2, e2, fitme, esf2, nbin, del2] = sfrmat3(1, 1, [0.3, 0.6, 0.1], sfr_patch);

dat = computed sfr data

oename optional name of oecf LUT file containing 3xn or 1xn array. The camera optoelectronic conversion function (OECF). 伽马校正部分用于校正非线性光输出（与之对应的是显示器的信号输入）。光输入强度和电信号输出强度之间的关系称为相机光电转换函数（OECF），为图像捕获设备提供伽马校正曲线形状。

Size or Shape of "a" (img_patch)

matlab


[nlin npix ncol] = size(a)

python


(ncol, nlin, npix) = a.shape

ncol=1

rotatev2

Code

matlab


nn=3;
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);

[a, nlin, npix, rflag] = rotatev2(a); %based on data values

python


nn=3
testv = abs(np.mean(a[mm,-nn-1,:])-np.mean(a[mm,nn-1,:]))
testh = abs(np.mean(a[mm,:,-nn-1])-np.mean(a[mm,:,nn-1]))
rflag=0
if testv>testh:
		rflag=1
		a=np.rot90(a, 1)

a, nlin, npix, rflag = rotatev2(a, nlin, npix) # based on data values

If the difference between 2 rows > 2 columns: rotate to vertical

Edge contrast

tleft, tright, test

Loop for each color

The first estimate of edge slope & offset

Code

matlab


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;
% clear c
fitme(color,1:2) = findedge(loc(color,:), nlin);

python


if color == 0:
		c = deriv1(a, nlin, npix, fil1)
else:
		c = deriv1(a[color, :, :], nlin, npix, fil1)
# compute centroid for derivative array for each line in ROI. NOTE WINDOW array 'win'
loc = np.zeros((ncol, nlin))
for n in range(nlin):
# -0.5 shift for FIR phase
		loc[color, n] = centroid(c[n, 0:npix].transpose()*win1[:, 0]) - 0.5
fitme[color, 0:2] = findedge(loc[color, :], nlin)

deriv1 (Compute derivative)? (27 May 2020 updated to use 'same' conv option)

matlab


temp = conv(fil, a(i,:));
b(i, nn:npix) = temp(nn:npix);    %ignore edge effects, preserve size
b(i, nn-1) = b(i, nn);

c = deriv1(a(:,:,color), nlin, npix, fil1);

python


temp=np.convolve(fil1, a[i,:])
c[i, nn-1:npix]=temp[nn-1:npix] #ignore edge effects, preserve size
c[i, nn-2]=c[i, nn-1]

c = deriv1(a, nlin, npix, fil1)

Apply convolution for each row of the image matrix.

💡

The "derivative" it calculated is quite likely! Whereas, I don't handle how to calculate the derivative of a discrete sequence!

centroid (Compute centroid for derivative array for each line in ROI)

matlab


n = 1:length(x);
sumx = sum(x);
if sumx < 1e-4;
		loc = 0;
else loc = sum(n*x)/sumx;

loc(color, n) = centroid( c(n, 1:npix )'.*win1) - 0.5; % -0.5 shift for FIR phase

python


n=range(len(x))
n=np.asarray(n)
n=n+1
sumx=np.sum(x)
if sumx<0.0001:
		loc=0
else:
		loc=np.sum(np.dot(n, x))/sumx # 
    return loc

loc[color, n] = centroid(c[n, 0:npix].transpose()*win1[:, 0]) - 0.5

💡

Different from ISO document.

Each row (r) of the edge spread image is an estimate of the camera edge spread function (ESF). Each of these ESFs is differentiated to form its discrete line spread function (LSF). The position of the centroid (C) of each r LSF is determined along the continuous variable x, where x has the range (1, X) [see Formula (D.3)].

💡

It really works! If you don't believe it, use x=[3, 2, 3] to make a test.

💡

Though it is in operation, I have no idea about the weird algorithm.

The final estimate of edge slope and offset

Code

matlab


fitme(color,1:2) = findedge(loc(color,:), nlin);
place = zeros(nlin,1);
for n=1:nlin;
	place(n) =fitme(color,1)*n + fitme(color,2);
	win2 = ahamming(npix, place(n));
	loc(color, n) = centroid( c(n, 1:npix )'.*win2) -0.5;
end
fitme(color,1:2) = findedge(loc(color,:), nlin);

sfrmat3.m

python


fitme[color, 0:2] = findedge(loc[color, :], nlin)
place = np.zeros((nlin, 1))
for n in range(nlin):
		place[n] = fitme[color, 1]+fitme[color, 0]*(n) # !!! (n+1) !!!
		win2 = ahamming(npix, place[n])
		loc[color, n] = centroid(c[n, 0:npix].transpose()*win2[:, 0])-0.5
fitme[color, 0:2] = findedge(loc[color, :], nlin)

sfrmat3.py

findedge 线性拟合(Done in the first estimate)

matlab


index=[0:nlin-1];
[slope int] = polyfit(index, cent, 1);

[slope, int] = findedge(cent, nlin)

slope从高次幂开始

cent=slope(1)*index + slope(2)

place(n) =fitme(color,1)*n + fitme(color,2)

python


index=range(nlin)
slope = np.polyfit(index, cent, 1)

findedge(cent, nlin)

slope从高次幂开始

cent=slope[0]*index + slope[1]

place(n) =fitme[color,0]*n + fitme[color,1]

x = fitme(color, 1) y + fitme(color, 2) % so the slope is the inverse of the one that you may expect

hamming: place= mid = midpoint (maximum) of the window function

If mid = (n+1)/2 then the usual symmetric Hamming

centroid (Mentioned above)

💡

根据centroid线性拟合；根据线性拟合再计算centroid；根据线性的centroid再计算slope

Edge slope judgment

Code

matlab


vslope = fitme(1,1);
slope_deg= 180*atan(abs(vslope))/pi;
%disp(['Edge angle: ',num2str(slope_deg, 3),' degrees'])
if slope_deg < 3.5 % Unit: Degree
    % beep, warndlg(['High slope warning ',num2str(slope_deg,3),' degrees'], 'Watch it!')
    disp(' ** WARNING: High edge slope.');
end

python


vslope = fitme[0, 0]
slope_deg = 180*math.atan(abs(vslope))/math.pi
if slope_deg < 3.5: # Unit: Degree
		print(" ** WARNING: High edge slope.")

Correction (Derivative correction)

Code (fir2fix)

matlab


del2=0;
if oldflag ~= 1;
    %Correct sampling inverval for sampling parallel to edge
    delfac = cos(atan(vslope));
    del = del*delfac;
    del2 = del/nbin;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nn =   floor(npix *nbin);
mtf =  zeros(nn, ncol);
nn2 =  floor(nn/2) + 1;

if oldflag ~=1;
    %    disp('Derivative correction')
    dcorr = fir2fix(nn2, 3);    % dcorr corrects SFR for response of FIR filter
end

python


del2 = 0
if oldflag != 1:
		# Correct sampling inverval for sampling parallel to edge
		delfac = math.cos(math.atan(vslope))
		del1 = del1*delfac
		del2 = del1/nbin

nn = math.floor(npix*nbin)
mtf = np.zeros((nn, ncol))
n2 = math.floor(nn/2)+1

if oldflag != 1:
		# % dcorr corrects SFR for response of FIR filter
		dcorr = fir2fix(nn2, 3)

fir2fix

matlab


m=3;
correct(i) = abs((pi*i*m/(2*(n+1))) / sin(pi*i*m/(2*(n+1))));
if correct(i) > 10;  % Note limiting the correction to the range [1, 10]
   correct(i) = 10;

[correct] = fir2fix(n, m); dcorr = fir2fix(nn2, 3)

python


m=3
correct[i]=np.abs((math.pi*(i+1)*m/(2*(n+1)))/math.sin(math.pi*(i+1)*m/(2*(n+1))))
if correct[i]>10:
		correct[i]=10

fir2fix(n, m)

D ( k) is the correction for the frequency response of the discrete derivative used to derive the point spread function from the edge spread function,

Project (shift the image data along the edge direction to the edge of the ROI)

Code

matlab


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);
% 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;
        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;
for i = 0:nn-1;
    point(i+1) = barray(2, i+start)/ barray(1, i+start);
end

% project and bin data in 4x sampled array [point, status] = project(a(:,:,color), loc(color, 1), fitme(color,1), nbin); [point, status] = project(bb, loc, slope, fac)

python


nn=npix*fac
slope=1/slope
offset= int(Decimal(str(fac*(0-(nlin-1)/slope))).quantize(Decimal('0'), rounding=ROUND_HALF_UP))
del1=abs(offset)
if offset>0:
		offset=0
barray=np.zeros((2, nn+del1+100))
# Projection and binning
for n in range(npix):
		for m in range(nlin):
				x=n-1+1
				y=m-1+1
				ling=math.ceil((x-y/slope)*fac)+1-offset
				barray[0, ling-1]+=1
				barray[1, ling-1]=barray[1, ling-1]+bb[m, n]
point=np.zeros((nn, 1))
start=1+int(Decimal(str(0.5*del1)).quantize(Decimal('0'), rounding=ROUND_HALF_UP))
# Check for zero counts
nz=0
for i in range(start, start+nn):
		if barray[0, i]==0:
				nz+=1
				status=0
		if i==1:
				barray[0, i]=barray[0, i+1]
		else:
				barray[0, i]=(barray[0, i-1]+barray[1, i+1])/2
for i in range(nn):
    point[i]=barray[1, i+start-1]/barray[0, i+start-1]

point, status = project(a, loc[color, 0], fitme[color, 0], nbin)

project 到线性拟合的方向上

Projects the data in array bb along the direction defined by $npix = (\frac{1}{slope})*nlin.$ Used by sfrmat11 and sfrmat2 functions. Data is accumulated in 'bins' that have a width $\cfrac{1}{fac}$ pixel.

Compute the one-dimensional derivative of the edge-spread function

Code

matlab


% compute first derivative via FIR (1x3) filter fil
c = deriv1(point', 1, nn, fil2);
c = c';

python


# compute first derivative via FIR (1x3) filter fil
c = deriv1(point.transpose(), 1, nn, fil2)
c = c.transpose()

Shift to center the line spread function (LSF) vector

Code

matlab


mid = centroid(c);
temp = cent(c, round(mid));              % shift array so it is centered
c = temp;

[b] = cent(a, center)

python


mid = centroid(c)
temp = cent(c, int(Decimal(str(mid)).quantize(Decimal('0'), rounding=ROUND_HALF_UP)))  # % shift array so it is centered
c = temp

cent(a, center) ... return b

cent: shift of one-dimensional array, so that

a(center)

is located at

b(round(\cfrac{n+1}{2}))

(round(\cfrac{n+1}{2}))

是横坐标的中心，centroid 也返回横坐标，但与纵坐标数值有关

Apply a Hamming window to the LSF vector

Code

matlab


% apply window (symmetric Hamming)
c = win.*c;

python


# apply window (symmetric Hamming)
c = win[:, 0]*c[:, 0]

w_{Hm}(n)=[0.54-0.46\cos(\cfrac{2\pi n}{N-1})]R_N(n)

Compute the discrete Fourier transform (DFT) of the windowed, binned LSF vector; Compute the modulus of the complex DFT array

matlab


temp = abs(fft(c, nn));

python


temp = np.abs(np.fft.fft(c[:, ], nn))

Normalize the modulus vector by the zero-frequency value (first element of the array) to obtain the E-SFR

matlab


mtf(1:nn2, color) = temp(1:nn2)/temp(1);

python


mtf[:nn2, color] = temp[:nn2]/temp[0]

Correct the E-SFR for the discrete derivative response

Code

matlab


if oldflag ~=1;
		mtf(1:nn2, color) = mtf(1:nn2, color).*dcorr;
end

python


if oldflag != 1:
		mtf[:nn2, color] = mtf[:nn2, color]*dcorr[:, 0]

Report E-SFR result

Code

matlab


del=1;
%Correct sampling inverval for sampling parallel to edge
delfac = cos(atan(vslope));
del = del*delfac;
nn =   floor(npix *nbin);
for n=1:nn;
    freq(n) = nbin*(n-1)/(del*nn);
end;

dat = zeros(nn2out, ncol+1);
for i=1:nn2out;
    dat(i,:) = [freq(i), mtf(i,:)];
end

patch_results(chan, :) = interp1(dat2(:, 1), dat2(:, 2), freqs);

python


del1=1
delfac = math.cos(math.atan(vslope))
del1 = del1*delfac
freq = np.zeros((nn, 1))
for n in range(nn):
		freq[n] = nbin*(n-1+1)/(del1*nn) # !!! (n-1+1) !!!

dat = np.zeros((nn2out, ncol+1))
for i in range(nn2out):
		dat[i, :] = np.asarray([freq[i, 0], mtf[i, :]], dtype=np.double)

patch_results[chan, :]=np.interp(freqs, dat2[:, 0], dat2[:, 1])

dat=[ freq, mtf ]

Derived metrics

Construct the square_results matrix from patch_results

matlab


patch_results(chan, :) = interp1(dat2(:, 1), dat2(:, 2), freqs);

square_results(patch,:,:) = patch_results;

dat(i,:) = [freq(i), mtf(i,:)]

python


patch_results[chan, :]=np.interp(freqs, dat2[:, 0], dat2[:, 1])

square_results[0, patch, :]=patch_results[:, 0].T # 0.25
square_results[1, patch, :]=patch_results[:, 1].T # 0.5

dat[i, :] = np.asarray([freq[i, 0], mtf[i, :]], dtype=np.double)

根据特殊算出的采样频率 freq 作为样本点（横坐标），对 mtf 值（对应值）进行插值，在

freq=[0.25, 0.5]

查询

Screenshots of MATLAB & Python

patch_result: each color channel

$square\_result: each patch name ("sub\_patch\_names": ["t", "b", "r", "l"])$

Condense and assess patch results

matlab


h_idx = contains(patch_names, 'l') | contains(patch_names, 'r');    % index of edges for horizontal MTF
v_idx = contains(patch_names, 't') | contains(patch_names, 'b');    % index of edges for vertical MTF
sfr = mean(square_results(:, :, f))';

sfr_h = mean(square_results(h_idx, :, f))';
sfr_v = mean(square_results(v_idx, :, f))';

python


sfr=np.mean(square_results[f, :, :], axis = 0).T  # axis=0，计算每一列的均值

sfr_h=np.mean(square_results[f, 2:4, :], axis=0).T  # axis=0，计算每一列的均值
sfr_v=np.mean(square_results[f, 0:2, :], axis=0).T  # axis=0，计算每一列的均值

"sub_patch_names": ["t", "b", "r", "l"]

💡

Each SFR data of a color channel is equal to the mean of t, b, r, l

Calculate tilt and falloff

tilt

matlab


% Get all outer positon data
for o = 1:length(outer_name)    % for each out square
    h_data = results.data.(strcat(outer_name{o}, "_", freq_output, '_h'));
    v_data = results.data.(strcat(outer_name{o}, "_", freq_output, '_v'));
    outer_data(:, o) = mean([h_data v_data], 2);    % average all edges of square
end
outer_data = mean(outer_data(gr_index | gb_index, :));  % average green channels

tilt = max(outer_data) - min(outer_data);

python


# Get all outer position data
outer_data0=np.zeros((channels, 4))
results_data2=np.asarray(results_data2)
results_data2=np.nan_to_num(results_data2, nan=0)
outer_data0[:, 0]=np.mean(np.hstack((results_data2[:, 20].reshape(4, 1), results_data2[:, 21].reshape(4, 1))), axis=1) # ne_05
outer_data0[:, 1]=np.mean(np.hstack((results_data2[:, 24].reshape(4, 1), results_data2[:, 25].reshape(4, 1))), axis=1) # nw_05
outer_data0[:, 2]=np.mean(np.hstack((results_data2[:, 28].reshape(4, 1), results_data2[:, 29].reshape(4, 1))), axis=1) # sw_05
outer_data0[:, 3]=np.mean(np.hstack((results_data2[:, 32].reshape(4, 1), results_data2[:, 33].reshape(4, 1))), axis=1) # se_05
results_data2=results_data2.tolist()
outer_data=np.mean(np.vstack((outer_data0[0, :], outer_data0[3, :])), axis=0) # 压缩行，对各列求均值

tilt=np.max(outer_data)-np.min(outer_data)

Screenshots of Python & MATLAB

for outer_name: outer 的 05Data h & v，列拼接

for outer_name: 对行求 mean，形成 outer_data

for outer_name: outer_data 取 Gr 和 Gb 行，行拼接

求列 mean

$tilt = max(outer\_data) - min(outer\_data)$ ;

falloff

matlab


center_data = results.data.(strcat(center_name, "_", freq_output));
center_data = mean(center_data(gr_index | gb_index));   % average green channels

falloff = center_data - mean(outer_data);

python


results_data[0].append(sfr[0])
results_data[1].append(sfr[1])
results_data[2].append(sfr[2])
results_data[3].append(sfr[3])

center_data=(results_data[0][1]+results_data[-1][1])/2

falloff=center_data-np.mean(outer_data)

fetch 05Data Gr & Gb from center

calculate the mean

falloff = center_data - mean(outer_data);