机器视觉原生算法入门

本文以边缘检测为例，提供一个机器视觉原生算法的入门案例。效果如下图，左边是源图片，右边是检测结果：

      

基本思路：逐行扫描探测边缘；对每行的边缘位置做直线拟合；使用方差剔除缺损点。

边缘探测算法：使用类似【2, 2, 2, -2, -2, -2】的卷积算子，在边缘位置，卷积值最大，利用最大值两边的次大值作比例均衡，提高精度。

废话不多说了，上代码：

class CDetectVHLine
{
public:
    CDetectVHLine(BYTE* image, int width, int height, int stride,
        int left, int right, int bottom, int top,
        int halfEdge, bool horizontal, bool blackToWhite)
    {
        m_image = image;
        m_width = width;
        m_height = height;
        m_stride = stride;

        m_left = left;
        m_right = right;
        m_bottom = bottom;
        m_top = top;

        m_halfEdge = halfEdge;
        m_horizontal = horizontal;
        m_blackToWhite = blackToWhite;

        // 生成检测算子
        m_factor.resize(m_halfEdge);
        for (int i = 0; i < m_halfEdge; ++i){
            m_factor[i] = 2;
        }

        // 卷积区间要去除边缘
        if (m_horizontal) {
            m_bottom += m_halfEdge;
            m_top -= m_halfEdge;
        }
        else {
            m_left += m_halfEdge;
            m_right -= m_halfEdge;
        }
    }
    ~CDetectVHLine(){}

    bool Detect()
    {
        if (m_horizontal){
            vector<double>    edges;                    // 保存计算结果，每个扫描线的边缘Y坐标
            edges.resize(m_right - m_left + 1);

            // 逐列扫描边缘
            for (int i = m_left, idx = 0; i <= m_right; ++i, idx++){
                edges[idx] = scanLineV(i, m_bottom, m_top);
            }
            if (!fitLine(m_left, edges, m_k, m_b)){
                return false;
            }
            m_x1 = m_left;
            m_y1 = m_k * m_x1 + m_b;
            m_x2 = m_right + 1;
            m_y2 = m_k * m_x2 + m_b;
        }
        else{
            vector<double> edges;
            edges.resize(m_top - m_bottom + 1);

            // 逐行扫描边缘
            for (int i = m_bottom, idx = 0; i <= m_top; ++i, idx++){
                edges[idx] = scanLineH(i, m_left, m_right);
            }
            if (!fitLine(m_bottom, edges, m_k, m_b)){
                return false;
            }
            m_y1 = m_bottom;
            m_x1 = m_k * m_y1 + m_b;
            m_y2 = m_top + 1;
            m_x2 = m_k* m_y2 + m_b;
        }

        return true;
    }

    // 检测结果
    double m_x1, m_y1, m_x2, m_y2;
    double m_k, m_b;

private:
    BYTE*    m_image;
    int        m_width, m_height;                            // 图像尺寸
    int        m_stride;                                    // 扫描行字节数
    int        m_left, m_right, m_bottom, m_top;            // 探测区域
    int        m_halfEdge;                                    // 边缘宽度的一半：像素
    bool    m_horizontal;                                // true：水平线；false：竖直线
    bool    m_blackToWhite;                                // true：左黑右白/下黑上白；false：左白右黑/下白上黑
    vector<int>    m_factor;                                // 对称边缘检测算子正半部分

    const    int    c_productThresh = 100;                    // 边缘卷积阈值，小于该值则不认为是边缘

    // 水平方向扫描检测边缘
    double scanLineH(int y, int left, int right)
    {
        int maxL = 0, maxM = 0, maxR = 0, maxX;
        int currL = 0, currM = 0, currR = 0;
        BYTE* pData = m_image + m_stride * y + left;
        for (int i = left; i <= right; ++i){
            currR = m_blackToWhite ? sumProductH(pData) : -sumProductH(pData);
            if (currM >= currL && currM >= currR){            // 判定极值点
                if (currM > maxM){
                    maxL = currL; maxM = currM; maxR = currR; maxX = i - 1;
                }
            }
            currL = currM; currM = currR;
            pData++;
        }

        if (maxM < c_productThresh){
            return -1;
        }

        return maxX + (double)(maxM - maxL) / (maxM - maxL + maxM - maxR) - 0.5;
    }

