織布缺陷——斷針缺陷檢測

本文針對織布生產過程中，由於斷針造成的織布缺陷圖像，進行檢測,如下圖1所示.

    

1. 圖像質量增強

　　由於織布原圖的圖像整體偏暗，不利於缺陷部分的檢測，考慮改善圖像質量，而圖像質量增強的算法有很多，本文借鑒Retinex的圖像增強算法(SSR, MSR, MSRCR)的實現，效果如圖2所示.

 

 

參考:https://blog.csdn.net/ajianyingxiaoqinghan/article/details/71435098

2. 低通濾波/邊緣檢測

 觀測圖像2中增強后的圖像，雖然織布瑕疵部分，相對正常區域有明顯的細長暗條紋痕跡，但是如果直接使用邊緣檢測的方式，來定位瑕疵部分的邊緣，很難找到明顯的邊界條紋（受織布正常紋理部分的影響,如圖3的對比效果圖所示）.

 這里用低通濾波器，實現對織布紋理背景的過濾，以提供織布瑕疵缺陷的高頻細節部分，利於后續的邊緣檢測效果.

   低通濾波算法，參考：https://blog.csdn.net/weixin_40647819/article/details/80600918?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.channel_param&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-2.channel_param

  邊緣檢測算法，參考：https://blog.csdn.net/dieju8330/article/details/82814529?utm_medium=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.channel_param&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromMachineLearnPai2-1.channel_param

3. 形態學處理

對邊緣檢測后的圖像，進行形態學的膨脹處理，就可以增大缺陷區域的輪廓，以便於利用opencv中找輪廓的方式，定位瑕疵部分的位置.

4. 定位瑕疵區域

     對形態學處理后的圖像，利用opencv中的findContours()函數,通過輪廓的寬高、面積大小等條件設置，來定位織布的瑕疵輪廓區域，具體實現參見代碼部分：

5. 實現代碼部分

MSRCR.h

#pragma once

#ifndef _MSRCR_H_
#define _MSRCR_H_

#include "cv.h"
#include "highgui.h"
#include "opencv2\opencv.hpp"
#include <math.h>

// 該處使用USE_EXACT_SIGMA，則使用自定義的濾波算法；
// 不使用USE_EXACT_SIGMA，則使用OpenCV自帶的高斯濾波算法；
//#define USE_EXACT_SIGMA

using namespace cv;
using namespace std;

class Msrcr
{
private:
#define pc(image, x, y, c) image->imageData[(image->widthStep * y) + (image->nChannels * x) + c]
#define INT_PREC 1024.0
#define INT_PREC_BITS 10
    inline double int2double(int x) { return (double)x / INT_PREC; }
    inline int double2int(double x) { return (int)(x * INT_PREC + 0.5); }
    inline int int2smallint(int x) { return (x >> INT_PREC_BITS); }
    inline int int2bigint(int x) { return (x << INT_PREC_BITS); }
public:
    vector<double> CreateKernel(double sigma);
    vector<int> CreateFastKernel(double sigma);
    void FilterGaussian(IplImage* img, double sigma);
    void FilterGaussian(Mat src, Mat &dst, double sigma);
    void FastFilter(IplImage *img, double sigma);
    void FastFilter(Mat src, Mat &dst, double sigma);
    void Retinex(IplImage *img, double sigma, int gain = 128, int offset = 128);
    void Retinex(Mat src, Mat &dst, double sigma, int gain = 128, int offset = 128);
    void MultiScaleRetinex(IplImage *img, vector<double> weights, vector<double> sigmas, int gain = 128, int offset = 128);
    void MultiScaleRetinex(Mat src, Mat &dst, vector<double> weights, vector<double> sigmas, int gain = 128, int offset = 128);
    void MultiScaleRetinexCR(IplImage *img, vector<double> weights, vector<double> sigmas, int gain = 128, int offset = 128,
        double restoration_factor = 6, double color_gain = 2);
    void MultiScaleRetinexCR(Mat src, Mat &dst, vector<double> weights, vector<double> sigmas, int gain = 128, int offset = 128,
        double restoration_factor = 6, double color_gain = 2);
};

#endif

 

MSRCR.cpp

#include "MSRCR.h"

/*===========================================================
* 函數: CreateKernel
* 說明：創建一個標准化的一維高斯核；
* 參數：
*   double sigma: 高斯核標准偏差
* 返回值：
*   double*: 長度為((6*sigma)/2) * 2 + 1的double數組
* 注：
*   調用者需要刪除該內核；
* --------------------------------------------------------
*  Summary:
*  Creates a normalized 1 dimensional gaussian kernel.
*
*  Arguments:
*  double sigma - the standard deviation of the gaussian kernel.
*
*  Returns:
*  vector<double> - an vector of values of length ((6*sigma)/2) * 2 + 1.
*  Note:
*  Caller is responsable for deleting the kernel.
=============================================================
*/
vector<double> Msrcr::CreateKernel(double sigma)
{
    int i, x, filter_size;
    vector<double> filter;
    double sum;

    // set sigma's upline
    // 為sigma設定上限
    if (sigma > 300) sigma = 300;

    // get needed filter size (enforce oddness)
    // 獲取需要的濾波尺寸，且強制為奇數；
    filter_size = (int)floor(sigma * 6) / 2;
    filter_size = filter_size * 2 + 1;

    // Calculate exponential
    // 計算指數
    sum = 0;
    for (i = 0; i < filter_size; i++)
    {
        double tmpValue;
        x = i - (filter_size / 2);
        tmpValue = exp(-(x*x) / (2 * sigma*sigma));
        filter.push_back(tmpValue);

        sum += tmpValue;
    }

    // Normalize
    // 歸一化計算
    for (i = 0, x; i < filter_size; i++)
        filter[i] /= sum;

    return filter;
}

