基于运动补偿的前景检测算法

这段代码实现了基于运动补偿的前景检测算法。

主要功能包括：

• 运动补偿模块：使用基于网格的 KLT 特征跟踪算法计算两帧之间的运动，然后通过单应性变换实现帧间运动补偿。
• 前景检测模块：结合两帧运动补偿结果，通过帧间差分计算前景掩码。
• 异常处理：添加了图像加载检查和异常捕获，提高了代码的健壮性。
• 路径处理：自动创建保存目录，避免因目录不存在导致的错误。

使用时需要提供三帧连续图像：两个参考帧和当前帧。代码会计算出前景掩码并保存为图像文件。

import cv2
import numpy as np
import os
import sys


def motion_compensate(frame1, frame2):
    """
    使用基于网格的KLT特征跟踪实现两帧之间的运动补偿
    
    参数:
    frame1: 前一帧图像(BGR格式)
    frame2: 当前帧图像(BGR格式)
    
    返回:
    compensated: 运动补偿后的图像
    mask: 补偿区域的掩码
    avg_dst: 平均运动距离
    motion_x: x方向平均运动量
    motion_y: y方向平均运动量
    homography_matrix: 单应性变换矩阵
    """
    # 设置LK光流参数
    lk_params = dict(winSize=(15, 15), maxLevel=3,
                     criteria=(cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 30, 0.003))

    # 图像预处理和网格点生成
    width = frame2.shape[1]
    height = frame2.shape[0]
    scale = 2  # 放大图像以获得更精确的跟踪
    
    # 调整图像大小以提高特征点检测精度
    frame1_grid = cv2.resize(frame1, (960 * scale, 540 * scale), dst=None, interpolation=cv2.INTER_CUBIC)
    frame2_grid = cv2.resize(frame2, (960 * scale, 540 * scale), dst=None, interpolation=cv2.INTER_CUBIC)

    width_grid = frame2_grid.shape[1]
    height_grid = frame2_grid.shape[0]
    gridSizeW = 32 * 2  # 网格宽度
    gridSizeH = 24 * 2  # 网格高度
    
    # 生成网格点作为特征点
    p1 = []
    grid_numW = int(width_grid / gridSizeW - 1)
    grid_numH = int(height_grid / gridSizeH - 1)
    for i in range(grid_numW):
        for j in range(grid_numH):
            # 将点放置在每个网格中心
            point = (np.float32(i * gridSizeW + gridSizeW / 2.0), np.float32(j * gridSizeH + gridSizeH / 2.0))
            p1.append(point)

    p1 = np.array(p1)
    pts_num = grid_numW * grid_numH
    pts_prev = p1.reshape(pts_num, 1, 2)

    # 计算光流
    pts_cur, st, err = cv2.calcOpticalFlowPyrLK(frame1_grid, frame2_grid, pts_prev, None, **lk_params)

    # 选择跟踪成功的点
    good_new = pts_cur[st == 1]  # 当前帧中的跟踪点
    good_old = pts_prev[st == 1]  # 前一帧中的跟踪点

    # 计算运动距离和位移
    motion_distance = []
    translate_x = []
    translate_y = []
    for i, (new, old) in enumerate(zip(good_new, good_old)):
        a, b = new.ravel()
        c, d = old.ravel()
        motion_distance0 = np.sqrt((a - c) * (a - c) + (b - d) * (b - d))

        # 过滤异常大的运动值
        if motion_distance0 > 50:
            continue

        translate_x0 = a - c
        translate_y0 = b - d

        motion_distance.append(motion_distance0)
        translate_x.append(translate_x0)
        translate_y.append(translate_y0)

    motion_dist = np.array(motion_distance)
    motion_x = np.mean(np.array(translate_x)) if translate_x else 0
    motion_y = np.mean(np.array(translate_y)) if translate_y else 0

    avg_dst = np.mean(motion_dist) if motion_dist.size > 0 else 0

    # 计算单应性变换矩阵
    if len(good_old) < 15:
        # 点太少时使用近似恒等变换
        homography_matrix = np.array([[0.999, 0, 0], [0, 0.999, 0], [0, 0, 1]])
    else:
        # 使用RANSAC算法估计单应性矩阵
        homography_matrix, status = cv2.findHomography(good_new, good_old, cv2.RANSAC, 3.0)

    # 应用单应性变换进行运动补偿
    compensated = cv2.warpPerspective(frame1, homography_matrix, (width, height), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)

    # 生成掩码以指示变换区域
    vertex = np.array([[0, 0], [width, 0], [width, height], [0, height]], dtype=np.float32).reshape(-1, 1, 2)
    homo_inv = np.linalg.inv(homography_matrix)
    vertex_trans = cv2.perspectiveTransform(vertex, homo_inv)
    vertex_transformed = np.array(vertex_trans, dtype=np.int32).reshape(1, 4, 2)
    im = np.zeros(frame1.shape[:2], dtype='uint8')
    cv2.polylines(im, vertex_transformed, 1, 255)
    cv2.fillPoly(im, vertex_transformed, 255)
    mask = 255 - im

    return compensated, mask, avg_dst, motion_x, motion_y, homography_matrix


