图像处理基础 第4课/共35课

第04课：形态学操作

📖 课程概述

本课学习形态学操作的核心原理与实现。我们将从数学基础出发，通过代码实践真正理解每一个算法的运作机制。

📑 本课目录

1. 形态学基础概念2. 腐蚀与膨胀3. 开运算与闭运算4. 形态学梯度5. 顶帽与黑帽6. 结构元素设计7. 形态学实战应用

1. 形态学基础概念

数学形态学基于集合论，用结构元素探测图像中的形状特征。

膨胀: A ⊕ B = {z | (B̂)z ∩ A ≠ ∅}
腐蚀: A ⊖ B = {z | Bz ⊆ A}
import cv2, numpy as np
shapes = {
    '矩形': cv2.getStructuringElement(cv2.MORPH_RECT, (5,5)),
    '椭圆': cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5)),
    '十字': cv2.getStructuringElement(cv2.MORPH_CROSS, (5,5)),
}
for name, se in shapes.items():
    print(f"{name}SE(5x5): 非零={np.count_nonzero(se)}")

2. 腐蚀与膨胀

腐蚀: (A ⊖ B)(x,y) = min(i,j)∈B A(x+i, y+j)
膨胀: (A ⊕ B)(x,y) = max(i,j)∈B A(x+i, y+j)
import cv2, numpy as np
img = np.zeros((256,256), dtype=np.uint8)
cv2.circle(img, (100,100), 50, 255, -1)
cv2.rectangle(img, (30,180), (80,230), 255, -1)
noise = np.random.random(img.shape)
img_noisy = img.copy()
img_noisy[noise < 0.03] = 255
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5))
eroded = cv2.erode(img_noisy, kernel, iterations=1)
dilated = cv2.dilate(img_noisy, kernel, iterations=1)
print(f"原始: {np.count_nonzero(img_noisy)}, 腐蚀: {np.count_nonzero(eroded)}, 膨胀: {np.count_nonzero(dilated)}")
cv2.imwrite('/var/www/ttl/cv/l04_erode_dilate.png', np.hstack([img_noisy, eroded, dilated]))

3. 开运算与闭运算

开运算: A ˆ B = (A ⊖ B) ⊕ B (先腐蚀再膨胀→去除小白点)
闭运算: A • B = (A ⊕ B) ⊖ B (先膨胀再腐蚀→填充小黑洞)
import cv2, numpy as np
img = np.zeros((256,256), dtype=np.uint8)
cv2.circle(img, (100,100), 60, 255, -1)
np.random.seed(42)
img_noisy = img.copy()
img_noisy[np.random.random(img.shape) < 0.02] = 255
fg = img > 0
img_noisy[fg & (np.random.random(img.shape) < 0.02)] = 0
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5))
opened = cv2.morphologyEx(img_noisy, cv2.MORPH_OPEN, kernel)
closed = cv2.morphologyEx(img_noisy, cv2.MORPH_CLOSE, kernel)
for name, r in [('原始',img_noisy),('开运算',opened),('闭运算',closed)]:
    white = (r[~fg] > 0).sum()
    black = (r[fg] == 0).sum()
    print(f"{name}: 白点={white}, 黑洞={black}")
cv2.imwrite('/var/www/ttl/cv/l04_open_close.png', np.hstack([img_noisy, opened, closed]))

4. 形态学梯度

基本梯度: G = A ⊕ B - A ⊖ B (膨胀-腐蚀=边缘)
内部梯度: Gin = A - A ⊖ B (原图-腐蚀=内边缘)
import cv2, numpy as np
img = np.zeros((256,256), dtype=np.uint8)
cv2.circle(img, (128,128), 80, 255, -1)
cv2.circle(img, (128,128), 40, 0, -1)
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (5,5))
gradient = cv2.morphologyEx(img, cv2.MORPH_GRADIENT, kernel)
internal = img - cv2.erode(img, kernel)
print(f"基本梯度: {np.count_nonzero(gradient)}, 内部: {np.count_nonzero(internal)}")
cv2.imwrite('/var/www/ttl/cv/l04_gradient.png', np.hstack([img, gradient, internal]))

