text:
function [ f,resL2,qualMeasOut] = AwPCSD(proj,geo,angles,maxiter,varargin)
%AwPCSD solves the reconstruction problem using adaptive-weighted
%projection-controlled steepest descent method
%
%   AwPCSD(PROJ,GEO,ALPHA,NITER) solves the reconstruction problem using
%   the projection data PROJ taken over ALPHA angles, corresponding to the
%   geometry described in GEO, using NITER iterations.
%
%   AwPCSD(PROJ,GEO,ALPHA,NITER,OPT,VAL,...) uses options and values for solving. The
%   possible options in OPT are:
%
%
%   'lambda':      Sets the value of the hyperparameter for the SART iterations.
%                  Default is 1
%
%   'lambdared':   Reduction of lambda.Every iteration
%                  lambda=lambdared*lambda. Default is 0.99
%
%       'init':    Describes diferent initialization techniques.
%                   •  'none'     : Initializes the image to zeros(default)

%                   •  'FDK'      : intializes image to FDK reconstrucition
%
%   'TViter':      Defines the amount of TV iterations performed per SART
%                  iteration. Default is 20
%
%   'maxL2err'     Maximum L2 error to accept an image as valid. This
%                  parameter is crucial for the algorithm, determines at
%                  what point an image should not be updated further.
%                  Default is 20% of the FDK L2 norm.
%   'Verbose'      1 or 0. Default is 1. Gives information about the
%                  progress of the algorithm.
%   'delta'        Defines Parameter to control the amount of smoothing
%                  for pixels at the edges. A large 'delta' is not able to
%                  differentiate image gradients at different pixels. A
%                  small 'delta' give low weights to almost every pixels,
%                  making the algorithm inefficient in removing noise or
%                  straking artifacts. Default is -0.00055
% 'redundancy_weighting': true or false. Default is true. Applies data
%                         redundancy weighting to projections in the update step
%                         (relevant for offset detector geometry)
%  'groundTruth'  an image as grounf truth, to be used if quality measures
%                 are requested, to plot their change w.r.t. this known
%                 data.
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
% This file is part of the TIGRE Toolbox
%
% Copyright (c) 2015, University of Bath and
%                     CERN-European Organization for Nuclear Research
%                     All rights reserved.
%
% License:            Open Source under BSD.
%                     See the full license at
%                     https://github.com/CERN/TIGRE/blob/master/LICENSE
%
% Contact:            tigre.toolbox@gmail.com
% Codes:              https://github.com/CERN/TIGRE/
% Coded by:           Ander Biguri and Manasavee Lohvithee
%--------------------------------------------------------------------------
%% parse inputs
[beta,beta_red,f,ng,verbose,epsilon,delta,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,varargin);
measurequality=~isempty(QualMeasOpts) | ~any(isnan(gt(:)));
if ~any(isnan(gt(:)))
    QualMeasOpts{end+1}='error_norm';
    res_prev=gt;
    clear gt
end
if nargout<3 && measurequality
    warning("Image metrics requested but none catched as output. Call the algorithm with 3 outputs to store them")
    measurequality=false;
end
qualMeasOut=zeros(length(QualMeasOpts),niter);

if nargout>1
    computeL2=true;
else
    computeL2=false;
end
resL2=zeros(1,niter);

% [alphablocks,orig_index]=order_subsets(angles,blocksize,OrderStrategy);
%
% angles=cell2mat(alphablocks);
% index_angles=cell2mat(orig_index);

% does detector rotation exists?
if ~isfield(geo,'rotDetector')
    geo.rotDetector=[0;0;0];
end

%% Create weigthing matrices for the SART step
% the reason we do this, instead of calling the SART fucntion is not to
% recompute the weigths every AwASD-POCS iteration, thus effectively doubling
% the computational time
% Projection weigth, W
W=computeW(geo,angles,gpuids);

% Back-Projection weigth, V
V=computeV(geo,angles,num2cell(angles),num2cell(1:length(angles)),'gpuids',gpuids);

