text:
function [out] = GLCM_Features1(glcmin,pairs)
%
% GLCM_Features1 helps to calculate the features from the different GLCMs
% that are input to the function. The GLCMs are stored in a i x j x n
% matrix, where n is the number of GLCMs calculated usually due to the
% different orientation and displacements used in the algorithm. Usually
% the values i and j are equal to 'NumLevels' parameter of the GLCM
% computing function graycomatrix(). Note that matlab quantization values
% belong to the set {1,..., NumLevels} and not from {0,...,(NumLevels-1)}
% as provided in some references
% http://www.mathworks.com/access/helpdesk/help/toolbox/images/graycomatrix
% .html
%
% Although there is a function graycoprops() in Matlab Image Processing
% Toolbox that computes four parameters Contrast, Correlation, Energy,
% and Homogeneity. The paper by Haralick suggests a few more parameters
% that are also computed here. The code is not fully vectorized and hence
% is not an efficient implementation but it is easy to add new features
% based on the GLCM using this code. Takes care of 3 dimensional glcms
% (multiple glcms in a single 3D array)
%
% If you find that the values obtained are different from what you expect
% or if you think there is a different formula that needs to be used
% from the ones used in this code please let me know.
% A few questions which I have are listed in the link
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/239608
%
% I plan to submit a vectorized version of the code later and provide
% updates based on replies to the above link and this initial code.
%
% Features computed
% Autocorrelation: [2]                      (out.autoc)
% Contrast: matlab/[1,2]                    (out.contr)
% Correlation: matlab                       (out.corrm)
% Correlation: [1,2]                        (out.corrp)
% Cluster Prominence: [2]                   (out.cprom)
% Cluster Shade: [2]                        (out.cshad)
% Dissimilarity: [2]                        (out.dissi)
% Energy: matlab / [1,2]                    (out.energ)
% Entropy: [2]                              (out.entro)
% Homogeneity: matlab                       (out.homom)
% Homogeneity: [2]                          (out.homop)
% Maximum probability: [2]                  (out.maxpr)
% Sum of sqaures: Variance [1]              (out.sosvh)
% Sum average [1]                           (out.savgh)
% Sum variance [1]                          (out.svarh)
% Sum entropy [1]                           (out.senth)
% Difference variance [1]                   (out.dvarh)
% Difference entropy [1]                    (out.denth)
% Information measure of correlation1 [1]   (out.inf1h)
% Informaiton measure of correlation2 [1]   (out.inf2h)
% Inverse difference (INV) is homom [3]     (out.homom)
% Inverse difference normalized (INN) [3]   (out.indnc)
% Inverse difference moment normalized [3]  (out.idmnc)
%
% The maximal correlation coefficient was not calculated due to
% computational instability
% http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
%
% Formulae from MATLAB site (some look different from
% the paper by Haralick but are equivalent and give same results)
% Example formulae:
% Contrast = sum_i(sum_j(  (i-j)^2 * p(i,j) ) ) (same in matlab/paper)
% Correlation = sum_i( sum_j( (i - u_i)(j - u_j)p(i,j)/(s_i.s_j) ) ) (m)
% Correlation = sum_i( sum_j( ((ij)p(i,j) - u_x.u_y) / (s_x.s_y) ) ) (p[2])
% Energy = sum_i( sum_j( p(i,j)^2 ) )           (same in matlab/paper)
% Homogeneity = sum_i( sum_j( p(i,j) / (1 + |i-j|) ) ) (as in matlab)
% Homogeneity = sum_i( sum_j( p(i,j) / (1 + (i-j)^2) ) ) (as in paper)
%
% Where:
% u_i = u_x = sum_i( sum_j( i.p(i,j) ) ) (in paper [2])
% u_j = u_y = sum_i( sum_j( j.p(i,j) ) ) (in paper [2])
% s_i = s_x = sum_i( sum_j( (i - u_x)^2.p(i,j) ) ) (in paper [2])
% s_j = s_y = sum_i( sum_j( (j - u_y)^2.p(i,j) ) ) (in paper [2])
%
%
% Normalize the glcm:
% Compute the sum of all the values in each glcm in the array and divide
% each element by it sum
%
% Haralick uses 'Symmetric' = true in computing the glcm
% There is no Symmetric flag in the Matlab version I use hence
% I add the diagonally opposite pairs to obtain the Haralick glcm
% Here it is assumed that the diagonally opposite orientations are paired
% one after the other in the matrix
% If the above assumption is true with respect to the input glcm then
% setting the flag 'pairs' to 1 will compute the final glcms that would result
% by setting 'Symmetric' to true. If your glcm is computed using the
% Matlab version with 'Symmetric' flag you can set the flag 'pairs' to 0
%
% References:
% 1. R. M. Haralick, K. Shanmugam, and I. Dinstein, Textural Features of
% Image Classification, IEEE Transactions on Systems, Man and Cybernetics,
% vol. SMC-3, no. 6, Nov. 1973
% 2. L. Soh and C. Tsatsoulis, Texture Analysis of SAR Sea Ice Imagery
% Using Gray Level Co-Occurrence Matrices, IEEE Transactions on Geoscience
% and Remote Sensing, vol. 37, no. 2, March 1999.
% 3. D A. Clausi, An analysis of co-occurrence texture statistics as a
% function of grey level quantization, Can. J. Remote Sensing, vol. 28, no.
% 1, pp. 45-62, 2002
% 4. http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
%
%
% Example:
%
% Usage is similar to graycoprops() but needs extra parameter 'pairs' apart
% from the GLCM as input
% I = imread('circuit.tif');
% GLCM2 = graycomatrix(I,'Offset',[2 0;0 2]);
% stats = GLCM_features1(GLCM2,0)
% The output is a structure containing all the parameters for the different
% GLCMs
%
% [Avinash Uppuluri: avinash_uv@yahoo.com: Last modified: 11/20/08]

% If 'pairs' not entered: set pairs to 0
if ((nargin > 2) || (nargin == 0))
   error('Too many or too few input arguments. Enter GLCM and pairs.');
elseif ( (nargin == 2) )
    if ((size(glcmin,1) <= 1) || (size(glcmin,2) <= 1))
       error('The GLCM should be a 2-D or 3-D matrix.');
    elseif ( size(glcmin,1) ~= size(glcmin,2) )
        error('Each GLCM should be square with NumLevels rows and NumLevels cols');
    end
elseif (nargin == 1) % only GLCM is entered
    pairs = 0; % default is numbers and input 1 for percentage
    if ((size(glcmin,1) <= 1) || (size(glcmin,2) <= 1))
       error('The GLCM should be a 2-D or 3-D matrix.');
    elseif ( size(glcmin,1) ~= size(glcmin,2) )
       error('Each GLCM should be square with NumLevels rows and NumLevels cols');
    end
end

format long e
if (pairs == 1)
    newn = 1;
    for nglcm = 1:2:size(glcmin,3)
        glcm(:,:,newn)  = glcmin(:,:,nglcm) + glcmin(:,:,nglcm+1);
        newn = newn + 1;
    end
elseif (pairs == 0)
    glcm = glcmin;
end

size_glcm_1 = size(glcm,1);
size_glcm_2 = size(glcm,2);
size_glcm_3 = size(glcm,3);

% checked
out.autoc = zeros(1,size_glcm_3); % Autocorrelation: [2]
out.contr = zeros(1,size_glcm_3); % Contrast: matlab/[1,2]
out.corrm = zeros(1,size_glcm_3); % Correlation: matlab
out.corrp = zeros(1,size_glcm_3); % Correlation: [1,2]
out.cprom = zeros(1,size_glcm_3); % Cluster Prominence: [2]
out.cshad = zeros(1,size_glcm_3); % Cluster Shade: [2]
out.dissi = zeros(1,size_glcm_3); % Dissimilarity: [2]
out.energ = zeros(1,size_glcm_3); % Energy: matlab / [1,2]
out.entro = zeros(1,size_glcm_3); % Entropy: [2]
out.homom = zeros(1,size_glcm_3); % Homogeneity: matlab
out.homop = zeros(1,size_glcm_3); % Homogeneity: [2]
out.maxpr = zeros(1,size_glcm_3); % Maximum probability: [2]

out.sosvh = zeros(1,size_glcm_3); % Sum of sqaures: Variance [1]
out.savgh = zeros(1,size_glcm_3); % Sum average [1]
out.svarh = zeros(1,size_glcm_3); % Sum variance [1]
out.senth = zeros(1,size_glcm_3); % Sum entropy [1]
out.dvarh = zeros(1,size_glcm_3); % Difference variance [4]
%out.dvarh2 = zeros(1,size_glcm_3); % Difference variance [1]
out.denth = zeros(1,size_glcm_3); % Difference entropy [1]
out.inf1h = zeros(1,size_glcm_3); % Information measure of correlation1 [1]
out.inf2h = zeros(1,size_glcm_3); % Informaiton measure of correlation2 [1]
%out.mxcch = zeros(1,size_glcm_3);% maximal correlation coefficient [1]
%out.invdc = zeros(1,size_glcm_3);% Inverse difference (INV) is homom [3]
out.indnc = zeros(1,size_glcm_3); % Inverse difference normalized (INN) [3]
out.idmnc = zeros(1,size_glcm_3); % Inverse difference moment normalized [3]

% correlation with alternate definition of u and s
%out.corrm2 = zeros(1,size_glcm_3); % Correlation: matlab
%out.corrp2 = zeros(1,size_glcm_3); % Correlation: [1,2]

glcm_sum  = zeros(size_glcm_3,1);
glcm_mean = zeros(size_glcm_3,1);
glcm_var  = zeros(size_glcm_3,1);

% http://www.fp.ucalgary.ca/mhallbey/glcm_mean.htm confuses the range of
% i and j used in calculating the means and standard deviations.
% As of now I am not sure if the range of i and j should be [1:Ng] or
% [0:Ng-1]. I am working on obtaining the values of mean and std that get
% the values of correlation that are provided by matlab.
u_x = zeros(size_glcm_3,1);
u_y = zeros(size_glcm_3,1);
s_x = zeros(size_glcm_3,1);
s_y = zeros(size_glcm_3,1);

% % alternate values of u and s
% u_x2 = zeros(size_glcm_3,1);
% u_y2 = zeros(size_glcm_3,1);
% s_x2 = zeros(size_glcm_3,1);
% s_y2 = zeros(size_glcm_3,1);

% checked p_x p_y p_xplusy p_xminusy
p_x = zeros(size_glcm_1,size_glcm_3); % Ng x #glcms[1]
p_y = zeros(size_glcm_2,size_glcm_3); % Ng x #glcms[1]
p_xplusy = zeros((size_glcm_1*2 - 1),size_glcm_3); %[1]
p_xminusy = zeros((size_glcm_1),size_glcm_3); %[1]
% checked hxy hxy1 hxy2 hx hy
hxy  = zeros(size_glcm_3,1);
hxy1 = zeros(size_glcm_3,1);
hx   = zeros(size_glcm_3,1);
hy   = zeros(size_glcm_3,1);
hxy2 = zeros(size_glcm_3,1);

%Q    = zeros(size(glcm));

