特征与匹配第10课/共35课

第10课：光流估计

📖 课程概述

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

📑 本课目录

1. 光流基本原理 2. Lucas-Kanade方法 3. Horn-Schunck方法 4. 金字塔光流 5. Farneback稠密光流 6. 光流可视化 7. 运动分析实战

1. 光流基本原理

光流是像素在连续帧之间的表观运动，基于亮度恒常假设。

光流约束方程: I_xu + I_yv + I_t = 0
即: ∇I · v + I_t = 0
1个方程2个未知数→孔径问题

import cv2, numpy as np
f1 = np.zeros((200,200), dtype=np.uint8)
cv2.circle(f1, (80,100), 20, 200, -1)
f2 = np.zeros((200,200), dtype=np.uint8)
cv2.circle(f2, (85,95), 20, 200, -1)  # 右移5,上移5
Ix = cv2.Sobel(f1, cv2.CV_64F, 1, 0, ksize=3)*0.5 + cv2.Sobel(f2, cv2.CV_64F, 1, 0, ksize=3)*0.5
Iy = cv2.Sobel(f1, cv2.CV_64F, 0, 1, ksize=3)*0.5 + cv2.Sobel(f2, cv2.CV_64F, 0, 1, ksize=3)*0.5
It = f2.astype(np.float64) - f1.astype(np.float64)
mask = (f1 > 100) | (f2 > 100)
residual = Ix[mask]*5 + Iy[mask]*(-5) + It[mask]
print(f"光流约束验证(u=5,v=-5): 残差均值={residual.mean():.4f}")

2. Lucas-Kanade方法

LK方法: [A]v = b, v = (A^TA)^-1A^Tb
A^TA = Harris矩阵!

import cv2, numpy as np
f1 = np.zeros((300,300), dtype=np.uint8)
cv2.circle(f1, (100,150), 30, 200, -1)
cv2.rectangle(f1, (200,80), (260,140), 150, -1)
dx, dy = 8, 5
M = np.float32([[1,0,dx],[0,1,dy]])
f2 = cv2.warpAffine(f1, M, (300,300))
corners = cv2.goodFeaturesToTrack(f1, maxCorners=50, qualityLevel=0.01, minDistance=10)
if corners is not None:
    next_pts, status, _ = cv2.calcOpticalFlowPyrLK(f1, f2, corners, None)
    good = status.ravel() == 1
    flow = next_pts[good] - corners[good]
    print(f"LK光流: 追踪{np.sum(good)}/{len(corners)}点")
    print(f"估计: dx={flow[:,0].mean():.2f}, dy={flow[:,1].mean():.2f}")
    print(f"真实: dx={dx}, dy={dy}")
vis = cv2.cvtColor(f2, cv2.COLOR_GRAY2BGR)
if corners is not None and next_pts is not None:
    for i in range(len(corners)):
        if status[i]:
            p1 = corners[i].ravel().astype(int)
            p2 = next_pts[i].ravel().astype(int)
            cv2.arrowedLine(vis, tuple(p1), tuple(p2), (0,255,0), 2)
cv2.imwrite('/var/www/ttl/cv/l10_lk.png', vis)

3. Horn-Schunck方法

HS能量: E = ∫∫ (I_xu+I_yv+I_t)² + α²(|∇u|²+|∇v|²) dxdy

import cv2, numpy as np
f1 = np.zeros((150,150), dtype=np.uint8)
cv2.circle(f1, (60,75), 20, 200, -1)
M = np.float32([[1,0,5],[0,1,3]])
f2 = cv2.warpAffine(f1, M, (150,150))
I1, I2 = f1.astype(np.float64)/255, f2.astype(np.float64)/255
Ix = cv2.Sobel(I1,cv2.CV_64F,1,0,ksize=3)*0.5 + cv2.Sobel(I2,cv2.CV_64F,1,0,ksize=3)*0.5
Iy = cv2.Sobel(I1,cv2.CV_64F,0,1,ksize=3)*0.5 + cv2.Sobel(I2,cv2.CV_64F,0,1,ksize=3)*0.5
It = I2 - I1
u, v = np.zeros_like(Ix), np.zeros_like(Ix)
kernel = np.array([[0,1/4,0],[1/4,0,1/4],[0,1/4,0]])
for _ in range(200):
    ua, va = cv2.filter2D(u,-1,kernel), cv2.filter2D(v,-1,kernel)
    P = Ix*ua + Iy*va + It
    D = 1.0 + Ix**2 + Iy**2
    u, v = ua - Ix*P/D, va - Iy*P/D
