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                 free(buff_sob1);
 161                 return sdev;
 162         }
 163
 164         static double compute_ti(const unsigned char *ref, const unsigned char *cmp, const unsigned int& x, const unsigned int& y){
 165                 int numelem = x*y;
 166                 int *buff_tmp = (int *)malloc(numelem*sizeof(int));
 167
 168                 if(ref==NULL){
 169                         return 0;
 170                 }
 171
 172                 // difference between 2nd frame and 1st one to buff_tmp
 173                 for (int i=0; i<numelem; i++) {
 174                         buff_tmp[i] = (int)cmp[3*i]-(int)ref[3*i];
 175                 }
 176
 177                 // find and return the standard deviation of buff_tmp
 178                 double sdev = std_dev(buff_tmp, numelem);
 179                 free(buff_tmp);
 180                 return sdev;
 181         }
 182
 183         static inline double r_0_1(const double& d) {
 184                 return std::max(0.0, std::min(1.0, d));
 185         }
 186
 187         static void rgb_2_hsi(unsigned char *p, const int& sz) {
 188                 /*
 189                 I = (1/3) *(R+G+B)
 190                 S = 1 - ( (3/(R+G+B)) * min(R,G,B))
 191                 H = cos^-1( ((1/2) * ((R-G) + (R - B))) / ((R-G)^2 + (R-B)*(G-B)) ^(1/2) )
 192                 H = cos^-1 ( (((R-G)+(R-B))/2)/ (sqrt((R-G)^2 + (R-B)*(G-B) )))
 193                 */
 194                 const static double PI = 3.14159265;
 195                 for(int j =0; j < sz; j += 3) {
 196                         const double    r = p[j+0]/255.0,
 197                                         g = p[j+1]/255.0,
 198                                         b = p[j+2]/255.0,
 199                                         i = r_0_1((r+g+b)/3),
 200                                         s = r_0_1(1.0 - (3.0/(r+g+b) * std::min(r, std::min(g, b)))),
 201                                         h = r_0_1(acos(0.5*(r-g + r-b) / sqrt((r-g)*(r-g) + (r-b)*(g-b)))/PI);
 202                         p[j+0] = 255.0*h + 0.5;
 203                         p[j+1] = 255.0*s + 0.5;
 204                         p[j+2] = 255.0*i + 0.5;
 205                 }
 206         }
 207
 208         static void rgb_2_YCbCr(unsigned char *p, const int& sz) {
 209                 // http://en.wikipedia.org/wiki/YCbCr
 210                 for(int j =0; j < sz; j += 3) {
 211                         const double    r = p[j+0]/255.0,
 212                                         g = p[j+1]/255.0,
 213                                         b = p[j+2]/255.0;
 214                         p[j+0] = 16 + (65.481*r + 128.553*g + 24.966*b);
 215                         p[j+1] = 128 + (-37.797*r - 74.203*g + 112.0*b);
 216                         p[j+2] = 128 + (112.0*r - 93.786*g - 12.214*b);
 217                 }
 218         }
 219
 220         static void rgb_2_Y(unsigned char *p, const int& sz) {
 221                 // http://en.wikipedia.org/wiki/YCbCr
 222                 for(int j =0; j < sz; j += 3) {
 223                         const double    r = p[j+0]/255.0,
 224                                         g = p[j+1]/255.0,
 225                                         b = p[j+2]/255.0;
 226                         p[j+0] = 16 + (65.481*r + 128.553*g + 24.966*b);
 227                         p[j+1] = 0.0;
 228                         p[j+2] = 0.0;
 229                 }
 230         }
 231
 232         static int getCPUcount(void) {
 233                 std::ifstream           cpuinfo("/proc/cpuinfo");
 234                 std::string             line;
 235                 std::set<std::string>   IDs;
 236                 while (cpuinfo){
 237                         std::getline(cpuinfo,line);
 238                         if (line.empty())
 239                                 continue;
 240                         if (line.find("processor") != 0)
 241                                 continue;
 242                         IDs.insert(line);
 243                 }
 244                 if (IDs.empty()) return 1;
 245                 return IDs.size();
 246         }
 247
 248         static mt::ThreadPool   __stats_tp(getCPUcount());
 249
 250         class psnr_job : public mt::ThreadPool::Job {
 251                 const unsigned char     *_ref,
 252                                         *_cmp;
 253                 const unsigned int      _sz;
 254                 double                  &_res;
 255         public:
 256                 psnr_job(const unsigned char *ref, const unsigned char *cmp, const unsigned int& sz, double& res) :
 257                 _ref(ref), _cmp(cmp), _sz(sz), _res(res) {
 258                 }
 259
 260                 virtual void run(void) {
 261                         _res = compute_psnr(_ref, _cmp, _sz);
 262                 }
 263         };
 264
 265         static void get_psnr_tp(const VUCHAR& ref, const std::vector<bool>& v_ok, const std::vector<VUCHAR>& streams, std::vector<double>& res) {
 266                 const unsigned int                      sz = v_ok.size();
 267                 std::vector<shared_ptr<psnr_job> >      v_jobs;
 268                 for(unsigned int i =0; i < sz; ++i) {
 269                         if (v_ok[i]) {
 270                                 v_jobs.push_back(new psnr_job(&ref[0], &(streams[i][0]), ref.size(), res[i]));
 271                                 __stats_tp.add(v_jobs.rbegin()->get());
 272                         } else res[i] = 0.0;
 273                 }
 274                 //wait for all
 275                 for(std::vector<shared_ptr<psnr_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 276                         (*it)->wait();
 277                         (*it) = 0;
 278                 }
 279         }
 280
 281         class ssim_job : public mt::ThreadPool::Job {
 282                 const unsigned char     *_ref,
 283                                         *_cmp;
 284                 const unsigned int      _x,
 285                                         _y,
 286                                         _b_sz;
 287                 double                  &_res;
 288         public:
 289                 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) :
 290                 _ref(ref), _cmp(cmp), _x(x), _y(y), _b_sz(b_sz), _res(res) {
 291                 }
 292
 293                 virtual void run(void) {
 294                         _res = compute_ssim(_ref, _cmp, _x, _y, _b_sz);
 295                 }
 296         };
 297
 298         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) {
 299                 const unsigned int                      sz = v_ok.size();
 300                 std::vector<shared_ptr<ssim_job> >      v_jobs;
 301                 for(unsigned int i =0; i < sz; ++i) {
 302                         if (v_ok[i]) {
 303                                 v_jobs.push_back(new ssim_job(&ref[0], &(streams[i][0]), x, y, b_sz, res[i]));
 304                                 __stats_tp.add(v_jobs.rbegin()->get());
 305                         } else res[i] = 0.0;
 306                 }
 307                 //wait for all
 308                 for(std::vector<shared_ptr<ssim_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 309                         (*it)->wait();
 310                         (*it) = 0;
 311                 }
 312         }
 313
 314         class siti_job : public mt::ThreadPool::Job {
 315                         const unsigned char     *_ref,
 316                                                 *_cmp;
 317                         const unsigned int      _x,
 318                                                 _y;
 319                         std::pair<double,double>                        &_res;
 320                 public:
 321                         siti_job(const unsigned char *ref, const unsigned char *cmp, const unsigned int& x, const unsigned int& y, std::pair<double,double>& res) :
 322                         _ref(ref), _cmp(cmp), _x(x), _y(y), _res(res) {
 323                         }
 324
 325                         virtual void run(void) {
 326                                 _res.first = compute_si(_ref, _cmp, _x, _y);
 327                                 _res.second= compute_ti(_ref, _cmp, _x, _y);
 328                         }
 329                 };
 330
 331         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) {
 332                         const unsigned int                      sz = v_ok.size();
 333                         std::vector<shared_ptr<siti_job> >      v_jobs;
 334                         for(unsigned int i =0; i < sz; ++i) {
 335                                 if (v_ok[i]) {
 336                                         v_jobs.push_back(new siti_job(&((*ref)[0]), &(streams[i][0]), x, y, res[i]));
 337                                         __stats_tp.add(v_jobs.rbegin()->get());
 338                                 } else{
 339                                         res[i].first = 0.0;
 340                                         res[i].second=0.0;
 341                                 }
 342                         }
 343                         //wait for all
 344                         for(std::vector<shared_ptr<siti_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 345                                 (*it)->wait();
 346                                 (*it) = 0;
 347                         }
 348                 }
 349
 350         class hsi_job : public mt::ThreadPool::Job {
 351                 unsigned char           *_buf;
 352                 const unsigned int      _sz;
 353         public:
 354                 hsi_job(unsigned char *buf, const unsigned int& sz) :
 355                 _buf(buf), _sz(sz) {
 356                 }
 357
 358                 virtual void run(void) {
 359                         rgb_2_hsi(_buf, _sz);
 360                 }
 361         };
 362
 363         class YCbCr_job : public mt::ThreadPool::Job {
 364                 unsigned char           *_buf;
 365                 const unsigned int      _sz;
 366         public:
 367                 YCbCr_job(unsigned char *buf, const unsigned int& sz) :
 368                 _buf(buf), _sz(sz) {
 369                 }
 370
 371                 virtual void run(void) {
 372                         rgb_2_YCbCr(_buf, _sz);
 373                 }
 374         };
 375
 376         class Y_job : public mt::ThreadPool::Job {
 377                 unsigned char           *_buf;
 378                 const unsigned int      _sz;
 379         public:
 380                 Y_job(unsigned char *buf, const unsigned int& sz) :
 381                 _buf(buf), _sz(sz) {
 382                 }
 383
 384                 virtual void run(void) {
 385                         rgb_2_Y(_buf, _sz);
 386                 }
 387         };
 388
 389         static void rgb_2_hsi_tp(VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 390                 const unsigned int                      sz = v_ok.size();
 391                 std::vector<shared_ptr<hsi_job> >       v_jobs;
 392                 v_jobs.push_back(new hsi_job(&ref[0], ref.size()));
 393                 __stats_tp.add(v_jobs.rbegin()->get());
 394                 for(unsigned int i =0; i < sz; ++i) {
 395                         if (v_ok[i]) {
 396                                 v_jobs.push_back(new hsi_job(&(streams[i][0]), streams[i].size()));
 397                                 __stats_tp.add(v_jobs.rbegin()->get());
 398                         }
 399                 }
 400                 //wait for all
 401                 for(std::vector<shared_ptr<hsi_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 402                         (*it)->wait();
 403                         (*it) = 0;
 404                 }
 405         }
 406
 407         static void rgb_2_YCbCr_tp(VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 408                 const unsigned int                      sz = v_ok.size();
 409                 std::vector<shared_ptr<YCbCr_job> >     v_jobs;
 410                 v_jobs.push_back(new YCbCr_job(&ref[0], ref.size()));
 411                 __stats_tp.add(v_jobs.rbegin()->get());
 412                 for(unsigned int i =0; i < sz; ++i) {
 413                         if (v_ok[i]) {
 414                                 v_jobs.push_back(new YCbCr_job(&(streams[i][0]), streams[i].size()));
 415                                 __stats_tp.add(v_jobs.rbegin()->get());
 416                         }
 417                 }
 418                 //wait for all
 419                 for(std::vector<shared_ptr<YCbCr_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 420                         (*it)->wait();
 421                         (*it) = 0;
 422                 }
 423         }
 424
 425         static void rgb_2_Y_tp(VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 426                 const unsigned int                      sz = v_ok.size();
 427                 std::vector<shared_ptr<Y_job> >         v_jobs;
 428                 v_jobs.push_back(new Y_job(&ref[0], ref.size()));
 429                 __stats_tp.add(v_jobs.rbegin()->get());
 430                 for(unsigned int i =0; i < sz; ++i) {
 431                         if (v_ok[i]) {
 432                                 v_jobs.push_back(new Y_job(&(streams[i][0]), streams[i].size()));
 433                                 __stats_tp.add(v_jobs.rbegin()->get());
 434                         }
 435                 }
 436                 //wait for all
 437                 for(std::vector<shared_ptr<Y_job> >::iterator it = v_jobs.begin(); it != v_jobs.end(); ++it) {
 438                         (*it)->wait();
 439                         (*it) = 0;
 440                 }
 441         }
 442
 443         class psnr : public s_base {
 444                 std::string     _colorspace;
 445         protected:
 446                 void print(const int& ref_frame, const std::vector<double>& v_res) {
 447                         _ostr << ref_frame << ',';
 448                         for(int i = 0; i < _n_streams; ++i)
 449                                 _ostr << v_res[i] << ',';
 450                         _ostr << std::endl;
 451                 }
 452
 453                 void process_colorspace(VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 454                         if (_colorspace == "hsi") {
 455                                 rgb_2_hsi_tp(ref, v_ok, streams);
 456                         } else if (_colorspace == "ycbcr") {
 457                                 rgb_2_YCbCr_tp(ref, v_ok, streams);
 458                         } else if (_colorspace == "y") {
 459                                 rgb_2_Y_tp(ref, v_ok, streams);
 460                         }
 461                 }
 462         public:
 463                 psnr(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 464                 s_base(n_streams, i_width, i_height, ostr), _colorspace("rgb") {
 465                 }
 466
 467                 virtual void set_parameter(const char* p_name, const char *p_value) {
 468                         const std::string s_name(p_name);
 469                         if (s_name == "colorspace") {
 470                                 _colorspace = p_value;
 471                                 if (_colorspace != "rgb" && _colorspace != "hsi"
 472                                 && _colorspace != "ycbcr" && _colorspace != "y")
 473                                         throw std::runtime_error("Invalid colorspace passed to analyzer");
 474                         }
 475                 }
 476
 477                 virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 478                         if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 479                         // process colorspace
 480                         process_colorspace(ref, v_ok, streams);
 481                         //
 482                         std::vector<double>     v_res(_n_streams);
 483                         get_psnr_tp(ref, v_ok, streams, v_res);
 484                         //
 485                         print(ref_frame, v_res);
 486                 }
 487         };
 488
 489         class avg_psnr : public psnr {
 490                 int                     _fpa,
 491                                         _accum_f,
 492                                         _last_frame;
 493                 std::vector<double>     _accum_v;
 494         protected:
 495                 void print(const int& ref_frame) {
 496                         // check if we have to print infor or not
 497                         if (_accum_f < _fpa) return;
 498                         // else print data
 499                         for(int i=0; i < _n_streams; ++i)
 500                                 _accum_v[i] /= _accum_f;
 501                         psnr::print(ref_frame, _accum_v);
 502                         // clear the vector and set _accum_f to 0
 503                         _accum_v.clear();
 504                         _accum_v.resize(_n_streams);
 505                         _accum_f = 0;
 506                 }
 507         public:
 508                 avg_psnr(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 509                 psnr(n_streams, i_width, i_height, ostr), _fpa(1), _accum_f(0), _last_frame(-1), _accum_v(n_streams) {
 510                 }
 511
 512                 virtual void set_parameter(const char* p_name, const char *p_value) {
 513                         psnr::set_parameter(p_name, p_value);
 514                         //
 515                         const std::string s_name(p_name);
 516                         if (s_name == "fpa") {
 517                                 const int fpa = atoi(p_value);
 518                                 if (fpa > 0) _fpa = fpa;
 519                         }
 520                 }
 521
 522                 virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 523                         if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 524                         // set last frame
 525                         _last_frame = ref_frame;
 526                         // process colorspace
 527                         process_colorspace(ref, v_ok, streams);
 528                         // compute the psnr
 529                         std::vector<double>     v_res(_n_streams);
 530                         get_psnr_tp(ref, v_ok, streams, v_res);
 531                         // accumulate for each
 532                         for(int i = 0; i < _n_streams; ++i) {
 533                                 if (v_ok[i]) {
 534                                          _accum_v[i] += v_res[i];
 535                                 } else _accum_v[i] += 0.0;
 536                         }
 537                         ++_accum_f;
 538                         // try to print data
 539                         print(ref_frame);
 540                 }
 541
 542                 virtual ~avg_psnr() {
 543                         // on exit check if we have some data
 544                         if (_accum_f > 0) {
 545                                 _ostr << _last_frame << ',';
 546                                 for(int i = 0; i < _n_streams; ++i) {
 547                                         _ostr << _accum_v[i]/_accum_f << ',';
 548                                 }
 549                                 _ostr << std::endl;
 550                         }
 551                 }
 552         };
 553
 554         class ssim : public s_base {
 555         protected:
 556                 int     _blocksize;
 557
 558                 void print(const int& ref_frame, const std::vector<double>& v_res) {
 559                         _ostr << ref_frame << ',';
 560                         for(int i = 0; i < _n_streams; ++i)
 561                                 _ostr << v_res[i] << ',';
 562                         _ostr << std::endl;
 563                 }
 564         public:
 565                 ssim(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 566                 s_base(n_streams, i_width, i_height, ostr), _blocksize(8) {
 567                 }
 568
 569                 virtual void set_parameter(const char* p_name, const char *p_value) {
 570                         const std::string s_name(p_name);
 571                         if (s_name == "blocksize") {
 572                                 const int blocksize = atoi(p_value);
 573                                 if (blocksize > 0) _blocksize = blocksize;
 574                         }
 575                 }
 576
 577                 virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 578                         if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 579                         // convert to Y colorspace
 580                         rgb_2_Y_tp(ref, v_ok, streams);
 581                         //
 582                         std::vector<double>     v_res(_n_streams);
 583                         get_ssim_tp(ref, v_ok, streams, v_res, _i_width, _i_height, _blocksize);
 584                         //
 585                         print(ref_frame, v_res);
 586                 }
 587         };
 588
 589         class avg_ssim : public ssim {
 590                 int                     _fpa,
 591                                         _accum_f,
 592                                         _last_frame;
 593                 std::vector<double>     _accum_v;
 594         protected:
 595                 void print(const int& ref_frame) {
 596                         // check if we have to print infor or not
 597                         if (_accum_f < _fpa) return;
 598                         // else print data
 599                         for(int i=0; i < _n_streams; ++i)
 600                                 _accum_v[i] /= _accum_f;
 601                         ssim::print(ref_frame, _accum_v);
 602                         // clear the vector and set _accum_f to 0
 603                         _accum_v.clear();
 604                         _accum_v.resize(_n_streams);
 605                         _accum_f = 0;
 606                 }
 607         public:
 608                 avg_ssim(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 609                 ssim(n_streams, i_width, i_height, ostr), _fpa(1), _accum_f(0), _last_frame(-1), _accum_v(n_streams) {
 610                 }
 611
 612                 virtual void set_parameter(const char* p_name, const char *p_value) {
 613                         ssim::set_parameter(p_name, p_value);
 614                         //
 615                         const std::string s_name(p_name);
 616                         if (s_name == "fpa") {
 617                                 const int fpa = atoi(p_value);
 618                                 if (fpa > 0) _fpa = fpa;
 619                         }
 620                 }
 621
 622                 virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 623                         if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 624                         // set last frame
 625                         _last_frame = ref_frame;
 626                         // convert to Y colorspace
 627                         rgb_2_Y_tp(ref, v_ok, streams);
 628                         //
 629                         std::vector<double>     v_res(_n_streams);
 630                         get_ssim_tp(ref, v_ok, streams, v_res, _i_width, _i_height, _blocksize);
 631                         // accumulate for each
 632                         for(int i = 0; i < _n_streams; ++i) {
 633                                 if (v_ok[i]) {
 634                                          _accum_v[i] += v_res[i];
 635                                 } else _accum_v[i] += 0.0;
 636                         }
 637                         ++_accum_f;
 638                         // try to print data
 639                         print(ref_frame);
 640                 }
 641
 642                 virtual ~avg_ssim() {
 643                         // on exit check if we have some data
 644                         if (_accum_f > 0) {
 645                                 _ostr << _last_frame << ',';
 646                                 for(int i = 0; i < _n_streams; ++i) {
 647                                         _ostr << _accum_v[i]/_accum_f << ',';
 648                                 }
 649                                 _ostr << std::endl;
 650                         }
 651                 }
 652         };
 653
 654         class siti : public s_base {
 655                 protected:
 656                         VUCHAR* _prev;
 657                         void print(const int& ref_frame, const std::vector<std::pair<double,double> >& v_res) {
 658                                 _ostr << ref_frame << ',';
 659                                 for(int i = 0; i < _n_streams; ++i)
 660                                         _ostr << v_res[i].first << ',' << v_res[i].second << ",";
 661                                 _ostr << std::endl;
 662                         }
 663                 public:
 664                         siti(const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) :
 665                         s_base(n_streams, i_width, i_height, ostr){
 666                                 _prev=NULL;
 667                         }
 668
 669                         virtual void process(const int& ref_frame, VUCHAR& ref, const std::vector<bool>& v_ok, std::vector<VUCHAR>& streams) {
 670                                 if (v_ok.size() != streams.size() || v_ok.size() != _n_streams) throw std::runtime_error("Invalid data size passed to analyzer");
 671                                 // convert to Y colorspace
 672                                 rgb_2_Y_tp(ref, v_ok, streams);
 673                                 //
 674                                 if(_prev==NULL){
 675                                         _prev=new std::vector<unsigned char> (streams[0]);
 676                                 }
 677                                 std::vector<std::pair<double,double> >  v_res(_n_streams);
 678                                 get_siti_tp(_prev, v_ok, streams, v_res, _i_width, _i_height);
 679                                 delete(_prev);
 680                                 _prev=new std::vector<unsigned char> (streams[0]);
 681                                 //
 682                                 print(ref_frame, v_res);
 683                         }
 684                 };
 685 }
 686
 687 stats::s_base* stats::get_analyzer(const char* id, const int& n_streams, const int& i_width, const int& i_height, std::ostream& ostr) {
 688         const std::string       s_id(id);
 689         if (s_id == "psnr") return new psnr(n_streams, i_width, i_height, ostr);
 690         else if (s_id == "avg_psnr") return new avg_psnr(n_streams, i_width, i_height, ostr);
 691         else if (s_id == "ssim") return new ssim(n_streams, i_width, i_height, ostr);
 692         else if (s_id == "avg_ssim") return new avg_ssim(n_streams, i_width, i_height, ostr);
 693         else if (s_id == "siti") return new siti(n_streams, i_width, i_height, ostr);
 694         throw std::runtime_error("Invalid analyzer id");
 695 }