import cv2
import numpy as np
import matplotlib.pyplot as plt

def otsu_threshold_builtin(image_path):
    """
    使用OpenCV内置函数实现Otsu阈值
    :param image_path: 输入图像路径
    :return: 原图和Otsu二值化后的图像
    """
    # 读取图像
    img = cv2.imread(image_path)
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    
    # 应用Otsu阈值
    _, binary = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    
    # 获取计算出的阈值
    threshold_value = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)[0]
    
    return img, cv2.cvtColor(binary, cv2.COLOR_GRAY2BGR), threshold_value

def manual_otsu_threshold(image_path):
    """
    手动实现Otsu阈值算法
    :param image_path: 输入图像路径
    :return: 二值化后的图像和计算出的阈值
    """
    # 读取图像
    img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
    
    # 计算直方图
    hist, _ = np.histogram(img.flatten(), 256, [0, 256])
    
    # 归一化直方图
    hist_norm = hist.astype(float) / hist.sum()
    
    # 初始化变量
    max_variance = 0
    optimal_threshold = 0
    
    # 遰度级总数
    total_levels = len(hist_norm)
    
    # 计算图像的总体平均灰度
    total_mean = sum(i * hist_norm[i] for i in range(total_levels))
    
    # 初始化背景概率和平均灰度
    omega_0 = 0  # 背景像素比例
    mu_0 = 0     # 背景平均灰度
    
    # 遰度级遍历
    for t in range(total_levels):
        # 更新背景概率
        omega_0 += hist_norm[t]
        
        # 更新背景平均灰度
        if omega_0 > 0:
            mu_0 += t * hist_norm[t]
            mu_0_t = mu_0 / omega_0 if omega_0 > 0 else 0
        else:
            mu_0_t = 0
        
        # 计算前景概率和平均灰度
        omega_1 = 1.0 - omega_0
        if omega_1 > 0:
            mu_1_t = (total_mean - omega_0 * mu_0_t) / omega_1
        else:
            mu_1_t = 0
        
        # 计算类间方差
        between_class_variance = omega_0 * omega_1 * (mu_0_t - mu_1_t) ** 2
        
        # 更新最优阈值
        if between_class_variance > max_variance:
            max_variance = between_class_variance
            optimal_threshold = t
    
    # 应用计算出的阈值进行二值化
    binary = np.zeros_like(img)
    binary[img > optimal_threshold] = 255
    binary[img <= optimal_threshold] = 0
    
    return cv2.cvtColor(img, cv2.COLOR_GRAY2BGR), cv2.cvtColor(binary, cv2.COLOR_GRAY2BGR), optimal_threshold

def otsu_with_histogram(image_path):
    """
    显示Otsu算法的直方图和阈值选择过程
    :param image_path: 输入图像路径
    :return: 原图、二值化结果和相关信息
    """
    # 读取图像
    img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
    
    # 计算直方图
    hist, bins = np.histogram(img.flatten(), 256, [0, 256])
    
    # 归一化直方图
    hist_norm = hist.astype(float) / hist.sum()
    
    # 计算Otsu阈值
    total_mean = sum(i * hist_norm[i] for i in range(256))
    omega_0 = 0
    mu_0 = 0
    max_variance = 0
    optimal_threshold = 0
    
    variances = []  # 存储每个阈值的类间方差
    
    for t in range(256):
        omega_0 += hist_norm[t]
        if omega_0 > 0:
            mu_0 += t * hist_norm[t]
            mu_0_t = mu_0 / omega_0 if omega_0 > 0 else 0
        else:
            mu_0_t = 0
        
        omega_1 = 1.0 - omega_0
        if omega_1 > 0:
            mu_1_t = (total_mean - omega_0 * mu_0_t) / omega_1
        else:
            mu_1_t = 0
        
        between_class_variance = omega_0 * omega_1 * (mu_0_t - mu_1_t) ** 2
        variances.append(between_class_variance)
        
        if between_class_variance > max_variance:
            max_variance = between_class_variance
            optimal_threshold = t
    
    # 应用阈值
    binary = np.zeros_like(img)
    binary[img > optimal_threshold] = 255
    binary[img <= optimal_threshold] = 0
    
    return cv2.cvtColor(img, cv2.COLOR_GRAY2BGR), cv2.cvtColor(binary, cv2.COLOR_GRAY2BGR), optimal_threshold, variances

