Add mcc infer uint8

modules/cudaimgproc/include/opencv2/cudaimgproc.hpp

@@ -148,6 +148,16 @@ The methods support arbitrary permutations of the original channels, including r
  */
 CV_EXPORTS void swapChannels(InputOutputArray image, const int dstOrder[4], Stream& stream = Stream::Null());
 
+
+/** @brief calibrate input CV_8UC3 image (GpuMat, BGR format) and return new CV_8UC3 image (GpuMat, BGR format)
+
+@param src - input GpuMat with type CV_8UC3, BGR format
+@param ccm - input CCM with type CV_64F and size 3x3
+@param dst - output GpuMat with type CV_8UC3, BGR format
+The methods support arbitrary permutations of the original channels, including replication.
+ */
+CV_EXPORTS_W void calibrateImageF32C3(InputArray src, Mat ccm, OutputArray dst, Stream& stream = Stream::Null());
+
 /** @brief Routines for correcting image color gamma.
 
 @param src Source image (3- or 4-channel 8 bit).

modules/cudaimgproc/src/cuda/hist.cu

@@ -52,6 +52,49 @@ using namespace cv::cuda::device;
 
 namespace hist
 {
+
+__device__ __forceinline__ float fromLFuncEW(const float& x) {
+    const float gamma2 = 2.3999999999999999f;
+    const float alpha = 1.0549999999999999f;
+    const float beta = 0.0030399346397784318f;
+    const float phi = 12.923210180787855f;
+    if (x > beta)
+        return alpha * pow(x, 1 / gamma2) - (alpha - 1);
+    else if (x >= -beta)
+        return x * phi;
+    else
+        return -(alpha * pow(-x, 1 / gamma2) - (alpha - 1));
+}
+
+__global__ void calibrateKernalF32C3(PtrStepSz<float3> input, PtrStepSz<float3> output, PtrStepSz<float3> ccm) {
+    int x = blockDim.x * blockIdx.x + threadIdx.x;
+    int y = blockDim.y * blockIdx.y + threadIdx.y;
+    int cols = input.cols;
+    int rows = input.rows;
+    const float gamma1 = 2.2000000000000002f;
+
+    if (x < cols && y < rows) {
+        const float3 in = input(y, x);
+        const float3 gamma = make_float3(pow(in.x / 255.f, gamma1), pow(in.y / 255.f, gamma1), pow(in.z / 255.f, gamma1));
+        const float3 mult = make_float3(gamma.x*ccm(0, 0).x + gamma.y*ccm(0, 0).y + gamma.z*ccm(0, 0).z,
+                                        gamma.x*ccm(1, 0).x + gamma.y*ccm(1, 0).y + gamma.z*ccm(1, 0).z,
+                                        gamma.x*ccm(2, 0).x + gamma.y*ccm(2, 0).y + gamma.z*ccm(2, 0).z);
+        output(y, x) = make_float3(fromLFuncEW(mult.x)*255.f, fromLFuncEW(mult.y)*255.f, fromLFuncEW(mult.z)*255.f);
+    }
+}
+
+void calibrateImageF32C3(PtrStepSz<float3> src, PtrStepSz<float3> dst, PtrStepSz<float3> ccm, cudaStream_t stream) {
+    const dim3 block(32, 8);
+    //const dim3 grid(divUp(src.rows, block.y));
+    const dim3 grid(divUp(src.cols, block.x), divUp(src.rows, block.y));
+    calibrateKernalF32C3<<<grid, block, 0, stream>>>(src, dst, ccm);
+
+    cudaSafeCall( cudaGetLastError() );
+
+    if (stream == 0)
+        cudaSafeCall( cudaDeviceSynchronize() );
+}
+
     template<bool fourByteAligned>
     __global__ void histogram256Kernel(const uchar* src, int cols, int rows, size_t step, int* hist, const int offsetX = 0)
     {

modules/cudaimgproc/src/histogram.cpp

@@ -70,13 +70,44 @@ namespace hist
 {
     void histogram256(PtrStepSzb src, int* hist, const int offsetX, cudaStream_t stream);
     void histogram256(PtrStepSzb src, PtrStepSzb mask, int* hist, const int offsetX, cudaStream_t stream);
+    void calibrateImageF32C3(PtrStepSz<float3> src, PtrStepSz<float3> dst, PtrStepSz<float3> ccm, cudaStream_t stream);
 }
 
 void cv::cuda::calcHist(InputArray _src, OutputArray _hist, Stream& stream)
 {
     calcHist(_src, cv::cuda::GpuMat(), _hist, stream);
 }
 
+void cv::cuda::calibrateImageF32C3(InputArray _src, Mat ccm, OutputArray _dst, Stream& _stream) {
+    GpuMat src = _src.getGpuMat();
+    CV_Assert(src.type() == CV_8UC3);
+    CV_Assert(ccm.type() == CV_64F && ccm.rows == 3 && ccm.cols == 3);
+    cudaStream_t stream = StreamAccessor::getStream(_stream);
+    NppStreamHandler h(stream);
+
+    Mat ccm_float(3, 3, CV_32F);
+    ccm.convertTo(ccm_float, CV_32FC1);
+    cv::flip(ccm_float, ccm_float, -1);
+    cv::transpose(ccm_float, ccm_float);
+    ccm_float = ccm_float.reshape(3, 3);
+    GpuMat ccm_gpu(3, 1, CV_32FC3);
+    ccm_gpu.upload(ccm_float);
+
+
+    GpuMat float_tmp(src.size(), CV_32FC3);
+    src.convertTo(float_tmp, CV_32FC3, _stream);
+    hist::calibrateImageF32C3(float_tmp, float_tmp, ccm_gpu, stream);
+    cudaSafeCall( cudaDeviceSynchronize() );
+
+    _dst.create(src.size(), src.type());
+    GpuMat dst = _dst.getGpuMat();
+
+
+    float_tmp.convertTo(dst, CV_8UC3, _stream);
+    if (stream == 0)
+        cudaSafeCall( cudaDeviceSynchronize() );
+}
+
 void cv::cuda::calcHist(InputArray _src, InputArray _mask, OutputArray _hist, Stream& stream)
 {
     GpuMat src = _src.getGpuMat();

modules/cudaimgproc/test/test_histogram.cpp

@@ -264,6 +264,37 @@ INSTANTIATE_TEST_CASE_P(CUDA_ImgProc, EqualizeHist, testing::Combine(
     ALL_DEVICES,
     DIFFERENT_SIZES));
 
+TEST(EqualizeHist, infer)
+{
+    string path = cvtest::findDataFile("../cv/mcc/mcc_ccm_test.jpg"); // mcc_ccm_test
+    Mat src = imread(path, IMREAD_COLOR);
+    //Mat src = randomMat(Size(4, 4), CV_8UC3);
+
+
+    Mat ccm = (Mat_<double>(3, 3) <<
+        1.571657000185009, -0.280924394459268, 0.08453922109770311,
+        -0.02194530839037822, 1.377608706481052, -0.1803987912579378,
+        -0.02843995357976759, 0.3033531498800233, 1.471473866992415);
+
+    const double gamma1 = 2.2000000000000002;
+    const double gamma2 = 2.3999999999999999;
+    const double alpha = 1.0549999999999999;
+    const double beta = 0.0030399346397784318;
+    const double phi = 12.923210180787855;
+
+    cv::cuda::GpuMat dst;
+    cv::cuda::calibrateImageF32C3(loadMat(src), ccm, dst);
+
+    Mat tmp;
+    dst.download(tmp);
+    //imshow("src", src);
+    //imshow("dst", tmp);
+    //waitKey(0);
+    path = cvtest::findDataFile("../cv/mcc/mcc_ccm_test_res.png");
+    Mat gold_img = imread(path);
+    EXPECT_MAT_NEAR(gold_img, tmp, 1);
+}
+
 TEST(EqualizeHistIssue, Issue18035)
 {
     std::vector<std::string> imgPaths;

modules/mcc/include/opencv2/mcc/ccm.hpp

@@ -507,6 +507,8 @@ class CV_EXPORTS_W ColorCorrectionModel
     */
     CV_WRAP Mat infer(const Mat& img, bool islinear = false);
 
