1 数据

  • 数据属性特征
    • 数据值类别
      • Nominal (labels)
      • Ordinal (ordered)
      • Interval (location of zero arbitrary)
      • Ratio (zero fixed)
      • 注意区分区间型和比值型，后者可以定义乘除法
    • 属性间的距离
      • 比值型下的各种定义
  • 数据预处理
    • 缺失值
    • 噪声值：回归分析、离群点分析
    • 可视化数据清洗
    • 数据整合：多数据源
  • 数据存储
    • 文件形式：CSV、HTML、XML
    • 数据库
    • 数据仓库 data warehouse

2.0 Image 图像

  • 从图像到数字图像
    • 天文、医学……
    • 采样 sampling
    • 量化 quantization
  • 分辨率
    • Radiometric resolution 量化的大小
    • Geometric resolution 通常意义下的分辨率
    • Image resolution 图片尺寸
  • 存储形式
    • 数字图像表达：\(f(X)=Y\)
    • Pixel
    • Coordinates
    • Gray level
  • 亮度和对比度
    • 视觉：brightness adaptation level
    • Weber ratio
    • Contrast
  • 坐标系/医学影像坐标系
    • 医学影像中不同的坐标系
    • 图像坐标-世界坐标的变换
    • \([w_x,w_y,w_z,1]'=M_{I2W}\times [i,j,k,1]'\)
    • 右侧是图像坐标系中的坐标，注意 M 是一个 \(3\times 4\) 矩阵，另加的一个维度用于平移。

2.1 Pixel operation 图像像素操作

  • Matrix persentation 以矩阵的形式对于图像进行表示
  • Neighbors of a pixel
  • Distance between pixels
  • Elementwise ops
  • Pixel ops
    • Single pixel op 如下面的灰度值变换
    • Neighbothood op 基于邻域的操作，很多算法的思想都有设计，如局部的直方图算法、插值等
    • Geometric spatial transformation 空间变换，在 9.0 中讲到
    • Image intensity as r.v.

3.1 Intensity transform 灰度值变换

  • Contrast enhancement

  • Negative 灰度值取反

  • Logarithmic

  • Power-law transformation 通过控制系数 γ 来调节凸显的部分，例如 γ 越大，灰度值高的区域对比度越高

  • Contrast stretching 自定义变换曲线 \(T(r)\)

  • Intensity-Level Slicing 凸显灰度值在一定区间内的部分，其余部分置零或是原样输出

3.2 Histogram-based transform

  • Histogram/nD joint histogram

Histogram equalization

  • 主要的想法是：灰度值集中的某个区域内，导致那部分的图像没有区分度，因此考虑做一个转换使得 Histogram 均匀分布；自然想到了👇的累积分布的转化（叫什么名字来着，统计好像学过）：

  • 连续情况是很理想的，但在离散状态下，若过多的像素点灰度值为同一个值，效果可能还是不够好

Histogram matching HM

  • 均衡化是对于全局进行的操作，然而有时我们仅仅需要突出某一区间内的灰度值细节，这就是 HM 的思想，将直方图再变换为所需的分布形式

  • 其中 T 是把原直方图均衡化，G 是把目标直方图均衡化；这里借由中间的均匀分布进行二次转换

Local Histogram processing

  • 一些场景下对于整张图片来进行 HM 是不合适的

  • 局部直方图均衡化：对于每一个关注的像素点 \(p_i\) ，根据其附近的区域 \(RoI_i\) 的直方图，计算转换 \(T_i\) ，对于点 \(p_i\) 进行更新；
  • 显然，我们对于每个点直接计算其 RoI 的直方图分布是低效的；也即需要 Efficient update of local histogram，比如我们把点向右移动一个像素点，可以发现其 RoI 的的变化仅仅是左右的两列，我们利用这些变化的像素点对于 Histogram 进行更新。

4.0 Basic Thresholding

  • 之间的直方图变换是为了增强对比度，而这里的 Threshold 可以看成是一种特殊的变换：二值化
  • 我们设定一个区间，在这个区间内的像素点灰度值变换为 1，其余部分变换为 0
  • 下面介绍了一些二值化算法（fix threshold），即仅确定通过在数轴上设立一个点分成两部分

Basic Global Thresholding (BGT)

  • A simple iterative threshold method, referred to as the Basic Global Thresholding (BGT) algorithm

  • 其实还挺迷的，想法是，根据这个 threshold 分割出来的直方图，左右两部分的两个均值的平均值正好等于我们设定的 threshold，迭代计算即可；不知道其背后有什么思想

OTSU Algorithm 大津算法

  • 这个的想法比较简单：想要使得两类像素点之间的类间方差尽可能大

  • 定义了几个量 ω 和 μ 分别是灰度值在一定区间内的频率和均值；m 是直接用 \(p(i)\) 计算的结果，没有进行归一化，因此有 \(m=\omega\mu\) 的关系（下面的符号上进行了简化）
  • 我们的目标是组间方差 \(\delta_{bet}^2\) ，好像和机器学习里的定义是一致的？以某一个的频率为权重对 \((\mu_{i}-\mu)^2\) 进行求和；
  • 计算的时候，我们希望尽量减少需要计算的量；这里第一行把 \(\mu\) 的形式带进去，上下同乘；第二三行进行变换，把 front 类的变量消去；可以看到最终的结果仅包含了 \(\mu, \omega_0, m_0\) 三个需要计算的量.