/*===========================================================
* 函數: CreateFastKernel
* 說明：創建一個近似浮點的整數類型（左移8bits）的快速高斯核；
* 參數：
*   double sigma: 高斯核標准偏差
* 返回值：
*   double*: 長度為((6*sigma)/2) * 2 + 1的int數組
* 注：
*   調用者需要刪除該內核；
* --------------------------------------------------------
* Summary:
* Creates a faster gaussian kernal using integers that
* approximate floating point (leftshifted by 8 bits)
*
* Arguments:
* double sigma - the standard deviation of the gaussian kernel.
*
* Returns:
* vector<int> - an vector of values of length ((6*sigma)/2) * 2 + 1.
*
* Note:
* Caller is responsable for deleting the kernel.
=============================================================
*/
vector<int> Msrcr::CreateFastKernel(double sigma)
{
    vector<double> fp_kernel;
    vector<int> kernel;
    int i, filter_size;

    // Reject unreasonable demands
    // 設置上限
    if (sigma > 300) sigma = 300;

    // get needed filter size (enforce oddness)
    // 獲取需要的濾波尺寸，且強制為奇數；
    filter_size = (int)floor(sigma * 6) / 2;
    filter_size = filter_size * 2 + 1;

    // Create Kernel
    // 創建內核
    fp_kernel = CreateKernel(sigma);

    // Change kernel's data type from double to int
    // double內核轉為int型
    for (i = 0; i < filter_size; i++)
    {
        int tmpValue;
        tmpValue = double2int(fp_kernel[i]);
        kernel.push_back(tmpValue);
    }

    return kernel;
}

/*===========================================================
* 函數：FilterGaussian
* 說明：通過內核計算高斯卷積，內核由sigma值得到，且在內核兩端值相等；
* 參數：
*   img: 被濾波的IplImage*類型圖像
*   double sigma: 高斯核標准偏差
* --------------------------------------------------------
* Summary:
* Performs a gaussian convolution for a value of sigma that is equal
* in both directions.
*
* Arguments:
* IplImage* img - the image to be filtered in place.
* double sigma - the standard deviation of the gaussian kernel to use.
=============================================================
*/
void Msrcr::FilterGaussian(IplImage* img, double sigma)
{
    int i, j, k, source, filter_size;
    vector<int> kernel;
    IplImage* temp;
    int v1, v2, v3;

    // Reject unreasonable demands
    // 設置上限
    if (sigma > 300) sigma = 300;

    // get needed filter size (enforce oddness)
    // 獲取需要的濾波尺寸，且強制為奇數；
    filter_size = (int)floor(sigma * 6) / 2;
    filter_size = filter_size * 2 + 1;

    // Create Kernel
    // 創建內核
    kernel = CreateFastKernel(sigma);

    temp = cvCreateImage(cvSize(img->width, img->height), img->depth, img->nChannels);

    // filter x axis
    // X軸濾波
    for (j = 0; j < temp->height; j++)
    {
        for (i = 0; i < temp->width; i++)
        {
            // inner loop has been unrolled
            // 內層循環已經展開
            v1 = v2 = v3 = 0;
            for (k = 0; k < filter_size; k++)
            {
                source = i + filter_size / 2 - k;

                if (source < 0) source *= -1;
                if (source > img->width - 1) source = 2 * (img->width - 1) - source;

                v1 += kernel[k] * (unsigned char)pc(img, source, j, 0);
                if (img->nChannels == 1) continue;
                v2 += kernel[k] * (unsigned char)pc(img, source, j, 1);
                v3 += kernel[k] * (unsigned char)pc(img, source, j, 2);
            }

            // set value and move on
            pc(temp, i, j, 0) = (char)int2smallint(v1);
            if (img->nChannels == 1) continue;
            pc(temp, i, j, 1) = (char)int2smallint(v2);
            pc(temp, i, j, 2) = (char)int2smallint(v3);

        }
    }

    // filter y axis
    // Y軸濾波
    for (j = 0; j < img->height; j++)
    {
        for (i = 0; i < img->width; i++)
        {
            v1 = v2 = v3 = 0;
            for (k = 0; k < filter_size; k++)
            {
                source = j + filter_size / 2 - k;

                if (source < 0) source *= -1;
                if (source > temp->height - 1) source = 2 * (temp->height - 1) - source;

                v1 += kernel[k] * (unsigned char)pc(temp, i, source, 0);
                if (img->nChannels == 1) continue;
                v2 += kernel[k] * (unsigned char)pc(temp, i, source, 1);
                v3 += kernel[k] * (unsigned char)pc(temp, i, source, 2);
            }

            // set value and move on
            pc(img, i, j, 0) = (char)int2smallint(v1);
            if (img->nChannels == 1) continue;
            pc(img, i, j, 1) = (char)int2smallint(v2);
            pc(img, i, j, 2) = (char)int2smallint(v3);

        }
    }

    cvReleaseImage(&temp);
}

