图像分析与特征提取

1. 边缘检测与分析

Canny边缘检测详解

import cv2
import numpy as np
import matplotlib.pyplot as plt
from skimage import feature, measure, filters
from scipy import ndimage

def comprehensive_edge_detection(img):
    """全面的边缘检测对比"""
  
    # 转为灰度图
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
  
    # 不同的边缘检测方法
  
    # 1. Canny边缘检测
    canny_50_150 = cv2.Canny(gray, 50, 150)
    canny_100_200 = cv2.Canny(gray, 100, 200)
  
    # 2. Sobel边缘检测
    sobel_x = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=3)
    sobel_y = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=3)
    sobel_mag = np.sqrt(sobel_x**2 + sobel_y**2)
    sobel_mag = np.uint8(sobel_mag / sobel_mag.max() * 255)
  
    # 3. Laplacian边缘检测
    laplacian = cv2.Laplacian(gray, cv2.CV_64F)
    laplacian = np.uint8(np.absolute(laplacian))
  
    # 4. Scharr算子
    scharr_x = cv2.Scharr(gray, cv2.CV_64F, 1, 0)
    scharr_y = cv2.Scharr(gray, cv2.CV_64F, 0, 1)
    scharr_mag = np.sqrt(scharr_x**2 + scharr_y**2)
    scharr_mag = np.uint8(scharr_mag / scharr_mag.max() * 255)
  
    # 显示结果
    images = [gray, canny_50_150, canny_100_200, sobel_mag, laplacian, scharr_mag]
    titles = ['原图', 'Canny(50,150)', 'Canny(100,200)', 'Sobel', 'Laplacian', 'Scharr']
  
    fig, axes = plt.subplots(2, 3, figsize=(15, 10))
    axes = axes.ravel()
  
    for i, (image, title) in enumerate(zip(images, titles)):
        axes[i].imshow(image, cmap='gray')
        axes[i].set_title(title)
        axes[i].axis('off')
  
    plt.tight_layout()
    plt.show()
  
    return canny_50_150, sobel_mag, laplacian

def adaptive_canny(img, sigma=0.33):
    """自适应Canny边缘检测"""
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
  
    # 计算图像的中值
    v = np.median(gray)
  
    # 应用自适应阈值
    lower = int(max(0, (1.0 - sigma) * v))
    upper = int(min(255, (1.0 + sigma) * v))
  
    # Canny边缘检测
    edged = cv2.Canny(gray, lower, upper)
  
    return edged, lower, upper

# 使用示例
# edges = comprehensive_edge_detection(img)
# adaptive_edges, low, high = adaptive_canny(img)

边缘特征分析