Entropy Method

  • 最大化划分出来的两部分灰度值分布的熵之和

  • 这里用 b,w 的符号有点混淆，反正就是前景和背景部分各自的熵；
  • 其实和大津算法很相似，后者是以组内/组间方差方差为准则，希望最大化组间方差，即最小化组内方差；而这里的熵算法则是最大化了两个分组内部的熵，同样也是希望是的组内的变异尽可能小。

Adaptive Thresholding

  • also called local (or regional) thresholding

  • 显然，全局的 threshold 会出现之前 HE 一样的问题，考虑局部处理

  • 还是一个简单的快速更新均值的公式（一维形式）；另外需要对于边缘进行细致的处理；
  • 可以采用简单的和 RoI 均值的比较设定阈值，公式的常数 C 可能是基于经验的；也可以用如大津算法等自动设定阈值；

5.0 Interpolation 插值

  • Four concepts
    • Single-pixel operation
    • Neighborhood ops
    • Geometric spatial transformation
    • Image intensity as random variables
    • 单像素的操作如 2.1 节的灰度值变换；邻域操作如平均模糊化，以及这里的图像插值；
  • 采样和量化
    • Sampling：digitizing the coordinate values
    • Quantization: digitizing the amplitude values
    • 如前所述采样和量化有其分辨率
  • Interpolation is the process of using known data to estimate values ar unknown locations
    • Image interpolation is an image resampling methods for image resizing (shrinking and zooming)
      • up-sampling and down-sampling
  • Gray-level interpolation methods for images

    • Nearest neighbot interpolation 找最近的那个点即可

    • Bilinear interpolation 线性插值，例如，在二维情况下，插值结果是四周各个点的平均，权重的其对角的矩形的面积

  • Cubic polynomial interpolation 上面用了周围的四个点，这里用四周的 16 个点
  • 空间标量场插值
    • 图像描述可能不止矩阵等形式，对于一个空间标量场来说，其插值的思想和在正常的笛卡尔系下是一样的，无非是利用连续性的假设，利用周围点对于目标位置加权求和，权重随距离增加而衰减

  • 当然，没有在上面的线性插值那么标准的算法了，例如这里的三角面片来说的，这里定义了两个方向，用了这样的形式定义了权重；注意到这里应该是和矩阵方式的存储一样各条边的长度为 1（尽管现实上不是等长的）；注意到，不同的两根轴的选取所得的结果有所不同；
  • 当然，例如像 cell 为六边形标量场，可以直接以目标点到各定点之间的距离，其倒数的 e 次方为权重进行衰减，考虑到发散情况加上一个平滑项 \(({1\over r_i}+\epsilon)^e\) 。

5.1 Filter

  • Spatial filter 空间滤波

    • Filter: passing, modifying or rejecting specific frequency
    • Contrast & Enhancement
    • Spatial filtering: Linear or Nonlinear

  • Convolution & Spatial correlation

  • 有下面的一些性质

Smoothing (noise reduction)

  • 平滑操作的背景是为了去噪，例如 X-ray image of a circuit board 中的 salt-and pepper noise 椒盐噪声

  • 最简单的如 Averaging filter，选取卷积核为取平均的大小为 3 的矩阵；

  • 再如 Gaussian Filter 基于二维高斯分布，取不同方差的高斯核，kernel 的大小也会不同； \[ w(s,t)=G(s,t)=Ke^{-{s^2+t^2\over 2\sigma^2}} \] 一般可取三倍方差范围内的点，例如取 \(\sigma=7\)，则 kernel size=21

  • 再如 Order-Statistics filter，如取 Median，这是 nonliner 的；在图像原本没有噪声的情况下，用前两者会有模糊化的效果，而这个相对来说好一些。

Edge detectation

  • 不同类型的边缘，前面的两个是快速和缓慢变化的边缘，最后一个可以理解成一条线/缝；它们的一些一阶导数和二阶导数呈现右边的形态

  • Image gradient：因此考虑借助图像梯度的信息进行边缘检测；在这样的离散状态下，可以通过不同的方式（只要合理）来定义梯度，都可以描述成一个 kernel

  • 这里用单独一个方向的导数矩阵变换可以得到一个方向上的边缘（实现的时候可以不用 kernel 的形式直接算出来即可）；分别计算两个方向的导数，汇总成梯度值图像可以实现多方向上的边缘检

