import cv2
import numpy as np
from matplotlib import pyplot as plt
import copy

img = cv2.imread(r'.\lena.jpg')
image = cv2.imread(r'.\rice.jpg')
img_chess = cv2.imread(r'.\chessboard.png')

cv2.namedWindow("lena", 0)
cv2.resizeWindow("lena", 640, 480)
cv2.imshow("lena", img)

# 题目一，均值滤波、高斯滤波、中值滤波
mean_img = cv2.blur(img, (101, 101))
cv2.namedWindow("mean", 0)
cv2.resizeWindow("mean", 640, 480)
cv2.imshow("mean", mean_img)

gaussian_img = cv2.GaussianBlur(img, (101, 101), 0)
cv2.namedWindow("gaussian", 0)
cv2.resizeWindow("gaussian", 640, 480)
cv2.imshow("gaussian", gaussian_img)

median_img = cv2.medianBlur(img, 101)
cv2.namedWindow("median", 0)
cv2.resizeWindow("median", 640, 480)
cv2.imshow("median", median_img)

# Sobel算子、Canny算子检测

# Sobel
x = cv2.Sobel(img, cv2.CV_16S, 1, 0)
y = cv2.Sobel(img, cv2.CV_16S, 0, 1)
absX = cv2.convertScaleAbs(x)  # 转回uint8
absY = cv2.convertScaleAbs(y)

sobel_img = cv2.addWeighted(absX, 0.5, absY, 0.5, 0)

cv2.namedWindow("sobel", 0)
cv2.resizeWindow("sobel", 640, 480)
cv2.imshow("sobel", sobel_img)

# Canny
canny_img = cv2.Canny(img, 100, 100)
cv2.namedWindow("canny", 0)
cv2.resizeWindow("canny", 640, 480)
cv2.imshow("canny", canny_img)

# 题目三，计算灰度直方图，使用大津算法进行分割
cv2.namedWindow("rice", 0)
cv2.resizeWindow("rice", 640, 480)
cv2.imshow("rice", image)

hist = cv2.calcHist([image], [0], None, [256], [0, 255])
plt.plot(hist)
plt.show()

gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
temp = cv2.threshold(gray, 0, 0xff, cv2.THRESH_OTSU)
otsu = temp[1]
cv2.namedWindow("otsu", 0)
cv2.resizeWindow("otsu", 640, 480)
cv2.imshow("otsu", otsu)

# 题目四，使用米粒图像，分割得到各米粒，首先计算各区域(米粒)的面积、长度等信息，进一步计算面积、长度的均值及方差，分析落在3sigma范围内米粒的数量
element = cv2.getStructuringElement(cv2.MORPH_CROSS, (3, 3))
otsu = cv2.morphologyEx(otsu, cv2.MORPH_OPEN, element)