mask = (f1>100)|(f2>100)
if np.any(mask):
    print(f"HS alpha=1: u={u[mask].mean():.2f}, v={v[mask].mean():.2f} (真实5,3)")

4. 金字塔光流

import cv2, numpy as np
f1 = np.zeros((300,300), dtype=np.uint8)
cv2.circle(f1, (80,150), 30, 200, -1)
dx, dy = 30, 15  # 大位移
M = np.float32([[1,0,dx],[0,1,dy]])
f2 = cv2.warpAffine(f1, M, (300,300))
corners = cv2.goodFeaturesToTrack(f1, maxCorners=30, qualityLevel=0.01, minDistance=10)
if corners is not None:
    p1, s1, _ = cv2.calcOpticalFlowPyrLK(f1, f2, corners, None, maxLevel=0)
    p3, s3, _ = cv2.calcOpticalFlowPyrLK(f1, f2, corners, None, maxLevel=3)
    f1_ok, f3_ok = s1.ravel()==1, s3.ravel()==1
    if np.any(f1_ok):
        m1 = (p1[f1_ok]-corners[f1_ok]).mean(axis=0)
        print(f"单层LK: dx={m1[0]:.1f}, dy={m1[1]:.1f}")
    if np.any(f3_ok):
        m3 = (p3[f3_ok]-corners[f3_ok]).mean(axis=0)
        print(f"3层金字塔: dx={m3[0]:.1f}, dy={m3[1]:.1f}")
    print(f"真实: dx={dx}, dy={dy}")
    print(f"单层追踪: {np.sum(f1_ok)}/{len(corners)}, 金字塔: {np.sum(f3_ok)}/{len(corners)}")

5. Farneback稠密光流

import cv2, numpy as np
f1 = np.zeros((200,200), dtype=np.uint8)
cv2.circle(f1, (60,100), 25, 200, -1)
cv2.rectangle(f1, (130,70), (170,130), 150, -1)
dx, dy = 10, 5
M = np.float32([[1,0,dx],[0,1,dy]])
f2 = cv2.warpAffine(f1, M, (200,200))
flow = cv2.calcOpticalFlowFarneback(f1, f2, None, 0.5, 3, 15, 3, 5, 1.2, 0)
mask = (f1>50)|(f2>50)
if np.any(mask):
    print(f"Farneback: u={flow[:,:,0][mask].mean():.2f}, v={flow[:,:,1][mask].mean():.2f}")
    print(f"真实: u={dx}, v={dy}")

6. 光流可视化

import cv2, numpy as np
f1 = np.zeros((200,200), dtype=np.uint8)
cv2.circle(f1, (70,100), 30, 200, -1)
cv2.rectangle(f1, (140,60), (180,140), 150, -1)
M = np.float32([[1,0,8],[0,1,4]])
f2 = cv2.warpAffine(f1, M, (200,200))
flow = cv2.calcOpticalFlowFarneback(f1, f2, None, 0.5, 3, 15, 3, 5, 1.2, 0)
hsv = np.zeros((200,200,3), dtype=np.uint8)
mag, ang = cv2.cartToPolar(flow[:,:,0], flow[:,:,1])
hsv[:,:,0] = ang*180/np.pi/2
hsv[:,:,1] = 255
hsv[:,:,2] = cv2.normalize(mag, None, 0, 255, cv2.NORM_MINMAX)
flow_bgr = cv2.cvtColor(hsv, cv2.COLOR_HSV2BGR)
vis = cv2.cvtColor(f2, cv2.COLOR_GRAY2BGR)
for y in range(0,200,10):
    for x in range(0,200,10):
        fx, fy = flow[y,x]
        if fx**2+fy**2 > 0.5:
            cv2.arrowedLine(vis, (x,y), (int(x+fx*3),int(y+fy*3)), (0,255,0), 1)
cv2.imwrite('/var/www/ttl/cv/l10_flow_vis.png', np.hstack([flow_bgr, vis]))
print(f"最大幅值: {mag.max():.2f}")

7. 运动分析实战