/*===========================================================
* 函數：FilterGaussian
* 說明：通過內核計算高斯卷積，內核由sigma值得到，且在內核兩端值相等；
* 參數：
*   Mat src: 輸入圖像
*   Mat &dst: 輸出圖像
*   double sigma: 高斯核標准偏差
* --------------------------------------------------------
* Summary:
* Performs a gaussian convolution for a value of sigma that is equal
* in both directions.
*
* Arguments:
* Mat src - Input Image.
* Mat &dst - Output Image.
* double sigma - the standard deviation of the gaussian kernel to use.
=============================================================
*/
void Msrcr::FilterGaussian(Mat src, Mat &dst, double sigma)
{
    IplImage tmp_ipl;
    tmp_ipl = IplImage(src);
    FilterGaussian(&tmp_ipl, sigma);
    dst = cvarrToMat(&tmp_ipl);

    //dst = Mat(&tmp_ipl);
}

/*===========================================================
* 函數：FastFilter
* 說明：給出任意大小的sigma值，都可以通過使用圖像金字塔與可分離濾波器計算高斯卷積；
* 參數：
*   IplImage *img: 被濾波的圖像
*   double sigma: 高斯核標准偏差
* --------------------------------------------------------
* Summary:
* Performs gaussian convolution of any size sigma very fast by using
* both image pyramids and seperable filters.  Recursion is used.
*
* Arguments:
* img - an IplImage to be filtered in place.
* double sigma - the standard deviation of the gaussian kernel to use.
=============================================================
*/
void Msrcr::FastFilter(IplImage *img, double sigma)
{
    int filter_size;

    // Reject unreasonable demands
    // 設置上限
    if (sigma > 300) sigma = 300;

    // get needed filter size (enforce oddness)
    // 獲取需要的濾波尺寸，且強制為奇數；
    filter_size = (int)floor(sigma * 6) / 2;
    filter_size = filter_size * 2 + 1;

    // If 3 sigma is less than a pixel, why bother (ie sigma < 2/3)
    // 如果3 * sigma小於一個像素，則直接退出
    if (filter_size < 3) return;

    // Filter, or downsample and recurse
    // 處理方式：(1) 濾波  (2) 高斯光滑處理  (3) 遞歸處理濾波器大小
    if (filter_size < 10) {

#ifdef USE_EXACT_SIGMA
        FilterGaussian(img, sigma);
#else
        cvSmooth(img, img, CV_GAUSSIAN, filter_size, filter_size);
#endif

    }
    else
    {
        if (img->width < 2 || img->height < 2) return;
        IplImage* sub_img = cvCreateImage(cvSize(img->width / 2, img->height / 2), img->depth, img->nChannels);
        cvPyrDown(img, sub_img);
        FastFilter(sub_img, sigma / 2.0);
        cvResize(sub_img, img, CV_INTER_LINEAR);
        cvReleaseImage(&sub_img);
    }
}

/*===========================================================
* 函數：FastFilter
* 說明：給出任意大小的sigma值，都可以通過使用圖像金字塔與可分離濾波器計算高斯卷積；
* 參數：
*   Mat src: 輸入圖像
*   Mat &dst: 輸出圖像
*   double sigma: 高斯核標准偏差
* --------------------------------------------------------
* Summary:
* Performs gaussian convolution of any size sigma very fast by using
* both image pyramids and seperable filters.  Recursion is used.
*
* Arguments:
* Mat src - Input Image.
* Mat &dst - Output Image.
* double sigma - the standard deviation of the gaussian kernel to use.
=============================================================
*/
void Msrcr::FastFilter(Mat src, Mat &dst, double sigma)
{
    IplImage tmp_ipl;
    tmp_ipl = IplImage(src);
    FastFilter(&tmp_ipl, sigma);
    dst = cvarrToMat(&tmp_ipl);
    //dst = Mat(&tmp_ipl);
}

/*===========================================================
* 函數：Retinex
* 說明：單通道SSR方法，基礎Retinex復原算法。原圖像和被濾波的圖像需要被轉換到
*   對數域，並做減運算；
* 參數：
*   IplImage *img: 被濾波的圖像
*   double sigma: 高斯核標准偏差
*   int gain: 圖像像素值改變范圍的增益
*   int offset: 圖像像素值改變范圍的偏移量
* --------------------------------------------------------
* Summary:
* Basic retinex restoration. The image and a filtered image are converted
* to the log domain and subtracted.
*
* Arguments:
* img - an IplImage to be enhanced in place.
* sigma - the standard deviation of the gaussian kernal used to filter.
* gain - the factor by which to scale the image back into visable range.
* offset - an offset similar to the gain.
=============================================================
*/
void Msrcr::Retinex(IplImage *img, double sigma, int gain, int offset)
{
    IplImage *A, *fA, *fB, *fC;

    // Initialize temp images
    // 初始化緩存圖像
    fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
    fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
    fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);

    // Compute log image
    // 計算對數圖像
    cvConvert(img, fA);
    cvLog(fA, fB);

    // Compute log of blured image
    // 計算濾波后模糊圖像的對數圖像
    A = cvCloneImage(img);
    FastFilter(A, sigma);
    cvConvert(A, fA);
    cvLog(fA, fC);

    // Compute difference
    // 計算兩圖像之差
    cvSub(fB, fC, fA);

    // Restore
    // 恢復圖像
    cvConvertScale(fA, img, gain, offset);

    // Release temp images
    // 釋放緩存圖像
    cvReleaseImage(&A);
    cvReleaseImage(&fA);
    cvReleaseImage(&fB);
    cvReleaseImage(&fC);

}

