text:
function [s, err_mse, iter_time]=greed_omp_pinv(x,A,m,varargin)
% greed_omp_pinv: Orthogonal Matching Pursuit algorithm based on matlab
% pinv solution. Use not recommended, provided for reference only!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Usage
% [s, err_mse, iter_time]=greed_omp_pinv(x,P,m,'option_name','option_value')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Input
%   Mandatory:
%               x   Observation vector to be decomposed
%               P   Either:
%                       1) An nxm matrix (n must be dimension of x)
%                       2) A function handle (type "help function_format"
%                          for more information)
%                          Also requires specification of P_trans option.
%                       3) An object handle (type "help object_format" for
%                          more information)
%               m   length of s
%
%   Possible additional options:
%   (specify as many as you want using 'option_name','option_value' pairs)
%   See below for explanation of options:
%__________________________________________________________________________
%   option_name    |     available option_values                | default
%--------------------------------------------------------------------------
%   stopCrit       | M, corr, mse, mse_change                   | M
%   stopTol        | number (see below)                         | n/4
%   P_trans        | function_handle (see below)                |
%   maxIter        | positive integer (see below)               | n
%   verbose        | true, false                                | false
%   start_val      | vector of length m                         | zeros
%
%   Available stopping criteria :
%               M           -   Extracts exactly M = stopTol elements.
%               corr        -   Stops when maximum correlation between
%                               residual and atoms is below stopTol value.
%               mse         -   Stops when mean squared error of residual
%                               is below stopTol value.
%               mse_change  -   Stops when the change in the mean squared
%                               error falls below stopTol value.
%
%   stopTol: Value for stopping criterion.
%
%   P_trans: If P is a function handle, then P_trans has to be specified and
%            must be a function handle.
%
%   maxIter: Maximum of allowed iterations.
%
%   verbose: Logical value to allow algorithm progress to be displayed.
%
%   start_val: Allows algorithms to start from partial solution.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outputs
%    s              Solution vector
%    err_mse        Vector containing mse of approximation error for each
%                   iteration
%    iter_time      Vector containing times for each iteration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Description
%   greed_omp_pinv performs a greedy signal decomposition.
%   In each iteration a new element is selected depending on the inner
%   product between the current residual and columns in P.
%   The non-zero elements of s are approximated by orthogonally projecting
%   x onto the selected elements in each iteration.
%   This implementation uses matlab's pinv command to calculate the required
%   projection. This is very slow and only works for small matrices. It is
%   therefore only provided for reference.
%
% See Also
%   greed_qr, greed_omp_chol, greed_omp_cg, greed_omp_cgp,
%   greed_omp_linsolve, greed_gp, greed_nomp
%
% Copyright (c) 2007 Thomas Blumensath
%
% The University of Edinburgh
% Email: thomas.blumensath@ed.ac.uk
% Comments and bug reports welcome
%
% This file is part of sparsity Version 0.1
% Created: April 2007
%
% Part of this toolbox was developed with the support of EPSRC Grant
% D000246/1
%
% Please read COPYRIGHT.m for terms and conditions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    Default values and initialisation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n1 n2]=size(x);
if n2 == 1
    n=n1;
elseif n1 == 1
    x=x';
    n=n2;
else
   display('x must be a vector.');
   return
end

sigsize     = x'*x/n;
initial_given=0;
err_mse     = [];
iter_time   = [];
STOPCRIT    = 'M';
STOPTOL     = ceil(n/4);
MAXITER     = n;
verbose     = false;
s_initial   = zeros(m,1);
vectnfact   = ones(m,1);

if verbose
   display('Initialising...')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Output variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch nargout
    case 3
        comp_err=true;
        comp_time=true;
    case 2
        comp_err=true;
        comp_time=false;
    case 1
        comp_err=false;
        comp_time=false;
    case 0
        error('Please assign output variable.')
    otherwise
        error('Too many output arguments specified')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       Look through options
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Put option into nice format
Options={};
OS=nargin-3;
c=1;
for i=1:OS
    if isa(varargin{i},'cell')
        CellSize=length(varargin{i});
        ThisCell=varargin{i};
        for j=1:CellSize
            Options{c}=ThisCell{j};
            c=c+1;
        end
    else
        Options{c}=varargin{i};
        c=c+1;
    end
end
OS=length(Options);
if rem(OS,2)
   error('Something is wrong with argument name and argument value pairs.')
end

for i=1:2:OS
   switch Options{i}
        case {'stopCrit'}
            if (strmatch(Options{i+1},{'M'; 'corr'; 'mse'; 'mse_change'},'exact'));
                STOPCRIT    = Options{i+1};
            else error('stopCrit must be char string [M, corr, mse, mse_change]. Exiting.'); end
        case {'stopTol'}
            if isa(Options{i+1},'numeric') ; STOPTOL     = Options{i+1};
            else error('stopTol must be number. Exiting.'); end
        case {'P_trans'}
            if isa(Options{i+1},'function_handle'); Pt = Options{i+1};
            else error('P_trans must be function _handle. Exiting.'); end
        case {'maxIter'}
            if isa(Options{i+1},'numeric'); MAXITER     = Options{i+1};
            else error('maxIter must be a number. Exiting.'); end
        case {'verbose'}
            if isa(Options{i+1},'logical'); verbose     = Options{i+1};
            else error('verbose must be a logical. Exiting.'); end
        case {'vecNormFac'}
            if isa(Options{i+1},'numeric')& length(Options{i+1}) == m , vectnfact = Options{i+1};
            else error('verbose must be a logical. Exiting.'); end
        case {'start_val'}
            if isa(Options{i+1},'numeric') & length(Options{i+1}) == m ;
                s_initial     = Options{i+1};
                initial_given=1;
            else error('start_val must be a vector of length m. Exiting.'); end
        otherwise
            error('Unrecognised option. Exiting.')
   end
end

if strcmp(STOPCRIT,'M')
    maxM=STOPTOL;
else
    maxM=MAXITER;
end

if nargout >=2
    err_mse = zeros(maxM,1);
end
if nargout ==3
    iter_time = zeros(maxM,1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Make P and Pt functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if          isa(A,'float')      P =@(z) A*z;  Pt =@(z) A'*z;
elseif      isobject(A)         P =@(z) A*z;  Pt =@(z) A'*z;
elseif      isa(A,'function_handle')
    try
        if          isa(Pt,'function_handle'); P=A;
        else        error('If P is a function handle, Pt also needs to be a function handle. Exiting.'); end
    catch error('If P is a function handle, Pt needs to be specified. Exiting.'); end
else        error('P is of unsupported type. Use matrix, function_handle or object. Exiting.'); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Do we start from zero or not?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if initial_given ==1;
    IN          = find(s_initial);
    s=zeros(m,1);
    if isa(A,'function_handle') || isobject(A)
         Pmat=zeros(n,length(IN));
         for i=1:length(IN)
             mask=zeros(m,1);
             mask(IN(i))=1;
             Pmat(:,i)=P(mask);
         end
         s(IN)=pinv(Pmat)*x;
     else
         s(IN)=pinv(A(:,IN))*x;
     end
    Residual    = x-P(s);
    oldERR      = Residual'*Residual/n;

else
    IN          = [];
    Residual    = x;
    s           = s_initial;
    sigsize     = x'*x/n;
    oldERR      = sigsize;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 Random Check to see if dictionary is normalised
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%         mask=zeros(m,1);
%         mask(ceil(rand*m))=1;
%         nP=norm(P(mask));
%         if abs(1-nP)>1e-3;
%             display('Dictionary appears not to have unit norm columns.')
%         end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Main algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if verbose
   display('Main iterations...')
end
tic
t=0;
DR=Pt(Residual).*vectnfact;
done = 0;
iter=1;
while ~done
     DR(IN)=0;
     [v I]=max(abs(DR));
     IN=[IN I];
     if isa(A,'function_handle') || isobject(A)
         Pmat=zeros(n,length(IN));
         for i=1:length(IN)
             mask=zeros(m,1);
             mask(IN(i))=1;
             Pmat(:,i)=P(mask);
         end
         s(IN)=pinv(Pmat)*x;
     else
         s(IN)=pinv(A(:,IN))*x;
     end
         Residual=x-P(s);
         DR=Pt(Residual).*vectnfact;

    ERR=Residual'*Residual/n;
     if comp_err
         err_mse(iter)=ERR;
     end

     if comp_time
         iter_time(iter)=toc;
     end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        Are we done yet?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     if strcmp(STOPCRIT,'M')
         if iter >= STOPTOL
             done =1;
         elseif verbose && toc-t>10
            display(sprintf('Iteration %i. --- %i iterations to go',iter ,STOPTOL-iter))
            t=toc;
         end
    elseif strcmp(STOPCRIT,'mse')
         if comp_err
            if err_mse(iter)<STOPTOL;
                done = 1;
            elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse',iter ,err_mse(iter)))
                t=toc;
            end
         else
             if ERR<STOPTOL;
                done = 1;
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse',iter ,ERR))
                t=toc;
             end
         end
     elseif strcmp(STOPCRIT,'mse_change') && iter >=2
         if comp_err && iter >=2
              if ((err_mse(iter-1)-err_mse(iter))/sigsize <STOPTOL);
                done = 1;
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse change',iter ,(err_mse(iter-1)-err_mse(iter))/sigsize ))
                t=toc;
             end
         else
             if ((oldERR - ERR)/sigsize < STOPTOL);
                done = 1;
             elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i mse change',iter ,(oldERR - ERR)/sigsize))
                t=toc;
             end
         end
     elseif strcmp(STOPCRIT,'corr')
          if max(abs(DR)) < STOPTOL;
             done = 1;
          elseif verbose && toc-t>10
                display(sprintf('Iteration %i. --- %i corr',iter ,max(abs(DR))))
                t=toc;
          end
     end

    % Also stop if residual gets too small or maxIter reached
     if comp_err
         if err_mse(iter)<1e-16
             display('Stopping. Exact signal representation found!')
             done=1;
         end
     else

         if iter>1
             if ERR<1e-16
                 display('Stopping. Exact signal representation found!')
                 done=1;
             end
         end
     end

     if iter >= MAXITER
         display('Stopping. Maximum number of iterations reached!')
         done = 1;
     end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    If not done, take another round
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     if ~done
        iter=iter+1;
        oldERR=ERR;
     end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  Only return as many elements as iterations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargout >=2
    err_mse = err_mse(1:iter);
end
if nargout ==3
    iter_time = iter_time(1:iter);
end

if verbose
   display('Done')
end

% Change history
%
% 8 of Februray: Algo does no longer stop if dictionary is not normaliesd.

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
function params = makeRoughBars(params,id)
% makeRoughBars - make rough bar mapping stimulus
%
% params = makeRoughBars(stimparams,stimulus_id);
%
% Hack: copied almost straight from stimulus code. Little
% flexibility - good enough for now.
%
% 2007/09 SOD: wrote it.

if notDefined('params');     error('Need params'); end;
if notDefined('id');         id = 1;                   end;

fprintf(1,'[%s]:Processing:...',mfilename);drawnow;

fixMov = -5;
eyeMov = [20 10 20./180*pi]./2;
X = params.analysis.X;
Y = params.analysis.Y;

stimlimits = [14 14/600*800];

% make images
images=zeros(numel(X),21);
imii  = 1;

%--- fix: left up
X2 = X-fixMov;
Y2 = Y-fixMov;
% stim: right
X3 = X2 - eyeMov(1);
Y3 = Y2;
t  = atan2 (Y3, X3);
mask = t>=-pi/2+eyeMov(3) & t<=pi/2-eyeMov(3) &...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii) = mask(:);
imii = imii+2;

% stim: bottom
X3 = X2;
Y3 = Y2 - eyeMov(2);
t  = atan2 (Y3, X3);
mask = t>=0+eyeMov(3) & t<=pi-eyeMov(3) &...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii) = mask(:);
imii = imii+2;

%-- fix: right up
X2 = X+fixMov;
Y2 = Y-fixMov;
% stim: left
X3 = X2 + eyeMov(1);
Y3 = Y2;
t  = atan2 (Y3, X3);
mask = (t<=-pi/2-eyeMov(3) | t>=pi/2+eyeMov(3)) &...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii)= mask(:);
imii = imii+2;

% stim: bottom
X3 = X2;
Y3 = Y2 - eyeMov(2);
t  = atan2 (Y3, X3);
mask = t>=0+eyeMov(3) & t<=pi-eyeMov(3)&...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii) = mask(:);
imii=imii+2;

%--- fix: center
X2 = X;
Y2 = Y;
% vertical
mask = abs(X2)<=eyeMov(1)./2 & ...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii) =  mask(:);
imii=imii+2;
% horizontal
mask = abs(Y2)<=eyeMov(2)./2 &...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii) =  mask(:);
imii=imii+2;