Sharpening 锐化

  • （上面的边缘检测似乎只是用到了一阶梯度信息？）这里的锐化是想要凸显出图像的细节特征。一个典型的算法是 Laplacian filter 拉普拉斯滤波，用到了二阶梯度的信息，对于「边缘检测」一节中的（一维）边缘图像全部二阶求导；例如对于上一节中的 Step edge，原本直接只要一个单峰的一阶梯度图像，经过求导后会变成有一个正、紧接着一个负的两个小峰；若我们将其取反叠加到原图像上，那么这个边缘会变得更加「明显」。对于其他的情况可以有类似的结论。

  • 这是一个例子，可以看到 Laplacian 变换之后提取到了边缘信息，叠加之后月球的细节部分更为清晰。

  • 这里用上面的一种梯度公式推导出了 Laplacian kernel，形式上比较好理解；用不同的导数定义可以得到多种形式的 kernel。

  • Unsharp masking and highboost filtering：除了上面的借助二阶梯度的算法，我们可以从另外的角度来进行锐化——之前讲到的 smoothing 可以消去图像中尖锐的信息，现在我们取两者之差即可得到「残差图」，包含了原本图像中尖锐部分的信息；将其叠加在源图像上也可以达到锐化效果，称为 highboost 算法

  • 总结

    Laplacian: \(f-w\nabla^2f\)

    Highboost: \(f+k(f-\bar f)\)

    另外，如果我们采用之前的定义方式，那么 \[ \nabla^2f=f(-1,0)+f(1,0)+f(0,-1)+f(0,1)+f(0,0)-5f(0,0)=-5(f-\bar f) \] 于是，两者在一定条件下是等价的。

7.0 频域变换概念

关于傅立叶变换的部分，之前课上的时候查阅了一些材料，整理了单独的一篇笔记，见 https://www.cnblogs.com/easonshi/p/12676491.html 。因此在下面简单罗列不加以说明。

  • 傅立叶级数 Fourier Series

  • 傅立叶变换 Fourier Transform

  • 傅立叶逆变换

    • 统一利用欧拉公式写成指数形式，注意到傅立叶级数的本质是将一个函数分解到一个函数空间上，其中的系数是对应的表示向量；不断增加周期，原本离散的系数/向量就变成了一个连续的函数，即为「傅立叶变换」，将原本的时空域函数 \(f(t)\) 变换到了频率域上的表示 \(F(u)\)，而原本的级数公式就变为傅立叶逆变换的形式，两者组成一个傅立叶变换对。

傅立叶变换对

  • 频域和时空域（图像域）

  • 正弦、余弦项的频率由 u 决定，称傅里叶变换域为频率域；

  • 对于图像，t 成为时空域

  • 频谱 frequency spectrum: \(|F(u)|\)

  • 相位 phase angle: \(\theta(u)=arctan[I(u)/R(u)]\)

  • 上面用的是指数形式，虽然简洁但变换后每一个点的函数值都是复数；我们还是回到经典的基函数为常值函数、正弦、余弦的形式，并且可以看到 \(c_n\) 的实部和虚部 \(a_n, b_n\) 分别对应了周期相同的正余弦函数上的系数；在此基础上我们定义了频谱和相位。

    下图展示了两张图像的频谱图和相位图

常用的变换

  • 首先介绍了 冲击函数 impulse \(\delta(t)\)，即仅某点有取值，积分为 1；
  • 其具有 sifting property，即和某个函数作积分可以「筛选」出其某个点的函数值；

  • 离散状态下表现为仅在一点取值为 1；冲击序列 Impulse train，即一系列的在 \(k\Delta T\) 处的单位冲击 \[ s_{\Delta T}(t)=\sum_{-\infty}^\infty \delta(t-\Delta T) \]
  1. box function 的傅立叶变换，是一个 sinc 函数形式；证明的话利用欧拉公式直接推

  1. 常值函数 \(f(t)=1\) 的傅立叶变换是一个连续冲击函数，可利用逆变换的唯一性证明
  2. 冲击函数的傅立叶变换就相当于作用 shifting property；冲击串的傅立叶变换还是一个冲击串

卷积定理

  • 两个函数的卷积 \[ f(t)*h(t)=\int f(s)h(t-s)ds \]

  • 一维状态下的卷积定理 \[ f(t)*h(t)\Leftrightarrow H(u)F(u)\\ f(t)h(t)\Leftrightarrow H(u)*F(u) \] 即两函数卷积的傅立叶变换等于其对应的变换后函数的乘积；乘积的傅立叶变换等于量函数变换后的卷积。证明的话用到一次变元，在后面二维情况下进行了证明这里略去。

7.1离散频域变换

Sampling & Reconstruction

  • 采用上面提到的冲击串进行间隔为 \(\Delta T\) 的采样，于是采样得到的结果是下面的形式，注意到其可以表示为原函数和采样函数的乘积；

对于采样的结果进行傅立叶变换，利用到卷积定理和上面冲击串的卷积形式，可以得到如下结果

  • 注意到，对于这一结果，由于对于 n 从负无穷到正无穷进行了求和，所以是一个周期为 \(1/\Delta T\) 的周期函数；
    • 由于这一特性，假设原函数的傅立叶变换结果有上下限，若我们的采样率过小，可能会出现下面的这张图中的混淆现象

  • Sampling theorem （取样定理）

    带限函数：\(F(u), u\in[-u_{max}, u_{max}]\)

    奈奎斯特（Nyquist）采样率：\(2u_{max}\)

    采样间隔：\({1\over \Delta T}>2u_{max}\)

    • 在这样的情况下，我们可以根据变换后的结果重建出原信号而不损失信息
    • 恢复过程中，我们抽取出 \([-u_{max}, u_{max}]\) 区间内的部分并乘以 \(\Delta T\) 消除系数，就得到了原始的连续函数傅立叶变换结果；要从频域空间恢复时空域信号，施加逆变换即可；

