图像去色算法matlab到C++/opencv的移植

是对（Real-time Contrast Preserving Decolorization，作者 Cewu Lu Li Xu Jiaya Jia）论文配套的matlab代码一直到C++；并测试与其他算法的效果及时间对比

平台：qt5.12+MSVC2015（vs2015的编译器）+opencv4.2；

相关项目：http://cadik.posvete.cz/color_to_gray_overview/

测试图片为1200*768为28ms；512*512为21ms；因时间仓促对移植的代码没有优化，未来有时间可以优化代码，以减少运行时间；后边我会把核心代码粘贴上，希望感兴趣的一块交流学习。

测试：1、opencv自带的opencv自带的rgb2gray（左下方的灰度图时，时间a）；2、opencv中的decolor方法（右上方的灰度图，时间b）http://www.cse.cuhk.edu.hk/~leojia/projects/color2gray/。3、移植的灰度图（右下方的灰度图，时间c）；

测试结果：1、对于openv自带的rgb2gray（最常用的）速度最快但某些时候失真严重。因为它直接以Gray = R*0.299 + G*0.587 + B*0.114公式进行线性叠加导致灰度图像某些时候失真，如左下方的灰度图；

 2、opencv中的decolor方法是源于 香港中文大学 计算机科学与工程系 的论文 Contrast Preserving Decolorization（OpenCV 3.0以上版本）；效果好但速度最慢对于1200*768需要400ms（右上方的灰度）；

3、本文移植到opencv下灰度图转换，也是源于香港中文大学论文Real-time Contrast Preserving Decolorization

（chrome-extension://ibllepbpahcoppkjjllbabhnigcbffpi/http://www.cse.cuhk.edu.hk/leojia/papers/siga12t_color2gray.pdf）；原论文提供matlab代码；在下边会给出；他的测试结果（右下方的灰度图）时间c也能能够满足实时需求； 

测试demo链接 https://download.csdn.net/download/niuyuanye/12108794     

 

 

matlab  源码

function  img  = rtcprgb2gray(im)
   
  %%  Proprocessing 
 %im=imread('flower.jpg');
 [n,m,ch] = size(im); 
 sigma = 0.05;
 %W = wei();
 W = [    0         0    1.0000
         0    0.1000    0.9000
         0    0.2000    0.8000
         0    0.3000    0.7000
         0    0.4000    0.6000
         0    0.5000    0.5000
         0    0.6000    0.4000
         0    0.7000    0.3000
         0    0.8000    0.2000
         0    0.9000    0.1000
         0    1.0000         0
    0.1000         0    0.9000
    0.1000    0.1000    0.8000
    0.1000    0.2000    0.7000
    0.1000    0.3000    0.6000
    0.1000    0.4000    0.5000
    0.1000    0.5000    0.4000
    0.1000    0.6000    0.3000
    0.1000    0.7000    0.2000
    0.1000    0.8000    0.1000
    0.1000    0.9000         0
    0.2000         0    0.8000
    0.2000    0.1000    0.7000
    0.2000    0.2000    0.6000
    0.2000    0.3000    0.5000
    0.2000    0.4000    0.4000
    0.2000    0.5000    0.3000
    0.2000    0.6000    0.2000
    0.2000    0.7000    0.1000
    0.2000    0.8000         0
    0.3000         0    0.7000
    0.3000    0.1000    0.6000
    0.3000    0.2000    0.5000
    0.3000    0.3000    0.4000
    0.3000    0.4000    0.3000
    0.3000    0.5000    0.2000
    0.3000    0.6000    0.1000
    0.3000    0.7000    0.0000
    0.4000         0    0.6000
    0.4000    0.1000    0.5000
    0.4000    0.2000    0.4000
    0.4000    0.3000    0.3000
    0.4000    0.4000    0.2000
    0.4000    0.5000    0.1000
    0.4000    0.6000    0.0000
    0.5000         0    0.5000
    0.5000    0.1000    0.4000
    0.5000    0.2000    0.3000
    0.5000    0.3000    0.2000
    0.5000    0.4000    0.1000
    0.5000    0.5000         0
    0.6000         0    0.4000
    0.6000    0.1000    0.3000
    0.6000    0.2000    0.2000
    0.6000    0.3000    0.1000
    0.6000    0.4000    0.0000
    0.7000         0    0.3000
    0.7000    0.1000    0.2000
    0.7000    0.2000    0.1000
    0.7000    0.3000    0.0000
    0.8000         0    0.2000
    0.8000    0.1000    0.1000
    0.8000    0.2000    0.0000
    0.9000         0    0.1000
    0.9000    0.1000    0.0000
    1.0000         0         0];
 
   
  
