EE103(2)：图像线性处理

不使用opencv的库函数，完全运用numpy和线性代数知识，我们能在多大程度上处理图像呢

多通道转单通道

• 可以使用(r,g,b)=(1/3,1/3,1/3)来转换
• 也可以使用(r,g,b)=(0.299,0.587,0.114)来转换

原图（pixiv id：490219）：

两种转换方式结果对比

亮度、对比度

$y=ax+[(1-a)avg(x)+b]\cdot 1$

• x为像素矩阵
• 1为全1矩阵
• a越大，对比度越高
• b越大亮度越高

不同的a（0、0.4、0.8、1.2、1.6、2）情况下图像的情况（不同b情况下略）：

变换

• 水平缩放：对于m*n的原图像，右乘k*m的矩阵A可以水平缩放为k*n的图像。其中，A的纵列代表了缩放过程中对原图像横排每个像素的加权，A的横行代表了最终图像横排每个位的生成，原图像的纵列代表了最终图像纵列每个位的生成。在缩放过程中，最终图像横排每个位的生成次数应该和原图像不同，并且每次生成时应该对原图像横排不同像素进行加权

  • 临近插值放大一倍的矩阵A：$\left[\begin{array}{ccccccc}1& 1& 0& 0& ...& 0& 0\\ 0& 0& 1& 1& ...& 0& 0\\ 0& 0& 0& 0& ...& 0& 0\\ ⋮& & & & & & \\ 0& 0& 0& 0& ...& 1& 1\end{array}\right]$
  • 线性插值放大一倍的矩阵A：$\left[\begin{array}{ccccccc}1& 0.5& 0& 0& ...& 0& 0\\ 0& 0.5& 1& 0.5& ...& 0& 0\\ 0& 0& 0& 0.5& ...& 0& 0\\ ⋮& & & & & & \\ 0& 0& 0& 0& ...& 0.5& 1\end{array}\right]$

• 垂直缩放：对于m*n的原图像，左乘n*k的矩阵A可以垂直缩放为m*k的图像。其中，A的横行代表了缩放过程中对原图像纵列每个像素的加权，A的纵列代表了最终图像纵列每个位的生成，原图像的横行代表了最终图像横行每个位的生成。在缩放过程中，最终图像纵列每个位的生成次数应该和原图像不同，并且每次生成时应该对原图像纵列不同像素进行加权

  • 临近插值放大一倍的矩阵A：$\left[\begin{array}{ccccc}1& 0& 0& ...& 0\\ 1& 0& 0& ...& 0\\ 0& 1& 0& ...& 0\\ 0& 1& 0& ...& 0\\ ⋮& & & & \\ 0& 0& 0& ...& 1\\ 0& 0& 0& ...& 1\end{array}\right]$
  • 线性插值放大一倍的矩阵A：$\left[\begin{array}{ccccc}1& 0& 0& ...& 0\\ 0.5& 0.5& 0& ...& 0\\ 0& 1& 0& ...& 0\\ 0& 0.5& 0.5& ...& 0\\ ⋮& & & & \\ 0& 0& 0& ...& 0.5\\ 0& 0& 0& ...& 1\end{array}\right]$

• 旋转、翻转：矩阵旋转与转置

• 偏移：切片与连接

常用图层混合模式

设像素矩阵$x$位于像素矩阵$y$之上且范围比y小，坐标$i,j\in x坐标范围$

• 正常：$y_{ij}=x_{ij}$

• 变暗系

  • 变暗： y i j = m i n { y i j , x i j } y_{ij}=min\{y_{ij},x_{ij}\} yij​=min{yij​,xij​}，分通道比较替换
  • 正片叠底：$y_{ij}=y_{ij}\ast x_{ij}/255$
  • 颜色加深：$y_{ij}=(y_{ij}+x_{ij}-255)\ast 255/x_{ij}$
  • 线性加深：$y_{ij}=y_{ij}+x_{ij}-255$
  • 深色： y i j = m i n { y i j , x i j } y_{ij}=min\{y_{ij},x_{ij}\} yij​=min{yij​,xij​}，比较全通道数值

