StSD_lc(distorted_left, reference_left, distorted_right,...

function [quality] = StSD_lc(distorted_left, reference_left, distorted_right,...
    reference_right, num_frames, rows, cols)
%========================================================================
%
% Calculation of StSD image quality measure
%
% StSD is Stereoscopic Structural Distortion metric taking into account 
% asymmetric structural distortions, asymmetric blur and information    
% information content level in to accout to provide a measure of percieved
% distortion of stereoscopic video
%
% Copyright(c) 2012 Varuna De Silva & Hemantha Kodikararachchi 
% All Rights Reserved
%
% E-mail: varunax{AT}gmail.com
%
%----------------------------------------------------------------------
%
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is hereby
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm for calculating the StSD
% metric between two videos. Please cite the following paper:
%
% [1] Varuna De Silva, Hemantha Kodikararachchi and Ahmet Kondoz, 
%     "Towards an Impairment Metric for Stereoscopic Video: A 
%     Full-Reference Video Quality Metric to Assess Compressed Stereoscopic
%     Video", IEEE Transactions on Image Processing, June 2013
%
% This paper is accepted for publication in  IEEE transactions 
% on Image Processing. 
%
% Kindly report any suggestions or corrections to varunax{AT}gmail[dot]com
%
%----------------------------------------------------------------------
%
% Input : (1) distorted_files: the distorted video sequence in YUV format
%         (2) reference_files: the reference video sequence in YUV format
%         The inputs should be in YUV 420 format 
%         (3) num_frames: the number of frames to process
%         (4) rows: the number of rows in yuv file
%         (5) cols: the number of cols in yuv file
% However, it is assumed that the videos are watched on a full HD screen.
%
% Output: (1) quality: the The StSD measure.
%
%
%   [quality] = StSD_lc(distorted_left.yuv, reference_left.yuv, 
%   distorted_right.yuv,reference_right.yuv,100, 1080, 960)
%
% See the results:
%  
%
%========================================================================

       
%Process the video frames one after other

%scale the number of frames to a multiple of 10 for memory management
%issues
step = 10; %Reduce this number to n if you need StSD every n frames
%Increase the step size when running for larger number of frames. 
nf = mod(num_frames,step);

siml = zeros((num_frames-nf)/step,step);
simr = zeros((num_frames-nf)/step,step);
blr_l = zeros((num_frames-nf)/10,10);
blr_r = zeros((num_frames-nf)/10,10);

h = fspecial('sobel'); 
hb = waitbar(0,'Please wait...');
progress = 0;
for j = 1:(num_frames-nf)/step;
    [Y1l, ~, ~] = yuv_import(distorted_left,[cols rows],step,(j-1)*step);
    [Y2l, ~, ~] = yuv_import(reference_left,[cols rows],step,(j-1)*step);
    
    [Y1r, ~, ~] = yuv_import(distorted_right,[cols rows],step,(j-1)*step);
    [Y2r, ~, ~] = yuv_import(reference_right,[cols rows],step,(j-1)*step);
    
    for s = 1:1:step,
        
        B = Y1l{s};
        A = Y2l{s};
        D = Y1r{s};
        C = Y2r{s};
        
        %Structural distortion
        [siml(j,s)] = structure_frame(A,B);
        [simr(j,s)] = structure_frame(C,D);
        
        %Blur
        Edge_left_ref1 = imfilter(uint8(A),h);
        Edge_left_prc1 = imfilter(uint8(B),h);
        Edge_right_ref1 = imfilter(uint8(C),h);
        Edge_right_prc1 = imfilter(uint8(D),h);
        Edge_left_ref2 = imfilter(uint8(A),h');
        Edge_left_prc2 = imfilter(uint8(B),h');
        Edge_right_ref2 = imfilter(uint8(C),h');
        Edge_right_prc2 = imfilter(uint8(D),h');
        Edge_left_ref = abs(Edge_left_ref1)+ abs(Edge_left_ref2);
        Edge_left_prc = abs(Edge_left_prc1)+ abs(Edge_left_prc2);
        Edge_right_ref = abs(Edge_right_ref1)+ abs(Edge_right_ref2);
        Edge_right_prc = abs(Edge_right_prc1)+ abs(Edge_right_prc2);
        
        Th1 = mean2([Edge_left_ref, Edge_right_ref]);
        Th2 = std2([Edge_left_ref, Edge_right_ref])/2;
        
        diff_edge_l = abs(Edge_left_ref(Edge_left_ref>Th1)-Edge_left_prc(Edge_left_ref>Th1));
        diff_edge_r = abs(Edge_right_ref(Edge_right_ref>Th1)-Edge_right_prc(Edge_right_ref>Th1));
        
        blr_l(j,s) = sum(diff_edge_l((diff_edge_l-Th2)>0))/(rows*cols);
        blr_r(j,s) = sum(diff_edge_r((diff_edge_r-Th2)>0))/(rows*cols);
        
        %SI & TI calculation
        %should be there by now!
        progress = (j-1)*step+s;
        waitbar(progress / (num_frames-nf))
    end
end

close(hb);

%Averaging of results
mean_sl = sum(siml( ~(isinf(siml)|isnan(siml))))/(sum(sum( ~(isinf(siml)|isnan(siml))))); %mean2(siml); 
mean_sr = sum(simr( ~(isinf(simr)|isnan(simr))))/(sum(sum( ~(isinf(simr)|isnan(simr)))));% mean2(simr);

mean_blrl = sum(blr_l( ~(isinf(blr_l)|isnan(blr_l))))/(sum(sum( ~(isinf(blr_l)|isnan(blr_l))))); %mean2(blr_l);
mean_blrr = sum(blr_r( ~(isinf(blr_r)|isnan(blr_r))))/(sum(sum( ~(isinf(blr_r)|isnan(blr_r))))); %mean2(blr_r);

%Aggregation of results
%Structrue: add the structural distortions
s2 = real(mean_sl + mean_sr);
k2=0.7343; r2=15.7784; t2=0.1397;         %------------------From Training
D = k2./(1+exp(-r2*(s2-t2))); %Eq. (14)

%Blur: to aid reduce prediction error
if (D<=0.8 & D>=0.45)                         %------------------From Training
    B = min(mean_blrl, mean_blrr);
else
    B = 0;
end

pe = -0.0730+0.0085*B; %-----------------From Training

quality = D+pe;

end

function [d_cod_sim] = structure_frame(ref,proc)
%This function intakes the 1/2 width and 1/2height
ref = imresize(ref,[540 960]);
proc = imresize(proc,[540 960]);

%compute similariy and difference vectors
%downsample images
dn_ref = imresize(ref,[273 481]); %21 rows, 37
dn_prc= imresize(proc,[273 481]);

%process every 13by13 block
S = zeros(21,37);

for i = 1:21
    for j = 1:37
        block_ref = dn_ref((i-1)*13+1:(i-1)*13+13,(j-1)*13+1:(j-1)*13+13);
        block_prc = dn_prc((i-1)*13+1:(i-1)*13+13,(j-1)*13+1:(j-1)*13+13);
        A = cov(block_ref(:),block_prc(:));
        S(i,j) = (A(1,2)+25)/(mean2((block_ref-mean2(block_ref)).^2)+25);
        
    end
end

%compute feature values
s_m = trimmean(S(:),20);

s_delta = s_m - mean(S(S<quantile(S(:),0.20)));

d_sim = 1 - s_m +1.5*s_delta;

d_cod_sim = stransformf(d_sim,0.07,0.1,2);

end

function s = stransformf(x, p_x, p_y, q)
d = (1-p_y);
c = 2*q/d;
b = q*p_x/p_y;
a = p_y/(p_x^b);

if x<=p_x
    s = a*x^b;
else
    s = 2*d*(1/(1+exp(-c*(x - p_x)))-1/2)  + p_y;
end

end