%%  Global and Local Contrast Computing
ims = imresize(im, round(64/sqrt(n*m)*[n,m]),'nearest');
ims=im2double(ims)
R = ims(:,:,1);G = ims(:,:,2);B = ims(:,:,3);
imV = [R(:),G(:),B(:)];
%defaultStream = RandStream.getDefaultStream; 
%savedState = defaultStream.State;
t1 = randperm(size(imV,1));
Pg = [imV - imV(t1,:)];
 
ims = imresize(ims, 0.5 ,'nearest');
Rx = ims(:,1:end-1,1) - ims(:,2:end,1);
Gx = ims(:,1:end-1,2) - ims(:,2:end,2);
Bx = ims(:,1:end-1,3) - ims(:,2:end,3);
 
Ry = ims(1:end-1,:,1) - ims(2:end,:,1);
Gy = ims(1:end-1,:,2) - ims(2:end,:,2);
By = ims(1:end-1,:,3) - ims(2:end,:,3);
Pl = [[Rx(:),Gx(:),Bx(:)];[Ry(:),Gy(:),By(:)]];

P = [Pg;Pl ];
 

det = sqrt(sum(P.^2,2))/1.41;
 
% P( (det < 0.05),:) = []; 
% det( (det < 0.05)) = [];

detM =  repmat(det,[1,size(W,1)]);
L = P*W';
%% Energy optimization 
%bbbb=exp(- ( detM ).^2/sigma.^2);
U = log(exp(- (L + detM ).^2/sigma.^2) + exp(- (L- detM).^2/sigma.^2));
Es = mean(U);
%% Output
[NULLval,bw] = max(Es); 
img = imlincomb(W(bw,1),im(:,:,1) , W(bw,2),im(:,:,2) ,  W(bw,3),im(:,:,3));
%imshow(img);

end

c++ opencv 移植源码 