+    CV_WRAP Mat infer_uint8(Mat& img, bool islinear = false);
+
     class Impl;
 private:
     std::shared_ptr<Impl> p;

modules/mcc/src/ccm.cpp

@@ -299,6 +299,38 @@ Mat ColorCorrectionModel::infer(const Mat& img, bool islinear)
     return p->cs.fromL(img_ccm);
 }
 
+Mat ColorCorrectionModel::infer_uint8(Mat& img, bool islinear)
+{
+    CV_Assert(img.type() == CV_8UC3 || img.type() == CV_8UC1 && !img.empty());
+    if (!p->ccm.data)
+        CV_Error(Error::StsBadArg, "No CCM values!" );
+    CV_Assert(p->ccm_type == CCM_TYPE::CCM_3x3);
+    CV_Assert(typeid(*p->linear).name() == typeid(LinearGamma).name());
+    CV_Assert(typeid(p->cs).name() == typeid(sRGB_).name());
+    double gamma = dynamic_cast<LinearGamma*>(p->linear.get())->gamma;
+
+    Mat res_img = LUTGammaCorrection(img, gamma, img);
+
+    Mat ccm_bgr;
+    p->ccm.convertTo(ccm_bgr, CV_32F);
+    flip(ccm_bgr, ccm_bgr, -1);
+
+    Mat float_img;
+    res_img.reshape(1, res_img.rows * res_img.cols).convertTo(float_img, CV_32F);
+
+    Mat mult_res = float_img * ccm_bgr;
+    mult_res = mult_res.reshape(img.channels(), img.rows);
+    mult_res.convertTo(res_img, CV_8UC3);
+    if (islinear)
+        return res_img;
+    sRGB_& color_space = dynamic_cast<sRGB_&>(p->cs);
+    res_img = LUT_EW(res_img, color_space.gamma, color_space.alpha, res_img);
+
+    //mult_res = fromLFuncEW(mult_res, color_space.gamma, color_space.alpha, color_space.beta, color_space.phi);
+    //mult_res.convertTo(res_img, CV_8UC3);
+    return res_img;
+}
+
 void ColorCorrectionModel::Impl::getColor(CONST_COLOR constcolor)
 {
     dst = (GetColor::getColor(constcolor));

modules/mcc/src/utils.cpp

@@ -40,6 +40,69 @@ Mat gammaCorrection(const Mat& src, const double& gamma)
     return elementWise(src, [gamma](double element) -> double { return gammaCorrection_(element, gamma); });
 }
 