/*===========================================================
* 函數：Retinex
* 說明：單通道SSR方法，基礎Retinex復原算法。原圖像和被濾波的圖像需要被轉換到
*   對數域，並做減運算；
* 參數：
*   Mat src: 輸入圖像
*   Mat &dst: 輸出圖像
*   double sigma: 高斯核標准偏差
*   int gain: 圖像像素值改變范圍的增益
*   int offset: 圖像像素值改變范圍的偏移量
* --------------------------------------------------------
* Summary:
* Basic retinex restoration. The image and a filtered image are converted
* to the log domain and subtracted.
*
* Arguments:
* Mat src - Input Image.
* Mat &dst - Output Image.
* sigma - the standard deviation of the gaussian kernal used to filter.
* gain - the factor by which to scale the image back into visable range.
* offset - an offset similar to the gain.
=============================================================
*/
void Msrcr::Retinex(Mat src, Mat &dst, double sigma, int gain, int offset)
{
    IplImage tmp_ipl;
    tmp_ipl = IplImage(src);
    Retinex(&tmp_ipl, sigma, gain, offset);

    dst = cvarrToMat(&tmp_ipl);
    //dst = Mat(&tmp_ipl);
}

/*===========================================================
* 函數：MultiScaleRetinex
* 說明：多通道MSR算法。原圖像和一系列被濾波的圖像轉換到對數域，並與帶權重的原圖像做減運算。
* 通常情況下，三個權重范圍選擇低、中、高標准偏差；
*
* 參數：
*   IplImage *img: 被濾波的圖像
*   vector<double> weights: 通道權重
*   vector<double> sigmas: 高斯核標准偏差
*   int gain: 圖像像素值改變范圍的增益
*   int offset: 圖像像素值改變范圍的偏移量
* --------------------------------------------------------
* Summary:
* Multiscale retinex restoration.  The image and a set of filtered images are
* converted to the log domain and subtracted from the original with some set
* of weights. Typicaly called with three equaly weighted scales of fine,
* medium and wide standard deviations.
*
* Arguments:
* IplImage* img - an IplImage to be enhanced in place.
* vector<double> weights - Weights of channels
* vector<double> sigma - the standard deviation of the gaussian kernal used to filter.
* int gain - the factor by which to scale the image back into visable range.
* int offset - an offset similar to the gain.
=============================================================
*/
void Msrcr::MultiScaleRetinex(IplImage *img, vector<double> weights, vector<double> sigmas, int gain, int offset)
{
    int i;
    double weight;
    int scales = sigmas.size();
    IplImage *A, *fA, *fB, *fC;

    // Initialize temp images
    fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
    fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
    fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);


    // Compute log image
    cvConvert(img, fA);
    cvLog(fA, fB);

    // Normalize according to given weights
    for (i = 0, weight = 0; i < scales; i++)
        weight += weights[i];

    if (weight != 1.0) cvScale(fB, fB, weight);

    // Filter at each scale
    for (i = 0; i < scales; i++)
    {
        A = cvCloneImage(img);
        double tmp = sigmas[i];
        FastFilter(A, tmp);

        cvConvert(A, fA);
        cvLog(fA, fC);
        cvReleaseImage(&A);

        // Compute weighted difference
        cvScale(fC, fC, weights[i]);
        cvSub(fB, fC, fB);
    }

    // Restore
    cvConvertScale(fB, img, gain, offset);

    // Release temp images
    cvReleaseImage(&fA);
    cvReleaseImage(&fB);
    cvReleaseImage(&fC);
}

/*===========================================================
* 函數：MultiScaleRetinex
* 說明：多通道MSR算法。原圖像和一系列被濾波的圖像轉換到對數域，並與帶權重的原圖像做減運算。
* 通常情況下，三個權重范圍選擇低、中、高標准偏差；
*
* 參數：
*   Mat src: 輸入圖像
*   Mat &dst: 輸出圖像
*   vector<double> weights: 通道權重
*   vector<double> sigmas: 高斯核標准偏差
*   int gain: 圖像像素值改變范圍的增益
*   int offset: 圖像像素值改變范圍的偏移量
* --------------------------------------------------------
* Summary:
* Multiscale retinex restoration.  The image and a set of filtered images are
* converted to the log domain and subtracted from the original with some set
* of weights. Typicaly called with three equaly weighted scales of fine,
* medium and wide standard deviations.
*
* Arguments:
* Mat src - Input Image.
* Mat &dst - Output Image.
* vector<double> weights - Weights of channels
* vector<double> sigmas - the standard deviation of the gaussian kernal used to filter.
* int gain - the factor by which to scale the image back into visable range.
* int offset - an offset similar to the gain.
=============================================================
*/
void Msrcr::MultiScaleRetinex(Mat src, Mat &dst, vector<double> weights, vector<double> sigmas, int gain, int offset)
{
    IplImage tmp_ipl;
    tmp_ipl = IplImage(src);
    MultiScaleRetinex(&tmp_ipl, weights, sigmas, gain, offset);

    dst = cvarrToMat(&tmp_ipl);

    //dst = Mat(&tmp_ipl);
}

