text:
function r8vec_indexed_heap_d_extract_test ( )

%*****************************************************************************80
%
%% R8VEC_INDEXED_HEAP_D_EXTRACT_TEST tests R8VEC_INDEXED_HEAP_D_EXTRACT.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
  m = 20;
  n_max = 20;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8VEC_INDEXED_HEAP_D_EXTRACT_TEST\n' );
  fprintf ( 1, '  For an indexed R8VEC,\n' );
  fprintf ( 1, '  R8VEC_INDEXED_HEAP_D_EXTRACT extracts the maximum value;\n' );
%
%  Set the data array.  To keep things easy, we will use the indicator vector.
%
  a = r8vec_indicator1 ( m );
%
%  The index array will initially be a random subset of the numbers 1 to M,
%  in random order.
%
  n = 5;
  indx(1:11,1) = [ 9, 2, 8, 14, 5, 7, 15, 1, 19, 20, 3 ]';

  r8vec_print ( m, a, '  The data vector:' );
  i4vec_print ( n, indx, '  The index vector:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A(INDX):\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %14.6f\n', i, a(indx(i)) );
  end
%
%  Create the descending heap.
%
  indx = r8vec_indexed_heap_d ( n, a, indx );

  i4vec_print ( n, indx, '  The index vector after heaping:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A(INDX) after heaping:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %14.6f\n', i, a(indx(i)) );
  end