5. 顶帽与黑帽

顶帽: T = A - (A ˆ B) (原图-开运算=亮结构)
黑帽: Bh = (A • B) - A (闭运算-原图=暗结构)
import cv2, numpy as np
img = np.zeros((256,256), dtype=np.uint8)
for i in range(256):
    for j in range(256):
        img[i,j] = int(128 + 60*np.sin(i/40) + 40*np.cos(j/50))
cv2.putText(img, 'CV', (80,150), cv2.FONT_HERSHEY_SIMPLEX, 3, 30, 4)
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (51,51))
tophat = cv2.morphologyEx(img, cv2.MORPH_TOPHAT, kernel)
blackhat = cv2.morphologyEx(img, cv2.MORPH_BLACKHAT, kernel)
corrected = img - tophat + blackhat
print(f"原图std: {img.std():.1f}, 校正后std: {corrected.std():.1f}")
cv2.imwrite('/var/www/ttl/cv/l04_tophat.png', np.hstack([img, tophat, blackhat, corrected]))

6. 结构元素设计

import cv2, numpy as np
img = np.zeros((256,256), dtype=np.uint8)
cv2.line(img, (20,128), (236,128), 255, 2)  # 水平线
cv2.line(img, (128,20), (128,236), 255, 2)  # 垂直线
se_h = cv2.getStructuringElement(cv2.MORPH_RECT, (15,3))  # 水平SE
se_v = cv2.getStructuringElement(cv2.MORPH_RECT, (3,15))  # 垂直SE
h_open = cv2.morphologyEx(img, cv2.MORPH_OPEN, se_h)
v_open = cv2.morphologyEx(img, cv2.MORPH_OPEN, se_v)
print(f"水平SE保留水平线: {np.count_nonzero(h_open[126:130,:])}, 垂直线: {np.count_nonzero(h_open[:,126:130])}")
print(f"垂直SE保留垂直线: {np.count_nonzero(v_open[:,126:130])}, 水平线: {np.count_nonzero(v_open[126:130,:])}")
cv2.imwrite('/var/www/ttl/cv/l04_se.png', np.hstack([img, h_open, v_open]))

7. 形态学实战应用

import cv2, numpy as np
img = np.zeros((300,500), dtype=np.uint8)
cv2.rectangle(img, (100,100), (400,200), 255, -1)
cv2.rectangle(img, (100,100), (400,200), 0, 3)
for i, ch in enumerate('ABCDEFGH'):
    cv2.putText(img, ch, (115+i*35, 165), cv2.FONT_HERSHEY_SIMPLEX, 0.8, 0, 2)
img[np.random.random(img.shape) < 0.005] = 255

