A 依据颜色的k均值聚类

def _1rgb_kmeans(picname='', cutnum=50, clus=5):
    '''
    彩色图像按颜色k-means聚类.
    依赖:from scipy.cluster import vq
        from scipy.misc import imresize #This function is only available if Python Imaging Library (PIL) is installed.
        from PIL import Image
    :param picname: 图像名,字符串
    :param cutnum: 切成cut*cut个块
    :param clus: 分几类
    :return: 无,显示图像
    '''

    from scipy.misc import imresize
    import  scipy.cluster.vq    # 包导入只会在第一次导入时执行
    gridsize = cutnum   # 多少格
    im_a = np.array(Image.open(picname))
    x_num = int(im_a.shape[0] / gridsize)  # 行方向上的像素数
    y_num = int(im_a.shape[1] / gridsize)

    grids = []
    for x in range(gridsize):
        for y in range(gridsize):
            R = np.mean(im_a[x*x_num: (x+1)*x_num, y*y_num: (y+1)*y_num, 0])   # 区域内的像素均值,红色
            G = np.mean(im_a[x * x_num:(x + 1) * x_num, y * y_num:(y + 1) * y_num, 1])
            B = np.mean(im_a[x * x_num:(x + 1) * x_num, y * y_num:(y + 1) * y_num, 2])
            grids.append([R,G,B])
    grids = np.array(grids, 'f')    # 变为一溜数组

    centroids, _ = vq.kmeans(grids, clus)   # (数据,分几类)->中心点
    lable, dist = vq.vq(grids, centroids)   # 为原数据根据中心点标签
    lable_im_a = lable.reshape(gridsize, gridsize)
    lable_im_a = imresize(lable_im_a, im_a.shape[0:2], interp='nearest')#(矩阵,大小)

    plt.close('all')    #一个习惯
    plt.figure(1)
    plt.imshow(lable_im_a)  # 从图片矩阵中显示
    plt.show()


B 图片求导

def _2deriv_im():
    sigma = 2
    im_a = np.array(Image.open('big.jpg').convert('L'))
    im_x = np.zeros(im_a.shape)
    filters.gaussian_filter(im_a, (sigma,sigma), (1,1), im_x)
    # (1,1):    order : int or sequence of ints, optional
    #         The order of the filter along each axis is given as a sequence
    #         of integers, or as a single number.  An order of 0 corresponds
    #         to convolution with a Gaussian kernel. A positive order
    #         corresponds to convolution with that derivative of a Gaussian.
    # 可以试试(0,1)或(1,0),效果会不同
    # filters.gaussian_filter有两种等价写法:output可以放在参数里也可放在
    # 等号左边函数的返回值。注意要求input与output不重名
    plt.figure()
    plt.imshow(im_x)
    plt.show()


C ROF去噪算法