%
%  Extract the first 5 largest elements.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Now extract the maximum several times.\n' );
  fprintf ( 1, '\n' );

  for i = 1 : 5
    [ indx_extract, n, indx ] = r8vec_indexed_heap_d_extract ( n, a, indx );
    fprintf ( 1, '  Extracting maximum element A(%d) = %f\n', ...
      indx_extract, a(indx_extract) );
  end

  r8vec_print ( m, a, '  The data vector after extractions:' );
  i4vec_print ( n, indx, '  The index vector after extractions:' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  A(INDX) after extractions:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : n
    fprintf ( 1, '  %4d  %14.6f\n', i, a(indx(i)) );
  end

  return
end


text:
function header=make_image_http_header(html,found)
if(found)
    header.HTTP='200 OK';
else
    header.HTTP='404 Not Found';
end
header.Date='Tue, 12 Oct 2010 09:19:05 GMT';
header.Server='Matlab Webserver';
header.LastModified='Last-Modified: Thu, 21 Jun 2007 14:56:37 GMT';
header.AcceptRanges='Accept-Ranges: bytes';
header.ContentLength=num2str(length(html));
header.ETag='"948921-15ae-c0dbf340"';
%header.KeepAlive='Keep-Alive: timeout=15, max=100';
header.ContentType='image/png';

author: jaderberg

repo: python-matlab-bridge


text:
% dc_remove remove the DC offset of the time domain signal

function [z] = dc_remove(x,a);

Num = [a -a];
Dem = [1 -a];
z = filter(Num,Dem,x);

author: singaxiong

repo: SignalGraph


text:
%  Figure 7.21      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% Script to generate fig. 7.21.
clf;
num=1;
den=[1 0 0 0 0];
rlocus(num,den);
title('Fig. 7.21 Symmetric root locus for satellite system')
grid;

author: Sable

repo: mcbench-benchmarks


text:
function check_unimodalGrad

%% test 1

cs1 = crystalSymmetry('-3m1');
cs2 = crystalSymmetry('m-3m');

center = orientation.rand(100,cs1,cs2);

odf = unimodalODF(center);

ori = orientation.rand(1000,odf.CS,odf.SS);

g1 = odf.grad(ori(:),'check','delta',0.001*degree);
g2 = odf.grad(ori(:));

if max(norm(g1-g2)./(1+norm(g1))) < 1e-1
  disp(' Unimoal gradient test passed');
else
  disp(' Unimoal gradient test failed');
end

%% test 2
odf = SantaFe;

ori = orientation.rand(1000,odf.CS,odf.SS);

g1 = odf.grad(ori(:),'check','delta',0.001*degree);
g2 = odf.grad(ori(:));

if max(norm(g1-g2)./norm(g1)) < 1e-2
  disp(' Unimoal gradient test passed');
else
  disp(' Unimoal gradient test failed');
end

end

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [success,origScanNum] = segmentationWrapper(cerrPath,fullSessionPath,...
    containerPath, algorithm, skipMaskExport)
% function success =heart(cerrPath,segResultCERRPath,fullSessionPath,deepLabContainerPath)
%
% This function serves as a wrapper for the all segmentation models
%
% INPUT:
% cerrPath - path to the original CERR file to be segmented
% fullSessionPath - path to write temporary segmentation metadata.
% deepLabContainerPath - path to the MR Prostate DeepLab V3+ container on the
% algorithm - name of the algorithm to run
% skipMaskExport (optional) - Set to false if model requires segmentation
%                             masks as input.Default: true.
%
%
%
% RKP, 5/21/2019
% RKP, 9/11/19 Updates for compatibility with training pipeline

if ~exist('skipMaskExport','var')
    skipMaskExport = true;
end

%Copy config file to session dir
configFilePath = fullfile(getCERRPath,'ModelImplementationLibrary',...
    'SegmentationModels', 'ModelConfigurations', [algorithm, '_config.json']);
copyfile(configFilePath,fullSessionPath);

%% Export data to fmt reqd by model
cmdFlag = 'singcontainer';
planCfiles = dir(fullfile(cerrPath,'*.mat'));
for p=1:length(planCfiles)

    % Load plan
    planCfiles(p).name
    fileNam = fullfile(planCfiles(p).folder,planCfiles(p).name);
    planC = loadPlanC(fileNam, tempdir);
    planC = quality_assure_planC(fileNam,planC);

    %Pre-process data and export to model input fmt
    [~,command,userOptS,~,scanNumV,planC] = ...
        prepDataForSeg(planC,fullSessionPath,algorithm,cmdFlag,...
        containerPath,[],skipMaskExport);

    %Save updated planC file
    tic
    save_planC(planC,[],'PASSED',fileNam);
    toc

end

%% Execute the container
disp('Running container....');
%Print command to stdout
disp(command);
tic
status = system(command);
toc

% Run container app to get git_hash
gitHash = 'unavailable';
[~,hashChk] = system(['singularity exec ' containerPath ' ls /scif/apps | grep get_hash'],'-echo');
if ~isempty(hashChk)
    [~,gitHash] = system(['singularity run --app get_hash ' containerPath],'-echo');
end
roiDescrpt = '';
if isfield(userOptS,'roiGenerationDescription')
    roiDescrpt = userOptS.roiGenerationDescription;
end
roiDescrpt = [roiDescrpt, '  __git_hash:',gitHash];
userOptS.roiGenerationDescription = roiDescrpt;

%% Stack output mask files
fprintf('\nRreading output masks...');
tic
outFmt = userOptS.modelInputFormat;
passedScanDim = userOptS.passedScanDim;
outC = stackDLMaskFiles(fullSessionPath,outFmt,passedScanDim);
toc

%% Import segmented mask to planC
fprintf('\nImporting to CERR...\n');
tic
identifierS = userOptS.structAssocScan.identifier;
labelPath = fullfile(fullSessionPath,'outputLabelMap');
if ~isempty(fieldnames(userOptS.structAssocScan.identifier))
    origScanNum = getScanNumFromIdentifiers(identifierS,planC);
else
    origScanNum = 1; %Assoc with first scan by default
end
outScanNum = scanNumV(origScanNum);
userOptS.scan(outScanNum) = userOptS.scan(origScanNum);
userOptS.scan(outScanNum).origScan = origScanNum;
success = joinH5CERR(cerrPath,outC{1},labelPath,outScanNum,userOptS); %Updated
toc

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
classdef FMeasureMerger < handle
    properties(Constant)
        USE_ORIGINAL_IDEA=false;
    end

    methods(Static)
        function [bestIds, fMeasure]=GetBest(matchIdPerRow, ...
                mergeIdPerRow, matchId, mergeIds, adaptiveBins, transpose)
            if isempty(adaptiveBins)
                [matchSize, fMeasures4MergeIds,mergeSizes,truePositives]...
                    =FMeasureMerger.Compute(matchIdPerRow, ...
                    mergeIdPerRow, matchId, mergeIds, transpose);
                if FMeasureMerger.USE_ORIGINAL_IDEA
                    [bestIds, fMeasure]=...
                        FMeasureMerger.FindBestUsingOriginalIdea(...
                        matchSize, mergeIds, fMeasures4MergeIds, mergeSizes, ...
                        truePositives);
                else
                    [bestIds, fMeasure]=FMeasureMerger.FindBest(...
                    matchSize, mergeIds, fMeasures4MergeIds, mergeSizes, ...
                    truePositives);
                end
            else
                [matchChoices, fMeasures4MergeIds, allMergeChoices] ...
                    =FMeasureMerger.ComputeBins(matchIdPerRow, ...
                    mergeIdPerRow, matchId, mergeIds, adaptiveBins, transpose);
                [bestIds, fMeasure]=FMeasureMerger.FindBestBins(...
                    mergeIds, fMeasures4MergeIds, matchChoices, ...
                    allMergeChoices, adaptiveBins, transpose);
            end
        end

        function [matchChoices, fMeasures4MergeIds, allMergeChoices]=...
                ComputeBins(idPerRow1, idPerRow2, ...
                matchId, mergeIds, adaptiveBins, transpose)
            if ~transpose
                matchChoices=MatBasics.LookForIds(idPerRow1, matchId);
            else
                matchChoices=MatBasics.LookForIds(idPerRow2, matchId);
            end
            nMergeIds=length(mergeIds);
            allMergeChoices=cell(1, nMergeIds);
            fMeasures4MergeIds=zeros(1,nMergeIds);
            for i=1:nMergeIds
                if ~transpose
                    mergeChoices=MatBasics.LookForIds(idPerRow2, mergeIds(i));
                    fMeasures4MergeIds(i)=adaptiveBins.fMeasure(matchChoices, mergeChoices);
                else
                    mergeChoices=MatBasics.LookForIds(idPerRow1, mergeIds(i));
                    fMeasures4MergeIds(i)=adaptiveBins.fMeasure(mergeChoices, matchChoices);
                end
                allMergeChoices{i}=mergeChoices;
                if isnan(fMeasures4MergeIds(i))
                    fMeasures4MergeIds(i)=0;
                end
            end
        end

        function [bestIds, fMeasure]=FindBestBins(mergeIds,...
                fMeasures4MergeIds, matchChoices, allMergeChoices, ...
                adaptiveBins, transpose)
            [maxF, idx]=max(fMeasures4MergeIds);
            bestIds=mergeIds(idx);
            fMeasure=-100;% force at least ONE merger
            nMergeCandidates=length(fMeasures4MergeIds);
            nextFms=zeros(1, nMergeCandidates);
            mergeChoices=allMergeChoices{idx};
            mergeIds(idx)=nan;
            while ~all(isnan(mergeIds))
                for i=1:nMergeCandidates
                    if isnan(mergeIds(i))
                        nextFms(i)=-1;
                    else
                        if ~transpose
                            nextFms(i)=adaptiveBins.fMeasure(...
                                matchChoices, mergeChoices|allMergeChoices{i});
                        else
                            nextFms(i)=adaptiveBins.fMeasure(...
                                mergeChoices|allMergeChoices{i}, matchChoices);
                        end
                    end
                end
                [f, idx]=max(nextFms);
                if f>fMeasure
                    fMeasure=f;
                    mergeChoices=mergeChoices|allMergeChoices{idx};
                    bestIds(end+1)=mergeIds(idx);
                else
                    break; % nothing changed!
                end
                mergeIds(idx)=nan;
            end
            if QfHiDM.F_MEASURE_MERGE_FAST==-1
                FMeasureMerger.NoteIfNoMergerImprovement(fMeasure, maxF, ...
                    bestIds(1));
            end
            bestIds=sort(bestIds);
        end

        function [matchSize, fMeasures4MergeIds, mergeSizes, truePositives]=...
                Compute(idPerRow1, idPerRow2, matchId, mergeIds, transpose)
            if ~transpose
                matchChoices=MatBasics.LookForIds(idPerRow1, matchId);
            else
                matchChoices=MatBasics.LookForIds(idPerRow2, matchId);
            end
            matchSize=sum(matchChoices);
            nMergeIds=length(mergeIds);
            fMeasures4MergeIds=zeros(1,nMergeIds);
            truePositives=zeros(1,nMergeIds);
            mergeSizes=zeros(1,nMergeIds);
            for i=1:nMergeIds
                if ~transpose
                    mergeChoices=MatBasics.LookForIds(idPerRow2, mergeIds(i));
                else
                    mergeChoices=MatBasics.LookForIds(idPerRow1, mergeIds(i));
                end
                mergeSizes(i)=sum(mergeChoices);
                truePositives(i)=sum(matchChoices&mergeChoices);
                p=truePositives(i)/mergeSizes(i);
                r=truePositives(i)/matchSize;
                fMeasures4MergeIds(i)=(2*p*r)/(p+r);
                if isnan(fMeasures4MergeIds(i))
                    fMeasures4MergeIds(i)=0;
                end
            end
        end

%This method was necessary over the
%FindBestUsingOriginalIdea when it was found that
%the next highest F-measure does not ALWAYS
%lead to the best F-measure after merging.
%This was found with the Neutrophil scenario
%when trying to find first best merger of 2 after
%starting with a merge candidate with max F-measure
%that is based on a precision of 212/584 and matchSize of 300
%and the next 2 candidates in decreasing F-measure
%order is precision 8/105 and 2/54.
%The best F-measure when combining the 2nd candidate
%is 2/54 and not 8/105

        function [bestIds, fMeasure]=FindBest(matchSize, mergeIds, ...
            fMeasures4MergeIds, mergeSizes, truePositives)
            [maxF, idx]=max(fMeasures4MergeIds);
            bestIds=mergeIds(idx);
            fMeasure=-100;% force at least ONE merger
            truePositive=truePositives(idx);
            mergeSize=mergeSizes(idx);
            truePositives(idx)=nan;
            while ~all(isnan(truePositives))
                tp=truePositive+truePositives;
                p=tp./(mergeSizes+mergeSize);
                r=tp./matchSize;
                [f, idx]=max((2*p.*r)./(p+r));
                if f>fMeasure
                    fMeasure=f;
                    truePositive=truePositive+truePositives(idx);
                    mergeSize=mergeSize+mergeSizes(idx);
                    bestIds(end+1)=mergeIds(idx);
                else
                    break; % nothing changed!
                end
                truePositives(idx)=nan;
            end
            if QfHiDM.F_MEASURE_MERGE_FAST==-1
                FMeasureMerger.NoteIfNoMergerImprovement(fMeasure, maxF,...
                    bestIds(1));
            end
            bestIds=sort(bestIds);
        end

        function [bestIds, fMeasure]=FindBestUsingOriginalIdea(matchSize,...
                mergeIds, fMeasures4MergeIds, mergeSizes, truePositives)
            [fmBest2Worst, I]=sort(fMeasures4MergeIds, 'descend');
            idx=I(1);
            bestIds=mergeIds(idx);
            fMeasure=-100;% force at least ONE merger
            truePositive=truePositives(idx);
            mergeSize=mergeSizes(idx);
            truePositives(idx)=nan;
            tp=truePositive+truePositives;
            p=tp./(mergeSizes+mergeSize);
            r=tp./matchSize;
            [f, idx]=max((2*p.*r)./(p+r));
            if f>fMeasure
                fMeasure=f;
                truePositive=truePositive+truePositives(idx);
                truePositives(idx)=nan;
                mergeSize=mergeSize+mergeSizes(idx);
                bestIds(end+1)=mergeIds(idx);
                nMergeCandidates=length(fMeasures4MergeIds);
                for i=2:nMergeCandidates
                    if fMeasure<.5 && fmBest2Worst(i)<.5
                        %NOT going to get better!
                        break;
                    end
                    idx=I(i);
                    if isnan(truePositives(idx))
                        continue;
                    end
                    c=mergeSizes(idx);
                    tp=truePositive+truePositives(idx);
                    p=tp/(c+mergeSize);
                    r=tp/matchSize;
                    f=(2*p*r)/(p+r);
                    if f>fMeasure
                        fMeasure=f;
                        truePositive=tp;
                        mergeSize=c+mergeSize;
                        bestIds(end+1)=mergeIds(idx);
                    end
                end
            end

            if QfHiDM.F_MEASURE_MERGE_FAST==-1
                FMeasureMerger.NoteIfNoMergerImprovement(fMeasure,...
                    fmBest2Worst(1), bestIds(1));
            end
            bestIds=sort(bestIds);
        end

        function NoteIfNoMergerImprovement(fMeasure, maxF, bestId)
            if fMeasure<maxF
                fprintf(['No improvement by merging:  Best merger (fm=%s) is LESS than '...
                    'ID=%d (fm=%s)!\n'], ...
                    String.encodeRounded(fMeasure, 3), bestId, ...
                    String.encodeRounded(maxF(1),3));
            end
        end

        function [dif, newF, oldF]=DetectChange(matchSize, ...
                 currentTruePositives, currentMergeSize, newTruePositives, newMergeSize)
             p=currentTruePositives/currentMergeSize;
             r=currentTruePositives/matchSize;
             oldF=(2*p*r)/(p+r);
             p=newTruePositives/currentMergeSize;
             r=newTruePositives/matchSize;
             individualF=(2*p*r)/(p+r);
             tp=currentTruePositives+newTruePositives;
             p=tp/(newMergeSize+currentMergeSize);
             r=tp/matchSize;
             newF=(2*p*r)/(p+r);
             dif=newF-oldF;
             if dif<0
                 word='WORSER';
             else
                 word='better';
             end
             fprintf('%s F-measures:  old=%s, individual=%s, merged=%s, change=%s\n',...
                 word, ...
                 String.encodeRounded(oldF, 3), ...
                 String.encodeRounded(individualF, 3), ...
                 String.encodeRounded(newF, 3), ...
                 String.encodeRounded(dif, 3));
        end
    end
end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [ y, j, ierror ] = yj_check_common ( y, j )

%*****************************************************************************80
%
%% YJ_CHECK_COMMON checks a Common YJ date.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    25 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input/output, integer Y, J, the YJ date.
%
%    Output, integer IERROR, is 0 if no serious error was found in
%    the date, and 1 otherwise.
%

%
%  Check the year.
%
  [ y, ierror ] = y_check_common ( y );

  if ( ierror ~= 0 )
    return
  end
%
%  Make sure J is not too small or too big.
%
  [ y, j ] = j_borrow_common ( y, j );

  [ y, j ] = j_carry_common ( y, j );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%  Figure 7.39      Feedback Control of Dynamic Systems, 5e
%                        Franklin, Powell, Emami
%
% script to generate Fig. 7.39
clf;
np=1;
dp=[1 0 0];
nc=[1 0.619];
dc=[1 6.41];
num=conv(np,nc);
den=conv(dp,dc);
rlocus(num,den);
axis([-10 2 -4.5 4.5]);
grid;
title('Fig.7.39  Root locus for reduced-order controller for 1/s^2');


text:
function plotExcitationPhase(varargin)
% Plots the excitation force phase for each hydro structure's bodies in
% the heave, surge and pitch degrees of freedom.
%
% Usage:
% ``plotExcitationPhase(hydro, hydro2, hydro3, ...)``
%
% Parameters
% ----------
%     varargin : struct(s)
%         The hydroData structure(s) created by the other BEMIO functions.
%         One or more may be input.
%

if isempty(varargin)
    error(['plotExcitationPhase: No arguments passed. Include one or more hydro ' ...
        'structures when calling: plotExcitationPhase(hydro1, hydro2, ...)']);
end

B=1;  % Wave heading index
figHandle = figure('Position',[950,300,975,521]);
titleString = ['Excitation Force Phase: $$\phi_i(\omega,\theta)$$'];
subtitleString = {'Surge','Heave','Pitch'};
xString = {'$$\omega (rad/s)$$','$$\omega (rad/s)$$','$$\omega (rad/s)$$'};
yString = {['$$\phi_1(\omega,\theta$$',' = ',num2str(varargin{1}.theta(B)),'$$^{\circ})$$'],...
    ['$$\phi_3(\omega,\theta$$',' = ',num2str(varargin{1}.theta(B)),'$$^{\circ}$$)'],...
    ['$$\phi_5(\omega,\theta$$',' = ',num2str(varargin{1}.theta(B)),'$$^{\circ}$$)']};

notes = {''};

numHydro = length(varargin);
for ii = 1:numHydro
    numBod = varargin{ii}.Nb;
    tmp1 = strcat('X',num2str(ii));
    X.(tmp1) = varargin{ii}.w;
    tmp2 = strcat('Y',num2str(ii));
    a = 0;
    for i = 1:numBod
        m = varargin{ii}.dof(i);
        Y.(tmp2)(1,i,:) = squeeze(varargin{ii}.ex_ph(a+1,B,:));
        Y.(tmp2)(2,i,:) = squeeze(varargin{ii}.ex_ph(a+3,B,:));
        Y.(tmp2)(3,i,:) = squeeze(varargin{ii}.ex_ph(a+5,B,:));
        legendStrings{i,ii} = [varargin{ii}.body{i}];
        a = a + m;
    end
end

formatPlot(figHandle,titleString,subtitleString,xString,yString,X,Y,legendStrings,notes)
saveas(figHandle,'Excitation_Phase.png');

end

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
function ompspeedtest
%OMPSPEEDTEST Test the speed of the OMP functions.
%  OMPSPEEDTEST invokes the three operation modes of OMP and compares
%  their speeds. The function automatically selects the number of signals
%  for the test based on the speed of the system.
%
%  To run the test, type OMPSPEEDTEST from the Matlab prompt.
%
%  See also OMPDEMO.

%  Ron Rubinstein
%  Computer Science Department
%  Technion, Haifa 32000 Israel
%  ronrubin@cs
%
%  April 2009

% random dictionary %

n = 500;
L = 1000;
T = 20;

D = randn(n,L);
D = D*diag(1./sqrt(sum(D.*D)));    % normalize the dictionary

% select signal number according to computer speed %

x = randn(n,5);
tic; omp(D,x,[],T,'messages',-1); t=toc;
signum = ceil(20/(t/5));     % approximately 20 seconds of OMP-Cholesky

% generate random signals %

X = randn(n,signum);

% run OMP  %

printf('\nRunning OMP-Cholesky...');
tic; omp(D,X,[],T,'messages',4); t1=toc;

printf('\nRunning Batch-OMP...');
tic; omp(D,X,D'*D,T,'messages',1); t2=toc;

printf('\nRunning Batch-OMP with D''*X specified...');
tic; omp(D'*X,D'*D,T,'messages',1); t3=toc;

% display summary  %

printf('\n\nSpeed summary for %d signals, dictionary size %d x %d:\n', signum, n, L);
printf('Call syntax        Algorithm               Total time');
printf('--------------------------------------------------------');
printf('OMP(D,X,[],T)      OMP-Cholesky            %5.2f seconds', t1);
printf('OMP(D,X,G,T)       Batch-OMP               %5.2f seconds', t2);
printf('OMP(DtX,G,T)       Batch-OMP with D''*X     %5.2f seconds\n', t3);

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function showresult(node,elem,u,viewangle)
%% SHOWRESULT display the mesh and the solution
%
%  showresult(node,elem,u,viewangle) displays the mesh and the solution in
%  one figure. The left one is the mesh, the middle
%  one is the contour of the solution, and the right one is the graph of
%  the function. The last viewangle is used to adjust the view angle of the
%  graph of the function.
%
%  Example:
%     f = inline('sin(2*pi*x).*cos(2*pi*y)');
%     node = [0,0; 1,0; 1,1; 0,1];
%     elem = [2,3,1; 4,1,3];
%     for k = 1:4
%         [node,elem] = uniformrefine(node,elem);
%     end
%     u = f(node(:,1),node(:,2));
%     showresult(node,elem,u,[-62,58]);
%
% See also showrate, showmesh, showsolution
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

if (length(u) == size(elem,1)) || (length(u) == size(node,1))
    % show mesh
    set(gcf,'Units','normal');
    set(gcf,'Position',[0.25,0.25,0.6,0.25]);
    if size(elem,1) < 6e4
        subplot(1,3,1);
        showmesh(node,elem);
        pause(0.05)
    else
        subplot(1,3,1);
        title('The mesh is too dense to display')
    end
    % show solution
    subplot(1,3,2);
    showsolution(node,elem,u,2);
    colorbar;
    pause(0.05)
    subplot(1,3,3);
    showsolution(node,elem,u);
    if nargin>3
        view(viewangle);
    end
    pause(0.05)
else
    showmesh(node,elem);
end

repo: Aero-matlab


text:
function c=comp_uwpfbt(f,wtNodes,rangeLoc,nodesUps,scaling,interscaling)
%COMP_UWPFBT Compute Undecimated Wavelet Packet Filterbank Tree
%   Usage:  c=comp_uwpfbt(f,wtNodes,nodesUps);
%
%   Input parameters:
%         f        : Input data as L*W array.
%         wtNodes  : Filterbank tree nodes (elementary filterbanks) in
%                    BF order. Cell array of structures of length *nodeNo*.
%         nodesUps : Filters upsampling factor of each node. Array of
%                    length *nodeNo*.
%
%   Output parameters:
%         c        : Coefficients stored in L*M*W array.
%

% Pre-allocated output
[L, W] = size(f);
M = sum(cellfun(@(wtEl) numel(wtEl.h),wtNodes));
c = zeros(L,M,W,assert_classname(f,wtNodes{1}.h{1}.h));

% Convenience input reshape
ca = reshape(f,size(f,1),1,size(f,2));
cOutRunIdx = 1;
cInRunIdxs = [1];

interscalingfac = 1;
if strcmp('intscale',interscaling)
    interscalingfac = 1/2;
elseif strcmp('intsqrt',interscaling)
    interscalingfac = 1/sqrt(2);
end

% For each node in tree in the BF order...
for jj=1:numel(wtNodes)
   % Node filters subs. factors
   a = wtNodes{jj}.a;

   % Optionally scale the filters
   h = comp_filterbankscale(wtNodes{jj}.h(:),a(:),scaling);

   % Node filters to a matrix
   % hMat = cell2mat(cellfun(@(hEl) conj(flipud(hEl.h(:))),h','UniformOutput',0));
   hMat = cell2mat(cellfun(@(hEl) hEl.h(:),h','UniformOutput',0));

   % Node filters initial skips
   % hOffet = cellfun(@(hEl) 1-numel(hEl.h)-hEl.offset,wtNodes{jj}.h);
   hOffet = cellfun(@(hEl) hEl.offset, wtNodes{jj}.h);
   % Number of filters of the current node
   filtNo = size(hMat,2);
   % Zero index position of the upsampled filters.
   offset = nodesUps(jj).*(hOffet);

   % Run filterbank
   c(:,cOutRunIdx:cOutRunIdx + filtNo-1,:)=...
      comp_atrousfilterbank_td(squeeze(ca(:,1,:)),hMat,nodesUps(jj),offset);

   % Bookkeeping
   outRange = cOutRunIdx:cOutRunIdx+filtNo-1;
   outRange(rangeLoc{jj}) = [];
   cInRunIdxs = [cInRunIdxs(2:end),outRange];

   cOutRunIdx = cOutRunIdx + filtNo;

   % Prepare input for the next iteration
   if ~isempty(cInRunIdxs)
      c(:,cInRunIdxs(1),:) = c(:,cInRunIdxs(1),:)*interscalingfac;
      ca = c(:,cInRunIdxs(1),:);
   end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
%% This file is used to generate the simulation data of the Schwarzian-Kroteweg-de Vires equation.
% Coded By: K
% Last Updated: 2019/06/17
%% Clear all
clc;clear;close all;

%% Create a folder
addpath('Functions')
%
[fld_status, fld_msg, fld_msgID]=mkdir('Datas');
addpath('Datas')
%
[fld_status, fld_msg, fld_msgID]=mkdir('TempFunctions');
addpath('TempFunctions')
%% Define the parameters for the SKdV equation
% Define the time span
dt=0.01;T=20;
tspan=0:dt:T;

% Define the spatial domain
L=50;

% Define the discretization point
n=256;

% Get the index of each point
x2=linspace(-L/2,L/2,n+1);
x=x2(1:n);
dx=x(2)-x(1);

% Get the frequency
k=(2*pi/L)*[0:n/2-1 -n/2:-1].';
k2=fftshift(k);

% Get the initial value
a1=0.5;
u1=(a1/2*(sech(0.5*sqrt(1)*(x+10))).^2 ).';

a2=0.5;
u2=(a2/2*(sech(0.5*sqrt(1)*(x-10))).^2 ).';

u0=u1+u2;

% Transfer the initial point to spectram domain
ut0=fft(u0);

% Define parameter for system
% Omega: This term should be positive, and it will determine the
% dissapation of the system
Omega = 0.1;

alpha = 1;

% g0: This term should be positive, and it determines the gain of the
% system.
g0 =[0;0.0001;0.001;0.01;0.1;1];

%% Loop through the g0
for kkk=1:size(g0,1)
    % Solve for the PDE
    tic
    opts = odeset('RelTol',1e-12,'AbsTol',1e-13);
    [time,utsol]=ode45(@(time,u0)mkdv_rhs(time,u0,k,x,Omega,g0(kkk,1),alpha),tspan,ut0,opts);
    toc

    % Back to time domain
    for j=1:length(tspan)
        usol(:,j)=ifft(utsol(j,1:n).');
    end

    u=zeros(size(usol));
    ut=zeros(size(usol));
    ux=zeros(size(usol));
    uxx=zeros(size(usol));
    uxxx=zeros(size(usol));

    % Use the previous utsol to get the derivative
    for j=1:length(tspan)
        [~,ut(:,j),u(:,j),ux(:,j),uxx(:,j),uxxx(:,j)]=mkdv_rhs(0,utsol(j,1:n)',k,x,Omega,g0(kkk,1),alpha);
    end

    u=u';ut=ut';ux=ux';uxx=uxx';uxxx=uxxx';

    name=strcat('Datas/','SimulationData_V_',num2str(kkk),'.mat');
    save(name,'u','ut','ux','uxx','uxxx','x','tspan','dt','T','L','n');

end

%% Plot the result
figure(1)
waterfall(x,tspan,real(u))
colormap([0 0 0]);
view(42,55)

figure(2)
surfl(x,tspan,real(u))
colormap(gray)
shading interp
view(42,55)

figure(3)
pcolor(x,tspan,real(u))
shading interp

author: dynamicslab

repo: SINDy-PI


text:
function title = p28_title ( )

%*****************************************************************************80
%
%% P28_TITLE returns the title of problem P28.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'Problem 28, Enright and Pryce #F3';

  return
end


text:
function fatrix2_block_free(arg)
%printm 'freeing fatrix2_block object static memory'
for mm=1:length(arg.blocks)
	try
		free(blocks{mm})
	catch
	end
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function keyinfo
%% KEYINFO catch all the keyboard infomation
% The function requires no input arguments and return a value in terms of
% a column vector [3 x 1]
%
% Finally, I wish you will enjoy this tiny function, and any comment is
% appreciated
%
% Name: Fu-Sung Wang
% E-mail: fusung.wang@gmail.com
% Organization: www.terasoft.com.tw


text:
function [MW] = hpb2MW(hpb)
% Convert power from boiler horsepower to megawatts.
% Chad A. Greene 2012
MW = hpb*9809.5e-6;


text:
%% Structured grid
% This is a sample code showing how to use ConstructProjector2D().
% This example is going to setup a transformation matrix based on forth
% order polynomials in 2D.
% Pay attention that no data is needed to construct the interpolator. Only
% the coordinates of the points on source and destination grid.
clear;clc;
%% Initializing Part I
disp('- Initializing')
xMin=0;
xMax=2*pi;
yMin=0;
yMax=2*pi;

nx=40;
ny=40;

nPoly=4; % We gonna fit a fourth order polynomial,
         % i.e. sum_{n,m=0}^{n,m=nPoly} x^n*y^m
nInterp=30; % nPoly=4 requires 25 points. However,
            % we are going to use 30 points. This would be
            % the least-square fit of the surface.
            % Note that: min(nInterp)=(nPoly+1)^2. Otherwise it won't work.
%% Turning off the warnings
% Depending on the parameter that you choose the matrices might be nearly
% singular. In this example, despite being nearly singular, they are not
% affecting the final interpolation. Therefore, I set the warning to off.
% But generally do not turn it off and if you got this warning check if it
% is going to affect your results or not.
disp('- Turning the warnings off')
warning('OFF','MATLAB:nearlySingularMatrix');

%% Generating the Source Grid
disp('- Generating the source grid.')
[xn,yn]=meshgrid(linspace(xMin,xMax,nx),linspace(yMin,yMax,ny));

%% Finding the cell centers
disp('- Generating the destination grid')
xc=(xn(1:ny-1,1:nx-1)+xn(2:ny,1:nx-1)+xn(1:ny-1,2:nx)+xn(2:ny,2:nx))*0.25;
yc=(yn(1:ny-1,1:nx-1)+yn(2:ny,1:nx-1)+yn(1:ny-1,2:nx)+yn(2:ny,2:nx))*0.25;

figure
surface(xn,yn,zeros(size(xn)));
hold on
plot(xn,yn,'k.');
plot(xc,yc,'b.');
axis tight;
axis square;
title('Source Grid, black dots, & the destination grid, blue dots.', ...
      'FontName','Arial','FontSize',12,'FontWeight','Bold');
%% Constructing the Interpolant
% Note: that we do not need the data on the source grid to create the
% interpolant.
disp('- Constructing the interpolant')
P=ConstructProjector2D(xn(:),yn(:),xc(:),yc(:),nPoly,nInterp);

%% Generarting some Data
disp('- Generating some data and interpolating')
F1=@(x,y) (sin(sqrt(x.^2+y.^2)));
zn=F1(xn,yn);
zc_interp=reshape(P*zn(:),size(xc));
zc_Analytic=F1(xc,yc);
RMSE=sqrt(mean((zc_Analytic(:)-zc_interp(:)).^2));

figure
surface(xc,yc,zc_interp,'EdgeColor','none');
title(['nPoly: ' num2str(nPoly) ', RMSE= ' num2str(RMSE)]);
axis tight

%% Generating some more data
disp('- Generating multiple data field on the source grid and interpolating them.')
F2=@(x,y) (sin(x).*cos(y));
F3=@(x,y) (exp(-sqrt(x.^2+y.^2)));
F4=@(x,y,x0,y0) (exp(-sqrt((x-x0).^2+(y-y0).^2)));

zn(:,:,1)=F1(xn,yn);
zn(:,:,2)=F2(xn,yn);
zn(:,:,3)=F3(xn,yn);
zn(:,:,4)=F4(xn,yn,mean(mean(xn)),mean(mean(yn)));

zc_Analytic(:,:,1)=F1(xc,yc);
zc_Analytic(:,:,2)=F2(xc,yc);
zc_Analytic(:,:,3)=F3(xc,yc);
zc_Analytic(:,:,4)=F4(xc,yc,mean(mean(xn)),mean(mean(yn)));

% Now interpolating
% Note that the same interpolant is used for all data fields and they can
% be all interpolated with one sparse matrix multiplication.
zc_interp=P*reshape(zn,nx*ny,4); % There are 4 data fields.
zc_interp=reshape(zc_interp,(nx-1),(ny-1),4); % these two commands can be combined in one.
                                              % They were separated for clarity.

% calculating the RMSE and plotting
RMSE=zeros(4,1);
figure
for i=1:4
  RMSE(i)=sqrt(mean((reshape(zc_Analytic(:,:,i),(nx-1)*(ny-1),1)-reshape(zc_interp(:,:,i),(nx-1)*(ny-1),1)).^2));
  subplot(2,2,i);
  surface(xc,yc,squeeze(zc_interp(:,:,i)),'EdgeColor','none');
  title(['F' num2str(i) ', nPoly: ' num2str(nPoly) ', RMSE:' num2str(RMSE(i))])
  axis tight
  axis square
end


text:
% Sparse covariance estimation for Gaussian variables
% Joëlle Skaf - 04/24/08
% (a figure is generated)
%
% Suppose y \in\reals^n is a Gaussian random variable with zero mean and
% covariance matrix R = \Expect(yy^T), with sparse inverse S = R^{-1}
% (S_ij = 0 means that y_i and y_j are conditionally independent).
% We want to estimate the covariance matrix R based on N independent
% samples y1,...,yN drawn from the distribution, and using prior knowledge
% that S is sparse
% A good heuristic for estimating R is to solve the problem
%           maximize    logdet(S) - tr(SY) - lambda*sum(sum(abs(S)))
%           subject to  S >= 0
% where Y is the sample covariance of y1,...,yN, and lambda is a sparsity
% parameter to be chosen or tuned.
% A figure showing the sparsity (number of nonzeros) of S versus lambda
% is generated.

% Input data
rand('state',0); %#ok
randn('state',0); %#ok
n = 10;
N = 100;
Strue = sprandsym(n,0.5,0.01,1);
Strue(abs(Strue(:))<=1e-4) = 0;
nnz_true = nnz(Strue);
R = inv(full(Strue));
y_sample = sqrtm(R)*randn(n,N);
Y = cov(y_sample');
Nlambda = 20;
lambda = logspace(-2, 3, Nlambda);
nnz_i = zeros(1,Nlambda);

for i=1:Nlambda
    disp(['i = ' num2str(i) ', lambda(i) = ' num2str(lambda(i))]);
    % Maximum likelihood estimate of R^{-1}
    cvx_begin sdp quiet
        variable S(n,n) symmetric
        maximize log_det(S) - trace(S*Y) - lambda(i)*sum(sum(abs(S)))
        S >= 0; %#ok
    cvx_end
    nnz_i(i) = sum(abs(S(:))>1e-4);
end

figure;
semilogx(lambda, nnz_i);
hold on;
semilogx(lambda, nnz_true*ones(1,Nlambda),'r');
xlabel('\lambda');
legend('nonzeros in S', 'nonzeros in R^{-1}');

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function [bestError bestMap] = compare_labels(correctLabels, computedLabels)

% compare_labels
%
%   Given two labelings of a set of data, find a map between the
%   labelings that minimizes the number of samples that they assign
%   different labels. This version assumes the correct labeling has a small
%   number of groups (<= 5), and does an exhaustive search to find
%   the best map. If the correct number of groups is larger, the result
%   from relabel_samples (which will likely be incorrect) is run.
%
% Inputs:
%   correctLabels - the true labels of some data.
%   computedLabels - the labels computed by some algorithm for that data.
%
% Outputs:
%   bestError - the misclassification rate obtained by this relabelling.
%   bestMap - the map that converts correct labels to computed labels
%
% Sep. '07  Shankar Rao -- srrao@uiuc.edu

% Copyright 2007, University of Illinois. All rights reserved.

correctLabels = correctLabels(:)';
computedLabels = computedLabels(:)';

groupCount = max(correctLabels);
computedGroupCount = max(computedLabels);
sampleCount = length(correctLabels);

maxGroupCount = max(groupCount, computedGroupCount);

if groupCount <= 5
    bestError = 1;
    bestMap = [1:groupCount];
    candidateLabels = nchoosek(1:maxGroupCount, groupCount);
    permCount = factorial(groupCount);
    for candidateIndex = 1:size(candidateLabels,1)
        maps = perms(candidateLabels(candidateIndex,:));
        for permIndex = 1:permCount
            map = maps(permIndex,:);
            currentError = mean(map(correctLabels) ~= computedLabels);
            if currentError < bestError
                bestError = currentError;
                bestMap = map;
            end
        end
    end
else
    [bestMap, bestError] = relabel_samples(correctLabels, computedLabels, ones(1,groupCount));
end

author: SuTanTank

repo: VideoStitchingViaShakinessRemoving


text:
function [ node_coord, node_att, node_marker ] = triangle_node_data_example ( ...
  node_num, node_dim, node_att_num, node_marker_num )

%*****************************************************************************80
%
%% TRIANGLE_NODE_DATA_EXAMPLE returns the node information for the example.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    11 October 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NODE_NUM, the number of nodes.
%
%    Input, integer NODE_DIM, the spatial dimension.
%
%    Input, integer NODE_ATT_NUM, number of node attributes
%    listed on each node record.
%
%    Input, integer NODE_MARKER_NUM, 1 if every node record
%    includes a final boundary marker value.
%
%    Output, real NODE_COORD(NODE_DIM,NODE_NUM), the nodal
%    coordinates.
%
%    Output, real NODE_ATT(NODE_ATT_NUM,NODE_NUM), the nodal
%    attributes.
%
%    Output, integer NODE_MARKER(NODE_MARKER_NUM,NODE_NUM),
%    the node markers.
%
  node_coord = [ ...
    0.0, 0.0; ...
    1.0, 0.0; ...
    2.0, 0.0; ...
    3.0, 0.0; ...
    4.0, 0.0; ...
    0.0, 1.0; ...
    1.0, 1.0; ...
    2.0, 1.0; ...
    3.0, 1.0; ...
    4.0, 1.0; ...
    0.0, 2.0; ...
    1.0, 2.0; ...
    2.0, 2.0; ...
    3.0, 2.0; ...
    4.0, 2.0; ...
    0.0, 3.0; ...
    1.0, 3.0; ...
    2.0, 3.0; ...
    0.0, 4.0; ...
    1.0, 4.0; ...
    2.0, 4.0 ]';

  node_att = [];

  node_marker = [ ...
    1; 1; 1; 1; 1; 1; 0; 0; 0; 1; ...
    1; 0; 0; 1; 1; 1; 0; 1; 1; 1; ...
    1 ]';

  return
end


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXERCISE: WRITE THE LAGRANGE FORMULATION FOR A 2DOF PLANAR ROBOT ARM
%   fext are expressed in the base reference system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) 2012, by Arturo Gil Aparicio
%
% This file is part of ARTE (A Robotics Toolbox for Education).
%
% ARTE is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% ARTE is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with ARTE.  If not, see <http://www.gnu.org/licenses/>.
function tau = exercise_lagrange_2dofplanar(robot, q, qd, qdd, g, fext)
a = eval(robot.DH.a);
m=robot.dynamics.masses(1);

%DEFINE M

%DEFINE V

%DEFINE G

%compute tau!
%tau = M*Qdd + V + G

repo: ARTE


text:
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  calc_VIDEVAL_feat.m
%
%  This function extract the VIDEVAL feature vector (60-dim)
%
%
%  Input:
%           test_video:             Path to the test video file (YUV420 format)
%           width, height:          Resolution of the YUV video
%           framerate:              number of frames per second
%
%  Output:
%           VIDEVAL_all_features:   Resulting VIDEVAL feature vector (60-dim)
%
% Note: this is a light version of VIDEVAL, which has been spatially and
%       temporally sped up.
%
%  - Resized frames to 720p if larger than 720p keeping aspect ratio
%  - Sampled at most 6 frames for each 1-second block

%%
function [VIDEVAL_all_features] = calc_VIDEVAL_feats_light(test_video, ...
        width, height, framerate, max_reso, frs_per_blk)
    % setup speededup params
%     max_reso = 480;
%     frs_per_blk = 3;

    % subsampling ratio
    ratio = max_reso / min(width, height);
    fr_step = max(floor(framerate / frs_per_blk),1);

    % Try to open test_video; if cannot, return
    test_file = fopen(test_video,'r');
    if test_file == -1
        fprintf('Test YUV file not found.');
        VIDEVAL_all_features = [];
        return;
    end
    % Open test video file
    fseek(test_file, 0, 1);
    file_length = ftell(test_file);
    fprintf('Video file size: %d bytes (%d frames)\n',file_length, ...
            floor(file_length/width/height/1.5));
    % get frame number
    frame_start = 1;
    frame_end = (floor(file_length/width/height/1.5)-3);
    first_frame_loaded = 0;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Loop through all the frames in the frame_range to compute frame features
    %
    fprintf('Computing frame features on frames %d..%d\n', ...
            frame_start, frame_end);

    frame_features_all = [];
    for i = frame_start:fr_step:frame_end
        % Read frames i-i, i and i+1 (note that frame_start must be > 0)
        if first_frame_loaded
            % prev_YUV_frame = next_YUV_frame;
            prev_YUV_frame = YUVread(test_file,[width height],i-1);
            this_YUV_frame = YUVread(test_file,[width height],i);
            next_YUV_frame = YUVread(test_file,[width height],i+1);
            if ratio < 1
                this_YUV_frame = imresize(this_YUV_frame, ratio);
                next_YUV_frame = imresize(next_YUV_frame, ratio);
            end
        else
            prev_YUV_frame = YUVread(test_file,[width height],i-1);
            this_YUV_frame = YUVread(test_file,[width height],i);
            next_YUV_frame = YUVread(test_file,[width height],i+1);
            first_frame_loaded = 1;
            if ratio < 1
                prev_YUV_frame = imresize(prev_YUV_frame, ratio);
                this_YUV_frame = imresize(this_YUV_frame, ratio);
                next_YUV_frame = imresize(next_YUV_frame, ratio);
            end
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Compute BRISQUE features # 5(std), 9, 23, 29
        % features_all(1:4) #4
        brisque_feats = compute_brisque_features(this_YUV_frame(:,:,1));

        % compute GM-LOG features # 2(std), 4(std), 8, 18(std), 20, 23, 34(std), 40(both)
        % features_all(5:12) #8
        gmlog_feats = compute_gmlog_features(this_YUV_frame(:,:,1));

        % compute HIGRADE features # 3, 13, 14, 16, 21(both), 28, 30, 33
        % features_all(13:20) #8
        higrade_feats = compute_higrade_features(this_YUV_frame(:,:,1));

        % compute FRIQUEE_lUMA features # 5(std), 8(std), 27(std), 32(std), 58, 62(std), 66, 68(std), 70(both), 71, 72(std)
        % features_all(21:31) #11
        friquee_luma_feats = compute_friquee_luma_features(this_YUV_frame(:,:,1));

        % compute FRIQUEE_CHRO features # 4  7(std) 15   24   26   47   60   62   72   73   79
        % features_all(32:42) #11
        friquee_chroma_feats = compute_friquee_chroma_features(this_YUV_frame(:,:,:));

        % compute FRIQUEE LMS features # 52
        % features_all(43:43) # 1
        friquee_lms_feats = compute_friquee_lms_features(this_YUV_frame(:,:,:));

        % compute TLVQM LCF features # 1 5(both) 12 16 18 19(std) 22(std)
        % features_all(44:50) # 7
        tlvqm_lcf_feats = compute_tlvqm_lcf_features(this_YUV_frame./255, ...
                                                     prev_YUV_frame./255, ...
                                                     next_YUV_frame./255);
        % concat features
        frame_features_all = [frame_features_all; [brisque_feats, gmlog_feats, higrade_feats, ...
            friquee_luma_feats, friquee_chroma_feats, friquee_lms_feats, tlvqm_lcf_feats]];
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Compute mean and std pooling
    %
    std_feats = [];
    avg_feats = [];
    std_pool_idx = [1, 5, 6, 8, 11, 12, 17, 21, 22, 23, 24, 26, 28, 29, 31, ...
        33, 45, 49, 50];
    avg_pool_idx = [2:4, 7, 9, 10, 12, 13:20, 25, 27, 29, 30, 32, 34:48];
    n_temp_vecs = length(frame_features_all(:,1));
%     half_framerate = floor(framerate/2);

    fprintf('Mean and standard deviation pooling of frame-wise features...\n');
    for j = 1:frs_per_blk:n_temp_vecs - frs_per_blk
        j_start = j;
        j_end = min(n_temp_vecs, j + frs_per_blk);
        % mean and std pooling
        std_feats = [std_feats; nanstd(frame_features_all(j_start:j_end, std_pool_idx))];
        avg_feats = [avg_feats; nanmean(frame_features_all(j_start:j_end, avg_pool_idx))];
    end

    % compute TLVQM HCF features # 1 2 17 22 24 30
    tlvqm_hcf_feats = [];
    fprintf('Computing TLVQM-HCF features at 1fps...\n');
    for fr=floor(framerate/2):framerate:frame_end
        YUV_frame = YUVread(test_file,[width height],fr-1);
        if ratio < 1
            YUV_frame = imresize(YUV_frame, ratio);
        end
        ftrs = compute_tlvqm_hcf_features(YUV_frame./255);
        tlvqm_hcf_feats = [tlvqm_hcf_feats; ftrs];
    end

    % Combine feature vectors
    VIDEVAL_all_features = [nanmean(avg_feats, 1)   ...
                    nanmean(std_feats, 1) ...
                    nanmean(tlvqm_hcf_feats, 1)];

    fclose(test_file);
    fprintf('Extracting VIDEVAL features for %s finished!\n', test_video);
end

%% compute_brisque_features
function feat = compute_brisque_features(imgray)
% Compute BRISQUE features # 5(std), 9, 23, 29
    window = fspecial('gaussian',7,7/6);
    window = window/sum(sum(window));
    feat = [];
    %% scale1 5, 9
    mu = filter2(window, imgray, 'same');
    mu_sq = mu.*mu;
    sigma = sqrt(abs(filter2(window, imgray.*imgray, 'same') - mu_sq));
    imstruct = (imgray-mu)./(sigma+1);
    shifts = [0 1; 1 0];
    for itr_shift =1:2
        shifted_imstruct = circshift(imstruct,shifts(itr_shift,:));
        pair = imstruct(:).*shifted_imstruct(:);
        [~, leftstd, ~] = estimateaggdparam(pair);
        feat =[feat leftstd^2];
    end

    %% scale2 5, 11
    imgray = imresize(imgray,0.5);
    mu = filter2(window, imgray, 'same');
    mu_sq = mu.*mu;
    sigma = sqrt(abs(filter2(window, imgray.*imgray, 'same') - mu_sq));
    imstruct = (imgray-mu)./(sigma+1);
    shifts = [0 1; 1 1];
    for itr_shift =1:2
        shifted_imstruct = circshift(imstruct,shifts(itr_shift,:));
        pair = imstruct(:).*shifted_imstruct(:);
        [alpha, leftstd, ~] = estimateaggdparam(pair);
        if itr_shift == 1
            feat = [feat leftstd^2];
        else
            feat = [feat alpha];
        end
    end
end

%% compute_gmlog_features
function feat = compute_gmlog_features(imgray)
% compute GM-LOG features # 2(std), 4(std), 8, 18(std), 20, 23, 34(std), 40(both)
    sigma = 0.5;
    [gx,gy] = gaussian_derivative(imgray,sigma);
    grad_im = sqrt(gx.^2+gy.^2);

    window2 = fspecial('log', 2*ceil(3*sigma)+1, sigma);
    window2 =  window2/sum(abs(window2(:)));
    log_im = abs(filter2(window2, imgray, 'same'));

    ratio = 2.5; % default value 2.5 is the average ratio of GM to LOG on LIVE database
    grad_im = abs(grad_im/ratio);

    %Normalization
    c0 = 4*0.05;
    sigmaN = 2*sigma;
    window1 = fspecial('gaussian',2*ceil(3*sigmaN)+1, sigmaN);
    window1 = window1/sum(window1(:));
    Nmap = sqrt(filter2(window1,mean(cat(3,grad_im,log_im).^2,3),'same'))+c0;
    grad_im = (grad_im)./Nmap;
    log_im = (log_im)./Nmap;
    % remove the borders, which may be the wrong results of a convolution
    % operation
    h = ceil(3*sigmaN);
    grad_im = abs(grad_im(h:end-h+1,h:end-h+1,:));
    log_im = abs(log_im(h:end-h+1,h:end-h+1));

    ctrs{1} = 1:10;ctrs{2} = 1:10;
    % histogram computation
    step1 = 0.20;
    step2 = 0.20;
    grad_qun = ceil(grad_im/step1);
    log_im_qun = ceil(log_im/step2);

    N1 = hist3([grad_qun(:),log_im_qun(:)],ctrs);
    N1 = N1/sum(N1(:));
    NG = sum(N1,2); NL = sum(N1,1);

    alpha1 = 0.0001;
    % condition probability: Grad conditioned on LOG
    cp_GL = N1./(repmat(NL,size(N1,1),1)+alpha1);
    cp_GL_H=  sum(cp_GL,2)';
    cp_GL_H = cp_GL_H/sum(cp_GL_H);
    % condition probability: LOG conditioned on Grad
    cp_LG = N1./(repmat(NG,1,size(N1,2))+alpha1);
    cp_LG_H = sum(cp_LG,1);
    cp_LG_H = cp_LG_H/(sum(cp_LG_H));

    out = [NG', NL, cp_GL_H,cp_LG_H];

    feat = out([2, 4, 8, 18, 20, 23, 34, 40]);
end

function [gx,gy] = gaussian_derivative(imd,sigma)
    window1 = fspecial('gaussian',2*ceil(3*sigma)+1+2, sigma);
    winx = window1(2:end-1,2:end-1)-window1(2:end-1,3:end);winx = winx/sum(abs(winx(:)));
    winy = window1(2:end-1,2:end-1)-window1(3:end,2:end-1);winy = winy/sum(abs(winy(:)));
    gx = filter2(winx,imd,'same');
    gy = filter2(winy,imd,'same');
end

%% compute_higrade_features
function feat = compute_higrade_features(imgray)
% compute HIGRADE features # 3, 13, 14, 16, 21(both), 28, 30, 33
[Gmag, ~] = imgradient(imgray);
window = fspecial('gaussian',7,7/6);
window = window/sum(sum(window));
scale_num = 2;
feat = [];
for scale = 1:scale_num
    if scale == 1
        im = Gmag;
%         ggdpara = [];
        shifts = [0 1;1 0 ; 1 1; -1 1];

        mu = filter2(window, im, 'same');
        mu_sq = mu.*mu;
        sigma = sqrt(abs(filter2(window, im.*im, 'same') - mu_sq));
        structdis = (im-mu)./(sigma+1);
        [sigma_tmp, ~] = gaussian_para_esti(structdis(:));
        feat = [feat sigma_tmp]; %3
        structdis = log(abs(structdis) + .1);
        for itr_shift = 1:4
            structdis_shift_tmp = circshift(structdis, shifts(itr_shift,:));
            structdis_diff = structdis - structdis_shift_tmp;
            structdis_shift(itr_shift) = {[structdis_shift_tmp]};
        end
        structdis_diff = structdis + structdis_shift{1,3} - structdis_shift{1,1} - structdis_shift{1,2};
        [sigma_tmp_0, alpha_tmp_0] = gaussian_para_esti(structdis_diff(:));
        feat = [feat  sigma_tmp_0 alpha_tmp_0]; % 13 14
        win_tmp_1 = [0 1 0; -1 0 -1; 0 1 0;];
        structdis_diff_1 = filter2(win_tmp_1, structdis, 'same');
        [~, alpha_tmp_1] = gaussian_para_esti(structdis_diff_1(:));%16
        feat = [feat alpha_tmp_1];
    else
        im = Gmag;
        shifts = [0 1;1 0 ; 1 1; -1 1];

        mu = filter2(window, im, 'same');
        mu_sq = mu.*mu;
        sigma = sqrt(abs(filter2(window, im.*im, 'same') - mu_sq));
        structdis = (im-mu)./(sigma+1);
        [sigma_tmp, ~] = gaussian_para_esti(structdis(:));
        feat = [feat sigma_tmp]; %3

        for itr_shift = 1:4
            structdis_shift_tmp = circshift(structdis, shifts(itr_shift,:));
            if itr_shift == 3
                structdis_diff = structdis - structdis_shift_tmp;
                [~, alpha_tmp] = gaussian_para_esti(structdis_diff(:));
                feat = [feat alpha_tmp]; % 10
            elseif itr_shift == 4
                structdis_diff = structdis - structdis_shift_tmp;
                [~, alpha_tmp] = gaussian_para_esti(structdis_diff(:));
                feat = [feat alpha_tmp]; % 12
            end
            structdis_shift(itr_shift) = {[structdis_shift_tmp]};
        end
        win_tmp_1 = [0 1 0; -1 0 -1; 0 1 0;];
        structdis_diff_1 = filter2(win_tmp_1, structdis, 'same');
        [sigma_tmp_1, ~] = gaussian_para_esti(structdis_diff_1(:));%15
        feat = [feat sigma_tmp_1];
    end
    Gmag = imresize(Gmag, 0.5);
end
end

%% compute_friquee_luma_features
function feat = compute_friquee_luma_features(imgray)
% # 5(std), 8(std), 27(std), 32(std), 58, 62(std), 66, 68(std), 70(both), 71, 72(std)
% 1:31 32:62

    scalenum=2;
    imGray1=imgray;

    % The array that aggregates all the features from different feature
    % maps.
    feat = [];
    %% 1 scale 31 features x2 = 62
    %% scale1: 5, 8, 27
    [structdis,~] = divisiveNormalization(imgray);
    shifts = [ 0 1;1 0];
    for itr_shift =1:2
        % Construct the product neighborhood map.
        shifted_structdis = circshift(structdis,shifts(itr_shift,:));
        pair = structdis(:).*shifted_structdis(:); % Element wise product.
        if itr_shift == 1
            feat = [feat skewness(pair(:))]; %5
        else
                    % Fit an AGGD and extract its parameters.
            [alpha, leftstd, rightstd] = estimateaggdparam(pair);
            const = (sqrt(gamma(1/alpha))/sqrt(gamma(3/alpha)));
            meanparam = (rightstd-leftstd)*(gamma(2/alpha)/gamma(1/alpha))*const;
            feat = [feat meanparam]; %8
        end
    end
    %Extract the sigma field from the image.
    sigmaMap = computeSigmaMap(imgray);

    % Compute kurtosis, skewness and average of the sigma field
    feat = [feat mean2(sigmaMap)];%27

    %% sclae2: 32,58,62 - 31 = 1, 27,31
    imgray = imresize(imgray,0.5);
    shifts = [ 0 1];
    for itr_shift =1:1
        % Construct the product neighborhood map.
        shifted_structdis = circshift(structdis,shifts(itr_shift,:));
        pair = structdis(:).*shifted_structdis(:); % Element wise product.

        % Fit an AGGD and extract its parameters.
        [alpha, ~, ~] = estimateaggdparam(pair);
%         const = (sqrt(gamma(1/alpha))/sqrt(gamma(3/alpha)));
%         meanparam = (rightstd-leftstd)*(gamma(2/alpha)/gamma(1/alpha))*const;
        feat = [feat alpha];%% 1
    end
     %Extract the sigma field from the image.
    sigmaMap = computeSigmaMap(imgray);

    % Compute kurtosis, skewness and average of the sigma field
    feat = [feat mean2(sigmaMap)];%27

      % Apply DNT
    [imDivNorm,~] = divisiveNormalization(imgray);
    feat = [feat skewness(imDivNorm(:))]; %31

    %%
    imgray=imGray1;
    %% Original Scales

    % 63:64
    % Features from the Yellow color channel map.
%     yFeat=yellowColorChannelMap(rgb);
%     feat = [feat yFeat];

    % 65:70 --> 63:68
    % Features from the Difference of Gaussian (DoG) of Sigma Map

    %Get the sigma map of the image.
    sigmaMap = computeSigmaMap(imgray);

    % Construct two Gaussian windows where sig2 = 1.5*sig1;
    k=1.5;

    window1 = fspecial('gaussian',7,7/6);
    window1 = window1/sum(sum(window1));

    window2 = fspecial('gaussian',7,7*k/6);
    window2 = window2/sum(sum(window2));

    %Compute DOG of the sigmaMap
    DoGSigma= filter2((window1-window2), sigmaMap, 'same');

    %Compute the sigmaMap of DoGSigma
    DoGSigma1 = computeSigmaMap(DoGSigma);

    %DNT of DOG of sigma.
    divNormDoGSigma = divisiveNormalization(DoGSigma);

    %shape, variance, skewness.
    [~,lsigma] = estimateggdparam(divNormDoGSigma(:));

    % Apply DNT on DoGSigma1
    divNormDoGSigma1 = divisiveNormalization(DoGSigma1);

    feat = [feat lsigma kurtosis(divNormDoGSigma(:)) kurtosis(divNormDoGSigma1(:))];

    % 71:74 --> 69:72
    % Laplacian of the Luminance Map
    addpath(genpath(fullfile('include', 'matlabPyrTools')));

    % Extract the first laplacian of the image.
    [pyr pind] = buildLpyr(imgray,4);
    res =  pyrBand(pyr, pind, 1);

    % Fit a GGD to the laplacian image.
    [alpha, sigma] = estimateggdparam(res(:));

    feat = [feat alpha sigma];
%     feat = feat([5,8,27,32,58,62,64,66,68,69,70]);
end

%% compute_friquee_chroma_features
function feat = compute_friquee_chroma_features(imyuv)
%% # 4  7(std) 15   24   26   47   60   62   72   73   79
    rgb = ycbcr2rgb(uint8(imyuv));
%     lab=rgb2lab(rgb);
    colorTransform = makecform('srgb2lab');
    lab = applycform(rgb, colorTransform);

    % Get the A and B components of the LAB color space.
    A = double(lab(:,:,2));
    B = double(lab(:,:,3));

    % Compute the chroma map.
    imChroma = sqrt(A.*A + B.*B);

    % Initializations
    scalenum=2;
    imChroma1=imChroma;

    % The array that aggregates all the features from different feature
    % maps.
    feat = [];

    %% scale1: 4 7 15 24
    % Applying divisive normalization operation on the given gray scale image.
    [structdis,~] = divisiveNormalization(imChroma);

    % The four neighborhood maps that would be constructed [1]
    shifts = [ 0 1;1 0 ; 1 1; -1 1];
    for itr_shift =1:4
        % Construct the product neighborhood map.
        shifted_structdis = circshift(structdis,shifts(itr_shift,:));
        pair = structdis(:).*shifted_structdis(:); % Element wise product.

        if itr_shift ==1
        % Fit an AGGD and extract its parameters.
        [~, ~, rightstd] = estimateaggdparam(pair);
        feat = [feat rightstd^2];%4
        elseif itr_shift ==2
        % Fit an AGGD and extract its parameters.
        [alpha, ~, ~] = estimateaggdparam(pair);
        feat = [feat alpha];%7
        elseif itr_shift ==3
        [~, leftstd, ~] = estimateaggdparam(pair);
        feat = [feat  leftstd^2];%15
        elseif itr_shift ==4
          feat = [feat kurtosis(pair(:))];%24
        end
    end
    sigmaMap = computeSigmaMap(imChroma);
    feat = [feat skewness(sigmaMap(:))];%26
    %% scale2:   47   60   62 - 35 = 12    25    27
    imChroma = imresize(imChroma,0.5);

    % Applying divisive normalization operation on the given gray scale image.
    [structdis,~] = divisiveNormalization(imChroma);

    shifts = [1 0];

    % Construct the product neighborhood map.
    shifted_structdis = circshift(structdis,shifts);
    pair = structdis(:).*shifted_structdis(:); % Element wise product.

    feat = [feat kurtosis(pair(:))]; %12

    sigmaMap = computeSigmaMap(imChroma);
    feat = [feat kurtosis(sigmaMap(:)) mean2(sigmaMap)]; %25 27
    % Some features are computedcccccccc at multiple scales %% 1 scale 35 features x2 =
    % 70

    imChroma=imChroma1;

    %Get the sigma map of the image.
    sigmaMap = computeSigmaMap(imChroma);
        % Construct two Gaussian windows where sig2 = 1.5*sig1;
    k=1.5;

    window1 = fspecial('gaussian',7,7/6);
    window1 = window1/sum(sum(window1));

    window2 = fspecial('gaussian',7,7*k/6);
    window2 = window2/sum(sum(window2));
     %Compute DOG of the sigmaMap
    DoGSigma= filter2((window1-window2), sigmaMap, 'same');
    %DNT of DOG of sigma.
    divNormDoGSigma = divisiveNormalization(DoGSigma);
    %shape, variance, skewness.
    [~,lsigma] = estimateggdparam(divNormDoGSigma(:));
    feat = [feat lsigma skewness(divNormDoGSigma(:))];% 72 73

    % Laplacian of the Chroma Map
    addpath(genpath(fullfile('include', 'matlabPyrTools')));

    [pyr pind] = buildLpyr(imChroma,4);
    res =  pyrBand(pyr, pind, 1);
    feat = [feat skewness(res(:))];%79
end

%% compute_friquee_lms_features
function feat = compute_friquee_lms_features(imyuv)
    % compute FRIQUEE LMS features # 52
    rgb = ycbcr2rgb(uint8(imyuv));
    lms = (colorspace('RGB->CAT02 LMS',rgb));
     % Get the M and S components from the LMS color space.
%     LM = double(lms(:,:,2));
    LS = double(lms(:,:,3));
    %FRIQUEE features from the LMS color space.
%     fM = friqueeMS(LM);%1:32
%     fS = friqueeMS(LS);%33:64
    LS = imresize(LS,0.5);
    sigmaMap = computeSigmaMap(LS);
    dnSigmaFeat = debiasedNormalizedFeats(sigmaMap);%8:11 19:20:21
%     featOpp = lmsColorOpponentFeats(lms); % Color opponency features.
    feat = dnSigmaFeat(2);
end

%% compute_tlvqm_lcf_features
function features = compute_tlvqm_lcf_features(this_fr, prev_fr, next_fr)
 % compute TLVQM LCF features # 1 5(both) 12 16 18 19(std) 22(std)
    [height,width,~] = size(this_fr);
    % Initialize parameters
    bl_size = floor(width/40);
    src_win = floor(width/40);

    % The following computations are done with reduced resolution
    this_fr = imresize(this_fr,0.5);
    prev_fr = imresize(prev_fr,0.5);
    next_fr = imresize(next_fr,0.5);
    [height,width,~] = size(this_fr);
    this_Y = this_fr(:,:,1);
    prev_Y = prev_fr(:,:,1);
    next_Y = next_fr(:,:,1);
    this_fr = ycbcr2rgb(this_fr);

    % Apply Sobel filter to the frames
    H = [-1 -2 -1; 0 0 0; 1 2 1]./8;
    sob_h_this = imfilter(this_Y,H');
    sob_v_this = imfilter(this_Y,H);

    % Reset edge pixels in the Sobeled frames
    sob_h_this(1:4,1:width)=0;
    sob_h_this(height-3:height,1:width)=0;
    sob_h_this(1:height,1:4)=0;
    sob_h_this(1:height,width-3:width)=0;
    sob_v_this(1:4,1:width)=0;
    sob_v_this(height-3:height,1:width)=0;
    sob_v_this(1:height,1:4)=0;
    sob_v_this(1:height,width-3:width)=0;

    sob_tot = sqrt(sob_v_this.^2+sob_h_this.^2);
    sob_h_prev = imfilter(prev_Y,H');
    sob_v_prev = imfilter(prev_Y,H);
    sob_h_next = imfilter(next_Y,H');
    sob_v_next = imfilter(next_Y,H);

    H1 = [1 1 1 1 1;1 1 1 1 1;-2 -2 0 1 1;-2 -2 -2 1 1;-2 -2 -2 1 1]./32;
    H2 = [-2 -2 -2 1 1;-2 -2 -2 1 1;-2 -2 0 1 1;1 1 1 1 1;1 1 1 1 1]./32;
    H3 = [1 1 -2 -2 -2;1 1 -2 -2 -2;1 1 0 -2 -2;1 1 1 1 1;1 1 1 1 1]./32;
    H4 = [1 1 1 1 1;1 1 1 1 1;1 1 0 -2 -2;1 1 -2 -2 -2;1 1 -2 -2 -2]./32;

    corner_avg(:,:,1) = abs(imfilter(this_Y, H1));
    corner_avg(:,:,2) = abs(imfilter(this_Y, H2));
    corner_avg(:,:,3) = abs(imfilter(this_Y, H3));
    corner_avg(:,:,4) = abs(imfilter(this_Y, H4));
    corner_max = max(corner_avg,[],3);
    corner_this = corner_max-min(corner_avg,[],3);

    mot_threshold = 0.01;

%     cor_max = sort(corner_max(:),'ascend');
%     glob_blockiness = 0;
%     if std2(cor_max(1:floor(0.99*end)))>0
%         glob_blockiness = 0.5*((mean(cor_max(1:floor(0.99*end)))/ ...
%                           std2(cor_max(1:floor(0.99*end))))^2);
%     end

    % Reset edge pixels in the corner point filtered frame
    corner_this(1:src_win+3,1:width)=0;
    corner_this(height-src_win-2:height,1:width)=0;
    corner_this(1:height,1:src_win+3)=0;
    corner_this(1:height,width-src_win-2:width)=0;

    corner_this_copy = corner_this(:);
    key_pix = zeros((height-6)*(width-6),2);
    n_key_pix = 0;

    im_y_vec = mod(0:width*height, height)+1;
    im_x_vec = floor((0:width*height-1)/height)+1;
    sob_this_cp = corner_this_copy(corner_this_copy>mot_threshold);
    im_y_vec = im_y_vec(corner_this_copy>mot_threshold);
    im_x_vec = im_x_vec(corner_this_copy>mot_threshold);

    % In the following loop, find the key pixels
    [mx,idx] = max(sob_this_cp);
    if ~isempty(idx)
        while mx>mot_threshold
            i = im_y_vec(idx(1));
            j = im_x_vec(idx(1));

            n_key_pix = n_key_pix + 1;
            key_pix(n_key_pix,:) = [i j];

            idx_remove = find(im_y_vec>=i-floor(bl_size) & ...
                              im_y_vec<=i+floor(bl_size) & ...
                              im_x_vec>=j-floor(bl_size) & ...
                              im_x_vec<=j+floor(bl_size));
            sob_this_cp(idx_remove)=[];
            im_y_vec(idx_remove)=[];
            im_x_vec(idx_remove)=[];

            [mx,idx] = max(sob_this_cp);
        end
    end
    key_pix=key_pix(1:n_key_pix,:);

    non_mot_area = ones(height, width);

    num_mot_points = 0;
    max_mot_points = (height/bl_size)*(width/bl_size);

    %tic
    % In the following loop, find the motion vectors for each key pixel
    motion_vec = [];

    distance_matrix = ones(2*src_win+1);
    for i=1:2*src_win+1
        for j=1:2*src_win+1
            distance_matrix(i,j) = ...
                sqrt((1+src_win-i).^2+(1+src_win-j).^2)/sqrt(2*src_win^2);
        end
    end
    distances = distance_matrix(:);

    uncertain = 0;

    % Loop through the key pixels
    for z = 1:n_key_pix

        tar_y = key_pix(z,1);
        tar_x = key_pix(z,2);
        match_y_bw = tar_y;
        match_x_bw = tar_x;
        match_y_fw = tar_y;
        match_x_fw = tar_x;

        surr_win_v_prev = sob_v_prev(tar_y-src_win-2:tar_y+src_win+2, ...
                                     tar_x-src_win-2:tar_x+src_win+2);
        surr_win_h_prev = sob_h_prev(tar_y-src_win-2:tar_y+src_win+2, ...
                                     tar_x-src_win-2:tar_x+src_win+2);
        diff_win_prev = (sob_v_this(tar_y, tar_x)-surr_win_v_prev).^2 + ...
                        (sob_h_this(tar_y, tar_x)-surr_win_h_prev).^2;

        surr_win_v_next = sob_v_next(tar_y-src_win-2:tar_y+src_win+2, ...
                                     tar_x-src_win-2:tar_x+src_win+2);
        surr_win_h_next = sob_h_next(tar_y-src_win-2:tar_y+src_win+2, ...
                                     tar_x-src_win-2:tar_x+src_win+2);
        diff_win_next = (sob_v_this(tar_y, tar_x)-surr_win_v_next).^2 + ...
                        (sob_h_this(tar_y, tar_x)-surr_win_h_next).^2;

        for i=-1:1
            for j=-1:1
                if i~=0 || j~=0
                    diff_win_prev(3:end-2,3:end-2) = ...
                        diff_win_prev(3:end-2,3:end-2) + ...
                        (sob_v_this(tar_y+i, tar_x+j)- ...
                          surr_win_v_prev(3+i:end-2+i,3+j:end-2+j)).^2+ ...
                        (sob_h_this(tar_y+i, tar_x+j)- ...
                          surr_win_h_prev(3+i:end-2+i,3+j:end-2+j)).^2;
                    diff_win_next(3:end-2,3:end-2) = ...
                        diff_win_next(3:end-2,3:end-2) + ...
                        (sob_v_this(tar_y+i, tar_x+j)- ...
                          surr_win_v_next(3+i:end-2+i,3+j:end-2+j)).^2+...
                        (sob_h_this(tar_y+i, tar_x+j)- ...
                        surr_win_h_next(3+i:end-2+i,3+j:end-2+j)).^2;
                end
            end
        end
        diff_win_prev = diff_win_prev(3:end-2,3:end-2);
        diff_win_next = diff_win_next(3:end-2,3:end-2);

        orig_diff_bw = diff_win_prev(1+src_win,1+src_win);
        orig_diff_fw = diff_win_next(1+src_win,1+src_win);

        diff_bw = diff_win_prev(1+src_win,1+src_win);
        if orig_diff_bw>0.005
            [sorted,idx] = sort(diff_win_prev(:),'ascend');
            min_diff = orig_diff_bw;
            if length(sorted)>=2
                if sorted(1)<=0.8*sorted(2) || ...
                   distances(idx(1))<distances(idx(2))
                    min_diff = sorted(1);
                else
                    [idx,~] = find(0.8.*diff_win_prev(:)<=sorted(1));
                    [~,idx2] = sort(distances(idx),'ascend');
                    if diff_win_next(idx(idx2(1)))<1.1*sorted(1)
                        min_diff = diff_win_prev(idx(idx2(1)));
                    elseif sorted(1)<diff_bw*0.9
                        min_diff = sorted(1);
                    end
                    uncertain = uncertain + 1;
                end
                if min_diff*1.01<orig_diff_bw
                    [y,x] = find(diff_win_prev==min_diff);
                    match_y_bw = tar_y+y(1)-src_win-1;
                    match_x_bw = tar_x+x(1)-src_win-1;
                    diff_bw = diff_win_prev(y(1),x(1));
                end
            end
        end

        diff_fw = diff_win_next(1+src_win,1+src_win);
        if orig_diff_fw>0.005
            [sorted,idx] = sort(diff_win_next(:),'ascend');
            min_diff = orig_diff_fw;
            if length(sorted)>=2
                if sorted(1)<0.8*sorted(2) || ...
                   distances(idx(1))<distances(idx(2))
                    min_diff = sorted(1);
                else
                    [idx,~] = find(0.8.*diff_win_next(:)<=sorted(1));
                    [~,idx2] = sort(distances(idx),'ascend');
                    if diff_win_next(idx(idx2(1)))<1.1*sorted(1)
                        min_diff = diff_win_next(idx(idx2(1)));
                    elseif sorted(1)<diff_fw*0.9
                        min_diff = sorted(1);
                    end
                    uncertain = uncertain + 1;
                end
                if min_diff*1.01<orig_diff_fw
                    [y,x] = find(diff_win_next==min_diff);
                    match_y_fw = tar_y+y(1)-src_win-1;
                    match_x_fw = tar_x+x(1)-src_win-1;
                    diff_fw = diff_win_next(y(1),x(1));
                end
            end
        end

        % Add motion vector to the list of motion vectors
        if (orig_diff_bw > diff_bw*1.01 && ...
                (tar_y ~= match_y_bw || tar_x ~= match_x_bw)) || ...
           (orig_diff_fw > diff_fw*1.01 && ...
                (tar_y ~= match_y_fw || tar_x ~= match_x_fw))

            non_mot_area(max(1,tar_y-bl_size):min(height,tar_y+bl_size),...
                max(1,tar_x-bl_size):min(width,tar_x+bl_size))=0;
            non_mot_area(max(1,match_y_bw-bl_size): ...
                min(height,match_y_bw+bl_size),...
                max(1,match_x_bw-bl_size):...
                min(width,match_x_bw+bl_size)) = 0;
            non_mot_area(max(1,match_y_fw-bl_size):...
                min(height,match_y_fw+bl_size),...
                max(1,match_x_fw-bl_size):...
                min(width,match_x_fw+bl_size)) = 0;
        end

        num_mot_points = num_mot_points + 1;
        motion_vec = [motion_vec; ...
                      tar_y-match_y_bw tar_x-match_x_bw ...
                      match_y_fw-tar_y match_x_fw-tar_x ...
                      tar_y tar_x ...
                      orig_diff_bw diff_bw ...
                      orig_diff_fw diff_fw];
    end
    %toc
 % Compute motion point related statistics
%     motion_uncertainty = 0.5*uncertain/max_mot_points;
%     motion_density = 0;
    motion_intensity = 0;
%     std_mot_intensity = 0;
%     avg_mot_pos = 0;
%     avg_mot_sprd = 0;
%     mot_pred_acc = 0;
    mot_y = 0.5;
%     mot_x = 0.5;
    jerkiness = 0;
%     jerk_cons = 0;
    motion_vec_bg = [];
    num_bg_mot_points = 0;
    if num_mot_points>0
%         motion_density = num_mot_points/(width*height/bl_size^2);
        mot_intensity_vec = sqrt(((motion_vec(:,1)./src_win).^2 + ...
                                  (motion_vec(:,2)./src_win).^2 + ...
                                  (motion_vec(:,3)./src_win).^2 + ...
                                  (motion_vec(:,4)./src_win).^2)./4.0);
        sum_mot_int = sum(mot_intensity_vec);
        motion_intensity = (sum(mot_intensity_vec)/max_mot_points)^0.25;
%         std_mot_intensity = std(mot_intensity_vec);

        if sum_mot_int>0
            % Compute motion position in relation with the screen midpoint
%             avg_motp_y = sum(mot_intensity_vec.*motion_vec(:,5))/...
%                            sum_mot_int;
%             std_motp_y = sqrt(sum(mot_intensity_vec.*...
%                            (motion_vec(:,5)-avg_motp_y).^2)/sum_mot_int);
%             avg_mot_pos_y = (avg_motp_y-height/2)/(height/2);
%             sprd_mot_pos_y = std_motp_y/height;
%             avg_motp_x = sum(mot_intensity_vec.*motion_vec(:,6))/...
%                            sum_mot_int;
%             std_motp_x = sqrt(sum(mot_intensity_vec.*...
%                            (motion_vec(:,6)-avg_motp_x).^2)/sum_mot_int);
%             avg_mot_pos_x = (avg_motp_x-width/2)/(width/2);
%             sprd_mot_pos_x = std_motp_x/width;
%
%             avg_mot_pos = sqrt(avg_mot_pos_y^2+avg_mot_pos_x^2);
%             avg_mot_sprd = sqrt(sprd_mot_pos_y^2+sprd_mot_pos_x^2);

            % Mean motion along x and y axis
            mot_y = mean(0.25.*(motion_vec(:,1)+motion_vec(:,3))./ ...
                      src_win+0.5);
%             mot_x = mean(0.25.*(motion_vec(:,2)+motion_vec(:,4))./ ...
%                       src_win+0.5);

            % Average motion prediction improvement
%             mot_pred_acc_bw = mean(motion_vec(:,7)-motion_vec(:,8));
%             mot_pred_acc_fw = mean(motion_vec(:,9)-motion_vec(:,10));
%             mot_pred_acc = 0.5*(mot_pred_acc_bw+mot_pred_acc_fw).^0.5;

            % Motion jerkiness
            mot_y_diff = 0.5.*(motion_vec(:,1)'-motion_vec(:,3)')./src_win;
            mot_x_diff = 0.5.*(motion_vec(:,2)'-motion_vec(:,4)')./src_win;
            mot_diff = sqrt(mot_y_diff.^2+mot_x_diff.^2);
            jerkiness = mean(mot_diff.^0.5);
%             jerk_cons = std(mot_diff.^0.5);
        end

        avg_mot_x = mean(0.5.*motion_vec(:,2)+0.5.*motion_vec(:,4));
        avg_mot_y = mean(0.5.*motion_vec(:,1)+0.5.*motion_vec(:,3));
        std_mot_x = std(0.5.*motion_vec(:,2)+0.5.*motion_vec(:,4));
        std_mot_y = std(0.5.*motion_vec(:,1)+0.5.*motion_vec(:,3));

        for z=1:num_mot_points
            mot_x_this = 0.5*motion_vec(z,2)+0.5*motion_vec(z,4);
            mot_y_this = 0.5*motion_vec(z,1)+0.5*motion_vec(z,3);
            if mot_x_this > avg_mot_x-std_mot_x && ...
               mot_x_this < avg_mot_x+std_mot_x && ...
               mot_y_this > avg_mot_y-std_mot_y && ...
               mot_y_this < avg_mot_y+std_mot_y

                num_bg_mot_points = num_bg_mot_points + 1;
                motion_vec_bg = [motion_vec_bg; motion_vec(z,:)];
            end
        end
    end
 % Compute motion point related statistics
%     egomotion_density = 0;
%     egomotion_intensity = 0;
    std_egomot_intensity = 0;
%     avg_egomot_pos = 0;
%     avg_egomot_sprd = 0;
%     egomot_pred_acc = 0;
%     mot_y_bg = 0.5;
%     mot_x_bg = 0.5;
    if num_bg_mot_points>0
%         egomotion_density = num_bg_mot_points/(width*height/bl_size^2);
        bg_mot_intensity_vec = sqrt(((motion_vec_bg(:,1)./src_win).^2 + ...
                                     (motion_vec_bg(:,2)./src_win).^2 + ...
                                     (motion_vec_bg(:,3)./src_win).^2 + ...
                                     (motion_vec_bg(:,4)./src_win).^2)  ...
                                      ./4.0);
%         sum_bg_mot_int = sum(bg_mot_intensity_vec);
%         egomotion_intensity = (sum(bg_mot_intensity_vec)/...
%                                 max_mot_points)^0.25;
        std_egomot_intensity = std(bg_mot_intensity_vec);

        % Compute motion position in relation with the screen midpoint
%         if sum_bg_mot_int>0
%             avg_motp_y = sum(bg_mot_intensity_vec.*motion_vec_bg(:,5))/...
%                            sum_bg_mot_int;
%             std_motp_y = sqrt(sum(bg_mot_intensity_vec.*...
%                            (motion_vec_bg(:,5)-avg_motp_y).^2)/...
%                               sum_bg_mot_int);
%             avg_mot_pos_y = (avg_motp_y-height/2)/(height/2);
%             sprd_mot_pos_y = std_motp_y/height;
%             avg_motp_x = sum(bg_mot_intensity_vec.*motion_vec_bg(:,6))/...
%                            sum_bg_mot_int;
%             std_motp_x = sqrt(sum(bg_mot_intensity_vec.*...
%                            (motion_vec_bg(:,6)-avg_motp_x).^2)/...
%                            sum_bg_mot_int);
%             avg_mot_pos_x = (avg_motp_x-width/2)/(width/2);
%             sprd_mot_pos_x = std_motp_x/width;
%
%             avg_egomot_pos = sqrt(avg_mot_pos_y^2+avg_mot_pos_x^2);
%             avg_egomot_sprd = sqrt(sprd_mot_pos_y^2+sprd_mot_pos_x^2);

            % Average egomotion prediction improvement
%             mot_pred_acc_bw = mean(motion_vec_bg(:,7)-motion_vec_bg(:,8));
%             mot_pred_acc_fw = mean(motion_vec_bg(:,9)-motion_vec_bg(:,10));
%             egomot_pred_acc = 0.5*(mot_pred_acc_bw+mot_pred_acc_fw).^0.5;

%             mot_y_bg = mean(0.25.*(motion_vec_bg(:,1)+...
%                                    motion_vec_bg(:,3))./src_win+0.5);
%             mot_x_bg = mean(0.25.*(motion_vec_bg(:,2)+...
%                                    motion_vec_bg(:,4))./src_win+0.5);
%         end
    end

%     mot_size = sum(sum(1-non_mot_area));
    non_mot_size = sum(sum(non_mot_area));
%  static_area_flicker = 0;
    static_area_flicker_std = 0;
    if non_mot_size>0
        % Sum of the pixel differences in the static area
        static_area_flicker_bw = sum(non_mot_area(:) .* ...
                                 abs(this_Y(:)-prev_Y(:)))/non_mot_size;
        static_area_flicker_fw = sum(non_mot_area(:) .* ...
                                 abs(this_Y(:)-next_Y(:)))/non_mot_size;
        static_area_flicker = 0.5*(static_area_flicker_bw + ...
                                   static_area_flicker_fw);
        % Variance of pixel differences in the static area
        st_diff_bw = abs(this_Y(:)-prev_Y(:));
        st_diff_fw = abs(this_Y(:)-next_Y(:));
        static_area_flicker_std = sum(non_mot_area(:)' .* ...
                                  abs(max([st_diff_bw'; st_diff_fw']) - ...
                                  static_area_flicker))/non_mot_size;
    end

    % Spatial activity in the static area
    si = std2(sob_tot).^0.25;

    %[blur glob_blockiness si]

    % Temporal activity standard deviation in the static area
    ti_prev = mean(abs(this_Y(:)-prev_Y(:)));
    ti_next = mean(abs(this_Y(:)-next_Y(:)));
    ti_mean = mean([ti_prev ti_next]).^0.25;

    % Normalize static area size
%     mot_size = mot_size / (width*height);
    features = [motion_intensity,  ...
                std_egomot_intensity, ...
                si, ...
                jerkiness, ...
                ti_mean, ...
                mot_y, ...
                static_area_flicker_std];
end

%% compute_tlvqm_hcf_features
function features = compute_tlvqm_hcf_features(imyuv)
% compute TLVQM HCF features # 1 2    17    22    24    30
% Initializations
    mono_image = imyuv(:,:,1);
    image = ycbcr2rgb(imyuv);
    lab_image = rgb2lab(image);
    [height,width,depth] = size(image);

    % Make Sobeled image
    mask = zeros(height,width);
    mask(2:end-1,2:end-1)=1;
    H = [1 2 1; 0 0 0; -1 -2 -1]./8;

    % Make Sobeled image in CIELAB color space
    sob_image_lab_x = (imfilter(lab_image(:,:,1)./100.0,H).^2 + ...
                       imfilter(lab_image(:,:,2)./50.0,H).^2 + ...
                       imfilter(lab_image(:,:,3)./50.0,H).^2).*mask;
    sob_image_lab_y = (imfilter(lab_image(:,:,1)./100.0,H').^2 + ...
                       imfilter(lab_image(:,:,2)./50.0,H').^2 + ...
                       imfilter(lab_image(:,:,3)./50.0,H').^2).*mask;
    sob_image = sqrt(sob_image_lab_x+sob_image_lab_y);

    % Compute fetures for different feature groups
%     [a,b,sat_image1] = compute_saturation(mono_image,1);
% %     sat_bright = [a];
%     [a,b,sat_image2] = compute_saturation(mono_image,0);
% %     sat_dark = [a b];
%     sat_image = max(sat_image1, sat_image2);
%     saturation_ftr = [sat_bright sat_dark]; % 5:8
%     saturation_ftr = sat_bright;
    spatial_ftr = spatial_activity_features(sob_image, []);% 1:4
% %     noisiness_ftr = noise_features(mono_image, sat_image, lab_image); % 9:11
%     blockiness_ftr = blockiness_features(sob_image_lab_x.^0.5, ...
%                                          sob_image_lab_y.^0.5); % 12:14
    contrast_color_ftr = contrast_chroma_features(lab_image, []); % 15:18
%     dct_ftr = dct_features(mono_image); % 19:21
    sharpness_ftr = sharpness_features(sob_image); % 22:30

    % Make the HC feature vector
    features = [spatial_ftr   ...
                contrast_color_ftr   ...
                sharpness_ftr];
end

% This function computes the saturation (bright or dark)
function [len,num,segs] = compute_saturation(image, is_bright)

    [height,width] = size(image);

    lens = [];
    num = 0;

    segs = zeros(height,width);

    if (is_bright==1 && max(max(image))>0.9) || ...
       (is_bright==0 && min(min(image))<0.1)

        segs = seg_loop(image,segs,3,3,0.05, is_bright);
        for i=1:max(max(segs))
            len = length(find(segs==i));
            if len<50
                segs(find(segs==i))=0;
            else
                lens = [lens len];
                num = num + 1;
            end
        end
        segs(find(segs>0))=1;
    end

    len = sum(lens)/(width*height);
    if num > 0
        num = len / num;
    end

end

% This function is used for segmentation by measure_saturation
function segim = seg_loop(image, segs, wh, ww, interval, is_bright)

    [height,width] = size(image);

    segim = segs;

    maxi = max(max(image));
    mini = min(min(image));

    for i=1:height-wh+1
        for j=1:width-ww+1
            if (is_bright == 1 && ...
              min(min(image(i:i+wh-1,j:j+ww-1)))>maxi-interval) || ...
              (is_bright == 0 && ...
              max(max(image(i:i+wh-1,j:j+ww-1)))<mini+interval)

                maxsg = max(max(segim(i:i+wh-1,j:j+ww-1)));
                if maxsg>0
                    segs_temp = reshape(segim(i:i+wh-1,j:j+ww-1),wh*ww,1);
                    minsg=min(segs_temp(find(segs_temp>0)));
                    segim(i:i+wh-1,j:j+ww-1)=minsg;
                    if minsg<maxsg
                        segim(find(segim==maxsg))=minsg;
                    end
                else
                    segim(i:i+wh-1,j:j+ww-1)=max(max(segim(:,:)))+1;
                end
            end
        end
    end

end

% This function is used to compute noise related features
function out = noise_features(mono_image, sat_im, lab_image)

    [height,width] = size(mono_image);

    new_im = zeros(height, width, 3);

    nonsat_pix = 0;
    noise_pix = 0;
    noise_int = [];

    % Loop through pixels to find noise pixels
    for i=5:height-4
        for j=5:width-4
            if sat_im(i,j)==0
                surr_pix = mono_image(i-2:i+2,j-2:j+2);
                surr_pix = surr_pix(:);
                surr_pix = [surr_pix(1:12); surr_pix(14:25)];
                if (mono_image(i,j)>max(surr_pix) || ...
                    mono_image(i,j)<min(surr_pix))
                    surr_pix = mono_image(i-4:i+4,j-4:j+4);
                    if std(surr_pix)<0.05
                        new_im(i,j,2) = 1;
                        pix_diff = sqrt( ...
                            (mean(lab_image(i-3:i+3,j-3:j+3,1))-...
                                 lab_image(i,j,1)).^2 + ...
                            (mean(lab_image(i-3:i+3,j-3:j+3,2))-...
                                 lab_image(i,j,2)).^2 + ...
                            (mean(lab_image(i-3:i+3,j-3:j+3,3))-...
                                  lab_image(i,j,3)).^2);
                        noise_int = [noise_int pix_diff/100];
                        noise_pix = noise_pix + 1;
                    end
                end
                nonsat_pix = nonsat_pix + 1;
            end
        end
    end

    a = 0;
    b = 0;
    c = 0;

    if nonsat_pix > 0 && noise_pix > 0
        % noise density
        a = noise_pix / nonsat_pix;
        b = mean(noise_int);
        c = std(noise_int);
    end

    out = [a b c];

end

% This function is used to compute spatial activity features
function out = spatial_activity_features(sobel_image, ~)

%     [height,width] = size(sobel_image);
%
%     sob_dists = zeros(1,height*width);
%     sob_dists2 = zeros(height*width,2);
%     sob_str = zeros(1,height*width);
%     sumstr = 0;
%
%     n = 0;
%     for i=1:height
%         for j=1:width
%             if sat_image(i,j)==0
%                 if sobel_image(i,j)<0.01
%                     sobel_image(i,j)=0;
%                 end
%                 sumstr = sumstr + sobel_image(i,j);
%                 if sobel_image(i,j) > 0
%                     n = n + 1;
%                     sob_str(n) = sobel_image(i,j);
%                     sob_dists(n) = sqrt((i/height-0.5)^2+(j/width-0.5)^2);
%                     sob_dists2(n,1) = i/height-0.5;
%                     sob_dists2(n,2) = j/width-0.5;
%                 end
%             end
%         end
%     end
%
%     sob_str = sob_str(1:n);
%     sob_dists = sob_dists(1:n);
%     sob_dists2 = sob_dists2(1:n,:);

%     a = 0;
%     b = 0;
%     c = 0;
%     d = 0;

%     if ~isempty(sob_str)>0
        a = mean(mean(sobel_image));
        b = std2(sobel_image);
%         d = w_std(sob_dists, sob_str);
%         mean_y = sum(sob_str'.*sob_dists2(:,1))/sum(sob_str);
%         mean_x = sum(sob_str'.*sob_dists2(:,2))/sum(sob_str);
%         c = sqrt(mean_y^2+mean_x^2);
%     end

%     out = [a b c d];
    out = [ a b];

end

% Function for "weighted standard deviation", used by function
% measure_spatial_activity
function res = w_std(input, weights)

    wg_n = sum(weights);
    wg_input = input.*weights;
    wg_mean = mean(input.*weights);

    res = sqrt(sum((wg_input-wg_mean).^2)/wg_n);
end

% This function is used to compute blockiness index
function blockiness = blockiness_features(sob_y, sob_x)

    [height,width] = size(sob_y);

    hor_tot = zeros(1,height-4);
    ver_tot = zeros(1,width-4);

    for i=3:height-2
        hor_tot(i)=mean(sob_y(i,:)-sob_x(i,:));
    end
    for j=3:width-2
        ver_tot(j)=mean(sob_x(:,j)-sob_y(:,j));
    end

    % compute autocorrelations
    autocr_hor = zeros(1,23);
    autocr_ver = zeros(1,23);
    for i=0:22
        autocr_hor(i+1) = sum(hor_tot(1:end-i).*hor_tot(1+i:end));
        autocr_ver(i+1) = sum(ver_tot(1:end-i).*ver_tot(1+i:end));
    end

    % Find the highest local maximum (other than 0)
    localpeaks = 0;
%     peakdist = 0;
    max_hor = 0;
    max_ver = 0;
    min_hor = autocr_hor(1);
    min_ver = autocr_ver(1);
    max_hor_diff = 0;
    max_ver_diff = 0;
    for i=2:22
        if autocr_hor(i)>max(autocr_hor(i-1),autocr_hor(i+1))
            localpeaks = localpeaks+1/42;
        end
        if autocr_hor(i)<min(autocr_hor(i-1),autocr_hor(i+1)) && ...
                autocr_hor(i)<min_hor
            min_hor = autocr_hor(i);
        elseif autocr_hor(i)>max(autocr_hor(i-1),autocr_hor(i+1)) && ...
                autocr_hor(i)-min_hor>max_hor_diff
            max_hor = autocr_hor(i);
            max_hor_diff = max_hor-min_hor;
%             peakdist = (i-1)/21;
        end
        if autocr_ver(i)>max(autocr_ver(i-1),autocr_ver(i+1))
            localpeaks = localpeaks + 1/42;
        end
        if autocr_ver(i)<min(autocr_ver(i-1),autocr_ver(i+1)) && ...
                autocr_ver(i)<min_ver
            min_ver = autocr_ver(i);
        elseif autocr_ver(i)>max(autocr_ver(i-1),autocr_ver(i+1)) && ...
                autocr_ver(i)-min_ver>max_ver_diff
            max_ver = autocr_ver(i);
            max_ver_diff = max_ver-min_ver;
%             peakdist = (i-1)/21;
        end
    end

    a = 0;
    if autocr_hor(1)>0 && autocr_ver(1)>0
        if max_hor>0 && max_ver>0
            a = max((max_hor_diff/autocr_hor(1)), ...
                             (max_ver_diff/autocr_ver(1)))^0.5;
        elseif max_hor>0
            a = (max_hor_diff/autocr_hor(1))^0.5;
        elseif max_ver>0
            a = (max_ver_diff/autocr_ver(1))^0.5;
        end
    end

%     b = peakdist;
%     c = localpeaks;
    blockiness = a;
end

% This function is used to compute contrast and chroma related features
function out = contrast_chroma_features(lab_image, ~)

%     a=0;
%     b=0;
    c=0;
%     d=0;

%     [height,width,depth] = size(lab_image);
%     yuv_int = floor(lab_image(:,:,1));

    %sat_image = sat_image(:);
%     yuv_int2 = yuv_int(sat_image(:)==0);
%     cumu_err = 0;
%     cumu_tar = 0;
%     if ~isempty(yuv_int2)
%         for i=0:100
%             cumu_tar = cumu_tar + 1/100;
%             cumu_err = cumu_err + (sum(yuv_int2<=i)/length(yuv_int2) - ...
%                                    cumu_tar)/100;
%         end
% %         a = (cumu_err+1.0)/2.0;
%         b = 0.5*(1-cumu_err);
%     else
% %         a = 1;
%         b = sum(sum(lab_image(:,:,1)))/50;
%     end
    c = sqrt(mean(mean((lab_image(:,:,2)./50).^2 + ...
         (lab_image(:,:,3)./50).^2)));
%     d = 0;
%     if std2(lab_image(:,:,1))>0
%         d = 0.01*(std2(lab_image(:,:,2))+std2(lab_image(:,:,3)));
%     end

    out = [c];
end

% This function is used to compute dct derived features
function out = dct_features(im)

    % Input is monochrome image
    [height,width] = size(im);

    out_im = abs(dct2(im)).^.5;

    area1 = imresize(out_im(1:floor(height/2),1:floor(width/2)),0.25);
    area2 = imresize(out_im(1:floor(height/2),...
                            width:-1:width-floor(width/2)+1),0.25);
    area3 = imresize(out_im(height:-1:height-floor(height/2)+1,...
                            1:floor(width/2)),0.25);
    area4 = imresize(out_im(height:-1:height-floor(height/2)+1,...
                            width:-1:width-floor(width/2)+1),0.25);
    a = max(0,max(corr(area1(:),area2(:)),corr(area1(:),area3(:))));
%     b = 0;
%     if mean(area1)>0
%         b = mean(area4)/mean(area1);
%     end
%     c = 0;
%     if max(mean(area2),mean(area3))>0
%         c = min(mean(area2),mean(area3))/max(mean(area2),mean(area3));
%     end

    out = [a];

end

% This function is used to compute sharpness related features
function out = sharpness_features(im)

    [~, width] = size(im);

    % Full HD video could be downsized
    if width>1280
        im = imresize(im,0.5);
    end

    H = [-1 -2 -1; 0 0 0; 1 2 1]./8;
    im_s_h = imfilter(im,H');
    im_s_v = imfilter(im,H);
    im_s = sqrt(im_s_h.^2+im_s_v.^2);

    [height,width] = size(im_s_h);
    bl_size = 16;
    conv_list = [];

    blur_im = zeros(height,width);
    edge_strong = [];
    edge_all = [];

    conv_cube = [];
    blurvals = [];

    n_blks = 0;

    conv_val_tot = zeros(17);
    for y=floor(bl_size/2):bl_size:height-ceil(3*bl_size/2)
        for x=floor(bl_size/2):bl_size:width-ceil(3*bl_size/2)

            n_blks = n_blks + 1;

            conv_val = zeros(17);
            for i=0:6
                for j=0:6
                    if i==0 || j==0 || i==j
                        weight_h = 1;
                        weight_v = 1;
                        if i~=0 || j~=0
                            weight_h = abs(i)/(abs(i)+abs(j));
                            weight_v = abs(j)/(abs(i)+abs(j));
                        end
                        diff_h = (im_s_h(y+i:y+bl_size+i,   ...
                                         x+j:x+bl_size+j).* ...
                                  im_s_h(y:y+bl_size,       ...
                                         x:x+bl_size));
                        diff_v = (im_s_v(y+i:y+bl_size+i,   ...
                                         x+j:x+bl_size+j).* ...
                                  im_s_v(y:y+bl_size,       ...
                                         x:x+bl_size));
                        conv_val(i+9,j+9) = weight_h*(mean(diff_h(:)))+ ...
                                            weight_v*(mean(diff_v(:)));
                    end
                end
            end
            blur_im(y:y+bl_size-1,x:x+bl_size-1)=0.5;
            edge_all =  [edge_all conv_val(9,9)];
            if conv_val(9,9)>0.0001
                edge_strong =  [edge_strong conv_val(9,9)];
                conv_val=conv_val./conv_val(9,9);
                conv_val_tot = conv_val_tot + conv_val;

                new_conv_v = [];
                for i=1:6
                    new_conv_v = [new_conv_v sum(sum(conv_val(9-i:9+i,...
                                                            9-i:9+i)))- ...
                                             sum(sum(conv_val(10-i:8+i, ...
                                                            10-i:8+i)))];
                end
                if new_conv_v(1)>0
                    new_conv_v=new_conv_v./new_conv_v(1);
                end

                conv_list = [conv_list; new_conv_v];
                conv_cube(:,:,1)=conv_val;
                blurvals = [blurvals std2(im_s(y:y+bl_size, x:x+bl_size))];

                blur_im(y:y+bl_size-1,x:x+bl_size-1) = ...
                                  0.5 + mean(new_conv_v(2:6))/5;
            end
        end
    end

    % Find the sharpest blocks
    blurs_sharp = [];
    blurs_blur = [];
    if ~isempty(blurvals)
        for i=1:length(blurvals)
            if blurvals(i)>mean(blurvals)
                conv_val_tot = + conv_val_tot + conv_cube(:,:,1);
                blurs_sharp = [blurs_sharp blurvals(i)];
            else
                blurs_blur = [blurs_blur blurvals(i)];
            end
        end
    end

    n_sharps = length(blurs_sharp)/n_blks;
    n_blurs = length(blurs_blur)/n_blks;
    mean_sharps = 0;
    mean_blurs = 0;
    if ~isempty(blurs_sharp)
        mean_sharps = mean(blurs_sharp);
    end
    if ~isempty(blurs_blur)
        mean_blurs = mean(blurs_blur);
    end

%     if conv_val_tot(9,9)>0
%         conv_val_tot=conv_val_tot./conv_val_tot(9,9);
%     end

    new_conv_v = zeros(1,9);
    if ~isempty(edge_strong)>0
        if length(conv_list(:,1))>1
            new_conv_v = mean(conv_list);
        else
            new_conv_v = conv_list;
        end
    end

    % find local min and/or local max
%     localmin=0;
%     localmindist=0;
%     localmax=0;
%     localmaxdist=0;
%     for i=9:14
%         for j=9:14
%             if (i~=9 && j==9) || (j~=9 && i==9) || (i==j && i>9)
%                 conv_val_comp=conv_val_tot(i-1:i+1,j-1:j+1);
%                 conv_val_comp=conv_val_comp(:);
%                 if i==9
%                     conv_val_comp=conv_val_comp([4 6]);
%                 elseif j==9
%                     conv_val_comp=conv_val_comp([2 8]);
%                 else
%                     conv_val_comp=conv_val_comp([1 9]);
%                 end
%                 if conv_val_tot(i,j)>max(conv_val_comp) && ...
%                         conv_val_tot(i,j)>localmax
%                     localmax = conv_val_tot(i,j);
% %                     localmaxdist = sqrt((i-9)^2+(j-9)^2);
%                 elseif conv_val_tot(i,j)<min(conv_val_comp) && ...
%                         conv_val_tot(i,j)<localmin
%                     localmin = conv_val_tot(i,j);
% %                     localmindist = sqrt((i-9)^2+(j-9)^2);
%                 end
%             end
%         end
%     end

%     out = [mean(new_conv_v(2:6)) mean(new_conv_v(2:4)) new_conv_v(2) ...
%            localmaxdist/5 localmindist/5 ...
%            n_sharps n_blurs mean_sharps mean_blurs];

    out = [mean(new_conv_v(2:6)) new_conv_v(2) mean_blurs];
end

% Read one frame from YUV file
function YUV = YUVread(f,dim,frnum)

    % This function reads a frame #frnum (0..n-1) from YUV file into an
    % 3D array with Y, U and V components

    fseek(f,dim(1)*dim(2)*1.5*frnum,'bof');

    % Read Y-component
    Y=fread(f,dim(1)*dim(2),'uchar');
    if length(Y)<dim(1)*dim(2)
        YUV = [];
        return;
    end
    Y=cast(reshape(Y,dim(1),dim(2)),'double');

    % Read U-component
    U=fread(f,dim(1)*dim(2)/4,'uchar');
    if length(U)<dim(1)*dim(2)/4
        YUV = [];
        return;
    end
    U=cast(reshape(U,dim(1)/2,dim(2)/2),'double');
    U=imresize(U,2.0);

    % Read V-component
    V=fread(f,dim(1)*dim(2)/4,'uchar');
    if length(V)<dim(1)*dim(2)/4
        YUV = [];
        return;
    end
    V=cast(reshape(V,dim(1)/2,dim(2)/2),'double');
    V=imresize(V,2.0);

    % Combine Y, U, and V
    YUV(:,:,1)=Y';
    YUV(:,:,2)=U';
    YUV(:,:,3)=V';
end

author: vztu

repo: VIDEVAL

save_path: github-repos/MATLAB/vztu-VIDEVAL


text:
%CODEGENERATOR.GENSLBLOCKFKINE Generate Simulink block for forward kinematics
%
% cGen.genslblockfkine() generates a robot-specific Simulink block to compute
% forward kinematics.
%
% Notes::
% - Is called by CodeGenerator.genfkine if cGen has active flag genslblock.
% - The Simulink blocks are generated and stored in a robot specific block
%   library cGen.slib in the directory cGen.basepath.
% - Blocks are created for intermediate transforms T0, T1 etc. as well.
%
% Author::
%  Joern Malzahn, (joern.malzahn@tu-dortmund.de)
%
% See also CodeGenerator.CodeGenerator, CodeGenerator.genfkine.

% Copyright (C) 2012-2014, by Joern Malzahn
%
% This file is part of The Robotics Toolbox for Matlab (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com
%
% The code generation module emerged during the work on a project funded by
% the German Research Foundation (DFG, BE1569/7-1). The authors gratefully
% acknowledge the financial support.

function genslblockfkine(CGen)

%% Open or create block library
bdclose('all')                                                              % avoid problems with previously loaded libraries
load_system('simulink');
if ~(exist([CGen.slibpath,simulinkext]) == 2)                                  % Create new block library if none exists
 CGen.createnewblocklibrary;
end
open_system(CGen.slibpath);
set_param(CGen.slib,'lock','off');

q = CGen.rob.gencoords;

%% Forward kinematics up to tool center point
CGen.logmsg([datestr(now),'\tGenerating forward kinematics Simulink block up to the end-effector frame']);
symname = 'fkine';
fname = fullfile(CGen.sympath,[symname,'.mat']);

if exist(fname,'file')
    tmpStruct = load(fname);
else
    error ('genslblockfkine:SymbolicsNotFound','Save symbolic expressions to disk first!')
end

blockaddress = [CGen.slib,'/',symname];          % treat intermediate transformations separately
if ~isempty(find_system(CGen.slib,'SearchDepth',1,'Name',symname))                    % Delete previously generated block
    delete_block(blockaddress);
    save_system;
end

symexpr2slblock(blockaddress,tmpStruct.(symname).T,'vars',{q});

CGen.logmsg('\t%s\n',' done!');

%% Individual joint forward kinematics
CGen.logmsg([datestr(now),'\tGenerating forward kinematics Simulink block up to joint']);
for iJoints=1:CGen.rob.n

    CGen.logmsg(' %i ',iJoints);
    symname = ['T0_',num2str(iJoints)];
    fname = fullfile(CGen.sympath,[symname,'.mat']);

    tmpStruct = struct;
    tmpStruct = load(fname);

    funFileName = fullfile(CGen.robjpath,[symname,'.m']);
    q = CGen.rob.gencoords;

    blockaddress = [CGen.slib,'/',symname];          % treat intermediate transformations separately
    if doesblockexist(CGen.slib,symname)
        delete_block(blockaddress);
        save_system;
    end

    symexpr2slblock(blockaddress,tmpStruct.(symname).T,'vars',{q});

end
CGen.logmsg('\t%s\n',' done!');

%% Cleanup
% Arrange blocks
distributeblocks(CGen.slib);

% Lock, save and close library
set_param(CGen.slib,'lock','on');
save_system(CGen.slib,CGen.slibpath);
close_system(CGen.slib);

end

author: petercorke

repo: robotics-toolbox-matlab