    // 垂直方向扫描检测边缘
    double scanLineV(int x, int bottom, int top)
    {
        int maxB = 0, maxM = 0, maxT = 0, maxY;
        int currB = 0, currM = 0, currT = 0;
        BYTE* pData = m_image + m_stride * bottom + x;
        for (int i = bottom; i <= top; ++i){
            currT = m_blackToWhite ? sumProductV(pData) : -sumProductV(pData);
            if (currM >= currB && currM >= currT){            // 判定极值点
                if (currM > maxM){
                    maxB= currB; maxM = currM; maxT = currT; maxY = i - 1;
                }
            }
            currB = currM; currM = currT;
            pData += m_stride;
        }

        if (maxM < c_productThresh){
            return -1;
        }

        return maxY + (double)(maxM - maxB) / (maxM - maxB + maxM - maxT) - 0.5;
    }

    // 计算水平扫描卷积值，pData对应算子中心
    int    sumProductH(const BYTE* pData)
    {
        int product = 0;
        const BYTE* pDataN = pData - 1;
        for (size_t i = 0; i < m_factor.size(); ++i){
            product += (m_factor[i] * pData[0] - m_factor[i] * pDataN[0]);
            pData++;
            pDataN--;
        }

        return product;
    }

    // 计算垂直扫描卷积值，pData对应算子中心
    int    sumProductV(const BYTE* pData)
    {
        int product = 0;
        const BYTE* pDataN = pData - m_stride;
        for (size_t i = 0; i < m_factor.size(); ++i){
            product += (m_factor[i] * pData[0] - m_factor[i] * pDataN[0]);
            pData += m_stride;
            pDataN -= m_stride;
        }

        return product;
    }

    // 直线拟合
    bool fitLine(int start, vector<double> &points, double &k, double &b)
    {
        // 统计有效点数目
        int validCount = 0;
        for (size_t i = 0; i < points.size(); ++i){
            if (points[i] > 0){
                validCount++;
            }
        }
        if (validCount < 2)
            return false;

        // 根据方差迭代
        const int iterCount = 3;    // 迭代次数
        for (int iterator = 0; iterator < iterCount; ++iterator){
            double sigmaX = 0, sigmaY = 0, sigmaX2 = 0, sigmaXY = 0;
            double X = start - 0.5;
            for (size_t i = 0; i < points.size(); ++i){
                X += 1;
                if (points[i] < 0){
                    continue;
                }
                sigmaX += X;
                sigmaY += points[i];
                sigmaXY += X * points[i];
                sigmaX2 += X * X;
            }
            double denominator = (validCount * sigmaX2 - sigmaX * sigmaX);
            k = (validCount * sigmaXY - sigmaX * sigmaY) / denominator;
            b = (sigmaX2 * sigmaY - sigmaX * sigmaXY) / denominator;

            // 方差
            double e = 0;
            X = start - 0.5;
            for (size_t i = 0; i < points.size(); ++i){
                X += 1;
                if (points[i] < 0){
                    continue;
                }
                e += (points[i] - k * X - b) * (points[i] - k * X - b);
            }
            e /= validCount;
            if (e < 2){
                break;
            }
            e = sqrt(e);

            // 剔除误差过大的点
            X = start - 0.5;
            for (size_t i = 0; i < points.size(); ++i){
                X += 1;
                if (points[i] < 0){
                    continue;
                }
                if (abs(points[i] - k * X - b) > e){
                    points[i] = -1;
                    validCount--;
                }
            }
            if (validCount < 2){
                break;
            }
        }

        return true;
    }
};