DFT: discrete freq/Fourier transform

  • 在图像处理中，我们用的显然是离散状态下的傅立叶变换，即取一个周期内的 M 个样本点；注意到在离散状态下需要多出来一个系数以满足相互变换的一致性；

DFT for 2D

  • 2D impulse & sifting property 同样介绍了二维状态下的冲击函数
  • 2D continueous FT pair

2D sampling & sampling theorem 采样定理

  • 在二维状态下，冲击串 \(s_{\Delta T\Delta Z}(t,z)-\sum\sum \delta(t-m\Delta T, z-n\Delta Z)\) ，若取样间隔满足 \[ \Delta T<{1\over 2u_{max}}\\ \Delta Z<{1\over 2v_{max}} \] 则连续带限函数可由其一组样本基本无误地恢复。

  • DFT for 2D/3D，注意到同连续形式相比，同样多出来了一个系数

一些性质

没有细讲，不过在后面的处理当中必要的部分是将频谱图平移到中心，即这里的对于时空域图像做 \((-1)^{x+y}\) 变换

  • 频率域的频谱和相位角

下面可视化地分析了傅立叶变换后得到的结果中，频谱图和相位图的一些性质；我们记下面的一个 box 函数图像为 alpha，图 2 是中心化后的频谱图，图 3 进行了 Log 变换，图 4 是相位图；分别对于 alpha 做平移和卷积，傅立叶变换的结果也相应展示如下，在频谱上的结果是比较直观的，而我们很难从相位图上直观地观察到原图的一些信息；

然而，事实上，（在图像重建上）相位所带的信息要比频谱更为重要。这里记👇的图 1 为 A，图 2 为其相位图，图 3 是仅保留其相位，经过傅立叶逆变换之后重建的图像，可以看到一些边缘的信息；图 4 是仅保留其频谱信息重建出来的结果；更有意思的最后两张图，图 5 用了 A 的相位信息，而频谱没有置零而用了上面的 alpha 的，最后的结果可以模糊地看到小孩；图 6 用了 A 的频谱和 alpha 的相位，最后重建的结果从形态上更像 alpha。

二维状态下的离散卷积定理

要注意同样是在卷积定理中前面的一个系数

2 证明二维变量的离散傅里叶变换的卷积定理即：

\[ \begin{gathered} f(x,y)*h(x,y) <==> F(u,v)H(u,v)\\ f(x,y)h(x,y) <==> {1\over MN}F(u,v)*H(u,v) \end{gathered} \]

其中, \(*\) 表示卷积运算。

根据卷积和傅立叶变换公式： \[ \begin{gathered} F(u,v) = \sum_{x,y}f(x,y)e^{-j2\pi(ux/M+vy/N)}\\ f(x,y)*h(x,y) = \sum_{m,n}f(m,n)h(x-m,y-n) \end{gathered} \]

方便起见省略了加和的上标。参考一维情况，先证第一条

\[ \begin{aligned} \mathbf{F}\{f * g\}(u, v) &=\sum_{x, y}\left(\sum_{m, n} f(m, n) h(x-m, y-n)\right) e^{-j 2 \pi(u x / M+v y / N)} \\ &=\sum_{m, n} f(m, n) e^{-j 2 \pi(u m / M+v n / N)} \sum_{x, y} h(x-m, y-n) e^{-j 2 \pi(u(x-m) / M+v(y-n) / N)} \\ &=\sum_{m, n} f(m, n) e^{-j 2 \pi(u m / M+v n / N)} \sum_{i, j} h(i, j) e^{-j 2 \pi(u i / M+v j / N)} \\ &=F(u, v) H(u, v) \end{aligned} \]

即 \[ f(x,y)*h(x,y) <==> F(u,v)H(u,v) \]

下通过逆 FT 公式证第二条

\[ \begin{aligned} \mathbf F^{-1}\{\frac{1}{MN}F*H\}(x,y)&=\frac{1}{M N} \frac{1}{M N} \sum_{u, v}\left(\sum_{m, n} F(m, n) H(u-m, v-n)\right) e^{j 2 \pi(u x / M+v y / N)} \\ &=\frac{1}{M N} \sum_{m, n} F(m, n) e^{j 2 \pi(m x / M+n y / N)} \frac{1}{M N} \sum_{u, v} H(u-m, v-n) e^{j 2 \pi(x(u-m) / M+y(v-n) / N)} \\ &=\frac{1}{M N} \sum_{m, n} F(m, n) e^{j 2 \pi(m x / M+n y / N)} \frac{1}{M N} \sum_{i, j} H(i, j) e^{j 2 \pi(x i / M+y j / N)} \\ &=f(x, y) h(x, y) \end{aligned} \]

即 \[ f(x,y)h(x,y) <==>{1\over MN} F(u,v)*H(u,v) \]

8.1 Filters in frequency domain

  • 在时空域中，我们之前讲过 Filter 操作，如简单的平滑操作、Laplacian operator 等

  • Basics of filtering in frequency domain

    • 观察上面的二维傅立叶变换公式，我们可以看到：

    当 u=v=0 的时候，傅立叶变换得到的实际上是整个图像的平均灰度值； 当 u,v 较小的时候，低频，对应的是图像中灰度值变换比较缓慢的信息；

    • 例如，对于👇的这张图来说，我们可以看到频谱图的对角线方向有两个比较明显的高频分量，对应了原图像中的斜向的边缘（灰度值变化大）的信息。

  • 频域滤波公式