/*===========================================================
* 函數：MultiScaleRetinexCR
* 說明：MSRCR算法，MSR算法加上顏色修復。原圖像和一系列被濾波的圖像轉換到對數域，並與帶權重的原圖像做減運算。
* 通常情況下，三個權重范圍選擇低、中、高標准偏差；之后，顏色修復權重應用於每個顏色通道中；
*
* 參數：
*   IplImage *img: 被濾波的圖像
*   double sigma: 高斯核標准偏差
*   int gain: 圖像像素值改變范圍的增益
*   int offset: 圖像像素值改變范圍的偏移量
*   double restoration_factor: 控制顏色修復的非線性
*   double color_gain: 控制顏色修復增益
* --------------------------------------------------------
* Summary:
* Multiscale retinex restoration with color restoration.  The image and a set of
* filtered images are converted to the log domain and subtracted from the
* original with some set of weights. Typicaly called with three equaly weighted
* scales of fine, medium and wide standard deviations. A color restoration weight
* is then applied to each color channel.
*
* Arguments:
* IplImage *img - an IplImage to be enhanced in place.
* double sigma - the standard deviation of the gaussian kernal used to filter.
* int gain - the factor by which to scale the image back into visable range.
* int offset - an offset similar to the gain.
* double restoration_factor - controls the non-linearaty of the color restoration.
* double color_gain - controls the color restoration gain.
=============================================================
*/
void Msrcr::MultiScaleRetinexCR(IplImage *img, vector<double> weights, vector<double> sigmas,
    int gain, int offset, double restoration_factor, double color_gain)
{
    int i;
    double weight;
    int scales = sigmas.size();
    IplImage *A, *B, *C, *fA, *fB, *fC, *fsA, *fsB, *fsC, *fsD, *fsE, *fsF;

    // Initialize temp images
    // 初始化緩存圖像
    fA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
    fB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
    fC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, img->nChannels);
    fsA = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);
    fsB = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);
    fsC = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);
    fsD = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);
    fsE = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);
    fsF = cvCreateImage(cvSize(img->width, img->height), IPL_DEPTH_32F, 1);

    // Compute log image
    // 計算對數圖像
    cvConvert(img, fB);
    cvLog(fB, fA);

    // Normalize according to given weights
    // 依照權重歸一化
    for (i = 0, weight = 0; i < scales; i++)
        weight += weights[i];

    if (weight != 1.0) cvScale(fA, fA, weight);

    // Filter at each scale
    // 各尺度上進行濾波操作
    for (i = 0; i < scales; i++) {
        A = cvCloneImage(img);
        FastFilter(A, sigmas[i]);

        cvConvert(A, fB);
        cvLog(fB, fC);
        cvReleaseImage(&A);

        // Compute weighted difference
        // 計算權重后兩圖像之差
        cvScale(fC, fC, weights[i]);
        cvSub(fA, fC, fA);
    }

    // Color restoration
    // 顏色修復
    if (img->nChannels > 1) {
        A = cvCreateImage(cvSize(img->width, img->height), img->depth, 1);
        B = cvCreateImage(cvSize(img->width, img->height), img->depth, 1);
        C = cvCreateImage(cvSize(img->width, img->height), img->depth, 1);

        // Divide image into channels, convert and store sum
        // 將圖像分割為若干通道，類型轉換為浮點型，並存儲通道數據之和
        cvSplit(img, A, B, C, NULL);
        cvConvert(A, fsA);
        cvConvert(B, fsB);
        cvConvert(C, fsC);

        cvReleaseImage(&A);
        cvReleaseImage(&B);
        cvReleaseImage(&C);

        // Sum components
        // 求和
        cvAdd(fsA, fsB, fsD);
        cvAdd(fsD, fsC, fsD);

        // Normalize weights
        // 帶權重矩陣歸一化
        cvDiv(fsA, fsD, fsA, restoration_factor);
        cvDiv(fsB, fsD, fsB, restoration_factor);
        cvDiv(fsC, fsD, fsC, restoration_factor);

        cvConvertScale(fsA, fsA, 1, 1);
        cvConvertScale(fsB, fsB, 1, 1);
        cvConvertScale(fsC, fsC, 1, 1);

        // Log weights
        // 帶權重矩陣求對數
        cvLog(fsA, fsA);
        cvLog(fsB, fsB);
        cvLog(fsC, fsC);

        // Divide retinex image, weight accordingly and recombine
        // 將Retinex圖像切分為三個數組，按照權重和顏色增益重新組合
        cvSplit(fA, fsD, fsE, fsF, NULL);

        cvMul(fsD, fsA, fsD, color_gain);
        cvMul(fsE, fsB, fsE, color_gain);
        cvMul(fsF, fsC, fsF, color_gain);

        cvMerge(fsD, fsE, fsF, NULL, fA);
    }

    // Restore
    // 恢復圖像
    cvConvertScale(fA, img, gain, offset);

    // Release temp images
    // 釋放緩存圖像
    cvReleaseImage(&fA);
    cvReleaseImage(&fB);
    cvReleaseImage(&fC);
    cvReleaseImage(&fsA);
    cvReleaseImage(&fsB);
    cvReleaseImage(&fsC);
    cvReleaseImage(&fsD);
    cvReleaseImage(&fsE);
    cvReleaseImage(&fsF);
}

