function score = jpeg_quality_score(img)

%========================================================================
%
%Copyright (c) 2002 The University of Texas at Austin
%All Rights Reserved.
% 
%This program is free software; you can redistribute it and/or modify
%it under the terms of the GNU General Public License as published by
%the Free Software Foundation; either version 2 of the License, or
%(at your option) any later version.
% 
%This program is distributed in the hope that it will be useful,
%but WITHOUT ANY WARRANTY; without even the implied warranty of
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%GNU General Public License for more details.
% 
%The GNU Public License is available in the file LICENSE, or you
%can write to the Free Software Foundation, Inc., 59 Temple Place -
%Suite 330, Boston, MA 02111-1307, USA, or you can find it on the
%World Wide Web at http://www.fsf.org.
%
%Author  : Zhou Wang 
%Version : 1.0
% 
%The authors are with the Laboratory for Image and Video Engineering
%(LIVE), Department of Electrical and Computer Engineering, The
%University of Texas at Austin, Austin, TX.
%
%Kindly report any suggestions or corrections to zhouwang@ieee.org
%
%========================================================================
%
%This is an implementation of the algorithm for calculating the quality
%score of JPEG compressed images proposed by Zhou Wang, Hamid R. Sheikh
%and Alan C. Bovik. Please refer to the paper: Zhou Wang, Hamid R. Sheikh
%and Alan C. Bovik, "No-Reference Perceptual Quality Assessment of JPEG
%Compressed Images," submitted to IEEE International Conference on Image
%Processing, Sept. 2002.
%
%You can change this program as you like and use it anywhere, but please
%refer to its original source (cite our paper and our web page at
%http://anchovy.ece.utexas.edu/~zwang/research/nr_jpeg_quality/index.html).
%
%Input : A test 8bits/pixel grayscale image loaded in a 2-D array
%Output: A quality score of the image. The score typically has a value
%        between 1 and 10 (10 represents the best quality, 1 the worst).
%
%Usage:
%
%1. Load the image, for example
%
%   image = imread(‘testimage.jpg‘); 
%
%2. Call this function to calculate the quality score:
%
%   Quality_Score = jpeg_quality_score(image)
%
%========================================================================
% img=imread(‘/Users/anitafang/Desktop/testimg/test0.jpg‘);

if (nargin > 1)
    score = -1;
    return;
end

%  [M,N] = size(img);
%  if (M < 16 | N < 16)
%      score = -2;
%      return;
%  end
 
x1= rgb2gray(img);
x = double(x1);
[M,N] = size(x);

% Feature Extraction:

% 1. horizontal features

d_h = x(:, 2:N) - x(:, 1:(N-1));
% [m, n]=size(d_h);
%  disp(d_h);
%  fprintf(‘img width %d,and height %d,length %d, dims %d\n‘,M,N,length(img),ndims(img));
%  fprintf(‘d_h width %d,and height %d, length %d, dims %d\n‘,m,n,length(d_h),ndims(d_h));



B_h = mean2(abs(d_h(:, 8:8:8*(floor(N/8)-1))));

A_h = (8*mean2(abs(d_h)) - B_h)/7;

sig_h = sign(d_h);
left_sig = sig_h(:, 1:(N-2));
right_sig = sig_h(:, 2:(N-1));
Z_h = mean2((left_sig.*right_sig)<0);

%  fprintf(‘B_h:%f,A_h:%f,Z_h:%f,\n‘,B_h,A_h,Z_h);

% 2. vertical features

d_v = x(2:M, :) - x(1:(M-1), :);

B_v = mean2(abs(d_v(8:8:8*(floor(M/8)-1), :)));

A_v = (8*mean2(abs(d_v)) - B_v)/7;

sig_v = sign(d_v);
up_sig = sig_v(1:(M-2), :);
down_sig = sig_v(2:(M-1), :);
Z_v = mean2((up_sig.*down_sig)<0);

% 3. combined features

B = (B_h + B_v)/2;
A = (A_h + A_v)/2;
Z = (Z_h + Z_v)/2;

% Quality Prediction

alpha = -245.8909;
beta = 261.9373;
gamma1 = -239.8886;
gamma2 = 160.1664;
gamma3 = 64.2859;

score = alpha + beta*(B.^(gamma1/10000))*(A.^(gamma2/10000))*(Z.^(gamma3/10000));

function main()

allNum = 0;
blurNum = 0;

threshold = 7.7;


for  k = 0:2000
    try
        image = imread([‘/Users/anitafang/Desktop/testimg/‘,‘test‘,num2str(k),‘.jpg‘]);
        %image = imread([‘/Users/user/Desktop/test/ye_blur/‘,num2str(k),‘.jpg‘]);
        quality_score = jpeg_quality_score(image);
        fprintf(‘Quality_Score:%d.jpg  %f\n‘,k,quality_score);
        allNum = allNum + 1;
        
        if quality_score < threshold
           blurNum = blurNum + 1;
           savePath = [‘/Users/anitafang/Desktop/blur/‘,num2str(k),‘.jpg‘];
           imwrite(image,savePath);
        end
%         break
    catch err
        %throw(err);
    end
end

% fprintf(‘relevance:%f\n‘,blurNum/allNum)

//
//  jpegquality.hpp
//  SDM_Train
//
//  Created by anitafang on 2017/6/27.
//  Copyright ? Anita，fang. All rights reserved.
//

#ifndef jpegquality_hpp
#define jpegquality_hpp

#include <stdio.h>
#include <vector>
#include <iostream>
#include <time.h>
#include <fstream>
#include <math.h> 
#include <cmath>

#include "opencv2/opencv.hpp"
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/objdetect/objdetect.hpp"

#endif /* jpegquality_hpp */

using namespace std;
using namespace cv;

class JPEGQ{
public:
    JPEGQ();
    void h_feature(Mat x,int M,int N);
    void v_feature(Mat x,int M,int N);
    void combine_feature();
    double qual_predict();
    
private:
    
    double B_h,A_h,Z_h;
    double B_v,A_v,Z_v;
    double B,A,Z;
};

//
//  jpegquality.cpp
//  SDM_Train
//
//  Created by anitafang on 2017/6/27.
//  Copyright ? 2017年 antia.fang  All rights reserved.
//

#include "jpegquality.hpp"


using namespace std;
using namespace cv;


JPEGQ::JPEGQ(){
    B_h = 0;
    A_h = 0;
    Z_h = 0;
    B_v = 0;
    A_v = 0;
    Z_v = 0;
    B = 0;
    A = 0;
    Z = 0;
}
  //% 1. horizontal features
void JPEGQ::h_feature(Mat x,int M,int N){
   
    Mat d_h(M,N-1,CV_64F, Scalar(0,0,255));
    
    for(int i=0;i<M;i++){
        for(int j=0;j<N-1;j++){
            d_h.at<double_t>(i,j)=x.at<double_t>(i,j+1)-x.at<double_t>(i,j);
        }
    }
  //  cout << "d_h = " << d_h << "\n";
    int DEL= floor(N/8);
    Mat d1_h(M,DEL,CV_64F, Scalar(0,0,255));
    
    for (int i=0;i<M;i++){
        for(int j=0;j<DEL;j++){
            d1_h.at<double_t>(i,j)= abs(d_h.at<double_t>(i,8*j+7));
        }
    }
    
    //求矩阵像素的的平均值
    B_h = mean(d1_h)[0];
    A_h = (8*mean(abs(d_h))[0] - B_h)/7;
    
    Mat sig_h(M,N-2,CV_64F, Scalar(0,0,255));
  
    for(int i=0;i<M;i++){
        for(int j=0;j<N-1;j++){
            if(d_h.at<double_t>(i,j) < 0){
                sig_h.at<double_t>(i,j) = -1;
            }
            else if(d_h.at<double_t>(i,j) > 0){
                sig_h.at<double_t>(i,j) = 1;
            }
            else {
                sig_h.at<double_t>(i,j) = 0;
            }
        }
    }
    
    Mat left_sig(M,N-2,CV_64F,Scalar(0,0,255));
    Mat right_sig(M,N-2,CV_64F,Scalar(0,0,255));
    
    for(int i=0;i<M;i++){
        for(int j=0;j<N-2;j++){
            left_sig.at<double_t>(i,j)=sig_h.at<double_t>(i,j);
            right_sig.at<double_t>(i,j)=sig_h.at<double_t>(i,j+1);
        }
    }
    
    
    
    //double Z_h = mean2((left_sig.*right_sig)<0);
    Mat multi_sig(M,N-2,CV_64F, Scalar(0,0,255));
    for(int i=0;i<M;i++){
        for(int j=0;j<N-2;j++){
            double temp1=left_sig.at<double_t>(i,j)* right_sig.at<double_t>(i,j);
            if(temp1<0){
                multi_sig.at<double_t>(i,j)= 1;
            }
            else {
                multi_sig.at<double_t>(i,j)= 0;
            }
            
        }
    }
    Z_h =mean(multi_sig)[0];
    
   //  cout <<"B_h:  "<< B_h<<"A_h:  "<< A_h<<"Z_h:  "<< Z_h << endl;
    
}

    
 //   % 2. vertical features
void JPEGQ::v_feature(Mat x,int M,int N){
    
    Mat d_v(M-1,N,CV_64F, Scalar(0,0,255));
    
    for(int i=0;i<M-1;i++){
        for(int j=0;j<N;j++){
            
            d_v.at<double_t>(i,j)=x.at<double_t>(i+1,j)-x.at<double_t>(i,j);
            
        }
    }
    
    int DELV= floor(M/8);
    Mat d1_v(DELV,N,CV_64F, Scalar(0,0,255));
    
    for (int i=0;i<DELV;i++){
        for(int j=0;j<N;j++){
            
            d1_v.at<double_t>(i,j)= abs(d_v.at<double_t>(8*i+7,j));
        }
    }
    
    //求矩阵像素的的平均值
    
    B_v=mean(d1_v)[0];
    A_v = (8*mean(abs(d_v))[0] - B_v)/7;
    
    
    Mat sig_v(M-1,N,CV_64F, Scalar(0,0,255));
    
    for(int i=0;i<M-1;i++){
        
        for(int j=0;j<N;j++){
            
            if(d_v.at<double_t>(i,j)<0)
            {  sig_v.at<double_t>(i,j)=-1;  }
            
            else if(d_v.at<double_t>(i,j) >0)
            {  sig_v.at<double_t>(i,j) = 1;  }
            
            else {  sig_v.at<double_t>(i,j) = 0;}
            
        }
    }
    
    Mat up_sig(M-2,N,CV_64F, Scalar(0,0,255));
    Mat down_sig(M-2,N,CV_64F, Scalar(0,0,255));
    
    for(int i=0;i<M-2;i++){
        for(int j=0;j<N;j++){
            
            up_sig.at<double_t>(i,j)=sig_v.at<double_t>(i,j);
            down_sig.at<double_t>(i,j)=sig_v.at<double_t>(i+1,j);
            
        }
    }
    
    //double Z_h = mean2((left_sig.*right_sig)<0);
    Mat vmulti_sig(M-2,N,CV_64F, Scalar(0,0,255));
    for(int i=0;i<M-2;i++){
        for(int j=0;j<N;j++){
            double temp2=up_sig.at<double_t>(i,j)* down_sig.at<double_t>(i,j);
            if(temp2<0)
            { vmulti_sig.at<double_t>(i,j)= 1;}
            else { vmulti_sig.at<double_t>(i,j)= 0; }
            
        }
    }
    Z_v =mean(vmulti_sig)[0];

    
}


//% 3. combined features
void JPEGQ::combine_feature(){
    
    B = (B_h + B_v)/2;
    A = (A_h + A_v)/2;
    Z = (Z_h + Z_v)/2;
    
}


//% Quality Prediction
double JPEGQ::qual_predict(){
    
    double alpha = -245.8909;
    double beta = 261.9373;
    double gamma1 = -239.8886;
    double gamma2 = 160.1664;
    double gamma3 = 64.2859;
    
    double score = alpha + beta*(pow(B,gamma1/10000)*pow(A,gamma2/10000)*pow(Z,gamma3/10000));
    
   
    return score;
}

//
//  main.cpp
//  jpg_quality
//
//  Created by anitafang on 2017/6/28.
//  Copyright ? 2017年 anitafang. All rights reserved.
//
#include <iostream>
#include "jpegquality.hpp"


using namespace std;
using namespace cv;


int main(int argc, const char * argv[]) {
    
    char filename[100];
    double threshold = 7.7;
    int num =12;
   // int *pia = new int[num] ();  // 每个元素初始化为0
    
    
    for(int k=0;k<num;k++){
        
        sprintf(filename,"/Users/anitafang/Desktop/VIP/xcode-demo/test-img/testimg/test%d.jpg",k);
    
        Mat x1=imread(filename,IMREAD_GRAYSCALE);
        
        int M=x1.rows;
        int N=x1.cols;
        
        Mat x;
        x1.convertTo(x, CV_64F);//转换成浮点运算
        
        JPEGQ *jpegq = new JPEGQ();
        jpegq->h_feature(x,M, N);
        jpegq->v_feature(x,M, N);
        jpegq->combine_feature();
        double score= jpegq->qual_predict();
        
        if(score<threshold){
            cout<<"this is a blur image"<<endl;
        }
    
        cout<<"the image :"<<k<<" "<<" score is :"<<score<<endl;
        
    }

    
    return 0;
}