def FD_mask(lastFrame1, lastFrame2, currentFrame, save_path='data/mask.jpg'):
    """
    使用两帧运动补偿计算前景掩码
    
    参数:
    lastFrame1: 第一参考帧(BGR格式)
    lastFrame2: 第二参考帧(BGR格式)
    currentFrame: 当前帧(BGR格式)
    save_path: 结果掩码保存路径
    """
    # 图像预处理:高斯模糊和灰度转换
    lastFrame1 = cv2.GaussianBlur(lastFrame1, (11, 11), 0)
    lastFrame1 = cv2.cvtColor(lastFrame1, cv2.COLOR_BGR2GRAY)

    lastFrame2 = cv2.GaussianBlur(lastFrame2, (11, 11), 0)
    lastFrame2 = cv2.cvtColor(lastFrame2, cv2.COLOR_BGR2GRAY)

    currentFrame = cv2.GaussianBlur(currentFrame, (11, 11), 0)
    currentFrame = cv2.cvtColor(currentFrame, cv2.COLOR_BGR2GRAY)

    # 计算第一参考帧到第二参考帧的运动补偿
    img_compensate1, mask1, avg_dist1, motion_x1, motion_y1, homo_matrix = motion_compensate(lastFrame1, lastFrame2)
    frameDiff1 = cv2.absdiff(lastFrame2, img_compensate1)

    # 计算当前帧到第二参考帧的运动补偿
    img_compensate2, mask2, avg_dist2, motion_x2, motion_y2, homo_matrix2 = motion_compensate(currentFrame, lastFrame2)
    frameDiff2 = cv2.absdiff(lastFrame2, img_compensate2)

    # 融合两个差分结果
    frameDiff = (frameDiff1 + frameDiff2) / 2

    # 保存结果
    os.makedirs(os.path.dirname(save_path), exist_ok=True)
    cv2.imwrite(save_path, frameDiff)

    print(f'前景掩码已保存至: {save_path}')
    return frameDiff


if __name__ == "__main__":
    # 示例:加载三帧图像并计算前景掩码
    # 请确保这些图像存在，或者修改为您自己的图像路径
    try:
        lastFrame1 = cv2.imread('data/Test_images/images/phantom05_0600.jpg')
        lastFrame3 = cv2.imread('data/Test_images/images/phantom05_0602.jpg')
        currentFrame = cv2.imread('data/Test_images/images/phantom05_0604.jpg')

        if lastFrame1 is None or lastFrame3 is None or currentFrame is None:
            print("错误: 无法加载图像，请检查文件路径!")
        else:
            FD_mask(lastFrame1, lastFrame3, currentFrame)
    except Exception as e:
        print(f"程序执行出错: {e}")    

#include <opencv2/opencv.hpp>
#include <vector>
#include <iostream>
#include <numeric>
#include <cmath>

using namespace cv;
using namespace std;

struct MotionCompensationResult {
    Mat compensated;
    Mat mask;
    float avg_dst;
    float motion_x;
    float motion_y;
    Mat homography_matrix;
};