/*===========================================================
* 函數：MultiScaleRetinexCR
* 說明：MSRCR算法，MSR算法加上顏色修復。原圖像和一系列被濾波的圖像轉換到對數域，並與帶權重的原圖像做減運算。
* 通常情況下，三個權重范圍選擇低、中、高標准偏差；之后，顏色修復權重應用於每個顏色通道中；
*
* 參數：
*   Mat src: 輸入圖像
*   Mat &dst: 輸出圖像
*   double sigma: 高斯核標准偏差
*   int gain: 圖像像素值改變范圍的增益
*   int offset: 圖像像素值改變范圍的偏移量
*   double restoration_factor: 控制顏色修復的非線性
*   double color_gain: 控制顏色修復增益
* --------------------------------------------------------
* Summary:
* Multiscale retinex restoration with color restoration.  The image and a set of
* filtered images are converted to the log domain and subtracted from the
* original with some set of weights. Typicaly called with three equaly weighted
* scales of fine, medium and wide standard deviations. A color restoration weight
* is then applied to each color channel.
*
* Arguments:
* Mat src - Input Image.
* Mat &dst - Output Image.
* double sigma - the standard deviation of the gaussian kernal used to filter.
* int gain - the factor by which to scale the image back into visable range.
* int offset - an offset similar to the gain.
* double restoration_factor - controls the non-linearaty of the color restoration.
* double color_gain - controls the color restoration gain.
=============================================================
*/
void Msrcr::MultiScaleRetinexCR(Mat src, Mat &dst, vector<double> weights, vector<double> sigmas,
    int gain, int offset, double restoration_factor, double color_gain)
{
    IplImage tmp_ipl;
    tmp_ipl = IplImage(src);
    MultiScaleRetinexCR(&tmp_ipl, weights, sigmas, gain, offset, restoration_factor, color_gain);

    dst = cvarrToMat(&tmp_ipl);

    //dst = Mat(&tmp_ipl);
}

 

 Main.cpp

#include "cv.h"
#include "highgui.h"
#include "opencv2\opencv.hpp"

#include "MSRCR.h"

using namespace cv;
using namespace std;

//*****************高斯低通濾波器***********************
Mat gaussianlbrf(Mat scr, float sigma)
{
    Mat gaussianBlur(scr.size(), CV_32FC1); //，CV_32FC1
    float d0 = 2 * sigma*sigma;//高斯函數參數，越小，頻率高斯濾波器越窄，濾除高頻成分越多，圖像就越平滑
    for (int i = 0; i<scr.rows; i++)
    {
        for (int j = 0; j<scr.cols; j++)
        {
            float d = pow(float(i - scr.rows / 2), 2) + pow(float(j - scr.cols / 2), 2);//分子,計算pow必須為float型
            gaussianBlur.at<float>(i, j) = expf(-d / d0);//expf為以e為底求冪（必須為float型）
        }
    }
    //imshow("高斯低通濾波器", gaussianBlur);
    return gaussianBlur;
}