k1 = cv2.getStructuringElement(cv2.MORPH_RECT, (25,5))
closed = cv2.morphologyEx(img, cv2.MORPH_CLOSE, k1)
k2 = cv2.getStructuringElement(cv2.MORPH_RECT, (15,3))
opened = cv2.morphologyEx(closed, cv2.MORPH_OPEN, k2)
k3 = cv2.getStructuringElement(cv2.MORPH_RECT, (20,10))
dilated = cv2.dilate(opened, k3)
contours, _ = cv2.findContours(dilated, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
print(f"检测到 {len(contours)} 个区域")
for cnt in contours:
    x,y,w,h = cv2.boundingRect(cnt)
    aspect = w / max(h, 1)
    if 2 < aspect < 6 and cv2.contourArea(cnt) > 2000:
        print(f"  疑似车牌! 位置=({x},{y}), {w}x{h}, 宽高比={aspect:.2f}")
cv2.imwrite('/var/www/ttl/cv/l04_plate.png', np.hstack([img, closed, opened, dilated]))

🔑 形态学操作核心要点

📝 课后练习

  1. 实现本课所有算法并对比效果
  2. 在不同参数下测试算法的鲁棒性
  3. 将本课方法应用到自己的数据集
  4. 分析算法的计算复杂度和优化方向
延伸阅读：推荐阅读《Digital Image Processing》(Gonzalez)相关章节，以及OpenCV官方文档中的详细API说明。实际项目中，建议先在简单数据上验证算法，再迁移到复杂场景。

📊 方法对比总结

方法优点缺点适用场景
方法A简单高效精度有限快速原型
方法B精度高计算量大离线处理
方法C平衡精度与速度参数调优复杂实际应用

📐 关键公式速查

本课涉及的核心数学公式汇总，方便快速参考：

🔬 深入理解：形态学运算的数学基础

形态学运算建立在集合论和格论的基础上，具有严格的数学性质：

对偶性: (A^c ⊖ B) = (A ⊕ B)^c (腐蚀和膨胀互为对偶)
开闭对偶: (A ∘ B)^c = A^c • B^c
单调性: A ⊆ B ⟹ A⊕C ⊆ B⊕C
扩展性: A ⊆ A⊕B (膨胀扩展)
反扩展性: A⊖B ⊆ A (腐蚀收缩)

🔬 粒度测量(Granulometry)

通过逐步增大结构元素并统计剩余前景面积，可以分析图像中物体的大小分布：

import cv2, numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# 粒度测量
img = np.zeros((256,256), dtype=np.uint8)
np.random.seed(42)
for _ in range(30):
    x, y = np.random.randint(20,236,2)
    r = np.random.randint(5,25)
    cv2.circle(img, (x,y), r, 200, -1)

sizes = range(1, 30, 2)
areas = []
for s in sizes:
    kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (s,s))
    opened = cv2.morphologyEx(img, cv2.MORPH_OPEN, kernel)
    areas.append(np.count_nonzero(opened))

plt.figure(figsize=(8,4))
plt.plot(sizes, areas, 'o-')
plt.xlabel('SE Size'); plt.ylabel('Remaining Area')
plt.title('Granulometry (Size Distribution)')
plt.savefig('/var/www/ttl/cv/l04_granulometry.png', dpi=100)
plt.close()
print("粒度测量: 面积随SE增大而下降，下降速率反映大小分布")

形态学重建(Morphological Reconstruction)

迭代膨胀直到稳定，比单次开/闭运算更强大：

💻 形态学完整Pipeline代码

以下是一个完整的形态学处理管道，从数据准备到结果评估，包含所有关键步骤和参数调优建议：

import cv2
import numpy as np
import time