MotionCompensationResult motion_compensate(const Mat& frame1, const Mat& frame2) {
    // KLT 跟踪参数
    TermCriteria term_criteria(TermCriteria::EPS | TermCriteria::COUNT, 30, 0.003);
    Size win_size(15, 15);
    int max_level = 3;

    // 图像缩放
    int scale = 2;
    Mat frame1_grid, frame2_grid;
    resize(frame1, frame1_grid, Size(960 * scale, 540 * scale), 0, 0, INTER_CUBIC);
    resize(frame2, frame2_grid, Size(960 * scale, 540 * scale), 0, 0, INTER_CUBIC);

    // 创建网格点
    int gridSizeW = 32 * 2;
    int gridSizeH = 24 * 2;
    int grid_numW = static_cast<int>(frame2_grid.cols / gridSizeW - 1);
    int grid_numH = static_cast<int>(frame2_grid.rows / gridSizeH - 1);
    
    vector<Point2f> p1;
    for (int i = 0; i < grid_numW; i++) {
        for (int j = 0; j < grid_numH; j++) {
            p1.push_back(Point2f(i * gridSizeW + gridSizeW / 2.0f, j * gridSizeH + gridSizeH / 2.0f));
        }
    }

    int pts_num = grid_numW * grid_numH;
    Mat pts_prev = Mat(p1).reshape(2, pts_num);

    // 计算光流
    vector<Point2f> pts_cur;
    vector<uchar> status;
    vector<float> err;
    calcOpticalFlowPyrLK(frame1_grid, frame2_grid, pts_prev, pts_cur, status, err, win_size, max_level, term_criteria);

    // 筛选好点
    vector<Point2f> good_new, good_old;
    for (size_t i = 0; i < status.size(); i++) {
        if (status[i]) {
            good_new.push_back(pts_cur[i]);
            good_old.push_back(p1[i]);
        }
    }

    // 计算运动距离和位移
    vector<float> motion_distance;
    vector<float> translate_x, translate_y;
    for (size_t i = 0; i < good_new.size(); i++) {
        float dx = good_new[i].x - good_old[i].x;
        float dy = good_new[i].y - good_old[i].y;
        float dist = sqrt(dx * dx + dy * dy);

        if (dist > 50) continue;

        motion_distance.push_back(dist);
        translate_x.push_back(dx);
        translate_y.push_back(dy);
    }

    // 计算平均值
    float avg_dst = 0, motion_x = 0, motion_y = 0;
    if (!motion_distance.empty()) {
        avg_dst = accumulate(motion_distance.begin(), motion_distance.end(), 0.0f) / motion_distance.size();
    }
    if (!translate_x.empty()) {
        motion_x = accumulate(translate_x.begin(), translate_x.end(), 0.0f) / translate_x.size();
        motion_y = accumulate(translate_y.begin(), translate_y.end(), 0.0f) / translate_y.size();
    }

    // 计算单应性矩阵
    Mat homography_matrix;
    if (good_old.size() < 15) {
        homography_matrix = (Mat_<double>(3, 3) << 0.999, 0, 0, 0, 0.999, 0, 0, 0, 1);
    } else {
        homography_matrix = findHomography(good_new, good_old, RANSAC, 3.0);
    }

    // 运动补偿
    Mat compensated;
    warpPerspective(frame1, compensated, homography_matrix, Size(frame1.cols, frame1.rows), INTER_LINEAR + WARP_INVERSE_MAP);

    // 计算掩膜
    vector<Point2f> vertex;
    vertex.push_back(Point2f(0, 0));
    vertex.push_back(Point2f(frame1.cols, 0));
    vertex.push_back(Point2f(frame1.cols, frame1.rows));
    vertex.push_back(Point2f(0, frame1.rows));
    
    Mat vertex_mat = Mat(vertex).reshape(2);
    Mat homo_inv = homography_matrix.inv();
    
    vector<Point2f> vertex_trans;
    perspectiveTransform(vertex_mat, vertex_trans, homo_inv);
    
    vector<Point> vertex_transformed;
    for (const auto& pt : vertex_trans) {
        vertex_transformed.push_back(Point(static_cast<int>(pt.x), static_cast<int>(pt.y)));
    }
    
    Mat mask = Mat::zeros(frame1.size(), CV_8UC1);
    vector<vector<Point>> contours;
    contours.push_back(vertex_transformed);
    polylines(mask, contours, true, Scalar(255), 1);
    fillPoly(mask, contours, Scalar(255));
    mask = 255 - mask;

    return {compensated, mask, avg_dst, motion_x, motion_y, homography_matrix};
}

void FD_mask(const Mat& lastFrame1, const Mat& lastFrame2, const Mat& currentFrame, const string& save_path = "mask.jpg") {
    // 图像预处理
    Mat lastFrame1_gray, lastFrame2_gray, currentFrame_gray;
    GaussianBlur(lastFrame1, lastFrame1_gray, Size(11, 11), 0);
    GaussianBlur(lastFrame2, lastFrame2_gray, Size(11, 11), 0);
    GaussianBlur(currentFrame, currentFrame_gray, Size(11, 11), 0);
    
    cvtColor(lastFrame1_gray, lastFrame1_gray, COLOR_BGR2GRAY);
    cvtColor(lastFrame2_gray, lastFrame2_gray, COLOR_BGR2GRAY);
    cvtColor(currentFrame_gray, currentFrame_gray, COLOR_BGR2GRAY);

    // 第一组运动补偿
    auto result1 = motion_compensate(lastFrame1_gray, lastFrame2_gray);
    Mat frameDiff1;
    absdiff(lastFrame2_gray, result1.compensated, frameDiff1);

    // 第二组运动补偿
    auto result2 = motion_compensate(currentFrame_gray, lastFrame2_gray);
    Mat frameDiff2;
    absdiff(lastFrame2_gray, result2.compensated, frameDiff2);

    // 计算最终差分
    Mat frameDiff;
    frameDiff1.convertTo(frameDiff1, CV_32F);
    frameDiff2.convertTo(frameDiff2, CV_32F);
    frameDiff = (frameDiff1 + frameDiff2) / 2;
    frameDiff.convertTo(frameDiff, CV_8U);

    // 保存结果
    imwrite(save_path, frameDiff);
    cout << "done!" << endl;
}

int main() {
    // 读取图像
    Mat lastFrame1 = imread("data/Test_images/images/phantom05_0600.jpg");
    Mat lastFrame3 = imread("data/Test_images/images/phantom05_0602.jpg");
    Mat currentFrame = imread("data/Test_images/images/phantom05_0604.jpg");
    
    /*
    Mat lastFrame1 = imread("250514_430.bmp");
    Mat lastFrame3 = imread("250514_463.bmp");
    Mat currentFrame = imread("250514_490.bmp");
    */

    // 检查图像是否成功加载
    if (lastFrame1.empty() || lastFrame3.empty() || currentFrame.empty()) {
        cout << "无法加载图像!" << endl;
        return -1;
    }

    // 执行帧间差分
    FD_mask(lastFrame1, lastFrame3, currentFrame);

    return 0;
}