import cv2, numpy as np
f1 = np.zeros((300,300), dtype=np.uint8)
cv2.rectangle(f1, (0,250), (300,300), 80, -1)
cv2.circle(f1, (80,200), 20, 200, -1)
f2 = f1.copy()
cv2.circle(f2, (100,195), 20, 200, -1)
cv2.rectangle(f1, (200,180), (240,220), 150, -1)
cv2.rectangle(f2, (170,185), (210,225), 150, -1)
diff = cv2.absdiff(f1, f2)
_, binary = cv2.threshold(diff, 20, 255, cv2.THRESH_BINARY)
flow = cv2.calcOpticalFlowFarneback(f1, f2, None, 0.5, 3, 15, 3, 5, 1.2, 0)
contours, _ = cv2.findContours(binary, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
print(f"检测到 {len(contours)} 个运动区域")
vis = cv2.cvtColor(f2, cv2.COLOR_GRAY2BGR)
for cnt in contours:
    x,y,w,h = cv2.boundingRect(cnt)
    cv2.rectangle(vis, (x,y), (x+w,y+h), (0,255,0), 2)
    region_flow = flow[y:y+h, x:x+w]
    mu, mv = region_flow[:,:,0].mean(), region_flow[:,:,1].mean()
    cv2.arrowedLine(vis, (x+w//2,y+h//2), (int(x+w//2+mu*5),int(y+h//2+mv*5)), (0,0,255), 2)
    print(f"  区域({x},{y}): 运动=({mu:.1f}, {mv:.1f})")
cv2.imwrite('/var/www/ttl/cv/l10_motion.png', vis)

🔑 光流估计核心要点

理解数学原理是正确应用算法的基础
不同方法各有优缺点，需要根据场景选择
实际工程中往往需要组合多种方法
性能优化：考虑实时性、内存占用和精度权衡

📝 课后练习

实现本课所有算法并对比效果
在不同参数下测试算法的鲁棒性
将本课方法应用到自己的数据集
分析算法的计算复杂度和优化方向

延伸阅读：推荐阅读《Digital Image Processing》(Gonzalez)相关章节，以及OpenCV官方文档中的详细API说明。实际项目中，建议先在简单数据上验证算法，再迁移到复杂场景。

📊 方法对比总结

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

📐 关键公式速查

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

卷积定理: 空间域卷积 = 频域乘积
梯度计算: 一阶差分近似偏导数
卷积参数量: C_out × (C_in × K² + 1)
感受野: RF(n) = RF(n-1) + (K-1) × ∏stride

🔬 光流估计进阶内容

本节深入探讨光流估计的进阶话题和工程实践中的关键考虑。

工业实践要点

算法选型：根据精度、速度、资源约束选择合适的方法
数据质量：高质量数据比复杂算法更重要
部署优化：模型压缩、量化、TensorRT加速
监控维护：线上模型需要持续监控性能退化

⚠️ 常见陷阱与解决方案

过拟合：数据增强+正则化+早停
类别不平衡：重采样+Focal Loss
推理速度慢：模型剪枝+量化+批处理优化
长尾分布：数据重平衡+罕见类别增强

📊 性能优化策略

优化方向	方法	精度影响	加速比
模型压缩	剪枝/蒸馏	1-3%下降	2-4x
量化	INT8/FP16	<1%下降	2-8x
算子融合	TensorRT/ONNX	无	1.5-3x
批处理	增大batch size	无	线性

💻 光流完整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为优秀")

✅ 实机验证

光流约束方程验证
LK光流小运动准确
金字塔LK解决大运动

⬅️ 第9课 🏠 目录第11课 ➡️

🏆

光流追踪者

你已经掌握了光流估计的核心知识，继续前进！

⬅️ 第9课 🏠 目录第11课 ➡️

第10课：光流估计

📖 课程概述

📑 本课目录

1. 光流基本原理

2. Lucas-Kanade方法

3. Horn-Schunck方法

4. 金字塔光流

5. Farneback稠密光流

6. 光流可视化

7. 运动分析实战

🔑 光流估计核心要点

📝 课后练习

📊 方法对比总结

📐 关键公式速查

🔬 光流估计进阶内容

工业实践要点

最新研究进展

⚠️ 常见陷阱与解决方案

📊 性能优化策略

💻 光流完整Pipeline代码

🔧 参数调优指南

📊 完整评估指标

✅ 实机验证

光流追踪者