static double calc_ssim()

in libvmaf/src/feature/integer_ssim.c [88:197]

• 110 lines of code
• 11 McCabe index (conditional complexity)

static double calc_ssim(const unsigned char *_src,int _systride,
 const unsigned char *_dst,int _dystride,double _par,int depth,int _w,int _h){
  ssim_moments  *line_buf;
  ssim_moments **lines;
  double         ssim;
  double         ssimw;
  unsigned      *hkernel;
  int            hkernel_sz;
  int            hkernel_offs;
  unsigned      *vkernel;
  int            vkernel_sz;
  int            vkernel_offs;
  int            log_line_sz;
  int            line_sz;
  int            line_mask;
  int            x;
  int            y;
  int            samplemax;
  samplemax = (1 << depth) - 1;
  vkernel_sz=gaussian_filter_init(&vkernel,1.5,5);
  vkernel_offs=vkernel_sz>>1;
  for(line_sz=1,log_line_sz=0;line_sz<vkernel_sz;line_sz<<=1,log_line_sz++);
  line_mask=line_sz-1;
  lines=(ssim_moments **)malloc(line_sz*sizeof(*lines));
  lines[0]=line_buf=(ssim_moments *)malloc(line_sz*_w*sizeof(*line_buf));
  for(y=1;y<line_sz;y++)lines[y]=lines[y-1]+_w;
  hkernel_sz=gaussian_filter_init(&hkernel,1.5,5);
  hkernel_offs=hkernel_sz>>1;
  ssim=0;
  ssimw=0;
  for(y=0;y<_h+vkernel_offs;y++){
    ssim_moments *buf;
    int           k;
    int           k_min;
    int           k_max;
    if(y<_h){
      buf=lines[y&line_mask];
      for(x=0;x<_w;x++){
        ssim_moments m;
        memset(&m,0,sizeof(m));
        k_min=hkernel_offs-x<=0?0:hkernel_offs-x;
        k_max=x+hkernel_offs-_w+1<=0?
         hkernel_sz:hkernel_sz-(x+hkernel_offs-_w+1);
        for(k=k_min;k<k_max;k++){
          signed s;
          signed d;
          signed window;
          if (depth > 8) {
            s = _src[(x-hkernel_offs+k)*2] +
             (_src[(x-hkernel_offs+k)*2 + 1] << 8);
            d = _dst[(x-hkernel_offs+k)*2] +
             (_dst[(x-hkernel_offs+k)*2 + 1] << 8);
          } else {
            s=_src[(x-hkernel_offs+k)];
            d=_dst[(x-hkernel_offs+k)];
          }
          window=hkernel[k];
          m.mux+=window*s;
          m.muy+=window*d;
          m.x2+=window*s*s;
          m.xy+=window*s*d;
          m.y2+=window*d*d;
          m.w+=window;
        }
        *(buf+x)=*&m;
      }
      _src+=_systride;
      _dst+=_dystride;
    }
    if(y>=vkernel_offs){
      k_min=vkernel_sz-y-1<=0?0:vkernel_sz-y-1;
      k_max=y+1-_h<=0?vkernel_sz:vkernel_sz-(y+1-_h);
      for(x=0;x<_w;x++){
        ssim_moments m;
        double       c1;
        double       c2;
        double       mx2;
        double       mxy;
        double       my2;
        double       w;
        memset(&m,0,sizeof(m));
        for(k=k_min;k<k_max;k++){
          signed window;
          buf = lines[(y + 1 - vkernel_sz + k) & line_mask] + x;
          window=vkernel[k];
          m.mux+=window*buf->mux;
          m.muy+=window*buf->muy;
          m.x2+=window*buf->x2;
          m.xy+=window*buf->xy;
          m.y2+=window*buf->y2;
          m.w+=window*buf->w;
        }
        w=m.w;
        c1=samplemax*samplemax*SSIM_K1*w*w;
        c2=samplemax*samplemax*SSIM_K2*w*w;
        mx2=m.mux*(double)m.mux;
        mxy=m.mux*(double)m.muy;
        my2=m.muy*(double)m.muy;
        ssim+=m.w*(2*mxy+c1)*(c2+2*(m.xy*w-mxy))/
         ((mx2+my2+c1)*(m.x2*w-mx2+m.y2*w-my2+c2));
        ssimw+=m.w;
      }
    }
  }
  free(line_buf);
  free(lines);
  free(vkernel);
  free(hkernel);
  return ssim/ssimw;
}