//*****************頻率域濾波*******************
Mat freqfilt(Mat scr, Mat blur)
{
    //***********************DFT*******************
    Mat plane[] = { scr, Mat::zeros(scr.size() , CV_32FC1) }; //創建通道，存儲dft后的實部與虛部（CV_32F，必須為單通道數）
    Mat complexIm;
    merge(plane, 2, complexIm);//合並通道 （把兩個矩陣合並為一個2通道的Mat類容器）
    dft(complexIm, complexIm);//進行傅立葉變換，結果保存在自身  

                              //***************中心化********************
    split(complexIm, plane);//分離通道（數組分離）

    int cx = plane[0].cols / 2; int cy = plane[0].rows / 2;//以下的操作是移動圖像  (零頻移到中心)
    Mat part1_r(plane[0], Rect(0, 0, cx, cy));  //元素坐標表示為(cx,cy)
    Mat part2_r(plane[0], Rect(cx, 0, cx, cy));
    Mat part3_r(plane[0], Rect(0, cy, cx, cy));
    Mat part4_r(plane[0], Rect(cx, cy, cx, cy));

    Mat temp;
    part1_r.copyTo(temp);  //左上與右下交換位置(實部)
    part4_r.copyTo(part1_r);
    temp.copyTo(part4_r);

    part2_r.copyTo(temp);  //右上與左下交換位置(實部)
    part3_r.copyTo(part2_r);
    temp.copyTo(part3_r);

    Mat part1_i(plane[1], Rect(0, 0, cx, cy));  //元素坐標(cx,cy)
    Mat part2_i(plane[1], Rect(cx, 0, cx, cy));
    Mat part3_i(plane[1], Rect(0, cy, cx, cy));
    Mat part4_i(plane[1], Rect(cx, cy, cx, cy));

    part1_i.copyTo(temp);  //左上與右下交換位置(虛部)
    part4_i.copyTo(part1_i);
    temp.copyTo(part4_i);

    part2_i.copyTo(temp);  //右上與左下交換位置(虛部)
    part3_i.copyTo(part2_i);
    temp.copyTo(part3_i);

    //*****************濾波器函數與DFT結果的乘積****************
    Mat blur_r, blur_i, BLUR;
    multiply(plane[0], blur, blur_r); //濾波（實部與濾波器模板對應元素相乘）
    multiply(plane[1], blur, blur_i);//濾波（虛部與濾波器模板對應元素相乘）
    Mat plane1[] = { blur_r, blur_i };
    merge(plane1, 2, BLUR);//實部與虛部合並

                           //*********************得到原圖頻譜圖***********************************
    magnitude(plane[0], plane[1], plane[0]);//獲取幅度圖像，0通道為實部通道，1為虛部，因為二維傅立葉變換結果是復數  
    plane[0] += Scalar::all(1);  //傅立葉變o換后的圖片不好分析，進行對數處理，結果比較好看 
    log(plane[0], plane[0]);    // float型的灰度空間為[0，1])
    normalize(plane[0], plane[0], 1, 0, CV_MINMAX);  //歸一化便於顯示
    //imshow("原圖像頻譜圖", plane[0]);


    //******************IDFT*******************************
    /*
    Mat part111(BLUR,Rect(0,0,cx,cy));  //元素坐標(cx,cy)
    Mat part222(BLUR,Rect(cx,0,cx,cy));
    Mat part333(BLUR,Rect(0,cy,cx,cy));
    Mat part444(BLUR,Rect(cx,cy,cx,cy));
    part111.copyTo(temp);  //左上與右下交換位置(虛部)
    part444.copyTo(part111);
    temp.copyTo(part444);

    part222.copyTo(temp);  //右上與左下交換位置
    part333.copyTo(part222);
    temp.copyTo(part333);
    */

    idft(BLUR, BLUR);    //idft結果也為復數
    split(BLUR, plane);//分離通道，主要獲取通道
    magnitude(plane[0], plane[1], plane[0]);  //求幅值(模)
    normalize(plane[0], plane[0], 1, 0, CV_MINMAX);  //歸一化便於顯示
    return plane[0];//返回參數
}

void img_processing(Msrcr &msrcr, Mat &img, Mat &outImg, vector<double> sigema, vector<double> weight)
{
    
    int edge = 10;
    cv::Rect ROI(edge, edge, img.cols - edge, img.rows - edge);

    // ************************1. 圖像增強
    Mat img_msrcr;
    msrcr.MultiScaleRetinexCR(img, img_msrcr, weight, sigema, 128, 128);

    // ********************** 2. 低通濾波圖
    cv::Mat img_gray;
    cv::cvtColor(img_msrcr, img_gray, cv::COLOR_BGR2GRAY);
    
    int w = getOptimalDFTSize(img_gray.cols);  //獲取進行dtf的最優尺寸
    int h = getOptimalDFTSize(img_gray.rows); //獲取進行dtf的最優尺寸

    Mat padded;
    copyMakeBorder(img_gray, padded, 0, h - img_gray.rows, 0, w - img_gray.cols, BORDER_CONSTANT, Scalar::all(0));  //邊界填充
    padded.convertTo(padded, CV_32FC1); //將圖像轉換為flaot型

    Mat gaussianKernel = gaussianlbrf(padded, 50);//高斯低通濾波器
    Mat out = freqfilt(padded, gaussianKernel);//頻率域濾波

    Mat imageF_8UC3;
    out.convertTo(imageF_8UC3, CV_8UC3, 255);


    // *********************3. Canny邊緣檢測
    Mat img_gaussian, img_canny;
    GaussianBlur(imageF_8UC3, img_gaussian, Size(5, 5), 0, 0, BORDER_DEFAULT); //高斯濾波

    int edgeThresh = 15;
    Canny(img_gaussian, img_canny, edgeThresh, edgeThresh*4, 3); //Canny檢測

    cv::Mat img_canny_cut = cv::Mat(img_canny, ROI);
    cv::Mat img_canny_roi = img_canny_cut.clone();


    // ***********************4. 形態學處理
    Mat element, img_dilate;
    element = getStructuringElement(MORPH_RECT, Size(5, 5));
    dilate(img_canny_roi, img_dilate, element, Point(-1, -1), 2, BORDER_CONSTANT);

    Mat img_bin;
    img_bin = img_dilate.clone();
    std::vector<std::vector<cv::Point> > contours;
    findContours(img_bin, contours, CV_RETR_EXTERNAL, CV_CHAIN_APPROX_NONE); //找輪廓,注意第5個參數為CV_RETR_EXTERNAL，只檢索外框  

    if (contours.size() <= 0 || contours.size() > 999999)
    {
        printf("obj_segmentSecond 函數：coutours的大小超出限制...");
        return;
    }

    int width_less = 0.2*w;

    double least_area = 1000;

    std::vector<cv::Rect> boundRect(contours.size());    // 繪制最小外接矩形集合
    std::vector<cv::RotatedRect> box(contours.size());    // 定義最小外接矩形集合

    Mat img_copy;
    img_copy = img.clone();

    cv::Mat img_copy_cut = cv::Mat(img_copy, ROI);
    cv::Mat img_copy_roi = img_copy_cut.clone();
    
    for (int i = 0; i < contours.size(); i++) 
    {
        double area = contourArea(contours[i]);

        if (area < least_area)
            continue;

        box[i] = cv::minAreaRect(cv::Mat(contours[i])); // 計算每個輪廓的最小外接矩形
        boundRect[i] = cv::boundingRect(cv::Mat(contours[i]));

        cv::Point2f rect[4];
        box[i].points(rect);  // 把最小外接矩形四個端點復制給rect數組

        //// 小於分塊的邊界長度

        // 確定豎_方向織布缺陷
        if (boundRect[i].width < boundRect[i].height) 
        {

            int box_X = boundRect[i].x;
            int box_Y = boundRect[i].y;

            // 判斷缺陷窗口的坐標X是否在邊界禁止范圍內
            if (box_X > img_copy_roi.cols - 50 || box_Y >img_copy_roi.rows - 15)
            {
                printf("坐標X在邊界禁止范圍內... \n");
                continue;
            }

            cv::circle(img_copy_roi, cv::Point(box[i].center.x, box[i].center.y),
                5, cv::Scalar(0, 255, 0), -1, 8);   // 繪制最小外接矩形的中心點
            /*cv::rectangle(img_copy, cv::Point(boundRect[i].x, boundRect[i].y),
            cv::Point(boundRect[i].x + boundRect[i].width, boundRect[i].y + boundRect[i].height), cv::Scalar(0, 255, 0), 2, 8);
            */

            for (int j = 0; j<4; j++)
            {
                cv::line(img_copy_roi, rect[j], rect[(j + 1) % 4], cv::Scalar(255), 2, 8);  //繪制最小外接矩形每條邊
                cv::line(img_copy_roi, rect[j], rect[(j + 1) % 4], cv::Scalar(0, 255, 0), 2, 8); //繪制最小外接矩形每條邊
            }
        }
        else  // 確定橫_方向織布缺陷
        {

            int box_X = boundRect[i].x;
            int box_Y = boundRect[i].y;

            // 判斷缺陷窗口的坐標X是否在邊界禁止范圍內
            if (box_X > img_copy_roi.cols - 50 || box_Y >img_copy_roi.rows - 15)
            {
                printf("坐標X在邊界禁止范圍內... \n");
                continue;
            }

            if (boundRect[i].width < 0.3*w)
                continue;

            cv::circle(img_copy_roi, cv::Point(box[i].center.x, box[i].center.y),
                5, cv::Scalar(0, 255, 0), -1, 8);   // 繪制最小外接矩形的中心點
            /*cv::rectangle(img_copy, cv::Point(boundRect[i].x, boundRect[i].y),
            cv::Point(boundRect[i].x + boundRect[i].width, boundRect[i].y + boundRect[i].height), cv::Scalar(0, 255, 0), 2, 8);
            */
            for (int j = 0; j<4; j++)
            {
                cv::line(img_copy_roi, rect[j], rect[(j + 1) % 4], cv::Scalar(255), 2, 8);  //繪制最小外接矩形每條邊
                cv::line(img_copy_roi, rect[j], rect[(j + 1) % 4], cv::Scalar(255, 0, 0), 2, 8); //繪制最小外接矩形每條邊
            }
        
        }
    }


    outImg = img_copy_roi.clone();

}