+Mat LUTGammaCorrection(const Mat& src, const double& gamma, Mat dst) {
+    Mat lookUpTable(1, 256, CV_8U);
+    uchar* p = lookUpTable.ptr();
+    for(int i = 0; i < 256; i++)
+        p[i] = saturate_cast<uchar>(pow(i / 255.0, gamma) * 255.0);
+    if (dst.empty() || !dst.isContinuous() || dst.total() != src.total() || dst.type() != src.type())
+        dst = Mat(src.rows, src.cols, src.type());
+    LUT(src, lookUpTable, dst);
+    return dst;
+}
+
+// used fromLFuncEW from sRGB_
+Mat LUT_EW(const Mat& src, const double gamma, const double alpha, Mat dst) {
+    Mat lookUpTable(1, 256, CV_8U);
+    uchar* p = lookUpTable.ptr();
+    for(int i = 0; i < 256; i++) {
+        double x = i / 255.;
+        double res = 0.;
+        if (x > 0.)
+            res = alpha * pow(x, 1 / gamma) - (alpha - 1);
+        // other conditions from original fromLFuncEW cannot be executed
+        p[i] = saturate_cast<uchar>(res * 255.0);
+    }
+    if (dst.empty() || !dst.isContinuous() || dst.total() != src.total() || dst.type() != src.type())
+        dst = Mat(src.rows, src.cols, src.type());
+    LUT(src, lookUpTable, dst);
+    return dst;
+}
+
+static inline float fromLFuncEW_core(const float x, const float gamma, const float alpha, const float beta, const float phi)
+{
+    if (x > beta)
+        return alpha * pow(x, 1 / gamma) - (alpha - 1);
+    else if (x >= -beta)
+        return x * phi;
+    else
+        return -(alpha * pow(-x, 1 / gamma) - (alpha - 1));
+}
+
+Mat fromLFuncEW(const Mat& src, const float gamma, const float alpha, const float beta, const float phi, Mat dst)
+{
+    if (dst.empty() || !dst.isContinuous() || dst.total() != src.total() || dst.type() != src.type())
+        dst = src.clone();
+    const int channel = src.channels();
+    if (src.isContinuous()) {
+        const int num_elements = (int)src.total()*channel;
+        const float *psrc = (float*)src.data;
+        float *pdst = (float*)dst.data;
+
+        const int batch = 128;
+        const int N = (num_elements / batch) + ((num_elements % batch) > 0);
+        parallel_for_(Range(0, N),[&](const Range& range)
+        {
+            const int start = range.start * batch;
+            const int end = std::min(range.end*batch, num_elements);
+            for (int i = start; i < end; i++)
+                pdst[i] = fromLFuncEW_core(psrc[i]/255.f, gamma, alpha, beta, phi)*255.f;
+        });
+        return dst;
+    }
+    return dst;
+}
+
 Mat maskCopyTo(const Mat& src, const Mat& mask)
 {
     Mat dst(countNonZero(mask), 1, src.type());

modules/mcc/src/utils.hpp

@@ -45,6 +45,13 @@ double gammaCorrection_(const double& element, const double& gamma);
  */
 Mat gammaCorrection(const Mat& src, const double& gamma);
 
+
+Mat LUTGammaCorrection(const Mat& src, const double& gamma, Mat dst=Mat());
+
+Mat LUT_EW(const Mat& src, double gamma, double alpha, Mat dst=Mat());
+
+Mat fromLFuncEW(const Mat& src, const float gamma, const float alpha, const float beta, const float phi, Mat dst=Mat());
+
 /** @brief maskCopyTo a function to delete unsatisfied elementwise.
     @param src the input array, type of Mat.
     @param mask operation mask that used to choose satisfided elementwise.

modules/mcc/test/test_mcc.cpp

@@ -177,5 +177,41 @@ TEST(CV_mcc_ccm_test, infer)
 }
 
 
+TEST(CV_mcc_ccm_test, infer_uint8)
+{
+    string path = cvtest::findDataFile("mcc/mcc_ccm_test.jpg"); // mcc_ccm_test
+    Mat img = imread(path, IMREAD_COLOR);
+    // read gold calibrate img
+    path = cvtest::findDataFile("mcc/mcc_ccm_test_res.png");
+    Mat gold_img = imread(path);
+
+    // read gold chartsRGB
+    path = cvtest::findDataFile("mcc/mcc_ccm_test.yml");
+    FileStorage fs(path, FileStorage::READ);
+    Mat chartsRGB;
+    FileNode node = fs["chartsRGB"];
+    node >> chartsRGB;
+    fs.release();
+
+    // compute CCM
+    ColorCorrectionModel model(chartsRGB.col(1).clone().reshape(3, chartsRGB.rows/3) / 255., COLORCHECKER_Macbeth);
+    model.run();
+    std::cout << model.getCCM() << endl;
+
+    // compute calibrate image
+    Mat calibratedImage = model.infer_uint8(img);
+    imshow("calibratedImage", calibratedImage);
+    waitKey(0);
+    Mat diff;
+    absdiff(gold_img, calibratedImage, diff);
+    Scalar m = mean(diff);
+    cout << m << endl;
+    imshow("diff", diff);
+    waitKey(0);
+    // check calibrated image
+    EXPECT_MAT_NEAR(gold_img, calibratedImage, 11);
+}
+
+
 } // namespace
 } // namespace opencv_test