写出边缘跟踪的c++代码（不使用OpenCV与MATLAB）

边缘跟踪是一种经典的计算机视觉算法，可以用于检测图像中物体的轮廓。下面是一个简单的C++实现，使用Sobel算子进行边缘检测，再利用非极大值抑制和双阈值处理得到二值化的边缘图像，最后进行边缘跟踪。

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

// Sobel算子
int sobel_x[3][3] = {
    {-1, 0, 1},
    {-2, 0, 2},
    {-1, 0, 1}
};

int sobel_y[3][3] = {
    {-1, -2, -1},
    {0, 0, 0},
    {1, 2, 1}
};

// 非极大值抑制
void non_maximum_suppression(vector<vector<int>>& grad, vector<vector<double>>& theta) {
    int rows = grad.size();
    int cols = grad[0].size();
    for (int i = 1; i < rows - 1; i++) {
        for (int j = 1; j < cols - 1; j++) {
            double angle = theta[i][j] * 180 / M_PI;
            if ((angle >= -22.5 && angle < 22.5) || (angle >= 157.5 && angle <= -157.5)) {
                if (grad[i][j] < grad[i][j - 1] || grad[i][j] < grad[i][j + 1]) {
                    grad[i][j] = 0;
                }
            } else if ((angle >= 22.5 && angle < 67.5) || (angle >= -157.5 && angle < -112.5)) {
                if (grad[i][j] < grad[i - 1][j - 1] || grad[i][j] < grad[i + 1][j + 1]) {
                    grad[i][j] = 0;
                }
            } else if ((angle >= 67.5 && angle < 112.5) || (angle >= -112.5 && angle < -67.5)) {
                if (grad[i][j] < grad[i - 1][j] || grad[i][j] < grad[i + 1][j]) {
                    grad[i][j] = 0;
                }
            } else if ((angle >= 112.5 && angle <= 157.5) || (angle >= -67.5 && angle <= -22.5)) {
                if (grad[i][j] < grad[i - 1][j + 1] || grad[i][j] < grad[i + 1][j - 1]) {
                    grad[i][j] = 0;
                }
            }
        }
    }
}

// 双阈值处理
void double_thresholding(vector<vector<int>>& grad, int T1, int T2) {
    int rows = grad.size();
    int cols = grad[0].size();
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            if (grad[i][j] < T1) {
                grad[i][j] = 0;
            } else if (grad[i][j] > T2) {
                grad[i][j] = 255;
            } else {
                grad[i][j] = 128;
            }
        }
    }
}

// 边缘跟踪
void edge_tracking(vector<vector<int>>& grad, vector<vector<int>>& edges) {
    int rows = grad.size();
    int cols = grad[0].size();
    int dx[8] = {-1, 0, 1, -1, 1, -1, 0, 1};
    int dy[8] = {-1, -1, -1, 0, 0, 1, 1, 1};
    for (int i = 1; i < rows - 1; i++) {
        for (int j = 1; j < cols - 1; j++) {
            if (grad[i][j] == 255) {
                edges[i][j] = grad[i][j];
                for (int k = 0; k < 8; k++) {
                    int x = i + dx[k];
                    int y = j + dy[k];
                    if (grad[x][y] == 128) {
                        edges[x][y] = 255;
                    }
                }
            }
        }
    }
}

// 边缘检测
void detect_edges(vector<vector<int>>& image, vector<vector<int>>& edges, int T1, int T2) {
    int rows = image.size();
    int cols = image[0].size();
    vector<vector<int>> grad(rows, vector<int>(cols));
    vector<vector<double>> theta(rows, vector<double>(cols));
    // Sobel算子
    for (int i = 1; i < rows - 1; i++) {
        for (int j = 1; j < cols - 1; j++) {
            int gx = 0, gy = 0;
            for (int k = -1; k <= 1; k++) {
                for (int l = -1; l <= 1; l++) {
                    gx += image[i + k][j + l] * sobel_x[k + 1][l + 1];
                    gy += image[i + k][j + l] * sobel_y[k + 1][l + 1];
                }
            }
            grad[i][j] = sqrt(gx * gx + gy * gy);
            theta[i][j] = atan2(gy, gx);
        }
    }
    // 非极大值抑制
    non_maximum_suppression(grad, theta);
    // 双阈值处理
    double_thresholding(grad, T1, T2);
    // 边缘跟踪
    edge_tracking(grad, edges);
}

int main() {
    int T1 = 50; // 低阈值
    int T2 = 100; // 高阈值
    vector<vector<int>> image = {
        {0, 0, 0, 0, 0, 0, 0, 0},
        {0, 100, 100, 100, 100, 100, 100, 0},
        {0, 100, 0, 0, 0, 0, 100, 0},
        {0, 100, 0, 100, 100, 0, 100, 0},
        {0, 100, 0, 100, 100, 0, 100, 0},
        {0, 100, 0, 0, 0, 0, 100, 0},
        {0, 100, 100, 100, 100, 100, 100, 0},
        {0, 0, 0, 0, 0, 0, 0, 0}
    };
    int rows = image.size();
    int cols = image[0].size();
    vector<vector<int>> edges(rows, vector<int>(cols));
    detect_edges(image, edges, T1, T2);
    // 输出结果
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            cout << edges[i][j] << " ";
        }
        cout << endl;
    }
    return 0;
}