def edge_analysis(img):
    """边缘特征分析"""
  
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    edges = cv2.Canny(gray, 50, 150)
  
    # 1. 边缘密度
    edge_density = np.sum(edges > 0) / edges.size
  
    # 2. 边缘方向直方图
    sobel_x = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=3)
    sobel_y = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=3)
  
    # 计算梯度方向
    gradient_direction = np.arctan2(sobel_y, sobel_x) * 180 / np.pi
    gradient_direction = (gradient_direction + 180) % 180  # 0-180度
  
    # 创建方向直方图
    hist_bins = 18  # 每10度一个bin
    direction_hist, _ = np.histogram(gradient_direction[edges > 0], 
                                   bins=hist_bins, range=[0, 180])
  
    # 3. 边缘长度统计
    contours, _ = cv2.findContours(edges, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    edge_lengths = [cv2.arcLength(contour, False) for contour in contours]
  
    # 可视化结果
    fig, axes = plt.subplots(2, 2, figsize=(15, 10))
  
    # 原图和边缘
    axes[0,0].imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
    axes[0,0].set_title('原图')
    axes[0,0].axis('off')
  
    axes[0,1].imshow(edges, cmap='gray')
    axes[0,1].set_title(f'边缘图 (密度: {edge_density:.3f})')
    axes[0,1].axis('off')
  
    # 方向直方图
    angles = np.arange(0, 180, 180/hist_bins)
    axes[1,0].bar(angles, direction_hist, width=8)
    axes[1,0].set_title('边缘方向直方图')
    axes[1,0].set_xlabel('角度 (度)')
    axes[1,0].set_ylabel('频次')
  
    # 边缘长度分布
    if edge_lengths:
        axes[1,1].hist(edge_lengths, bins=20, alpha=0.7)
        axes[1,1].set_title('边缘长度分布')
        axes[1,1].set_xlabel('长度')
        axes[1,1].set_ylabel('频次')
  
    plt.tight_layout()
    plt.show()
  
    return {
        'edge_density': edge_density,
        'direction_histogram': direction_hist,
        'edge_lengths': edge_lengths,
        'mean_edge_length': np.mean(edge_lengths) if edge_lengths else 0
    }

# 使用示例
# edge_features = edge_analysis(img)

2. 纹理分析

灰度共生矩阵(GLCM)

def glcm_features(img, distances=[1], angles=[0, 45, 90, 135]):
    """基于灰度共生矩阵的纹理特征提取"""
  
    # 转换为灰度图并降低灰度级别以减少计算量
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
  
    # 降低灰度级别到64级
    gray_reduced = (gray / 4).astype(np.uint8)
  
    # 计算GLCM特征（需要安装scikit-image）
    try:
        from skimage.feature import greycomatrix, greycoprops
    
        # 计算GLCM矩阵
        glcm = greycomatrix(gray_reduced, distances=distances, 
                           angles=np.radians(angles), levels=64, symmetric=True, normed=True)
    
        # 提取纹理特征
        contrast = greycoprops(glcm, 'contrast')
        dissimilarity = greycoprops(glcm, 'dissimilarity') 
        homogeneity = greycoprops(glcm, 'homogeneity')
        energy = greycoprops(glcm, 'energy')
        correlation = greycoprops(glcm, 'correlation')
    
        # 计算各个方向的平均值
        features = {
            'contrast': np.mean(contrast),
            'dissimilarity': np.mean(dissimilarity),
            'homogeneity': np.mean(homogeneity),
            'energy': np.mean(energy),
            'correlation': np.mean(correlation)
        }
    
        return features, glcm
    
    except ImportError:
        print("需要安装scikit-image库: pip install scikit-image")
        return None, None

def local_binary_pattern(img, radius=3, n_points=24):
    """局部二值模式纹理特征"""
  
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
  
    try:
        from skimage import feature
    
        # 计算LBP
        lbp = feature.local_binary_pattern(gray, n_points, radius, method='uniform')
    
        # 计算LBP直方图
        lbp_hist, _ = np.histogram(lbp.ravel(), bins=n_points + 2, 
                                  range=(0, n_points + 2), density=True)
    
        # 可视化
        fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
        axes[0].imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
        axes[0].set_title('原图')
        axes[0].axis('off')
    
        axes[1].imshow(lbp, cmap='gray')
        axes[1].set_title(f'LBP (R={radius}, P={n_points})')
        axes[1].axis('off')
    
        axes[2].bar(range(len(lbp_hist)), lbp_hist)
        axes[2].set_title('LBP直方图')
        axes[2].set_xlabel('LBP值')
        axes[2].set_ylabel('频率')
    
        plt.tight_layout()
        plt.show()
    
        return lbp, lbp_hist
    
    except ImportError:
        print("需要安装scikit-image库")
        return None, None

def gabor_filters(img):
    """Gabor滤波器纹理分析"""
  
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
  
    # 不同参数的Gabor滤波器
    gabor_responses = []
  
    # 不同方向和频率的Gabor滤波器
    for theta in [0, 45, 90, 135]:  # 方向
        for frequency in [0.1, 0.3, 0.5]:  # 频率
            # 创建Gabor滤波器
            kernel_real = cv2.getGaborKernel((21, 21), 5, np.radians(theta), 
                                           2*np.pi*frequency, 0.5, 0, ktype=cv2.CV_32F)
        
            # 应用滤波器
            response = cv2.filter2D(gray, cv2.CV_8UC3, kernel_real)
            gabor_responses.append(response)
  
    # 显示Gabor响应
    fig, axes = plt.subplots(4, 3, figsize=(12, 16))
    axes = axes.ravel()
  
    for i, response in enumerate(gabor_responses):
        axes[i].imshow(response, cmap='gray')
        theta = [0, 45, 90, 135][i // 3]
        freq = [0.1, 0.3, 0.5][i % 3]
        axes[i].set_title(f'θ={theta}°, f={freq}')
        axes[i].axis('off')
  
    plt.tight_layout()
    plt.show()
  
    return gabor_responses

# 使用示例
# glcm_feats, glcm_matrix = glcm_features(img)
# lbp_img, lbp_hist = local_binary_pattern(img)
# gabor_resp = gabor_filters(img)

3. 统计特征提取

颜色特征

def color_features(img):
    """颜色特征提取"""
  
    # 转换到不同颜色空间
    hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
    lab = cv2.cvtColor(img, cv2.COLOR_BGR2LAB)
  
    # 1. 颜色直方图
    hist_b = cv2.calcHist([img], [0], None, [256], [0, 256])
    hist_g = cv2.calcHist([img], [1], None, [256], [0, 256])
    hist_r = cv2.calcHist([img], [2], None, [256], [0, 256])
  
    # 2. 统计矩
    def calculate_moments(channel):
        mean = np.mean(channel)
        std = np.std(channel)
        skewness = np.mean(((channel - mean) / std) ** 3)
        kurtosis = np.mean(((channel - mean) / std) ** 4) - 3
        return mean, std, skewness, kurtosis
  
    # BGR通道统计矩
    b_moments = calculate_moments(img[:,:,0])
    g_moments = calculate_moments(img[:,:,1])
    r_moments = calculate_moments(img[:,:,2])
  
    # HSV通道统计矩
    h_moments = calculate_moments(hsv[:,:,0])
    s_moments = calculate_moments(hsv[:,:,1])
    v_moments = calculate_moments(hsv[:,:,2])
  
    # 3. 主要颜色提取（K-means聚类）
    data = img.reshape((-1, 3))
    data = np.float32(data)
  
    criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 20, 1.0)
    k = 5
    _, labels, centers = cv2.kmeans(data, k, None, criteria, 10, cv2.KMEANS_RANDOM_CENTERS)
    centers = np.uint8(centers)
  
    # 计算每个聚类的比例
    unique, counts = np.unique(labels, return_counts=True)
    percentages = counts / len(labels) * 100
  
    # 可视化结果
    fig, axes = plt.subplots(2, 3, figsize=(18, 12))
  
    # 原图
    axes[0,0].imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
    axes[0,0].set_title('原图')
    axes[0,0].axis('off')
  
    # BGR直方图
    axes[0,1].plot(hist_b, color='blue', alpha=0.7, label='Blue')
    axes[0,1].plot(hist_g, color='green', alpha=0.7, label='Green')
    axes[0,1].plot(hist_r, color='red', alpha=0.7, label='Red')
    axes[0,1].set_title('BGR直方图')
    axes[0,1].legend()
  
    # HSV图像
    axes[0,2].imshow(cv2.cvtColor(hsv, cv2.COLOR_HSV2RGB))
    axes[0,2].set_title('HSV图像')
    axes[0,2].axis('off')
  
    # 主要颜色
    color_bar = np.ones((100, 500, 3), dtype=np.uint8)
    x_start = 0
    for i, (center, percentage) in enumerate(zip(centers, percentages)):
        x_end = x_start + int(500 * percentage / 100)
        color_bar[:, x_start:x_end] = center[::-1]  # BGR to RGB
        x_start = x_end
  
    axes[1,0].imshow(color_bar)
    axes[1,0].set_title('主要颜色分布')
    axes[1,0].axis('off')
  
    # 颜色统计信息
    stats_text = f"BGR均值: {np.mean(img, axis=(0,1)):.1f}\n"
    stats_text += f"BGR标准差: {np.std(img, axis=(0,1)):.1f}\n"
    stats_text += f"HSV均值: {np.mean(hsv, axis=(0,1)):.1f}"
  
    axes[1,1].text(0.1, 0.5, stats_text, fontsize=12, 
                   verticalalignment='center', transform=axes[1,1].transAxes)
    axes[1,1].set_title('颜色统计')
    axes[1,1].axis('off')
  
    # 颜色饼图
    colors_rgb = centers / 255.0
    axes[1,2].pie(percentages, colors=colors_rgb, autopct='%1.1f%%')
    axes[1,2].set_title('颜色比例')
  
    plt.tight_layout()
    plt.show()
  
    return {
        'bgr_moments': (b_moments, g_moments, r_moments),
        'hsv_moments': (h_moments, s_moments, v_moments),
        'dominant_colors': centers,
        'color_percentages': percentages,
        'histograms': (hist_b, hist_g, hist_r)
    }

# 使用示例
# color_feats = color_features(img)

形状特征

def shape_features(img):
    """形状特征提取"""
  
    # 转换为灰度图并二值化
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
  
    # 使用Otsu自动阈值
    _, binary = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
  
    # 查找轮廓
    contours, _ = cv2.findContours(binary, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
  
    # 选择最大的轮廓
    if not contours:
        return None
  
    largest_contour = max(contours, key=cv2.contourArea)
  
    # 计算形状特征
    features = {}
  
    # 1. 基本几何特征
    features['area'] = cv2.contourArea(largest_contour)
    features['perimeter'] = cv2.arcLength(largest_contour, True)
  
    # 2. 紧致度（圆形度）
    features['compactness'] = (features['perimeter'] ** 2) / (4 * np.pi * features['area'])
  
    # 3. 边界矩形特征
    x, y, w, h = cv2.boundingRect(largest_contour)
    features['aspect_ratio'] = float(w) / h
    features['extent'] = features['area'] / (w * h)
  
    # 4. 最小外接矩形
    rect = cv2.minAreaRect(largest_contour)
    features['min_area_rect'] = rect
    features['orientation'] = rect[2]
  
    # 5. 凸包特征
    hull = cv2.convexHull(largest_contour)
    hull_area = cv2.contourArea(hull)
    features['solidity'] = features['area'] / hull_area
  
    # 6. 椭圆拟合
    if len(largest_contour) >= 5:
        ellipse = cv2.fitEllipse(largest_contour)
        features['ellipse'] = ellipse
        features['ellipse_area'] = np.pi * ellipse[1][0] * ellipse[1][1] / 4
        features['ellipse_ratio'] = features['area'] / features['ellipse_area']
  
    # 7. 中心矩和Hu矩
    moments = cv2.moments(largest_contour)
    if moments['m00'] != 0:
        features['centroid'] = (int(moments['m10'] / moments['m00']), 
                               int(moments['m01'] / moments['m00']))
    
        # Hu矩（形状描述符）
        hu_moments = cv2.HuMoments(moments)
        features['hu_moments'] = hu_moments.flatten()
  
    # 8. 傅里叶描述符
    contour_complex = np.empty(largest_contour.shape[:-1], dtype=complex)
    contour_complex.real = largest_contour[:,0,0]
    contour_complex.imag = largest_contour[:,0,1]
    fourier_result = np.fft.fft(contour_complex)
    features['fourier_descriptors'] = np.abs(fourier_result)[:20]  # 前20个系数
  
    # 可视化结果
    result_img = img.copy()
  
    # 绘制轮廓
    cv2.drawContours(result_img, [largest_contour], -1, (0, 255, 0), 2)
  
    # 绘制边界矩形
    cv2.rectangle(result_img, (x, y), (x+w, y+h), (255, 0, 0), 2)
  
    # 绘制最小外接矩形
    box = cv2.boxPoints(rect)
    box = np.int0(box)
    cv2.drawContours(result_img, [box], 0, (0, 0, 255), 2)
  
    # 绘制凸包
    cv2.drawContours(result_img, [hull], -1, (255, 255, 0), 2)
  
    # 绘制质心
    if 'centroid' in features:
        cv2.circle(result_img, features['centroid'], 5, (255, 0, 255), -1)
  
    # 绘制椭圆拟合
    if 'ellipse' in features:
        cv2.ellipse(result_img, ellipse, (0, 255, 255), 2)
  
    # 显示结果
    fig, axes = plt.subplots(1, 3, figsize=(18, 6))
  
    axes[0].imshow(cv2.cvtColor(img, cv2.COLOR_BGR2RGB))
    axes[0].set_title('原图')
    axes[0].axis('off')
  
    axes[1].imshow(binary, cmap='gray')
    axes[1].set_title('二值图')
    axes[1].axis('off')
  
    axes[2].imshow(cv2.cvtColor(result_img, cv2.COLOR_BGR2RGB))
    axes[2].set_title('形状分析')
    axes[2].axis('off')
  
    plt.tight_layout()
    plt.show()
  
    # 打印主要特征
    print(f"面积: {features['area']:.2f}")
    print(f"周长: {features['perimeter']:.2f}")
    print(f"紧致度: {features['compactness']:.3f}")
    print(f"长宽比: {features['aspect_ratio']:.3f}")
    print(f"充实度: {features['extent']:.3f}")
    print(f"凸性: {features['solidity']:.3f}")
  
    return features

# 使用示例
# shape_feats = shape_features(img)

4. 频域分析

傅里叶变换分析

def frequency_domain_analysis(img):
    """频域分析"""
  
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
  
    # 1. 二维傅里叶变换
    f_transform = np.fft.fft2(gray)
    f_shift = np.fft.fftshift(f_transform)
  
    # 计算幅度谱和相位谱
    magnitude_spectrum = np.log(np.abs(f_shift) + 1)
    phase_spectrum = np.angle(f_shift)
  
    # 2. 频域滤波
    rows, cols = gray.shape
    crow, ccol = rows // 2, cols // 2
  
    # 低通滤波器（去除高频噪声）
    low_pass = f_shift.copy()
    low_pass[crow-30:crow+30, ccol-30:ccol+30] = 0
    f_ishift_low = np.fft.ifftshift(low_pass)
    img_back_low = np.fft.ifft2(f_ishift_low)
    img_back_low = np.abs(img_back_low)
  
    # 高通滤波器（增强边缘）
    high_pass = f_shift.copy()
    high_pass[crow-30:crow+30, ccol-30:ccol+30] = 0
    f_ishift_high = np.fft.ifftshift(high_pass)
    img_back_high = np.fft.ifft2(f_ishift_high)
    img_back_high = np.abs(img_back_high)
  
    # 3. 频域特征
    # 能量分布
    energy_spectrum = np.abs(f_shift) ** 2
    total_energy = np.sum(energy_spectrum)
  
    # 计算不同频率范围的能量比例
    center_energy = np.sum(energy_spectrum[crow-50:crow+50, ccol-50:ccol+50])
    low_freq_ratio = center_energy / total_energy
  
    # 频域熵
    normalized_spectrum = energy_spectrum / total_energy
    normalized_spectrum = normalized_spectrum[normalized_spectrum > 0]
    spectral_entropy = -np.sum(normalized_spectrum * np.log2(normalized_spectrum))
  
    # 可视化
    fig, axes = plt.subplots(2, 4, figsize=(20, 10))
  
    # 原图
    axes[0,0].imshow(gray, cmap='gray')
    axes[0,0].set_title('原图')
    axes[0,0].axis('off')
  
    # 幅度谱
    axes[0,1].imshow(magnitude_spectrum, cmap='gray')
    axes[0,1].set_title('幅度谱 (对数)')
    axes[0,1].axis('off')
  
    # 相位谱
    axes[0,2].imshow(phase_spectrum, cmap='gray')
    axes[0,2].set_title('相位谱')
    axes[0,2].axis('off')
  
    # 能量谱
    axes[0,3].imshow(np.log(energy_spectrum + 1), cmap='hot')
    axes[0,3].set_title('能量谱')
    axes[0,3].axis('off')
  
    # 低通滤波结果
    axes[1,0].imshow(img_back_low, cmap='gray')
    axes[1,0].set_title('低通滤波')
    axes[1,0].axis('off')
  
    # 高通滤波结果
    axes[1,1].imshow(img_back_high, cmap='gray')
    axes[1,1].set_title('高通滤波')
    axes[1,1].axis('off')
  
    # 径向频率分布
    y, x = np.ogrid[:rows, :cols]
    center = np.array([crow, ccol])
    radius = np.sqrt((x - center[1])**2 + (y - center[0])**2)
  
    # 计算径向平均
    max_radius = int(np.min([crow, ccol, rows-crow, cols-ccol]))
    radial_profile = []
    for r in range(max_radius):
        mask = (radius >= r) & (radius < r+1)
        radial_profile.append(np.mean(energy_spectrum[mask]))
  
    axes[1,2].plot(radial_profile)
    axes[1,2].set_title('径向频率分布')
    axes[1,2].set_xlabel('半径')
    axes[1,2].set_ylabel('平均能量')
  
    # 统计信息
    stats_text = f"总能量: {total_energy:.2e}\n"
    stats_text += f"低频能量比: {low_freq_ratio:.3f}\n"
    stats_text += f"频域熵: {spectral_entropy:.3f}"
  
    axes[1,3].text(0.1, 0.5, stats_text, fontsize=12,
                   verticalalignment='center', transform=axes[1,3].transAxes)
    axes[1,3].set_title('频域统计')
    axes[1,3].axis('off')
  
    plt.tight_layout()
    plt.show()
  
    return {
        'magnitude_spectrum': magnitude_spectrum,
        'phase_spectrum': phase_spectrum,
        'energy_spectrum': energy_spectrum,
        'low_freq_ratio': low_freq_ratio,
        'spectral_entropy': spectral_entropy,
        'radial_profile': radial_profile
    }

# 使用示例
# freq_features = frequency_domain_analysis(img)

5. 多尺度分析

小波变换

def wavelet_analysis(img):
    """小波变换分析（需要PyWavelets库）"""
  
    gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
  
    try:
        import pywt
    
        # 1. 单级小波分解
        coeffs = pywt.dwt2(gray, 'db4')
        cA, (cH, cV, cD) = coeffs
    
        # 2. 多级小波分解
        coeffs_multilevel = pywt.wavedec2(gray, 'db4', level=3)
    
        # 3. 小波特征提取
        # 各子带的能量
        energy_LL = np.sum(cA ** 2)
        energy_LH = np.sum(cH ** 2)
        energy_HL = np.sum(cV ** 2)
        energy_HH = np.sum(cD ** 2)
    
        total_energy = energy_LL + energy_LH + energy_HL + energy_HH
    
        # 归一化能量
        energy_ratios = {
            'LL': energy_LL / total_energy,
            'LH': energy_LH / total_energy,
            'HL': energy_HL / total_energy,
            'HH': energy_HH / total_energy
        }
    
        # 4. 纹理能量度量
        def texture_energy(coeff):
            return np.mean(coeff ** 2)
    
        texture_energies = {
            'LL': texture_energy(cA),
            'LH': texture_energy(cH),
            'HL': texture_energy(cV),
            'HH': texture_energy(cD)
        }
    
        # 可视化
        fig, axes = plt.subplots(2, 3, figsize=(18, 12))
    
        # 原图
        axes[0,0].imshow(gray, cmap='gray')
        axes[0,0].set_title('原图')
        axes[0,0].axis('off')
    
        # 低频分量（近似）
        axes[0,1].imshow(cA, cmap='gray')
        axes[0,1].set_title(f'LL (近似) - 能量比: {energy_ratios["LL"]:.3f}')
        axes[0,1].axis('off')
    
        # 水平细节
        axes[0,2].imshow(cH, cmap='gray')
        axes[0,2].set_title(f'LH (水平细节) - 能量比: {energy_ratios["LH"]:.3f}')
        axes[0,2].axis('off')
    
        # 垂直细节
        axes[1,0].imshow(cV, cmap='gray')
        axes[1,0].set_title(f'HL (垂直细节) - 能量比: {energy_ratios["HL"]:.3f}')
        axes[1,0].axis('off')
    
        # 对角细节
        axes[1,1].imshow(cD, cmap='gray')
        axes[1,1].set_title(f'HH (对角细节) - 能量比: {energy_ratios["HH"]:.3f}')
        axes[1,1].axis('off')
    
        # 重构图像
        reconstructed = pywt.idwt2(coeffs, 'db4')
        axes[1,2].imshow(reconstructed, cmap='gray')
        axes[1,2].set_title('重构图像')
        axes[1,2].axis('off')
    
        plt.tight_layout()
        plt.show()
    
        # 多级分解可视化
        fig, axes = plt.subplots(1, 4, figsize=(20, 5))
    
        for i, coeff in enumerate(coeffs_multilevel[:4]):
            if i == 0:
                axes[i].imshow(coeff, cmap='gray')
                axes[i].set_title(f'Level {len(coeffs_multilevel)-1-i} - 近似')
            else:
                # 对于细节系数，显示其中一个子带
                axes[i].imshow(coeff[0], cmap='gray')
                axes[i].set_title(f'Level {len(coeffs_multilevel)-1-i} - 细节')
            axes[i].axis('off')
    
        plt.tight_layout()
        plt.show()
    
        return {
            'energy_ratios': energy_ratios,
            'texture_energies': texture_energies,
            'coefficients': coeffs,
            'multilevel_coeffs': coeffs_multilevel
        }
    
    except ImportError:
        print("需要安装PyWavelets库: pip install PyWavelets")
        return None

# 使用示例
# wavelet_feats = wavelet_analysis(img)

6. 综合特征提取系统

特征向量构建

class ImageFeatureExtractor:
    """综合图像特征提取器"""
  
    def __init__(self):
        self.features = {}
  
    def extract_all_features(self, img):
        """提取所有类型的特征"""
    
        print("正在提取特征...")
    
        # 1. 颜色特征
        print("- 颜色特征")
        color_feats = self.extract_color_features(img)
    
        # 2. 纹理特征
        print("- 纹理特征")
        texture_feats = self.extract_texture_features(img)
    
        # 3. 形状特征
        print("- 形状特征")
        shape_feats = self.extract_shape_features(img)
    
        # 4. 统计特征
        print("- 统计特征")
        stat_feats = self.extract_statistical_features(img)
    
        # 组合所有特征
        all_features = {
            'color': color_feats,
            'texture': texture_feats,
            'shape': shape_feats,
            'statistical': stat_feats
        }
    
        # 创建特征向量
        feature_vector = self.create_feature_vector(all_features)
    
        print(f"特征提取完成，特征向量维度: {len(feature_vector)}")
    
        return all_features, feature_vector
  
    def extract_color_features(self, img):
        """提取颜色特征"""
        hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)
    
        # 颜色统计矩
        color_moments = []
        for i in range(3):
            channel = img[:,:,i]
            moments = [
                np.mean(channel),
                np.std(channel),
                np.mean((channel - np.mean(channel))**3) / (np.std(channel)**3),  # 偏斜度
            ]
            color_moments.extend(moments)
    
        # HSV统计矩
        hsv_moments = []
        for i in range(3):
            channel = hsv[:,:,i]
            hsv_moments.extend([np.mean(channel), np.std(channel)])
    
        return {
            'bgr_moments': color_moments,
            'hsv_moments': hsv_moments
        }
  
    def extract_texture_features(self, img):
        """提取纹理特征"""
        gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    
        # 简化的GLCM特征
        # 计算灰度差分统计
        diff_horizontal = np.diff(gray, axis=1)
        diff_vertical = np.diff(gray, axis=0)
    
        texture_features = [
            np.mean(np.abs(diff_horizontal)),
            np.std(diff_horizontal),
            np.mean(np.abs(diff_vertical)),
            np.std(diff_vertical)
        ]
    
        return {'glcm_simple': texture_features}
  
    def extract_shape_features(self, img):
        """提取形状特征"""
        gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        _, binary = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
    
        contours, _ = cv2.findContours(binary, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    
        if not contours:
            return {'basic_shape': [0, 0, 0, 0]}
    
        largest_contour = max(contours, key=cv2.contourArea)
    
        # 基本形状特征
        area = cv2.contourArea(largest_contour)
        perimeter = cv2.arcLength(largest_contour, True)
    
        if perimeter > 0:
            compactness = (perimeter ** 2) / (4 * np.pi * area) if area > 0 else 0
        else:
            compactness = 0
    
        # 边界矩形
        x, y, w, h = cv2.boundingRect(largest_contour)
        aspect_ratio = float(w) / h if h > 0 else 0
    
        return {
            'basic_shape': [area, perimeter, compactness, aspect_ratio]
        }
  
    def extract_statistical_features(self, img):
        """提取统计特征"""
        gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    
        # 基本统计量
        stats = [
            np.mean(gray),
            np.std(gray),
            np.min(gray),
            np.max(gray),
            np.median(gray)
        ]
    
        # 直方图统计
        hist = cv2.calcHist([gray], [0], None, [256], [0, 256]).flatten()
        hist_normalized = hist / np.sum(hist)
    
        # 直方图熵
        hist_entropy = -np.sum(hist_normalized * np.log2(hist_normalized + 1e-10))
    
        stats.append(hist_entropy)
    
        return {'basic_stats': stats}
  
    def create_feature_vector(self, all_features):
        """创建特征向量"""
        feature_vector = []
    
        # 展平所有特征
        for category, features in all_features.items():
            for feature_name, feature_values in features.items():
                if isinstance(feature_values, list):
                    feature_vector.extend(feature_values)
                else:
                    feature_vector.append(feature_values)
    
        return np.array(feature_vector)
  
    def compare_images(self, img1, img2):
        """比较两个图像的特征相似性"""
    
        _, features1 = self.extract_all_features(img1)
        _, features2 = self.extract_all_features(img2)
    
        # 欧氏距离
        euclidean_dist = np.linalg.norm(features1 - features2)
    
        # 余弦相似度
        cosine_sim = np.dot(features1, features2) / (np.linalg.norm(features1) * np.linalg.norm(features2))
    
        # 相关系数
        correlation = np.corrcoef(features1, features2)[0, 1]
    
        return {
            'euclidean_distance': euclidean_dist,
            'cosine_similarity': cosine_sim,
            'correlation': correlation
        }

# 使用示例
# extractor = ImageFeatureExtractor()
# features, feature_vector = extractor.extract_all_features(img)
# similarity = extractor.compare_images(img1, img2)

总结

图像分析与特征提取涵盖了高级的图像理解技术：

核心技术

• 边缘分析：Canny、Sobel、Laplacian等算子，边缘密度和方向分析
• 纹理分析：GLCM、LBP、Gabor滤波器等纹理描述方法
• 颜色特征：颜色矩、直方图、主要颜色提取
• 形状特征：几何特征、Hu矩、傅里叶描述符

高级分析

• 频域分析：傅里叶变换、频域滤波、频谱特征
• 多尺度分析：小波变换、多尺度分解
• 统计特征：统计矩、熵、分布特征

特征应用

• 特征向量构建：多种特征融合
• 相似性度量：欧氏距离、余弦相似度、相关系数
• 图像检索：基于内容的图像检索
• 模式识别：图像分类和识别