sjero.net Git - qpsnr/blob - src/stats.cpp

   1 /*
   2 *       qpsnr (C) 2010 E. Oriani, ema <AT> fastwebnet <DOT> it
   3 *
   4 *       This file is part of qpsnr.
   5 *
   6 *       qpsnr is free software: you can redistribute it and/or modify
   7 *       it under the terms of the GNU General Public License as published by
   8 *       the Free Software Foundation, either version 3 of the License, or
   9 *       (at your option) any later version.
  10 *
  11 *       qpsnr is distributed in the hope that it will be useful,
  12 *       but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 *       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 *       GNU General Public License for more details.
  15 *
  16 *       You should have received a copy of the GNU General Public License
  17 *       along with qpsnr.  If not, see <http://www.gnu.org/licenses/>.
  18 */
  19
  20 #include "stats.h"
  21 #include "mt.h"
  22 #include "shared_ptr.h"
  23 #include "settings.h"
  24 #include <stdio.h>
  25 #include <climits>
  26 #include <string.h>
  27 #include <cmath>
  28 #include <string>
  29 #include <stdexcept>
  30 #include <set>
  31 #include <fstream>
  32 #include <algorithm>
  33 #include <cstdlib>
  34
  35 // define these classes just locally
  36 namespace stats {
  37         static double compute_psnr(const unsigned char *ref, const unsigned char *cmp, const unsigned int& sz) {
  38                 double mse = 0.0;
  39                 for(unsigned int i = 0; i < sz; ++i) {
  40                         const int       diff = ref[i]-cmp[i];
  41                         mse += (diff*diff);
  42                 }
  43                 mse /= (double)sz;
  44                 if (0.0 == mse) mse = 1e-10;
  45                 return 10.0*log10(65025.0/mse);
  46         }
  47
  48         static double compute_ssim(const unsigned char *ref, const unsigned char *cmp, const unsigned int& x, const unsigned int& y, const unsigned int& b_sz) {
  49                 // we return the average of all the blocks
  50                 const unsigned int      x_bl_num = x/b_sz,
  51                                         y_bl_num = y/b_sz;
  52                 if (!x_bl_num || !y_bl_num) return 0.0;
  53                 std::vector<double>     ssim_accum;
  54                 // for each block do it
  55                 for(int yB = 0; yB < y_bl_num; ++yB)
  56                         for(int xB = 0; xB < x_bl_num; ++xB) {
  57                                 const unsigned int      base_offset = xB*b_sz + yB*b_sz*x;
  58                                 double                  ref_acc = 0.0,
  59                                                         ref_acc_2 = 0.0,
  60                                                         cmp_acc = 0.0,
  61                                                         cmp_acc_2 = 0.0,
  62                                                         ref_cmp_acc = 0.0;
  63                                 for(int j = 0; j < b_sz; ++j)
  64                                         for(int i = 0; i < b_sz; ++i) {
  65                                                 // we have to multiply by 3, colorplanes are Y Cb Cr, we need
  66                                                 // only Y component
  67                                                 const unsigned char     c_ref = ref[3*(base_offset + j*x + i)],
  68                                                                         c_cmp = cmp[3*(base_offset + j*x + i)];
  69                                                 ref_acc += c_ref;
  70                                                 ref_acc_2 += (c_ref*c_ref);
  71                                                 cmp_acc += c_cmp;
  72                                                 cmp_acc_2 += (c_cmp*c_cmp);
  73                                                 ref_cmp_acc += (c_ref*c_cmp);
  74                                         }
  75                                 // now finally get the ssim for this block
  76                                 // http://en.wikipedia.org/wiki/SSIM
  77                                 // http://en.wikipedia.org/wiki/Variance
  78                                 // http://en.wikipedia.org/wiki/Covariance
  79                                 const double    n_samples = (b_sz*b_sz),
  80                                                 ref_avg = ref_acc/n_samples,
  81                                                 ref_var = ref_acc_2/n_samples - (ref_avg*ref_avg),
  82                                                 cmp_avg = cmp_acc/n_samples,
  83                                                 cmp_var = cmp_acc_2/n_samples - (cmp_avg*cmp_avg),
  84                                                 ref_cmp_cov = ref_cmp_acc/n_samples - (ref_avg*cmp_avg),
  85                                                 c1 = 6.5025, // (0.01*255.0)^2
  86                                                 c2 = 58.5225, // (0.03*255)^2
  87                                                 ssim_num = (2.0*ref_avg*cmp_avg + c1)*(2.0*ref_cmp_cov + c2),
  88                                                 ssim_den = (ref_avg*ref_avg + cmp_avg*cmp_avg + c1)*(ref_var + cmp_var + c2),
  89                                                 ssim = ssim_num/ssim_den;
  90                                 ssim_accum.push_back(ssim);
  91                         }
  92
  93                 double  avg = 0.0;
  94                 for(std::vector<double>::const_iterator it = ssim_accum.begin(); it != ssim_accum.end(); ++it)
  95                         avg += *it;
  96                 return avg/ssim_accum.size();
  97         }
  98
  99         double std_dev(int *buffer, int numelem)
 100         {
 101                 // find sum and mean
 102                 long sum = 0;
 103                 for (int i=0; i<numelem; i++) sum+=buffer[i];
 104                 int mean = sum/numelem;
 105
 106                 // sum of the squares
 107                 long long sum_squares = 0;
 108                 for (int i=0; i<numelem; i++)
 109                         sum_squares += (long long)pow((float)(buffer[i]-mean),2);
 110
 111                 // return the square root
 112                 return sqrt((double)(sum_squares/numelem));
 113         }
 114
 115         void conv_sobel(int *buffer, int width, int height)
 116         {
 117                 int mask1[9] = {-1,-2,-1,0,0,0,1,2,1}; // TODO right values
 118                 int mask2[9] = {-1,0,1,-2,0,2,-1,0,1};
 119
 120                 // convolution
 121                 int buffsize = width*height;
 122                 int *buffer_tmp = (int *)malloc(buffsize*sizeof(int));
 123                 memcpy(buffer_tmp, buffer, buffsize*sizeof(int));
 124                 for (int l=1; l<height-1; l++) {
 125                 for (int c=1; c<width-1; c++) {
 126
 127                         int p1 = (l-1)*width+(c-1);     // pos first pixel
 128                         int p2 = l*width+(c-1);         // pos 2nd line first pixel
 129                         int p3 = (l+1)*width+(c-1);     // pos 3nd line first pixel
 130
 131                         // apply the mask and set the result
 132                         long gv = 0, gh = 0;
 133                         gv += buffer[p1]*mask1[0] + buffer[p1+1]*mask1[1] + buffer[p1+2]*mask1[2];
 134                         gv += buffer[p2]*mask1[3] + buffer[p2+1]*mask1[4] + buffer[p2+2]*mask1[5];
 135                         gv += buffer[p3]*mask1[6] + buffer[p3+1]*mask1[7] + buffer[p3+2]*mask1[8];
 136                         gh += buffer[p1]*mask2[0] + buffer[p1+1]*mask2[1] + buffer[p1+2]*mask2[2];
 137                         gh += buffer[p2]*mask2[3] + buffer[p2+1]*mask2[4] + buffer[p2+2]*mask2[5];
 138                         gh += buffer[p3]*mask2[6] + buffer[p3+1]*mask2[7] + buffer[p3+2]*mask2[8];
 139
 140                         buffer_tmp[l*width+c] = sqrt(pow((double)gv,2)+pow((double)gh,2));
 141                 }
 142                 }
 143                 memcpy(buffer, buffer_tmp, buffsize*sizeof(int));
 144                 free(buffer_tmp);
 145         }
 146
 147         static double compute_si(const unsigned char *ref, const unsigned char *cmp, const unsigned int& x, const unsigned int& y){
 148                 int sizebuff = x*y;
 149                 int *buff_sob1 = (int *)malloc(sizebuff*sizeof(int));
 150                 char c;
 151
 152                 // create new buffer and apply both sobel filters
 153                 for (int i=0; i<sizebuff; i++)
 154                         buff_sob1[i] = cmp[3*i];
 155
 156                 conv_sobel(buff_sob1, x, y);
 157
 158                 // get the standard deviation
 159                 double sdev = std_dev(buff_sob1, sizebuff);
 160                 return sdev;
 161         }
 162
 163         static double compute_ti(const unsigned char *ref, const unsigned char *cmp, const unsigned int& x, const unsigned int& y){
 164                 int numelem = x*y;
 165                 int *buff_tmp = (int *)malloc(numelem*sizeof(int));
 166
 167                 if(ref==NULL){
 168                         return 0;
 169                 }
 170
 171                 // difference between 2nd frame and 1st one to buff_tmp
 172                 for (int i=0; i<numelem; i++) {
 173                         buff_tmp[i] = (int)cmp[3*i]-(int)ref[3*i];
 174                 }
 175
 176                 // find and return the standard deviation of buff_tmp
 177                 double sdev = std_dev(buff_tmp, numelem);
 178                 free(buff_tmp);
 179                 return sdev;
 180         }
 181
 182         static inline double r_0_1(const double& d) {
 183                 return std::max(0.0, std::min(1.0, d));
 184         }
 185
 186         static void rgb_2_hsi(unsigned char *p, const int& sz) {
 187                 /*
 188                 I = (1/3) *(R+G+B)
 189                 S = 1 - ( (3/(R+G+B)) * min(R,G,B))
 190                 H = cos^-1( ((1/2) * ((R-G) + (R - B))) / ((R-G)^2 + (R-B)*(G-B)) ^(1/2) )
 191                 H = cos^-1 ( (((R-G)+(R-B))/2)/ (sqrt((R-G)^2 + (R-B)*(G-B) )))
 192                 */
 193                 const static double PI = 3.14159265;
 194                 for(int j =0; j < sz; j += 3) {
 195                         const double    r = p[j+0]/255.0,
 196                                         g = p[j+1]/255.0,
 197                                         b = p[j+2]/255.0,
 198                                         i = r_0_1((r+g+b)/3),
 199                                         s = r_0_1(1.0 - (3.0/(r+g+b) * std::min(r, std::min(g, b)))),
 200                                         h = r_0_1(acos(0.5*(r-g + r-b) / sqrt((r-g)*(r-g) + (r-b)*(g-b)))/PI);
 201                         p[j+0] = 255.0*h + 0.5;
 202                         p[j+1] = 255.0*s + 0.5;
 203                         p[j+2] = 255.0*i + 0.5;
 204                 }
 205         }
 206
 207         static void rgb_2_YCbCr(unsigned char *p, const int& sz) {
 208                 // http://en.wikipedia.org/wiki/YCbCr
 209                 for(int j =0; j < sz; j += 3) {
 210                         const double    r = p[j+0]/255.0,
 211                                         g = p[j+1]/255.0,
 212                                         b = p[j+2]/255.0;
 213                         p[j+0] = 16 + (65.481*r + 128.553*g + 24.966*b);
 214                         p[j+1] = 128 + (-37.797*r - 74.203*g + 112.0*b);
 215                         p[j+2] = 128 + (112.0*r - 93.786*g - 12.214*b);
 216                 }
 217         }
 218
 219         static void rgb_2_Y(unsigned char *p, const int& sz) {
 220                 // http://en.wikipedia.org/wiki/YCbCr
 221                 for(int j =0; j < sz; j += 3) {
 222                         const double    r = p[j+0]/255.0,
 223                                         g = p[j+1]/255.0,
 224                                         b = p[j+2]/255.0;
 225                         p[j+0] = 16 + (65.481*r + 128.553*g + 24.966*b);
 226                         p[j+1] = 0.0;
 227                         p[j+2] = 0.0;
 228                 }
 229         }
 230
 231         static int getCPUcount(void) {
 232                 std::ifstream           cpuinfo("/proc/cpuinfo");
 233                 std::string             line;
 234                 std::set<std::string>   IDs;
 235                 while (cpuinfo){
 236                         std::getline(cpuinfo,line);
 237                         if (line.empty())
 238                                 continue;
 239                         if (line.find("processor") != 0)
 240                                 continue;
 241                         IDs.insert(line);
 242                 }
 243                 if (IDs.empty()) return 1;
 244                 return IDs.size();
 245         }
 246
 247         static mt::ThreadPool   __stats_tp(getCPUcount());
 248
 249         class psnr_job : public mt::ThreadPool::Job {
 250                 const unsigned char     *_ref,
 251                                         *_cmp;
 252                 const unsigned int      _sz;
 253                 double                  &_res;
 254         public:
 255                 psnr_job(const unsigned char *ref, const unsigned char *cmp, const unsigned int& sz, double& res) :
 256                 _ref(ref), _cmp(cmp), _sz(sz), _res(res) {
 257                 }
 258
 259                 virtual void run(void) {
 260                         _res = compute_psnr(_ref, _cmp, _sz);
 261                 }
 262         };
 263
 264         static void get_psnr_tp(const VUCHAR& ref, const std::vector<bool>& v_ok, const std::vector<VUCHAR>& streams, std::vector<double>& res) {
 265                 const unsigned int                      sz = v_ok.size();
 266                 std::vector<shared_ptr<psnr_job> >      v_jobs;
 267                 for(unsigned int i =0; i < sz; ++i) {
 268                         if (v_ok[i]) {
 269                                 v_jobs.push_back(new psnr_job(&ref[0], &(streams[i][0]), ref.size(), res[i]));
 270                                 __stats_tp.add(v_jobs.rbegin()->get());
 271                         } else res[i] = 0.0;
 272                 }
 273                 //wait for all
 274                 for(std::vector<shared_ptr<psnr_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 275                         (*it)->wait();
 276                         (*it) = 0;
 277                 }
 278         }
 279
 280         class ssim_job : public mt::ThreadPool::Job {
 281                 const unsigned char     *_ref,
 282                                         *_cmp;
 283                 const unsigned int      _x,
 284                                         _y,
 285                                         _b_sz;
 286                 double                  &_res;
 287         public:
 288                 ssim_job(const unsigned char *ref, const unsigned char *cmp, const unsigned int& x, const unsigned int& y, const unsigned int b_sz, double& res) :
 289                 _ref(ref), _cmp(cmp), _x(x), _y(y), _b_sz(b_sz), _res(res) {
 290                 }
 291
 292                 virtual void run(void) {
 293                         _res = compute_ssim(_ref, _cmp, _x, _y, _b_sz);
 294                 }
 295         };
 296
 297         static void get_ssim_tp(const VUCHAR& ref, const std::vector<bool>& v_ok, const std::vector<VUCHAR>& streams, std::vector<double>& res, const unsigned int& x, const unsigned int& y, const unsigned int& b_sz) {
 298                 const unsigned int                      sz = v_ok.size();
 299                 std::vector<shared_ptr<ssim_job> >      v_jobs;
 300                 for(unsigned int i =0; i < sz; ++i) {
 301                         if (v_ok[i]) {
 302                                 v_jobs.push_back(new ssim_job(&ref[0], &(streams[i][0]), x, y, b_sz, res[i]));
 303                                 __stats_tp.add(v_jobs.rbegin()->get());
 304                         } else res[i] = 0.0;
 305                 }
 306                 //wait for all
 307                 for(std::vector<shared_ptr<ssim_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 308                         (*it)->wait();
 309                         (*it) = 0;
 310                 }
 311         }
 312
 313         class siti_job : public mt::ThreadPool::Job {
 314                         const unsigned char     *_ref,
 315                                                 *_cmp;
 316                         const unsigned int      _x,
 317                                                 _y;
 318                         std::pair<double,double>                        &_res;
 319                 public:
 320                         siti_job(const unsigned char *ref, const unsigned char *cmp, const unsigned int& x, const unsigned int& y, std::pair<double,double>& res) :
 321                         _ref(ref), _cmp(cmp), _x(x), _y(y), _res(res) {
 322                         }
 323
 324                         virtual void run(void) {
 325                                 _res.first = compute_si(_ref, _cmp, _x, _y);
 326                                 _res.second= compute_ti(_ref, _cmp, _x, _y);
 327                         }
 328                 };
 329
 330         static void get_siti_tp(const VUCHAR* ref, const std::vector<bool>& v_ok, const std::vector<VUCHAR>& streams, std::vector<std::pair<double,double> >& res, const unsigned int& x, const unsigned int& y) {
 331                         const unsigned int                      sz = v_ok.size();
 332                         std::vector<shared_ptr<siti_job> >      v_jobs;
 333                         for(unsigned int i =0; i < sz; ++i) {
 334                                 if (v_ok[i]) {
 335                                         v_jobs.push_back(new siti_job(&((*ref)[0]), &(streams[i][0]), x, y, res[i]));
 336                                         __stats_tp.add(v_jobs.rbegin()->get());
 337                                 } else{
 338                                         res[i].first = 0.0;
 339                                         res[i].second=0.0;
 340                                 }
 341                         }
 342                         //wait for all
 343                         for(std::vector<shared_ptr<siti_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 344                                 (*it)->wait();
 345                                 (*it) = 0;
 346                         }
 347                 }
 348
 349         class hsi_job : public mt::ThreadPool::Job {
 350                 unsigned char           *_buf;
 351                 const unsigned int      _sz;
 352         public:
 353                 hsi_job(unsigned char *buf, const unsigned int& sz) :
 354                 _buf(buf), _sz(sz) {
 355                 }
 356
 357                 virtual void run(void) {
 358                         rgb_2_hsi(_buf, _sz);
 359                 }
 360         };
 361
 362         class YCbCr_job : public mt::ThreadPool::Job {
 363                 unsigned char           *_buf;
 364                 const unsigned int      _sz;
 365         public:
 366                 YCbCr_job(unsigned char *buf, const unsigned int& sz) :
 367                 _buf(buf), _sz(sz) {
 368                 }
 369
 370                 virtual void run(void) {
 371                         rgb_2_YCbCr(_buf, _sz);
 372                 }
 373         };
 374
 375         class Y_job : public mt::ThreadPool::Job {
 376                 unsigned char           *_buf;
 377                 const unsigned int      _sz;
 378         public:
 379                 Y_job(unsigned char *buf, const unsigned int& sz) :
 380                 _buf(buf), _sz(sz) {
 381                 }
 382
 383                 virtual void run(void) {
 384                         rgb_2_Y(_buf, _sz);
 385                 }
 386         };
 387
 388         static void rgb_2_hsi_tp(VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 389                 const unsigned int                      sz = v_ok.size();
 390                 std::vector<shared_ptr<hsi_job> >       v_jobs;
 391                 v_jobs.push_back(new hsi_job(&ref[0], ref.size()));
 392                 __stats_tp.add(v_jobs.rbegin()->get());
 393                 for(unsigned int i =0; i < sz; ++i) {
 394                         if (v_ok[i]) {
 395                                 v_jobs.push_back(new hsi_job(&(streams[i][0]), streams[i].size()));
 396                                 __stats_tp.add(v_jobs.rbegin()->get());
 397                         }
 398                 }
 399                 //wait for all
 400                 for(std::vector<shared_ptr<hsi_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 401                         (*it)->wait();
 402                         (*it) = 0;
 403                 }
 404         }
 405
 406         static void rgb_2_YCbCr_tp(VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 407                 const unsigned int                      sz = v_ok.size();
 408                 std::vector<shared_ptr<YCbCr_job> >     v_jobs;
 409                 v_jobs.push_back(new YCbCr_job(&ref[0], ref.size()));
 410                 __stats_tp.add(v_jobs.rbegin()->get());
 411                 for(unsigned int i =0; i < sz; ++i) {
 412                         if (v_ok[i]) {
 413                                 v_jobs.push_back(new YCbCr_job(&(streams[i][0]), streams[i].size()));
 414                                 __stats_tp.add(v_jobs.rbegin()->get());
 415                         }
 416                 }
 417                 //wait for all
 418                 for(std::vector<shared_ptr<YCbCr_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 419                         (*it)->wait();
 420                         (*it) = 0;
 421                 }
 422         }
 423
 424         static void rgb_2_Y_tp(VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 425                 const unsigned int                      sz = v_ok.size();
 426                 std::vector<shared_ptr<Y_job> >         v_jobs;
 427                 v_jobs.push_back(new Y_job(&ref[0], ref.size()));
 428                 __stats_tp.add(v_jobs.rbegin()->get());
 429                 for(unsigned int i =0; i < sz; ++i) {
 430                         if (v_ok[i]) {
 431                                 v_jobs.push_back(new Y_job(&(streams[i][0]), streams[i].size()));
 432                                 __stats_tp.add(v_jobs.rbegin()->get());
 433                         }
 434                 }
 435                 //wait for all
 436                 for(std::vector<shared_ptr<Y_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 437                         (*it)->wait();
 438                         (*it) = 0;
 439                 }
 440         }
 441
 442         class psnr : public s_base {
 443                 std::string     _colorspace;
 444         protected:
 445                 void print(const int& ref_frame, const std::vector<double>& v_res) {
 446                         _ostr << ref_frame << ',';
 447                         for(int i = 0; i < _n_streams; ++i)
 448                                 _ostr << v_res[i] << ',';
 449                         _ostr << std::endl;
 450                 }
 451
 452                 void process_colorspace(VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 453                         if (_colorspace == "hsi") {
 454                                 rgb_2_hsi_tp(ref, v_ok, streams);
 455                         } else if (_colorspace == "ycbcr") {
 456                                 rgb_2_YCbCr_tp(ref, v_ok, streams);
 457                         } else if (_colorspace == "y") {
 458                                 rgb_2_Y_tp(ref, v_ok, streams);
 459                         }
 460                 }
 461         public:
 462                 psnr(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 463                 s_base(n_streams, i_width, i_height, ostr), _colorspace("rgb") {
 464                 }
 465
 466                 virtual void set_parameter(const char* p_name, const char *p_value) {
 467                         const std::string s_name(p_name);
 468                         if (s_name == "colorspace") {
 469                                 _colorspace = p_value;
 470                                 if (_colorspace != "rgb" && _colorspace != "hsi"
 471                                 && _colorspace != "ycbcr" && _colorspace != "y")
 472                                         throw std::runtime_error("Invalid colorspace passed to analyzer");
 473                         }
 474                 }
 475
 476                 virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 477                         if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 478                         // process colorspace
 479                         process_colorspace(ref, v_ok, streams);
 480                         //
 481                         std::vector<double>     v_res(_n_streams);
 482                         get_psnr_tp(ref, v_ok, streams, v_res);
 483                         //
 484                         print(ref_frame, v_res);
 485                 }
 486         };
 487
 488         class avg_psnr : public psnr {
 489                 int                     _fpa,
 490                                         _accum_f,
 491                                         _last_frame;
 492                 std::vector<double>     _accum_v;
 493         protected:
 494                 void print(const int& ref_frame) {
 495                         // check if we have to print infor or not
 496                         if (_accum_f < _fpa) return;
 497                         // else print data
 498                         for(int i=0; i < _n_streams; ++i)
 499                                 _accum_v[i] /= _accum_f;
 500                         psnr::print(ref_frame, _accum_v);
 501                         // clear the vector and set _accum_f to 0
 502                         _accum_v.clear();
 503                         _accum_v.resize(_n_streams);
 504                         _accum_f = 0;
 505                 }
 506         public:
 507                 avg_psnr(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 508                 psnr(n_streams, i_width, i_height, ostr), _fpa(1), _accum_f(0), _last_frame(-1), _accum_v(n_streams) {
 509                 }
 510
 511                 virtual void set_parameter(const char* p_name, const char *p_value) {
 512                         psnr::set_parameter(p_name, p_value);
 513                         //
 514                         const std::string s_name(p_name);
 515                         if (s_name == "fpa") {
 516                                 const int fpa = atoi(p_value);
 517                                 if (fpa > 0) _fpa = fpa;
 518                         }
 519                 }
 520
 521                 virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 522                         if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 523                         // set last frame
 524                         _last_frame = ref_frame;
 525                         // process colorspace
 526                         process_colorspace(ref, v_ok, streams);
 527                         // compute the psnr
 528                         std::vector<double>     v_res(_n_streams);
 529                         get_psnr_tp(ref, v_ok, streams, v_res);
 530                         // accumulate for each
 531                         for(int i = 0; i < _n_streams; ++i) {
 532                                 if (v_ok[i]) {
 533                                          _accum_v[i] += v_res[i];
 534                                 } else _accum_v[i] += 0.0;
 535                         }
 536                         ++_accum_f;
 537                         // try to print data
 538                         print(ref_frame);
 539                 }
 540
 541                 virtual ~avg_psnr() {
 542                         // on exit check if we have some data
 543                         if (_accum_f > 0) {
 544                                 _ostr << _last_frame << ',';
 545                                 for(int i = 0; i < _n_streams; ++i) {
 546                                         _ostr << _accum_v[i]/_accum_f << ',';
 547                                 }
 548                                 _ostr << std::endl;
 549                         }
 550                 }
 551         };
 552
 553         class ssim : public s_base {
 554         protected:
 555                 int     _blocksize;
 556
 557                 void print(const int& ref_frame, const std::vector<double>& v_res) {
 558                         _ostr << ref_frame << ',';
 559                         for(int i = 0; i < _n_streams; ++i)
 560                                 _ostr << v_res[i] << ',';
 561                         _ostr << std::endl;
 562                 }
 563         public:
 564                 ssim(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 565                 s_base(n_streams, i_width, i_height, ostr), _blocksize(8) {
 566                 }
 567
 568                 virtual void set_parameter(const char* p_name, const char *p_value) {
 569                         const std::string s_name(p_name);
 570                         if (s_name == "blocksize") {
 571                                 const int blocksize = atoi(p_value);
 572                                 if (blocksize > 0) _blocksize = blocksize;
 573                         }
 574                 }
 575
 576                 virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 577                         if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 578                         // convert to Y colorspace
 579                         rgb_2_Y_tp(ref, v_ok, streams);
 580                         //
 581                         std::vector<double>     v_res(_n_streams);
 582                         get_ssim_tp(ref, v_ok, streams, v_res, _i_width, _i_height, _blocksize);
 583                         //
 584                         print(ref_frame, v_res);
 585                 }
 586         };
 587
 588         class avg_ssim : public ssim {
 589                 int                     _fpa,
 590                                         _accum_f,
 591                                         _last_frame;
 592                 std::vector<double>     _accum_v;
 593         protected:
 594                 void print(const int& ref_frame) {
 595                         // check if we have to print infor or not
 596                         if (_accum_f < _fpa) return;
 597                         // else print data
 598                         for(int i=0; i < _n_streams; ++i)
 599                                 _accum_v[i] /= _accum_f;
 600                         ssim::print(ref_frame, _accum_v);
 601                         // clear the vector and set _accum_f to 0
 602                         _accum_v.clear();
 603                         _accum_v.resize(_n_streams);
 604                         _accum_f = 0;
 605                 }
 606         public:
 607                 avg_ssim(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 608                 ssim(n_streams, i_width, i_height, ostr), _fpa(1), _accum_f(0), _last_frame(-1), _accum_v(n_streams) {
 609                 }
 610
 611                 virtual void set_parameter(const char* p_name, const char *p_value) {
 612                         ssim::set_parameter(p_name, p_value);
 613                         //
 614                         const std::string s_name(p_name);
 615                         if (s_name == "fpa") {
 616                                 const int fpa = atoi(p_value);
 617                                 if (fpa > 0) _fpa = fpa;
 618                         }
 619                 }
 620
 621                 virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 622                         if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 623                         // set last frame
 624                         _last_frame = ref_frame;
 625                         // convert to Y colorspace
 626                         rgb_2_Y_tp(ref, v_ok, streams);
 627                         //
 628                         std::vector<double>     v_res(_n_streams);
 629                         get_ssim_tp(ref, v_ok, streams, v_res, _i_width, _i_height, _blocksize);
 630                         // accumulate for each
 631                         for(int i = 0; i < _n_streams; ++i) {
 632                                 if (v_ok[i]) {
 633                                          _accum_v[i] += v_res[i];
 634                                 } else _accum_v[i] += 0.0;
 635                         }
 636                         ++_accum_f;
 637                         // try to print data
 638                         print(ref_frame);
 639                 }
 640
 641                 virtual ~avg_ssim() {
 642                         // on exit check if we have some data
 643                         if (_accum_f > 0) {
 644                                 _ostr << _last_frame << ',';
 645                                 for(int i = 0; i < _n_streams; ++i) {
 646                                         _ostr << _accum_v[i]/_accum_f << ',';
 647                                 }
 648                                 _ostr << std::endl;
 649                         }
 650                 }
 651         };
 652
 653         class siti : public s_base {
 654                 protected:
 655                         VUCHAR* _prev;
 656                         void print(const int& ref_frame, const std::vector<std::pair<double,double> >& v_res) {
 657                                 _ostr << ref_frame << ',';
 658                                 for(int i = 0; i < _n_streams; ++i)
 659                                         _ostr << v_res[i].first << ',' << v_res[i].second << ",";
 660                                 _ostr << std::endl;
 661                         }
 662                 public:
 663                         siti(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 664                         s_base(n_streams, i_width, i_height, ostr){
 665                                 _prev=NULL;
 666                         }
 667
 668                         virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 669                                 if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 670                                 // convert to Y colorspace
 671                                 rgb_2_Y_tp(ref, v_ok, streams);
 672                                 //
 673                                 if(_prev==NULL){
 674                                         _prev=new std::vector<unsigned char> (streams[0]);
 675                                 }
 676                                 std::vector<std::pair<double,double> >  v_res(_n_streams);
 677                                 get_siti_tp(_prev, v_ok, streams, v_res, _i_width, _i_height);
 678                                 delete(_prev);
 679                                 _prev=new std::vector<unsigned char> (streams[0]);
 680                                 //
 681                                 print(ref_frame, v_res);
 682                         }
 683                 };
 684 }
 685
 686 stats::s_base* stats::get_analyzer(const char* id, const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) {
 687         const std::string       s_id(id);
 688         if (s_id == "psnr") return new psnr(n_streams, i_width, i_height, ostr);
 689         else if (s_id == "avg_psnr") return new avg_psnr(n_streams, i_width, i_height, ostr);
 690         else if (s_id == "ssim") return new ssim(n_streams, i_width, i_height, ostr);
 691         else if (s_id == "avg_ssim") return new avg_ssim(n_streams, i_width, i_height, ostr);
 692         else if (s_id == "siti") return new siti(n_streams, i_width, i_height, ostr);
 693         throw std::runtime_error("Invalid analyzer id");
 694 }