opencv第五次作业（边界线链码，灰度共生矩阵，角点特征，sift）

在这里插入代码片
```#!/usr/bin/env python
# coding=utf-8

import cv2
import numpy as np
import easygui as g
import math
np.set_printoptions(suppress=True)


def func(x,y,sigma=1):
    return 100*(1/(2*np.pi*sigma))*np.exp(-((x-1)**2+(y-1)**2)/(2.0*sigma**2))   #创建高斯函数，该函数中心为（1,1），所以创建3*3比较合适

def convolution(kernal_size, image):
    rows = image.shape[0]
    cols = image.shape[1]
    image_new = np.zeros((rows+2,cols+2))
    image_new[1:rows+1,1:cols+1] = image.copy()
    # gray_data = image.copy()
    img_new = []
    for i in range(1,rows+1):   #减去3的原因是自己定义的模板都是3*3
        line = []             #line记录着每一行的信息
        for j in range(1,cols+1):
            a = image_new[i-1:i+2,j-1:j+2]
            a=np.array(a)
            line.append(np.sum(np.multiply(kernal_size, a)))   #np.multiply表示两个矩阵点乘
        img_new.append(line)   #记录着每一行卷积后的结果
    return np.array(img_new)

def non_maximum_suppression(hough):
    r,l  = hough.shape

    ## non maximum suppression
    for y in range(r):
        for x in range(l):
            if x-1>=0 and x+2 <= l and y-1 >=0 and y+2<=r:
                # get 8 nearest neighbor
                x1 = x-1
                x2 = x+2
                y1 = y-1
                y2 = y+2
                if np.max(hough[y1:y2, x1:x2]) == hough[y, x] and hough[y, x] != 0:
                    pass
                    # hough[y,x] = 255
                else:
                    hough[y, x] = 0

    return hough
def grad_x(h):
    delta_h = h.copy()
    a = int(h.shape[0])
    b = int(h.shape[1])

    for i in range(a):
        for j in range(b):
            if i-1>=0 and i+1<a and j-1>=0 and j+1<b:           #防止超出图像边界
                c = abs(int(h[i-1,j-1]) - int(h[i+1,j-1]) + 2*(int(h[i-1,j]) - int(h[i+1,j])) + int(h[i-1,j+1]) - int(h[i+1,j+1]))
#                print c
                if c>255:
#                    print c
                    c = 255         #限制我们的最大灰度值
                delta_h[i,j] = c
            else:
                delta_h[i,j] = 0
    print ('x方向的梯度：\n %s \n' %delta_h)
    return delta_h

##计算y方向的梯度的函数（水平方向锐化算子）
def grad_y(h):
    a = int(h.shape[0])
    b = int(h.shape[1])
    delta_h = h.copy()

    for i in range(a):
        for j in range(b):
            if i-1>=0 and i+1<a and j-1>=0 and j+1<b:
                c = abs(int(h[i-1,j-1]) - int(h[i-1,j+1]) + 2*(int(h[i,j-1]) - int(h[i,j+1])) + (int(h[i+1,j-1]) - int(h[i+1,j+1])))  #注意像素不能直接计算，需要转化为整型
#                print c
                if c > 255:
                    c = 255
                delta_h[i,j] = c
            else:
                delta_h[i,j] = 0
    print ('y方向的梯度：\n %s \n' %delta_h)
    return delta_h

def laplace(image):
    laplace = np.array([[1, 1, 1],
                        [1, -8, 1],
                        [1, 1, 1]])

    image_laplace = convolution(image=image, kernal_size=laplace)  # 进行高斯平滑后在进行我们的拉place变换
    image_laplace = image_laplace * (255 / image_laplace.max())
    image_laplace = convolution(image=image_laplace, kernal_size=laplace)
    return image_laplace

def glcm(arr, d_x, d_y, gray_level=16):
    '''计算并返回归一化后的灰度共生矩阵'''
    max_gray = arr.max()
    height, width = arr.shape
    arr = arr.astype(np.float64)  # 将uint8类型转换为float64，以免数据失真
    arr = arr * (gray_level - 1) // max_gray  # 若灰度级数大于gray_level，则将图像的灰度级缩小至gray_level，减小灰度共生矩阵的大小。量化后灰度值范围：0 ~ gray_level - 1，整张图片的灰度值大小缩小至0到16之间
    ret = np.zeros([gray_level, gray_level])    #ret矩阵大小为 16*16
    for j in range(height - abs(d_y)):
        for i in range(width - abs(d_x)):
            rows = arr[j][i].astype(int)                    #arr经过缩减尺度后依旧为np.float类型，有可能还是小数形式，我们这里将其转化为整数形式，这里rows指代的是ret映射表的x
            cols = arr[j + d_y][i + d_x].astype(int)        #arr经过缩减尺度后依旧为np.float类型，有可能还是小数形式，我们这里将其转化为整数形式 ,这里cols指代的是ret映射表的y
            ret[rows][cols] += 1
    # if d_x >= d_y:
    #     ret = ret / float(height * (width - 1))  # 归一化, 水平方向或垂直方向
    # else:
    #     ret = ret / float((height - 1) * (width - 1))  # 归一化, 45度或135度方向
    # rols_ret = ret.shape[0]
    # cols_ret = ret.shape[1]
    ret = ret /(height * width)
    return ret

def gause_image(image):             #高斯滤波，第一步也实现了灰度化
    suanzi = np.fromfunction(func,(3,3),sigma=2)  # 创建高斯函数，（1,1）为函数的中心，这里是生成3*3的标准差为2的高斯算子
    # 打开图像并转化成灰度图像
    # image = image[:, :, 0] * 0.11 + image[:, :, 1] * 0.59 + image[:, :, 2] * 0.3
    # image = image.astype(np.uint8)  # GRAY=0.3*R+0.59*G+0.11*B：
    image2 = convolution(image=image, kernal_size=suanzi)
    # 将结果灰度值转化到0-255
    image2 = (image2 / float(image2.max())) * 255
    return image2

def duibidu_number(ret):
    rows = ret.shape[0]
    cols  = ret.shape[1]
    sum = 0
    for i in range(rows):
        for j in range(cols):
           sum += (i-j)**2*ret[i][j]
    return sum
def max_number(ret):
    rows = ret.shape[0]
    cols = ret.shape[1]
    sum = 0
    for i in range(rows):
        for j in range(cols):
            if sum <= ret[i][j]:
                sum = ret[i][j]
    return sum
def er_jie_ju(ret):
    rows = ret.shape[0]
    cols = ret.shape[1]
    sum = 0
    for i in range(rows):
        for j in range(cols):
            sum += ret[i][j]**2
    return sum
def fancha_fen_matrix(ret):
    rows = ret.shape[0]
    cols = ret.shape[1]
    sum = 0
    for i in range(rows):
        for j in range(cols):
            sum +=(ret[i][j])/(1+(i-j)**2)
    return sum

def fun1(ret):
    dui_bi_du = duibidu_number(ret)
    max__number = max_number(ret)
    er_jieju_number = er_jie_ju(ret)
    fan=fancha_fen_matrix(ret)
    print('对比度为%f' % (dui_bi_du),end='\n')
    print('最大概率为 %f' % (max__number),end='\n')
    print('二阶矩 %f ' % (er_jieju_number),end='\n')
    print('反差分为%f' % (fan))


def line_length(dict):
    sum = 0
    for i in range(0,len(dict),2):
        x1 = dict[i][0]
        y1 = dict[i][1]
        if i+1 <= len(dict)-1:
            x2 = dict[i+1][0]
            y2 = dict[i+1][1]
            sum += math.sqrt((x1-x2)**2+(y1 -y2)**2)
    sum += math.sqrt((dict[0][0]-dict[len(dict)-1][0])**2+(dict[0][1]-dict[len(dict)-1][1])**2)
    return sum
def area_lianma(dict,img):                              #基于三角型的目标面积计算方法
    area = 0
    for i in range(0,len(dict),2):
        x1 = dict[i][0]
        y1 = dict[i][1]
        if i + 1 <= len(dict) - 1:
            x2 = dict[i + 1][0]
            y2 = dict[i + 1][1]
            area += (x2 * y1 - x1 * y2)/2               #面积求导公式
    rols = img.shape[0]
    cols = img.shape[1]
    area_all = rols * cols
    gray_point = 0
    white_point = 0
    for i in range(rols):
        for j in range(cols):
            if img[i,j] == 255:
                white_point +=1
            else:
                img[i,j] == 0
                gray_point +=1
    # print('利用像素点算出来的白色面积为 %d ' % white_point)
    #
    # print('总面积为%d' % (area_all))
    # print('外界的面积 %d ' % (area))
    return area_all + area


def lianma(img):
# img = cv2.imread("E:\\360Downloads\\lianma.jpg")
    [img_h, img_w, img_channel] = img.shape
    trace = []
    start_x = 0
    start_y = 0

    img1 = img[:,:,0] * 0.11 + img[:,:,1] * 0.59 + img[:,:,2] * 0.3
    img1 = img1.astype(np.uint8)
    gray = img1.copy()
    img3 = img[:,:,0] * 0.11 + img[:,:,1] * 0.59 + img[:,:,2] * 0.3
    img3 = img3.astype(np.uint8)
    # for h in range(img_h):
    #     for w in range(img_w):
    #         if (gray[h, w] > 128):
    #             gray[h, w] = 255
    #         else:
    #             gray[h, w] = 0
    gray[gray>128] = 255
    gray[gray<128] =0


    # python 跳出多重循环

    class getoutofloop(Exception): pass


    try:
        for h in range(10,img_h - 2):
            for w in range(10,img_w - 2):
                if gray[h, w] == 255:                   #先根据灰度值为255找到我们的一个起始点
                    start_x = w
                    start_y = h
                    raise getoutofloop
    except getoutofloop:
        pass
    # start_x = 110
    # start_y = 15

    print("Start Point (%d %d)" % (start_x, start_y))
    trace.append([start_x, start_y])

    # 8邻域 顺时针方向搜索
    neighbor = [ [ 1, 1 ], [ 1, 0 ], [1,-1],[0,-1],[-1,-1],[-1,0],[-1,1],[0,1]]  #从右下角45度的地方开始按照顺时针的方向来依次遍历
    # neighbor = [[0,1],[],[],[],[],[],[],[],[]]
    neighbor_len = len(neighbor)

    # 先从当前点的左上方开始，
    # 如果左上方也是黑点(边界点)：
    #         搜索方向逆时针旋转90 i-=2
    # 否则：
    #         搜索方向顺时针旋转45 i+=1
    i = 0
    cur_x = start_x + neighbor[i][0]                #当前点的坐标位置 = 起始点 + 邻域坐标
    cur_y = start_y + neighbor[i][1]

    is_contour_point = 0

    try:
        while not ((cur_x == start_x) and (cur_y == start_y)):                  #如果当前点的坐标与起始点坐标相同就跳出循环
            is_contour_point = 0
            while is_contour_point == 0:
                # neighbor_x = cur_x +
                if gray[cur_y, cur_x] == 0:                 #为黑色边界点（像素值为0）才让trace列表添加进去
                    is_contour_point = 1
                    trace.append([cur_x, cur_y])            #trace列表添加的是（x,y）
                    i -= 2                                  #像素值为0时逆时针旋转90度（i-2,i代表着我们8邻域当前序号）
                    if i < 0:
                        i += neighbor_len
                else:
                    i += 1                                  #像素值为255时顺时针旋转45度（i+1）
                    if i >= neighbor_len:
                        i -= neighbor_len
                # print(i)
                cur_x = cur_x + neighbor[i][0]
                cur_y = cur_y + neighbor[i][1]
    except:
        print("throw error")


    img3 = cv2.cvtColor(img3, cv2.COLOR_GRAY2RGB)
    for i in range(len(trace) - 1):
        cv2.line(img3, (trace[i][0], trace[i][1]), (trace[i + 1][0], trace[i + 1][1]), (0, 0, 255), 3)           #根据两点画一条线
        cv2.imshow("img", img3)
        cv2.waitKey(10)
    image2 = np.zeros((img.shape[0],img.shape[1]),np.uint8)
    image2 = cv2.cvtColor(image2,cv2.COLOR_GRAY2RGB)            #很重要，转化为三通道的图像
    for i in range(len(trace)):
        # if trace[i][0] <= image2.shape[0] and trace[i][1] <= image2.shape[1]:
        image2[trace[i][1],trace[i][0]] = [0,0,255]     #注意是image[y,x]的形式
    cv2.imshow('dot image',image2)
    cv2.rectangle(img3, (start_x-10, start_y-10), (start_x + 10, start_y + 10) , (0, 255, 0), 2)
    cv2.imshow("img", img)
    cv2.imshow('img2',img3)
    line = line_length(trace)
    g.msgbox('线的长度为 %d' % (line),title='线的长度')
    # print('length',line)
    lianma_area = area_lianma(trace,gray)
    g.msgbox('链码包围的面积为 %d' % (lianma_area),title='链码面积' )           #利用三角型的目标面积法来进行计算
    # print('area',lianma_area)
    cv2.waitKey(0)
    cv2.destroyWindow("img")

def harris(img):
    img2 = img.copy()
    img_gray = img2[:, :, 0] * 0.11 + img2[:, :, 1] * 0.59 + img2[:, :, 2] * 0.3
    img_gray = img_gray.astype(np.uint8)

    # 转换为灰度图
    gray = img_gray.copy()
    m = np.array(gray)
    # 计算x方向的梯度的函数（水平方向锐化算子）
    # delta_h = m
    img_laplace = laplace(gray)  # 先进行laplace变换

    dx = np.array(grad_x(gray))  # 再根据梯度函数来获取我们需要的x偏导与y偏导
    dy = np.array(grad_y(gray))

    A = dx * dx  # 获得Ix*Ix
    B = dy * dy  # 获得Iy * Iy
    C = dx * dy  # 获得Ixy

    A1 = A
    B1 = B
    C1 = C

    a = int(gray.shape[0])
    b = int(gray.shape[1])

    A1 = gause_image(A1)  # 高斯平滑处理
    B1 = gause_image(B1)
    C1 = gause_image(C1)

    a_rols = A1.shape[0]
    a_cols = A1.shape[1]

    R = np.zeros(gray.shape)                        #构造我们自己的 R = (AC-B^2) - α（A+C）^2
    for i in range(a):
        for j in range(b):
            M = [[A1[i, j], C1[i, j]], [C1[i, j], B1[i, j]]]
            R[i, j] = np.linalg.det(M) - 0.06 * (np.trace(M)) * (np.trace(M))       #这里的阿尔法我们设为0.06
    threshold = 2500
    R[R < threshold] = 0            #小于0 的R设为0

    R = non_maximum_suppression(R)  # 进行非极大值抑制

    print(R)
    rols = R.shape[0]
    cols = R.shape[1]
    # print(rols,cols)
    pix_number = []
    for i in range(rols):
        for j in range(cols):
            if R[i, j] >= threshold:
                pix_number.append((i, j))           #记录非零值的R函数的当前坐标，并用列表存储起来

    for i in (pix_number):
        # print(pix_number[i])
        cv2.circle(img2, (i[1], i[0]), 2, (0, 0, 255), 1)
    # gray = cv2.resize(gray,(rols,cols))
    cv2.imshow('harris_detection', img2)
    cv2.namedWindow('R', cv2.WINDOW_NORMAL)
    cv2.imshow('R', R)

    cv2.waitKey(0)
    cv2.destroyAllWindows()


if __name__ == '__main__':
    msg = "请输入您想要完成的任务（建议您第一步先打开图片）"
    title = '第五次作业'
    choice = ('打开图片', '退出')
    a = g.buttonbox(msg=msg, title=title, choices=choice)
    if a == '打开图片':
        filename = g.fileopenbox(msg="请打开一个jpg文件")
        img = cv2.imread(filename)
        msg1 = "选择您想要实现的功能"
        title1 = '第五次作业'
        choice1 = ('边界线链码', '灰度共生矩阵', '角点特征','SIFT','显示原图','重新选择图片','退出')
        q = 1
        while q:
            b = g.buttonbox(msg=msg1, title=title1, choices=choice1)
            if b == '边界线链码':
               lianma(img)
            elif b == '灰度共生矩阵':
                img_gray = img[:,:,0] * 0.11+img[:,:,1] *0.59 + img[:,:,2] *0.3
                img_gray = img_gray.astype(np.uint8)
                glcm_0 = glcm(img_gray, 1, 0)  # 水平方向
                fun1(glcm_0)
                print(glcm_0,
                      end='\n\n**************************************************************************************\n\n\n')
                glcm_1 = glcm(img_gray, 0, 1)  # 垂直方向
                fun1(glcm_1)
                print(glcm_1,
                      end='\n\n**************************************************************************************\n\n\n')
                glcm_2 = glcm(img_gray, 1, 1)  # 45度方向
                fun1(glcm_2)
                print(glcm_2,
                      end='\n\n**************************************************************************************\n\n\n')
                glcm_3 = glcm(img_gray, -1, -1)  # 135度方向
                fun1(glcm_3)
                print(glcm_3,
                      end='\n\n**************************************************************************************\n\n\n')

            elif b == '角点特征':
               harris(img)


            elif b == 'SIFT':
                pass
            elif b == '显示原图':
                cv2.imshow('原图',img)
            elif b == '重新选择图片':
                filename = g.fileopenbox(msg="请打开一个jpg文件")
                img = cv2.imread(filename)
            else:
                q = 0