在图像域内，我们的滤波操作是采用了卷积操作；而根据卷积定理，相应的我们可以在频率域进行两张图像的点积操作，再将其变换到图像域，即上面的公式；

  • 下面直观的展示了低通和高通滤波器；低通滤波器可以达到 smooth 的效果；高通滤波器可以进行锐化；第三张图直接上图 2 的高通滤波器上加了一个常数，修正了最终处理结果的亮度；

另外，若直接对原图像进行操作，可能出现 缠绕错误 wraparound error ，因此需要将进行 zero padding

频率域步骤五步骤

  1. 给定大小为 \(M\times N\) 的图像 \(f(x,y)\) ，填充成大小为 \(P\times G\) 的图像 \(f_P(x,y)\) ，一般可取为原图像大小的 2 倍；
  2. 以 \((-1)^{x+y}\) 乘以 \(f_p\)，计算 DFT，得到中心化 \(F(u,v)\) ；
  3. 生成实对称中心化滤波函数 \(H(u,v)\)，\(G(u,v)=H(u,x)F(u,v)\) ；
  4. 计算 IDFT \(g_p(x,y)={real[F^{-1}(G(u,v))]}(-1)^{x+y}\) ，逆变换，取其实数部分，在进行去中心化操作；
  5. \(g_p\) 左上象限提取 \(M\times N\) 图像 \(g(x,y)\) ；

上述过程如下所示，图 1 为原图；图 2 进行了 padding；图 3 进行中心化；图 4 为傅立叶变换结果；应用图 5 的低通滤波器后得到图 6 的结果；逆变换得到图 7，再去除之前施加的中心化算子，取左上角的那部分图像得到图 8

Image smoothing

  • Ideal lowpass filter，理想低通滤波，即一个类似 box function 的圆柱

  • Ringing effect，直接用的话会出现 Ringing effect 效果；类似于时空域里的 box function 经过变换后得到 sinc 函数的特征，在频率域内的这种函数经过逆变换也会有类似的效果；

  • Gaussian lowpass filter GLPF，常用的是高斯低通滤波 \[ H(u,v)=e^{-D^2(u,b)/2\sigma^2} \] 设置不同的方差，可以得到不同的滤波效果；

Image sharpening

  • 高通滤波器，简单的可取 \(H_{HP}(u,v)=1-H_{LP}(u,v)\)
  • 于是有理想高通滤波，高斯高通滤波
  • 对于频率域中的 Laplacian 算子和 Highboost 滤波算子，在频率域中分别有如下的等效算子

选择性滤波

  • Bandreject and bandpass filters 带阻/带通滤波器

  • Notch filters 陷波滤波器

    用一个例子来进行说明，我们得到的图片 1 可能出现一些条纹，这在时空域显然是很难去除的；我们转换为频率域得到图 2，横向条纹对应的频谱图分量类似于图 3；因此设计图 4 的滤波器；滤波结果如图 5 所示；图 6 是对于图 3，也即横向条纹对应的频谱分量进行逆变换，直观地展示了这些杂纹；

    作业中用到了这个算法，不过找到对应的频谱分量真的有点难 Orz，可能是基于 Python 的手动实现的原因吧，这时候才意识到在图像处理领域中，交互式的设计有多重要；

9.0 Spatial Transform 图像变换

Linear transformation

  • Rigid transformation 刚性变换：Rigid transformation preserves the angels and distances within the model

    包括了平移和旋转

  • Similarity transformation 相似变换：Similarity transformation adds scaling {s}

    增加了伸缩

  • Affine transformation 仿射变换：Affine transformation applies a function between affine spaces which preserves points, straight lines and planes.

    增加了 Shear？

    可以用矩阵来表示；保持平行关系

Nonlinear transformation

  • Thin-Plate Spline 薄板样条
  • FFD Free-form deformation 自由形变变换
  • Locally affine 局部仿射