def _3rof_denoise(image, tolerance=0.1, step=0.125, tv_weight=100, U_init=None):
    '''
    rof去噪算法,目标为min{ ||I-U||**2 + 2*weight* J(U)}, J(U)为全变差
    :param image: 灰度图像矩阵
    :param tolerance: 停止循环阈值
    :param step:
    :param tv_weight:正则项J的权值
    :param U_init:如果提供初始值。否则初始值就是原矩阵
    :return:去噪后新图像和残余
    '''
    m,n = image.shape
    U = U_init if U_init is not None else image
    Px = image  #x分量
    Py = image
    error = 1

    while(error > tolerance):   # 梯度下降,直到小于阈值,实现min{}
        Uold = U
        gradUx = np.roll(U, -1, axis=1) - U # 原始变量的x梯度
        # ~这个式子计算了U的宽度向差分差分
        gradUy = np.roll(U, -1, axis=0) - U # roll:在某一或所有坐标轴上“滚动”,使得整体移位。
        #     >>> x2
        #     array([[0, 1, 2, 3, 4],
        #            [5, 6, 7, 8, 9]])
        #     >>> np.roll(x2, 1)
        #     array([[9, 0, 1, 2, 3],
        #            [4, 5, 6, 7, 8]])
        #     >>> np.roll(x2, 1, axis=0)
        #     array([[5, 6, 7, 8, 9],
        #            [0, 1, 2, 3, 4]])
        #     >>> np.roll(x2, 1, axis=1)
        #     array([[4, 0, 1, 2, 3],
        #            [9, 5, 6, 7, 8]])

        PxNew = Px + (step / tv_weight) * gradUx
        PyNew = Py + (step / tv_weight) * gradUy
        normNew = np.maximum(1,np.sqrt(PxNew**2 + PyNew**2))

        Px = PxNew / normNew    # 更新x分量
        Py = PyNew / normNew

        rPx = np.roll(Px, 1, axis=1)    # 更新原始x分量,坐标轴平移
        rPy = np.roll(Py, 1, axis=0)

        divP = (Px - rPx) + (Py - rPy)    #对偶域散度
        U = image + tv_weight * divP    # 更新输出变量
        error = np.linalg.norm(U - Uold) / np.sqrt(m*n)    # norm()衡量矩阵间差异

    return U, image-U


D 图像角点值计算与角点数统计

def _4compute_harris_response(ima, sig=3):
    '''
    计算一张图像(矩阵形式)的harris响应值
    :param ima:
    :param sig: 越大最后得到值越小
    :return:
    '''
    imx = np.zeros(im_a.shape)
    imy = np.zeros(im_a.shape)
    filters.gaussian_filter(im_a, (sig,sig), (0,1), imx)
    filters.gaussian_filter(im_a, (sig,sig), (1,0), imy)      # 计算y方向高斯导数

    Wxx = filters.gaussian_filter(imx*imx, sig)     # harris矩阵的各分量
    Wyy = filters.gaussian_filter(imy*imy, sig)
    Wxy = filters.gaussian_filter(imx*imy, sig)

    Wdet = Wxx * Wyy - Wxy**2       # 矩阵乘法
    Wtr = Wxx + Wyy     # 特征值和迹
    Wtr = Wtr+ 0.01        # 发现对大块相同的区域,Wxx Wyy会出现导数为0,Wtr中也有0,后面除法做不下去了。
                            # 万般无奈,试试+0.1吧……验证了一下,影响不大
    ans = Wdet / Wtr    # 矩阵除法
    return ans


def _5get_harris_points(harris_ima, min_dist=10, threshold=0.1):
    '''
    输入已经计算过harris响应值的图像矩阵,输出相应的角点
    :param harris_ima: 以harris响应值为值的图像矩阵
    :param min_dist: 分割角点和图像边界的最小像素数,角点的间隔范围的控制
    :param threshold: 候选点阈值,越小生成的角点越多
    :return: 角点(x,y)坐标列表
    '''
    corner_thre = harris_ima.max() * threshold
    harris_ima_ = (harris_ima > corner_thre) * 1  # 寻找高于阈值的候选点并0-1化

    coords = np.array(harris_ima_.nonzero()).T  # 得到这些点的坐标,coords第一维为点集序号,每个点的第0维存x、1为y。 @nonzero()
    candidate_val = [harris_ima[c[0], c[1]] for c in coords] # 得到这些点的角点值和坐标
    indexs = np.argsort(candidate_val)  # 按harris响应值排序,得到序列 @argsort

    allowed_loca = np.zeros(harris_ima.shape)
    allowed_loca[ min_dist:-min_dist, min_dist:-min_dist] = 1   # 标记可行点

    selected = []
    for i in indexs:
        if allowed_loca[ coords[i,0], coords[i,1] ] == 1:
            selected.append( (coords[i,1], coords[i,0]) )
            allowed_loca[ (coords[i,0] - min_dist):(coords[i,0] + min_dist),
                            (coords[i, 1] - min_dist):(coords[i, 1] + min_dist) ] = 0   #已定点周围标记为0
    return selected