if redundancy_weights
    % Data redundancy weighting, W_r implemented using Wang weighting
    % reference: https://iopscience.iop.org/article/10.1088/1361-6560/ac16bc

    num_frames = size(proj,3);
    W_r = redundancy_weighting(geo);
    W_r = repmat(W_r,[1,1,num_frames]);
    % disp('Size of redundancy weighting matrix');
    % disp(size(W_r));
    W = W.*W_r; % include redundancy weighting in W
end
%Initialize image.
%f=zeros(geo.nVoxel','single');

iter=0;
offOrigin=geo.offOrigin;
offDetector=geo.offDetector;
rotDetector=geo.rotDetector;
stop_criteria=0;
DSD=geo.DSD;
DSO=geo.DSO;
while ~stop_criteria %POCS
    % If quality is going to be measured, then we need to save previous image
    if measurequality && ~strcmp(QualMeasOpts,'error_norm')
        res_prev = f; % only store if necesary
    end
    f0=f;
    if (iter==0 && verbose==1);tic;end
    iter=iter+1;

    %Estimation error in the projection domain
    est_proj=Ax(f,geo,angles,'interpolated');
    delta_p=im3Dnorm(est_proj-proj,'L2');

    %Enforcing ART along all projections if squared delta_p > epsilon
    if (delta_p^2)>epsilon
        for jj=1:size(angles,2)
            if size(offOrigin,2)==size(angles,2)
                geo.offOrigin=offOrigin(:,jj);
            end
            if size(offDetector,2)==size(angles,2)
                geo.offDetector=offDetector(:,jj);
            end
            if size(rotDetector,2)==size(angles,2)
                geo.rotDetector=rotDetector(:,jj);
            end
            if size(DSD,2)==size(angles,2)
                geo.DSD=DSD(jj);
            end
            if size(DSO,2)==size(angles,2)
                geo.DSO=DSO(jj);
            end

            f=f+beta* bsxfun(@times,1./sum(V(:,:,jj),3),Atb(W(:,:,jj).*(proj(:,:,jj)-Ax(f,geo,angles(:,jj),'gpuids',gpuids)),geo,angles(:,jj),'gpuids',gpuids));

        end
    end

    %Non-negativity projection on all pixels
    if nonneg
        f=max(f,0);
    end

    geo.offDetector=offDetector;
    geo.offOrigin=offOrigin;

    if measurequality
        qualMeasOut(:,iter)=Measure_Quality(res_prev,f,QualMeasOpts);
    end

    % Compute L2 error of actual image. Ax-b
    dd=im3Dnorm(Ax(f,geo,angles,'gpuids',gpuids)-proj,'L2');

    % Compute change in the image after last SART iteration
    dp_vec=(f-f0);

    if iter==1
        step=1;
    else
        step=delta_p/delta_p_first;
    end
    f0=f;
    %  TV MINIMIZATION
    % =========================================================================
    %  Call GPU to minimize AwTV
    f=minimizeAwTV(f0,step,ng,delta,'gpuids',gpuids);    %   This is the MATLAB CODE, the functions are sill in the library, but CUDA is used nowadays
    %                                             for ii=1:ng
    %                                                 %delta=-0.00038 for thorax phantom
    %                                                 df=weighted_gradientTVnorm(f,delta);
    %                                                 df=df./im3Dnorm(df,'L2');
    %                                                 f=f-(step.*df);
    %                                             end
    % Compute change by TV min
    dg_vec=(f-f0);

    if iter==1
        delta_p_first=im3Dnorm((Ax(f0,geo,angles,'interpolated','gpuids',gpuids))-proj,'L2');
    end

    % Reduce SART step
    beta=beta*beta_red;

    % Check convergence criteria
    % ==========================================================================

    %Define c_alpha as in equation 21 in the journal
    c=dot(dg_vec(:),dp_vec(:))/(norm(dg_vec(:),2)*norm(dp_vec(:),2));
    %This c is examined to see if it is close to -1.0

    if (c<-0.99 && dd<=epsilon) || beta<0.005|| iter>maxiter
        if verbose
            disp(['Stopping criteria met']);
            disp(['   c    = ' num2str(c)]);
            disp(['   beta = ' num2str(beta)]);
            disp(['   iter = ' num2str(iter)]);
        end
        stop_criteria=true;
    end

    if computeL2
        geo.offOrigin=offOrigin;
        geo.offDetector=offDetector;
        errornow=im3Dnorm(proj-Ax(f,geo,angles,'gpuids',gpuids),'L2');                       % Compute error norm2 of b-Ax
        % If the error is not minimized.
        if  iter~=1 && errornow>errorL2(end)
            if verbose
                disp(['Convergence criteria met, exiting on iteration number:', num2str(iter)]);
            end
            return;
        end
        errorL2=[errorL2 errornow];
    end

    if (iter==1 && verbose==1)
        expected_time=toc*maxiter;
        disp('AwPCSD');
        disp(['Expected duration  :    ',secs2hms(expected_time)]);
        disp(['Expected finish time:    ',datestr(datetime('now')+seconds(expected_time))]);
        disp('');
    end

end
end

function [beta,beta_red,f0,ng,verbose,epsilon,delta,QualMeasOpts,nonneg,gpuids,redundancy_weights,gt]=parse_inputs(proj,geo,angles,argin)
opts=     {'lambda','lambda_red','init','tviter','verbose','maxl2err','delta','qualmeas','nonneg','gpuids','redundancy_weighting','groundtruth'};
defaults=ones(length(opts),1);
% Check inputs
nVarargs = length(argin);
if mod(nVarargs,2)
    error('TIGRE:AwPCSD:InvalidInput','Invalid number of inputs')
end

% check if option has been passed as input
for ii=1:2:nVarargs
    ind=find(ismember(opts,lower(argin{ii})));
    if ~isempty(ind)
        defaults(ind)=0;
    else
        error('TIGRE:AwPCSD:InvalidInput',['Optional parameter "' argin{ii} '" does not exist' ]);
    end
end

for ii=1:length(opts)
    opt=opts{ii};
    default=defaults(ii);
    % if one option isnot default, then extract value from input
    if default==0
        ind=double.empty(0,1);jj=1;
        while isempty(ind)
            ind=find(isequal(opt,lower(argin{jj})));
            jj=jj+1;
        end
        if isempty(ind)
            error('TIGRE:AwPCSD:InvalidInput',['Optional parameter "' argin{jj} '" does not exist' ]);
        end
        val=argin{jj};
    end
    % parse inputs
    switch opt
        % Verbose
        %  =========================================================================
        case 'verbose'
            if default
                verbose=1;
            else
                verbose=val;
            end
            if ~is2014bOrNewer
                warning('TIGRE:Verbose mode not available for older versions than MATLAB R2014b');
                verbose=false;
            end
            % Lambda
            %  =========================================================================
        case 'lambda'
            if default
                beta=1;
            else
                if length(val)>1 || ~isnumeric( val)
                    error('TIGRE:AwPCSD:InvalidInput','Invalid lambda')
                end
                beta=val;
            end
            % Lambda reduction
            %  =========================================================================
        case 'lambda_red'
            if default
                beta_red=0.99;
            else
                if length(val)>1 || ~isnumeric( val)
                    error('TIGRE:AwPCSD:InvalidInput','Invalid lambda')
                end
                beta_red=val;
            end
            % Initial image
            %  =========================================================================
        case 'init'
            if default || strcmp(val,'none')
                f0=zeros(geo.nVoxel','single');

            else
                if strcmp(val,'FDK')
                    f0=FDK(proj, geo, angles);
                else
                    error('TIGRE:AwPCSD:InvalidInput','Invalid init')
                end
            end
            % Number of iterations of TV
            %  =========================================================================
        case 'tviter'
            if default
                ng=20;
            else
                ng=val;
            end
            %  Maximum L2 error to have a "good image"
            %  =========================================================================
        case 'maxl2err'
            if default
                epsilon=im3Dnorm(FDK(proj,geo,angles))*0.2; %heuristic
            else
                epsilon=val;
            end
            %Parameter to control the amount of smoothing for pixels at the
            %edges
            %  =========================================================================
        case 'delta'
            if default
                delta=-0.005;
            else
                delta=val;
            end
            % Image Quality Measure
            %  =========================================================================
        case 'qualmeas'
            if default
                QualMeasOpts={};
            else
                if iscellstr(val)
                    QualMeasOpts=val;
                else
                    error('TIGRE:AwPCSD:InvalidInput','Invalid quality measurement parameters');
                end
            end
            %  Non negative
            %  =========================================================================
        case 'nonneg'
            if default
                nonneg=true;
            else
                nonneg=val;
            end
            %  GPU Ids
            %  =========================================================================
        case 'gpuids'
            if default
                gpuids = GpuIds();
            else
                gpuids = val;
            end
        case 'redundancy_weighting'
            if default
                redundancy_weights = true;
            else
                redundancy_weights = val;
            end
                case 'groundtruth'
            if default
                gt=nan;
            else
                gt=val;
            end
        otherwise
            error('TIGRE:AwPCSD:InvalidInput',['Invalid input name:', num2str(opt),'\n No such option in AwPCSD()']);

    end
end

end

author: CERN

repo: TIGRE

save_path: github-repos/MATLAB/CERN-TIGRE


text:
clear all; close all; clc;

% The matrix to be displayed.
data = spiral(100);

spx.graphics.display.matrix(data);

data = randn(100, 100);
options.color_map = 'gray';
spx.graphics.display.matrix(data, options);

author: indigits

repo: sparse-plex


text:
% Test for bezier curves

N = 1024;
cont = true;
hold on;
axis([0 1 0 1]);
C = [];
while cont
    [x,y,b] = ginput(1);
    cont = b==1;
    if cont
        plot(x,y, '.');
        axis([0 1 0 1]);
        C(:,end+1) = [x;y];
    end
end
compute_bezier_curve(C,N);
hold off;

author: gpeyre

repo: matlab-toolboxes


text:
function geometry_test043 ( )

%*****************************************************************************80
%
%% TEST043 tests SEGMENTS_DIST_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 December 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST043\n' );
  fprintf ( 1, '  SEGMENTS_DIST_3D computes the distance between\n' );
  fprintf ( 1, '    line segments in 3D.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Case   Computed    True\n' );
  fprintf ( 1, '\n' );
%
%  Case 1, parallel, not coincident.
%
%  LS1: (2,3,0) + t * (2,1,0) for t = 0 to 3.
%  LS2: (11,6,4) + t * (2,1,0) for t = 0 to 3.
%  Distance is 5.
%
  p1(1:3,1) = [  2.0; 3.0; 0.0 ];
  p2(1:3,1) = [  8.0; 6.0; 0.0 ];
  q1(1:3,1) = [ 11.0; 6.0; 4.0 ];
  q2(1:3,1) = [ 17.0; 9.0; 4.0 ];

  dist = segments_dist_3d ( p1, p2, q1, q2 );

  fprintf ( 1, '  1  %8f  5.0\n', dist );
%
%  Case 2, parallel, coincident, overlapping.
%
%  (1,2,3) + t * ( 1,-1,2)
%  LS1: t = 0 to t = 3.
%  Distance is 0.
%
  p1(1:3,1) = [  1.0;  2.0;  3.0 ];
  p2(1:3,1) = [  4.0; -1.0;  9.0 ];
  q1(1:3,1) = [  3.0;  0.0;  7.0 ];
  q2(1:3,1) = [  6.0; -3.0; 13.0 ];

  dist = segments_dist_3d ( p1, p2, q1, q2 );

  fprintf ( 2, '  1  %8f  0.0\n', dist );
%
%  Case 3, parallel, coincident, disjoint.
%
%  LS1: (3,4,5) + t * ( 2,2,1) for 0 <= t <= 2.
%  LS2: (3,4,5) + t * ( 2,2,1) for 3 <= t <= 5.
%  Distance = 3.
%
  p1(1:3,1) = [  3.0;  4.0;  5.0 ];
  p2(1:3,1) = [  7.0;  8.0;  7.0 ];
  q1(1:3,1) = [  9.0; 10.0;  8.0 ];
  q2(1:3,1) = [ 13.0; 14.0; 10.0 ];

  dist = segments_dist_3d ( p1, p2, q1, q2 );

  fprintf ( 1, '  3  %8f  3.0\n', dist );
%
%  Case 4, nonparallel, could intersect, and does intersect.
%
%  L1: (1,1,1) + t * (0,1,2)
%  L2: (0,2,3) + t * (1,0,0)
%  intersect at (1,2,3)
%  Distance is 0.
%
  p1(1:3,1) = [  1.0;  1.0;  1.0 ];
  p2(1:3,1) = [  1.0;  4.0;  7.0 ];
  q1(1:3,1) = [  0.0;  2.0;  3.0 ];
  q2(1:3,1) = [  5.0;  2.0;  3.0 ];

  dist = segments_dist_3d ( p1, p2, q1, q2 );

  fprintf ( 1, '  4  %8f  0.0\n', dist );
%
%  Case 5, nonparallel, could intersect, and does not intersect.
%
%  L1: (1,1,1) + t * (0,1,2)
%  L2: (0,2,3) + t * (1,0,0)
%  lines intersect at (1,2,3), line segments do not.
%  Distance is 1.0
%
  p1(1:3,1) = [  1.0;  1.0;  1.0 ];
  p2(1:3,1) = [  1.0;  4.0;  7.0 ];
  q1(1:3,1) = [  0.0;  2.0;  3.0 ];
  q2(1:3,1) = [ -5.0;  2.0;  3.0 ];

  dist = segments_dist_3d ( p1, p2, q1, q2 );

  fprintf ( 1, '  5  %8f  1.0\n', dist );
%
%  Case 6, nonparallel, can not intersect, "end-to-end".
%
%  L1: (2,2,1) + t * (0,1,2)  0 <= t <= 5
%  L2: (0,0,0) + t * (-1,-1,-1) 0 <= t <= 5
%  Distance is 3.
%
  p1(1:3,1) = [  2.0;  2.0;  1.0 ];
  p2(1:3,1) = [  2.0;  7.0;  11.0 ];
  q1(1:3,1) = [  0.0;  0.0;  0.0 ];
  q2(1:3,1) = [ -5.0; -5.0; -5.0 ];

  dist = segments_dist_3d ( p1, p2, q1, q2 );

  fprintf ( 1, '  6  %8f  3.0\n', dist );
%
%  Case 7, nonparallel, can not intersect, "end-to-mid".
%
%  L1: (1,1,1) + t * (0,1,2) 0 <= t <= 5
%  L2: (0,4,7) + t * (-1,0,0) 0 <= t <= 5
%  Distance is 1.
%
  p1(1:3,1) = [  1.0;  1.0;  1.0 ];
  p2(1:3,1) = [  1.0;  6.0; 11.0 ];
  q1(1:3,1) = [  0.0;  4.0;  7.0 ];
  q2(1:3,1) = [ -5.0;  4.0;  7.0 ];

  dist = segments_dist_3d ( p1, p2, q1, q2 );

  fprintf ( 1, '  7  %8f  1.0\n', dist );
%
%  Case 8, nonparallel, can not intersect, "mid-to-mid".
%
%  L1: (0,5,10) + t * (1,-1,0) 0 <= t <= 5
%  L2: (0,0,0) + t * (1,1,0) 0 <= t <= 6
%  Distance = 10.
%
  p1(1:3,1) = [  0.0;  5.0; 10.0 ];
  p2(1:3,1) = [  5.0;  0.0; 10.0 ];
  q1(1:3,1) = [  0.0;  0.0;  0.0 ];
  q2(1:3,1) = [  6.0;  6.0;  0.0 ];

  dist = segments_dist_3d ( p1, p2, q1, q2 );

  fprintf ( 1, '  8  %8f 10.0\n', dist );
%
%  Case 9, nonparallel, can not intersect, "mid-to-end".
%
%  L1: (-2,0,0) + t * (1,0,0) 0 <= t <= 12
%  L2: (-2,8,1) + t * (9,-4,-1) 0 <= t <= 1
%  Distance = 4.
%
  p1(1:3,1) = [ -2.0;  0.0;  0.0 ];
  p2(1:3,1) = [ 10.0;  0.0;  0.0 ];
  q1(1:3,1) = [ -2.0;  8.0;  1.0 ];
  q2(1:3,1) = [  7.0;  4.0;  0.0 ];

  dist = segments_dist_3d ( p1, p2, q1, q2 );

  fprintf ( 1, '  9  %8f  4.0\n', dist );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function f = cumsum(f, m, dim)
%CUMSUM   Indefinite integral of a TRIGTECH.
%   CUMSUM(F) is the indefinite integral of the TRIGTECH F, whose mean
%   is zero, with the constant of integration chosen so that F(-1) = 0.
%   If the mean of F is not zero then an error is thrown since the indefinite
%   integral would no longer be periodic.
%
%   CUMSUM(F, M) will compute the Mth definite integral with the constant of
%   integration chosen so that each intermediary integral evaluates to 0 at
%   -1.
%   Thus, CUMSUM(F, 2) is equivalent to CUMSUM(CUMSUM(F)).
%
%   CUMSUM(F, M, 2) will take the Mth cumulative sum over the columns F an
%   array-valued TRIGTECH.
%
% See also DIFF, SUM.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If the TRIGTECH G of length n is represented as
%       \sum_{k=-(n-1)/2}^{(n-1)/2} c_k exp(i*pi*kx)
% its integral is represented with a TRIGTECH of length n given by
%       \sum_{k=-(n-1)/2}^{(n-1)/2} b_k exp(i*pi*kx)
% where b_0 is determined from the constant of integration as
%       b_0 = \sum_{k=-(n-1)/2}^{(n-1)/2} (-1)^k/(i*pi*k) c_k;
% with c_0 := 0. The other coefficients are given by
%       b_k = c_k/(i*pi*k).
%
% If the TRIGTECH G of length n is represented as
%       \sum_{k=-n/2+1}^{n/2-1} c_k exp(i*pi*kx) + c(n/2)cos(n*pi/2x)
% then first set c(n) = 0.5*c(n) and define a = [0.5*c(n/2) c] so that we
% have the equivalent expansion:
%       \sum_{k=-n/2}^{n/2} a_k exp(i*pi*kx)
% The integral of this is represented with a TRIGTECH of length n+1 given by
%       \sum_{k=-n/2}^{n/2} b_k exp(i*pi*kx)
% where b_0 is determined from the constant of integration as
%       b_0 = \sum_{k=-n/2}^{n/2} (-1)^k/(i*pi*k) a_k;
% with a_0 := 0. The other coefficients are given by
%       b_k = a_k/(ik).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Trivial case of an empty TRIGTECH:
if ( isempty(f) )
    return
end

if ( nargin < 2 || isempty(m) )
    % Order of integration not passed in. Assume 1 by default:
    m = 1;
elseif ( m == 0 )
    % Nothing to do here!
    return
end

% Sum with respect to the continuous variable by default:
if ( nargin < 3 )
    dim = 1;
end

if ( dim == 1 )
    % Take difference across 1st dimension:
    f = cumsumContinuousDim(f, m);
else
    % Take difference across 2nd dimension:
    f = cumsumFiniteDim(f, m);
end

end

function f = cumsumContinuousDim(f, m)
% CUMSUM over the continuous dimension.

    % Initialize storage:
    c = f.coeffs; % Obtain Fourier coefficients {c_k}
    numCoeffs = size(c,1);

    fIsEven = mod(numCoeffs, 2) == 0;

    % index of constant coefficient
    if fIsEven
       ind = numCoeffs/2 + 1;
    else
       ind = (numCoeffs + 1)/2;
    end

    % Check that the mean of the TRIGtech is zero.  If it is not, then
    % throw an error.
    if ( any(abs(c(ind,:)) > 1e1*vscale(f)*eps) )
        error('CHEBFUN:TRIGTECH:cumsum:meanNotZero', ...
            ['Indefinite integrals are only possible for TRIGTECH objects '...
            'with zero mean.']);
    end

    % Force the mean to be exactly zero.
    if ( fIsEven )
        % Set coeff corresponding to the constant mode to zero:
        c(numCoeffs/2+1,:) = 0;
        % Expand the coefficients to be symmetric (see above discussion).
        c(1,:) = 0.5*c(1,:);
        c = [c; c(1,:)];
        highestDegree = numCoeffs/2;
    else
        c((numCoeffs+1)/2,:) = 0;
        highestDegree = (numCoeffs-1)/2;
    end

    % Loop for integration factor for each coefficient:
    sumIndicies = (-highestDegree:highestDegree).';
    integrationFactor = (-1i./sumIndicies/pi).^m;
    % Zero out the one corresponding to the constant term.
    integrationFactor(highestDegree+1) = 0;
    c = bsxfun(@times,c,integrationFactor);
    % If this is an odd order cumsum and there are an even number of
    % coefficients then zero out the cofficient corresponding to sin(N/2x)
    % term, since this will be zero on the Fourier grid.
    if ( (mod(m, 2) == 1) && fIsEven )
        c(1,:) = 0;
        c(numCoeffs+1,:) = 0;
    end

    % Fix the constant term.
    c(highestDegree+1,:) = -sum(bsxfun(@times,c,(-1).^sumIndicies));

    % If the original TRIGTECH had an even number of coefficients then
    % shrink the coefficent vector corresponding to its indefinite integral
    % back to its original size since it was increased by one above to make
    % the integration code slicker.
    if ( fIsEven )
        c = c(1:end-1,:);
    end

    % Recover values and attach to output:
    f.values = f.coeffs2vals(c);
    f.values(:,f.isReal) = real(f.values(:,f.isReal));

    f.coeffs = c;

    % Simplify (as suggested in Chebfun ticket #128)
    f = simplify(f);

    % Ensure f(-1) = 0:
    lval = get(f, 'lval');
    f.coeffs(1,:) = f.coeffs(1,:) - lval;
    f.values = bsxfun(@minus, f.values, lval);

end

function f = cumsumFiniteDim(f, m)
% CUMSUM over the finite dimension.

    for k = 1:m
        f.values = cumsum(f.values, 2);
        f.coeffs = cumsum(f.coeffs, 2);
    end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function displaytable(colname,varargin)
%% DISPLAYTABLE display data in table
%
% Example
%
%	load dispdata
%   colname = {'#nodes'     '|u_I-u_h|_1'      '||p-p_h||'};
%   displaytable(colname,err.N,[],err.uH1,'%0.5e',err.pL2,'%0.5e');
%
% Note that all data are column arrays.
%
% See also disptable
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

col = varargin{1};
rownum = size(col,1);
colnum = size(colname,2);
colwidth = zeros(1,colnum);
dispresult = char(zeros(rownum,6*colnum));
idx = 1;
for c = 1:nargin/2
    if isempty(varargin{2*c})
        currentstr = num2str(varargin{2*c-1});
    else
        currentstr = num2str(varargin{2*c-1},varargin{2*c});
    end
    colwidth(c) = size(currentstr,2);
    dispresult(:,idx:idx+colwidth(c)-1) = currentstr;
    idx = idx+colwidth(c)+3;
end
for c = 1:colnum
    spnum = colwidth(c) - length(colname{c})+3;
    if spnum >0
        if c>1
            colname{c} = [char(zeros(1,floor(spnum/2))) colname{c} ...
                          char(zeros(1,ceil(spnum/2)))];
        else
            colname{c} = [char(zeros(1,floor(spnum/3))) colname{c} ...
                          char(zeros(1,ceil(spnum/2)))];
        end
    end
end
% fprintf('\n');
disp(cell2mat(colname));
fprintf('\n');
disp(dispresult);
fprintf('\n');

repo: ifem


text:
%  When you load any ANALYZE or NIfTI file with 'load_nii.m', and view
%  it with 'view_nii.m', you may find that the image is L-R flipped.
%  This is because of the confusion of radiological and neurological
%  convention in the medical image before NIfTI format is adopted. You
%  can find more details from:
%
%  http://www.rotman-baycrest.on.ca/~jimmy/UseANALYZE.htm
%
%  Sometime, people even want to convert RAS (standard orientation) back
%  to LAS orientation to satisfy the legend programs or processes. This
%  program is only written for those purpose. So PLEASE BE VERY CAUTIOUS
%  WHEN USING THIS 'FLIP_LR.M' PROGRAM.
%
%  With 'flip_lr.m', you can convert any ANALYZE or NIfTI (no matter
%  3D or 4D) file to a flipped NIfTI file. This is implemented simply
%  by flipping the affine matrix in the NIfTI header. Since the L-R
%  orientation is determined there, so the image will be flipped.
%
%  Usage: flip_lr(original_fn, flipped_fn, [old_RGB],[tolerance],[preferredForm])
%
%  original_fn  -  filename of the original ANALYZE or NIfTI (3D or 4D) file
%
%  flipped_fn  -  filename of the L-R flipped NIfTI file
%
%  old_RGB (optional)  -  a scale number to tell difference of new RGB24
%	from old RGB24. New RGB24 uses RGB triple sequentially for each
%	voxel, like [R1 G1 B1 R2 G2 B2 ...]. Analyze 6.0 from AnalyzeDirect
%	uses old RGB24, in a way like [R1 R2 ... G1 G2 ... B1 B2 ...] for
%	each slices. If the image that you view is garbled, try to set
%	old_RGB variable to 1 and try again, because it could be in
%	old RGB24. It will be set to 0, if it is default or empty.
%
%  tolerance (optional) - distortion allowed for non-orthogonal rotation
%	or shearing in NIfTI affine matrix. It will be set to 0.1 (10%),
%	if it is default or empty.
%
%  preferredForm (optional)  -  selects which transformation from voxels
%	to RAS coordinates; values are s,q,S,Q.  Lower case s,q indicate
%	"prefer sform or qform, but use others if preferred not present".
%	Upper case indicate the program is forced to use the specificied
%	tranform or fail loading.  'preferredForm' will be 's', if it is
%	default or empty.	- Jeff Gunter
%
%  Example: flip_lr('avg152T1_LR_nifti.nii', 'flipped_lr.nii');
%           flip_lr('avg152T1_RL_nifti.nii', 'flipped_rl.nii');
%
%  You will find that 'avg152T1_LR_nifti.nii' and 'avg152T1_RL_nifti.nii'
%  are the same, and 'flipped_lr.nii' and 'flipped_rl.nii' are also the
%  the same, but they are L-R flipped from 'avg152T1_*'.
%
%  NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function flip_lr(original_fn, flipped_fn, old_RGB, tolerance, preferredForm)

   if ~exist('original_fn','var') | ~exist('flipped_fn','var')
      error('Usage: flip_lr(original_fn, flipped_fn, [old_RGB],[tolerance])');
   end

   if ~exist('old_RGB','var') | isempty(old_RGB)
      old_RGB = 0;
   end

   if ~exist('tolerance','var') | isempty(tolerance)
      tolerance = 0.1;
   end

   if ~exist('preferredForm','var') | isempty(preferredForm)
      preferredForm= 's';				% Jeff
   end

   nii = load_nii(original_fn, [], [], [], [], old_RGB, tolerance, preferredForm);
   M = diag(nii.hdr.dime.pixdim(2:5));
   M(1:3,4) = -M(1:3,1:3)*(nii.hdr.hist.originator(1:3)-1)';
   M(1,:) = -1*M(1,:);
   nii.hdr.hist.sform_code = 1;
   nii.hdr.hist.srow_x = M(1,:);
   nii.hdr.hist.srow_y = M(2,:);
   nii.hdr.hist.srow_z = M(3,:);
   save_nii(nii, flipped_fn);

   return;					% flip_lr

author: CoSMoMVPA

repo: CoSMoMVPA

save_path: github-repos/MATLAB/CoSMoMVPA-CoSMoMVPA


text:
function feature = mk_mfcc(x,use_logE, fs)
if nargin<3
    fs = 8000;
end
if nargin<2
    use_logE = 0;
end

feature = fbank2mfcc(wav2fbank(x, fs));
if use_logE == 1
    feature(:,13) = comp_logE(x);
end
feature(:,14:26) = comp_delta(feature,3);
feature(:,27:39) = comp_delta(feature(:,14:26),2);

author: singaxiong

repo: SignalGraph