seg = copy.deepcopy(otsu)
bin, cnts, hier = cv2.findContours(seg, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
count = 0
area_array = np.array([])  # 米粒面积数组
line_array = np.array([])  # 米粒长度数组
for i in range(len(cnts), 0, -1):
    c = cnts[i-1]
    area = cv2.contourArea(c)
    if area < 10:
        continue
    count = count + 1
    area_array = np.append(area_array, area)

    minAreaRect = cv2.minAreaRect(c)  # 精确米粒边界
    x, y, z, p = cv2.boxPoints(minAreaRect)  # 得到米粒的4个坐标
    # 计算米粒的长度
    l1 = ((x[0] - y[0]) ** 2 + (x[1] - y[1]) ** 2) ** 0.5  # 计算米粒边界长度
    l2 = ((x[0] - p[0]) ** 2 + (x[1] - p[1]) ** 2) ** 0.5
    if l1 > l2:
        line = l1
    else:
        line = l2
    print("米粒", i, "面积:", area, "长度:", round(line, 2))
    line_array = np.append(line_array, line)
    cv2.line(image, tuple(x), tuple(y), (255, 0, 0), 1)
    cv2.line(image, tuple(y), tuple(z), (255, 0, 0), 1)
    cv2.line(image, tuple(z), tuple(p), (255, 0, 0), 1)
    cv2.line(image, tuple(p), tuple(x), (255, 0, 0), 1)
    cv2.putText(image, str(count), (x[0], x[1]), cv2.FONT_HERSHEY_PLAIN, 0.5, (0, 0xff, 0))

print("------------------------------")

print("米粒总数量: ", len(area_array))
print("------------------------------")

area_all = area_array.sum()
print("总面积:", round(area_all, 3))
average_area = area_all / count
print("平均面积:", round(average_area, 3))
std = area_array.std()
print("面积标准差:", round(std, 3))
area1 = average_area - std * 1.5
area2 = average_area + std * 1.5
print("面积3sigma的取值范围:", round(area1, 3), "--", round(area2, 3))
count1 = 0
for i in area_array:
    if area1 < i < area2:
        count1 = count1 + 1
print("米粒面积在3sigma内的数量为:", count1)

print("------------------------------")

line_all = line_array.sum()
print("总长度:", round(line_all, 3))
average_line = line_all / count
print("平均长度:", round(average_line, 3))
std_line = line_array.std()
print("长度标准差:", round(std_line, 3))
line1 = average_line - std_line * 1.5
line2 = average_line + std_line * 1.5
print("长度3sigma的取值范围:", round(line1, 3), "--", round(line2, 3))
count2 = 0
for i in line_array:
    if line2 > i > line1:
        count2 = count2 + 1
print("米粒长度在3sigma内的数量为:", count2)
cv2.imshow("source", image)
cv2.imshow("otsu2", otsu)

# 题目五，使用棋盘格及自选风景图像，分别使用SIFT、FAST及ORB算子检测角点;使用Harris角点检测算子检测棋盘格
gray = cv2.cvtColor(img_chess, cv2.COLOR_BGR2GRAY)
cv2.namedWindow("origin", 0)
cv2.resizeWindow("origin", 640, 480)
cv2.imshow('origin', gray)

# SIFT
image_sift = img_chess.copy()
detector = cv2.xfeatures2d.SIFT_create()
keypoints = detector.detect(gray, None)
cv2.drawKeypoints(gray, keypoints, image_sift, (255, 0, 0), cv2.DRAW_MATCHES_FLAGS_DEFAULT)
cv2.namedWindow("SIFT", 0)
cv2.resizeWindow("SIFT", 640, 480)
cv2.imshow('SIFT', image_sift)

# FAST
image_fast = img_chess.copy()
fast = cv2.FastFeatureDetector_create(40, True, cv2.FAST_FEATURE_DETECTOR_TYPE_9_16)
keypoints1 = fast.detect(image_fast, None)
# 在图像上绘制关键点
cv2.drawKeypoints(image_fast, keypoints1, image_fast, (255, 0, 255), cv2.DRAW_MATCHES_FLAGS_DEFAULT)
print("Threshold: ", fast.getThreshold())
print("nonmaxSuppression: ", fast.getNonmaxSuppression())
print("neighborhood: ", fast.getType())
print("Total Keypoints with nonmaxSuppression: ", len(keypoints1))
cv2.namedWindow("FAST", 0)
cv2.resizeWindow("FAST", 1200, 1200)
cv2.imshow('FAST', image_fast)

# ORB
image_orb = img_chess.copy()
orb = cv2.ORB_create()
kp = orb.detect(image_orb, None)
kp, des = orb.compute(image_orb, kp)
cv2.drawKeypoints(image_orb, kp, image_orb, (0, 255, 0), cv2.DRAW_MATCHES_FLAGS_DEFAULT)
cv2.namedWindow("ORB", 0)
cv2.resizeWindow("ORB", 1200, 1200)
cv2.imshow('ORB', image_orb)

# Harris
image_harris = img_chess.copy()
gray_img = cv2.cvtColor(image_harris, cv2.COLOR_BGR2GRAY)
gray_img = np.float32(gray_img)

# 对图像执行harris
Harris_detector = cv2.cornerHarris(gray_img, 100, 19, 0.04)

# 腐蚀harris结果
dst = cv2.dilate(Harris_detector, None)

# 设置阈值
thres = 0.01 * dst.max()

image_harris[dst > thres] = [0, 0, 255]

cv2.namedWindow("Harris", 0)
cv2.resizeWindow("Harris", 1200, 1200)
cv2.imshow('Harris', image_harris)

cv2.waitKey()
cv2.destroyAllWindows()