function par = vfd_clippar_loop_psnr(test_structs, clip_structs, feature_base_dir)
% VFD_CLIPPAR_LOOP_PSNR
%  Computes the psnr_vfd parameter of the VQM_VFD model. This function
%  takes variable (1) 'test_structs' (of the same format as GTests), which
%  describes each video test and the location of the associated video, and
%  variable (2) 'clip_structs' (of the same format as GClips), which
%  specifies the set of video clips, and (3) the feature directory
%  'feature_base_dir', which specifies the directory location of the
%  features. The returned par will be an ITS parameter structure that
%  contains the parameters for all the processed clips in clip_structs. The
%  sort_clips_by function is used to sort the clips before parameter
%  calculation, and the original clips are skipped if they are encountered.
% SYNTAX
%  par = vfd_clippar_loop_psnr(test_structs, clip_structs, feature_base_dir)
% DESCRIPTION
%  This script uses the features generated by vfd_feature_loop_mse to
%  compute the psnr_vfd described in NTIA TM-11-475. In particular, the
%  feature feature_Y_FR_vfd_0.4deg_0.2s_rms must reside in the
%  feature_base_dir. This feature is converted to PSNR using a peak signal
%  level of 255. PSNR is clipped at an upper level of 48 dB.
%

peak_signal = 255;  % Defines the peak signal level for PSNR
psnr_limit = 48;  % limit on PSNR max value, otherwise PSNR could correlate poorly

%  Define the features used for the block parameters.
%  This psnr_vfd parameter uses 0.4deg_0.2s ST blocks.
f_psnr = 'feature_Y_FR_vfd_0.4deg_0.2s_rms';

% Define the par_name for the returned parameter structure.
par.par_name = {'psnr_vfd'};

% Loop through all clips, sorted alphabetically
offsets = sort_clips_by('none', clip_structs, test_structs);
clip_structs = clip_structs(offsets);
ccnt = 1;  % processed clip counter for returned par array
for loop = 1:length(offsets)
    
    % Skip if an original clip was included.
    if strcmpi(clip_structs(loop).hrc{1},'original'),
        continue;
    end
    
    % Load rmse features
    data = zeros(1,0);
    datao = zeros(1,0);
    name = sprintf('%s/%s/%s_%s_%s.mat', feature_base_dir, f_psnr, ...
        clip_structs(loop).test{1}, ...
        clip_structs(loop).scene{1}, ...
        clip_structs(loop).hrc{1});
    if ~exist(name, 'file'),
        fprintf('Skipping clip %s:%s(%s), RMSE feature filename does not exist.\n', clip_structs(loop).test{1}, clip_structs(loop).scene{1}, clip_structs(loop).hrc{1});
        continue;
    end
    load( name );
    
    % Also skip this clip if no data exists (datao is a 3D array of zeros
    % in this case so no need to check if that exists).
    [proc_r,proc_c,proc_t] = size(data);
    if proc_r*proc_c*proc_t == 0,
        fprintf('Skipping clip %s:%s(%s), no processed RMSE features.\n', clip_structs(loop).test{1}, clip_structs(loop).scene{1}, clip_structs(loop).hrc{1});
        continue;
    end
    this_par = data;
    clear data datao;
    
    % Convert RMSE to MSE of each block by squaring.  Then compute mean
    % over all space and time (i.e., the MSE of the video sequence).  Then
    % convert to PSNR.
    this_par = this_par.^2;  % convert rmse to mse
    
    % Calculate mean over ST
    this_par = st_collapse('mean', this_par, '3D');
    
    % Convert to PSNR
    if (this_par ~= 0)  % Might be zero MSE in rare cases
        this_par = 10*log10(peak_signal*peak_signal/this_par);
        this_par = min(this_par, psnr_limit);
    else
        this_par = psnr_limit;
    end
    
    % Fill the clip_structs information into the returned parameter.
    par.clip_name{ccnt} = sprintf('%s_%s_%s', clip_structs(loop).test{1}, ...
        clip_structs(loop).scene{1}, clip_structs(loop).hrc{1});
    par.inlsa_mos(ccnt) = clip_structs(loop).inlsa_mos;
    par.mos(ccnt) = clip_structs(loop).mos;
    par.data(1, ccnt) = this_par;
    
    % Update clip counter.
    ccnt = ccnt + 1;
    
end


