brisque

gistfile1.txt

#include "brisque.h"


//function definitions
static void AGGDfit(const cv::Mat & structdis, double& lsigma_best, double& rsigma_best, double& gamma_best)
{
    cv::Mat_<double> Imarr(structdis);
    long int poscount=0, negcount=0;
    double possqsum=0, negsqsum=0, abssum=0;
    for (int i=0;i<structdis.rows;i++)
    {
        for (int j =0; j<structdis.cols; j++)
        {
            double pt = Imarr(i, j);
            if (pt > 0)
            {
                poscount++;
                possqsum += pt*pt;
                abssum += pt;
            }
            else if (pt < 0)
            {
                negcount++;
                negsqsum += pt*pt;
                abssum -= pt;
            }
        }
    }
    lsigma_best = pow(negsqsum/negcount, 0.5); //1st output parameter set
    rsigma_best = pow(possqsum/poscount, 0.5); //2nd output parameter set

    double gammahat = lsigma_best/rsigma_best;
    long int totalcount = structdis.total();
    double rhat = pow(abssum/totalcount, static_cast<double>(2))/((negsqsum + possqsum)/totalcount);
    double rhatnorm = rhat*(pow(gammahat,3) +1)*(gammahat+1)/pow(pow(gammahat,2)+1,2);

    double prevgamma = 0;
    double prevdiff = 1e10;
    float sampling = 0.001;
    for (float gam=0.2; gam<10; gam+=sampling) //possible to coarsen sampling to quicken the code, with some loss of accuracy
    {
        double r_gam = tgamma(2/gam)*tgamma(2/gam)/(tgamma(1/gam)*tgamma(3/gam));
        double diff = abs(r_gam-rhatnorm);
        if(diff> prevdiff) break;
        prevdiff = diff;
        prevgamma = gam;
    }
    gamma_best = prevgamma;
}


//function definitions
void ComputeBrisqueFeature(const cv::Mat &orig, vector<double>& featurevector)
{
    cv::Mat orig_bw, orig_bw_int;
    cv::cvtColor(orig, orig_bw_int, CV_RGB2GRAY);
    orig_bw_int.convertTo(orig_bw, CV_64F, 1.0/255);
    //orig_bw now contains the grayscale image normalized to the range 0,1

    int scalenum = 2;
    for (int itr_scale = 1; itr_scale<=scalenum; itr_scale++)
    {
        cv::Size siz(orig_bw.cols/pow((double)2,itr_scale-1), orig_bw.rows/pow((double)2,itr_scale-1));
        cv::Mat imdist_scaled; resize(orig_bw, imdist_scaled, siz, 0, 0, cv::INTER_CUBIC);

        cv::Mat structdis;
        //compute mu and mu squared
        cv::Mat mu,mu_sq;
        cv::GaussianBlur( imdist_scaled, mu, cv::Size(7,7), CV_GAUSSIAN, 1.16666 );
        cv::multiply(mu, mu, mu_sq);

        //compute sigma
        cv::Mat sigma;
        cv::multiply(imdist_scaled, imdist_scaled, sigma);
        cv::GaussianBlur(sigma, sigma, cv::Size(7,7), CV_GAUSSIAN, 1.16666 );
        cv::subtract(sigma, mu_sq, sigma);
        cv::pow(sigma, 0.5, sigma);

        //compute structdis = (x-mu)/sigma
        cv::add(sigma, cv::Scalar(1.0/255), sigma);
        cv::subtract(imdist_scaled, mu, structdis);
        cv::divide(structdis, sigma, structdis);

        //Compute AGGD fit
        double lsigma_best, rsigma_best, gamma_best;
        AGGDfit(structdis, lsigma_best, rsigma_best, gamma_best);
        featurevector.push_back(gamma_best);
        featurevector.push_back((lsigma_best*lsigma_best + rsigma_best*rsigma_best)/2);

        cv::Rect bounds(0, 0, structdis.cols, structdis.rows);
        //Compute paired product images
        int shifts[4][2]={{0,1},{1,0},{1,1},{0,0}};
        for (int itr_shift=0; itr_shift<4; itr_shift++)
        {
            int *rshift = shifts[itr_shift];
            cv::Rect r(rshift[0], rshift[1], structdis.cols - rshift[0], structdis.rows - rshift[1]);

            cv::Mat shiftarr(structdis.size(), structdis.type(), 0.0f);
            shiftarr(r & bounds) = 1;
            //computing correlation
            cv::Mat shifted;
            cv::multiply(structdis, shiftarr, shifted);

            AGGDfit(shifted, lsigma_best, rsigma_best, gamma_best);

            double constant = sqrt(tgamma(1/gamma_best))/sqrt(tgamma(3/gamma_best));
            double meanparam = (rsigma_best-lsigma_best)*(tgamma(2/gamma_best)/tgamma(1/gamma_best))*constant;

            featurevector.push_back(gamma_best);
            featurevector.push_back(meanparam);
            featurevector.push_back(pow(lsigma_best,2));
            featurevector.push_back(pow(rsigma_best,2));
        }
    }
}