Possibly wrong formula implementation in BRISQUE's MSCN algorithm

ndujar commented 4 years ago

System information (version)

OpenCV => 4.2
Operating System / Platform => all
Compiler => N/A

Detailed description

I have carefully gone through the transformations implemented in the qualitybrisque.cpp file (lines 147 to 165) to produce the MSCN image (structdis) and I believe there might be a mistake in the way it is conceived, particularly in the obtainment of the denominator.

Steps to reproduce

qualitybrisque.cpp: lines 145-165:

        // calculating MSCN coefficients
        // compute mu (local mean)
        brisque_mat_type mu;
        cv::GaussianBlur(imdist_scaled, mu, cv::Size(7, 7), 7. / 6., 0., cv::BORDER_REPLICATE );

        brisque_mat_type mu_sq;
        cv::pow(mu, double(2.0), mu_sq);

Here the computation of mu_sq leads to the wrong result...

        //compute sigma (local sigma)
        brisque_mat_type sigma;// (imdist_scaled.size(), CV_64FC1, 1);
        cv::multiply(imdist_scaled, imdist_scaled, sigma);
        cv::GaussianBlur(sigma, sigma, cv::Size(7, 7), 7./6., 0., cv::BORDER_REPLICATE );

..because it is not the same the square of the difference as the difference of the square

        cv::subtract(sigma, mu_sq, sigma);
        cv::pow(sigma, double(0.5), sigma);
        cv::add(sigma, Scalar(1.0 / 255), sigma); // to avoid DivideByZero Error

        brisque_mat_type structdis;// (imdist_scaled.size(), CV_64FC1, 1);
        cv::subtract(imdist_scaled, mu, structdis);
        cv::divide(structdis, sigma, structdis);  // structdis is MSCN image

Proposed solution

brisque_mat_type mu;//  (imdist_scaled.size(), CV_64FC1, 1);
cv::GaussianBlur(imdist_scaled, mu, cv::Size(7, 7), 7. / 6., 0., cv::BORDER_REPLICATE );

brisque_mat_type gamma;
cv::subtract(imdist_scaled, mu, gamma);

//compute sigma (local sigma)
brisque_mat_type sigma;// (imdist_scaled.size(), CV_64FC1, 1);
cv::GaussianBlur(gamma, sigma, cv::Size(7, 7), 7./6., 0., cv::BORDER_REPLICATE );
cv::pow(sigma, double(0.5), sigma);
cv::add(sigma, cv::Scalar(1.0 / 255), sigma); // to avoid DivideByZero Error

brisque_mat_type structdis;// (imdist_scaled.size(), CV_64FC1, 1);
cv::divide(gamma, sigma, structdis);  // structdis is MSCN image

Issue submission checklist

[x] I report the issue, it's not a question
[x] I checked the problem with documentation, FAQ, open issues, answers.opencv.org, Stack Overflow, etc and have not found solution
[x] I updated to latest OpenCV version and the issue is still there
[x] There is reproducer code and related data files: videos, images, onnx, etc

alalek commented 4 years ago

/cc @krshrimali @clunietp Could you please take a look on this?

krshrimali commented 4 years ago

Thanks @ndujar for reporting this here. I'll be able to look into this, this weekend. Will update back with my views. Hope it sounds good, @alalek @ndujar. Thanks!

clunietp commented 4 years ago

Hi @ndujar

Here is the original author's code, available from https://live.ece.utexas.edu/research/Quality/index_algorithms.htm

          //compute mu and mu squared
    IplImage* mu = cvCreateImage(cvGetSize(imdist_scaled), IPL_DEPTH_64F, 1);
    cvSmooth( imdist_scaled, mu, CV_GAUSSIAN, 7, 7, 1.16666 );
    IplImage* mu_sq = cvCreateImage(cvGetSize(imdist_scaled), IPL_DEPTH_64F, 1);
    cvMul(mu, mu, mu_sq);

    //compute sigma
    IplImage* sigma = cvCreateImage(cvGetSize(imdist_scaled), IPL_DEPTH_64F, 1);
    cvMul(imdist_scaled, imdist_scaled, sigma);
    cvSmooth(sigma, sigma, CV_GAUSSIAN, 7, 7, 1.16666 );
    cvSub(sigma, mu_sq, sigma);
    cvPow(sigma, sigma, 0.5);

    //compute structdis = (x-mu)/sigma
    cvAddS(sigma, cvScalar(1.0/255), sigma);
    IplImage* structdis = cvCreateImage(cvGetSize(imdist_scaled), IPL_DEPTH_64F, 1);
    cvSub(imdist_scaled, mu, structdis);
    cvDiv(structdis, sigma, structdis);

Matlab version:

window = fspecial('gaussian',7,7/6);
...  
mu            = filter2(window, imdist, 'same');
mu_sq         = mu.*mu;
sigma         = sqrt(abs(filter2(window, imdist.*imdist, 'same') - mu_sq));
structdis     = (imdist-mu)./(sigma+1);

Are you suggesting the current OpenCV implementation is wrong, or the original author's implementation is wrong? It is possible we missed something when porting to OpenCV.

Thanks

ndujar commented 4 years ago

Hi @clunietp

Indeed, the original implementation seems to be mistaken :/

However, the code says otherwise, as they seem to be applying the Gaussian filter to both the squares of the I and mu separately, then computing the pixelwise squared root. Actually, that explains why they need to use the abs() function, when in reality it shouldn't be needed because quadratic products are always positive.

clunietp commented 4 years ago

Interesting. I wonder how your proposed method compares to the original method. From what I remember, the current implementation did not score nearly as well on TID2008 as is documented in the paper, and I wonder if this may be the cause. Ref: Table VIII in https://www.live.ece.utexas.edu/publications/2012/TIP%20BRISQUE.pdf

Are you (or anyone) able to test the current impl vs the proposed changes on TID2008? I'd be interested in the result. Unfortunately I don't really have a lot of time right now. Eval code is here: https://github.com/opencv/opencv_contrib/blob/master/modules/quality/samples/brisque_eval_tid2008.cpp

If the revised method does improve performance on TID2008 and brings it closer to Mittal's original work, then I would agree with changing the current implementation.

opencv / opencv_contrib