Possibly wrong formula implementation in BRISQUE's MSCN algorithm

[ndujar]

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

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

[alalek]

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

[krshrimali]

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]

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]

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]

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.

[clunietp]

Hi @ndujar I finally had some time to take a look at this. I implemented your proposed solution, retrained the model, and obtained the following results on TID2008, where SROCC is the Spearman rank order correlation coefficient:

Original MSCN SROCC: -0.836977
Proposed MSCN SROCC: -0.565237

The MOS scores in TID2008 have a linear association with rated image quality, while there is an inverse association in the LIVE DB R2 dataset, explaining why the Spearman score is negative. With the proposed MSCN implementation, the performance of BRISQUE is notably worse according to the above results.

This makes me wonder if perhaps the paper is incorrect, and the MATLAB/C++ code was correct all along. I'm happy to post source if you'd like to review the implementation of the proposed algorithm.