数字图像处理MATLAB学习笔记（三）

数字图像处理MATLAB学习笔记（三）

Image Restoration and Reconstruction 图像复原

1. A Model of the Image Degradation/Restoration Process

• Image Degradation Process：对图像进行退化和噪声添加处理。
  • Input：Imgae：记为$f(x,y)$
  • Degradation function：记为$H[f(x,y)]$
    • 如果H是线性、空间不变的过程：
      • $H[f(x,y)]=h(x,y)\odot f(x,y)$
      • 其对应的频域下的Fourier Transformer是：$G(u,v)=H(u,v)F(u,v)+N(u,v)$
  • Nosie：记为$\eta (x,y)$
  • Output：degraded image：记为$g(x,y)=H[f(x,y)]+\eta (x,y)$
• Image Restoration Process：对drgraded image进行“复原”。实际上，是通过对g进行最大估计，使得输出尽可能接近于原输入。
  • Input：degraded image：记为$g(x,y)$
  • Output：restored image：记为$\hat{f}(x,y)$

两个函数：

• Optical Transfer Function(OTF)：光传递函数，Fourier Analysis中的degradation function，即$H(u,v)$
• Point Spread Function(PSF)：点扩散函数，空间域下degradation function的空间大小，即$h(x,y)$

matlab中实现两个函数之间转换的方法是：function otf2psf和function psf2otf。

在退化函数的线性和空间不变时，退化函数实际上是卷积过程，有时候被称为“PSF与图像卷积”；同理，restoration process被称为反卷积。

2. Noise Models

2.1 Adding Noise to Image with Function imnoise

内置函数imnoise主要是用于给图像添加噪声干扰，语法为

g = imnoise(f, type, parameters)

参数信息：

• f：输入图像
• type：添加噪声类型，主要有：gaussian、lacalvar、salt & pepper、speckle、poisson
• parameters：根据type输入参数

有一下几种类型：

g = imnoise(f, 'gaussian', m , var) % 均值m、方差var的高斯噪声。默认m=0，var=0.01
g = imnoise(f, 'localvar', V)       % 均值0、局部方差V的高斯噪声。V和f大小相同，包含每个点的理想方差值
g = imnoise(f, 'localvar', image_intensity, var) % 均值0、局部方差var(f的灰度值函数)的高斯噪声。var和image_intensity是大小相同的向量，且image_intensity包含归一化灰度值(range [0,1])。
g = imnoise(f, 'salt & pepper', d)  % 密度d的噪声。有d*numel(f)个像素点会收到噪音污染，默认密度0.05
g = imnoise(f, 'speckle', var)      % 方程式g=f+n.*f的均值n为0、方差为var的均匀分布的随机噪声。默认var为0.04
g = imnoise(f, 'poisson')           % 泊松分布噪声。

2.2 Generating Spatial Random Noise with a Specified Distribution使用给定分布产生空间随机噪声

Probability Distributed Function(PDF)：概率密度函数

Cumulative Distributed FUnction(CDF)：累计分布函数

概率论中一个重要的结果：对于(0,1)区间均匀分布的随机变量w，可以通过求解以下方程得到具有规定CDF的F的随机变量z：$z=F^{-1}(w)$，等价为寻找方程$w=F(z)$的解。

Example：Using uniform random numbers to generate random numbers with a specified distribution