for k = 1:size_glcm_3 % number glcms

    glcm_sum(k) = sum(sum(glcm(:,:,k)));
    glcm(:,:,k) = glcm(:,:,k)./glcm_sum(k); % Normalize each glcm
    glcm_mean(k) = mean2(glcm(:,:,k)); % compute mean after norm
    glcm_var(k)  = (std2(glcm(:,:,k)))^2;

    for i = 1:size_glcm_1

        for j = 1:size_glcm_2

            out.contr(k) = out.contr(k) + (abs(i - j))^2.*glcm(i,j,k);
            out.dissi(k) = out.dissi(k) + (abs(i - j)*glcm(i,j,k));
            out.energ(k) = out.energ(k) + (glcm(i,j,k).^2);
            out.entro(k) = out.entro(k) - (glcm(i,j,k)*log(glcm(i,j,k) + eps));
            out.homom(k) = out.homom(k) + (glcm(i,j,k)/( 1 + abs(i-j) ));
            out.homop(k) = out.homop(k) + (glcm(i,j,k)/( 1 + (i - j)^2));
            % [1] explains sum of squares variance with a mean value;
            % the exact definition for mean has not been provided in
            % the reference: I use the mean of the entire normalized glcm
            out.sosvh(k) = out.sosvh(k) + glcm(i,j,k)*((i - glcm_mean(k))^2);

            %out.invdc(k) = out.homom(k);
            out.indnc(k) = out.indnc(k) + (glcm(i,j,k)/( 1 + (abs(i-j)/size_glcm_1) ));
            out.idmnc(k) = out.idmnc(k) + (glcm(i,j,k)/( 1 + ((i - j)/size_glcm_1)^2));
            u_x(k)          = u_x(k) + (i)*glcm(i,j,k); % changed 10/26/08
            u_y(k)          = u_y(k) + (j)*glcm(i,j,k); % changed 10/26/08
            % code requires that Nx = Ny
            % the values of the grey levels range from 1 to (Ng)
        end

    end
    out.maxpr(k) = max(max(glcm(:,:,k)));
end
% glcms have been normalized:
% The contrast has been computed for each glcm in the 3D matrix
% (tested) gives similar results to the matlab function

for k = 1:size_glcm_3

    for i = 1:size_glcm_1

        for j = 1:size_glcm_2
            p_x(i,k) = p_x(i,k) + glcm(i,j,k);
            p_y(i,k) = p_y(i,k) + glcm(j,i,k); % taking i for j and j for i
            if (ismember((i + j),[2:2*size_glcm_1]))
                p_xplusy((i+j)-1,k) = p_xplusy((i+j)-1,k) + glcm(i,j,k);
            end
            if (ismember(abs(i-j),[0:(size_glcm_1-1)]))
                p_xminusy((abs(i-j))+1,k) = p_xminusy((abs(i-j))+1,k) +...
                    glcm(i,j,k);
            end
        end
    end

%     % consider u_x and u_y and s_x and s_y as means and standard deviations
%     % of p_x and p_y
%     u_x2(k) = mean(p_x(:,k));
%     u_y2(k) = mean(p_y(:,k));
%     s_x2(k) = std(p_x(:,k));
%     s_y2(k) = std(p_y(:,k));

end

% marginal probabilities are now available [1]
% p_xminusy has +1 in index for matlab (no 0 index)
% computing sum average, sum variance and sum entropy:
for k = 1:(size_glcm_3)

    for i = 1:(2*(size_glcm_1)-1)
        out.savgh(k) = out.savgh(k) + (i+1)*p_xplusy(i,k);
        % the summation for savgh is for i from 2 to 2*Ng hence (i+1)
        out.senth(k) = out.senth(k) - (p_xplusy(i,k)*log(p_xplusy(i,k) + eps));
    end

end
% compute sum variance with the help of sum entropy
for k = 1:(size_glcm_3)

    for i = 1:(2*(size_glcm_1)-1)
        out.svarh(k) = out.svarh(k) + (((i+1) - out.senth(k))^2)*p_xplusy(i,k);
        % the summation for savgh is for i from 2 to 2*Ng hence (i+1)
    end

end
% compute difference variance, difference entropy,
for k = 1:size_glcm_3
% out.dvarh2(k) = var(p_xminusy(:,k));
% but using the formula in
% http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
% we have for dvarh
    for i = 0:(size_glcm_1-1)
        out.denth(k) = out.denth(k) - (p_xminusy(i+1,k)*log(p_xminusy(i+1,k) + eps));
        out.dvarh(k) = out.dvarh(k) + (i^2)*p_xminusy(i+1,k);
    end
end

% compute information measure of correlation(1,2) [1]
for k = 1:size_glcm_3
    hxy(k) = out.entro(k);
    for i = 1:size_glcm_1

        for j = 1:size_glcm_2
            hxy1(k) = hxy1(k) - (glcm(i,j,k)*log(p_x(i,k)*p_y(j,k) + eps));
            hxy2(k) = hxy2(k) - (p_x(i,k)*p_y(j,k)*log(p_x(i,k)*p_y(j,k) + eps));
%             for Qind = 1:(size_glcm_1)
%                 Q(i,j,k) = Q(i,j,k) +...
%                     ( glcm(i,Qind,k)*glcm(j,Qind,k) / (p_x(i,k)*p_y(Qind,k)) );
%             end
        end
        hx(k) = hx(k) - (p_x(i,k)*log(p_x(i,k) + eps));
        hy(k) = hy(k) - (p_y(i,k)*log(p_y(i,k) + eps));
    end
    out.inf1h(k) = ( hxy(k) - hxy1(k) ) / ( max([hx(k),hy(k)]) );
    out.inf2h(k) = ( 1 - exp( -2*( hxy2(k) - hxy(k) ) ) )^0.5;
%     eig_Q(k,:)   = eig(Q(:,:,k));
%     sort_eig(k,:)= sort(eig_Q(k,:),'descend');
%     out.mxcch(k) = sort_eig(k,2)^0.5;
% The maximal correlation coefficient was not calculated due to
% computational instability
% http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
end

corm = zeros(size_glcm_3,1);
corp = zeros(size_glcm_3,1);
% using http://www.fp.ucalgary.ca/mhallbey/glcm_variance.htm for s_x s_y
for k = 1:size_glcm_3
    for i = 1:size_glcm_1
        for j = 1:size_glcm_2
            s_x(k)  = s_x(k)  + (((i) - u_x(k))^2)*glcm(i,j,k);
            s_y(k)  = s_y(k)  + (((j) - u_y(k))^2)*glcm(i,j,k);
            corp(k) = corp(k) + ((i)*(j)*glcm(i,j,k));
            corm(k) = corm(k) + (((i) - u_x(k))*((j) - u_y(k))*glcm(i,j,k));
            out.cprom(k) = out.cprom(k) + (((i + j - u_x(k) - u_y(k))^4)*...
                glcm(i,j,k));
            out.cshad(k) = out.cshad(k) + (((i + j - u_x(k) - u_y(k))^3)*...
                glcm(i,j,k));
        end
    end
    % using http://www.fp.ucalgary.ca/mhallbey/glcm_variance.htm for s_x
    % s_y : This solves the difference in value of correlation and might be
    % the right value of standard deviations required
    % According to this website there is a typo in [2] which provides
    % values of variance instead of the standard deviation hence a square
    % root is required as done below:
    s_x(k) = s_x(k) ^ 0.5;
    s_y(k) = s_y(k) ^ 0.5;
    out.autoc(k) = corp(k);
    out.corrp(k) = (corp(k) - u_x(k)*u_y(k))/(s_x(k)*s_y(k));
    out.corrm(k) = corm(k) / (s_x(k)*s_y(k));
%     % alternate values of u and s
%     out.corrp2(k) = (corp(k) - u_x2(k)*u_y2(k))/(s_x2(k)*s_y2(k));
%     out.corrm2(k) = corm(k) / (s_x2(k)*s_y2(k));
end
% Here the formula in the paper out.corrp and the formula in matlab
% out.corrm are equivalent as confirmed by the similar results obtained

% % The papers have a slightly different formular for Contrast
% % I have tested here to find this formula in the papers provides the
% % same results as the formula provided by the matlab function for
% % Contrast (Hence this part has been commented)
% out.contrp = zeros(size_glcm_3,1);
% contp = 0;
% Ng = size_glcm_1;
% for k = 1:size_glcm_3
%     for n = 0:(Ng-1)
%         for i = 1:Ng
%             for j = 1:Ng
%                 if (abs(i-j) == n)
%                     contp = contp + glcm(i,j,k);
%                 end
%             end
%         end
%         out.contrp(k) = out.contrp(k) + n^2*contp;
%         contp = 0;
%     end
%
% end

%       GLCM Features (Soh, 1999; Haralick, 1973; Clausi 2002)
%           f1. Uniformity / Energy / Angular Second Moment (done)
%           f2. Entropy (done)
%           f3. Dissimilarity (done)
%           f4. Contrast / Inertia (done)
%           f5. Inverse difference
%           f6. correlation
%           f7. Homogeneity / Inverse difference moment
%           f8. Autocorrelation
%           f9. Cluster Shade
%          f10. Cluster Prominence
%          f11. Maximum probability
%          f12. Sum of Squares
%          f13. Sum Average
%          f14. Sum Variance
%          f15. Sum Entropy
%          f16. Difference variance
%          f17. Difference entropy
%          f18. Information measures of correlation (1)
%          f19. Information measures of correlation (2)
%          f20. Maximal correlation coefficient
%          f21. Inverse difference normalized (INN)
%          f22. Inverse difference moment normalized (IDN)


text:
function T = spm_type(x, arg)
% Translate data type specifiers between SPM & MATLAB representations
% FORMAT T = spm_type(x, arg)
% x    - specifier
% T    - type
% arg  - optional string argument, can be:
%         - 'maxval'  - return maximum allowed value.
%         - 'minval'  - return minimum allowed value.
%         - 'nanrep'  - return 1 if there is a NaN representation.
%         - 'bits'    - return the number of bits per voxel.
%         - 'intt'    - return 1 if values rounded to nearest integer.
%__________________________________________________________________________
%
% Format specifiers are based on NIFTI-1.
% If the input is a number then the corresponding MATLAB string is
% returned by default.
% If the input is a string then the appropriate TYPE is returned.
% However, if the optional arg argument is supplied then other
% information will be returned instead.
%
% With no arguments, a list of data types is returned.
%__________________________________________________________________________
% Copyright (C) 1996-2014 Wellcome Trust Centre for Neuroimaging

% John Ashburner & Andrew Holmes
% $Id: spm_type.m 5925 2014-03-20 16:47:44Z guillaume $

prec   = {'uint8','int16','int32','float32','float64','int8','uint16','uint32'};
conv   = {@uint8,@int16,@int32,@single,@double,@int8,@uint16,@uint32};
types  = [    2      4      8   16   64   256    512    768];
maxval = [2^8-1 2^15-1 2^31-1  Inf  Inf 2^7-1 2^16-1 2^32-1];
minval = [    0  -2^15  -2^31 -Inf -Inf  -2^7      0      0];
nanrep = [    0      0      0    1    1     0      0      0];
bits   = [    8     16     32   32   64     8     16     32];
intt   = [    1      1      1    0    0     1      1      1];