• 变亮系

  • 变亮： y i j = m a x { y i j , x i j } y_{ij}=max\{y_{ij},x_{ij}\} yij​=max{yij​,xij​}，分通道比较替换
  • 滤色：$y_{ij}=255-(255-y_{ij})\ast (255-x_{ij})/255$
  • 颜色减淡：$y_{ij}=y_{ij}+x_{ij}\ast y_{ij}/(255-y_{ij})$
  • 线性减淡：$y_{ij}=y_{ij}+x_{ij}$
  • 浅色： y i j = m a x { y i j , x i j } y_{ij}=max\{y_{ij},x_{ij}\} yij​=max{yij​,xij​}，比较全通道数值

• 叠加：根据y的数值来决定使用正片叠底还是滤色，尽量保留y的颜色$\{\begin{array}{l}{y}_{ij}=x_{ij}\ast y_{ij}/128y_{ij}\le 128\\ y_{ij}=255-(255-x_{ij})\ast (255-y_{ij})/128y_{ij}>128\end{array}$

• 照射系

  • 柔光：对x数值高的部分保留较好，$\{\begin{array}{l}{y}_{ij}=y_{ij}+(2x_{ij}-255)(y_{ij}-y_{ij}^2/255)/255x_{ij}\le 128\\ y_{ij}=y_{ij}+(2x_{ij}-255)(\sqrt{y_{ij}/255}-y_{ij})/255x_{ij}>128\end{array}$
  • 强光：根据x的数值来决定使用正片叠底还是滤色，尽量保留x的颜色$\{\begin{array}{l}{y}_{ij}=x_{ij}\ast y_{ij}/128x_{ij}\le 128\\ y_{ij}=255-(255-x_{ij})\ast (255-y_{ij})/128x_{ij}>128\end{array}$
  • 亮光：尽量保留x的颜色，x亮度高时图像将被降低对比度并且变亮，x亮度低时图像会被提高对比度并且变暗$\{\begin{array}{l}{y}_{ij}=255-(255-y_{ij})/(2\ast x_{ij})x_{ij}\le 128\\ y_{ij}=y_{ij}/(2\ast (255-x_{ij}))\ast 255x_{ij}>128\end{array}$
  • 线性光：根据x的数值来决定使用线性加深还是线性减淡，尽量保留x的颜色$y_{ij}=2\ast x_{ij}+y_{ij}-255$

模糊

生成一个$M\ast N$的模糊用卷积核x，与原图像y进行卷积运算，即$y_{ij}=\sum _{m=0}^M\sum _{n=0}^N{x}_{mn}{y}_{i-m,j-n}$

结果的每个ij位置都相当于原图像ij位置左上$M\ast N$范围内的像素加权和，因此可以实现模糊

高斯核

生成原理很简单，对于每个位置计算其与右下角横纵距离，代入二维高斯函数即可。sigma越大模糊程度越大

def gaussian(u, v, sigma):
    pi = 3.1416
    intensity = 1 / (2.0 * pi * sigma * sigma) * math.exp(- 1 / 2.0 * ((u ** 2) + (v ** 2)) / (sigma ** 2))
    return intensity

def gaussianKernel(r, sigma):
    kernel = np.zeros([r, r])
    for i in range(r):
        for j in range(r):
            kernel[i, j] = gaussian(r-i,r-j, sigma)
    kernel /= np.sum(np.sum(kernel))
    return kernel

运动模糊核

运动模糊的原理是一个方向上的图像叠加，给出运动模糊的角度和程度，需要生成一个权重和为1的核。对于向右下、左下、右上、左上运动，核分别越靠近左上、右上、左下、右下的权重最大，其中对角线上r个格子中有值

设首项等于公差，解方程$1=rd+\frac{r(r-1)}{2}d$，得$d=\frac{2}{r+r^2}$

生成核时，从根据不同的弧度情况，从左上、右上、左下、右下分别开始向着弧度对应的方向寻找r个格并在对应点计算等差数列对应项，如果该格已经有值，则在其旁边一格设置项

def mbkernel(agl,r):
    rad=math.radians(agl)
    w=math.ceil(r*np.abs(np.cos(rad)))
    h=math.ceil(r*np.abs(np.sin(rad)))
    if h==0:
        h+=1
    if w==0:
        w+=1
    d=2/(r+r**2)
    kernel=np.zeros((int(h),int(w)))
    for i in range(0,r):
        if np.sin(rad)>=0:
            y=math.floor(h-i*np.sin(rad)-1)
        else:
            y=math.floor(-i*np.sin(rad))
        if np.cos(rad)>=0:
            x=math.floor(i*np.cos(rad))
        else:
            x=math.floor(w+i*np.cos(rad)-1)
        if kernel[y][x]==0:
            kernel[y][x]=(i+1)*d
        else:
            if np.sin(rad)>=0:
                kernel[y+1][x]=(i+1)*d
            else:
                kernel[y-1][x]=(i+1)*d
    return kernel

卷积

opencv提供了卷积函数

dst=cv2.filter2D(src,-1,kernel=k,anchor=(-1,-1))	#anchor默认位于中间，深度为-1默认保持和原图像相同深度
dst=cv2.filter2D(src,-1,kernel=flip(k),anchor=(k.cols-1,k.rows-1))	#真正意义上的卷积，深度为-1默认保持和原图像相同深度

不使用opencv提供的卷积函数，同样可以手写卷积，但运行速度极慢

半径（纵轴）分别为5、15、30下，$\sigma$（横轴）分别为10、50、100下高斯模糊的结果

半径（纵轴）分别为15、50、100下，角度（横轴）分别为$0^o$、$45^o$、$135^o$下运动模糊的结果

差分

水平差分

对于m*n的原图像，右乘m*m的矩阵A可以进行水平差分。其中，A的纵列代表了缩放过程中对原图像横排每个像素的加权，A的横行代表了最终图像横排每个位的生成，原图像的纵列代表了最终图像纵列每个位的生成。在差分过程中，最终图像横排每个位的生成次数应该和原图像相同，并且每次生成时应该对原图像横排不同像素进行加权（某个像素为1，其下一个像素为-1）

• 矩阵A：$\left[\begin{array}{cccccc}1& 0& 0& ...& 0& 0\\ -1& 1& 0& ...& 0& 0\\ 0& -1& 1& ...& 0& 0\\ ⋮& & & & & \\ 0& 0& 0& ...& -1& 1\end{array}\right]$

垂直差分

对于m*n的原图像，左乘n*n的矩阵A可以进行垂直差分。其中，A的横行代表了缩放过程中对原图像纵列每个像素的加权，A的纵列代表了最终图像纵列每个位的生成，原图像的横行代表了最终图像横行每个位的生成。在缩放过程中，最终图像纵列每个位的生成次数应该和原图像不同，并且每次生成时应该对原图像纵列不同像素进行加权（某个像素为1，其下一个像素为-1）

• 矩阵A：$\left[\begin{array}{ccccc}1& -1& 0& ...& 0\\ 0& 1& -1& ...& 0\\ 0& 0& 1& ...& 0\\ ⋮& & & & \\ 0& 0& 0& ...& -1\\ 0& 0& 0& ...& 1\end{array}\right]$

图像复原

对我们的图像添加10%椒盐噪声

base = cv2.imread("base.jpg",cv2.IMREAD_GRAYSCALE)
noise=np.random.randint(0,100,base.shape)
base=np.where(noise<5,0,base)
base=np.where(noise>=95,255,base)

添加前：

添加后：

均值滤波

设半径为r，均值滤波核是一个$r\cdot r$，每个位置值为$\frac{1}{r\cdot r}$的方阵

kernel=np.zeros((r,r))
kernel+=1/r/r

分别使用半径为3、5、15的核滤波结果

中值滤波

设半径为r，中值滤波核在图像上滑动，每$r\cdot r$个位置统计中值

def med(inputs,R):
    H, W = inputs.shape
    result = np.zeros(inputs.shape)
    H, W = inputs.shape
    for r in range(0, H - R + 1):
        for c in range(0, W - R + 1):
            cur = inputs[r:r + R,c:c + R]
            val = np.median(cur)
            result[r, c] = val
    return result

分别使用半径为3、5、15的核滤波结果