class 形态学Pipeline:
    """完整的形态学处理管道"""
    
    def __init__(self, params=None):
        self.params = params or {}
        self.results = {}
    
    def preprocess(self, image):
        """预处理步骤"""
        # 1. 去噪
        if self.params.get('denoise', True):
            image = cv2.GaussianBlur(image, (3, 3), 1.0)
        
        # 2. 对比度增强
        if self.params.get('enhance', False):
            if len(image.shape) == 2:
                image = cv2.equalizeHist(image)
            else:
                lab = cv2.cvtColor(image, cv2.COLOR_BGR2LAB)
                lab[:,:,0] = cv2.equalizeHist(lab[:,:,0])
                image = cv2.cvtColor(lab, cv2.COLOR_LAB2BGR)
        
        # 3. 尺寸调整
        max_size = self.params.get('max_size', 1024)
        h, w = image.shape[:2]
        if max(h, w) > max_size:
            scale = max_size / max(h, w)
            image = cv2.resize(image, (int(w*scale), int(h*scale)))
        
        return image
    
    def process(self, image):
        """核心处理步骤"""
        t0 = time.time()
        preprocessed = self.preprocess(image)
        t1 = time.time()
        
        # 主处理逻辑
        result = self._main_process(preprocessed)
        t2 = time.time()
        
        self.results['timing'] = {
            'preprocess': (t1-t0)*1000,
            'process': (t2-t1)*1000,
            'total': (t2-t0)*1000
        }
        return result
    
    def _main_process(self, image):
        """主处理逻辑（子类可重写）"""
        return image
    
    def evaluate(self, result, ground_truth=None):
        """评估处理结果"""
        metrics = {}
        if ground_truth is not None:
            if len(result.shape) == 2:
                mse = np.mean((result.astype(float) - ground_truth.astype(float))**2)
                metrics['mse'] = mse
                metrics['psnr'] = 10 * np.log10(255**2 / (mse + 1e-10))
        metrics['timing'] = self.results.get('timing', {})
        return metrics
    
    def visualize(self, image, result):
        """可视化对比"""
        if len(image.shape) == 2:
            image = cv2.cvtColor(image, cv2.COLOR_GRAY2BGR)
        if len(result.shape) == 2:
            result = cv2.cvtColor(result, cv2.COLOR_GRAY2BGR)
        h1, w1 = image.shape[:2]
        h2, w2 = result.shape[:2]
        h, w = max(h1,h2), w1+w2+10
        vis = np.zeros((h, w, 3), dtype=np.uint8)
        vis[:h1, :w1] = image
        vis[:h2, w1+10:] = result
        cv2.putText(vis, 'Input', (10, 25), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0,255,0), 2)
        cv2.putText(vis, 'Output', (w1+20, 25), cv2.FONT_HERSHEY_SIMPLEX, 0.7, (0,255,0), 2)
        return vis


# 使用示例
pipeline = 形态学Pipeline({'denoise': True, 'enhance': False, 'max_size': 512})
test_img = np.random.randint(0, 256, (256, 256, 3), dtype=np.uint8)
result = pipeline.process(test_img)
metrics = pipeline.evaluate(result)
print(f"处理时间: {metrics['timing'].get('total', 0):.1f}ms")

🔧 参数调优指南

不同参数对结果的影响：

建议从默认参数开始，根据结果逐步调整。先在少量数据上快速迭代，确定参数后再全量处理。

生产环境建议：1) 添加输入验证和异常处理 2) 记录每次处理的参数和结果 3) 设置性能监控和告警 4) 定期在测试集上评估质量 5) 保持模型/算法版本管理

📊 完整评估指标

import numpy as np

class Metrics:
    @staticmethod
    def mse(img1, img2):
        return np.mean((img1.astype(float) - img2.astype(float))**2)
    
    @staticmethod
    def psnr(img1, img2, max_val=255):
        mse_val = Metrics.mse(img1, img2)
        if mse_val == 0: return float('inf')
        return 10 * np.log10(max_val**2 / mse_val)
    
    @staticmethod
    def iou(mask1, mask2):
        inter = np.logical_and(mask1 > 0, mask2 > 0).sum()
        union = np.logical_or(mask1 > 0, mask2 > 0).sum()
        return inter / (union + 1e-10)
    
    @staticmethod
    def f1_score(pred, gt, threshold=0.5):
        pred_bin = (pred > threshold)
        gt_bin = (gt > 0)
        tp = np.logical_and(pred_bin, gt_bin).sum()
        precision = tp / (pred_bin.sum() + 1e-10)
        recall = tp / (gt_bin.sum() + 1e-10)
        return 2 * precision * recall / (precision + recall + 1e-10)

# 示例
img1 = np.random.randint(0, 256, (100,100), dtype=np.uint8)
img2 = img1 + np.random.randint(-10, 10, (100,100))
img2 = np.clip(img2, 0, 255).astype(np.uint8)

print(f"MSE: {Metrics.mse(img1,img2):.2f}")
print(f"PSNR: {Metrics.psnr(img1,img2):.2f} dB")
print(f"说明: PSNR>30dB为良好, >40dB为优秀")

✅ 实机验证

⬅️ 第3课🏠 目录第5课 ➡️
🏆

形态学专家

你已经掌握了形态学操作的核心知识，继续前进！

⬅️ 第3课🏠 目录第5课 ➡️