if ~nargin
    T = types;
    return;
end

if ischar(x)
    sel = find(strcmpi(prec,deblank(x)));
else
    sel = find(ismember(types,x));
end
if nargin == 1
    if ischar(x)
        if isempty(sel), T = NaN;
        else T = types(sel); end
    else
        if isempty(sel), T = 'unknown';
        else T = char(prec(sel)); end
    end
elseif isempty(sel)
    T = NaN;
else
    switch lower(arg)
        case 'maxval', T = maxval(sel);
        case 'minval', T = minval(sel);
        case 'nanrep', T = nanrep(sel);
        case 'bits',   T = bits(sel);
        case 'intt',   T = intt(sel);
        case 'conv',   T = conv(sel);
        otherwise,     T = NaN;
    end
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function Callbacks_ideal_vocal_tract_GUI25(f,C)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x=C{1,1};
y=C{1,2};
a=C{1,3};
b=C{1,4};
u=C{1,5};
v=C{1,6};
m=C{1,7};
n=C{1,8};
lengthbutton=C{1,9};
widthbutton=C{1,10};
enterType=C{1,11};
enterString=C{1,12};
enterLabel=C{1,13};
noPanels=C{1,14};
noGraphicPanels=C{1,15};
noButtons=C{1,16};
labelDist=C{1,17};%distance that the label is below the button
noTitles=C{1,18};
buttonTextSize=C{1,19};
labelTextSize=C{1,20};
textboxFont=C{1,21};
textboxString=C{1,22};
textboxWeight=C{1,23};
textboxAngle=C{1,24};
labelHeight=C{1,25};
fileName=C{1,26};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%PANELS
for j=0:noPanels-1
uipanel('Parent',f,...
'Units','Normalized',...
'Position',[x(1+4*j) y(1+4*j) x(2+4*j)-x(1+4*j) y(3+4*j)-y(2+4*j)]);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%GRAPHIC PANELS
for i=0:noGraphicPanels-1
switch (i+1)
case 1
graphicPanel1 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 2
graphicPanel2 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 3
graphicPanel3 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 4
graphicPanel4 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 5
graphicPanel5 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
case 6
graphicPanel6 = axes('parent',f,...
'Units','Normalized',...
'Position',[a(1+4*i) b(1+4*i) a(2+4*i)-a(1+4*i) b(3+4*i)-b(2+4*i)],...
'GridLineStyle','--');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TITLE BOXES
for k=0:noTitles-1
switch (k+1)
case 1
titleBox1 = uicontrol('parent',f,...
'Units','Normalized',...
'Position',[u(1+4*k) v(1+4*k) u(2+4*k)-u(1+4*k) v(3+4*k)-v(2+4*k)],...
'Style','text',...
'FontSize',textboxFont{k+1},...
'String',textboxString(k+1),...
'FontWeight',textboxWeight{k+1},...
'FontAngle',textboxAngle{k+1});
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%BUTTONS
for i=0:(noButtons-1)
enterColor='w';
if strcmp(enterType{i+1},'pushbutton')==1 ||strcmp(enterType{i+1},'text')==1
enterColor='default';
end
if (strcmp(enterLabel{1,(i+1)},'')==0 &&...
        strcmp(enterLabel{1,(i+1)},'...')==0) %i.e. there is a label
%creating a label for some buttons
uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i)-labelDist-labelHeight(i+1) ...
(m(2+2*i)-m(1+2*i)) labelHeight(i+1)],...
'Style','text',...
'String',enterLabel{i+1},...
'FontSize', labelTextSize(i+1),...
'HorizontalAlignment','center');
end
switch (i+1)
case 1
button1=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button1Callback);
case 2
button2=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button2Callback);
case 3
button3=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button3Callback);
case 4
button4=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button4Callback);
case 5
button5=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button5Callback);
case 6
button6=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button6Callback);
case 7
button7=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button7Callback);
case 8
button8=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button8Callback);
case 9
button9=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button9Callback);
case 10
button10=uicontrol('Parent',f,...
'Units','Normalized',...
'Position',[m(1+2*i) n(1+2*i) (m(2+2*i)-m(1+2*i)) (n(2+2*i)-n(1+2*i))],...
'Style',enterType{i+1},...
'String',enterString{i+1},...
'FontSize', buttonTextSize(1+i),...
'BackgroundColor',enterColor,...
'HorizontalAlignment','center',...
'Callback',@button10Callback);
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%USER CODE FOR THE VARIABLES AND CALLBACKS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize Variables
vowel_index=1;
ipd=100;
Lm=40;
L=400;
ilog=1;
ftable=[270 2290 3010; 390 1990 2550; 530 1840 2480; 660 1720 2410;...
    520 1190 2390; 730 1090 2390; 570 840 2410; 490 1350 1690;...
    440 1020 2240; 300 870 2240; 300 870 2240; 500 1500 2500];
btable=[40 60 75];
fk=[270 2290 3010];
bk=[40 60 75];
M=3;
fs=10000;
hn=[];
g=[];
iglottal=0;
vowels=['IY'; 'IH'; 'EH'; 'AE'; 'AH'; 'AA'; 'AO'; 'ER'; 'UH'; 'UW'; 'OO'; 'Ne'];
filename=vowels(1,:);
alpha1=40;
alpha2=20;

% nfft: size of fft and ifft for conversion from time to frequency or
% frequency to time domain
    nfft=1000;

% Name GUI
    set(f,'Name','ideal_vocal_tract_GUI');

% CALLBACKS
% Callback for button1 -- vowel choice from popuplist
 function button1Callback(h,eventdata)
     vowel_index=get(button1,'val');
     fk=ftable(vowel_index,:);
     bk=btable;
     filename=vowels(vowel_index,:);
 end

% Callback for button2 -- ipd: pitch period in msec
 function button2Callback(h,eventdata)
     ipd=str2num(get(button2,'string'));
     if ipd < 2
         waitfor(msgbox('pitch period must lie between 2 and 15 msec'));
     elseif ipd > 15
         waitfor(msgbox('pitch period must lie between 2 and 15 msec'));
     end
     ipd=ipd*fs/1000;
 end

% Callback for button3 -- Lm: frame duration in msec, (editable)
 function button3Callback(h,eventdata)
     Lm=str2num(get(button3,'string'));
      if Lm < 0
         waitfor(msgbox('The frame duration cannot be negative'));
     end
     L=round(Lm*fs/1000);
 end

% Callback for button4 -- log/linear popuplist (1=log, 2=linear); option
% for plotting vowel and ideal frequency responses
 function button4Callback(h,eventdata)
     ilog=get(button4,'val');
 end

% Callback for button10 -- iglottal: use glottal puse (1) or skip glottal
% pulse (0)
    function button10Callback(h,eventdata)
        iglottal=get(button10,'val');
    end

% Callback for button8 -- alpha1 for glottal pulse
    function button8Callback(h,eventdata)
        alpha1=str2num(get(button8,'string'))/100;
        if alpha1 < 0
            waitfor(msgbox('alpha1 must lie between 0 and 100'));
        elseif alpha1 > 1
            waitfor(msgbox('alpha1 must lie between 0 and 100'));
        end
    end

% Callback for button9 -- alpha2 for glottal pulse
    function button9Callback(h,eventdata)
        alpha2=str2num(get(button9,'string'))/100;
    if alpha2 < 0
        waitfor(msgbox('alpha2 must lie between 0 and 100'));
    elseif alpha2 > 1
        waitfor(msgbox('alpha2 must lie between 0 and 100'));
     end
    end

% Callback for button5 -- Run ideal vocal tract responses
 function button5Callback(h,eventdata)

% eps: smallest spectral magnitude for plotting log magnitude spectra
    eps=1.e-3;

% form rosenberg glottal pulse, g; rise % from alpha1; fall % from alpha2
    button8Callback(h,eventdata);
    button9Callback(h,eventdata);
    button10Callback(h,eventdata);
    button2Callback(h,eventdata);
    button3Callback(h,eventdata);
    button4Callback(h,eventdata);

    n1=round(ipd*alpha1);
    n2=round(ipd*alpha2);
    g=[];
    x1=0:n1;
    g=[g 0.5*(1-cos(pi*x1/n1))];
    x2=1:n2;
    g=[g cos(pi*x2/(2*n2))];
    g=[g zeros(1,ipd-length(g))];

% create periodic voiced excitation
    e=zeros(1,nfft);
    e(1:ipd:nfft)=1;
    EM=zeros(1,nfft);
    EM(1:round(nfft/ipd):nfft)=1;

