python中的scipy库

Scipy建立在Numpy的基础之上，用于数值运算.具有很多搞笑操作,如数值积分、优化、统计、信号处理，以及图像处理功能。

自己在阅读书籍或代码的时候遇到的一些函数，记录在此，后面或有更新。

图像模糊处理

图像模糊的本质：图像I与一个高斯核$G_{\sigma}$进行卷积

$$I_{\sigma}=I*G_{\sigma}$$
$$G_{\sigma}=\frac{1}{2\pi{\sigma}^2}e^{-\frac{(x^2+y^2)}{2{\sigma}^2}}$$
滤波操作模块：scipy.ndimage.filters，使用快速一维分离方式计算卷积

from PIL import Image
from numpy import  *
from scipy.ndimage import filters

im = array(Image.open('sample.jpg').convert('L'))
im2 = filters.gaussian_filter(im, 5) # 标准差为5

若是模糊彩色图像，则对每个颜色通道进行高斯模糊

im = array(Image.open('smaple.jpg'))
im2 = zeros(im.shape)
for i in range(3):
    im2[;, ;, i] = filters.gaussian_filter(im[;, ;, i], 5)

im2 = uint8(im2) # 或者 im2 = array(im2, 'uint8')

图像导数

强度的变化可用灰度图像I（彩色图像则单独对每个通道计算导数）的x、y方向导数Ix、Iy描述。

图像的梯度向量$\nabla I = [I_x, I_y]^T$，其大小描述了图像强度变化的强弱：

$$|\nabla I|=\sqrt{I^2_x+I^2_y}$$

角度描述了图像中在每个像素上强度变化最大的方向，$[-\pi, \pi]$：

$$\alpha = arctan2(I_x, I_y)$$

图像导数通常可用卷积实现

$$I_x = I*G_{\sigma x}\space\space I_y = I*G_{\sigma y}$$

$G_{\sigma x},G_{\sigma y}$为$G_{\sigma}$不同方向上的导数

# 简单处理
from PIL import Image
from numpy import *
from scipy.ndimage import filters

im = array(Image.open('sample.jpg').convert('L'))
sigma = 5 # 标准差

imx = zeros(im.shape)
# 第二个参数：使用的标准差
# 第三个参数：指定对每个方向计算哪种类型的导数
# 第四个参数：保存输出的变量
filters.gaussian_filter(im, (sigma, sigma), (0,1), imx)

imy = zeros(im.shape)
filters.gaussian_filter(im, (sigma, sigma), (1,0), imy)

magnitude = sqrt(imx**2+imy**2)

读写.mat文件

使用scipy.io模块进行读取

• 保存一般数组变量

data = scipy.io.loadmat('test.mat')
# data对象包含一个字典，字典中的健对应保存在原始.mat文件中的变量名

# 保存时使用savemat()函数
data = {}
data['x'] = 12
scipy.io.savemat('test.mat', data)

• 以图像形式保存数组

# 将数组im直接保存为图像文件
from scipy.misc import imsave
imsave('test.jpg', im)

图像去噪

在去除噪声的同时，尽可能地保留图像细节和结构

这里采用ROF（Rudin-Osher-Fatemi）去噪模型，ROF处理后的图像更平滑，同时保持图像边缘和结构信息
ROF模型简介
一副（灰度）图像I的全变差$J(I)$（Total Variation， TV）定义为梯度范数之和。

$$(连续表示情况) J(\mathbf I) = \int{|\nabla I|d\mathbf x}$$
$$(离散表示情况) J(\mathbf I) = \Sigma_{\mathbf x}{|\nabla I|}$$
x为所有图像坐标，在Chambolle提出的ROF模型中，目标函数为 寻找降噪后的图像 U，使得下面表达式最小：
$$min_{U}\Vert{I-U}\Vert^2 + 2\lambda J(\mathbf U)$$
其中 $\Vert{I-U}\Vert$是去噪后图像与原始图像差异的度量，本质上是该模型使去噪后的图像像素值”平坦”变化，但在图像区域的边缘上，允许去噪后的图像像素值”跳跃”变化

# ROF模型实现 rof.py
#! -*- coding:utf-8 -*-
def denoise(im, U_init, tolerance=0.1, tau=0.125, tv_weight=100):
    """使用Chambolle在公式（11）中的计算步骤实现ROF
    input：含有噪声的输入图像（灰度图像）、U的初始值、TV正则项权值、步长、停止条件
    output：去噪和去除纹理后的图像、纹理残留
    """
    m,n = im.shape

    # 初始化
    U = U_init
    Px = im # 对偶域x分量
    Py = im # 对偶域y分量
    error = 1

    while(error < tolerance):
        Uold = U
        # 原始变量的梯度
        GradUx = roll(U, -1, axis=1)-U # 变量U梯度的x分量
        GradUy = roll(U, -1, axis=0)-U # 变量U梯度的y分量

        # 更新对偶变量
        Pxnew = Px + (tau/tv_weight)*GradUx
        Pynew = Py + (tau/tv_weight)*GradUy
        NormNew = maximum(1, sqrt(Pxnew**2 + Pynew**2))

        Px = Pxnew/NormNew # 更新x分量(对偶)
        Py = Pynew/NormNew # 更新y分量(对偶)

        # 更新原始变量
        RxPx = roll(Px, 1, axis=1) # 对x分量进行向右x轴平移
        RyPy = roll(Py, 1, axis=0) # 对y分量进行向右y轴平移

        DivP = (Px - RxPx) + (Py - RyPy) # 对偶域的散度
        U = im +tv_weight*DivP # 更新原始变量

        # 更新误差
        error = linalg.norm(U-Uold)/sqrt(n*m)

    return U, im-U # 去噪后图像及纹理残余

ROF测试

#! -*- coding:utf-8 -*-

from pylab import *
from numpy import  *
from numpy import random
from scipy.ndimage import filters
from scipy.misc import imsave
from rof import *

im = zeros((500,500))
im[100:400,100:400] = 128
im[200:300,200:300] = 255
im = im + 30*random.standard_normal((500,500))

U, I = rof.denoise(im, im)
G = filters.gaussian_filter(im, 10)

# 保存结果
imsave('synth_original.pdf',im)
imsave('synth_gaussian.pdf',G)
imsave('synth_rof.pdf',U)

# 绘图
figure()
gray()

subplot(1,3,1)
imshow(im)
axis('equal')
axis('off')

subplot(1,3,2)
imshow(G)
axis('equal')
axis('off')

subplot(1,3,3)
imshow(U)
axis('equal')
axis('off')

show()