OpenCV 的 Contrast Preserving Decolorization 源码解析

运行效果为：

 出乎我意料的是，不仅仅保留了对比度，居然还增强了图像的对比度（去雾，不过只适用于比较均匀的雾），不过运行的速度堪忧，500*500的图像都需要 1s 多！

经过 OpenMP 优化，执行时间减少了一半左右

该代码是源于 香港中文大学 计算机科学与工程系 的一篇论文 Contrast Preserving Decolorization

其代码已被收录到 OpenCV 的源码中，位于（注意，3.* 以上才有）

以下是在通读了论文《IEEE International Conference on Computational Photography(ICCP), 2012》之后，对原始代码做的详细注释（ OpenCV 原代码中有多处 Bug，我已 fix 掉）：

  

   

    

     

    
    

     

      void deColor(InputArray _src, OutputArray _dst, OutputArray _color_boost)
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     

              Mat I = _src.getMat();
     
    
   

    

     

    
    

     

              _dst.create(I.size(), CV_8UC1);
     
    
   

    

     

    
    

     

              Mat dst = _dst.getMat();
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              _color_boost.create(I.size(), CV_8UC3);
     
    
   

    

     

    
    

     

              Mat color_boost = _color_boost.getMat();
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        CV_Assert(!I.empty() && (I.channels() == 3));
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        // Parameter Setting
     
    
   

    

     

    
    

     
        const int maxIter = 15;
     
    
   

    

     

    
    

     
        const double tol = .0001;
     
    
   

    

     

    
    

     
        int iterCount = 0;
     
    
   

    

     

    
    

     
        double E = 0;
     
    
   

    

     

    
    

     
        double pre_E = std::numeric_limits<double>::infinity(); // 返回编译器允许的 double 型的正无穷大
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              Mat img;
     
    
   

    

     

    
    

     

              I.convertTo(img, CV_32FC3, 1.0 / 255.0); // 8UC3 -> 32FC3，归一化
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        // Initialization
     
    
   

    

     

    
    

     

              Decolor obj;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              vector <double> Cg;
     
    
   

    

     

    
    

     

              vector < vector <double> > polyGrad;
     
    
   

    

     

    
    

     

              vector <Vec3i> comb;
     
    
   

    

     

    
    

     

              vector <double> alf;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              obj.gradSystem(img, polyGrad, Cg, comb);
     
    
   

    

     

    
    

     

              obj.weakOrder(img, alf);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        // Solver
     
    
   

    

     

    
    

     
		
     
    
   

    

     

    
    

     

              Mat Mt = Mat(int(polyGrad.size()), int(polyGrad[0].size()), CV_32FC1);
     
    
   

    

     

    
    

     

              obj.weightUpdateMatrix(polyGrad, Cg, Mt);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              vector <double> wei;
     
    
   

    

     

    
    

     

              obj.weightInit(comb, wei);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

               main loop starting 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              vector <double> G_pos(alf.size());
     
    
   

    

     

    
    

     

              vector <double> G_neg(alf.size());
     
    
   

    

     

    
    

     

              vector <double> EXPsum(G_pos.size());
     
    
   

    

     

    
    

     

              vector <double> EXPterm(G_pos.size());
     
    
   

    

     

    
    

     

              vector <double> temp(polyGrad[0].size());
     
    
   

    

     

    
    

     

              vector <double> temp1(polyGrad[0].size());
     
    
   

    

     

    
    

     

              vector <double> temp2(EXPsum.size());
     
    
   

    

     

    
    

     

              vector <double> wei1(polyGrad.size());
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        while (sqrt(pow(E - pre_E, 2)) > tol)
     
    
   

    

     

    
    

     

              {
     
    
   

    

     

    
    

     

                  iterCount += 1;
     
    
   

    

     

    
    

     

                  pre_E = E;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
			// 梯度图大小
     
    
   

    

     

    
    

     
            for (size_t i = 0; i < polyGrad[0].size(); i++)
     
    
   

    

     

    
    

     

                  {
     
    
   

    

     

    
    

     
                double val = 0.0;
     
    
   

    

     

    
    

     
				// 基的个数，公式（10）
     
    
   

    

     

    
    

     
                for (size_t j = 0; j < polyGrad.size(); j++)
     
    
   

    

     

    
    

     

                          val = val + (polyGrad[j][i] * wei[j]);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
				// 公式（4） 的绝对值内部
     
    
   

    

     

    
    

     

                      temp[i] = val - Cg[i];
     
    
   

    

     

    
    

     

                      temp1[i] = val + Cg[i];
     
    
   

    

     

    
    

     

                  }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
			// 近似公式（4）
     
    
   

    

     

    
    

     
            for (size_t i = 0; i < alf.size(); i++)
     
    
   

    

     

    
    

     

                  {
     
    
   

    

     

    
    

     
                const double sqSigma = obj.sigma * obj.sigma;
     
    
   

    

     

    
    

     
                const double pos = ((1 + alf[i]) / 2) * exp(-1.0 * 0.5 * (temp[i] * temp[i]) / sqSigma);
     
    
   

    

     

    
    

     
                const double neg = ((1 - alf[i]) / 2) * exp(-1.0 * 0.5 * (temp1[i] * temp1[i]) / sqSigma);
     
    
   

    

     

    
    

     

                      G_pos[i] = pos;
     
    
   

    

     

    
    

     

                      G_neg[i] = neg;
     
    
   

    

     

    
    

     

                  }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
			// 近似公式（12）（14）
     
    
   

    

     

    
    

     
            for (size_t i = 0; i < G_pos.size(); i++)
     
    
   

    

     

    
    

     

                      EXPsum[i] = G_pos[i] + G_neg[i];
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
            for (size_t i = 0; i < EXPsum.size(); i++)
     
    
   

    

     

    
    

     

                      temp2[i] = (EXPsum[i] == 0) ? 1.0 : 0.0;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
            for (size_t i = 0; i < G_pos.size(); i++)
     
    
   

    

     

    
    

     

                      EXPterm[i] = (G_pos[i] - G_neg[i]) / (EXPsum[i] + temp2[i]);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
			// 更新权重
     
    
   

    

     

    
    

     
            for (int i = 0; i < int(polyGrad.size()); i++)
     
    
   

    

     

    
    

     

                  {
     
    
   

    

     

    
    

     
                double val1 = 0.0;
     
    
   

    

     

    
    

     
                for (int j = 0; j < int(polyGrad[0].size()); j++)
     
    
   

    

     

    
    

     

                      {
     
    
   

    

     

    
    

     

                          val1 = val1 + (Mt.at<float>(i, j) * EXPterm[j]);
     
    
   

    

     

    
    

     

                      }
     
    
   

    

     

    
    

     

                      wei1[i] = val1;
     
    
   

    

     

    
    

     

                  }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
			// 替换权重
     
    
   

    

     

    
    

     
            for (size_t i = 0; i < wei.size(); i++)
     
    
   

    

     

    
    

     

                      wei[i] = wei1[i];
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
			// 公式（11）
     
    
   

    

     

    
    

     

                  E = obj.energyCalcu(Cg, polyGrad, wei);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
            if (iterCount > maxIter)
     
    
   

    

     

    
    

     
                break;
     
    
   

    

     

    
    

     

              }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              Mat Gray = Mat::zeros(img.size(), CV_32FC1);
     
    
   

    

     

    
    

     

              obj.grayImContruct(wei, img, Gray);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              Gray.convertTo(dst, CV_8UC1, 255);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        /// Contrast Boosting /
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              Mat lab;
     
    
   

    

     

    
    

     
        cvtColor(I, lab, COLOR_BGR2Lab);
     
    
   

    

     

    
    

     

              vector <Mat> lab_channel;
     
    
   

    

     

    
    

     
        split(lab, lab_channel);
     
    
   

    

     

    
    

     

              dst.copyTo(lab_channel[0]); // 仅仅替换 L 通道
     
    
   

    

     

    
    

     
        merge(lab_channel, lab);
     
    
   

    

     

    
    

     
        cvtColor(lab, color_boost, COLOR_Lab2BGR);
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
    double Decolor::energyCalcu(const vector <double> &Cg, const vector < vector <double> > &polyGrad, const vector <double> &wei) const
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     
        const size_t size = polyGrad[0].size();
     
    
   

    

     

    
    

     

              vector <double> energy(size);
     
    
   

    

     

    
    

     

              vector <double> temp(size);
     
    
   

    

     

    
    

     

              vector <double> temp1(size);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
		// 公式（11）两个 exp {} 中的部分
     
    
   

    

     

    
    

     
		// 公式中的 l(x,y)
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < polyGrad[0].size(); i++)
     
    
   

    

     

    
    

     

              {
     
    
   

    

     

    
    

     
            double val = 0.0;
     
    
   

    

     

    
    

     
			// 公式中的 i
     
    
   

    

     

    
    

     
            for (size_t j = 0; j < polyGrad.size(); j++)
     
    
   

    

     

    
    

     

                      val = val + (polyGrad[j][i] * wei[j]);
     
    
   

    

     

    
    

     

                  temp[i] = val - Cg[i];
     
    
   

    

     

    
    

     

                  temp1[i] = val + Cg[i];
     
    
   

    

     

    
    

     

              }
     
    
   

    

     

    
    

     
		// 注意这里和公式（11）不同，左右两边的比重都为 1
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < polyGrad[0].size(); i++)
     
    
   

    

     

    
    

     

                  energy[i] = -1.0 * log(exp(-1.0 * pow(temp[i], 2) / sigma) + exp(-1.0 * pow(temp1[i], 2) / sigma));
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        double sum = 0.0;
     
    
   

    

     

    
    

     
		// 把整幅图的能量（代价）都相加起来
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < polyGrad[0].size(); i++)
     
    
   

    

     

    
    

     

                  sum += energy[i];
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
		// 平均
     
    
   

    

     

    
    

     
        return (sum / polyGrad[0].size());
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

          Decolor::Decolor()
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     

              kernelx = Mat(1, 2, CV_32FC1);
     
    
   

    

     

    
    

     

              kernely = Mat(2, 1, CV_32FC1);
     
    
   

    

     

    
    

     

              kernelx.at<float>(0, 0) = 1.0; // 1.， -1.
     
    
   

    

     

    
    

     

              kernelx.at<float>(0, 1) = -1.0;
     
    
   

    

     

    
    

     

              kernely.at<float>(0, 0) = 1.0;  // 1.; -1.
     
    
   

    

     

    
    

     

              kernely.at<float>(1, 0) = -1.0;
     
    
   

    

     

    
    

     

              order = 2; // degree
     
    
   

    

     

    
    

     

              sigma = 0.02f;
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
    vector<double> Decolor::product(const vector <Vec3i> &comb, const double initRGB[3])
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     

              vector <double> res(comb.size());
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
		// 从 comb 容器中，逐个取出 vec3 和 rgb 进行点乘
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < comb.size(); i++)
     
    
   

    

     

    
    

     

              {
     
    
   

    

     

    
    

     
            double dp = 0.0;
     
    
   

    

     

    
    

     
            for (int j = 0; j < 3; j++)
     
    
   

    

     

    
    

     

                      dp = dp + (comb[i][j] * initRGB[j]);
     
    
   

    

     

    
    

     

                  res[i] = dp;
     
    
   

    

     

    
    

     

              }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
		// 返回每个 vec3 点乘的结果
     
    
   

    

     

    
    

     
        return res;
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
    // 计算横向的单通道梯度图
     
    
   

    

     

    
    

     
    void Decolor::singleChannelGradx(const Mat &img, Mat &dest) const
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     
        const int w = img.size().width;
     
    
   

    

     

    
    

     
        // kernels.cols/2-1, kernelx.rows/2-1
     
    
   

    

     

    
    

     
        // 调整卷积图像的锚点，默认是 (-1,-1)
     
    
   

    

     

    
    

     
        // Σ kernel(x',y')*src(x+x'-anchor.x, y+y'-anchor.y)
     
    
   

    

     

    
    

     
        const Point anchor(kernelx.cols - kernelx.cols / 2 - 1, kernelx.rows - kernelx.rows / 2 - 1);
     
    
   

    

     

    
    

     
        filter2D(img, dest, -1, kernelx, anchor, 0.0, BORDER_CONSTANT); // 超出边界的地方用常数填充
     
    
   

    

     

    
    

     

              dest.col(w - 1) = 0.0;// 最右侧设置为 0
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
    // 计算纵向的单通道梯度图
     
    
   

    

     

    
    

     
    void Decolor::singleChannelGrady(const Mat &img, Mat &dest) const
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     
        const int h = img.size().height;
     
    
   

    

     

    
    

     
        const Point anchor(kernely.cols - kernely.cols / 2 - 1, kernely.rows - kernely.rows / 2 - 1);
     
    
   

    

     

    
    

     
        filter2D(img, dest, -1, kernely, anchor, 0.0, BORDER_CONSTANT);
     
    
   

    

     

    
    

     

              dest.row(h - 1) = 0.0;
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
	// 计算横向和纵向的梯度（转置）并合并到一个容器中
     
    
   

    

     

    
    

     
    void Decolor::gradvector(const Mat &img, vector <double> &grad) const
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     

              Mat dest;
     
    
   

    

     

    
    

     

              Mat dest1;
     
    
   

    

     

    
    

     
        singleChannelGradx(img, dest);
     
    
   

    

     

    
    

     
        singleChannelGrady(img, dest1);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        // 得到转置的单通道梯度图
     
    
   

    

     

    
    

     

              Mat d_trans = dest.t();
     
    
   

    

     

    
    

     

              Mat d1_trans = dest1.t();
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        const int height = d_trans.size().height;
     
    
   

    

     

    
    

     
        const int width = d_trans.size().width;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        // 把两张梯度图合并到一个容器当中（前：横向，后：纵向）
     
    
   

    

     

    
    

     
        // OpenCV 源码中此处有越界问题（难以置信）
     
    
   

    

     

    
    

     

              grad.resize(width * height * 2);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        for (int i = 0; i < height; i++)
     
    
   

    

     

    
    

     
            for (int j = 0; j < width; j++)
     
    
   

    

     

    
    

     

                      grad[i * width + j] = d_trans.at<float>(i, j);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        const int offset = width * height;
     
    
   

    

     

    
    

     
        for (int i = 0; i < height; i++)
     
    
   

    

     

    
    

     
            for (int j = 0; j < width; j++)
     
    
   

    

     

    
    

     

                      grad[offset + i * width + j] = d1_trans.at<float>(i, j);
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
	// 计算 CIELab 空间的颜色对比度
     
    
   

    

     

    
    

     
    void Decolor::colorGrad(const Mat &img, vector <double> &Cg) const
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     

              Mat lab;
     
    
   

    

     

    
    

     
        // 转换到 CIELab 颜色空间
     
    
   

    

     

    
    

     
        cvtColor(img, lab, COLOR_BGR2Lab);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              vector <Mat> lab_channel;
     
    
   

    

     

    
    

     
        split(lab, lab_channel);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              vector <double> ImL;
     
    
   

    

     

    
    

     

              vector <double> Ima;
     
    
   

    

     

    
    

     

              vector <double> Imb;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        // 分别计算 Lab 三个通道的梯度图
     
    
   

    

     

    
    

     
        gradvector(lab_channel[0], ImL);
     
    
   

    

     

    
    

     
        gradvector(lab_channel[1], Ima);
     
    
   

    

     

    
    

     
        gradvector(lab_channel[2], Imb);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              Cg.resize(ImL.size());
     
    
   

    

     

    
    

     
		// 计算颜色对比度（Color Contrast）
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < ImL.size(); i++)
     
    
   

    

     

    
    

     

              {
     
    
   

    

     

    
    

     
            const double res = sqrt(pow(ImL[i], 2) + pow(Ima[i], 2) + pow(Imb[i], 2)) / 100;
     
    
   

    

     

    
    

     

                  Cg[i] = res;
     
    
   

    

     

    
    

     

              }
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
    // 构造一个 vec3 张量，并添加到 comb 容器中
     
    
   

    

     

    
    

     
    void Decolor::addVector(vector <Vec3i> &comb, int &idx, int r, int g, int b)
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     

              comb.push_back(Vec3i(r, g, b));
     
    
   

    

     

    
    

     

              idx++;
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
	// 保存梯度容器到一个容器的容器
     
    
   

    

     

    
    

     
    void Decolor::addToVectorPoly(vector < vector <double> > &polyGrad, const vector <double> &curGrad, int &idx1)
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     

              polyGrad.push_back(curGrad);
     
    
   

    

     

    
    

     

              idx1++;
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
	// 和论文上公式不同
     
    
   

    

     

    
    

     
    void Decolor::weakOrder(const Mat &src, vector <double> &alf) const
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     
        const int h = src.size().height;
     
    
   

    

     

    
    

     
        const int w = src.size().width;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              cv::Mat img;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        if ((h + w) > 800)
     
    
   

    

     

    
    

     

              {
     
    
   

    

     

    
    

     
            const double sizefactor = double(800) / (h + w);
     
    
   

    

     

    
    

     
            resize(src, img, Size(cvRound(w * sizefactor), cvRound(h * sizefactor)));
     
    
   

    

     

    
    

     

              }
     
    
   

    

     

    
    

     
        else
     
    
   

    

     

    
    

     

                  img = src;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              Mat curIm = Mat(img.size(), CV_32FC1);
     
    
   

    

     

    
    

     

              vector <Mat> rgb_channel;
     
    
   

    

     

    
    

     
        split(img, rgb_channel);// 又对图像进行了分裂
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              vector <double> Rg, Gg, Bg;
     
    
   

    

     

    
    

     
        gradvector(rgb_channel[2], Rg); // 计算 RGB 三个通道的梯度图
     
    
   

    

     

    
    

     
        gradvector(rgb_channel[1], Gg);
     
    
   

    

     

    
    

     
        gradvector(rgb_channel[0], Bg);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              vector <double> t1(Rg.size()), t2(Rg.size()), t3(Rg.size());
     
    
   

    

     

    
    

     

              vector <double> tmp1(Rg.size()), tmp2(Rg.size()), tmp3(Rg.size());
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        const double level = .05;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
		// 和 level 、-level 进行比较
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < Rg.size(); i++)
     
    
   

    

     

    
    

     

              {
     
    
   

    

     

    
    

     

                  t1[i] = (Rg[i] > level) ? 1.0 : 0.0;
     
    
   

    

     

    
    

     

                  t2[i] = (Gg[i] > level) ? 1.0 : 0.0;
     
    
   

    

     

    
    

     

                  t3[i] = (Bg[i] > level) ? 1.0 : 0.0;
     
    
   

    

     

    
    

     

                  tmp1[i] = (Rg[i] < -1.0 * level) ? 1.0 : 0.0;
     
    
   

    

     

    
    

     

                  tmp2[i] = (Gg[i] < -1.0 * level) ? 1.0 : 0.0;
     
    
   

    

     

    
    

     

                  tmp3[i] = (Bg[i] < -1.0 * level) ? 1.0 : 0.0;
     
    
   

    

     

    
    

     

              }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              alf.resize(Rg.size());
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < Rg.size(); i++)
     
    
   

    

     

    
    

     

                  alf[i] = (t1[i] * t2[i] * t3[i]);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < Rg.size(); i++)
     
    
   

    

     

    
    

     

                  alf[i] -= tmp1[i] * tmp2[i] * tmp3[i];
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
	// 构造 polynomial space 每个基的梯度图，还有得到 CIELab 颜色空间的对比度，以及 polynomial space 的基
     
    
   

    

     

    
    

     
    void Decolor::gradSystem(const Mat &src, vector < vector < double > > &polyGrad,
     
    
   

    

     

    
    

     

       vector < double > &Cg, vector <Vec3i> &comb) const
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     
        int h = src.size().height;
     
    
   

    

     

    
    

     
        int w = src.size().width;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              cv::Mat img;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        // 如果宽高和超过一定大小则进行缩放
     
    
   

    

     

    
    

     
        if ((h + w) > 800)
     
    
   

    

     

    
    

     

              {
     
    
   

    

     

    
    

     
            const double sizefactor = double(800) / (h + w);
     
    
   

    

     

    
    

     
            resize(src, img, Size(cvRound(w * sizefactor), cvRound(h * sizefactor)));
     
    
   

    

     

    
    

     

              }
     
    
   

    

     

    
    

     
        else
     
    
   

    

     

    
    

     

                  img = src;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              h = img.size().height;
     
    
   

    

     

    
    

     

              w = img.size().width;
     
    
   

    

     

    
    

     
        colorGrad(img, Cg);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
		// 将一副图像映射到 polynomial space
     
    
   

    

     

    
    

     

              Mat curIm = Mat(img.size(), CV_32FC1);
     
    
   

    

     

    
    

     

              vector <Mat> rgb_channel;
     
    
   

    

     

    
    

     
        split(img, rgb_channel); // 得到 BGR 三个通道
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        int idx = 0, idx1 = 0;
     
    
   

    

     

    
    

     
        for (int r = 0; r <= order; r++)
     
    
   

    

     

    
    

     
            for (int g = 0; g <= order; g++)
     
    
   

    

     

    
    

     
                for (int b = 0; b <= order; b++)
     
    
   

    

     

    
    

     

                      {
     
    
   

    

     

    
    

     
                    if ((r + g + b) <= order && (r + g + b) > 0)
     
    
   

    

     

    
    

     

                          {
     
    
   

    

     

    
    

     
						// 保存 polynomial space 的基
     
    
   

    

     

    
    

     
                        addVector(comb, idx, r, g, b);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
						// 每个 polynomial space 的基（r, g, b） 都要对整张图像进行计算（w*h 大小）
     
    
   

    

     

    
    

     
                        for (int i = 0; i < h; i++)
     
    
   

    

     

    
    

     
                            for (int j = 0; j < w; j++)
     
    
   

    

     

    
    

     
								// 映射每一个 rgb 像素
     
    
   

    

     

    
    

     

                                      curIm.at<float>(i, j) =
     
    
   

    

     

    
    

     
                                    pow(rgb_channel[2].at<float>(i, j), r) * 
     
    
   

    

     

    
    

     
									pow(rgb_channel[1].at<float>(i, j), g) *
     
    
   

    

     

    
    

     
                                    pow(rgb_channel[0].at<float>(i, j), b);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

                              vector <double> curGrad;
     
    
   

    

     

    
    

     
                        gradvector(curIm, curGrad);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
						// 保存每个基计算的梯度图
     
    
   

    

     

    
    

     
                        addToVectorPoly(polyGrad, curGrad, idx1);
     
    
   

    

     

    
    

     

                          }
     
    
   

    

     

    
    

     

                      }
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
    void Decolor::weightUpdateMatrix(const vector < vector <double> > &poly, const vector <double> &Cg, Mat &X)
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     
		// 容器转为矩阵
     
    
   

    

     

    
    

     
        const int size = static_cast<int>(poly.size());
     
    
   

    

     

    
    

     
        const int size0 = static_cast<int>(poly[0].size());
     
    
   

    

     

    
    

     

              Mat P = Mat(size, size0, CV_32FC1);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        for (int i = 0; i < size; i++)
     
    
   

    

     

    
    

     
            for (int j = 0; j < size0; j++)
     
    
   

    

     

    
    

     

                      P.at<float>(i, j) = static_cast<float>(poly[i][j]);
     
    
   

    

     

    
    

     
		// 转置
     
    
   

    

     

    
    

     
        const Mat P_trans = P.t();
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              Mat B = Mat(size, size0, CV_32FC1);
     
    
   

    

     

    
    

     
        for (int i = 0; i < size; i++)
     
    
   

    

     

    
    

     

              {
     
    
   

    

     

    
    

     
            for (int j = 0, end = int(Cg.size()); j < end; j++)
     
    
   

    

     

    
    

     

                      B.at<float>(i, j) = static_cast<float>(poly[i][j] * Cg[j]);
     
    
   

    

     

    
    

     

              }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
		// 得到一个方阵，大小为 size*size
     
    
   

    

     

    
    

     

              Mat A = P * P_trans;
     
    
   

    

     

    
    

     
		// 求解线性方程组，这里的 X 应该是论文公式（14）的部分
     
    
   

    

     

    
    

     
		// cv::solve -- https://docs.opencv.org/3.1.0/d2/de8/group__core__array.html#ga12b43690dbd31fed96f213eefead2373
     
    
   

    

     

    
    

     
		// DECOMP_NORMAL -- https://docs.opencv.org/3.1.0/d2/de8/group__core__array.html#gaaf9ea5dcc392d5ae04eacb9920b9674c
     
    
   

    

     

    
    

     
        // 这意味着求解的方程是 src1 T⋅src1⋅dst = src1 T src2，而不是原始方程 src1⋅dst = src2
     
    
   

    

     

    
    

     
        solve(A, B, X, DECOMP_NORMAL);
     
    
   

    

     

    
    

     
        //--------------------- 
     
    
   

    

     

    
    

     
        // DECOMP_LU（LU分解）
     
    
   

    

     

    
    

     
        // http://blog.csdn.net/myhaspl/article/details/49450743
     
    
   

    

     

    
    

     
        // DECOMP_CHOLESKY（CHOLESKY分解）
     
    
   

    

     

    
    

     
        // http://blog.csdn.net/acdreamers/article/details/44656847
     
    
   

    

     

    
    

     
        // DECOMP_EIG（EIG分解）
     
    
   

    

     

    
    

     
        // DECOMP_SVD（SVD分解）
     
    
   

    

     

    
    

     
        // http://www.cnblogs.com/LeftNotEasy/archive/2011/01/19/svd-and-applications.html
     
    
   

    

     

    
    

     
        // DECOMP_QR（QR分解）
     
    
   

    

     

    
    

     
        // http://blog.sina.com.cn/s/blog_3f41287a0101ke2s.html
     
    
   

    

     

    
    

     
        //--------------------- 
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
    void Decolor::weightInit(const vector <Vec3i> &comb, vector <double> &wei)
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     
        double initRGB[3] = { .33, .33, .33 };
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
		// 通过 polynomial space 的基和 rgb 系数进行点乘，对 weight 进行初始化
     
    
   

    

     

    
    

     

              wei = product(comb, initRGB);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              vector <int> sum(comb.size());
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < comb.size(); i++)
     
    
   

    

     

    
    

     

                  sum[i] = (comb[i][0] + comb[i][1] + comb[i][2]);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
		// 除了 r,g,b 这三种基，其他基的权重初始化为 0
     
    
   

    

     

    
    

     
        for (size_t i = 0; i < sum.size(); i++)
     
    
   

    

     

    
    

     

              {
     
    
   

    

     

    
    

     
            if (sum[i] == 1)
     
    
   

    

     

    
    

     

                      wei[i] = wei[i] * double(1);
     
    
   

    

     

    
    

     
            else
     
    
   

    

     

    
    

     

                      wei[i] = wei[i] * double(0);
     
    
   

    

     

    
    

     

              }
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              sum.clear();
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

          }
     
    
   

    

     

    
    

     
    
     
    
   

    

     

    
    

     
    // 将迭代出各个基的权重和各个基相乘，并把结果累加到原灰度上
     
    
   

    

     

    
    

     
    void Decolor::grayImContruct(vector <double> &wei, const Mat &img, Mat &Gray) const
     
    
   

    

     

    
    

     

          {
     
    
   

    

     

    
    

     
        const int h = img.size().height;
     
    
   

    

     

    
    

     
        const int w = img.size().width;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              vector <Mat> rgb_channel;
     
    
   

    

     

    
    

     
        split(img, rgb_channel);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        int kk = 0;
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     
        for (int r = 0; r <= order; r++)
     
    
   

    

     

    
    

     
            for (int g = 0; g <= order; g++)
     
    
   

    

     

    
    

     
                for (int b = 0; b <= order; b++)
     
    
   

    

     

    
    

     
                    if ((r + g + b) <= order && (r + g + b) > 0)
     
    
   

    

     

    
    

     

                          {
     
    
   

    

     

    
    

     
                        for (int i = 0; i < h; i++)
     
    
   

    

     

    
    

     
                            for (int j = 0; j < w; j++)
     
    
   

    

     

    
    

     

                                      Gray.at<float>(i, j) = Gray.at<float>(i, j) +
     
    
   

    

     

    
    

     
                                                       static_cast<float>(wei[kk]) * pow(rgb_channel[2].at<float>(i, j), r) * pow(rgb_channel[1].at<float>(i, j), g) *
     
    
   

    

     

    
    

     
                                                       pow(rgb_channel[0].at<float>(i, j), b);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

                              kk = kk + 1; // 遍历各个基
     
    
   

    

     

    
    

     

                          }
     
    
   

    

     

    
    

     
        // 找出最值，并调整归一化便于显示
     
    
   

    

     

    
    

     
        double minval, maxval;
     
    
   

    

     

    
    

     
        minMaxLoc(Gray, &minval, &maxval);
     
    
   

    

     

    
    

     
 
     
    
   

    

     

    
    

     

              Gray -= minval;
     
    
   

    

     

    
    

     

              Gray /= maxval - minval;
     
    
   

    

     

    
    

     

          }
     
    
  
 

文章来源: panda1234lee.blog.csdn.net，作者：panda1234lee，版权归原作者所有，如需转载，请联系作者。

原文链接：panda1234lee.blog.csdn.net/article/details/87978829