%--- fix: left down
X2 = X-fixMov;
Y2 = Y+fixMov;
% stim: right
X3 = X2 - eyeMov(1);
Y3 = Y2;
t  = atan2 (Y3, X3);
mask = t>=-pi/2+eyeMov(3) & t<=pi/2-eyeMov(3)&...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii) =  mask(:);
imii=imii+2;
% stim: top
X3 = X2;
Y3 = Y2 + eyeMov(2);
t  = atan2 (Y3, X3);
mask = t>=-pi+eyeMov(3) & t<=0-eyeMov(3)&...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii) = mask(:);
imii=imii+2;

%-- fix: right down
X2 = X+fixMov;
Y2 = Y+fixMov;
% stim: left
X3 = X2 + eyeMov(1);
Y3 = Y2;
t  = atan2 (Y3, X3);
mask = (t<=-pi/2-eyeMov(3) | t>=pi/2+eyeMov(3)) &...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii) = mask(:);
imii=imii+2;

% stim: top
X3 = X2;
Y3 = Y2 + eyeMov(2);
t  = atan2 (Y3, X3);
mask = t>=-pi+eyeMov(3) & t<=0-eyeMov(3) &...
       abs(X)<=stimlimits(2) & ...
       abs(Y)<=stimlimits(1);
images(:,imii) = mask(:);
%imii=imii+2;

% make sequence
sequence = [5:8 9:12 13:16 21,21 17:20 9:12 3,4,1,2 21,21 15,16 13,14 11,12 9,10 7,8 5,6 21,21 19,20 17,18 11,12 9,10 1:4 21,21];
fixseq   = ceil(sequence./4);