FFD Free-form deformation 自由形变变换

  • 想法还是一个局部效应求和的过程；设定网格状的控制点 P，每个控制点可以发生一个形变 \(\Phi\)；这里的想法就是 \(X'=X+Q\)，其中的「形变」部分由 X 周围的一些控制点决定；

常见的是三次 B 样条核函数，一个类似于高斯分布的形式（忽略里面的导数曲线）；可以看到 B 样条函数只在 \([-2,2]\) 区间内取值；对于一个目标点 X，划出这个区间可以找到与其相邻的四个点，对着四个点的形变量 \(\Phi\) 以 X 到控制点的距离所对应的 B 样条核函数的值为权重（注意这四个点之和=1）求和；最后将形变量加到原来的 X 上即可；

Locally affine 局部仿射

  • 这里可以在多个区域内进行仿射仿射变换，作为看参考的区域，对于此外的像素点，按照其离这些参考区域的距离倒数的 e 次方衰减作为权重，把这些变换汇总起来；实现各个区域的局部仿射在算法上还是比较复杂的，作业里面用的是一些参考点而非区域。下面是一些效果图

Transform an image

  • Forward VS Backward (inverse)

  • 在具体实现的时候，若是前向进行，即对原图像上的每一个像素点，做 T 变换到新的图像，则可能出现其中很多的点没有被映射到的问题；

  • 因此，更多时候用的是反向图变换，对于目标图像的每一个像素点位置，填充原图像上 \(T^{-1}(x')\) 位置的颜色。注意变换之后需要用到插值。

  • 作业实现了局部仿射的逆变换，摘录报告中的表述如下

    具体而言，算法的输入包括人脸和猩猩脸两张图片，需要利用两者的相关特征，将人脸 \(X\) 风格化为猩猩的面部特征，记输出图像为 \(Y\)。考虑到图像存储的离散性，采用逆向变换，即若 \(X\) 到 \(Y\) 的变换为 \(Y=T(X)\)，则算法实现其逆变换 \(S(Y)=T^{-1}(Y)=X\) 。

    对于原图像的若干控制点 \(X_i\)，设其经过各自变换后新的坐标点为 \(Y_i=G_i(X_i)\) ，则基于局部仿射变换，对于图像中的任意点 \(X\)，其变换结果为个变换的加权平均，权重根据该点与各控制点的距离衰减。即 \[ T(X)=\sum_{i=1}^n w_i(X)G_i(X_i)\\ w_i(X) = {1/d_i(X)^e\over \sum 1/d_i(X)^e}\\ \] 其中 \(d_i(X)\) 表示 \(X\) 到各控制点 \(X_i\) 之间的距离，权重根据距离倒数的 \(e\) 次方衰减，并进行归一化。

    记 \(G_i(X_i)\) 的逆变换为 \(G'(Y_i)\) ，则逆变换过程可近似表为 \[ S(Y)=\sum_{i=1}^n w'_i(Y)G'_i(Y_i)\\ \tag{2.1} w'_i(Y)={1/d_i'(Y_i)^e\over \sum1/d_i'(Y_i)^e} \] 其中 \(d'_i(Y_i)\) 为 \(Y\) 到各控制点 \(Y_i\) 的距离。

10.0 颜色和彩色图像

  • 颜色、色彩空间

    • 光谱
    • 人眼：视杆细胞、视锥细胞
    • 基色：RGB、红黄蓝、青品红黄黑 CMYK（Cyan, Magenta）
    • 颜色视觉障碍
    • 色彩空间
      • RGB 加色法系统 采用笛卡尔坐标系定义颜色，用处如显示器
      • CMYK 减色法系统 白色减去某个颜色，多用于打印机
      • HSL/HSI 颜色空间
        • Hue 色调，人类认为的颜色
        • Saturation 饱和度，纯度，与灰色的距离
        • Lightness/Intense 亮度，从黑色到亮色
      • Lab颜色空间<L,a,b>
        • 感知上均匀分布
        • a, b用来近似“红/绿” 和“黄/蓝”通道
    • 色彩转换：颜色空间的转换 #Further reading

彩色图像处理

  • 根据不同的色彩空间，需要用不同的函数进行变换，例如要对降低一张图像的亮度，在 RGB 空间中对每一个通道进行收缩即可；仅是 CMY 减色法系统中，统一增加数值即可（加起来是黑色）；而在 CMYK 空间中，仅需要调节 K 通道；在 HSL 空间中，单独调节亮度通道即可。

  • 而要计算补色，在 RGB 空间内转置即可；而在 HSI 空间中，将亮度转置（白变黑），纯度不变，而对于色调，想象一个色环，我们要将其映射到「对面」去，或者说变换保持 1/2 的距离，如上图所示。

  • 注意到，一些算法对于不同的色彩空间需要有不同的设计；例如，对于基于 Histogram 的算法，分别在 RGB 通道上操作会有问题，可在 HSI 空间中进行，仅调节亮度 I。

  • 图像的 Filter 操作，下图展示了在 RGB 空间和 HSI 空间的结果：

10.1 Plot - basics

  • 视觉基础

    • 视觉：低阶（物理性质）、高阶（识别、分类等，属于认知能力重要组成部分）

    • 格式塔 Gestalt 理论（可以指导视觉设计）

      • 强调经验和行为的整体性

      感知的事情大于眼睛见到的事物

      • 贴近原则

      相似原则

      对称性原则

      经验原则

      连续原则

      闭合原则

      共势原则 common fate，方向同？

      好图原则 good figure

    • 视觉感知的相对性：相对性&绝对性，介绍了而一些常见的错视效应

  • 视觉编码

    • 标记和视觉通道

      • 标记：是数据属性到可视化元素映射（图形元素：点线面）
      • 视觉通道：是数据的值到标记的视觉表现属性的映射（位置、大小、形状、方向、色调、饱和度、亮度……）

    • 视觉通道的概念

      • 表现力：视觉通道在编码数据信息时，需要表达且仅表达数据的完整属性

      要求视觉通道准确编码数据包含的所有信息；在对数据进行编码的时候，需要尽量忠于原始数据

      • 不同类型的数据选取不同的通道进行可视化

  • 视觉通道的特性

    • 不同的视觉通道，人眼对其的感知能力是不同的；如尺寸/长度是先行感知的话，平面位置、颜色是超线性的，角度、形状等式亚线性的。

11.0 Plot

  • 统计图介绍

    • 部分占整体比例：饼状图

    • 比较与比例：柱状图为例

    • 堆叠图：相对于柱状图，可以同时加上时间维度

    • 趋势与模式：折线图

    • 散点图

    • 盒须图/箱线图

    • 其他…

  • 设计原则

    • 数据到可视化的直观映射：确定数据类型、合适编码方式、优先级
    • 视图选择与交互设计：滚动和缩放、颜色修改、更改映射方式
    • 信息密度--数据的筛选
    • 美学因素：坐标轴刻度、网络；标注、图例、标题；正确使用颜色
    • 颜色与透明度：颜色的混合（线性加和）
    • 可视化隐喻：拟物？
    • 动画与过渡：时间换取空间

11.1 Map

  • 地图和地图投影
    • 地理坐标表示：经纬度系统
    • 经纬度系统，实质上是用二维坐标描述了一个三维球体表面的位置；然而我们要可视化到二维的屏幕上，就涉及到如何选择投影；根据不同的优化原则，有不同的选择（等角度/正形投影、等面积、等距离/等方位角）
  • 著名的墨卡托投影，可以保证等角度。因此多用于各种航海地图，导航地图（百度地图，高德地图）；然而，由于不具备等面积特性，所以在远离赤道的区域，两个点之间的距离会远远超过真实值，容易给人带来混乱

  • 亚尔勃斯投影，又称等面积圆锥投影，由德国人亚尔勃斯于 1805 年发明，按照原则将地图投影到圆锥面上；由于其保持了等面积，常用于表现国家疆域面积大小的图示中

  • 等方位角 azimuth 投影，属于等距投影的一种，地图上任何一点沿着经度线到投影中原点的距离保持不变。

  • https://en.wikipedia.org/wiki/List_of_map_projections

12.0 地理数据可视化

  • 点数据可视化
  • 线数据：长度属性、连接关系
  • 区域数据可视化，简单介绍了一些
    • Choropleth 地图（等值区间地图），在区域内部均匀分布，填上一致的颜色；有时候的问题在于，区域内的数据和数据大小不一致；
    • 于是有了 Cartogram，克服了对空间实用的不合理性，它按照地理区域的属性值对各个区域进行适当的变形；有些不同的形变

12.1 地图地理数据可视化实现

初步讲了两类

  • Python Basemap+Matplotlib
  • Baidu Echart and PyEchart

12.2 空间标量场可视化介绍

  • 一维数据
  • 二维数据：颜色映射、高度映射、等值线映射 iso-contour
  • 三维标量场
    • 截面可视化：如三个正交平面或任意角度切片
    • 间接体绘制：等值面提取
    • 直接体绘制

13.0 空间标量场可视化 等值线/面的计算

  • 等值线 (isocontour)
    • 类似于线性插值的想法，若一条边两端分别大于和小于所关注的值，线性找到所需要的点；把一个正方形中边上至少两个等值线上的点连起来；
    • 分治思想，分出正方形单元；每个点和 C 比较；线性插值；连接线
  • 等值面 (isosurface)
    • 根据不同的点大小分布，等值线计算有 4 种不同的情况；而等值面的计算中，立方体共有 14 中不同的情况。注意这时候，一个 cell 里的分割面不是平面而是一些三角面片。
  • Marching Cubes Algorithm MC
    • 对于各点相较于 C 的大小，变成一个二值向量，映射到最终的分割方式；
    • 需要注意的是歧义性问题 Ambiguity，例如在二维情况下，对角线上两组点分别是大于/小于 C 的话，可能出现两种分割方式；在三维情况下，可能出现不匹配的问题，从而出现「洞」；解决算法 #further reading

14.0 直接体渲染

  • 直接计算最终可视化里的每一个像素，体绘制中所有体素对最后的图像亮度都有贡献；因此，其在计算量上的需求要高于面渲染

  • 基于光线投射的直接体绘制方法 Ray casting

    • 采样重建：体采样

    数据分类：体分类，传输函数设计

    光照计算：体光照模型

    光学积分：体积分

  • 采样的话，除了网格点之外，更为理想的自适应采样，即平缓均匀区域增大采样间隔

    • 对于采样率的要求：每个体素至少需要 2 个采样点（之前的定理）

  • 体分类和传输函数

    • 我们拿到的数据多为灰度值数据，因此我们需要决定每个体素吸收的值（不透明度α）和发射的部分颜色(RGB）
    • Transfer funciton: 一组定义了数据值及其相关属性，不透明度等视觉元素之间的映射关系的函数
    • 主要是其中不透明度函数的设计，原则是：高不透明度对应重要特征、感兴趣的区域；而颜色的设计可以帮助我们进行区域分割；
    • 我们可以利用体素的标量值、梯度模等信息构建传输函数，当然还可以加入一些其他的信息帮助设计；
    • 例如，用上面的两个属性构建二维直方图，进行 1. 映射规则设计；2. 光学属性设计。
    • 注意到传输函数的设计很大程度上决定了最终渲染的好坏，用户交互的传输函数是一个很好的方案。

  • 体光照模型

    • 面绘制中的 Phong 光照模型：最终的 Result color = Ambient color + Diffuse color + Specular color 即环境光、散射光和反射光的共同作用；

  • 从公式可以看到，最终的颜色三个因素都有影响；环境光仅和物体的属性相关；漫反射和物体法线与光线方向、物体颜色有关，而与视线方向无关；而镜面反射光与光照方向、物体法线方向、视线方向头有关，而与物体颜色无光，还有一个高光系数。

  • 这里需要计算法向量，类似之前的「边缘检测」可以表示为 spatial correlation 的形式；当然，更为常用的是中心差分，kernel 仅仅是概念上的存在；

  • 体积分

  • 这里是把视线/屏幕作为原点，连续情况下是进行积分，对于每一个点的颜色值，其权重为从该点发射的光线，经过了前面的所有体素阻隔之后的剩余部分，即式中指数上面对于不透明作积分，取负指数衰减；

  • 这是一种思路，而在下面的离散情况下（似乎和连续情况并非对应关系），积分变为求和，这里的视角从「不透明度」变为「透明度」，直接以连乘的透明度作为权重作用到发射点的颜色值上；

  • 工程上，显然直接用上式求 \(D_M\) 之前的所有点的一个体积分需要两次循环，复杂度 \(O(M^2)\) ；下面讲了两种简化的思路：

  • 一种是「从后向前积分」，采用 \(D_i=C_i+D_{i+1}(1-a_i)\) ，物理意义是从后向前不断有光线叠加，并透过「玻璃」吸收掉一部分，累加到最前面，这样的复杂度变为线性；

  • 更高效或许是「从前向后积分」 \[ D_i=D_{i-1}+C_i(1-A_i)\\ A_i=A_{i-1}+(1-A_{i-1})a_i=a_i+A_{i-1}(1-a_i) \] 这里维护了一个 \(A_i\)，是从前向后累计下来的对当前像素的「总吸收率」，我们按照从前向后的顺序进行加总；这里的好处是，当总吸收率较大的时候，我们就不需要再计算后面的体素了。下面是伪代码

  • 根据合成顺序的不同，体绘制方法分为
    • 图像空间扫描的体绘制方法（屏幕），想象屏幕发出视线
    • 物体空间扫描的体绘制方法（数据），这里是针对数据的每个像素点，叠加到屏幕上
  • 上面的光线投影法数据图像空间方法，其他的还有一些比较直观的算法

  • 这里的 First，即第一次遇到感兴趣的数值 C 时停下来，记录该点深度；也就是通过体绘制的方式所构成的等值面；

  • 取 Max，在医学影像中用到如最大密度投影法，进行血管造影/显像。

  • 下面的「splitting」属于物体空间扫描的体绘制方法，想法也很直观

总结如下

14.1 三维数据场渲染实现 VTK + NIFTY + Python可视化

  • 这里讲了 VTK http://www.vtk.org/，当然用的 Python 版本

  • 面渲染，用了 MC 算法

  • 体绘制

    直接调用的 vtkVolumeRayCastMapper 类，相关的设置挺流程化的，所以真正实现的时候最为关键的应该是传输函数的设计；实验的时候似乎环境光系数、漫反射系数、镜面反射系数等参数也挺重要的。

15.0 向量张量场数据可视化

  • 这一讲简单介绍了张量场的可视化技巧

张量场可视化

  • 向量场的可视化：包括图标法、几何法、纹理法、拓扑法等

图标法

  • 采用的图标主要有：线条、箭头、方向标志符，经典方法：hedge hogs（刺猬），用带方向线段来表示矢量场的一个点图标；
  • 此外，我们可以用其他的视觉通道表示向量场的其他信息，例如对于向量的方向，不同的方向赋予一些颜色

  • 这一算法的优点在于实现简单、直观、灵活；然而，缺点也很明显：
    • 可视混乱
    • 无法解释数据的内在连续性
    • 难以表达结构特征如涡流等

几何法

  • 几何法指采用不同类型的几何元素，如线、面、体模拟向量场的特征
  • 这里讲了三类曲线的可视化：流线、迹线、脉线
  • 流线 stream line 主要面向的是稳定的向量场，或许不稳定向量场的某一时刻；即在一个静止的向量场中，在某处放下一个点（不考虑质量？），其运动轨迹；
  • 迹线 path line 是在一个动态的向量场中，某处放置一个点，其运动轨迹；
  • 脉线 streak line 的意思是，在某一个点持续释放粒子，经过一段时间之后，各个粒子所连成的曲线（在稳定场中，迹线和脉线相同？）；因此这里和上面的两种「轨迹」不太一样，可以想象「烟圈」；
  • 设在 0-3 时刻向量场的方向分别如下所示，则流线迹线脉线区别如下

纹理法

  • 点噪音法（spot noise），按照局部流场方向对圆点变形
  • 线积分卷积（line integral convolution, LIC），在流线方向进行卷积操作

张量 Tensor 场可视化

  • 主要面向二阶（对称）张量场
  • 想法是对其作特征分解；根据三个特征值的大小关系，可以有 spherical, linear, planar 三种情况
  • 几何图标法，用球、椭球、线等方式可视化这些特征值

  • 几何：纤维追踪法 Fiber Tracking #further reading

多变量空间数据场可视化

  • 降维
  • 对于特定的应用，必须根据数据自身的特点为其选择和设计适用的特征抽取方法和可视化表达方式。
  • 例如，选取不同的视觉通道
  • 图标法可清晰编码多个变量属性，尤其是可表达方向信息