cv::Mat MainWindow::myRgb2Gray(cv::Mat image)
{
    cv::Mat IMAGE=image.clone();
    cv::Mat img_src=image.clone();

     double star=(double)cv::getTickCount();
     double scale=64/qSqrt(img_src.rows*img_src.cols);
     int h=round(scale*img_src.rows);
     int w=round(scale*img_src.cols);

     cv::Mat ims;
     cv::resize(img_src, img_src,cv::Size(w,h),cv::INTER_LINEAR);
     cv::cvtColor(img_src,img_src,cv::COLOR_BGR2RGB);

     img_src.convertTo(img_src,CV_64FC3);
     cv::normalize(img_src, ims, 0, 1, cv::NORM_MINMAX);




     std::vector<cv::Mat> vResultImage(3);
     std::vector<cv::Mat> vSrcImage;
    cv::split(ims, vSrcImage);

    for(int i=0;i<3;i++)
    {
//        Mat dst = data.reshape(0, 1);      //序列成行向量
//        Mat dst = data.reshape(0, 1).t();  //序列成列向量
        vResultImage[i]=vSrcImage[i].reshape(0,1).t();
    }
    cv::Mat img_re;
    cv::merge(vResultImage,img_re);
    cv::Mat img_Pg=cv::Mat::eye(img_re.rows,3,CV_64FC1);
    cv::Mat deM_Pg=cv::Mat::eye(img_re.rows,1,CV_64FC1);

    numbersList.clear();
    randData(ims.rows*ims.cols);
    for(int h=0;h<img_re.rows;h++)
    {
        double tmp=0;
        double tmp_dem=0;
        for(int c=0;c<3;c++)
        {
            tmp=img_re.at<cv::Vec3d>(h,0)[c]-img_re.at<cv::Vec3d>(numbersList[h],0)[c];
            img_Pg.at<double>(h,c)=tmp;
            tmp_dem+=tmp*tmp;
        }
        deM_Pg.at<double>(h,0)=qSqrt(tmp_dem)/1.41;
        //qDebug()<<("tmp:")<<(QString::number(h))<<(QString::number(tmp));
    }

    cv::resize(ims, ims,cv::Size(round(0.5*ims.cols), round(0.5*ims.rows)),cv::INTER_LINEAR);

    int num_xy=(ims.rows-1)*(ims.cols-1);
    cv::Mat ims_xy=cv::Mat::zeros(num_xy*2,3,CV_64FC1);
    cv::Mat deM_xy=cv::Mat::zeros(num_xy*2,1,CV_64FC1);
    //ims_x
    for(int h=0;h<ims.rows-1;h++)
    {
        int hh=h*(ims.cols-1);
        for(int w=0;w<ims.cols-1;w++)
        {
            double tmp=0;
            double tmp_x=0;
            double tmp_y=0;

            for(int c=0;c<3;c++)
            {
                tmp=ims.at<cv::Vec3d>(h,w)[c]-ims.at<cv::Vec3d>(h,w+1)[c];
                ims_xy.at<double>(hh+w,c)=tmp;
                tmp_x+=tmp*tmp;

                tmp=ims.at<cv::Vec3d>(h,w)[c]-ims.at<cv::Vec3d>(h+1,w)[c];
                ims_xy.at<double>(num_xy-1+hh+w,c)=tmp;
                tmp_y+=tmp*tmp;
            }
            deM_xy.at<double>(hh+w,0)=qSqrt(tmp_x)/1.41;
            deM_xy.at<double>(num_xy-1+hh+w,0)=qSqrt(tmp_y)/1.41;
        }
    }
    cv::Mat ims_P;
    cv::Mat ims_det;

    //图像竖向连接
    cv::vconcat(img_Pg, ims_xy,ims_P);
    ims_xy.release();
    cv::vconcat(deM_Pg, deM_xy,ims_det);
    deM_xy.release();

    cv::Mat ims_detM=ims_det.clone();
\
    cv::Mat ims_L;
    ims_L=ims_P*W;
    cv::Mat ims_U=cv::Mat::eye(ims_L.rows,ims_L.cols,CV_64FC1);
    for(int h=0;h<ims_detM.rows;h++)
    {
        double tmp=0;
        double add=0;
        double sub=0;
        for(int w=0;w<W.cols;w++)
        {
            add=ims_detM.at<double>(h,0)+ims_L.at<double>(h,w);
            sub=ims_detM.at<double>(h,0)-ims_L.at<double>(h,w);

            double sub_tmp=sub*sub/(0.05*0.05);
            double add_tmp=add*add/(0.05*0.05);
            tmp=qLn(qExp(-sub_tmp)+qExp(-add_tmp));
            ims_U.at<double>(h,w)=tmp;
        }
    }


    QList<double> MeanList;
    for(int w=0;w<ims_U.cols-1;w++)
    {
        double tmp=0;
        for(int h=0;h<ims_U.rows;h++)
        {
            tmp+=ims_U.at<double>(h,w);
        }
        tmp=tmp/ims_U.rows;
        MeanList.append(tmp);
    }

    int idx=0;
    double tmp=MeanList[0];
    for(int i=0;i<MeanList.count();i++)
    {
        if(MeanList[i]>tmp)
        {
            tmp=MeanList[i];
            idx=i;
        }
    }
    //qDebug()<<(QString::number(idx)+":")<<(QString::number(tmp));





    cv::Mat gray=cv::Mat::eye(IMAGE.rows,IMAGE.cols,CV_8UC1);
    for(int h=0;h<IMAGE.rows;h++)
    {
        //double tmp=0;
        for(int w=0;w<IMAGE.cols;w++)
        {
            double r=IMAGE.at<cv::Vec3b>(h,w)[0]*W.at<double>(2,idx);
            double g=IMAGE.at<cv::Vec3b>(h,w)[1]*W.at<double>(1,idx);
            double b=IMAGE.at<cv::Vec3b>(h,w)[2]*W.at<double>(0,idx);
            gray.at<int8_t>(h,w)=(int)(r+g+b);
        }
    }
    return gray;

}