sequence = reshape(sequence,2,numel(sequence)./2)';% pair up
sequence = sequence(:,1);
sequence = repmat(sequence',5,1);
sequence = sequence(:);

fixseq = reshape(fixseq,2,numel(fixseq)./2)';% pair up
fixseq = fixseq(:,1);
fixseq = repmat(fixseq,1,5);
fixseq = fixseq';
fixseq = fixseq(:);

% fill in fixseq with nearest values
while sum(fixseq==6),
    a = diff(fixseq==6);
    ii=find(a==1);
    fixseq(ii+1)= fixseq(ii);
    ii=find(a==-1);
    fixseq(ii) = fixseq(ii+1);
end

nuisance = zeros(numel(fixseq),10);
change = [0;diff(fixseq)~=0] -  [diff(fixseq)~=0;0];
counter =1;
for n=1:5,
    nuisance(:,counter) = fixseq==n & change ==1;
    if fixseq(1)==n, nuisance(1,counter) = 1; end;
    counter = counter + 1;
    nuisance(:,counter) = fixseq==n & change ==-1;
    if fixseq(end)==n, nuisance(end,counter) = 1; end;
    counter = counter + 1;
end

% all images
img = images(:,sequence);
preimg = img(:,1+end-params.stim(id).prescanDuration:end);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
doEyeCorrection = params.stim(id).stimStart(1)~=0;
doVelocityWeighting = false;
%velocityMax     = 200;
if doEyeCorrection,
    fprintf(1,'\n[%s]:Correcting for eye-movements.\n',mfilename);
else
    fprintf(1,'\n[%s]:NOT correcting for eye-movements.\n',mfilename);
end
% incorporate eyemovements
newimg = zeros(size(img));

% HACK: stored in non-used stimDir
scans = params.stim(id).stimStart;
scans = max(scans,1);
for scanid=scans,
    eyeTrackData = sprintf('%s/eyeMovements/scan%d_eyetrackdata.mat',pwd,scanid);
    e = eyetrack2deg_roughBars(eyeTrackData);
    if doVelocityWeighting,
        stepsize = mean(diff(unique(params.analysis.Y)));
        e.velocity = e.velocity/2.3548*stepsize/4;
    end
    if ~doEyeCorrection,
        % this provides a visual check of the motion correction, because
        % fixation is moved during the scan.
        e.eye = e.fixation;
        e.blinks = false(size(e.blinks));
    end

    % resample to stimulus sequence
    stimindex    = linspace(0,numel(e.blinks),numel(e.blinks)/(e.framerate.*params.stim(id).framePeriod)+1)+1;
    stimduration = e.framerate.*.2;
    X = reshape(params.analysis.X,[1 1]*sqrt(numel(params.analysis.X)));
    Y = reshape(params.analysis.Y,[1 1]*sqrt(numel(params.analysis.Y)));

    for n=1:size(img,2),
        if sum(img(:,n))~=0,
            tmp = zeros(size(X));
            zi  = tmp;
            tmp(:) = img(:,n);
            for n2=stimindex(n):stimindex(n)+stimduration-1,
                if ~e.blinks(n2)

                    % resample to 'corrected' position + give weight based
                    % on velocity.
                    newtmp = interp2(X,Y,tmp,X+e.eye(n2,1),Y+e.eye(n2,2),'*linear',0);
                    if doVelocityWeighting && doEyeCorrection,
                        % assume velocity = fwhm
                        if e.velocity(n2) > 1,
                            f = fspecial('gaussian',ceil(e.velocity(n2)*6),e.velocity(n2));
                            newtmp = imfilter(newtmp,f,'same');
                        end
                    end
                    zi = zi + newtmp;%.*w(n); % old way
                end
            end
            newimg(:,n) = newimg(:,n) + zi(:)./stimduration;
        end
    end
end
newimg = newimg./numel(scans);

% append prestimuli
params.stim(id).images = cat(2,preimg,newimg);
% store nuisance parameters convolved with the Hrf
params.stim(id).nuisance = [nuisance(1+end-params.stim(id).prescanDuration:end,:); nuisance];
params.stim(id).nuisance = filter(params.analysis.Hrf{id},1,params.stim(id).nuisance);
params.stim(id).nuisance = rmAverageTime(params.stim(id).nuisance, ...
                                         params.stim(id).nUniqueRep);
params.stim(id).nuisance = params.stim(id).nuisance(params.stim(id).prescanDuration+1:end,:);
params.stim(id).nuisance = params.stim(id).nuisance';
return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [state] = CliffLeftSensorRoomba(serPort);
%[state] = CliffLeftSensorRoomba(serPort)
% Specifies the state of the Cliff Left sensor
% Either triggered or not triggered

% By; Joel Esposito, US Naval Academy, 2011

%Initialize preliminary return values
state = nan;

try

%Flush Buffer
N = serPort.BytesAvailable();
while(N~=0)
fread(serPort,N);
N = serPort.BytesAvailable();
end

warning off
global td

fwrite(serPort, [142]);  fwrite(serPort,9);
CliffLft = dec2bin(fread(serPort, 1));
state = bin2dec(CliffLft(end));

pause(td)
catch
    disp('WARNING:  function did not terminate correctly.  Output may be unreliable.')
end

author: petercorke

repo: robotics-toolbox-matlab


text:
function [RET, REL, RETREL, ndx] = LMprecisionRecall(retrievedBoundingBox, annotation, objectname, mindetectablesize)
% This function returns one point in the precision-recall curve.
% Provided a set of detected bounding boxes, it indicates how may objects
% have been detected (RETREL).
%
% The standard measures for IR are recall and precision. Assuming that:
%
%    * RET is the set of all items the system has retrieved for a specific inquiry;
%    * REL is the set of relevant items for a specific inquiry;
%    * RETREL is the set of the retrieved relevant items
%
% then precision and recall measures are obtained as follows:
%
%    precision = RETREL / RET
%    recall = RETREL / REL

% Search the target object in the annotation
j = LMobjectindex(annotation, objectname);

% Do not consider deleted objects:
j = j(find(strcmp({annotation.object(j).deleted}, '0')));

Ninstances = length(j);

% Extract the bounding boxes for each target present in the image
BoundingBox = []; REL = 0;
for i = 1:Ninstances
    [X,Y] = getLMpolygon(annotation.object(j(i)).polygon);
    BoundingBox = [BoundingBox; min(X) max(X) min(Y) max(Y)];

    % If we detect an object that we thought was too small, we should not
    % penalize performance for this. So, we will consider that the only
    % relevant targets are the ones with a size larger than the minimal
    % detectable size.
    if max(Y)-min(Y)>mindetectablesize(1) & max(X)-min(X)>mindetectablesize(2) & max(Y)>0 & max(X)>0
        REL = REL + 1;
    end
end

RET = size(retrievedBoundingBox, 1);

if REL == 0
    RETREL = 0;
    return
end

cxO = mean(BoundingBox(:,1:2),2);
cyO = mean(BoundingBox(:,3:4),2);
DxO = diff(BoundingBox(:,1:2),1,2);
DyO = diff(BoundingBox(:,3:4),1,2);

cxR = mean(retrievedBoundingBox(:,1:2),2);
cyR = mean(retrievedBoundingBox(:,3:4),2);
DxR = diff(retrievedBoundingBox(:,1:2),1,2);
DyR = diff(retrievedBoundingBox(:,3:4),1,2);

ndx = [];
for i = 1:RET
    d = sqrt(((cxR(i) - cxO)./DxO).^2 + ((cyR(i) - cyO)./DyO).^2)<.5 ...
        & max(DxO/DxR(i), DxR(i)./DxO)<1.5 ...
        & max(DyO/DyR(i), DyR(i)./DyO)<1.5;
    n = find(d);
    if length(n)>0
        ndx(i) = n(1);
    else
        ndx(i) = 0;
    end
end
RETREL = sum(unique(ndx)>0); % each object can be detected only once

% If we detect an object that we thought was too small, we should not
% penalize performance for this. So, if we detected one object that was not
% considered before within the relevan set, then we will move it into the
% relevant set:
REL = max(REL, RETREL);

author: CSAILVision

repo: LabelMeToolbox


text:
function [coeffs,sigma2]=multilinRegress(zSamp,y,w,hasIntercept)
%%MULTILINREGRESS Perform multilinear regression. This means, given a set
%         of multidimensional points of the form [x1;x2;x3;...;xn],
%         determine the line/plane/hyperplane best fitting the points
%         such that
%         xn=coeffs(1)*x(2)+coeffs(2)*x(3)+....coeffs(n-1)*x(n-1)+coeffs(n)
%         if there is an intercept or
%         xn=coeffs(1)*x(2)+coeffs(2)*x(3)+....coeffs(n)*x(n).
%         if it is known there is a zero intercept.
%         The "best" fitting line/plane/hyperplane is the one that
%         minimizes
%         sum_{i=1}^N w(i)*(dot(coeffs(1:(n-1)),zSamp(1:(n-1),i))+coeffs(n)-zSamp(n,i))^2
%         if there is an intercept and y is not provided or that minimizes
%         sum_{i=1}^N w(i)*(dot(coeffs(1:(n-1)),zSamp(:,i))+coeffs(n)-y(i))^2
%         if y is provided and there is a y intercept. If there is no y
%         intercept, then the coeffs(n) term goes away in both instances.
%         The w(i) terms are all positive weights. For unweighted fusion,
%         w(i)=1.
%
%INPUTS: There are two parameterizations for zSamp and y. If y is an empty
%        matrix, then zSamp is a numDimsXnumPoints set of numPoints data
%        samples numDims>=2, as given above. On the other hand, if y is
%        provided, then y corresponds to zSamp(n,:) in the previous
%        parameterization (and can be a row or column vector) and zSamp thus
%        lacks the final row that it has in the initial parameterization.
%        Thus, zSamp will be a (numDims-1)XnumPoints matrix when y is
%        given.
%      w If weighted linear least squares optimization is supposed to be
%        performed, this is a length numPoints vector of values >0
%        containing the weight of each sample. If the regression is
%        performed with un unweighted cost function, this parameter can be
%        omitted or an empty matrix can be passed.
% hasIntercept A boolean value indicating whether the fit has a nonzero y
%        intercept. The default if omitted or an empty matrix is passed is
%        true.
%
%OUTPUTS: coeffs The numDimsX1 set of weights for the linear equation
%                described above. If hasIntercept=true, then the last
%                element is the additive constant.
%         sigma2 The residual mean square.  This is SSE/(numPoints-numDims)
%                where SSE is the vector of ordinary residuals e'*diag(w)*e
%                where e=(y(:)-x'*coeffs), where x is
%                [zSamp;ones(1,numPoints)] if an intercept is desired and y
%                is given, or x=zSamp if an intercept is not desired and y
%                is given. Similar expressions for when y is not given
%                hold, whereby y is replaced by the last row of zSamp. This
%                is the quantity being minimized in the optimization.
%
%Formulae for this type of regression are given in Appendix A.7.1 and A.7.3
%of [1] for 2D and 3D measurements without the weights. This function just
%continues the pattern implementing the results for an arbitrary number of
%dimensions, for weights and a modification allowing for there to be a
%known zero y intercept is also added. The derivation is just basic
%differential calculus. Note that this function will fail if the plane
%being estimated does not depend on the independent parameter. In 2D, this
%would correspond to a vertical line (e.g. x=4).
%
%EXAMPLE 1:
%Here is an example with a nonzero y intercept using both input formats.
% coeffs=[4;6;2];
% numPts=20;
% xVals=linspace(-4,4,numPts);
% yVals=linspace(-4,4,numPts);
% [x,y]=ndgrid(xVals,yVals);
% xyPts=[x(:).';y(:).'];
% zPts=sum(bsxfun(@times,coeffs(1:2),xyPts),1)+coeffs(3);
% zPts=zPts+0.1*randn(1,numPts^2); %Add noise.
% zSamp=[xyPts;zPts];
% [coeffsEst,sigma2]=multilinRegress(zSamp)
% [coeffsEstAlt,sigma2Alt]=multilinRegress(xyPts,zPts)
%One finds that coeffsEst is the same as coeffsEstAlt and is close to
%coeffs. A perfect fit (wihtin finite precision limitaitions) is achieved
%if there is zero noise.
%
%EXAMPLE 2:
%Here is an example when it is known that the y intercept is zero.
% coeffs=[4;6];
% numPts=20;
% xVals=linspace(-4,4,numPts);
% yVals=linspace(-4,4,numPts);
% [x,y]=ndgrid(xVals,yVals);
% xyPts=[x(:).';y(:).'];
% zPts=sum(bsxfun(@times,coeffs(1:2),xyPts),1);
% zPts=zPts+0.1*randn(1,numPts^2); %Add noise.
% zSamp=[xyPts;zPts];
% [coeffsEst,sigma2]=multilinRegress(zSamp,[],[],false)
% [coeffsEstAlt,sigma2Alt]=multilinRegress(xyPts,zPts,[],false)
%Again, the results are comparable to those in problem 1 with the two input
%formulations equal and close to coeffs.
%
%REFERENCES:
%[1] P. J. Schneider and D. H. Eberly, Geometric Tools for Computer
%    Graphics. Amsterdam: Morgan Kaufmann Publishers, 2003.
%
%April 2017 David F. Crouse, Naval Research Laboratory, Washington D.C.
%(UNCLASSIFIED) DISTRIBUTION STATEMENT A. Approved for public release.

numDims=size(zSamp,1);
numPoints=size(zSamp,2);

if(nargin<4||isempty(hasIntercept))
    hasIntercept=true;
end

if(nargin<3)
   w=[];
end

if(nargin>1&&~isempty(y))
    if(hasIntercept)
        x=[zSamp;ones(1,numPoints)];
        p=numDims+1;
    else
        x=zSamp;
        p=numDims;
    end
else
    if(hasIntercept)
        x=[zSamp(1:(numDims-1),:);ones(1,numPoints)];
        p=numDims;
    else
        x=zSamp(1:(numDims-1),:);
        p=numDims-1;
    end

    y=zSamp(numDims,:);
    y=y(:);
end

if(isempty(w))
    A=x*x';
    b=x*y(:);
else
    %In Matlab, x*eye(size(x,2))*x'-x*x' and
    %bsxfun(@times,x,ones(numPoints,1).')*x'-x*x' are not numerically zero,
    %which means that just using w with weights all equal to 1 will change
    %the last few bits of the result.

    w=w(:);
    A=bsxfun(@times,w.',x)*x';
    b=x*(w.*y(:));
end

coeffs=A\b;

if(nargout>1)
    if(nargin>1&&~isempty(y))
        e=y(:)-x'*coeffs;
    else
        e=vec(zSamp(numDims,:))-x'*coeffs;
    end

    if(isempty(w))
        %The residual sum of squares.
        sigma2=e'*e/(numPoints-p);
    else%The weighted residual sum of squares.
        sigma2=(w.*e)'*e/(numPoints-p);
    end
end
end

%LICENSE:
%
%The source code is in the public domain and not licensed or under
%copyright. The information and software may be used freely by the public.
%As required by 17 U.S.C. 403, third parties producing copyrighted works
%consisting predominantly of the material produced by U.S. government
%agencies must provide notice with such work(s) identifying the U.S.
%Government material incorporated and stating that such material is not
%subject to copyright protection.
%
%Derived works shall not identify themselves in a manner that implies an
%endorsement by or an affiliation with the Naval Research Laboratory.
%
%RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF THE
%SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY THE NAVAL
%RESEARCH LABORATORY FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE ACTIONS
%OF RECIPIENT IN THE USE OF THE SOFTWARE.

author: USNavalResearchLaboratory

repo: TrackerComponentLibrary


text:
function y = glm_fit_model(x,k,v);
%
% y = glm_fit_model(x,k,v);
%
% A model used to generate a 'synthetic'
% autocorrelation function which may estimate
% temporally-correlated noise for within-
% skull voxels: see Burock and Dale, HBM 2000
% for more info.
%
% This is just broken off as a separate function
% for use in calling matlab's fminsearch routine.
%
% original code by gb, 11/04
% updated by ras, 05/05
y = norm((1 - x(1))*(x(2).^(1:k)') - v)^2;

return


text:
function result = glmnetPredict(object,s)

%--------------------------------------------------------------------------
% glmnetPredict.m: make predictions from a "glmnet" object.
%--------------------------------------------------------------------------
%
% DESCRIPTION:
%    Similar to other predict methods, this functions predicts fitted
%    values, logits, coefficients and more from a fitted "glmnet" object.
%
% USAGE:
%    glmnetPredict(object)
%    glmnetPredict(object, type)
%    glmnetPredict(object, type, newx)
%    glmnetPredict(object, type, newx, s)
%
% INPUT ARGUMENTS:
% fit         Fitted "glmnet" model object.
% type        Type of prediction required. Type "link" gives the linear
%             predictors for "binomial" or "multinomial" models; for
%             "gaussian" models it gives the fitted values. Type "response"
%             gives the fitted probabilities for "binomial" or
%             "multinomial"; for "gaussian" type "response" is equivalent
%             to type "link". Type "coefficients" computes the coefficients
%             at the requested values for s. Note that for "binomial"
%             models, results are returned only for the class corresponding
%             to the second level of the factor response. Type "class"
%             applies only to "binomial" or "multinomial" models, and
%             produces the class label corresponding to the maximum
%             probability. Type "nonzero" returns a list of the indices of
%             the nonzero coefficients for each value of s.
% newx        Matrix of new values for x at which predictions are to be
%             made. Must be a matrix; This argument is not used for
%             type=c("coefficients","nonzero")
% s           Value(s) of the penalty parameter lambda at which predictions
%             are required. Default is the entire sequence used to create
%             the model.
%
% DETAILS:
%    The shape of the objects returned are different for "multinomial"
%    objects. glmnetCoef(fit, ...) is equivalent to glmnetPredict(fit, "coefficients", ...)
%
% LICENSE: GPL-2
%
% DATE: 14 Jul 2009
%
% AUTHORS:
%    Algorithm was designed by Jerome Friedman, Trevor Hastie and Rob Tibshirani
%    Fortran code was written by Jerome Friedman
%    R wrapper (from which the MATLAB wrapper was adapted) was written by Trevor Hasite
%    MATLAB wrapper was written and maintained by Hui Jiang, jiangh@stanford.edu
%    Department of Statistics, Stanford University, Stanford, California, USA.
%
% REFERENCES:
%    Friedman, J., Hastie, T. and Tibshirani, R. (2009)
%    Regularization Paths for Generalized Linear Models via Coordinate Descent.
%    Journal of Statistical Software, 33(1), 2010
%
% SEE ALSO:
%    glmnet, glmnetSet, glmnetPrint, glmnetPlot and glmnetCoef methods.
%
% EXAMPLES:
%    x=randn(100,20);
%    y=randn(100,1);
%    g2=randsample(2,100,true);
%    g4=randsample(4,100,true);
%    fit1=glmnet(x,y);
%    glmnetPredict(fit1,'link',x(1:5,:),[0.01,0.005]') % make predictions
%    glmnetPredict(fit1,'coefficients')
%    fit2=glmnet(x,g2,'binomial');
%    glmnetPredict(fit2, 'response', x(2:5,:))
%    glmnetPredict(fit2, 'nonzero')
%    fit3=glmnet(x,g4,'multinomial');
%    glmnetPredict(fit3, 'response', x(1:3,:), 0.01)
%
% DEVELOPMENT:
%    14 Jul 2009: Original version of glmnet.m written.
%    20 Oct 2009: Fixed a bug in bionomial response, pointed out by Ramon
%                 Casanov from Wake Forest University.
%    26 Jan 2010: Fixed a bug in multinomial link and class, pointed out by
%                 Peter Rijnbeek from Erasmus University.
%    23 Jun 2010: Fixed a bug in multinomial with s, pointed out by
%                 Robert Jacobsen from Aalborg University.

% if nargin < 2
%     type = 'link';
% end
%
% if nargin < 3
%     newx = [];
% end
%
% if nargin < 4
%     s = object.lambda;
% end

% if strcmp(object.class, 'elnet')
a0=transpose(object.a0);
nbeta=[a0; object.beta];
% if nargin == 4
lambda=object.lambda;
lamlist=lambda_interp(lambda,s);
nbeta=nbeta(:,lamlist.left).*repmat(lamlist.frac',size(nbeta,1),1) +nbeta(:,lamlist.right).*(1-repmat(lamlist.frac',size(nbeta,1),1));
% end
% if strcmp(type, 'coefficients')
result = nbeta;
% elseif strcmp(type, 'link')
%     result = [ones(size(newx,1),1), newx] * nbeta;
% elseif strcmp(type, 'response')
%     result = [ones(size(newx,1),1), newx] * nbeta;
% elseif strcmp(type, 'nonzero')
%     result = nonzeroCoef(nbeta(2:size(nbeta,1),:), true);
% else
%     error('Unrecognized type');
% end
% elseif strcmp(object.class, 'lognet')
%
%     a0=transpose(object.a0);
%     nbeta=[object.a0; object.beta];
%     if nargin == 4
%         lambda=object.lambda;
%         lamlist=lambda_interp(lambda,s);
%         nbeta=nbeta(:,lamlist.left).*repmat(lamlist.frac',size(nbeta,1),1) +nbeta(:,lamlist.right).*(1-repmat(lamlist.frac',size(nbeta,1),1));
%     end
%     %%%   remember that although the fortran lognet makes predictions
%     %%%   for the first class, we make predictions for the second class
%     %%%   to avoid confusion with 0/1 responses.
%     %%%   glmnet flipped the signs of the coefficients
%     if strcmp(type,'coefficients')
%         result = nbeta;
%     elseif strcmp(type,'nonzero')
%         result = nonzeroCoef(nbeta(2:size(nbeta,1),:), true);
%     else
%         nfit = [ones(size(newx,1),1), newx] * nbeta;
%
%         if strcmp(type,'response')
%             pp=exp(-nfit);
%             result = 1./(1+pp);
%         elseif strcmp(type,'link')
%             result = nfit;
%         elseif strcmp(type,'class')
%             result = (nfit > 0) * 2 + (nfit <= 0) * 1;
%         else
%             error('Unrecognized type');
%         end
%     end
% elseif strcmp(object.class, 'multnet')
%     a0=object.a0;
%     nbeta=object.beta;
%     nclass=size(a0,1);
%     nlambda=length(s);
%     if nargin == 4
%         lambda=object.lambda;
%         lamlist=lambda_interp(lambda,s);
%         for i=1:nclass
%             kbeta=[a0(i,:); nbeta{i}];
% %            kbeta=kbeta(:,lamlist.left)*lamlist.frac +kbeta(:,lamlist.right)*(1-lamlist.frac);
%             kbeta=kbeta(:,lamlist.left).*repmat(lamlist.frac',size(kbeta,1),1)+kbeta(:,lamlist.right).*(1-repmat(lamlist.frac',size(kbeta,1),1));
%             nbeta{i}=kbeta;
%         end
%     else
%         for i=1:nclass
%             nbeta{i} = [a0(i,:);nbeta{i}];
%         end
%     end
%     if strcmp(type, 'coefficients')
%         result = nbeta;
%     elseif strcmp(type, 'nonzero')
%         for i=1:nclass
%             result{i}=nonzeroCoef(nbeta{i}(2:size(nbeta{i},1),:),true);
%         end
%     else
%         npred=size(newx,1);
%         dp = zeros(nclass,nlambda,npred);
%         for i=1:nclass
%             fitk = [ones(size(newx,1),1), newx] * nbeta{i};
%             dp(i,:,:)=dp(i,:,:)+reshape(transpose(fitk),1,nlambda,npred);
%         end
%         if strcmp(type, 'response')
%             pp=exp(dp);
%             psum=sum(pp,1);
%             result = permute(pp./repmat(psum,nclass,1),[3,1,2]);
%         elseif strcmp(type, 'link')
%             result=permute(dp,[3,1,2]);
%         elseif strcmp(type, 'class')
%             dp=permute(dp,[3,1,2]);
%             result = [];
%             for i=1:size(dp,3)
%                 result = [result, softmax(dp(:,:,i))];
%             end
%         else
%             error('Unrecognized type');
%         end
%     end
% else
%     error('Unrecognized class');
% end

%-------------------------------------------------------------
% End private function glmnetPredict
%-------------------------------------------------------------

function result = lambda_interp(lambda,s)
% lambda is the index sequence that is produced by the model
% s is the new vector at which evaluations are required.
% the value is a vector of left and right indices, and a vector of fractions.
% the new values are interpolated bewteen the two using the fraction
% Note: lambda decreases. you take:
% sfrac*left+(1-sfrac*right)

  if length(lambda)==1 % degenerate case of only one lambda
    nums=length(s);
    left=ones(nums,1);
    right=left;
    sfrac=ones(nums,1);
else
      s(s > max(lambda)) = max(lambda);
      s(s < min(lambda)) = min(lambda);
      k=length(lambda);
      sfrac =(lambda(1)-s)/(lambda(1) - lambda(k));
      lambda = (lambda(1) - lambda)/(lambda(1) - lambda(k));
      coord = interp1(lambda, 1:length(lambda), sfrac);
      left = floor(coord);
      right = ceil(coord);
      sfrac=(sfrac-lambda(right))./(lambda(left) - lambda(right));
      sfrac(left==right)=1;
  end
  result.left = left;
  result.right = right;
  result.frac = sfrac;

%-------------------------------------------------------------
% End private function lambda_interp
%-------------------------------------------------------------
%
% function result = softmax(x, gap)
% if nargin < 2
%     gap = false;
% end
% d = size(x);
% maxdist = x(:, 1);
% pclass = repmat(1, d(1), 1);
% for i =2:d(2)
%     l = x(:, i) > maxdist;
%     pclass(l) = i;
%     maxdist(l) = x(l, i);
% end
% if gap
%     x = abs(maxdist - x);
%     x(1:d(1), pclass) = x * repmat(1, d(2));
%     gaps = pmin(x);
% end
% if gap
%     result = {pclass, gaps};
% else
%     result = pclass;
% end

%-------------------------------------------------------------
% End private function softmax
%-------------------------------------------------------------

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function b = frank_rhs ( m, k )

%*****************************************************************************80
%
%% FRANK_RHS returns the FRANK right hand side.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 November 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the row dimension.
%
%    Input, integer K, the column dimension ( should be 2).
%
%    Output, real B(M,K), the right hand side matrix.
%
  b = zeros ( m, k );

  b(1:m,1) = 1.0;

  b(1  ,2) = m * ( m + 1 ) / 2;
  for i = 2 : m
    b(i,2) = ( m + 1 - i ) * ( m + 4 - i ) / 2;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function frst=show_progress(cnt,ttl,frst);

persistent hit

%Initialise hit only when it is first declared as an empty array
if frst==0
    hit=zeros(20,1);
    fprintf('10%% 20%% 30%% 40%% 50%% 60%% 70%% 80%% 90%% 100%%\n');
end

prp=cnt/ttl;
if prp>=0.05 & hit(1)==0
    fprintf('||');
    hit(1)=1;
elseif prp>=0.1 & hit(2)==0
    fprintf('||');
    hit(2)=1;
elseif prp>=0.15 & hit(3)==0
    fprintf('||');
    hit(3)=1;
elseif prp>=0.2 & hit(4)==0
    fprintf('||');
    hit(4)=1;
elseif prp>=0.25 & hit(5)==0
    fprintf('||');
    hit(5)=1;
elseif prp>=0.3 & hit(6)==0
    fprintf('||');
    hit(6)=1;
elseif prp>=0.35 & hit(7)==0
    fprintf('||');
    hit(7)=1;
elseif prp>=0.4 & hit(8)==0
    fprintf('||');
    hit(8)=1;
elseif prp>=0.45 & hit(9)==0
    fprintf('||');
    hit(9)=1;
elseif prp>=0.5 & hit(10)==0
    fprintf('||');
    hit(10)=1;
elseif prp>=0.55 & hit(11)==0
    fprintf('||');
    hit(11)=1;
elseif prp>=0.6 & hit(12)==0
    fprintf('||');
    hit(12)=1;
elseif prp>=0.65 & hit(13)==0
    fprintf('||');
    hit(13)=1;
elseif prp>=0.7 & hit(14)==0
    fprintf('||');
    hit(14)=1;
elseif prp>=0.75 & hit(15)==0
    fprintf('||');
    hit(15)=1;
elseif prp>=0.8 & hit(16)==0
    fprintf('||');
    hit(16)=1;
elseif prp>=0.85 & hit(17)==0
    fprintf('||');
    hit(17)=1;
elseif prp>=0.9 & hit(18)==0
    fprintf('||');
    hit(18)=1;
elseif prp>=0.95 & hit(19)==0
    fprintf('||');
    hit(19)=1;
elseif prp>=0.99 & hit(20)==0
    fprintf('||\n');
    hit(20)=1;
end

frst=1;

author: yetianmed

repo: subcortex

save_path: github-repos/MATLAB/yetianmed-subcortex


text:
function p = identityLine(ax,varargin)
%Draw an identity line on the current axis
%
%   p = identityLine(ax,varargin)
%
% Example:
%   plot(1:10,randn(1,10),'o')
%   identityLine(gca);
%   identityLine(gca,'linecolor',[1 0 0],'linestyle',':');
%
% (c) Stanford VISTA Team

%% if notDefined('ax'), ax = gca; end
p = inputParser;

% First optional argument
p.addRequired('ax',@isgraphics);

% Line color and style
vFunc = @(x) isnumeric(x) && logical(length(x) == 3);
p.addParameter('linecolor',[.5 .5 .5],vFunc);
p.addParameter('linestyle','--',@ischar);

% Parse inputs
p.parse(ax,varargin{:});
ax = p.Results.ax;
linecolor = p.Results.linecolor;
linestyle = p.Results.linestyle;

%% Set up axes
% Minimum and maximum of axes
xlim = get(ax,'xlim');
ylim = get(ax,'ylim');
mn = min(xlim(1),ylim(1));
mx = max(xlim(2),ylim(2));

% Here's the line
p = line([mn mx],[mn mx],'color',linecolor,'linestyle',linestyle);

% Set line properties.  These probably want to come in as an argument
set(p,'linewidth',2);
grid on

return


text:
function slide37

	a = [0 0 0;
		3 0 0;
		3 1 0;
		0 1 0;
		0 0 1;
		3 0 1;
		3 1 1;
		0 1 1];

	b = [1 2 6 5;
		2 3 7 6;
		3 4 8 7;
		4 1 5 8;
		1 2 3 4;
		5 6 7 8];

	p1 = patch('faces',b,...
		'vertices',a,...
		'facecolor',[.5 .5 .5],...
		'edgecolor',[1,1,1],...
		'facealpha',0.5);

	view(3)
	axis([-3 7 -3 5 -3 5])
	grid on

	while true

		V = get(p1,'vertices');
		V(5:8,3) = V(5:8,3) + 0.01;

		set(p1,'vertices',V)
		drawnow

	end

end

author: Sable

repo: mcbench-benchmarks


text:
function Ref = RefSelect(Population,k)
% Reference solutions selection by RSEA strategy

%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Cheng He

    k      = min(k,length(Population));
    PopObj = Population.objs;
	[FrontNO,MaxFNO] = NDSort(PopObj,k);
    Next = find(FrontNO<=MaxFNO);
    Pmin = min(PopObj,[],1) + 1e-6;
    Pmax = max(PopObj,[],1);
    if Pmax > Pmin
        PopObj = (PopObj-repmat(Pmin,size(PopObj,1),1))./repmat(Pmax-Pmin,size(PopObj,1),1);
    end

    %% Environmental selection
    Choose = LastSelection(PopObj(Next,:),ismember(Next,find(FrontNO<MaxFNO)),ceil(sqrt(k)),k);
    Ref    = Population(Next(Choose));
end

function Choose = LastSelection(PopObj,Choose,div,k)
% Select part of the solutions based on the radar grid

    %% Identify the extreme solutions
	[~,Extreme] = min(sqrt(sum(PopObj.^2,2)).*sqrt(1-(1-pdist2(PopObj,ones(1,size(PopObj,2)),'cosine')).^2),[],1); %Calculate the extreme points based on PBI
    Choose      = Choose | ismember(1:size(PopObj,1),Extreme);

    %% Calculate the convergence of each solution
	Con = sum(PopObj.^1,2).^1;
    Con = Con./max(Con);

    %% Calculate the radar grid of each solution
    [Site,RLoc] = RadarGrid(PopObj,div);
    RDis        = pdist2(RLoc,RLoc);
    RDis(logical(eye(length(RDis)))) = inf;
    CrowdG      = zeros(1,max(Site));
    temp        = tabulate(Site(Choose));
    CrowdG(temp(:,1)) = temp(:,2);

    %% Select k solutions
    while sum(Choose) < k
        % Delete outline solutions
        remainS  = find(~Choose);
        remainG  = unique(Site(remainS));
        bestG    = CrowdG(remainG) == min(CrowdG(remainG));
        current  = remainS(ismember(Site(remainS),remainG(bestG)));
        fitness  = 0.1.*size(PopObj,2).*Con(current) - min(RDis(current,Choose),[],2); % - 0.1.* min(Dis(current,Choose),[],2);
        [~,best] = min(fitness);
        Choose(current(best))       = true;
        CrowdG(Site(current(best))) = CrowdG(Site(current(best))) + 1;
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
%1. Load the image, for example
    image        = imread('img.bmp');

%2. Call this function to calculate the quality score:
    qualityscore = SSEQ(image)

author: thomas-koehler

repo: SupER


text:
%CREATECONCENTRICSPHERESTESTSET  Creates test set
%
%     [samples, responses] = cv.createConcentricSpheresTestSet(nsamples, nfeatures, nclasses)
%
% ## Input
% * __nsamples__ returned samples count.
% * __nfeatures__ returned features count.
% * __nclasses__ number of classes.
%
% ## Output
% * __samples__ returned samples array `nsamples-by-nfeatures`.
% * __sampClasses__ corresponding classes labels `nsamples-by-1`.
%
% See also: cv.randMVNormal, cv.randGaussMixture
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function hh = quadmesh(quad,x,y,z,varargin)
%QUADMESH Quadrilateral mesh plot.
%   QUADMESH(QUAD,X,Y,Z,C) displays the quadrilaterals defined in the M-by-4
%   face matrix QUAD as a mesh.  A row of QUAD contains indexes into
%   the X,Y, and Z vertex vectors to define a single quadrilateral face.
%   The edge color is defined by the vector C.
%
%   QUADMESH(QUAD,X,Y,Z) uses C = Z, so color is proportional to surface
%   height.
%
%   QUADMESH(TRI,X,Y) displays the quadrilaterals in a 2-d plot.
%
%   H = QUADMESH(...) returns a handle to the displayed quadrilaterals.
%
%   QUADMESH(...,'param','value','param','value'...) allows additional
%   patch param/value pairs to be used when creating the patch object.
%
%   See also PATCH.
%
% Script code based on copyrighted code from mathworks for TRIMESH.
% Allan P. Engsig-Karup, apek@mek.dtu.dk.

ax = axescheck(varargin{:});
ax = newplot(ax);

if nargin == 3 || (nargin > 4 && ischar(z))
  d = tri(:,[1 2 3 4 1])';
  if nargin == 3
    h = plot(ax, x(d), y(d));
  else
    h = plot(ax, x(d), y(d),z,varargin{1},varargin{2:end});
  end
  if nargout == 1, hh = h; end
  return;
end

start = 1;
if nargin>4 && rem(nargin-4,2)==1
  c = varargin{1};
  start = 2;
elseif nargin<3
  error(id('NotEnoughInputs'),'Not enough input arguments');
else
  c = z;
end

if ischar(get(ax,'color')),
  fc = get(gcf,'Color');
else
  fc = get(ax,'color');
end

h = patch('faces',quad,'vertices',[x(:) y(:) z(:)],'facevertexcdata',c(:),...
	  'facecolor',fc,'edgecolor',get(ax,'defaultsurfacefacecolor'),...
	  'facelighting', 'none', 'edgelighting', 'flat',...
      'parent',ax,...
	  varargin{start:end});
if ~ishold(ax), view(ax,3), grid(ax,'on'), end
if nargout == 1, hh = h; end

function str = id(str)
str = ['MATLAB:quadmesh:' str];

return

% % Example: plot quadmesh
% close all, clear all, clc
%
% x = [0 1 2 2 1 1 0 0];
% y = [0 0 0 1 1 2 2 1];
% for i = 1 : length(x)
%     text(x(i)+0.1,y(i)+0.1,sprintf('%d',i))
%     hold on
% end
% quad = [1 2 5 8;
%     2 3 4 5;
%     8 5 6 7];
%
% plot(x,y,'k.-')
% axis off
%
% z = sin(pi*x).*cos(pi*y);
% c = z;
% figure
%
% ax=newplot;
% fc = get(gcf,'Color');
% h = patch('faces',quad,'vertices',[x(:) y(:) z(:)],'facevertexcdata',c(:),...
%     'facecolor',fc,'edgecolor',get(ax,'defaultsurfacefacecolor'),...
%     'facelighting', 'none', 'edgelighting', 'flat',...
%     'parent',ax);


text:
%% faceToVertexMeasure
% Below is a demonstration of the features of the |patchPathAngles| function

%%
clear; close all; clc;

%% Syntax
% |[indAngles]=faceToVertexMeasure(F,V,ind,isClosedLoop);|

%% Description
% The |faceToVertexMeasure| function converts data for faces to data on
% vertices through averaging.

%% Examples

%%
% Plot settings
markerSize=150;

%% Example 1: Convert face data to vertex data
[F,V]=geoSphere(2,1);
VF=patchCentre(F,V);
CF=VF(:,1);
[CV]=faceToVertexMeasure(F,V,CF);

%%

cFigure;
subplot(1,2,1); hold on;
title('Color data on faces')
gpatch(F,V,CF);
scatterV(VF,markerSize,CF,'filled');
colormap gjet; colorbar;
axisGeom;
camlight headlight;

subplot(1,2,2); hold on;
title('Converted data on vertices')
hp=gpatch(F,V,CV);
hp.FaceColor='Interp';
scatterV(V,markerSize,CV,'filled');
colormap gjet; colorbar;
axisGeom;
camlight headlight;

drawnow;

%% Example 2: Convert face data on a mixed mesh to vertex data

%%
% Create an example of a mixed mesh consisting of triangles and
% quadrilateral faces.

% Create a triangular mesh
[F,V]=geoSphere(2,1);

% Converting to a quadrilateral mesh
optionStruct.maxAngleDeviation=45*(pi/180);
optionStruct.selectionMethod='best';
optionStruct.triangleConvert=0;
optionStruct.fourConnectConvert=0;
[F,V]=tri2quadGroupSplit(F,V,optionStruct);

VF=patchCentre(F,V);
CF=VF;
for q=1:1:numel(VF)
    CF{q}=VF{q}(:,1);
end
[CV]=faceToVertexMeasure(F,V,CF);

%%

cFigure;
subplot(1,2,1); hold on;
title('Color data on faces')
gpatch(F,V,CF);
for q=1:1:numel(VF)
    scatterV(VF{q},markerSize,CF{q},'filled');
end
colormap gjet; colorbar;
axisGeom;
camlight headlight;

subplot(1,2,2); hold on;
title('Converted data on vertices')
hp=gpatch(F,V,CV);
% hp.FaceColor='Interp';
scatterV(V,markerSize,CV,'filled');
colormap gjet; colorbar;
axisGeom;
camlight headlight;
drawnow;

%% Example 3: Convert multi-dimensional face data (e.g. on a mixed mesh) to vertex data

N=patchNormal(F,V); %get face normals
NV=faceToVertexMeasure(F,V,N); %Convert to vertex normals

%%

VF=patchCentre(F,V); %Get face centres for plotting

cFigure;
subplot(1,2,1); hold on;
title('Vector data on faces')
gpatch(F,V,CF);
for q=1:1:numel(F)
    quiverVec(VF{q},N{q},0.25,'k');
end
colormap gjet; colorbar;
axisGeom;
camlight headlight;

subplot(1,2,2); hold on;
title('Converted data on vertices')
gpatch(F,V,CV);
quiverVec(V,NV,0.25,'k');
colormap gjet; colorbar;
axisGeom;
camlight headlight;
drawnow;

%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program 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 General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function u = asinh(a)
%ASINH        Slope inverse hyperbolic sine asinh(a)
%

% written  12/06/98     S.M. Rump
% modified 04/04/04     S.M. Rump  set round to nearest for safety
% modified 04/06/05     S.M. Rump  rounding unchanged
% modified 09/28/08     S.M. Rump  check for rounding to nearest improved
% modified 08/26/12     S.M. Rump  global variables removed
%

  e = 1e-30;
  if 1+e==1-e                           % fast check for rounding to nearest
    rndold = 0;
  else
    rndold = getround;
    setround(0)
  end

  INTLAB_SLOPE = getappdata(0,'INTLAB_SLOPE');

  u = a;

  u.r = asinh(a.r);
  indexc = 1:INTLAB_SLOPE.NUMVAR;
  indexr = 2:INTLAB_SLOPE.NUMVAR+1;
  Xxs = hull(a.r(:,indexc),a.r(:,indexr));
  Index = 1:size(a.r.inf,1);

  index = all( a.r.sup<=0 , 2);
  if any(index)
    aindex.r = a.r(index,:);
    aindex.s = a.s(index,:);
    u.s(index,:) = slopeconvexconcave('asinh','1./sqrt(1+sqr(%))',aindex,1);
    Index(index) = 0;
  end

  index = all( a.r.inf>=0 , 2);
  if any(index)
    aindex.r = a.r(index,:);
    aindex.s = a.s(index,:);
    u.s(index,:) = slopeconvexconcave('asinh','1./sqrt(1+sqr(%))',aindex,0);
    Index(index) = 0;
  end

  if any(Index)
    Index( Index==0 ) = [];
    u.s(Index,:) = a.s(Index,:) ./ sqrt( 1 + sqr(Xxs(Index)) );
  end

  if rndold
    setround(rndold)
  end

repo: OpenMAS


text:
% SYNTAX:
% [yAvgStd, yAvgStdErr] = hmrG_SubjAvgStd(ySubjAvg)
%
% UI NAME:
% Subj_Average_Standard_Deviation_and_Standard_Error
%
% DESCRIPTION:
% Calculate avearge HRF standard deviation and standard error of all subjects in a group.
%
% INPUTS:
% ySubjAvg:
%
% OUTPUTS:
% yAvgStdOut: the standard deviation across subjects.
% yAvgStdErrOut: the standard error across subjects.
%
% USAGE OPTIONS:
% Subj_Average_Standard_Deviation_on_Concentration_Data:  [dcAvgStd, dcAvgStdErr]  = hmrG_SubjAvgStd(dcAvgSubjs)
% Subj_Average_Standard_Deviation_on_Delta_OD_Data:       [dodAvgStd, dodAvgStd] = hmrG_SubjAvgStd(dodAvgSubjs)
%

function [yAvgStdOut, yAvgStdErrOut] = hmrG_SubjAvgStd(ySubjAvg)

yAvgStdOut = DataClass().empty();
yAvgStdErrOut = DataClass().empty();

if isempty(ySubjAvg)
    return;
end
for iBlk = 1:length(ySubjAvg{1})
    yAvgStdOut(iBlk) = DataClass(ySubjAvg{1});
    yAvgStdErrOut(iBlk) = DataClass(ySubjAvg{1});
    foo = ySubjAvg{1}(iBlk).GetDataTimeSeries();
    dts = zeros(size(foo,1), size(foo,2), length(ySubjAvg));
    for iRun = 1:length(ySubjAvg)
        dts(:,:,iRun) = ySubjAvg{iRun}(iBlk).GetDataTimeSeries();
    end
    yAvgStdOut(iBlk).SetDataTimeSeries(std(dts,0,3,'omitnan'));
    yAvgStdErrOut(iBlk).SetDataTimeSeries(std(dts,0,3,'omitnan')/sqrt(length(ySubjAvg)-1));
end

author: BUNPC


text:
function out = densityremap(x)
% DENSITYREMAP Standard set of parameters for density remap
%
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

    out = uint8( amplitudetodensity(x) );
end
% //////////////////////////////////////////
% /// CLASSIFICATION: UNCLASSIFIED       ///
% //////////////////////////////////////////

author: ngageoint


text:
%-------------------------------------------------------------------------------------------------------------------%
%
% IB2d is an Immersed Boundary Code (IB) for solving fully coupled non-linear
% 	fluid-structure interaction models. This version of the code is based off of
%	Peskin's Immersed Boundary Method Paper in Acta Numerica, 2002.
%
% Author: Nicholas A. Battista
% Email:  nick.battista@unc.edu
% Date Created: May 27th, 2015
% Institution: UNC-CH
%
% This code is capable of creating Lagrangian Structures using:
% 	1. Springs
% 	2. Beams (*torsional springs)
% 	3. Target Points
%	4. Muscle-Model (combined Force-Length-Velocity model, "HIll+(Length-Tension)")
%
% One is able to update those Lagrangian Structure parameters, e.g., spring constants, resting %%	lengths, etc
%
% There are a number of built in Examples, mostly used for teaching purposes.
%
% If you would like us %to add a specific muscle model, please let Nick (nick.battista@unc.edu) know.
%
%--------------------------------------------------------------------------------------------------------------------%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the HEART-TUBE-EXAMPLE geometry and prints associated input files
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function HeartTube()

%
% Grid Parameters (MAKE SURE MATCHES IN input2d !!!)
%
Nx = 128;        % # of Eulerian Grid Pts. in x-Direction (MUST BE EVEN!!!)
Ny = 128;        % # of Eulerian Grid Pts. in y-Direction (MUST BE EVEN!!!)
Lx = 5.0;        % Length of Eulerian Grid in x-Direction
Ly = 5.0;        % Length of Eulerian Grid in y-Direction
dx = Lx/Nx;      % Eulerian Grid-step in x-direction
dy = Ly/Ny;      % Eulerian Grid-step in y-direction

% Immersed Structure Geometric / Dynamic Parameters %
ds = Lx/(2*Nx);  % Lagrangian Spacing
d = 1.0;         % Diameter of the tube
L = 3.0;         % Length of the heart-tube
struct_name = 'HeartTube'; % Name for .vertex, .spring, etc files.

% Call function to construct geometry
[xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,L,d,Lx,Ly);
yLag = yLag - 1.25;

fprintf('\nUSER-INPUT FOR MODEL IN FITZHUGH_NAGUMO_1d:\n');
fprintf('numPtsAlongTube = %d\n',length(xLag)/2);
fprintf('ind_Start = %d\n',length(xLag)-2+1);
fprintf('ind_End = %d\n\n',length(xLag)-2+length(xLag)/2)

% Plot Geometry to test
plot(xLag(1:end/2),yLag(1:end/2),'r-'); hold on;
plot(xLag(end/2+1:end),yLag(end/2+1:end),'r-'); hold on;
plot(xLag,yLag,'*'); hold on;
xlabel('x'); ylabel('y');
axis([0 Lx 0 Ly]);

% Prints .vertex file!
print_Lagrangian_Vertices(xLag,yLag,struct_name);

% Prints .spring file!
k_Spring = 1e7;
print_Lagrangian_Springs(xLag,k_Spring,ds,d,struct_name);

% Prints .muscle file!
%LFO = d; SK = 0.3; a = 0.25; b = 4.0; Fmax = 1e5;
%print_Lagrangian_Muscles(xLag,LFO,SK,a,b,Fmax,struct_name)

% Prints .beam file!
k_Beam = 2.5e6; C = 0.0;
print_Lagrangian_Beams(xLag,k_Beam,C,struct_name);

% Prints .target file!
k_Target = 1e7;
print_Lagrangian_Target_Pts(xLag,k_Target,struct_name);

% Prints .concentration file!
kDiffusion = 5e-5;
Concentration = give_Me_Initial_Concentration(Lx,Ly/2,Nx,Ny/2,dx,dy);
print_Concentration_Info(Nx,Ny/2,Concentration,kDiffusion,struct_name);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints VERTEX points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Vertices(xLag,yLag,struct_name)

    N = length(xLag);

    vertex_fid = fopen([struct_name '.vertex'], 'w');

    fprintf(vertex_fid, '%d\n', N );

    %Loops over all Lagrangian Pts.
    for s = 1:N
        X_v = xLag(s);
        Y_v = yLag(s);
        fprintf(vertex_fid, '%1.16e %1.16e\n', X_v, Y_v);
    end

    fclose(vertex_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints Vertex points to a file called rubberband.vertex
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Target_Pts(xLag,k_Target,struct_name)

    N = length(xLag); % Total Number of Lagrangian Pts
    num = 1;          % Number of target points on each end

    target_fid = fopen([struct_name '.target'], 'w');

    fprintf(target_fid, '%d\n', 4*num );

    %Loops over all Lagrangian Pts.
    %for s = 1:N

    %Left Bottom Target Points
    for s=1:num
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    %Right Bottom Target Points
    for s=(N/2)-(num-1):N/2
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    %Left Top Target Points
    for s=(N/2+1):(N/2+1)+(num-1)
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    %Right Top Target Points
    for s=N-(num-1):N
        fprintf(target_fid, '%d %1.16e\n', s, k_Target);
    end

    %end

    fclose(target_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints BEAM (Torsional Spring) points to a file called rubberband.beam
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Beams(xLag,k_Beam,C,struct_name)

    % k_Beam: beam stiffness
    % C: beam curvature

    N = length(xLag); % NOTE: Total number of beams = Number of Total Lag Pts. - 2

    beam_fid = fopen([struct_name '.beam'], 'w');

    fprintf(beam_fid, '%d\n', N - 4 );

    %spring_force = kappa_spring*ds/(ds^2);

    %BEAMS BETWEEN VERTICES
    for s = 2:N-1
            if  s <= N/2-1
                % Bottom of tube
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1, k_Beam, C);
            elseif ( ( s >= N/2+2 ) && ( s <= N-1 ) )
                % Top of Tube
                fprintf(beam_fid, '%d %d %d %1.16e %1.16e\n',s-1, s, s+1,   k_Beam, C);
            end
    end
    fclose(beam_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints MUSCLE points to a file called "struct_name".muscle
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Muscles(xLag,LFO,SK,a,b,Fmax,struct_name)

    N = length(xLag); %Number of Lagrangian Pts. Total

    muscle_fid = fopen([struct_name '.muscle'], 'w');

    fprintf(muscle_fid, '%d\n', N/2-2 );

    %spring_force = kappa_spring*ds/(ds^2);

    %MUSCLES BETWEEN VERTICES
    for s = 2:N/2-1
            fprintf(muscle_fid, '%d %d %1.16e %1.16e %1.16e %1.16e %1.16e\n', s, s+N/2, LFO, SK, a,b,Fmax);
    end
    fclose(muscle_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints SPRING points to a file called rubberband.spring
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Lagrangian_Springs(xLag,k_Spring,ds_Rest,d,struct_name)

    N = length(xLag); %Number of Lagrangian Pts. Total

    spring_fid = fopen([struct_name '.spring'], 'w');

    fprintf(spring_fid, '%d\n', N-2 + N/2 ); %(N-2) btwn adajcent Lag. Pts, (N/2) btwn opposite sides of HT

    %spring_force = kappa_spring*ds/(ds^2);

    %SPRINGS BETWEEN VERTICES
    for s = 1:N
            if s < N/2
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            elseif ( ( s >= N/2+1 ) && ( s < N ) )
                %Case s=N
                fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+1, k_Spring, ds_Rest);
            end
    end

    for s = 1:N/2
        fprintf(spring_fid, '%d %d %1.16e %1.16e\n', s, s+N/2,   k_Spring/1e10, 0.0);
    end
    fclose(spring_fid);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: creates the Lagrangian structure geometry
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [xLag,yLag] = give_Me_Immsersed_Boundary_Geometry(ds,L,d,Lx,Ly)

% The immsersed structure is a straight heart-tube %
x1 = [-L/2:ds:L/2 L/2];       % Constructs x-Values for bottom of tube
y1 = -d/2*ones(1,length(x1)); % Constructs y-Values for bottom of tube

x2 = x1;                      % Constructs x-Values for top of tube
y2 = -y1;                     % Constructs y-Values for top of tube

x1 = x1 + Lx/2;               % Shift into correct box
x2 = x2 + Lx/2;                % Shift into correct box

y1 = y1 + Ly/2;               % Shift into correct box
y2 = y2 + Ly/2;               % Shift into correct box

xLag = [x1 x2];
yLag = [y1 y2];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: prints CONCENTRATION INFO to file called
%           'struct_name'.concentration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function print_Concentration_Info(Nx,Ny,C,kDiffusion,struct_name)

    con_fid = fopen([struct_name '.concentration'], 'w');

    fprintf(con_fid, '%d\n', kDiffusion );

    for i=1:Ny
        for j=1:Nx
            fprintf(con_fid, '%1.16e ', C(i,j) );
        end
        fprintf(con_fid,'\n');
    end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: gives initial concentration gradient inside channel
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C = give_Me_Initial_Concentration(Lx,Ly,Nx,Ny,dx,dy)

%WHERE OUTER TUBE LIES
%xMin = 0.15; xMax = 0.45;
%yMin = 0.85; yMax = 1.15;

xMin = 1.6; xMax = 2.2;
yMin = 0.8; yMax = 1.7;

xMid = (xMin+xMax)/2;
yMid = (yMin+yMax)/2;

xDiff = (xMax-xMin)/2;
yDiff = (yMax-yMin)/2;

x = 0:dx:Lx;
y = 0:dy:Ly;
inds = give_Me_Indices_To_Apply_Force(x,y,xMin,xMax,yMin,yMax);

C = zeros(Ny,Nx);
for i=1:length( inds(:,1) )
    xInd = inds(i,1);
    yInd = inds(i,2);
    xPt = x(xInd);
    yPt = y(yInd);
    %C(xInd,yInd ) = (-0.5/yDiff^2)*( (yPt-yMid) - yDiff )*( (yPt-yMid) + yDiff ) +  (-0.5/xDiff^2)*( (xPt-xMid) - xDiff )*( (xPt-xMid) + xDiff ); %1.0;
    C(yInd,xInd ) = (-1.0/xDiff^2)*( (xPt-xMid) - xDiff )*( (xPt-xMid) + xDiff ); %1.0;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: computes indices for placing initial concentration
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function inds = give_Me_Indices_To_Apply_Force(x,y,xMin,xMax,yMin,yMax)

j=1; noMinYet = 1;
while noMinYet

    if ( x(j) >= xMin )
        iX_min = j;
        noMinYet = 0;
    end
    j=j+1;
end

j=length(x); noMaxYet = 1;
while noMaxYet

    if ( x(j) <= xMax )
        iX_max = j;
        noMaxYet = 0;
    end
    j=j-1;
end

j=1; noMinYet = 1;
while noMinYet

    if ( y(j) >= yMin )
        iY_min = j;
        noMinYet = 0;
    end
    j=j+1;
end

j=length(y); noMaxYet = 1;
while noMaxYet

    if ( y(j) <= yMax )
        iY_max = j;
        noMaxYet = 0;
    end
    j=j-1;
end

iX_Vec = iX_min:1:iX_max;
iY_Vec = iY_min:1:iY_max;

n = 1;
for i=1:length(iX_Vec)
    for j=1:length(iY_Vec)
        inds(n,1) = iX_Vec(i);
        inds(n,2) = iY_Vec(j);
        n = n+1;
    end
end

author: nickabattista


text:
function hydro = radiationIRF(hydro,tEnd,nDt,nDw,wMin,wMax)
% Calculates the normalized radiation impulse response function. This is
% equivalent to the radiation IRF in the theory section normalized by
% :math:`\rho`:
%
%     :math:`\overline{K}_{r,i,j}(t) = {\frac{2}{\pi}}\intop_0^{\infty}{\frac{B_{i,j}(\omega)}{\rho}}\cos({\omega}t)d\omega`
%
% Default parameters can be used by inputting [].
% See ``WEC-Sim\examples\BEMIO`` for examples of usage.
%
% Parameters
% ----------
%     hydro : struct
%         Structure of hydro data
%
%     tEnd : float
%         Calculation range for the IRF, where the IRF is calculated from t
%         = 0 to tEnd, and the default is 100 s
%
%     nDt : float
%         Number of time steps in the IRF, the default is 1001
%
%     nDw : float
%         Number of frequency steps used in the IRF calculation
%         (hydrodynamic coefficients are interpolated to correspond), the
%         default is 1001
%
%     wMin : float
%         Minimum frequency to use in the IRF calculation, the default is
%         the minimum frequency from the BEM data
%
%     wMax : float
%         Maximum frequency to use in the IRF calculation, the default is
%         the maximum frequency from the BEM data
%
% Returns
% -------
%     hydro : struct
%         Structure of hydro data with radiation IRF
%

p = waitbar(0,'Calculating radiation IRFs...');  % Progress bar

% Set defaults if empty
if isempty(tEnd)==1;  tEnd = 100;           end
if isempty(nDt)==1;    nDt = 1001;            end
if isempty(nDw)==1;    nDw = 1001;            end
if isempty(wMin)==1;  wMin = min(hydro.w);  end
if isempty(wMax)==1;  wMax = max(hydro.w);  end

% Interpolate to the given t and w
t = linspace(0,tEnd,nDt);
w = linspace(wMin,wMax,nDw);
N = sum(hydro.dof) * sum(hydro.dof);

% Calculate the impulse response function for radiation
n = 0;
hydro.ra_K = nan(sum(hydro.dof), sum(hydro.dof), length(t));
for i = 1:sum(hydro.dof)
    for j = 1:sum(hydro.dof)
        ra_B = interp1(hydro.w,squeeze(hydro.B(i,j,:)),w);
        hydro.ra_K(i,j,:) = (2/pi)*trapz(w,ra_B.*(cos(w.*t(:)).*w), 2);
        n = n+1;
    end
    waitbar(n/N)
end

% Calculate the infinite frequency added mass
ra_Ainf_temp = zeros(length(hydro.w),1);                                    %Initialize the variable
[~,F] = size(hydro);                                                        %Last data set in
if strcmp(hydro(F).code,'WAMIT')==0
    for i = 1:sum(hydro.dof)
        for j = 1:sum(hydro.dof)
            ra_A            = squeeze(hydro.A(i,j,:));
            ra_K            = squeeze(hydro.ra_K(i,j,:));
            for k = 1:length(hydro.w)                                       %Calculate the infinite frequency added mass at each input frequency
                ra_Ainf_temp(k,1)  = ra_A(k) + (1./hydro.w(k))*trapz(t,ra_K.*sin(hydro.w(k).*t.'));
            end
            hydro.Ainf(i,j) = mean(ra_Ainf_temp);                           %Take the mean across the vector of infinite frequency added mass
        end
    end
end

hydro.ra_t = t;
hydro.ra_w = w;
close(p)

end

author: WEC-Sim

repo: WEC-Sim

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


text:
%%
%
function [LT] = lraNTFbgmodel(T,lraNTF_opts)
  %% 4D TENSOR DEMO
  % inputPath = 'dataset/demo.avi';
  % [video nFrame vidWidth vidHeight] = load_input(inputPath);
  % T = im2double(video); size(T)

  %% 3D TENSOR DEMO
  % load('C:/ABSLibrary/dataset/trafficdb/trafficdb/traffic_patches.mat');
  % V = imgdb{100};
  % V = im2double(V);
  %T = tensor(V); size(T)

  %% Factorize subtensors with PARAFAC/CP algorithm with R = 1
  sT = size(T);
  blksize = 8;
  tdim = numel(size(T));

  if(tdim == 4)
    d = zeros(sT(4),prod(sT(1:2))/blksize^2);
    Tlen = size(T,4);
  end
  if(tdim == 3)
    d = zeros(sT(3),prod(sT(1:2))/blksize^2);
    Tlen = size(T,3);
  end

  kblk = 0;
  xy = zeros(1,2);
  %tic;
  displog('Factorizing subtensors with Parafac/CP algorithm');
  for xu = 1:8:sT(1)-7
    for yu = 1:8:sT(2)-7
      kblk = kblk+1;

      if(tdim == 4)
        Tblk = T(xu:xu+blksize-1, yu:yu+blksize-1,:,:);
        % show_4dvideo(Tblk);
        Yblk = permute(tensor(Tblk),[4 1 2 3]);
      end

      if(tdim == 3)
        Tblk = T(xu:xu+blksize-1, yu:yu+blksize-1,:);
        Yblk = permute(tensor(Tblk),[3 1 2]);
        % show_3dvideo(Tblk);
      end

      % lraNTF_opts = struct('NumOfComp',1,'CPalg','CP_HALS_opts.mat',...
      %  'lra_rank',2,'maxiter',500,'nlssolver','hals','initmode','cp','lra',true);
      [Ycap] = lraNTF(Yblk,lraNTF_opts);
      d(:,kblk) = double(Ycap.U{1});
      xy(kblk,:) = [xu yu];

      %[X_als,Uinit,out] = cp_als(Yblk,1);
      % X_als2 = double(permute(tensor(X_als),[2 3 4 1])); size(X_als2)
      % show_4dvideo(X_als2);
    end
  end
  %toc

  %% Find stationary blocks and build background model
  %tic;
  displog('Finding stationary blocks and build background model');
  maxd = max(d); mind = min(d);
  thresh = 0.005;

  if(tdim == 4) Imbgr = zeros(sT(1:3)); end
  if(tdim == 3) Imbgr = zeros(sT(1:2)); end

  for k = 1:size(d,2)
    edges = [mind(k):thresh:maxd(k) maxd(k)+eps];
    [n,bin] = histc(d(:,k),edges);
    m = mode(bin);
    indbgr = find((d(:,k) >= edges(m)) & (d(:,k) <= edges(m+1)));

    if(tdim == 4)
      bgrblk = median(T(xy(k,1):xy(k,1)+blksize-1, xy(k,2):xy(k,2)+blksize-1,:,indbgr),4);
      Imbgr(xy(k,1):xy(k,1)+blksize-1, xy(k,2):xy(k,2)+blksize-1,:) = bgrblk;
    end

    if(tdim == 3)
      bgrblk = median(double(T(xy(k,1):xy(k,1)+blksize-1, xy(k,2):xy(k,2)+blksize-1,indbgr)),3);
      Imbgr(xy(k,1):xy(k,1)+blksize-1, xy(k,2):xy(k,2)+blksize-1) = bgrblk;
    end
  end
  %toc
  displog('Finished');
  %%
  LT = repmat(Imbgr,[1 1 Tlen]);
end

author: andrewssobral

repo: mtt

save_path: github-repos/MATLAB/andrewssobral-mtt


text:
function [sr_img_deconv, sr_img, sampled_img, sr_affix] = lmsSR_generateSRusingNUIv4weighted(lr_seq, warped_meshXN, warped_meshYN, upscaling, psf, lucy_iter, weighting_struct)
%
% LMSSR_GENERATESRUSINGNUIV4WEIGHTED Performs a Non-Uniform Interpolation (NUI) Super-Resolution (SR) approach with Spatial Quantization using Cubic Interpolation
%                                    and a Triple Sample Weighting. Supports GS/Y and YUV.
%    [sr_img_deconv, sr_img, sampled_img, sr_affix] = LMSSR_GENERATESRUSINGNUIV4WEIGHTED(lr_seq, warped_meshXN, warped_meshYN, upscaling, psf, lucy_iter, weighting_struct)
%
% Parameters: lr_seq             -   4-D matrix containing a LR YUV444 or GS/Y sequence (height,width,dye,frames)
%             warped_meshXN      -   Warped image mesh for X-coordinates
%             warped_meshYN      -   Warped image mesh for Y-coordinates
%             upscaling          -   1x2 vector containing the upscaling factors for Y- and X-coordinates
%             psf                -   Point spread function
%             lucy_iter          -   Number of iteration for the Lucy-Richardson deconvolution
%             weighting_struct   -   Struct containing the selected weighting mode + corresponding extra information:
%
%                mec_weight_flag   -   MEC Weighting ON/OFF (0: OFF. 1: ON, 2: ON [MECv2], 3: ON [MECv3], 4: ON [MECv4])
%                dis_weight_flag   -   Distance Weighting ON/OFF
%                qps_weight_flag   -   QP Weighting ON/OFF
%
%                mec_confMapN      -   3-D matrix containing confidence information in form of SSD values for each pixel and each auxiliary frame
%                mec_thr_scaler    -   Scalar value responsible for scaling the standard deviation threshold which cuts off SSD values that are too large [default: 2]
%                mec_weight_exp    -   Exponent for intensifying the weighting [default: 2]
%                mec_weighting     -   3-D matrix containing weights for each pixel depending on the quality of the MEC step
%
%                dis_rho           -   Distance Weighting Decay Factor [default: 0.7]
%                dis_scaler        -   Distance Weighting Scaling Factor [default: upscaling*10]
%                dis_weighting     -   3-D matrix containing weights for each pixel depending on the distance to the nearest pixel center
%
%                qps_map           -   3-D matrix containing the quantization parameters for each pixel and each frame including the reference frame
%                qps_rho           -   Quantization Weighting Decay Factor [default: 0.7]
%                qps_weighting     -   3-D matrix containing weights for each pixel depending on the quantization parameters
%
%                alpha             -   MEC Weighting Mixing Factor [default: 1]
%                beta              -   Distance Weighting Mixing Factor [default: 1]
%                gamma             -   Quantization Weighting Mixing Factor [default: 1]
%
% Author: Michel Bätz (LMS)
%
% See also: lmsSR_framework
%

sr_affix = 'sr2weighted000';

if weighting_struct.mec_weight_flag == 1,
    sr_affix(end-2) = '1';
elseif weighting_struct.mec_weight_flag == 2,
    sr_affix(end-2) = '2';
elseif weighting_struct.mec_weight_flag == 3,
    sr_affix(end-2) = '3';
elseif weighting_struct.mec_weight_flag == 4,
    sr_affix(end-2) = '4';
end

if weighting_struct.dis_weight_flag == 1,
    sr_affix(end-1) = '1';
elseif weighting_struct.dis_weight_flag == 2,
    sr_affix(end-1) = '2';
end

if weighting_struct.qps_weight_flag,
    sr_affix(end) = '1';
end

disp(['Weighting-Mode [',sr_affix(end-2:end),'] Selected..']);

alpha = weighting_struct.alpha;
beta  = weighting_struct.beta;
gamma = weighting_struct.gamma;

% Extract information from lr_seq
[lr_height,lr_width,dye,num_imgs] = size(lr_seq);

% Build Original X-Y-Position Grid
[lr_meshX,lr_meshY] = meshgrid(1:lr_width,1:lr_height);

% Build SR X-Y-Position Grid (+Padding)
[sr_meshX,sr_meshY] = meshgrid(1:lr_width*upscaling(2)+2*upscaling(2),1:lr_height*upscaling(1)+2*upscaling(1));

% Computing MEC weights from SSD values
if weighting_struct.mec_weight_flag == 1,
    weighting_struct.mec_weighting = lmsSR_computeWeightsMEC(weighting_struct.mec_confMapN,weighting_struct.mec_thr_scaler,weighting_struct.mec_weight_exp);
elseif weighting_struct.mec_weight_flag == 2,
    weighting_struct.mec_weighting = lmsSR_computeWeightsMECv2(weighting_struct.mec_confMapN,weighting_struct.mec_thr_scaler,weighting_struct.mec_weight_exp);
elseif weighting_struct.mec_weight_flag == 3,
    weighting_struct.mec_weighting = lmsSR_computeWeightsMECv3(weighting_struct.mec_confMapN,weighting_struct.mec_thr_scaler,weighting_struct.mec_weight_exp);
elseif weighting_struct.mec_weight_flag == 4,
    weighting_struct.mec_weighting = lmsSR_computeWeightsMECv4(weighting_struct.mec_confMapN,weighting_struct.mec_thr_scaler*10,weighting_struct.mec_weight_exp); % MAGI-ADD: Scaling with 10 necessary!
else
    weighting_struct.mec_weighting = ones(lr_height,lr_width,num_imgs);
end

% Computing QP weights from quantization parameters
if weighting_struct.qps_weight_flag, % MAGI-TODO: Work-in-Progress
    weighting_struct.qps_weighting = lmsSR_computeWeightsQP(weighting_struct.qps_map,weighting_struct.qps_rho);
else
    weighting_struct.qps_weighting = ones(lr_height,lr_width,num_imgs);
end

% Linearized Meshes (Direct Computation)
merged_meshX = zeros(lr_height,lr_width,num_imgs);
merged_meshY = zeros(lr_height,lr_width,num_imgs);

merged_meshX(:,:,1)     = lr_meshX;
merged_meshX(:,:,2:end) = warped_meshXN;

merged_meshY(:,:,1)     = lr_meshY;
merged_meshY(:,:,2:end) = warped_meshYN;

merged_lin_meshX = merged_meshX(~isnan(merged_meshX));
merged_lin_meshY = merged_meshY(~isnan(merged_meshY));

merged_valsY = lr_seq(:,:,1,:);
merged_lin_valsY = merged_valsY(~isnan(merged_meshX));
if dye == 3,
    merged_valsU = lr_seq(:,:,2,:);
    merged_valsV = lr_seq(:,:,3,:);
    merged_lin_valsU = merged_valsU(~isnan(merged_meshX));
    merged_lin_valsV = merged_valsV(~isnan(merged_meshX));
end

% Removing NaN-Positions of the QP/MEC Weighting Matrix and Linearizing It
merged_mec_weighting = weighting_struct.mec_weighting(~isnan(merged_meshX));
merged_qps_weighting = weighting_struct.qps_weighting(~isnan(merged_meshX));

num_pix = numel(merged_lin_meshX);

disp('... Throwing points into spatially quantized buckets ...');
%timeSR2bucket = tic;

% Averaging all samples inside the HR grid pixel patch
rounded_meshX = round(merged_lin_meshX*upscaling(2)) - upscaling(2) + 1;
rounded_meshY = round(merged_lin_meshY*upscaling(1)) - upscaling(1) + 1;

rounding_errX = abs(merged_lin_meshX - (rounded_meshX + upscaling(2) - 1)/upscaling(2));
rounding_errY = abs(merged_lin_meshY - (rounded_meshY + upscaling(1) - 1)/upscaling(1));

% Memory Cleansing
clear warped_meshXN warped_meshYN lr_seq lr_meshX lr_meshY merged_meshX merged_meshY merged_valsY merged_valsU merged_valsV merged_lin_meshX merged_lin_meshY mec_weighting

% Throwing away all mesh points that lie outside the SR image grid
outlier_vec = (rounded_meshX < 1) | (rounded_meshX > lr_width*upscaling(2)) | (rounded_meshY < 1) | (rounded_meshY > lr_height*upscaling(1));

rounded_meshX = rounded_meshX(~outlier_vec);
rounded_meshY = rounded_meshY(~outlier_vec);

rounding_errX = rounding_errX(~outlier_vec);
rounding_errY = rounding_errY(~outlier_vec);

% Distance Weighting Stuff
if weighting_struct.dis_weight_flag == 1,
    weighting_struct.dis_weighting = weighting_struct.dis_rho.^(sqrt((rounding_errY*weighting_struct.dis_scaler(1)).^2 + (rounding_errX*weighting_struct.dis_scaler(1)).^2));
elseif weighting_struct.dis_weight_flag == 2,
    weighting_struct.dis_weighting = weighting_struct.dis_rho.^(sqrt((rounding_errY*weighting_struct.dis_scaler(1)/2).^2 + (rounding_errX*weighting_struct.dis_scaler(1)/2).^2));
else
    weighting_struct.dis_weighting = ones(size(rounding_errX));
end

merged_lin_valsY = merged_lin_valsY(~outlier_vec);

tmp_sr_imgY = zeros(lr_height*upscaling(1),lr_width*upscaling(2)); % Creating bucketed image
if dye == 3,
    merged_lin_valsU = merged_lin_valsU(~outlier_vec);
    merged_lin_valsV = merged_lin_valsV(~outlier_vec);

    tmp_sr_imgU = zeros(lr_height*upscaling(1),lr_width*upscaling(2)); % Creating bucketed image
    tmp_sr_imgV = zeros(lr_height*upscaling(1),lr_width*upscaling(2)); % Creating bucketed image
end

merged_mec_weighting = merged_mec_weighting(~outlier_vec);
merged_qps_weighting = merged_qps_weighting(~outlier_vec);

counter_matrix = zeros(lr_height*upscaling(1),lr_width*upscaling(2)); % Matrix for counting the number of entries in each bucket
weighting_matrix = zeros(lr_height*upscaling(1),lr_width*upscaling(2)); % Matrix for the sum of MEC-DIS-QPS weights

num_pix = num_pix - sum(outlier_vec);

if dye == 3, % Color case
    for itera = 1:num_pix,
        % Weighted case
        tmp_sr_imgY(rounded_meshY(itera),rounded_meshX(itera))    = tmp_sr_imgY(rounded_meshY(itera),rounded_meshX(itera)) + merged_lin_valsY(itera)*alpha*merged_mec_weighting(itera)*beta*weighting_struct.dis_weighting(itera)*gamma*merged_qps_weighting(itera);
        tmp_sr_imgU(rounded_meshY(itera),rounded_meshX(itera))    = tmp_sr_imgU(rounded_meshY(itera),rounded_meshX(itera)) + merged_lin_valsU(itera)*alpha*merged_mec_weighting(itera)*beta*weighting_struct.dis_weighting(itera)*gamma*merged_qps_weighting(itera);
        tmp_sr_imgV(rounded_meshY(itera),rounded_meshX(itera))    = tmp_sr_imgV(rounded_meshY(itera),rounded_meshX(itera)) + merged_lin_valsV(itera)*alpha*merged_mec_weighting(itera)*beta*weighting_struct.dis_weighting(itera)*gamma*merged_qps_weighting(itera);

        counter_matrix(rounded_meshY(itera),rounded_meshX(itera))   = counter_matrix(rounded_meshY(itera),rounded_meshX(itera)) + 1;
        weighting_matrix(rounded_meshY(itera),rounded_meshX(itera)) = weighting_matrix(rounded_meshY(itera),rounded_meshX(itera)) + alpha*merged_mec_weighting(itera)*beta*weighting_struct.dis_weighting(itera)*gamma*merged_qps_weighting(itera);
    end

    % Weighted case
    tmp_sr_imgY = tmp_sr_imgY ./ weighting_matrix;
    tmp_sr_imgU = tmp_sr_imgU ./ weighting_matrix;
    tmp_sr_imgV = tmp_sr_imgV ./ weighting_matrix;
else % Grayscale case
    for itera = 1:num_pix,
        % Weighted case
        tmp_sr_imgY(rounded_meshY(itera),rounded_meshX(itera))    = tmp_sr_imgY(rounded_meshY(itera),rounded_meshX(itera)) + merged_lin_valsY(itera)*alpha*merged_mec_weighting(itera)*beta*weighting_struct.dis_weighting(itera)*gamma*merged_qps_weighting(itera);

        counter_matrix(rounded_meshY(itera),rounded_meshX(itera))   = counter_matrix(rounded_meshY(itera),rounded_meshX(itera)) + 1;
        weighting_matrix(rounded_meshY(itera),rounded_meshX(itera)) = weighting_matrix(rounded_meshY(itera),rounded_meshX(itera)) + alpha*merged_mec_weighting(itera)*beta*weighting_struct.dis_weighting(itera)*gamma*merged_qps_weighting(itera);
    end

    % Weighted case
    tmp_sr_imgY = tmp_sr_imgY ./ weighting_matrix;
end

disp('... Throwing points into spatially quantized buckets: Done.');
%toc(timeSR2bucket);

sampled_img = tmp_sr_imgY;
if dye == 3,
    sampled_img(:,:,2) = tmp_sr_imgU;
    sampled_img(:,:,3) = tmp_sr_imgV;
end

disp('... Interpolating data to the HR grid ...');
%timeSR2interp = tic;

% For Cubic Interpolation a slight Padding is required to avoid incorrectly interpolated values near the image borders
tmp_sr_imgY = padarray(tmp_sr_imgY,upscaling,'symmetric','both');
if dye == 3,
    tmp_sr_imgU = padarray(tmp_sr_imgU,upscaling,'symmetric','both');
    tmp_sr_imgV = padarray(tmp_sr_imgV,upscaling,'symmetric','both');
end

% Bringing the bucketed image into a linearized form for the griddata call
bucket_meshX = sr_meshX(~isnan(tmp_sr_imgY));
bucket_meshY = sr_meshY(~isnan(tmp_sr_imgY));

bucket_valsY  = tmp_sr_imgY(~isnan(tmp_sr_imgY));
if dye == 3,
    bucket_valsU  = tmp_sr_imgU(~isnan(tmp_sr_imgY));
    bucket_valsV  = tmp_sr_imgV(~isnan(tmp_sr_imgY));
end

% SR Reconstruction using Non-Uniform Interpolation (NUI)
% Griddata-Solution (ScatteredInterpolant is called from inside the griddata function)
tmp_Y = griddata(bucket_meshX,bucket_meshY,bucket_valsY,sr_meshX,sr_meshY,'cubic');%,'linear');
tmp_Y(isnan(tmp_Y)) = 0; % Temporary NaN treatment for Y
sr_img(:,:,1) = tmp_Y;
if dye == 3, % YUV444 Input
  tmp_U = griddata(bucket_meshX,bucket_meshY,bucket_valsU,sr_meshX,sr_meshY,'cubic');%,'linear');
  tmp_V = griddata(bucket_meshX,bucket_meshY,bucket_valsV,sr_meshX,sr_meshY,'cubic');%,'linear');
  tmp_U(isnan(tmp_U)) = 127/255; % Temporary NaN treatment for U
  tmp_V(isnan(tmp_V)) = 127/255; % Temporary NaN treatment for V
  sr_img(:,:,2) = tmp_U;
  sr_img(:,:,3) = tmp_V;
end
%sr_affix = [sr_affix,'lin'];

% Cropping the Padding
sr_img = sr_img(1+upscaling(1):end-upscaling(1),1+upscaling(2):end-upscaling(2),:);

disp('... Interpolating data to the HR grid: Done.');
%toc(timeSR2interp);

disp('... Deblurring SR image ...');
%timeSR2deblur = tic;

% Padding Image for the Deblurring Step
pad_size = 16; % This was 5 prior

sr_img_tmp = padarray(sr_img,[pad_size pad_size],'replicate','both');

% SR Deblurring Step (Deconvolution)
sr_img_deconv = deconvlucy(sr_img_tmp,psf,lucy_iter); % Lucy-Richardson deconvolution for deblurring (more iterations improves the result for checkerboard)

% Remove the Pading after the Deblurring Step
sr_img_deconv = sr_img_deconv(1+pad_size:end-pad_size,1+pad_size:end-pad_size,:);

disp('... Deblurring SR image: Done.');
%toc(timeSR2deblur);


text:
function [pos, vel] = vectrs (ra, dec, pmra, pmdec, parllx, rv)

% this function converts angular quantities related to a star's
% position and motion to vectors.

% input

%  ra     = right ascension in hours

%  dec    = declination in degrees

%  pmra   = proper motion in ra in milliarcseconds per year

%  pmdec  = proper motion in dec in milliarcseconds per year

%  parllx = parallax in milliarcseconds

%  rv     = radial velocity in kilometers/second

% output

%  pos = position vector, equatorial rectangular coordinates,
%           with respect to solar system barycenter, components in au

%  vel = velocity vector, equatorial rectangular coordinates,
%        with respect to solar system barycenter, components in au/day

% ported from NOVAS 3.0

%%%%%%%%%%%%%%%%%%%%%%%

seccon = 180.0d0 * 3600.0d0 / pi;

% speed of light in kilometers/second

c = 1.0d-3 * 2997924580.0d0;

% au in kilometers

aukm = 1.0d-3 * 499.0047838061d0 * c;

% if parallax is unknown, undetermined, or zero, set it to 1e-6
% milliarcsecond, corresponding to a distance of 1 gigaparsec

paralx = parllx;

if ( paralx <= 0.0d0 )
    paralx = 1.0d-6;
end

% convert right ascension, declination, and parallax to position
% vector in equatorial system with units of au

dist = 1.0d0 / sin (paralx * 1.0d-3 / seccon);

r = ra * 54000.0d0 / seccon;

d = dec * 3600.0d0 / seccon;

cra = cos(r);

sra = sin(r);

cdc = cos(d);

sdc = sin(d);

pos(1) = dist * cdc * cra;

pos(2) = dist * cdc * sra;

pos(3) = dist * sdc;

% compute doppler factor, which accounts for change in
% light travel time to star

k = 1.d0 / (1.0d0 - rv / c);

% convert proper motion and radial velocity to orthogonal components
% of motion with units of au/day

pmr = pmra  / (paralx * 365.25d0) * k;

pmd = pmdec / (paralx * 365.25d0) * k;

rvl = rv * 86400.0d0 / aukm       * k;

% transform motion vector to equatorial system

vel(1) = - pmr * sra - pmd * sdc * cra + rvl * cdc * cra;

vel(2) =   pmr * cra - pmd * sdc * sra + rvl * cdc * sra;

vel(3) =               pmd * cdc       + rvl * sdc;

author: goGPS-Project


text:
function res = spm_eeg_regressors_tfpower(S)
% Generate regressors from power in TF dataset
% S                     - input structure
% fields of S:
%    S.D                - M/EEG object
%
%    Additional parameters can be defined specific for each plugin
% Output:
%  res -
%   If no input is provided the plugin returns a cfg branch for itself
%
%   If input is provided the plugin returns
%______________________________________________________________________________________
% Copyright (C) 2014 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak

SVNrev = '$Rev: 6188 $';

if nargin == 0
    %--------------------------------------------------------------------------
    % TF power dataset
    %--------------------------------------------------------------------------
    Dtf        = cfg_files;
    Dtf.tag    = 'Dtf';
    Dtf.name   = 'TF power dataset name';
    Dtf.filter = 'mat';
    Dtf.num    = [1 1];
    Dtf.help   = {'Select the M/EEG mat file containing TF power data'};

    %--------------------------------------------------------------------------
    % timewin
    %--------------------------------------------------------------------------
    timewin         = cfg_entry;
    timewin.tag     = 'timewin';
    timewin.name    = 'Time window';
    timewin.help    = {['Start and stop of the time window [ms]. ' ...
        '(Used only for summarised epoched cases)']};
    timewin.strtype = 'r';
    timewin.num     = [1 2];
    timewin.val     = {[-Inf Inf]};

    %--------------------------------------------------------------------------
    % freqwin
    %--------------------------------------------------------------------------
    freqwin         = cfg_entry;
    freqwin.tag     = 'freqwin';
    freqwin.name    = 'Frequency window';
    freqwin.help    = {'Start and stop of the frequency window (Hz).'};
    freqwin.strtype = 'r';
    freqwin.num     = [1 2];
    freqwin.val     = {[-Inf Inf]};

    %--------------------------------------------------------------------------
    % average
    %--------------------------------------------------------------------------
    average = cfg_menu;
    average.tag = 'average';
    average.name = 'Average over frequency';
    average.labels = {'Yes', 'No'};
    average.val = {1};
    average.values = {1,0};
    average.help = {'Average over the frequency window to produce one regressor',...
        'or output each frequency as a separate regressor'};

    %--------------------------------------------------------------------------
    % standardize
    %--------------------------------------------------------------------------
    standardize = cfg_menu;
    standardize.tag = 'standardize';
    standardize.name = 'Standardize';
    standardize.labels = {'Yes', 'No'};
    standardize.val = {0};
    standardize.values = {1,0};
    standardize.help = {'Standardize (zscore) regressor values'};

    %--------------------------------------------------------------------------
    % regname
    %--------------------------------------------------------------------------
    regname         = cfg_entry;
    regname.tag     = 'regname';
    regname.name    = 'Regressor name';
    regname.help    = {'Specify the string to be used as regressor name.'};
    regname.strtype = 's';
    regname.num     = [1 Inf];
    regname.val     = {'TFpower'};

    tfpower = cfg_branch;
    tfpower.tag = 'tfpower';
    tfpower.name = 'Time-frequency power';
    tfpower.val = {Dtf, spm_cfg_eeg_channel_selector, timewin, freqwin, average, standardize, regname};

    res = tfpower;

    return
end

%-Startup
%--------------------------------------------------------------------------
spm('sFnBanner', mfilename, SVNrev);
spm('FigName','Time-frequency power regressors');

if ~isfield(S, 'timewin'),      S.timewin  = [-Inf Inf];    end
if ~isfield(S, 'freqwin'),      S.freqwin  = [-Inf Inf];    end

if iscell(S.Dtf)
    S.Dtf = char(S.Dtf);
end

Dtf  = spm_eeg_load(S.Dtf);
D    = spm_eeg_load(S.D);

if ~isequal(Dtf.transformtype, 'TF');
    error('Time-frequency power dataset is expected as input.')
end

freqind = Dtf.indfrequency(min(S.freqwin)):Dtf.indfrequency(max(S.freqwin)); %BW
if isempty(freqind) || any(isnan(freqind))
    error('Selected frequency window is invalid.');
end

chanind = setdiff(Dtf.selectchannels(spm_cfg_eeg_channel_selector(S.channels)), Dtf.badchannels);

if isempty(chanind)
    error('No channels were selected');
end

if isequal(D.type, 'continuous')

    if ~isequal(Dtf.type, 'continuous') || (D.time(1) < Dtf.time(1)) || (D.time(end)>Dtf.time(end))
        error('All times of the input dataset should be within the power dataset.');
    end

    data = spm_squeeze(mean(Dtf(chanind, freqind, :), 1), [1 3]);

    if S.average
        data = mean(data, 1);
    end

    if D.fsample ~= Dtf.fsample
        [data, alpha] = spm_timeseries_resample(data, D.fsample/Dtf.fsample);
    else
        alpha = 1;
    end

    start = round(alpha*Dtf.indsample(D.time(1)));

    data = data(:, start:(start+D.nsamples-1));

else
    if D.ntrials ~= Dtf.ntrials
        error('Trial numbers should be equal between input and power dataset.');
    end

    if S.summarise
        timeind = D.indsample(1e-3*(min(S.timewin))):D.indsample(1e-3*(max(S.timewin)));
        if isempty(timeind) || any(isnan(timeind))
            error('Selected time window is invalid.');
        end
    else
        timeind = 1:D.nsamples;
    end

    data = spm_squeeze(mean(Dtf(chanind, freqind, timeind, :), 1), 1);

    if S.average
        data = mean(data, 1);
    end

    if S.summarise
        data = spm_squeeze(mean(data, 2), 2);
    else
        data = reshape(data, size(data, 1), []);
    end
end

data = data';

if S.standardize
    data = (data - repmat(mean(data, 1), size(data, 1), 1))./repmat(std(data, 1), size(data, 1), 1);
end

res.R     = data;

if S.average
    res.names = {S.regname};
else
    freq = Dtf.frequencies(freqind);
    for i = 1:length(freq)
        res.names{i} = sprintf('%s_%.1fHz', S.regname, freq(i));
    end
end

spm('FigName','Time-frequency power regressors: done');

author: spm