假设一个Rayleigh CDF的随机数的形式如下：
$$F(z)=\{\begin{array}{ll}1-e^{-(z-a{)}^2/b}& \text{for}z\ge a\\ 0& \text{for}z<a\end{array}$$
当$b<0$时，我们可以找到方程的解为$1-e^{-(z-a{)}^2/b}=w$或$z=a+\sqrt{-b\ln (1-w)}$，即随机数产生器方程式。

在matlab中，我们通过如下表达式产生随机数组R：

>> R = a + sqrt(b * log(1 - rand(M, N)));

随机变量的PDF、CDF和Genrator：

Name	PDF	Mean and Variance	CDF	Generator
Uniform	$p(z)=\{\begin{array}{ll}\frac{1}{b-a}& \text{if}0\le z\le b\\ 0& \text{otherwise}\end{array}$	$\small{m=\frac{a+b}{2}}\$，${\small\sigma 2=\frac{(b-a)2}{12}}\$	$F(z)=\{\begin{array}{ll}0& z<a\\ \frac{z-a}{b-a}& a\le z\le b\\ 1& z>b\end{array}$	MATLAB function rand
Gaussian	$\small{p(z)=\frac{1}{\sqrt{2\pi }b}exp(-(z-a)2/2b2)}\$	$m=a$，${\sigma }^2=b^2$	$\small{F(z)=\int^z_{\infty}p(v)dv}\$	MATLAB function randn
Lognormal	$$\begin{gathered} p(z)=\frac{1}{\sqrt{2\pi }bz}exp(-(\ln (z)-a{)}^2/2b^2) \\ z>0 \end{gathered}$$	$m=e^{a+\frac{b^2}{2}}$，${\sigma }^2=[e^{b^2}-1]e^{2a+b^2}$	$\small{F(z)=\int^z_0p(v)dv}\$	$z=e^{bN(0,1)+a}$
Rayleigh	$p(z)=\{\begin{array}{ll}\frac{2}{b}(z-a)exp(-(z-a{)}^2/b)& z\ge a\\ 0& z<a\end{array}$	$m=a+\sqrt{\pi b/2}$，${\small\sigma^2=\frac{b(4-\pi)}{4}}\$	$F(z)=\{\begin{array}{ll}1-e^{-(z-a{)}^2/b}& z\ge a\\ 0& z<a\end{array}$	$z=a+\sqrt{-b\ln [1-U(0,1)]}$
Exponential	$p(z)=\{\begin{array}{ll}a{e}^{-az}& z\ge 0\\ 0& z<0\end{array}$	$\small{m=\frac{1}{a}}\$，${\small\sigma2=\frac{1}{a2}}\$	$F(z)=\{\begin{array}{ll}1-e^{-az}& z\ge 0\\ 0& z<0\end{array}$	$\small{z=-\frac{1}{a}\ln{[1-U(0,1)]}}\$
Erlang	$$\begin{gathered} p(z)=\frac{a^b{z}^{b-1}}{(b-1)!}{e}^{-az} \\ z\ge 0 \end{gathered}$$	$\small{m=\frac{b}{a}}\$，${\small\sigma2=\frac{b}{a2}}\$	$$\begin{gathered} F(z)=[1-e^{-az}{\sum }_{n=0}^{b-1}\frac{(az{)}^n}{n!}] \\ z\ge 0 \end{gathered}$$	$z=E_1+E_2+E_3+\dots$($E$是参数为a的exponential rand)
Salt & Pepper	$p(z)=\{\begin{array}{ll}{P}_p& z=0\text{(pepper)}\\ P_s& z=2^n-1\text{(salt)}\\ 1-(P_p+P_s)& z=k(0<k<2^n-1)\end{array}$	$$\begin{gathered} m=(0)P_p+k(1-P_p-P_s)+(2^n-1)P_s \\ {\sigma }^2=(0-m{)}^2{P}_p+(k-m{)}^2(1-P_p-P_s)+(2^n-1-m{)}^2{P}_s \end{gathered}$$	$F(z)=\{\begin{array}{ll}0& z<0\\ P_p& 0\le z<k\\ 1-P_s& k\le z<2^n-1\\ 1& 2^n-1\le z\end{array}$	MATLAB function rand with some additional logic

这里需要注意一下Salt&Pepper Noise Function，被污染的总的概率求解为$P=P_s+P_p$，也就是salt noise和pepper noise的概率和。椒盐噪音模型产生的图像实际上具有三个值：原值、0白、255黑。结合表中的式子，n即代表图像的存储数值bit大小，也就是说，当n为8bit时，$P_p$的概率为0，$P_s$的概率为255，$k$的概率为$1-P_p-P_s$。

Matlab function rand：随机生成一个大小为MxN的数组，并遵循(0,1)区间内的均匀分布。如果参数数量为1，默认生成MxM大小的数组；若为0，默认随机生成一个数。

A = rand(M, N)

Matlab function randn：随机生一个遵循均值为0、单位方差的Gaussian正态分布的MxN的数组。参数规则与rand一致

A = randn(M, N)

Matlab function find：数组内索引查询。和python里的find实际上差不多，返回的事数组的索引数组。这里传参可以是数组，也可以是某种逻辑规则。

output = find(param)

我们用以上函数，定义一个imnoise2函数

function R = imnoise2(type, varargin)
%IMNOISE2 Generates an array of random numbers with specified PDF
%   R = IMNOISE2 (TYPE, M, N, A, B) generates an array, R, of size 
%   M-by-N, whose elements are random numbers of the speci fied TYPE 
%   with parameters A and B. If only TYPE is included in the
%   input argument list, a single random number of the specified
%   TYPE and default parameters shown below is generated. If only
%   TYPE, M, and N are provided, the default parameters shown below
%   are used. If M = N = 1, IMNOISE2 generates a single random
%   number of the specified TYPE and parameters A and B.
% 
%   Valid values for TYPE and parameters A and B are :
%   
%   'uniform'  Uniform random numbers in the interval (A, B) .
%              The default values are (0, 1) .
%   'gaussian' Gaussian random numbers with mean A and standard
%              deviation B. The default values are A = 0,
%              B=1.
%   'salt & pepper' Salt and pepper numbers of amplitude 0 with 
%                   probability Pa = A, and amplitude 1 with
%                   probability Pb = B. The default values are Pa = 
%                   Pb = A = B = 0.05. Note that the noise has
%                   values 0 (with probability Pa = A) and 1 (with
%                   probability Pb = B), so scaling is necessary if
%                   values other than 0 and 1 are required. The
%                   noise matrix R is assigned three values. If
%                   R(x, y) = 0, the noise at (x, y) is pepper
%                   (black) . If R(x, y) = 1, the noise at (x, y) is
%                   salt (white) . If R(x, y) = 0.5, there is no
%                   noise assigned to coordinates (x, y) .
%   'lognormal' Lognormal numbers with offset A and shape
%               parameter B. The defaults are A = 1 and B =
%               0.25.
%   'rayleigh'  Rayleigh noise with parameters A and B. The 
%               default values are A = 0 and B = 1.
%   'exponential' Exponential random numbers with parameter A.
%                 The default is A = 1.
%    'erlang'   Erlang (gamma) random numbers with parameters A
%               and B. B must be a positive integer. The
%               defaults are A = 2 and B = 5. Erlang random
%               numbers are approximated as the sum of B
%               exponential r andom numbers.

% Set defaults
[M, N, a, b] = setDefaults(type, varargin{:});

% Begin processing. Use lower(type) to protect against input being
% capitalized
switch lower(type)
    case 'uniform'
        R = a + (b - a) * rand(M, N);
    case 'gaussian'
        R = a + b * randn(M, N);
    case 'salt & pepper'
        R = saltpepper(M, N, a, b);
    case 'lognormal'
        R = exp(b*randn(M, N) + a);
    case 'rayleigh'
        R = a + (-b * log(1 - rand(M, N))) .^ 0.5;
    case 'exponential'
        R = exponential(M, N, a);
    case 'erlang'
        R = erlang(M, N, a, b);
    otherwise
        error('Unknow distribution type');
end

%------------------------------------------------------------------
function R = saltpepper(M, N, a, b)
% Check to make sure that Pa + Pb is not > 1
if (a + b) > 1
    error('The sum Pa + Pb must not exceed 1.');
end
R(1:M, 1:N) = 0.5;
% Generate an M-by-N array of uniformly-distributed random numbers
% in the range (0,1). Then, Pa*(M*N) of them will have values <= a
% The coordinates of these points we call 0 (pepper noise).
% Similarly, Pb*(M*N) points will have values in the range > a & <=
% (a+b). These we call 1 (salt noise).
X = rand(M, N);
R(X <= a) = 0;
u = a + b;
R(X > a & X <= u) = 1;

%------------------------------------------------------------------
function R = exponential(M, N, a)
if a <= 0
    error('Parameter a must be positive for exponential type.');
end

k = -1/a;
R = k*log(1 - rand(M, N));

%------------------------------------------------------------------
function R = erlang(M, N, a, b)
if (b ~= round(b) || b <= 0)
    error('Param b must be a positive integer for Erlang');
end
k = -1/a;
R = zeros(M, N);
for j = 1:b
    R = R + k*log(1 - rand(M, N));
end

%------------------------------------------------------------------
function varargout = setDefaults(type, varargin)
varargout = varargin;
P = numel(varargin);
if P < 4
    % Set default b.
    varargout{4} = 1;
end
if P < 3
    % Set default a.
    varargout{3} = 0;
end
if P < 2
    % Set default N.
    varargout{2} = 1;
end
if P < 1
    % Set default M.
    varargout{1} = 1;
end
if (P <= 2)
    switch type
        case 'salt & pepper'
            % a = b = 0.05
            varargout{3} = 0.05;
            varargout{4} = 0.05;
        case 'lognormal'
            % a = 1, b = 0.25
            varargout{3} = 1;
            varargout{4} = 0.25;
        case 'exponential'
            % a = 1
            varargout{3} = 1;
        case 'erlang'
            % a = 2, b = 5
            varargout{3} = 2;
            varargout{4} = 5;
    end
end

我们通过随机生成各种类型的随机数100000个，绘制成图像如图所示：

这里的function hist是绘制柱状图的函数，其中一个参数***bins***是柱子的宽度。

2.3 Periodic Noise

一般的周期噪声指的是图像获取过程中的电器和电动机械产生的干扰。这里我们讨论一个2-D的正弦波周期噪声信号：
$$r(x,y)=Asin[2\pi u_0(x+B_x)/M+2\pi v_0(y+B_y)/N]$$
其中，$A$是振幅，$x=0,1,\dots ,M-1$、$y=0,1,\dots ,N-1$，$u_0$和$v_0$决定x、y方向上的正弦波频率，$B_x$和$B_y$是关于远点的相移。

其DFT为：
$$R(u,v)=j\frac{AMN}{2}[e^{-j2\pi (u_0{B}_x/M+v_0{B}_y/N)}\delta (u+u_0,v+v_0)-e^{j2\pi (u_0{B}_x/M+v_0{B}_y/N)}\delta (u-u_0,v-v_0)]$$
其中，$u=0,1,\dots ,M-1$、$v=0,1,\dots ,N-1$。实际上这是一对分别位于$(u+u_0,v+v_0)$和$(u-u_0,v-v_0)$的共轭复数，也就是说，第一项为0，除了$u=-u_0,v=-v_0$；第二项为0，除了$u=u_0,v=v_0$。

代码实现为：

function [r, R, S] = imnoise3(M, N, C, A, B)
%IMNOISE3 Generates periodic noise
%   [r, R, S] = IMNOISE3 (M, N, c, A, B) , generates a spatial
%   sinusoidal noise pattern, r, of size M-by-N, its Fourier
%   transform, R, and spectrum, S. The remaining parameters are:
% 
%   C is a K-by-2 matrix with K pairs of frequency domain
%   coordinates (u, v) that define the locations of impulses in the
%   frequency domain. The locations are with respect to the
%   frequency rectangle center at [floor (M/2) + 1, floor (N/2) + 1] ,
%   where the use of function floor is necessary to guarantee that
%   all values of (u, v) are integers, as required by all Fourier
%   formulations in the book. The impulse locations are specified as
%   integer increments with respect to the center. For example, if M
%   N = 512, then the center is at (257, 257) . To specify an
%   impulse at (280, 300) we specify the pair (23, 43) ; i.e.; 257 +
%   23 = 280, and 257 + 43 = 300. Only one pair of coordinates is
%   required for each impulse. The conjugate pairs are generated
%   automatically.
% 
%   A is a 1-by-K vector that contains the amplitude of each of the
%   K impulse pairs. If A is not included in the argument, the
%   default used is A = ONES(1, K) . B is then automatically set to
%   its default values (see next paragraph) . The value specified
%   for A(j) is associated with the coordinates in C(j, :)
% 
%   B is a K-by-2 matrix containing the Bx and By phase components
%   for each impulse pair. The default value for B is zeros (K, 2) .

% Process input parameters
K = size(C, 1);
if nargin < 4
    A = ones(1, K);
end
if nargin < 5
    B = zeros(K, 2);
end

% Generate R.
R = zeros(M, N);
for j = 1:k
    % Bases on the equation for R(u,v), we know that the first term
    % of R(u,v) associated with a sinusoid is 0 unless u = -u0 and
    % v = -v0
    u1 = floor(M/2) + 1 - C(j, 1);
    v1 = floor(N/2) + 1 - C(j, 2);
    R(u1, v1) = 1i*M*N*(A(j)/2) * exp(-1i*2*pi*(C(j,1)*B(j,1)/M ...
                                + C(j,2)*B(j,2)/N));
    % Conjugate. The second term is zero unless u = u0 and v = v0
    u2 = floor(M/2) + 1 + C(j, 1);
    v2 = floor(N/2) + 1 + C(j, 2);
    R(u2, v2) = -1i*M*N*(A(j)/2) * exp(1i*2*pi*(C(j,1)*B(j,1)/M ...
                                + C(j,2)*B(j,2)/N));
end

% Compute the spectrum and spatial sinusoidal pattern.
S = abs(R);
r = real(ifft2(ifftshift(R)));

Example: Using function imnoise3：生成几幅不同的噪音图像及其光谱图

>> C = [0 64; 0 128; 32 32; 64 0; 128 0; -32 32];
>> [r, R, S] = imnoise3(512, 512, C);
>> subplot(2,3,1);
>> imshow(S, []);
>> title('Spectrum of specified impulses')
>> subplot(2,3,4);
>> imshow(r, []);
>> title('Corresponding sine noise patter in the spatial domain')
>> C2 = [0 32; 0 64; 16 16; 32 0; 64 0; -16 16];
>> [r2, R2, S2] = imnoise3(512, 512, C2);
>> subplot(2,3,2);
>> imshow(S2, []);
>> title('Spectrum of specified impulses')
>> subplot(2,3,5);
>> imshow(r2, []);
>> title('Corresponding sine noise patter in the spatial domain')
>> C3 = [6 32; -2 2];
>> A = [1 5];
>> [rf, Rf, Sf] = imnoise3(512, 512, C3, A);
>> [re, Re, Se] = imnoise3(512, 512, C3);
>> subplot(2,3,3);
>> imshow(re, []);
>> title('Higher-frequenct sine wave')
>> subplot(2,3,6);
>> imshow(rf, []);
>> title('Lower-frequenct sine wave')

2.4 Estimating Noise Parameters

周期噪声估计的典型方法是傅里叶频谱。

对于spatial domain noise的估计可以将问题转变成为估计sample image的mean和variance，并用估计值求解对应PDF的参数a和b。

令$z_i$表示一幅图像的灰度级的离散随机变量，其中$i=0,1,\dots ,L-1$。并令$p(z_i)$表示相应的归一化直方图，也是灰度PDF的近似值。

直方图的中心矩(central moments)：${\mu }_n={\sum }_{i=0}^{L-1}(z_i-m{)}^{n}p(z_i)$

n是中心矩的阶(order)，m是均值：$m={\sum }_{i=0}^{L-1}{z}_{i}p(z_i)$

假设直方图是归一化的，它的所有分量之和为1，结合上式，可以得到${\mu }_0=1$且${\mu }_1=0$。

方差是二阶矩：${\mu }_2={\sum }_{i=0}^{L-1}(z_i-m{)}^{2}p(z_i)$

在MATLAB中，我们使用function statmoments计算mean和n-order的central moments，语法如下：其中p为直方图向量，n为计算的矩的数量。输出v是标准化的n阶矩，unv包含了v相同的矩，但是用原始值所在区间内的数据进行计算。

[v, unv] = statmomments(p, n)

源码如下：

function [v, unv] = statmoments(p, n)
%STATMOMENTS Computes statistical central moments of image histogram.
%   [W, UNV] = STATMOMENTS(P, N) computes up to the Nth statistical central
%   moment of a hi stogram whose components are in vector P. The length of
%   P must equal 256 or 65536.
% 
%   The program outputs a vector V with V (1) = mean, V (2) = variance, V
%   (3) = 3rd moment, ... V (N) = Nth central moment. The random variable
%   values are normalized to the range [0, 1], so all moments also are in
%   this range .
% 
%   The program also outputs a vector UNV containing the same moments as V,
%   but using un-normalized random variable values (e.g., 0 to 255 if
%   length(P) = 2~8) . For example, if length(P) = 256 and V(1) = 0.5, then
%   UNV(1) would have the value UNV(1) = 127.5 (half of the [0 255] range).

Lp = length(p);
if (Lp ~= 256) && (LP ~= 65536)
    error('P must be a 256- or 65536-element vector');
end
G = Lp - 1;

% Make sure the histogram has unit area, and convert it to a column vector.
p = p/sum(p);
p = p(:);

% Form a vector of all the possible values of the random variable
z = 0:G;

% Now normalize the z's to the range [0, 1].
z = z./G;

% The mean
m = z*p;

% Center random variables about the mean
z = z - m;

% Compute the central moments.
v = zeros(1, n);
v(1) = m;
for j = 2:n
    v(j) = (z.^j)*p;
end

if nargout > 1
    % Compute the uncentralized moments
    unv = zeros(1, n);
    unv(1) = m.*G;
    for j = 2:n
        unv(j) = ((z*G).^j)*p;
    end
end

噪声参数通常直接由带噪声的图像或一组图像来估计。选取一块尽量和背景一样无特征的区域，作为ROI(Region of Interest)。在MATLAB中，使用function roipoly，语法如下

B = roipoly(f, c, r)

其中，参数c和r是相应有序多边形的坐标（列先被指定）。输出B是大小与f相同的一个二值图像，ROI区域内值为1，其余区域值为0。

作用：roipoly通常被作为模版限制ROI的操作。

根据roipoly的性质，编写ROI的直方图绘制函数function histroi：

function [p, npix] = histroi(f, c, r)
%HISTROI Computes the histogram of an ROI in an image.
%   [P, NPIX] = HISTROI(F, C, R) computes the histogram, P, of a
%   polygonal region of interest (ROI) in image F. The polygonal 
%   region is defined by the column and row coordinates of its
%   vertices, which are specified (sequentially) in vectors C and R,
%   respectively. All pixels of F must be >= 0. Parameters NPIX is the
%   number of pixels in the polygonal region.

% Generate the binary mask image.
B = roipoly(f, c, r);

% Compute the histogram of the pixels in the ROI
p = imhist(f(B));

% Obtain the number of pixels in the ROI if requested in the output.
if nargout > 1
    npix = sum(B(:));
end

3 Restoration in the Presence of Noise Only-Spatial Filtering

如果省略degradation function，退化过程就变为：
$$g(x,y)=f(x,y)+\eta (x,y)$$
这样的话，降噪较好的方法就是Spatial Filter（😓好嘛，又回到了滤波器…）

3.1 Spatial Noise Filters

各种各样的Spatial Filters(m和n分别表示行和列的步长)：

Name	Equation	Comments
Arithmetic Mean	${\hat{f}(x,y)=\frac{1}{mn}\sum_{(s,t\in S_{xy})}g(s,t)}\$	MATLAB函数实现采用*w=fspecial(‘average’,[m,n])和f=imfilter(g,w)*实现
Geometric Mean	$\hat{f}(x,y)=[\prod _{(s,t)\in S_{xy}}g(s,t){]}^{\frac{1}{mn}}$	非线性，MATLAB实现使用gmean
Harmonic Mean调和均值	$${\hat{f}(x,y)=\frac{mn}{\displaystyle\sum_{(s,t)\in S_{xy}}g(s,t)^{-1}}}\$$	非线性，MATLAB实现使用harmean
Contraharmonic Mean反调和均值	${\hat{f}(x,y)=\frac{\displaystyle\sum_{(s,t)\in S_{xy}}g(s,t)^{Q+1}}{\displaystyle\sum_{(s,t)\in S_{xy}}g(s,t)^{Q}}}\$	非线性，MATLAB实现使用charmean
Median	${\hat{f}(x,y)=\mathop{\text{median}}{(s,t)\in S{xy}}{g(s,t)}}\$	MATLAB实现使用medfilt2(g,[m n],‘symmetric’)
Max	${\hat{f}(x,y)=\max_{(s,t)\in S_{xy}}{g(s,t)}}\$	MATLAB实现使用imdilate(g,ones(m,n))
Min	${\hat{f}(x,y)=\min_{(s,t)\in S_{xy}}{g(s,t)}}\$	MATLAB实现使用imerode(g,ones(m,n))
Midpoint	${\hat{f}(x,y)=\frac{1}{2}\bigg[\max_{(s,t)\in S_{xy}}{g(s,t)}+\min_{(s,t)\in S_{xy}}{g(s,t)}\bigg]}\$	max与min滤波器的0.5sum
Alpha-trimmed mean	${\hat{f}(x,y)=\frac{1}{mn-d}\sum_{(s,t)\in S_{xy}}g_r(s,t)}\$	删除d/2个最低和最高的像素点，$g_r$表示剩下的mn-d个点，使用alphatrim函数实现。

根据给中filter的定义，实现function spfilt：

function f = spfilt(g, type, varargin)
%SPFILT Performs linear and nonlinear spatial filtering
%   F = SPFILT (G, TYPE, M, N, PARAMETER) performs spatial filtering
%   оf image G using a TYPE filter of size M-by-N. Valid cal1s to
%   SPFILT are as follows :
% 
%     F = SPFILT (G, 'amean', M, N)        Arithmetic mean filtering.
%     F = SPFILT (G, 'gmean', M, N)        Geometric mean filtering.
%     F = SPFILT (G, 'hmean' , M, N)       Harmonic mean filtering.
%     F = SPFILT (G, 'chmean' , M, N, Q)   Contraharmonic mean
%                                          filtering of order Q. The
%                                          default Q is 1.5.
%     F = SPFILT (G, 'median', M, N)       Median filtering.
%     F = SPFILT (G, 'max', M, N)          Max filtering.
%     F = SPFILT (G, 'min' , M, N)         Min filtering.
%     F = SPFILT (G, 'midpoint' , M, N)    Midpoint filtering.
%     F = SPFILT (G, 'atrimmed' , M, N, D) Alpha-trimmed mean
%                                          filtering. Parameter D must
%                                          be a nonnegative even
%                                          integer; its default value
%                                          is 2.
% 
%   The default values when only G and TYPE are input are M = N = 3,
%   Q =1.5, and D=2.

[m, n, Q, d] = processInputs(varargin{:});

% Do the filtering
switch type
    case 'ameam'
        w = fspecial('average', [m n]);
        f = imfilter(g, w, 'replicate');
    case 'gmean'
        f = gmean(g, m, n);
    case 'hmean'
        f = harmean(g, m, n);
    case 'chmean'
        f = charmean(g, m, n, Q);
    case 'median'
        f = medfilt2(g, [m n], 'symmertic');
    case 'max'
        f = imdilate(g, ones(m, n));
    case 'min'
        f = imerode(g, ones(m, n));
    case 'midpoint'
        f1 = ordfilt2(g, 1, ones(m, n), 'symmetric');
        f2 = ordfilt2(g, m*n, ones(m, n), 'symmetric');
        f = imlincomb(0.5, f1, 0.5, f2);
    case 'atrimmed'
        f = alphatrim(g, m, n, d);
    otherwise
        error('Uknown filter type.');
end

%-----------------------------------------------------------
function f = gmean(g, m, n)
% Implements a geometric mean filter
[g, revertClass] = tofloat(g);
f = exp(imfilter(log(g), ones(m, n), 'replicate')).^(1 / m / n);
f = revertClass(f);

%-----------------------------------------------------------
function f = harmean(g, m, n)
% Implements a harmonic mean filter
[g, revertClass] = tofloat(g);
f = m * n / imfilter(1./(g + eps), ones(m, n), 'replicate');
f = revertClass(f);

%-----------------------------------------------------------
function f = charmean(g, m, n, q)
% Implements a contraharmonic mean filter
[g, revertClass] = tofloat(g);
f = imfilter(g.^(q+1), ones(m, n), 'replicate');
f = f ./ imfilter(g.^(q), ones(m, n), 'replicate');
f = revertClass(f);

%-----------------------------------------------------------
function f = alphatrim(g, m, n, d)
% Implements a alpha-trimmed mean filter.
if (d <= 0) || (d/2 ~= round(d/2))
    error('d must be a positive, even integer');
end
[g, revertClass] = tofloat(g);
f = imfilter(g, ones(m, n), 'symmetric');
for k = 1:d/2
    f = f - ordfilt2(g, k, ones(m, n), 'symmetric');
end
for k = (m*n - (d/2) + 1):m*n
    f = f - ordfilt2(g, k, ones(m, n), 'symmetric');
end
f = f / (m*n - d);
f = revertClass(f);

%-----------------------------------------------------------
function [m, n, Q, d] = processInputs(varargin)
m = 3;
n = 3;
Q = 1.5;
d = 2;
if nargin > 0
    m = varargin{1};
end
if nargin > 1
    n = varargin{2};
end
if nargin > 2
    Q = varargin{3};
    d = varargin{3};
end

Example: Using function spfilt

>> [M, N] = size(f);
>> R = imnoise2('salt & pepper', M, N, 0.1, 0);
>> gp = f;
>> gp(R == 0) = 0;
>> R = imnoise2('salt & pepper', M, N, 0, 0.1);
>> gs = f;
>> gs(R == 1) = 255;
>> fp = spfilt(gp, 'chmean', 3, 3, 1.5);
>> fs = spfilt(gs, 'chmean', 3, 3, -1.5);
>> fpmax = spfilt(gp, 'max', 3, 3);
>> fsmin = spfilt(gs, 'min', 3, 3);

结果如图所示：

3.2 Adaptive Spatial Filters

这里制定以下规则：
$$\begin{array}{rlrl}\text{Level}A:& \text{IF}{z}_{min}<z_{med}<z_{max}& & \text{go to Level}B\\ & \text{ELSE}& & \text{increases the window size}\\ & \text{IF}\text{window size}\le S_{max}& & \text{repeat Leval}A\\ & \text{ELSE}& & \text{output}{z}_{med}\\ \text{Leval}B:& \text{IF}{z}_{min}<z_{xy}<z_{max}& & \text{output}{z}_{xy}\\ & \text{ELSE}& & \text{output}{z}_{med}\end{array}$$
实现如下：

function f = adpmedian(g, Smax)
%ADPMEDIAN Perform adaptive median filtering
%   F = ADPMEDIAN(G, SMAX) performs adaptive median filtering of image G.
%   The median filter starts at size 3-by-3 and iterates up to size
%   SMAX-by-SMAX, SMAX must be an odd integer greater than 1.
% 
%   Rules:
%      
%        LevelA:  If zmin < zmed < zmax, go to LevelB
%                 Else increase the window size
%                 If window size <= Smax, repeat LevelA
%                 Else output Zmed
% 
%        LevelB:  If zmin < zxy < zmax, output zxy
%                 Else output zmed


% SMAX must be and odd, positive integer greater than 1.
if (Smax <= 1) || (Smax/2 == round(Smax/2)) || (Smax ~= round(Smax))
    error('SMAX must be and odd integer > 1.');
end

% Initial setup
f = g;
f(:) = 0;
alreadyProcessed = false(size(g));

% Begin filtering
for k = 3:2:Smax  % repeat LevelA
    zmin = ordfilt2(g, 1, ones(k, k), 'symmetric');
    zmax = ordfilt2(g, k*k, ones(k, k), 'symmetric');
    zmed = medfilt2(g, [k k], 'symmetric');
    
    processUsingLevelB = (zmed > zmin) & (zmax > zmed) & ...
        ~alreadyProcessed;
    zB = (g > zmin) & (zmax > g);
    outputZxy = processUsingLevelB & zB;
    outputZmed = processUsingLevelB & ~zB;
    f(outputZxy) = g(outputZxy);
    f(outputZmed) = g(outputZmed);
    
    alreadyProcessed = alreadyProcessed | processUsingLevelB;
    if all(alreadyProcessed(:))
        break;
    end
end

% Output zmed for any remaining unprocessed pixels. Note that this zmed was
% computed using a window of size Smax-by-Smax, which is the final value of
% k in the loop.
f(~alreadyProcessed) = zmed(~alreadyProcessed);

4 Periodic Noise Reduction Using Frequency Domain Filtering

具有Q notch pairs的Notchreject Filter的表达式为
$$H_{NR}(u,v)=\prod _{k=1}^Q{H}_k(u,v)H_{-k}(u,v)$$
其中，$H_k(u,v)$和$H_{-k}(u,v)$是分别以$(u_k,v_k)$和$(-u_k,-v_k)$为中心的HPF，这些中心是以频率矩形的中心(M/2, N/2)来确定的，所以，每个Filter的distance function有：
$$\begin{gathered} D_k(u,v)=[(u-M/2-u_k{)}^2+(v-N/2-v_k{)}^2{]}^{\frac{1}{2}} \\ {D}_{-k}(u,v)=[(u-M/2+u_k{)}^2+(v-N/2+v_k{)}^2{]}^{\frac{1}{2}} \end{gathered}$$
沿着频率轴开槽分量的陷波带阻滤波的some special cases，可以用之前自定义的function cnotch和function recnotch实现image restoration。

5 Modeling the Degradation Function

主要使用到的方法：生成点扩散函数PSF，利用复原算法对结果进行测试。在不知道任何PSF信息时，可以使用***blind convolution***获得。

在MATLAB中，使用function imfilter、function fspecial和noise functions等对PSF进行建模。

根据前面的section可以得知，Image Degradation Process的一个重要步骤是Imgae Blur。可以通过使用function fspecial实现：

PSF = fspecial('motion', len, theta);

其中，len默认9，theta默认0(水平)，表示object逆时针在对应theta角度方向上移动len个pixel。

之后，根据创造的PSF，使用imfilter生成一个退化后的图像。padding方法使用circular，以此减少边缘效应。

g = imfilter(f, PSF, 'circular');

最后直接进行与noise的叠加即可。

g = g + noise;

在测试阶段，为了保证测试的可比性，一般要选用same size的image和test pattern。这里是用function checkerboard生成test pattern

C = checkerboard(NP, M, N)

参数信息：

• NP：每个正方形一遍的像素数。默认10个pixels
• M、N：M是行数，N是列数。N省略时，默认为M。M、N都省略时，默认为8x8的正方形。

Example: Modeling a blurred, noisy image.

这里，我们生成一个逆时针在45度方向上行走7步的PSF，和一个均值为0、方差为0.001的高斯分布噪声。

>> f = checkerboard(8); % Image if of class double. A board image
>> PSF = fspecial('motion', 7, 45);
>> gb = imfilter(f, PSF, 'circular');
>> noise = imnoise2('Gaussian', size(f, 1), size(f, 2), 0, sqrt(0.001));
>> g = gb + noise;
>> subplot(2,2,1);imshow(pixeldup(f, 8), []);title('original image');
>> subplot(2,2,2);imshow(pixeldup(gb, 8), []);title('degraded image');
>> subplot(2,2,3);imshow(pixeldup(noise, 8), []);title('noise image');
>> subplot(2,2,4);imshow(pixeldup(g, 8), []);title('blurred image with degradation and noise ');

6 Direct Inverse Filter

Image Restoration的一个最简单的方法是忽略degradation process中的noise。因为非常容易得到如下表达式：
$$\hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}$$
然后通过获得$\hat{F}(u,v)$的IDFT进行估计。这个过程称为inverse filtering逆滤波。

如果考虑noise，表达式应为：
$$\hat{F}(u,v)=F(u,v)+\frac{N(u,v)}{H(u,v)}$$
但由于noise function是随机的，所以它的DFT是未知的。所以这种情况下很难求的精确的original image function。

注意：inverse filtering最典型的方法是形成第一个比例式子。并将frequency的范围尽可能限制接近频率原点。

7 Wiener Filter

Wiener Filter是一种Linear Filter，通过估计$\hat{f}$去最小化error function：
e 2 = E { ( f − f ^ ) 2 } {e^2=E\{(f-\hat{f})^2\}} e2=E{(f−f^​)2}
其中$E$是期望值算子，$f$是未退化的图像。频域中的表达式和需要满足的条件为：
$$\begin{array}{rl}\text{Expression: }& \hat{F}(u,v)=\left[\frac{1}{H(u,v)}\frac{|H(u,v{|}^2}{|H(u,v){|}^2+S_{\eta }(u,v)/S_f(u,v)}\right]G(u,v)\\ & \\ \text{Conditions: }& H(u,v)=\text{the degradation function}\\ & |H(u,v){|}^2=H^{\ast }(u,v)H(u,v)\\ & H^{\ast }(u,v)=\text{the complex conjugate of}H(u,v)\\ & S_{\eta }(u,v)=|N(u,v){|}^2=\text{the power spectrum of the noise}\\ & S_f(u,v)=|F(u,v){|}^2=\text{the power spectrum of the undegraded image}\end{array}$$
其中比例式$S_{\eta }(u,v/S_f(u,v))$也称为noise-to-signal power ratio信噪功率比。当noise的功率谱为0时，wiener filter也就变成了最常见的inverse filter逆滤波器。

Average noise power平均噪声功率和average image power平均图像功率：
$$\begin{array}{rl}\text{Average Noise Power(ANP): }& {\eta }_A=\frac{1}{MN}\sum _u\sum _v{S}_{\eta }(u,v)\\ \text{Average Image Power(AIP): }& f_A=\frac{1}{MN}\sum _u\sum _v{S}_f(u,v)\\ \text{Average N/I ratio(ANIR): }& R=\frac{{\eta }_A}{f_A}\end{array}$$
其中，上述三个表达式的值均为标量。且$R$通常被用做代替信噪功率比以此产生常量数组。

Wiener Filter在MATLAB中使用function deconvwnr实现，语法如下。

如果信噪功率比为0时，也就是inverse filter：

frest = deconvwnr(g, PSF)

如果信噪功率比已知，或是constant，或是array，都可以背deconvwnr接受，NSPR为标量。语法如下：

frest = deconvwnr(g, PSF, NSPR)

假设noise和undegraded image的自相关函数是已知的，即NACORR、FACORR，deconvwnr里使用$\eta$和$f$的自相关来代替这些function的spectrum of power：

frest = deconvwnr(g, PSF, NACORR, FACORR)

Ringing Effect：振铃效应，是由于在image restoration中选取了不适当的图像模型造成的，从而导致信息量(尤其是高频信息)的丢失。表现为输出图像的灰度剧烈变化处产生的震荡，就好像钟被敲击后产生的空气震荡。

这里引用function edgetaper，语法为：

J = edgetaper(I, PSF)

edgetaper利用PSF模糊处理了输入图像I的边缘。输出J为图像I和I的模糊版本的加权和。

Example Using function deconvwnr to restore a blurred, noisy image：

>> subplot(2,2,1);imshow(pixeldup(g, 8), []);title('blurred image');
>> frest1 = deconvwnr(g, PSF);
>> subplot(2,2,2);imshow(pixeldup(frest1, 8), []);title('Result of inverse filtering');
>> Sn = abs(fft2(noise)).^2; % noise power spectrum
>> nA = sum(Sn(:))/numel(noise); % noise average power
>> Sf = abs(fft2(f)).^2; % image power spectrum
>> nf = sum(Sf(:))/numel(f); % image average power
>> fA = sum(Sf(:))/numel(f); % image average power
>> R = nA/fA;
>> frest2 = deconvwnr(g, PSF, R);
>> subplot(2,2,3);imshow(pixeldup(frest2, 8), []);title('Result of Wiener filtering using a constant ratio');
>> NCORR = fftshift(real(ifft2(Sn)));
>> ICORR = fftshift(real(ifft2(Sf)));
>> frest3 = deconvwnr(g, PSF, NCORR, ICORR);
>> subplot(2,2,4);imshow(pixeldup(frest3, 8), []);title('Result of Wiener filtering using autocorrelation functions');

8 Constrained Least Squares(Regularized) Filter约束的最小二乘(规则化)滤波

回顾一些知识点：
KaTeX parse error: Expected 'EOF', got '&' at position 18: …text{二维离散卷积公式：}&̲h(x,y)\odot f(x…
为了方便函数推导，所以采用向量的方式表示。接下来，定义一下准则函数C：
KaTeX parse error: Expected 'EOF', got '&' at position 15: \text{准则函数：} &̲ C=\sum_{m=0}^{…
上述formulation中，只有$\gamma$和$||\eta |{|}^2$是未知的。但当$||\eta |{|}^2$已知时，可以通过迭代求出$\gamma$。

在MATLAB中，约束的最小二乘法的实现使用function deconverg：

frest = deconvreg(g, PSF, NOISEPOWER, RANGE)

其中，NOISEPOWER是噪声功率，与$||\eta |{|}^2$成比例，最好的估计范围是$MN[{\sigma }_{\eta }^2+{mean}_{\eta }^2]$；RANCE时值的范围，默认范围是$[10^{-9},10^9]$

Example Using function deconvreg to restore a blurred noisy image：

>> frest1 = deconvreg(g, PSF, 4);
>> frest2 = deconvreg(g, PSF, 0.4, [1e-7, 1e7]);
>> subplot(1,2,1);imshow(pixeldup(frest1, 8), []);title('NOISEPOWER=4');
>> subplot(1,2,2);imshow(pixeldup(frest2, 8), []);title('NOISEPOWER=0.4, RANGE=[1e-7,1e7]');

9 Iterative Nonlinear Restoration Using the Lucy-Richardson Algorithm

😄开始一个非线性的图像还原算法！！！

先说一下nonlinear的优点和缺点：

• 优点：理论基础容易建立，执行简单；复原滤波一旦确定，相应的solution就会通过filter应用到。

• 缺点：其行为常常不可预见；需要重要的计算资源。

有关L-R算法，起源于maximum-likelihood formulation(极大似然公式)，图像是由Poisson statistics(泊松分布)建模。当迭代收敛时，可以得到方程式：
$${\hat{f}}_{k+1}(x,y)={\hat{f}}_k(x,y)[h(-x,-y)\odot \frac{g(x,y)}{h(x,y)\odot {\hat{f}}_k(x,y)}]$$
我们对L-R算法进行一个简单的推导：
$$\begin{array}{rl}\text{连续贝叶斯定理：}& P(x|y)=\frac{P(y|x)P(x)}{\int P(y|x)P(x)dx}\\ \text{离散贝叶斯定理：}& P(x|y)=\frac{P(y|x)P(x)}{\sum P(y|x)P(x)}\\ \text{模糊图像最大概率：}& P(f|g)=\frac{P(g|f)P(f)}{P(g)}\\ \text{泊松分布模型公式：}& P(x)=\frac{u^x{e}^{-u}}{x!}\\ \text{泊松分布条件概率：}& p(g|f)=\prod _{(x,y)}\frac{[h(x,y)\odot f(x,y){]}^{g(x,y)}{e}^{-(h(x,y)\odot f(x,y))}}{g(x,y)!}\\ \text{取对数化简上述方程：}& \ln p(g|f)=\sum _{(x,y)}[g(x,y)\ln [h(x,y)\odot f(x,y)]-h(x,y)\odot f(x,y)-\ln g(x,y)!]\\ \text{对}f(x,y)\text{求导：}& \frac{\partial \ln p(g|f)}{\partial f(x,y)}=0\text{ 解得 }[\frac{g(x,y)}{h(x,y)\odot f(x,y)}\odot {h(x,y)}^T]=1\\ \text{化简方程解：}& f(x,y)=[\frac{g(x,y)}{h(x,y)\odot f(x,y)}\odot {h(x,y)}^T]f(x,y)\\ \text{加入迭代概念得：}& {\hat{f}}_{k+1}(x,y)={\hat{f}}_k(x,y)[h(-x,-y)\odot \frac{g(x,y)}{h(x,y)\odot {\hat{f}}_k(x,y)}]\end{array}$$
在MATLAB中使用function deconvlucy：

f = deconvlucy(g, PSF, NUMIT, DAMPAR, WEIGHT)

其中，NUMIT是迭代次数，默认为10；DAMPAR是结果图像与原始图像g之间的偏离阈值，当pixel偏离original value在阈值范围内时，结束迭代，默认0；WEIGHT是数组，即每个像素的权重，反映了像素的质量，当PSF的大小为NxN时，WEIGHT用到的零边界的宽度就是ceil(N/2)，默认为输入图像g的大小的单位数组。

Example Using function deconvlucy to restore a blurred noisy image：

>> g = checkerboard(8);
>> subplot(2,3,1);imshow(pixeldup(g, 8));title('Original image');
>> PSF = fspecial('gaussian', 7, 10); % 7x7的高斯点扩散函数
>> SD = 0.01;
>> g = imnoise(imfilter(g, PSF), 'gaussian', 0, SD^2); % 均值为0，标准差为0.01的高斯噪声
>> subplot(2,3,4);imshow(pixeldup(g, 8));title('Image blurred and corrupted by Gaussian');
>> DAMPAR = 10*SD;
>> LIM = ceil(size(PSF, 1)/2); % WEIGHT中边界为0的宽度
>> WEIGHT = zeros(size(g)); % 生成WEIGHT数组
>> WEIGHT(LIM + 1:end - LIM, LIM + 1:end - LIM) = 1; % 给图像非边界的点赋值，边界点为0可以剔除来自图像边界的像素点
>> NUMIT = 5; % 迭代次数
>> f5 = deconvlucy(g, PSF, NUMIT, DAMPAR, WEIGHT);
>> subplot(2,3,2);imshow(pixeldup(f5, 8));title('Result from L-R:iteration=5');
>> f5 = deconvlucy(g, PSF, 10, DAMPAR, WEIGHT);
>> subplot(2,3,5);imshow(pixeldup(f5, 8));title('Result from L-R:iteration=10');
>> f5 = deconvlucy(g, PSF, 20, DAMPAR, WEIGHT);
>> subplot(2,3,3);imshow(pixeldup(f5, 8));title('Result from L-R:iteration=20');
>> f5 = deconvlucy(g, PSF, 100, DAMPAR, WEIGHT);
>> subplot(2,3,6);imshow(pixeldup(f5, 8));title('Result from L-R:iteration=100');

10 Blind Deconvolution

简单粗暴：不考虑PSF点扩散函数的Deconvolution称为Blind Deconvolution。也就是说，Blind Deconvolution要估计PSF。

基本方法是Maximum-Likelihood Estimation，是对随机噪声污染的量进行估计时采用的最佳策略。

在MATLAB中，使用function deconvblind实现：

[fr, PSF] = deconvblind(g, INITPSF)

其中，INITPSF是点扩散函数的初始估计，PSF是计算得到的估计值，fr是利用估计的PSF复原的图像。

如果加上L-R算法和function deconvlucy中需要考虑的变量的话，function deconvblind的使用如下：

[fr, PSF] = deconvblnd(g, INITPSF, NUMIT, DAMPAR, WEIGHT)

Example Using deconvblind to estimate a PSF：

>> PSF = fspecial('gaussian',7,10);
>> subplot(221);imshow(pixeldup(PSF, 73), []);title('Original PSF');
>> SD = 0.01;
>> g = imnoise(imfilter(g, PSF), 'gaussian', 0, SD^2);
>> INITPSF = ones(size(PSF));
>> NUMIT = 5;
>> [g5, PSF5] = deconvblind(g, INITPSF, NUMIT, DAMPAR, WEIGHT);
>> subplot(222);imshow(pixeldup(PSF5, 73), []);title('iteration = 5');
>> [g5, PSF5] = deconvblind(g, INITPSF, 10, DAMPAR, WEIGHT);
>> subplot(223);imshow(pixeldup(PSF5, 73), []);title('iteration = 10');
>> [g5, PSF5] = deconvblind(g, INITPSF, 20, DAMPAR, WEIGHT);
>> subplot(224);imshow(pixeldup(PSF5, 73), []);title('iteration = 20');

11 Image Reconstruction from Projections来自投影的图像重建

11.1 Parallel-Beam Projections and the Radon Transform

我们知道，一条直线在笛卡尔坐标系下的表达式为$y=ax+b$，其法线表达式为$x\cos \theta +y\sin \theta =\rho$。

Parallel-Beam Projections可以通过一组直线建模，在投影的刨面上的任意一点$({\rho }_j,{\theta }_k)$都可以沿着$x\cos {\theta }_k+y\sin {\theta }_k=\rho$的射线求和得出，即求线积分：
$$g({\rho }_j,{\theta }_k)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)\delta (x\cos {\theta }_k+y\sin {\theta }_k-\rho )dxdy$$
这里$\delta$指impulse function脉冲函数。根据脉冲函数的特性$\delta (t)=\{\begin{array}{ll}\infty & t=0\\ 0& t\ne 0\end{array}$和${\int }_{-\infty }^{\infty }\delta (t)dt=1$，我们可以得到线积分等式右边$\delta$函数的传参为0，即$x\cos {\theta }_k+y\sin {\theta }_k-\rho =0$时，等式为0。

如果我们考虑到所有的法线长度和偏向角度数的话，上述等式可以化简为：
$$g(\rho ,\theta )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)\delta (x\cos \theta +y\sin \theta -\rho )dxdy$$
这个表达式给出沿xy平面内任意一条直线f(x,y)的投影(线积分)，也就是radon transform雷登变换。其离散形式写为：
$$g(\rho ,\theta )=\sum _{x=0}^M\sum _{y=0}^{N}f(x,y)\delta (x\cos \theta +y\sin \theta -\rho )$$
Back-projection：这里简单说一下反投影的步骤

1. 定义：$g(\rho ,{\theta }_k)$全部投影、$g({\rho }_j,{\theta }_k)$投影上的单个点、${\theta }_k$固定
2. 对于固定角度的单个点，反投影图像表达式：$f_{{\theta }_k}(x,y)=g(\rho ,{\theta }_k)=g(x\cos {\theta }_k+y\sin {\theta }_k,{\theta }_k)$
3. 对于所有角度的单个点，反投影图像表达式：$f_{\theta }(x,y)=g(\rho ,\theta )=g(x\cos \theta +y\sin \theta ,\theta )$
4. 通过对所有点反投影图像进行积分得到最终图像：$f(x,y)={\int }_0^{\pi }{f}_{\theta }(x,y)d\theta$。由于$[0,\pi ]$和$[\pi ,2\pi ]$间隔投影相同，所以积分范围只在半周上执行。
5. 对于离散点，最终图像的表达式为：${f(x,y)=\sum_{\theta=0}^\pi f_\theta(x,y)}\$

11.2 The Fourier Slice Theorem and Filtered Backprojections

关于$\rho$的$g(\rho ,\theta )$的一维傅里叶转换方程为：${G(\omega,\theta)=\int_{-\infty}^\infty g(\rho,\theta)e^{-j2\pi\omega\rho}d\rho}\$，其中${\omega}$是频率变量，${\theta}$是固定值。

二维变换的切片，即***傅里叶切片定理***为：$G(\omega ,\theta )=[F(u,v){]}_{u=\omega \cos \theta ,v=\omega \sin \theta }=F(\omega \cos \theta ,\omega \sin \theta )$

那么，逆傅里叶变换的公式为：${f(x,y)=\int_{-\infty}^{\infty\int_{-\infty}}\infty F(u,v)e^{j2\pi(ux+vy)}dudv}\$

我们结合投影平面上2-D Fourier Transform Slice中u和v的定义我们可以得到极坐标下2-D Fourier Transform Slice Expression为：
$$\begin{gathered} u=\omega \cos \theta \\ v=\omega \sin \theta \\ dudv=\omega d\omega d\theta \\ f(x,y)={\int }_0^{2\pi }{\int }_0^{\infty }F(\omega \cos \theta ,\omega \sin \theta )e^{j2\pi \omega (x\cos \theta ,y\sin \theta )}\omega d\omega d\theta \\ f(x,y)={\int }_0^{2\pi }{\int }_0^{\infty }G(\omega ,\theta )e^{j2\pi \omega (x\cos \theta ,y\sin \theta )}\omega d\omega d\theta \end{gathered}$$

将上述的积分分成两个表达式：一个是针对$\theta$的表达式，范围是$[0,\pi ]$；另一个用于$G(\omega ,\theta +\pi )=G(-\omega ,\theta )$的情况，范围是$[\pi ,2\pi ]$：
$$f(x,y)={\int }_0^{\pi }{\int }_{-\infty }^{\infty }|\omega |G(\omega ,\theta )e^{j2\pi \omega (x\cos \theta +y\sin \theta )}d\omega d\theta$$
对于变量$\omega$，$x\cos \theta +y\sin \theta$是一个常量，记做$\rho$。那么上面的积分就可以写为：
$$f(x,y)={\int }_0^{\pi }[{\int }_{-\infty }^{\infty }|\omega |G(\omega ,\theta )e^{j2\pi \omega \rho }d\omega {]}_{\rho =x\cos \theta +y\sin \theta }d\theta$$
实际上，这里的$|\omega |$是一个滤波函数。那么，完整的back-projection过程为：

1. 计算每一个投影(点)的一维傅里叶转换
2. 利用滤波函数$|\omega |$乘以每一个傅里叶转换
3. 得到步骤2中的滤波结果后做逆傅里叶转换操作
4. 对步骤3中逆傅里叶转换的结果求积分，得到最终图像$f(x,y)$

因为在步骤2中使用了滤波函数，所以image reconstruction也称为filtered back-projection。

11.3 Filter Implement

在filter back-projection中，滤波器$|\omega |$形状呈斜坡状，在连续状态下不可积。

但是，作为频率函数，很容易受到噪声影响，在空间域上还会有振铃效应产生。一种解决方法是用窗口函数乘以斜坡滤波器，使滤波器的拖尾呈渐变状态，在high-frequency处减少幅度。

我们常见的窗口函数有sine、cosine、Hamming、Hann。

sine窗口的传递函数为：
$$H_s(\omega )=\frac{sin(\pi \omega /2\Delta \omega K)}{(\pi \omega /2\Delta \omega K)}\omega =0,\pm \Delta \omega ,\pm 2\Delta \omega ,\dots ,\pm K\Delta \omega$$
其中K是频率间隔数(点数减1)。

同理，cosine窗口的传递函数为：
$$H_c(\omega )=\cos \frac{\pi \omega }{2\Delta \omega K}$$
Hamming和Hann窗口的传递函数为：

$$H(\omega )=c+(c-1)\cos \frac{2\pi \omega }{\Delta \omega K}$$
当c=0.54时，为Hamming窗口(点的末端有较少的偏移)；当c=0.5时，为Hann窗口(点的末端是0)。

11.4 Function radon

我们用MATLAB中的function radon为给定的二维矩阵生成一组平行射线投影，基本语法为：

R = radon(I, theta)

其中，I是一个二维矩阵，theta是角度值的一维矩阵，输出R中的投影矩阵数量等于角度值矩阵theta的数量。产生的投影要足够长，以便在射束旋转时跨越观察的宽度。

function radon的一般形式为：

[R, xp] = radon(I, theta)

其中，xp包含了前面我们提到的$\rho$值。

在CT算法中，我们还要用到另一个function phantom

P = phantom(def, n)

其中，def指定头部模型，n是行数和列数，默认为256。关于def，可用值有：

• Shepp-Logan：主要用于测试图像，但对比度很低。
• Modified Shepp-Logan：上面模型的变量，对比度有所提高。

Example Using function radon：

>> g1 = zeros(600, 600);
>> g1(100:500, 250:350) = 1;
>> g2 = phantom('Modified Shepp-Logan', 600);
>> subplot(221);imshow(g1);xlabel('(a)');title('original image');
>> subplot(223);imshow(g2);xlabel('(c)');title('Modified Shepp-Logan Model');
>> theta = 0:0.5:179.5;
>> [R1, xp1] = radon(g1, theta);
>> [R2, xp2] = radon(g2, theta);
>> R1 = flipud(R1');
>> R2 = flipud(R2');
>> subplot(222);imshow(R1, [], 'XData', xp1([1 end]), 'YData', [179.5 0]);title('Radon transforms of image-a');
>> axis xy
>> axis on
>> xlabel('\rho'), ylabel('\theta')
>> subplot(224);imshow(R2, [], 'XData', xp2([1 end]), 'YData', [179.5 0]);title('Radon transforms of image-c');
>> axis xy
>> axis on
>> xlabel('\rho'), ylabel('\theta')

11.5 Function iradon

这个函数从给定的取自不同角度的一组投影来重建一幅图像(切片)，主要用于计算逆雷登变换。基本语法为：

I = iradon(R, theta, interp, filter, frequency_scaling, output_size)

参数信息如下：

• R：反投影数据，列是从左到右以角度渐增的函数组织的一维反投影
• theta：投影捕捉的角度。可以是标量或向量。若为向量，角度间距是相等的；若为标量，就认为投影是在$theta=m\ast D\_theta$角度下获得的，其中$m=0,1,2,\dots ,size(R,2)-1$，标量默认大小为180/size(R, 2)
• iterp：重建图像的插值方法。主要有四种：nearest、linear、cubic、spline
• filter：滤波反投影计算中使用的滤波器。主要有六种：
  • Ram-Lak：一种斜坡滤波器。且为默认值。
  • Shepp-Logan：用sine window乘以Ram-Lak滤波器产生的filter
  • Cosine：用Cosine window乘以Ram-Lak滤波器产生的filter
  • Hamming：用Hamming window乘以Ram-Lak滤波器产生的filter
  • Hann：用Hann window乘以Ram-Lak滤波器产生的filter
  • None：不进行滤波
• frequency_scaling：标量，范围是(0,1)，默认为1。通过改变频率轴的比例来选择滤波器。在规一化频率下，所有大于frequency_scaling的频率置于0。
• output_size：确定重建图像的大小。默认由投影长度决定：2*floor(size(R,1)/(2*sqrt(2)))

Example Using function iradon：

>> subplot(241);imshow(g1);xlabel('(a)');title('original image');
>> subplot(245);imshow(g2);xlabel('(b)');title('Modified Shepp-Logan Model');
>> theta = 0:0.5:179.5;
>> R1 = radon(g1, theta);
>> R2 = radon(g2, theta);
>> f1 = iradon(R1, theta, 'none');
>> f2 = iradon(R2, theta, 'none');
>> subplot(242);imshow(f1, []);xlabel('(c)');title('Backprojection without filtering from a');
>> subplot(246);imshow(f2, []);xlabel('(d)');title('Backprojection without filtering from b');
>> f1_ram = iradon(R1, theta);
>> f2_ram = iradon(R2, theta);
>> subplot(243);imshow(f1_ram, []);xlabel('(e)');title('Backprojection with Ram-Lak filtering');
>> subplot(247);imshow(f2_ram, []);xlabel('(f)');title('Backprojection with Ram-Lak filtering');
>> f1_hamm = iradon(R1, theta, 'Hamming');
>> f2_hamm = iradon(R2, theta, 'Hamming');
>> subplot(244);imshow(f1_hamm, []);xlabel('(g)');title('Backprojection with Hamming filtering');
>> subplot(248);imshow(f2_hamm, []);xlabel('(h)');title('Backprojection with Hamming filtering');

11.6 Working with Fan-Beam Data

在MATLAB中，主要使用function fanbeam

B = fanbeam(g, D, param1, val1, param2, val2,...)

其中，g是包含被投射物体的图像；D是从扇形射束的顶点到旋转中心的距离，如果旋转中心是图像的中心，则满足$D=K\times {\text{sqrt(size(g,1)}}^2{\text{size(g,2)}}^2)/2$，K是个大于1的常数，一般在1.5到2之间取值。

这里简单说一下fanbeam使用的参数和值：

• FanRotationIncrement：以度进行度量的扇形射束投影的旋转角度增量，值为正的标量，默认为1。
• FanSensorGeometry：规定如何等间隔安排传感器的文本串，默认为arc，也可以为line。
• FanSensorSpacing：规定扇形射束中传感器间隔的、正的实标量，若为arc，说明值以度为角度间距；若为line，说明以线为间隔。默认为1。

Example Using function fanbeam：

>> g1 = zeros(600, 600);
>> g1(100:500, 250:350) = 1;
>> g2 = phantom('Modified Shepp-Logan', 600);
>> D = 1.5*hypot(size(g1,1), size(g1,2))/2;
>> B1_line = fanbeam(g1, D, 'FanSensorGeometry', 'line', 'FanSensorSpacing', 1, 'FanRotationIncrement', 0.5);
>> B1_line = flipud(B1_line');
>> B2_line = fanbeam(g2, D, 'FanSensorGeometry', 'line', 'FanSensorSpacing', 1, 'FanRotationIncrement', 0.5);
>> B2_line = flipud(B2_line');
>> B1_arc = fanbeam(g1, D, 'FanSensorGeometry', 'arc', 'FanSensorSpacing', .1, 'FanRotationIncrement', 0.5);
>> B2_arc = fanbeam(g2, D, 'FanSensorGeometry', 'arc', 'FanSensorSpacing', .1, 'FanRotationIncrement', 0.5);
>> subplot(221);imshow(B1_line, []);
>> subplot(222);imshow(B2_line, []);
>> subplot(223);imshow(B1_arc, []);
>> subplot(224);imshow(B2_arc, []);

这里也说一下inverse fanbeam，也就是function ifanbeam：

I = ifanbeam(B, D, ..., param1, val1, param2, val2, ...)

参数信息：

• FanCoverage：指定射束的旋转范围。默认为circle，即旋转一周；minimal指的是物体所需的最小范围

Example Working with function ifanbeam：

>> g = phantom('Modified Shepp-Logan', 600);
>> D = 1.5*hypot(size(g1,1), size(g1,2))/2;
>> B1 = fanbeam(g, D);
>> f1 = ifanbeam(B1, D);
>> subplot(131);imshow(f1,[]);
>> B2 = fanbeam(g, D, 'FanRotationIncrement', 0.5, 'FanSensorSpacing', 0.5);
>> f2 = ifanbeam(B2, D, 'FanRotationIncrement', 0.5, 'FanSensorSpacing', 0.5, 'Filter', 'Hamming');
>> subplot(132);imshow(f2,[]);
>> B3 = fanbeam(g, D, 'FanRotationIncrement', 0.5, 'FanSensorSpacing', 0.05);
>> f3 = ifanbeam(B3, D, 'FanRotationIncrement', 0.5, 'FanSensorSpacing', 0.05, 'Filter', 'Hamming');
>> subplot(133);imshow(f3,[]);

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-uEmWzRyq-1647138748562)(/Users/wanghe/Library/Application Support/typora-user-images/image-20211119234834837.png)]

这里补充一个扇形和平行投影互相转换的函数fan2para和para2fan：

P = fan2para(F, D, parm1, val1, parm2, val2, ...)

参数信息如下：

• ParallelCoverage：circle表示平行覆盖360度，halfcyle为默认值，表示平行覆盖180度。
• ParallelRotationIncrement：平行射束的角度增量。默认与扇形射束的角度增量一致。
• ParallelSensorSpacing：平行射束传感器间隔，以像素计量。默认是均匀的。

F = fan2para(P, D, parm1, val1, parm2, val2, ...)

参数信息如下：

• FanRotationIncrement：以度进行度量的扇形射束投影的旋转角度增量，值为正的标量，默认为1。若FanCoverage为circly，则FanRotationIncrement为360的倍数。默认与平行射束的角度相同。
• FanSensorSpacing：如果为arc或line，与fanbeam中的一致。默认为平行射束ParallelSensorSpacing指定的最小值。若FanSensorGeometry为arc时，FanSensorSpacing是$\frac{180}{\pi }\arcsin (\frac{\text{PSPACE}}{D})$，其中PSPACE是ParallelSensorSpacing的值；若FanSensorGeometry为line时，FanSensorSpacing是$D\arcsin (\frac{\text{PSPACE}}{D})$。