def multi_level_otsu(image_path, levels=3):
    """
    多级Otsu阈值（将图像分为多个等级）
    :param image_path: 输入图像路径
    :param levels: 分割等级数
    :return: 原图和多级分割结果
    """
    # 读取图像
    img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
    
    # 使用OpenCV的多级Otsu阈值
    if levels == 2:
        _, result = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    elif levels == 3:
        _, result = cv2.threshold(img, 0, 255, cv2.THRESH_TRIANGLE)
    else:
        # 对于更多等级，我们可以使用分层方法
        # 这里简化处理，使用均匀分割
        step = 255 // levels
        result = (img // step) * step
    
    return cv2.cvtColor(img, cv2.COLOR_GRAY2BGR), cv2.cvtColor(result, cv2.COLOR_GRAY2BGR)

def compare_thresholding_methods(image_path):
    """
    比较不同阈值方法的效果
    :param image_path: 输入图像路径
    :return: 原图和不同方法的二值化结果
    """
    # 读取图像
    img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
    
    # 不同的阈值方法
    methods = {
        'Global (127)': cv2.threshold(img, 127, 255, cv2.THRESH_BINARY)[1],
        'Otsu': cv2.threshold(img, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)[1],
        'Adaptive Mean': cv2.adaptiveThreshold(img, 255, cv2.ADAPTIVE_THRESH_MEAN_C, cv2.THRESH_BINARY, 11, 2),
        'Adaptive Gaussian': cv2.adaptiveThreshold(img, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, 11, 2),
        'Triangle': cv2.threshold(img, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_TRIANGLE)[1]
    }
    
    results = {}
    for name, result in methods.items():
        results[name] = cv2.cvtColor(result, cv2.COLOR_GRAY2BGR)
    
    return cv2.cvtColor(img, cv2.COLOR_GRAY2BGR), results

def otsu_performance_analysis(image_path):
    """
    Otsu算法性能分析
    :param image_path: 输入图像路径
    :return: 分析结果
    """
    # 读取图像
    img = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
    
    # 计算直方图
    hist, _ = np.histogram(img.flatten(), 256, [0, 256])
    hist_norm = hist.astype(float) / hist.sum()
    
    # 计算Otsu阈值
    total_mean = sum(i * hist_norm[i] for i in range(256))
    omega_0 = 0
    mu_0 = 0
    max_variance = 0
    optimal_threshold = 0
    
    for t in range(256):
        omega_0 += hist_norm[t]
        if omega_0 > 0:
            mu_0 += t * hist_norm[t]
            mu_0_t = mu_0 / omega_0 if omega_0 > 0 else 0
        else:
            mu_0_t = 0
        
        omega_1 = 1.0 - omega_0
        if omega_1 > 0:
            mu_1_t = (total_mean - omega_0 * mu_0_t) / omega_1
        else:
            mu_1_t = 0
        
        between_class_variance = omega_0 * omega_1 * (mu_0_t - mu_1_t) ** 2
        
        if between_class_variance > max_variance:
            max_variance = between_class_variance
            optimal_threshold = t
    
    # 计算分割质量指标
    foreground_pixels = np.sum(img > optimal_threshold)
    background_pixels = np.sum(img <= optimal_threshold)
    total_pixels = img.size
    
    foreground_ratio = foreground_pixels / total_pixels
    background_ratio = background_pixels / total_pixels
    
    # 应用阈值
    binary = np.zeros_like(img)
    binary[img > optimal_threshold] = 255
    binary[img <= optimal_threshold] = 0
    
    return {
        'optimal_threshold': optimal_threshold,
        'max_between_class_variance': max_variance,
        'foreground_ratio': foreground_ratio,
        'background_ratio': background_ratio,
        'original_image': cv2.cvtColor(img, cv2.COLOR_GRAY2BGR),
        'binary_image': cv2.cvtColor(binary, cv2.COLOR_GRAY2BGR)
    }

# 使用示例
if __name__ == "__main__":
    # 注意：需要提供实际的图像路径
    # img, result, threshold = otsu_threshold_builtin('image.jpg')
    # manual_img, manual_result, manual_threshold = manual_otsu_threshold('image.jpg')
    # hist_img, hist_result, hist_threshold, variances = otsu_with_histogram('image.jpg')
    # multi_img, multi_result = multi_level_otsu('image.jpg', 3)
    # comp_img, comp_results = compare_thresholding_methods('image.jpg')
    # analysis_result = otsu_performance_analysis('image.jpg')
    pass