-- M module extracted from ITU-T J.249 (01/2010)


function [s1,s2,s3] = block_statistic(y,vsize,hsize,varargin)
% BLOCK_STATISTIC
%  Extract feature from each spatial-temporal (S-T) region, producing one number 
%  for each S-T block.  Takes a block of perceptually filtered images 
%  and produces features.  
% SYNTAX
%  [s1] = block_statistic(y,vsize,hsize,'Stat1');
%  [s1,s2] = block_statistic(y,vsize,hsize,'Stat1','Stat2');
%  [s1,s2,s3] = block_statistic(y,vsize,hsize,'Stat1','Stat2','Stat3');
% DEFINITION
%  [s1] = block_statistic(y,'Stat1'); divides time-slice of images 'y' into
%  abutting Spatial-Temporal (S-T) regions that contain 'vsize' pixels
%  vertically and 'hsize' pixels horizontally.  Then, statstic 'Stat1' is
%  computed over each S-T region, and the results are returned in 's1'.
%  When called with the names of two statistics, two statistics are 
%  computed and returned, and so forth.  Available statistics are:
%  'mean' ==> compute the mean over each S-T region.
%  'std' ==> compute the standard deviation over each S-T region.
%  'rms' ==> compute the RMS over each S-T region.
%  'fraction' ==> compute fraction of pixels that are greater than or
%                 equal to 1.0
% REMARKS
%  Functionality tested pretty well.
 
want_mean = 0;
want_std = 0;
want_rms = 0;
want_fraction = 0;
 
for cnt = 1:nargin-3,
    if strcmp(lower(varargin(cnt)),'mean'),
        want_mean = 1;
        want_mean_at = cnt;
    elseif strcmp(lower(varargin(cnt)),'std'),
        want_std = 1;
        want_std_at = cnt;
    elseif strcmp(lower(varargin(cnt)),'rms'),
        want_rms = 1;
        want_rms_at = cnt;
    elseif strcmp(lower(varargin(cnt)),'fraction'),
        want_fraction = 1;
        want_fraction_at = cnt;
    end
end
 
if want_mean + want_std + want_rms + want_fraction ~= nargin - 3,
    error('block_statistic ''Stat'' not recognized or repeated.');
end
 
% check block size request.
[row,col,time] = size(y);
if mod(row,vsize) ~= 0,
    error('vertical size of block must evenly divide the SROI.');
end
 
if mod(col,hsize) ~= 0,
    error('horizontal size of block must evenly divide the SROI.');
end
 
if want_mean | want_std,
    temp = sum(y,3);  % sum over time
    temp = sum( reshape(temp,vsize,row*col/vsize)); % sum block vertically
    temp = reshape(temp,row/vsize,col)';
    temp = sum(reshape(temp,hsize,row*col/(vsize*hsize))); % sum block horizontally
    y_sum = reshape(temp,col/hsize,row/vsize)' ./ (hsize*vsize*time); % reshape
end
 
if want_std | want_rms,
    temp = sum(y.^2,3);
    temp = sum(reshape(temp,vsize,row*col/vsize));
    temp = reshape(temp,row/vsize,col)';
    temp = sum(reshape(temp,hsize,row*col/(vsize*hsize)));
    y_square = reshape(temp,col/hsize,row/vsize)' ./ (hsize*vsize*time);
end
 
if want_fraction,
    y(find(y >= 1.0)) = 1.0;
    y(find(y < 1.0)) = 0.0;
    
    temp = sum(y,3);  % sum over time
    temp = sum( reshape(temp,vsize,row*col/vsize)); % sum block vertically
    temp = reshape(temp,row/vsize,col)';
    temp = sum(reshape(temp,hsize,row*col/(vsize*hsize))); % sum block horizontally
    y_fraction = reshape(temp,col/hsize,row/vsize)' ./ (hsize*vsize*time); % reshape
end
 
 
if want_mean,
    switch want_mean_at,
        case 1,
            s1 = y_sum;
        case 2,
            s2 = y_sum;
        case 3,
            s3 = y_sum;
        otherwise
            error('Code defect')
    end
end
 
if want_std,
    y_std = sqrt( max(y_square - y_sum .^ 2,0));
    switch want_std_at,
        case 1,
            s1 = y_std;
        case 2,
            s2 = y_std;
        case 3,
            s3 = y_std;
        otherwise
            error('Code defect')
    end
end
 
if want_rms,
    y_rms = sqrt( y_square );
    switch want_rms_at,
        case 1,
            s1 = y_rms;
        case 2,
            s2 = y_rms;
        case 3,
            s3 = y_rms;
        otherwise
            error('Code defect')
    end
end
 
 
if want_fraction,
    switch want_fraction_at,
        case 1,
            s1 = y_fraction;
        case 2,
            s2 = y_fraction;
        case 3,
            s3 = y_fraction;
        otherwise
            error('Code defect')
    end
end