void read_video() 
{
    vector<double> sigema;
    vector<double> weight;

    for (int i = 0; i < 3; i++)
        weight.push_back(1. / 3);

    // 由於MSRCR.h中定義了宏USE_EXTRA_SIGMA，所以此處的vector<double> sigema並沒有什么意義
    sigema.push_back(30);
    sigema.push_back(150);
    sigema.push_back(300);

    Msrcr msrcr;


    // *****************************1. 測試單張圖像*****************************************
    //Mat img, dst;
    //img = imread("F:\\BaiduNetdiskDownload\\織物瑕疵檢測\\video_photo_fail\\570.jpg", 1);

    //img_processing(msrcr, img, dst, sigema, weight);

    //imshow("src", img);
    //imshow("dst", dst);
    //imwrite("F:\\BaiduNetdiskDownload\\織物瑕疵檢測\\frame_fail\\570_filter_50.jpg", dst);

    //waitKey(10000);
    //return;

    
    // *****************************2. 連續測試多幀圖像*****************************************
    string file_path, fileName, grayFile, saveFile_path, saveFileName;

    file_path = "F:\\BaiduNetdiskDownload\\織物瑕疵檢測\\video_photo_fail\\";
    saveFile_path = "F:\\BaiduNetdiskDownload\\織物瑕疵檢測\\frame_detect\\";

    string str1, str2;
    Mat srcImage, grayImage, outImg;
    for (int i = 1; i <= 100000; i++)
    {
        /*if (i < 2000)
            continue;*/

        printf("i =%d \n", i);

        stringstream ss1, ss2;

        ss1 << i;
        ss1 >> str1;

        //cout << str1 << endl;
        fileName = file_path + str1 + ".jpg";
        srcImage = imread(fileName, 1);
        
        img_processing(msrcr, srcImage, outImg, sigema, weight);
        saveFileName = saveFile_path + std::to_string(i) + ".bmp";

        imwrite(saveFileName, outImg);

        imshow("src", srcImage);
        imshow("dst", outImg);
    
        waitKey(1);

    }
    
}


int main()
{

    read_video();
    system("pause");
    
    return 0;
}