% glottal_pulse excitation based on iglottal (1: use glottal pulse, 0 skip
% glottal pulse
if (iglottal == 1)
    eg=conv(e,g);
else
    eg=e;
end
    ege(1:nfft)=0;
    ege(1:nfft)=eg(1:nfft);
    EGM=abs(fft(ege,nfft));

% convert to second order system format
    zk=exp(-2*pi*bk/fs);
    fkn=cos(2*pi*fk/fs);

% transform into frequency domain
    HR=[];
    HT=zeros(1,nfft);
    HP=(complex(ones(1,nfft),zeros(1,nfft)))';
    for k=1:M
        B=1-2*zk(k)*fkn(k)+zk(k)*zk(k);
        A=[1 -2*zk(k)*fkn(k) zk(k)*zk(k)];
        [H,W]=freqz(B,A,nfft,'whole',fs);
        HP=HP.*H;
        HL=(abs(H));
        HT=HT.*HL';
    end

% inverse transform to give vocal tract impulse response
    hn=ifft(HP,nfft);
    hmax=max(max(hn),-min(hn));
    for index=1:nfft
        if (max(max(hn(index:nfft)),-min(hn(index:nfft)))) < hmax/50
            nl=index;
            break;
        end
    end

% plot elements of composite response in the 6 graphics Panels
% plot ideal impulse excitation signal in graphics Panel 5 (top left panel)
    reset(graphicPanel5);
    axes(graphicPanel5);

    if (iglottal == 0)
        stem(0:nfft-1,e(1:nfft),'b','LineWidth',2),axis tight;
    else
        plot(0:nfft-1,eg(1:nfft),'b','LineWidth',2),axis tight;
    end
    xpp=['Time in Samples; fs=',num2str(fs),' samples/second'];
    xlabel(xpp),ylabel('Amplitude');grid on;legend('Excitation Pulses');

% plot vowel impulse response in graphics Panel 3 (middle left panel)
    reset(graphicPanel3);
    axes(graphicPanel3);

    plot(0:nl,hn(1:nl+1),'r','LineWidth',2);
    xlabel(xpp),ylabel('Amplitude');
    axis tight, grid on;legend('Vowel Vocal Tract Impulse Response');

% plot linear or log magnitude frequency response of ideal impulse
% excitation signal in graphics Panel 6 (top right panel); use lin/log
% option to decide whether to plot linear or log magnitude responses
    if (ilog == 1)
        reset(graphicPanel6);
        axes(graphicPanel6);
        plot(W(1:nfft/2+1),20*log10(EGM(1:nfft/2+1)+eps),'b','LineWidth',2);
        axis tight,xlabel('Frequency in Hz'),ylabel('Log Magnitude (dB)');
        legend('Excitation Pulses Log Spectrum');

% plot linear or log magnitude frequency response of ideal vocal tract
% impulse response in graphics Panel 4 (middle right panel); use lin/log
% option to decide whether to plot linear or log magnitude response
        reset(graphicPanel4);
        axes(graphicPanel4);
        plot(W(1:nfft/2+1),20*log10(abs(HP(1:nfft/2+1))),'r','LineWidth',2);
        axis tight,grid on, xlabel('Frequency in Hz');
        ylabel('Log Magnitude (dB)');legend('Vocal Tract Log Spectrum');
    else
        reset(graphicPanel6);
        axes(graphicPanel6);
        plot(W(1:nfft/2+1),EGM(1:nfft/2+1),'b','LineWidth',2);
        axis tight,xlabel('Frequency in Hz'),ylabel('Magnitude');
        legend('Excitation Pulses Linear Spectrum');

        reset(graphicPanel4);
        axes(graphicPanel4);
        plot(W(1:nfft/2+1),abs(HP(1:nfft/2+1)),'r','LineWidth',2);
        xlabel('Frequency in Hz'),ylabel('Magnitude'),axis tight,grid on;
        legend('Vocal Tract Linear Spectrum');
    end
    axis tight;

% take product of magnitudes, convolve time responses
    he=zeros(1,L+nfft-1);
    for index=1:ipd:L
        he(index:index+nfft-1)=he(index:index+nfft-1)+hn(1:nfft)';
    end
    EMV=EGM(1:nfft/2+1).*abs(HP(1:nfft/2+1))';

% plot ideal periodic vowel impulse response in graphics Panel 1 (lower
% left panel); note: no glottal pulse used in this simple ideal vocal tract
% model
    reset(graphicPanel1);
    axes(graphicPanel1);
    plot(0:L-1,he(1:L),'b','LineWidth',2);
    xlabel(xpp),ylabel('Amplitude');
    axis tight,grid on;legend('Vowel Waveform');

% plot linear or log magnitude frequency response of ideal synthetic vowel
% impulse response in graphics Panel 2 (bottom right panel); use lin/log
% option to decide whether to plot linear or log magnitude response
    if (ilog == 1)
        reset(graphicPanel2);
        axes(graphicPanel2);
        plot(W(1:nfft/2+1),20*log10(EMV(1:nfft/2+1)+eps),'r','LineWidth',2);
        xlabel('Frequency in Hz'),ylabel('Log Magnitude (dB)');
        axis tight,grid on;legend('Vowel Log Spectrum');
    else
        reset(graphicPanel2);
        axes(graphicPanel2);
        plot(W(1:nfft/2+1),EMV(1:nfft/2+1),'r','LineWidth',2);
        xlabel('Frequency in Hz'),ylabel('Magnitude');
        axis tight, grid on;legend('Vowel Linear Spectrum');
    end

% set titleBox1 information for display
    stitle=sprintf('Ideal Vocal Tract -- Vowel:%s, ipd:%d, L:%d, fs:%d, eps:%6.3f',filename,ipd,L,fs,eps);
    set(titleBox1,'String',stitle);
    set(titleBox1,'FontSize',20);
 end

% Callback for button7 -- Play Synthetic Vowel
    function button7Callback(h,eventdata)
% play back synthetic speech output
    e=zeros(1,fs);
    e(1:ipd:fs)=1;
    if (iglottal == 1)
        y1=conv(e,g);
        yout=conv(y1,hn);
    else
        yout=conv(e,hn);
    end
    % uiwait(msgbox('Play out synthetic speech','speech out'));
    soundsc(yout,fs);
    end

% Callback for button6 -- Close GUI
 function button6Callback(h,eventdata)
     close(gcf);
 end
end

author: Sable

repo: mcbench-benchmarks


text:
function dtiMrMeshScreenSave(handles,fname)
%
%  dtiMrMeshScreenSave(handles,fname)
%
%Author: Wandell, Dougherty
%Purpose:
%   Save whatever is in the dti mrMesh window to a file
%

rgb = mrmGet(handles.mrMesh,'screenshot')/255;
imwrite(rgb, fname);
disp(['Screenshot saved to ' fname '.']);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [ r, t ] = xy_to_polar ( xy )

%*****************************************************************************80
%
%% XY_TO_POLAR converts XY coordinates to polar coordinates.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2010
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real XY(2,1), the Cartesian coordinates.
%
%    Output, real R, T, the radius and angle (in radians).
%
  r = norm ( xy );

  if ( r == 0.0 )
    t = 0.0;
  else
    t = atan2 ( xy(2,1), xy(1,1) );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [XYZ, Z, M] = pr_get_spm_results
% Fetch SPM results and return as point list
% FORMAT [XYZ, Z, M] = pr_get_spm_results
%
% Outputs
% XYZ    - XYZ point list in voxels (empty if not found)
% Z      - values at points in XYZ
% M      - 4x4 voxel -> world transformation matrix
%__________________________________________________________________________

% Matthew Brett
% $Id: pr_get_spm_results.m 6623 2015-12-03 18:38:08Z guillaume $

errstr = '''Cannot find SPM results in workspace''';
[XYZ,Z,M] = deal([]);

have_res = evalin('base', 'exist(''xSPM'', ''var'')');
if ~have_res, return, end
xSPM = evalin('base', 'xSPM', ['error(' errstr ')']);
XYZ = xSPM.XYZ;
Z   = xSPM.Z;
M   = xSPM.M;

author: spm


text:
function [mesh1,mesh2] = mshSplit(mesh,X0,U)
%+========================================================================+
%|                                                                        |
%|                 OPENMSH - LIBRARY FOR MESH MANAGEMENT                  |
%|           openMsh is part of the GYPSILAB toolbox for Matlab           |
%|                                                                        |
%| COPYRIGHT : Matthieu Aussal (c) 2017-2018.                             |
%| PROPERTY  : Centre de Mathematiques Appliquees, Ecole polytechnique,   |
%| route de Saclay, 91128 Palaiseau, France. All rights reserved.         |
%| LICENCE   : This program is free software, distributed in the hope that|
%| it will be useful, but WITHOUT ANY WARRANTY. Natively, you can use,    |
%| redistribute and/or modify it under the terms of the GNU General Public|
%| License, as published by the Free Software Foundation (version 3 or    |
%| later,  http://www.gnu.org/licenses). For private use, dual licencing  |
%| is available, please contact us to activate a "pay for remove" option. |
%| CONTACT   : matthieu.aussal@polytechnique.edu                          |
%| WEBSITE   : www.cmap.polytechnique.fr/~aussal/gypsilab                 |
%|                                                                        |
%| Please acknowledge the gypsilab toolbox in programs or publications in |
%| which you use it.                                                      |
%|________________________________________________________________________|
%|   '&`   |                                                              |
%|    #    |   FILE       : mshSplit.m                                    |
%|    #    |   VERSION    : 0.51                                          |
%|   _#_   |   AUTHOR(S)  : Matthieu Aussal                               |
%|  ( # )  |   CREATION   : 14.03.2019                                    |
%|  / 0 \  |   LAST MODIF : 21.06.2019                                    |
%| ( === ) |   SYNOPSIS   : Split mesh with planar cut                    |
%|  `---'  |                (to be clearly improved)                      |
%+========================================================================+

% Normalize direction
N = U./norm(U);

% Split mesh
lambda = (mesh.ctr-ones(length(mesh),1)*X0) * N';
ind    = (lambda>0);
mesh1  = mesh.sub(ind);
mesh2  = mesh.sub(~ind);

% Extract free boundary
bound  = mesh1.bnd;
mu     = (bound.vtx-ones(size(bound.vtx,1),1)*X0) * N';

% Move upper boundary
[~,IA] = intersect(msh(mesh1.vtx),msh(bound.vtx));
mesh1.vtx(IA,:) = mesh1.vtx(IA,:) - mu*N;

% Move lower boundary
[~,IA] = intersect(msh(mesh2.vtx),msh(bound.vtx));
mesh2.vtx(IA,:) = mesh2.vtx(IA,:) - mu*N;
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function tf = all(I)

tf = all(I.Z(:));

author: GERSL

repo: CCDC

save_path: github-repos/MATLAB/GERSL-CCDC


text:
function v = ise(p,q,type)
%
% ise(p,q [,'type'])  -- estimate the integrated squared error between
%                        two densities p,q
%   type:
%    [double] -- use "epsilon-exact" product with this value for epsilon
%    'p','q'  -- use the samples at p (or at q)
%    'pq'     -- use both samples at p & q
%
if (nargin < 3) type = 0; end;
if (isa(type,'double'))  v = iseEpsilon(p,q,type);
else
  switch type
    % Three different monte-carlo estimates (different proposals)
    case {'p'}, x = getPoints(p); quo = evaluate(p,x);
                v = mean( (evaluate(p,x) - evaluate(q,x)).^2 ./ quo);
    case {'q'}, x = getPoints(q); quo = evaluate(q,x);
                v = mean( (evaluate(p,x) - evaluate(q,x)).^2 ./ quo);
    case {'pq'},x = [getPoints(p),getPoints(q)]; quo = .5*(evaluate(p,x)+evaluate(q,x));
                v = mean( (evaluate(p,x) - evaluate(q,x)).^2 ./ quo);
    % and one (possible) estimate from Mark Girolami
    case {'abs'},   eval1 = evaluate(p,p); eval2 = evaluate(q,p);
                    if (min(eval2) == 0) v = inf;
                    else v = getWeights(p) * (eval1 .* abs(log(eval1./eval2)))';
                    end;

  end;
end;

author: ShapeNet

repo: RenderForCNN


text:
function [grad, roiCoords, view] = checkSpatialGradMap(view, scan, roi, preserveCoords);
%
% [grad, roiCoords, view] = checkSpatialGradMap(view, [scan=1], [roi], [preserveCoords=0]);
%
% Get a spatial gradient map for detrending purposes, loading and/or
% computing in the view as needed. Returns a gradient matrix (size
% is dataSize(view)) for the selected scan.
%
% If an ROI is specified (see tc_roiStruct for the ways to specify), will
% also return a 3xN set of coordinates into the gradient map corresponding
% to the ROI, appropriate for the current view. The roiCoords are intended
% for use with er_preporcessTSeries.
%
% ras, 02/2007
if notDefined('view'),      view = getCurView;          end
if notDefined('scan'),      scan = view.curScan;        end
if notDefined('preserveCoords'), preserveCoords = 0;    end

grad = [];
roiCoords = [];

if ~isfield(view, 'spatialGrad') | isempty(view.spatialGrad)
    % load and/or compute
    mapPath = fullfile(dataDir(view), 'spatialGrad.mat');
    if ~exist(mapPath, 'file')  % offer to compute it
        q = ['The inomogeity correction flag is set to use the spatial ' ...
             'gradient map (inhomoCorrect=3). This map is not found. ' ...
             'Compute it now? '];
        resp = questdlg(q);
        if ~isequal(resp, 'Yes')
            error('Aborted--no spatial gradient map.')
        end
        view = computeSpatialGradient(view);
    end

    view = loadSpatialGradient(view);

    % ensure the loaded spatial grad map propagates back
    % to the base workspace (global variable):
    updateGlobal(view);
end

grad = view.spatialGrad{scan};

if exist('roi', 'var') & ~isempty(roi)
    % get map coordinates for the ROI
    roi = tc_roiStruct(view, roi);

    if preserveCoords==0
       subCoords = roiSubCoords(view, roi.coords);
    else
       subCoords = roi.coords;
    end
    [indices roiCoords] = roiIndices(view, subCoords, preserveCoords);

    % for volume/gray views, we want the indices and not coords
    % for the map (which is size 1 x nIndices). Mock up a 3xN
    % coords array, in which the columns reflec the indices into the map:
    if ismember(view.viewType, {'Volume' 'Gray'})
        roiCoords = ones(3, size(indices,2));
        roiCoords(2,:) = indices;
    end
end

return


text:
function pn = line_par_point_near_2d ( f, g, x0, y0, p )

%*****************************************************************************80
%
%% LINE_PAR_POINT_NEAR_2D: nearest point on parametric line to given point, 2D.
%
%  Discussion:
%
%    The parametric form of a line in 2D is:
%
%      X = X0 + F * T
%      Y = Y0 + G * T
%
%    We may normalize by choosing F*F + G*G = 1, and F nonnegative.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer, John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983,
%    ISBN: 0408012420.
%
%  Parameters:
%
%    Input, real F, G, X0, Y0, the parametric line parameters.
%
%    Input, real P(2,1), the point whose distance from the line is
%    to be measured.
%
%    Output, real PN(2,1), the point on the parametric line which
%    is nearest to P.
%
  t = ( f * ( p(1) - x0 ) + g * ( p(2) - y0 ) ) / ( f * f + g * g );

  pn(1,1) = x0 + t * f;
  pn(2,1) = y0 + t * g;

  return
end


text:
function addCERRtoPath(pathStr)
% function addCERRtoPath(pathStr)
%
% This function adds sub-directories from CERR distribution to Matlab/Octave path
% pathStr is the absolute path of CERR.
%
% APA, 7/30/2020

pathC = genpath(pathStr);
pathC = strsplit(pathC,';');
indGitV = cellfun(@isempty,strfind(pathC,{'.git'}));
pathC = pathC(indGitV);
addpath(strjoin(pathC,';'))

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function D = spm_eeg_avgfreq(S)
% Average a TF-dataset over frequency to get a time-domain dataset
% FORMAT D = spm_eeg_avgfreq(S)
%
% S        - input struct
%  fields of S:
%   D        - MEEG object or filename of M/EEG mat-file with epoched data
%   freqwin  - frequency window to average over [default: [-Inf, Inf]]
%   prefix   - prefix for the output file [default: 'P']
%
% Output:
% D        - MEEG object
%__________________________________________________________________________
% Copyright (C) 2012-2017 Wellcome Trust Centre for Neuroimaging

% Vladimir Litvak
% $Id: spm_eeg_avgfreq.m 7132 2017-07-10 16:22:58Z guillaume $

SVNrev = '$Rev: 7132 $';

%-Startup
%--------------------------------------------------------------------------
spm('FnBanner', mfilename, SVNrev);
spm('FigName','Average over frequency'); spm('Pointer','Watch');

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

D = spm_eeg_load(S.D);

if ~strncmpi(D.transformtype,'TF',2)
    error('This function only works on TF datasets.');
end

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

%-Generate new MEEG object with new files
%--------------------------------------------------------------------------
Dnew = clone(D, [S.prefix fname(D)], [D.nchannels D.nsamples D.ntrials]);

%-Averaged data
%--------------------------------------------------------------------------
spm_progress_bar('Init', D.ntrials, 'Trials processed');
if D.ntrials > 100, Ibar = floor(linspace(1, D.ntrials, 100));
else Ibar = 1:D.ntrials; end

for i = 1:D.ntrials

    Dnew(:, :, i) = spm_squeeze(mean(D(:, freqind, :, i), 2), 2);

    if any(Ibar == i), spm_progress_bar('Set', i); end
end

spm_progress_bar('Clear');

%-Save the new M/EEG dataset
%--------------------------------------------------------------------------
Dnew = Dnew.history(mfilename, S);
save(Dnew);

D = Dnew;

%-Cleanup
%--------------------------------------------------------------------------
fprintf('%-40s: %30s\n','Completed',spm('time'));                       %-#
spm('FigName','Average over frequency: done'); spm('Pointer','Arrow');


text:
% comp_log_Mel receives a time domain speech signal and produces its
% log Mel filterbank coefficients. It is the Matlab counterpart of the
% feature extraction program provided with AURORA2 database.
% Author: Xiao Xiong
% Created: 18 Jul, 2005
% Last Modified: 28 Jul, 2005
% Inputs:
%       x   1-D time domain signal
% Outputs:
%       hist_log_mel    the log Mel filterbank coefficients
%       hist_mel        the Mel filterbank coefficients
%       hist_abs        the spectral magnitude
%       logE            the log energy item needed in computing MFCC

function [hist_log_mel, hist_mel, hist_abs, hist_fft,hist_x, logE] = comp_log_Mel(x);

x(:);   % convert x to a column vector

FFT_length = 256;
logE_floor = -50;           % minimum value for log energy
Fs = 8000;              % for AURORA2 database, the sampling rate is 8k
N_mel = 23;             % for AURORA2 database, 23 mel filterbank are used

frame_size = Fs * 0.025;    % 25ms frame
frame_shift = Fs * 0.01;    % 10ms shift
frame_overlap = frame_size - frame_shift;

x_old = x;
% DC offset removing
x = DC_remove(x,0.999);
% pre-emphasis, boost the high frequency spectrum
A    = [1 -0.97];
x = filter(A,1,x);    % y(n) = x(n) - 0.97*x(n-1)
% number of blocks
N_block = floor((length(x)-frame_size)/frame_shift)+1;
% generate the mel window
mel_win = mel_window_FE(N_mel, FFT_length/2, Fs);

% produce the hamming windowm
window = hamming(frame_size);

    hist_x = zeros(frame_size,N_block);
    hist_fft = zeros(FFT_length,N_block);
    hist_abs = zeros(FFT_length/2+1,N_block);
    hist_mel = zeros(N_mel,N_block);
    hist_log_mel = zeros(N_mel,N_block);

for i = 1:N_block
%     if mod(i,100) == 0, disp(i); end;
    % step 1. framing
    start = (i-1)*frame_shift+1;
    last = min(length(x),(i-1)*frame_shift+frame_size);
    x_fr = x(start:last);
    if i<N_block, hist_x(:,i) = x_fr; end
    % step 2. calculate the log energy
    logE(i) = log(max(logE_floor, sum(x_old(start:last).^2)));
    % step3. windowing
    x_fr = x_fr(:).*window;
    % step 4. zero padding, if the number of elements in x_fr is less than
    % the length of FFT, append zeros to its end
    x_fr = [x_fr' zeros(1,FFT_length-length(x_fr))];
    % step 5, calculate the Fourier transform using Fast Fourier Transform
    X = fft(x_fr(:),FFT_length);
    % step 6. extract the magnitude of the spectral coefficients
    X_abs = abs(X(2:FFT_length/2+1));
    % step 7. mel-scale window wraping
    X_mel = X_abs'*mel_win;
    % step 8. take the natural logarithm
    X_log_mel = log(X_mel);

    hist_fft(:,i) = X;
    hist_abs(:,i) = abs(X(1:FFT_length/2+1));
    hist_mel(:,i) = X_mel;
    hist_log_mel(:,i) = X_log_mel;
end
%figure;
%imageFE(hist_log_mel');
a=1;

author: singaxiong

repo: SignalGraph


text:
function [UV,F, res, edge_norms] = ...
  create_regular_grid(xRes, yRes, xWrap, yWrap, near, far)
% Creates list of triangle vertex indices for a rectangular domain,
% optionally wrapping around in X/Y direction.
%
% Usage:
%   [UV,F,res,edge_norms] = create_regular_grid(xRes, yRes, xWrap, yWrap)
%
% Input:
%   xRes, yRes: number of points in X/Y direction
%   wrapX, wrapY: wrap around in X/Y direction
%   near, far: near and far should be fractions of one which control the
%              pinching of the domain at the center and sides
%
% Output:
%   F : mesh connectivity (triangles)
%   UV: UV coordinates in interval [0,1]x[0,1]
%   res: mesh resolution
%
% Example:
%  % Create and m by n cylinder
%  m = 10; n = 20;
%  [V,F] = create_regular_grid(m,n,1,0);
%  V = [sin(2*pi*V(:,1)) cos(2*pi*V(:,1)) (n-1)*2*pi/(m-1)*V(:,2)];
%  tsurf(F,V); axis equal;
%
%  % Quads:
%  Q = [F(1:2:end-1,[1 2]) F(2:2:end,[2 3])];
%
%  % Alternating diagonals:
%  Q(1:2:end,:) = Q(1:2:end,[2 3 4 1]);
%  F = [Q(:,[1 2 3]);Q(:,[1 3 4])];
%

if (nargin<2) yRes=xRes; end
if (nargin<3) xWrap=0; end
if (nargin<4) yWrap=0; end
if (nargin<5) overlap=0; end

%res = [yRes, xRes];
res_wrap = [yRes+yWrap, xRes+xWrap];

%xSpace = linspace(0,1,xRes+xWrap); if (xWrap) xSpace = xSpace(1:end-1); end
%ySpace = linspace(0,1,yRes+yWrap); if (yWrap) ySpace = ySpace(1:end-1); end
xSpace = linspace(0,1,xRes+xWrap);
ySpace = linspace(0,1,yRes+yWrap);

[X, Y] = meshgrid(xSpace, ySpace);
UV_wrap = [X(:), Y(:)];

% Must perform pinch before edge_norms are taken
if(exist('near') & exist('far'))
  if(near>0 & far>0)
  t = ( ...
      UV_wrap(:,1).*(UV_wrap(:,1)<0.5)+ ...
      (1-UV_wrap(:,1)).*(UV_wrap(:,1)>=0.5) ...
    )/0.5;
  t = 1-sin(t*pi/2+pi/2);
  UV_wrap(:,2) = ...
    far/2 + ...
    near*(UV_wrap(:,2)-0.5).*(1-t) + ...
    far*(UV_wrap(:,2)-0.5).*t;
  else
    %error('Pinch must be between 0 and 1');
  end
end

idx_wrap = reshape(1:prod(res_wrap), res_wrap);

v1_wrap = idx_wrap(1:end-1, 1:end-1); v1_wrap=v1_wrap(:)';
v2_wrap = idx_wrap(1:end-1, 2:end  ); v2_wrap=v2_wrap(:)';
v3_wrap = idx_wrap(2:end  , 1:end-1); v3_wrap=v3_wrap(:)';
v4_wrap = idx_wrap(2:end  , 2:end  ); v4_wrap=v4_wrap(:)';

F_wrap = [v1_wrap;v2_wrap;v3_wrap; v2_wrap;v4_wrap;v3_wrap];
F_wrap = reshape(F_wrap, [3, 2*length(v1_wrap)])';

% old way
% edges = [F_wrap(:,1) F_wrap(:,2); F_wrap(:,2) F_wrap(:,3); F_wrap(:,3) F_wrap(:,1)];
% edge_norms = sqrt(sum((UV_wrap(edges(:,1),:)-UV_wrap(edges(:,2),:)).^2,2));
% edge_norms = reshape(edge_norms,size(F_wrap,1),3);

% edges numbered same as opposite vertices
edge_norms = [ ...
  sqrt(sum((UV_wrap(F_wrap(:,2),:)-UV_wrap(F_wrap(:,3),:)).^2,2)) ...
  sqrt(sum((UV_wrap(F_wrap(:,3),:)-UV_wrap(F_wrap(:,1),:)).^2,2)) ...
  sqrt(sum((UV_wrap(F_wrap(:,1),:)-UV_wrap(F_wrap(:,2),:)).^2,2)) ...
  ];

% correct indices
res = [yRes,xRes];
idx = reshape(1:prod(res),res);
if (xWrap) idx = [idx, idx(:,1)]; end
if (yWrap) idx = [idx; idx(1,:)]; end
idx_flat = idx(:);

% this might not be neccessary, could just rebuild UV like before
UV = reshape(UV_wrap,[size(idx_wrap),2]);
UV = UV(1:end-yWrap,1:end-xWrap,:);
UV = reshape(UV,xRes*yRes,2);

F = [idx_flat(F_wrap(:,1)),idx_flat(F_wrap(:,2)),idx_flat(F_wrap(:,3))];

author: alecjacobson

repo: gptoolbox


text:
function md5=mocov_util_md5(filename)
% return the md5sum from the contents of a file
%
% md5=mocov_util_md5(filename)
%
% Inputs:
%   filename            name of file for which md5 must be computed
%
% Output:
%   md5                 md5 checksum of the contents of the file filename
%
% Notes:
%    - this function requires either the presence of a function named
%      'md5sum' (available on the GNU Octave platform), or the usage of a
%      unix-like platform
%
%
    md5_with_whitespace=md5_from_file(filename);
    md5=regexprep(md5_with_whitespace,'\s','');%

function md5=md5_from_file(fn)
    md5_processors={@md5_builtin,...
                    @hash_builtin,...
                    @md5sum_shell,...
                    @md5_shell};

    n=numel(md5_processors);
    for k=1:n
        md5_processor=md5_processors{k};
        [is_ok,md5]=md5_processor(fn);
        if is_ok
            return;
        end
    end

    error('Unable to compute md5 - no method available');

function [is_ok,md5]=md5_builtin(fn)
% supported in GNU Octave <= 4.6
    is_ok=has_builtin_function('md5sum');
    if is_ok
        md5=md5sum(fn);
    else
        md5=[];
    end

function [is_ok,md5]=hash_builtin(fn)
% supported in GNU Octave >= 4.4

    is_ok=has_builtin_function('hash') && has_builtin_function('fileread');
    if is_ok
        md5=hash('md5',fileread(fn));
    else
        md5=[];
    end

function tf=has_builtin_function(name)
    tf=exist(name,'builtin');

function [is_ok,md5]=run_unix(cmd)
% helper function
    is_ok=false;
    md5=[];

    if ispc()
        return;
    end

    [status,md5]=unix(cmd);
    is_ok=status==0;

function [is_ok,md5]=md5sum_shell(fn)
% supported on Unix platform

    cmd=sprintf('md5sum "%s"',fn);
    [is_ok,md5_with_fn]=run_unix(cmd);
    parts=regexp(md5_with_fn,'\s+','split');
    md5=parts{1};

function [is_ok,md5]=md5_shell(fn)
% supported on Unix platform

    cmd=sprintf('md5 -q "%s"',fn);
    [is_ok,md5]=run_unix(cmd);

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function realignunwarp = spm_cfg_realignunwarp
% SPM Configuration file for Realign & Unwarp
%__________________________________________________________________________
% Copyright (C) 2005-2016 Wellcome Trust Centre for Neuroimaging

% $Id: spm_cfg_realignunwarp.m 6952 2016-11-25 16:03:13Z guillaume $

%--------------------------------------------------------------------------
% scans Images
%--------------------------------------------------------------------------
scans         = cfg_files;
scans.tag     = 'scans';
scans.name    = 'Images';
scans.help    = {
    'Select scans for this session.'
    'In the coregistration step, the sessions are first realigned to each other, by aligning the first scan from each session to the first scan of the first session.  Then the images within each session are aligned to the first image of the session. The parameter estimation is performed this way because it is assumed (rightly or not) that there may be systematic differences in the images between sessions.'
    }';
scans.filter  = 'image';
scans.ufilter = '.*';
scans.num     = [1 Inf];
scans.preview = @(f) spm_check_registration(char(f));

%--------------------------------------------------------------------------
% pmscan Phase map (vdm* file)
%--------------------------------------------------------------------------
pmscan         = cfg_files;
pmscan.tag     = 'pmscan';
pmscan.name    = 'Phase map (vdm* file)';
pmscan.help    = {
    'Select pre-calculated phase map, or leave empty for no phase correction.'
    'The vdm* file is assumed to be already in alignment with the first scan of the first session.'
    }';
pmscan.filter  = 'image';
pmscan.ufilter = '^vdm5_.*';
pmscan.num     = [0 1];
pmscan.val     = {''};
pmscan.preview = @(f) spm_image('Display',char(f));

%--------------------------------------------------------------------------
% data Session
%--------------------------------------------------------------------------
data         = cfg_branch;
data.tag     = 'data';
data.name    = 'Session';
data.val     = {scans pmscan};
data.help    = {
    'Only add similar session data to a realign+unwarp branch, i.e., choose Data or Data+phase map for all sessions, but don''t use them interchangeably.'
    ''
    'In the coregistration step, the sessions are first realigned to each other, by aligning the first scan from each session to the first scan of the first session.  Then the images within each session are aligned to the first image of the session. The parameter estimation is performed this way because it is assumed (rightly or not) that there may be systematic differences in the images between sessions.'
    }';

%--------------------------------------------------------------------------
% generic Data
%--------------------------------------------------------------------------
generic         = cfg_repeat;
generic.tag     = 'generic';
generic.name    = 'Data';
generic.help    = {'Data sessions to unwarp.'};
generic.values  = {data};
generic.num     = [1 Inf];

%--------------------------------------------------------------------------
% quality Quality
%--------------------------------------------------------------------------
quality         = cfg_entry;
quality.tag     = 'quality';
quality.name    = 'Quality';
quality.help    = {
    'Quality versus speed trade-off.'
    'Highest quality (1) gives most precise results, whereas lower qualities gives faster realignment. The idea is that some voxels contribute little to the estimation of the realignment parameters. This parameter is involved in selecting the number of voxels that are used.'
    }';
quality.strtype = 'r';
quality.num     = [1 1];
quality.extras  = [0 1];
quality.def     = @(val)spm_get_defaults('realign.estimate.quality', val{:});

%--------------------------------------------------------------------------
% sep Separation
%--------------------------------------------------------------------------
sep         = cfg_entry;
sep.tag     = 'sep';
sep.name    = 'Separation';
sep.help    = {
    'The separation (in mm) between the points sampled in the reference image.'
    'Smaller sampling distances gives more accurate results, but will be slower.'};
sep.strtype = 'r';
sep.num     = [1 1];
sep.def     = @(val)spm_get_defaults('realign.estimate.sep', val{:});

%--------------------------------------------------------------------------
% fwhm Smoothing (FWHM)
%--------------------------------------------------------------------------
fwhm         = cfg_entry;
fwhm.tag     = 'fwhm';
fwhm.name    = 'Smoothing (FWHM)';
fwhm.help    = {
    'The FWHM of the Gaussian smoothing kernel (mm) applied to the images before estimating the realignment parameters.'
    ''
    '    * PET images typically use a 7 mm kernel.'
    ''
    '    * MRI images typically use a 5 mm kernel.'
    }';
fwhm.strtype = 'r';
fwhm.num     = [1 1];
fwhm.def     = @(val)spm_get_defaults('realign.estimate.fwhm', val{:});

%--------------------------------------------------------------------------
% rtm Num Passes
%--------------------------------------------------------------------------
rtm         = cfg_menu;
rtm.tag     = 'rtm';
rtm.name    = 'Num Passes';
rtm.help    = {
    'Register to first: Images are registered to the first image in the series. Register to mean: A two pass procedure is used in order to register the images to the mean of the images after the first realignment.'
    ''
    '    * PET images are typically registered to the mean.'
    ''
    '    * MRI images are typically registered to the first image.'
    }';
rtm.labels  = {
               'Register to first'
               'Register to mean'
}';
rtm.values  = {0 1};
rtm.def     = @(val)spm_get_defaults('unwarp.estimate.rtm', val{:});

%--------------------------------------------------------------------------
% einterp Interpolation
%--------------------------------------------------------------------------
einterp         = cfg_menu;
einterp.tag     = 'einterp';
einterp.name    = 'Interpolation';
einterp.help    = {
    'The method by which the images are sampled when estimating the optimum transformation.'
    'Higher degree interpolation methods provide the better interpolation, but they are slower because they use more neighbouring voxels /* \cite{thevenaz00a,unser93a,unser93b}*/.'
    }';
einterp.labels  = {
                   'Nearest neighbour'
                   'Trilinear'
                   '2nd Degree B-spline'
                   '3rd Degree B-Spline'
                   '4th Degree B-Spline'
                   '5th Degree B-Spline'
                   '6th Degree B-Spline'
                   '7th Degree B-Spline'
}';
einterp.values  = {0 1 2 3 4 5 6 7};
einterp.def     = @(val)spm_get_defaults('realign.estimate.interp', val{:});

%--------------------------------------------------------------------------
% ewrap Wrapping
%--------------------------------------------------------------------------
ewrap         = cfg_menu;
ewrap.tag     = 'ewrap';
ewrap.name    = 'Wrapping';
ewrap.help    = {
    'This indicates which directions in the volumes the values should wrap around in.'
    'These are typically:'
    '* No wrapping - for images that have already been spatially transformed.'
    '* Wrap in Y   - for (un-resliced) MRI where phase encoding is in the Y direction (voxel space).'
    }';
ewrap.labels = {
                'No wrap'
                'Wrap X'
                'Wrap Y'
                'Wrap X & Y'
                'Wrap Z '
                'Wrap X & Z'
                'Wrap Y & Z'
                'Wrap X, Y & Z'
}';
ewrap.values = {[0 0 0] [1 0 0] [0 1 0] [1 1 0] [0 0 1] [1 0 1] [0 1 1]...
                [1 1 1]};
ewrap.def    = @(val)spm_get_defaults('realign.estimate.wrap', val{:});

%--------------------------------------------------------------------------
% weight Weighting
%--------------------------------------------------------------------------
weight         = cfg_files;
weight.tag     = 'weight';
weight.name    = 'Weighting';
weight.val     = {''};
weight.help    = {
    'Optional weighting image to weight each voxel of the reference image differently when estimating the realignment parameters.'
    'The weights are proportional to the inverses of the standard deviations. For example, when there is a lot of extra-brain motion - e.g., during speech, or when there are serious artifacts in a particular region of the images.'
    }';
weight.filter  = 'image';
weight.ufilter = '.*';
weight.num     = [0 1];
weight.preview = @(f) spm_image('Display',char(f));

%--------------------------------------------------------------------------
% eoptions Estimation Options
%--------------------------------------------------------------------------
eoptions         = cfg_branch;
eoptions.tag     = 'eoptions';
eoptions.name    = 'Estimation Options';
eoptions.val     = {quality sep fwhm rtm einterp ewrap weight };
eoptions.help    = {'Various registration options that could be modified to improve the results. Whenever possible, the authors of SPM try to choose reasonable settings, but sometimes they can be improved.'};

%--------------------------------------------------------------------------
% basfcn Basis Functions
%--------------------------------------------------------------------------
basfcn         = cfg_menu;
basfcn.tag     = 'basfcn';
basfcn.name    = 'Basis Functions';
basfcn.help    = {
    'Number of basis functions to use for each dimension.'
    'If the third dimension is left out, the order for that dimension is calculated to yield a roughly equal spatial cut-off in all directions. Default: [12 12 *]'
    }';
basfcn.labels  = {
                  '8x8x*'
                  '10x10x*'
                  '12x12x*'
                  '14x14x*'
}';
basfcn.values  = {[8 8] [10 10] [12 12] [14 14]};
basfcn.def     = @(val)spm_get_defaults('unwarp.estimate.basfcn', val{:});

%--------------------------------------------------------------------------
% regorder Regularisation
%--------------------------------------------------------------------------
regorder         = cfg_menu;
regorder.tag     = 'regorder';
regorder.name    = 'Regularisation';
regorder.help    = {
    'Unwarp looks for the solution that maximises the likelihood (minimises the variance) while simultaneously maximising the smoothness of the estimated field (c.f. Lagrange multipliers). This parameter determines how to balance the compromise between these (i.e. the value of the multiplier). Test it on your own data (if you can be bothered) or go with the defaults. '
    ''
    'Regularisation of derivative fields is based on the regorder''th (spatial) derivative of the field. The choices are 0, 1, 2, or 3.  Default: 1'
    }';
regorder.labels = {
                   '0'
                   '1'
                   '2'
                   '3'
}';
regorder.values = {0 1 2 3};
regorder.def    = @(val)spm_get_defaults('unwarp.estimate.regorder', val{:});

%--------------------------------------------------------------------------
% lambda Reg. Factor
%--------------------------------------------------------------------------
lambda         = cfg_menu;
lambda.tag     = 'lambda';
lambda.name    = 'Reg. Factor';
lambda.help    = {'Regularisation factor. Default: Medium.'};
lambda.labels  = {
                 'A little'
                 'Medium'
                 'A lot'
}';
lambda.values  = {10000 100000 1000000};
lambda.def     = @(val)spm_get_defaults('unwarp.estimate.regwgt', val{:});

%--------------------------------------------------------------------------
% jm Jacobian deformations
%--------------------------------------------------------------------------
jm         = cfg_menu;
jm.tag     = 'jm';
jm.name    = 'Jacobian deformations';
jm.help    = {
    'Option to include Jacobian intensity modulation when estimating the fields.'
    '"Jacobian intensity modulation" refers to the dilution/concentration of intensity that ensue as a consequence of the distortions. Think of a semi-transparent coloured rubber sheet that you hold against a white background. If you stretch a part of the sheet (induce distortions) you will see the colour fading in that particular area. In theory it is a brilliant idea to include also these effects when estimating the field (see e.g. Andersson et al, NeuroImage 20:870-888). In practice for this specific problem it is NOT a good idea. Default: No'
    }';
jm.labels  = {
              'Yes'
              'No'
}';
jm.values  = {1 0};
jm.def     = @(val)spm_get_defaults('unwarp.estimate.jm', val{:});

%--------------------------------------------------------------------------
% fot First-order effects
%--------------------------------------------------------------------------
fot         = cfg_entry;
fot.tag     = 'fot';
fot.name    = 'First-order effects';
fot.help    = {
    'Vector of first order effects to model.'
    ''
    'Theoretically (ignoring effects of shimming) one would expect the field to depend only on subject out-of-plane rotations. Hence the default choice ("Pitch and Roll", i.e., [4 5]). Go with that unless you have very good reasons to do otherwise'
    ''
    'Movements to be modelled are referred to by number. 1= x translation; 2= y translation; 3= z translation 4 = x rotation,  5 = y rotation and 6 = z rotation.'
    ''
    'To model pitch & roll enter: [4 5]'
    ''
    'To model all movements enter: [1:6]'
    ''
    'Otherwise enter a customised set of movements to model'
    }';
fot.strtype = 'n';
fot.num     = [1 Inf];
fot.def     = @(val)spm_get_defaults('unwarp.estimate.foe', val{:});

%--------------------------------------------------------------------------
% sot Second-order effects
%--------------------------------------------------------------------------
sot         = cfg_entry;
sot.tag     = 'sot';
sot.name    = 'Second-order effects';
sot.help    = {
    'List of second order terms to model second derivatives of.'
    'This is entered as  a vector of movement parameters similar to first order effects, or leave blank for NONE.'
    ''
    'Movements to be modelled are referred to by number:'
    ''
    '1= x translation; 2= y translation; 3= z translation 4 = x rotation,  5 = y rotation and 6 = z rotation.'
    ''
    'To model the interaction of pitch & roll enter: [4 5]'
    ''
    'To model all movements enter: [1:6]'
    ''
    'The vector will be expanded into an n x 2 matrix of effects. For example [4 5] will be expanded to:'
    ''
    '[ 4 4'
    ''
    '  4 5'
    ''
    '  5 5 ]'
    }';
sot.strtype = 'n';
sot.num     = [Inf Inf];
sot.def     = @(val)spm_get_defaults('unwarp.estimate.soe', val{:});

%--------------------------------------------------------------------------
% uwfwhm Smoothing for unwarp (FWHM)
%--------------------------------------------------------------------------
uwfwhm         = cfg_entry;
uwfwhm.tag     = 'uwfwhm';
uwfwhm.name    = 'Smoothing for unwarp (FWHM)';
uwfwhm.help    = {'FWHM (mm) of smoothing filter applied to images prior to estimation of deformation fields.'};
uwfwhm.strtype = 'r';
uwfwhm.num     = [1 1];
uwfwhm.def     = @(val)spm_get_defaults('unwarp.estimate.fwhm', val{:});

%--------------------------------------------------------------------------
% rem Re-estimate movement params
%--------------------------------------------------------------------------
rem         = cfg_menu;
rem.tag     = 'rem';
rem.name    = 'Re-estimate movement params';
rem.help    = {'Re-estimation means that movement-parameters should be re-estimated at each unwarping iteration.'};
rem.labels  = {
              'Yes'
              'No'
}';
rem.values  = {1 0};
rem.def     = @(val)spm_get_defaults('unwarp.estimate.rem', val{:});

%--------------------------------------------------------------------------
% noi Number of Iterations
%--------------------------------------------------------------------------
noi         = cfg_entry;
noi.tag     = 'noi';
noi.name    = 'Number of Iterations';
noi.help    = {'Maximum number of iterations.'};
noi.strtype = 'n';
noi.num     = [1 1];
noi.def     = @(val)spm_get_defaults('unwarp.estimate.noi', val{:});

%--------------------------------------------------------------------------
% expround Taylor expansion point
%--------------------------------------------------------------------------
expround         = cfg_menu;
expround.tag     = 'expround';
expround.name    = 'Taylor expansion point';
expround.help    = {
    'Point in position space to perform Taylor-expansion around.'
    'Choices are (''First'', ''Last'' or ''Average''). ''Average'' should (in principle) give the best variance reduction. If a field-map acquired before the time-series is supplied then expansion around the ''First'' MIGHT give a slightly better average geometric fidelity.'
    };
expround.labels = {
                   'Average'
                   'First'
                   'Last'
}';
expround.values = {
                   'Average'
                   'First'
                   'Last'
}';
expround.def     = @(val)spm_get_defaults('unwarp.estimate.expround', val{:});

%--------------------------------------------------------------------------
% uweoptions Unwarp Estimation Options
%--------------------------------------------------------------------------
uweoptions      = cfg_branch;
uweoptions.tag  = 'uweoptions';
uweoptions.name = 'Unwarp Estimation Options';
uweoptions.val  = {basfcn regorder lambda jm fot sot uwfwhm rem noi expround};
uweoptions.help = {'Various registration & unwarping estimation options.'};

%--------------------------------------------------------------------------
% uwwhich Reslices images (unwarp)?
%--------------------------------------------------------------------------
uwwhich         = cfg_menu;
uwwhich.tag     = 'uwwhich';
uwwhich.name    = 'Resliced images (unwarp)?';
uwwhich.help    = {
    'Specify the images to reslice.'
    'All Images (1..n) '
    '  This reslices and unwarps all the images. '
    ' '
    'All Images + Mean Image '
    '  In addition to reslicing the images, it also creates a mean of the resliced images.'
    }';
uwwhich.labels = {
                  ' All Images (1..n)'
                  ' All Images + Mean Image'
}';
uwwhich.values = {[2 0] [2 1]};
uwwhich.def    = @(val)spm_get_defaults('realign.write.which', val{:});

%--------------------------------------------------------------------------
% rinterp Interpolation
%--------------------------------------------------------------------------
rinterp         = cfg_menu;
rinterp.tag     = 'rinterp';
rinterp.name    = 'Interpolation';
rinterp.help    = {
    'The method by which the images are sampled when being written in a different space.'
    'Nearest Neighbour is fastest, but not recommended for image realignment. Trilinear Interpolation is probably OK for PET, but not so suitable for fMRI because higher degree interpolation generally gives better results/* \cite{thevenaz00a,unser93a,unser93b}*/. Although higher degree methods provide better interpolation, but they are slower because they use more neighbouring voxels.'
    }';
rinterp.labels = {
                  'Nearest neighbour'
                  'Trilinear'
                  '2nd Degree B-spline '
                  '3rd Degree B-Spline'
                  '4th Degree B-Spline'
                  '5th Degree B-Spline '
                  '6th Degree B-Spline'
                  '7th Degree B-Spline'
}';
rinterp.values = {0 1 2 3 4 5 6 7};
rinterp.def    = @(val)spm_get_defaults('realign.write.interp', val{:});

%--------------------------------------------------------------------------
% wrap Wrapping
%--------------------------------------------------------------------------
wrap         = cfg_menu;
wrap.tag     = 'wrap';
wrap.name    = 'Wrapping';
wrap.help    = {
    'This indicates which directions in the volumes the values should wrap around in.'
    'These are typically:'
    '* No wrapping - for images that have already been spatially transformed.'
    '* Wrap in Y  - for (un-resliced) MRI where phase encoding is in the Y direction (voxel space).'
    }';
wrap.labels = {
               'No wrap'
               'Wrap X'
               'Wrap Y'
               'Wrap X & Y'
               'Wrap Z '
               'Wrap X & Z'
               'Wrap Y & Z'
               'Wrap X, Y & Z'
}';
wrap.values = {[0 0 0] [1 0 0] [0 1 0] [1 1 0] [0 0 1] [1 0 1] [0 1 1]...
               [1 1 1]};
wrap.def    = @(val)spm_get_defaults('realign.write.wrap', val{:});

%--------------------------------------------------------------------------
% mask Masking
%--------------------------------------------------------------------------
mask         = cfg_menu;
mask.tag     = 'mask';
mask.name    = 'Masking';
mask.help    = {'Because of subject motion, different images are likely to have different patterns of zeros from where it was not possible to sample data. With masking enabled, the program searches through the whole time series looking for voxels which need to be sampled from outside the original images. Where this occurs, that voxel is set to zero for the whole set of images (unless the image format can represent NaN, in which case NaNs are used where possible).'};
mask.labels  = {
                'Mask images'
                'Dont mask images'
}';
mask.values  = {1 0};
mask.def     = @(val)spm_get_defaults('realign.write.mask', val{:});

%--------------------------------------------------------------------------
% prefix Filename Prefix
%--------------------------------------------------------------------------
prefix         = cfg_entry;
prefix.tag     = 'prefix';
prefix.name    = 'Filename Prefix';
prefix.help    = {'Specify the string to be prepended to the filenames of the smoothed image file(s). Default prefix is ''u''.'};
prefix.strtype = 's';
prefix.num     = [1 Inf];
prefix.def     = @(val)spm_get_defaults('unwarp.write.prefix', val{:});

%--------------------------------------------------------------------------
% uwroptions Unwarp Reslicing Options
%--------------------------------------------------------------------------
uwroptions         = cfg_branch;
uwroptions.tag     = 'uwroptions';
uwroptions.name    = 'Unwarp Reslicing Options';
uwroptions.val     = {uwwhich rinterp wrap mask prefix};
uwroptions.help    = {'Various registration & unwarping estimation options.'};

%--------------------------------------------------------------------------
% realignunwarp Realign & Unwarp
%--------------------------------------------------------------------------
realignunwarp         = cfg_exbranch;
realignunwarp.tag     = 'realignunwarp';
realignunwarp.name    = 'Realign & Unwarp';
realignunwarp.val     = {generic eoptions uweoptions uwroptions};
realignunwarp.help    = {
    'Within-subject registration and unwarping of time series.'
    ''
    'The realignment part of this routine realigns a time-series of images acquired from the same subject using a least squares approach and a 6 parameter (rigid body) spatial transformation.  The first image in the list specified by the user is used as a reference to which all subsequent scans are realigned. The reference scan does not have to the the first chronologically and it may be wise to chose a "representative scan" in this role.'
    ''
    'The aim is primarily to remove movement artefact in fMRI and PET time-series (or more generally longitudinal studies). This affects the header of each of the input images which contains details about the voxel-to-world mapping. The details of the transformation are displayed in the results window as plots of translation and rotation. A set of realignment parameters are saved for each session, named rp_*.txt.'
    ''
    'In the coregistration step, the sessions are first realigned to each other, by aligning the first scan from each session to the first scan of the first session.  Then the images within each session are aligned to the first image of the session. The parameter estimation is performed this way because it is assumed (rightly or not) that there may be systematic differences in the images between sessions.'
    ''
    'See spm_uw_estimate.m for a detailed description of the implementation.'
    ''
    'Even after realignment there is considerable variance in fMRI time series that covary with, and is most probably caused by, subject movements/* \cite{ja_geometric}*/. It is also the case that this variance is typically large compared to experimentally induced variance. Anyone interested can include the estimated movement parameters as covariates in the design matrix, and take a look at an F-contrast encompassing those columns. It is quite dramatic. The result is loss of sensitivity, and if movements are correlated to task specificity. I.e. we may mistake movement induced variance for true activations. The problem is well known, and several solutions have been suggested. A quite pragmatic (and conservative) solution is to include the estimated movement parameters (and possibly squared) as covariates in the design matrix. Since we typically have loads of degrees of freedom in fMRI we can usually afford this. The problems occur when movements are correlated with the task, since the strategy above will discard "good" and "bad" variance alike (i.e. remove also "true" activations).'
    ''
    'The "covariate" strategy described above was predicated on a model where variance was assumed to be caused by "spin history" effects, but will work pretty much equally good/bad regardless of what the true underlying cause is. Others have assumed that the residual variance is caused mainly by errors introduced by the interpolation kernel in the resampling step of the realignment. One has tried to solve this through higher order resampling (huge Sinc kernels, or k-space resampling). Unwarp is based on a different hypothesis regarding the residual variance. EPI images are not particularly faithful reproductions of the object, and in particular there are severe geometric distortions in regions where there is an air-tissue interface (e.g. orbitofrontal cortex and the anterior medial temporal lobes). In these areas in particular the observed image is a severely warped version of reality, much like a funny mirror at a fair ground. When one moves in front of such a mirror ones image will distort in different ways and ones head may change from very elongated to seriously flattened. If we were to take digital snapshots of the reflection at these different positions it is rather obvious that realignment will not suffice to bring them into a common space.'
    ''
    'The situation is similar with EPI images, and an image collected for a given subject position will not be identical to that collected at another. We call this effect susceptibility-by-movement interaction. Unwarp is predicated on the assumption that the susceptibility-by-movement interaction is responsible for a sizable part of residual movement related variance.'
    ''
    'Assume that we know how the deformations change when the subject changes position (i.e. we know the derivatives of the deformations with respect to subject position). That means that for a given time series and a given set of subject movements we should be able to predict the "shape changes" in the object and the ensuing variance in the time series. It also means that, in principle, we should be able to formulate the inverse problem, i.e. given the observed variance (after realignment) and known (estimated) movements we should be able to estimate how deformations change with subject movement. We have made an attempt at formulating such an inverse model, and at solving for the "derivative fields". A deformation field can be thought of as little vectors at each position in space showing how that particular location has been deflected. A "derivative field" is then the rate of change of those vectors with respect to subject movement. Given these "derivative fields" we should be able to remove the variance caused by the susceptibility-by-movement interaction. Since the underlying model is so restricted we would also expect experimentally induced variance to be preserved. Our experiments have also shown this to be true.'
    ''
    'In theory it should be possible to estimate also the "static" deformation field, yielding an unwarped (to some true geometry) version of the time series. In practise that doesn''t really seem to work. Hence, the method deals only with residual movement related variance induced by the susceptibility-by-movement interaction. This means that the time-series will be undistorted to some "average distortion" state rather than to the true geometry. If one wants additionally to address the issue of anatomical fidelity one should combine Unwarp with a measured fieldmap.'
    ''
    'The description above can be thought of in terms of a Taylor expansion of the field as a function of subject movement. Unwarp alone will estimate the first (and optionally second, see below) order terms of this expansion. It cannot estimate the zeroth order term (the distortions common to all scans in the time series) since that doesn''t introduce (almost) any variance in the time series. The measured fieldmap takes the role of the zeroth order term. Refer to the FieldMap toolbox and the documents FieldMap.man and FieldMap_principles.man for a description of how to obtain fieldmaps in the format expected by Unwarp.'
    ''
    'If we think of the field as a function of subject movement it should in principle be a function of six variables since rigid body movement has six degrees of freedom. However, the physics of the problem tells us that the field should not depend on translations nor on rotation in a plane perpendicular to the magnetic flux. Hence it should in principle be sufficient to model the field as a function of out-of-plane rotations (i.e. pitch and roll). One can object to this in terms of the effects of shimming (object no longer immersed in a homogenous field) that introduces a dependence on all movement parameters. In addition SPM/Unwarp cannot really tell if the transversal slices it is being passed are really perpendicular to the flux or not. In practice it turns out thought that it is never (at least we haven''t seen any case) necessary to include more than Pitch and Roll. This is probably because the individual movement parameters are typically highly correlated anyway, which in turn is probably because most heads that we scan are attached to a neck around which rotations occur. On the subject of Taylor expansion we should mention that there is the option to use a second-order expansion (through the defaults) interface. This implies estimating also the rate-of-change w.r.t. to some movement parameter of the rate-of-change of the field w.r.t. some movement parameter (colloquially known as a second derivative). It can be quite interesting to watch (and it is amazing that it is possible) but rarely helpful/necessary.'
    ''
    'In the defaults there is also an option to include Jacobian intensity modulation when estimating the fields. "Jacobian intensity modulation" refers to the dilution/concentration of intensity that ensue as a consequence of the distortions. Think of a semi-transparent coloured rubber sheet that you hold against a white background. If you stretch a part of the sheet (induce distortions) you will see the colour fading in that particular area. In theory it is a brilliant idea to include also these effects when estimating the field (see e.g. Andersson et al, NeuroImage 20:870-888). In practice for this specific problem it is NOT a good idea.'
    ''
    'It should be noted that this is a method intended to correct data afflicted by a particular problem. If there is little movement in your data to begin with this method will do you little good. If on the other hand there is appreciable movement in your data (>1deg) it will remove some of that unwanted variance. If, in addition, movements are task related it will do so without removing all your "true" activations. The method attempts to minimise total (across the image volume) variance in the data set. It should be realised that while (for small movements) a rather limited portion of the total variance is removed, the susceptibility-by-movement interaction effects are quite localised to "problem" areas. Hence, for a subset of voxels in e.g. frontal-medial and orbitofrontal cortices and parts of the temporal lobes the reduction can be quite dramatic (>90). The advantages of using Unwarp will also depend strongly on the specifics of the scanner and sequence by which your data has been acquired. When using the latest generation scanners distortions are typically quite small, and distortion-by-movement interactions consequently even smaller. A small check list in terms of distortions is '
    'a) Fast gradients->short read-out time->small distortions '
    'b) Low field (i.e. <3T)->small field changes->small distortions '
    'c) Low res (64x64)->short read-out time->small distortions '
    'd) SENSE/SMASH->short read-out time->small distortions '
    'If you can tick off all points above chances are you have minimal distortions to begin with and Unwarp might not be of use to you.'
    }';
realignunwarp.prog     = @spm_run_realignunwarp;
realignunwarp.vout     = @vout_realignunwarp;
realignunwarp.modality = {'PET','FMRI'};

%==========================================================================
function dep = vout_realignunwarp(job)
for k=1:numel(job.data)
    cdep(1)            = cfg_dep;
    cdep(1).sname      = sprintf('Realignment Param File (Sess %d)', k);
    cdep(1).src_output = substruct('.','sess', '()',{k}, '.','rpfile');
    cdep(1).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});
    cdep(2)            = cfg_dep;
    cdep(2).sname      = sprintf('Unwarp Params File (Sess %d)', k);
    cdep(2).src_output = substruct('.','sess', '()',{k}, '.','dsfile');
    cdep(2).tgt_spec   = cfg_findspec({{'filter','any','strtype','e'}});
    if job.uwroptions.uwwhich(1) == 2
        cdep(3)            = cfg_dep;
        cdep(3).sname      = sprintf('Unwarped Images (Sess %d)', k);
        cdep(3).src_output = substruct('.','sess', '()',{k}, '.','uwrfiles');
        cdep(3).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
    end
    if k == 1
        dep = cdep;
    else
        dep = [dep cdep];
    end
end

if job.uwroptions.uwwhich(2)
    dep(end+1)          = cfg_dep;
    dep(end).sname      = 'Unwarped Mean Image';
    dep(end).src_output = substruct('.','meanuwr');
    dep(end).tgt_spec   = cfg_findspec({{'filter','image','strtype','e'}});
end