text:
function [sf,out2,out3] = symmetric_dirichlet_sym(in1,in2)
%SYMMETRIC_DIRICHLET_SYM
%    [SF,OUT2,OUT3] = SYMMETRIC_DIRICHLET_SYM(IN1,IN2)

%    This function was generated by the Symbolic Math Toolbox version 8.6.
%    24-Apr-2021 15:46:34

P1 = in2(:,1);
P2 = in2(:,2);
P3 = in2(:,3);
P4 = in2(:,4);
P5 = in2(:,5);
P6 = in2(:,6);
P7 = in2(:,7);
P8 = in2(:,8);
P9 = in2(:,9);
U1 = in1(:,1);
U2 = in1(:,2);
U3 = in1(:,3);
U4 = in1(:,4);
U5 = in1(:,5);
U6 = in1(:,6);
t2 = U1.*2.0;
t3 = U2.*2.0;
t4 = U3.*2.0;
t5 = U4.*2.0;
t6 = U5.*2.0;
t7 = U6.*2.0;
t8 = -P2;
t9 = -P3;
t10 = -P5;
t11 = -P6;
t12 = -P8;
t13 = -P9;
t15 = -U2;
t17 = -U3;
t20 = -U5;
t22 = -U6;
t24 = U1./2.0;
t25 = U2./2.0;
t26 = U3./2.0;
t27 = U4./2.0;
t28 = U5./2.0;
t29 = U6./2.0;
t14 = -t2;
t16 = -t3;
t18 = -t4;
t19 = -t5;
t21 = -t6;
t23 = -t7;
t30 = P1+t8;
t31 = P1+t9;
t32 = P4+t10;
t33 = P4+t11;
t34 = P7+t12;
t35 = P7+t13;
t36 = U1+t15;
t37 = U1+t17;
t38 = U4+t20;
t39 = U4+t22;
t40 = -t25;
t41 = -t26;
t42 = -t28;
t43 = -t29;
t44 = U2+U3+t14;
t45 = U5+U6+t19;
t46 = t30.^2;
t47 = t31.^2;
t48 = t32.^2;
t49 = t33.^2;
t50 = t34.^2;
t51 = t35.^2;
t52 = t2+t16;
t53 = t2+t18;
t54 = t5+t21;
t55 = t5+t23;
t56 = t36.^2;
t57 = t37.^2;
t58 = t38.^2;
t59 = t39.^2;
t66 = t24+t40;
t67 = t24+t41;
t68 = t25+t41;
t69 = t27+t42;
t70 = t27+t43;
t71 = t28+t43;
t72 = t30.*t31;
t73 = t32.*t33;
t74 = t34.*t35;
t75 = t36.*t37;
t76 = t38.*t39;
t77 = (t30.*t33)./2.0;
t78 = (t31.*t32)./2.0;
t79 = (t30.*t35)./2.0;
t80 = (t31.*t34)./2.0;
t81 = (t32.*t35)./2.0;
t82 = (t33.*t34)./2.0;
t83 = (t36.*t39)./2.0;
t84 = (t37.*t38)./2.0;
t60 = t46.*2.0;
t61 = t47.*2.0;
t62 = t48.*2.0;
t63 = t49.*2.0;
t64 = t50.*2.0;
t65 = t51.*2.0;
t85 = -t78;
t86 = -t80;
t87 = -t82;
t88 = -t84;
t89 = t56+t58;
t90 = t57+t59;
t91 = t46+t48+t50;
t92 = t47+t49+t51;
t93 = t75+t76;
t110 = t72+t73+t74;
t94 = t60+t62+t64;
t95 = t61+t63+t65;
t96 = t77+t85;
t97 = t79+t86;
t98 = t81+t87;
t99 = t83+t88;
t100 = t53.*t91;
t101 = t52.*t92;
t102 = t55.*t91;
t103 = t54.*t92;
t111 = t90.*t91;
t112 = t89.*t92;
t104 = t96.^2;
t105 = t97.^2;
t106 = t98.^2;
t107 = 1.0./t99.^2;
t108 = 1.0./t99.^3;
t113 = t94+t95;
t114 = t100+t101;
t115 = t102+t103;
t116 = t111+t112;
t109 = t107.^2;
t117 = t104+t105+t106;
t118 = 1.0./sqrt(t117);
t121 = t107.*t117;
t119 = (t91.*t118)./2.0;
t120 = (t92.*t118)./2.0;
t122 = t121+1.0;
t123 = (t94.*t118)./4.0;
t124 = (t95.*t118)./4.0;
t127 = (t100.*t118)./4.0;
t128 = (t101.*t118)./4.0;
t129 = (t102.*t118)./4.0;
t130 = (t103.*t118)./4.0;
t135 = t110.*t118;
t145 = (t113.*t118)./4.0;
t148 = (t114.*t118)./4.0;
t149 = (t115.*t118)./4.0;
t150 = (t116.*t118)./4.0;
t125 = -t123;
t126 = -t124;
t131 = -t127;
t132 = -t128;
t133 = -t129;
t134 = -t130;
t136 = -t135;
t137 = t135./2.0;
t147 = t93.*t135.*(-1.0./2.0);
t151 = t119.*t122;
t152 = t120.*t122;
t156 = t122.*t135.*(-1.0./2.0);
t173 = t66.*t108.*t117.*(t127-(t36.*t135)./2.0).*-2.0;
t174 = t67.*t108.*t117.*(t127-(t36.*t135)./2.0).*-2.0;
t175 = t68.*t108.*t117.*(t127-(t36.*t135)./2.0).*-2.0;
t176 = t66.*t108.*t117.*(t128-(t37.*t135)./2.0).*-2.0;
t177 = t67.*t108.*t117.*(t128-(t37.*t135)./2.0).*-2.0;
t178 = t68.*t108.*t117.*(t128-(t37.*t135)./2.0).*-2.0;
t179 = t70.*t108.*t117.*(t127-(t36.*t135)./2.0).*-2.0;
t180 = t71.*t108.*t117.*(t127-(t36.*t135)./2.0).*-2.0;
t181 = t69.*t108.*t117.*(t128-(t37.*t135)./2.0).*-2.0;
t182 = t71.*t108.*t117.*(t128-(t37.*t135)./2.0).*-2.0;
t183 = t67.*t108.*t117.*(t129-(t38.*t135)./2.0).*-2.0;
t184 = t68.*t108.*t117.*(t129-(t38.*t135)./2.0).*-2.0;
t185 = t66.*t108.*t117.*(t130-(t39.*t135)./2.0).*-2.0;
t186 = t68.*t108.*t117.*(t130-(t39.*t135)./2.0).*-2.0;
t187 = t69.*t108.*t117.*(t129-(t38.*t135)./2.0).*-2.0;
t188 = t70.*t108.*t117.*(t129-(t38.*t135)./2.0).*-2.0;
t189 = t71.*t108.*t117.*(t129-(t38.*t135)./2.0).*-2.0;
t190 = t69.*t108.*t117.*(t130-(t39.*t135)./2.0).*-2.0;
t191 = t70.*t108.*t117.*(t130-(t39.*t135)./2.0).*-2.0;
t192 = t71.*t108.*t117.*(t130-(t39.*t135)./2.0).*-2.0;
t193 = t67.*t108.*t117.*(t127-(t36.*t135)./2.0).*2.0;
t194 = t67.*t108.*t117.*(t128-(t37.*t135)./2.0).*2.0;
t195 = t71.*t108.*t117.*(t127-(t36.*t135)./2.0).*2.0;
t196 = t69.*t108.*t117.*(t128-(t37.*t135)./2.0).*2.0;
t197 = t71.*t108.*t117.*(t128-(t37.*t135)./2.0).*2.0;
t198 = t67.*t108.*t117.*(t129-(t38.*t135)./2.0).*2.0;
t199 = t69.*t108.*t117.*(t129-(t38.*t135)./2.0).*2.0;
t200 = t71.*t108.*t117.*(t129-(t38.*t135)./2.0).*2.0;
t201 = t69.*t108.*t117.*(t130-(t39.*t135)./2.0).*2.0;
t202 = t71.*t108.*t117.*(t130-(t39.*t135)./2.0).*2.0;
t203 = -t122.*(t135-t145);
t138 = -t137;
t139 = t36.*t137;
t140 = t37.*t137;
t141 = t38.*t137;
t142 = t39.*t137;
t143 = t44.*t137;
t144 = t45.*t137;
t146 = t93.*t137;
sf = -t122.*(t146-t150);
if nargout > 1
    t155 = t122.*t137;
    t157 = t125+t137;
    t158 = t126+t137;
    t163 = t136+t145;
    t166 = t147+t150;
    t153 = t119+t138;
    t154 = t120+t138;
    t159 = t131+t139;
    t160 = t132+t140;
    t161 = t133+t141;
    t162 = t134+t142;
    t164 = t143+t148;
    t165 = t144+t149;
    t171 = -t122.*(t123+t138);
    t172 = -t122.*(t124+t138);
    t219 = -t108.*t117.*(t146-t150);
    t220 = t108.*t117.*(t146-t150);
    out2 = [t122.*t164+t71.*t220.*2.0;t70.*t220.*-2.0-t122.*(t128-(t37.*t135)./2.0);t69.*t220.*2.0-t122.*(t127-(t36.*t135)./2.0);t122.*t165-t68.*t220.*2.0;t67.*t220.*2.0-t122.*(t130-(t39.*t135)./2.0);t66.*t220.*-2.0-t122.*(t129-(t38.*t135)./2.0)];
end
if nargout > 2
    t221 = t66.*t67.*t109.*t117.*(t146-t150).*-6.0;
    t222 = t66.*t68.*t109.*t117.*(t146-t150).*-6.0;
    t223 = t67.*t68.*t109.*t117.*(t146-t150).*-6.0;
    t224 = t66.*t69.*t109.*t117.*(t146-t150).*-6.0;
    t225 = t66.*t70.*t109.*t117.*(t146-t150).*-6.0;
    t226 = t67.*t69.*t109.*t117.*(t146-t150).*-6.0;
    t227 = t66.*t71.*t109.*t117.*(t146-t150).*-6.0;
    t228 = t67.*t70.*t109.*t117.*(t146-t150).*-6.0;
    t229 = t68.*t69.*t109.*t117.*(t146-t150).*-6.0;
    t230 = t67.*t71.*t109.*t117.*(t146-t150).*-6.0;
    t231 = t68.*t70.*t109.*t117.*(t146-t150).*-6.0;
    t232 = t68.*t71.*t109.*t117.*(t146-t150).*-6.0;
    t233 = t69.*t70.*t109.*t117.*(t146-t150).*-6.0;
    t234 = t69.*t71.*t109.*t117.*(t146-t150).*-6.0;
    t235 = t70.*t71.*t109.*t117.*(t146-t150).*-6.0;
    t236 = t66.*t67.*t109.*t117.*(t146-t150).*6.0;
    t237 = t67.*t68.*t109.*t117.*(t146-t150).*6.0;
    t238 = t66.*t69.*t109.*t117.*(t146-t150).*6.0;
    t239 = t66.*t71.*t109.*t117.*(t146-t150).*6.0;
    t240 = t67.*t70.*t109.*t117.*(t146-t150).*6.0;
    t241 = t68.*t69.*t109.*t117.*(t146-t150).*6.0;
    t242 = t68.*t71.*t109.*t117.*(t146-t150).*6.0;
    t243 = t69.*t70.*t109.*t117.*(t146-t150).*6.0;
    t244 = t70.*t71.*t109.*t117.*(t146-t150).*6.0;
    t167 = t122.*t153;
    t168 = t122.*t154;
    t204 = t66.*t108.*t117.*t164.*2.0;
    t205 = t67.*t108.*t117.*t164.*2.0;
    t206 = t68.*t108.*t117.*t164.*2.0;
    t207 = t69.*t108.*t117.*t164.*2.0;
    t208 = t70.*t108.*t117.*t164.*2.0;
    t209 = t66.*t108.*t117.*t165.*2.0;
    t210 = t67.*t108.*t117.*t165.*2.0;
    t211 = t69.*t108.*t117.*t165.*2.0;
    t212 = t70.*t108.*t117.*t165.*2.0;
    t213 = t71.*t108.*t117.*t165.*2.0;
    t245 = t173+t199+t238;
    t246 = t191+t194+t240;
    t248 = t156+t179+t196+t243;
    t249 = t156+t185+t198+t236;
    t250 = t193+t201+t219+t226;
    t251 = t176+t188+t220+t225;
    t169 = -t167;
    t170 = -t168;
    t214 = -t205;
    t215 = -t207;
    t216 = -t210;
    t217 = -t211;
    t218 = -t213;
    t252 = t178+t212+t219+t231;
    t253 = t200+t204+t219+t239;
    t247 = t206+t218+t242;
    t254 = t202+t214+t220+t230;
    t255 = t175+t217+t220+t241;
    out3 = [t203-t71.*t108.*t117.*t164.*4.0-t71.^2.*t109.*t117.*(t146-t150).*6.0;t172+t197+t208+t244;t171+t195+t215+t234;t247;t254;t253;t170+t197+t208+t244;t152-t70.^2.*t109.*t117.*(t146-t150).*6.0-t70.*t108.*t117.*(t128-(t37.*t135)./2.0).*4.0;t248;t252;t246;t251;t169+t195+t215+t234;t248;t151-t69.^2.*t109.*t117.*(t146-t150).*6.0+t69.*t108.*t117.*(t127-(t36.*t135)./2.0).*4.0;t255;t250;t245;t247;t252;t255;t203+t68.*t108.*t117.*t165.*4.0-t68.^2.*t109.*t117.*(t146-t150).*6.0;t172+t186+t216+t237;t171+t184+t209+t222;t254;t246;t250;t170+t186+t216+t237;t152-t67.^2.*t109.*t117.*(t146-t150).*6.0+t67.*t108.*t117.*(t130-(t39.*t135)./2.0).*4.0;t249;t253;t251;t245;t169+t184+t209+t222;t249;t151-t66.^2.*t109.*t117.*(t146-t150).*6.0-t66.*t108.*t117.*(t129-(t38.*t135)./2.0).*4.0];
end

author: alecjacobson

repo: gptoolbox


text:
function cmap = classpal

nClasses = 51;

hFig  = gcf;
hAxes = gca;

if strcmpi(hAxes.Visible, 'on')
    bgColor = hAxes.Color;
else
    bgColor = hFig.Color;
end

cmap = distinguishable_colors(nClasses, bgColor);

end

author: pglira


text:
function cube_exactness_test01 ( )

%*****************************************************************************80
%
%% CUBE_EXACTNESS_TEST01 tests product Gauss-Legendre rules for the Legendre 3D integral.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 August 2014
%
%  Author:
%
%    John Burkardt
%
  a(1:3) = -1.0;
  b(1:3) = +1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'CUBE_EXACTNESS_TEST01\n' );
  fprintf ( 1, '  Product Gauss-Legendre rules for the 3D Legendre integral.\n' );
  fprintf ( 1, '  Density function rho(x) = 1.\n' );
  fprintf ( 1, '  Region: -1 <= x <= +1.\n' );
  fprintf ( 1, '          -1 <= y <= +1.\n' );
  fprintf ( 1, '          -1 <= z <= +1.\n' );
  fprintf ( 1, '  Level: L\n' );
  fprintf ( 1, '  Exactness: 2*L+1\n' );
  fprintf ( 1, '  Order: N = (L+1)*(L+1)*(L+1)\n' );

  for l = 0 : 5

    nx = l + 1;
    ny = l + 1;
    nz = l + 1;
    n = nx * ny * nz;
    t = 2 * l + 1;

    [ x, y, z, w ] = legendre_3d_set ( a, b, nx, ny, nz );

    p_max = t + 1;
    legendre_3d_exactness ( a, b, n, x, y, z, w, p_max );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [image_type, principal_filename, secondary_filenames] = MimGuessFileType(file_path, image_filename, default_guess, reporting)
    % MimGuessFileType. Heuristically determines a file format
    %
    %
    %     Licence
    %     -------
    %     Part of the TD MIM Toolkit. https://github.com/tomdoel
    %     Author: Tom Doel, Copyright Tom Doel 2014.  www.tomdoel.com
    %     Distributed under the MIT licence. Please see website for details.
    %

    [~, name, ext] = fileparts(image_filename);
    image_filename_without_extension = fullfile(file_path, name);
    if strcmp(ext, '.mat')
        image_type = MimImageFileFormat.Matlab;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;

    elseif strcmp(ext, '.gipl')
        image_type = MimImageFileFormat.Gipl;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;

    elseif strcmp(ext, '.hdr')
        image_type = MimImageFileFormat.Analyze;
        principal_filename = {image_filename};
        secondary_filenames = {};
        if CoreDiskUtilities.FileExists(file_path, [name '.img'])
            secondary_filenames = [image_filename_without_extension '.img'];
        end
        return;

    elseif strcmp(ext, '.nii')
        image_type = MimImageFileFormat.Nifti;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;

    elseif strcmp(ext, '.img')
        hdr_filename = [name '.hdr'];
        nii_filename = [name '.nii'];
        if CoreDiskUtilities.FileExists(file_path, nii_filename)
            image_type = MimImageFileFormat.Nifti;
            principal_filename = {nii_filename};
            secondary_filenames = {image_filename};
            return;
        elseif CoreDiskUtilities.FileExists(file_path, hdr_filename)
            image_type = MimImageFileFormat.Analyze;
            principal_filename = {hdr_filename};
            secondary_filenames = {image_filename};
            return;
        end

    elseif strcmp(ext, '.isi')
        image_type = MimImageFileFormat.Isi;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;

    elseif strcmp(ext, '.v3d')
        image_type = MimImageFileFormat.V3d;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;

    elseif strcmp(ext, '.vmp')
        image_type = MimImageFileFormat.Vmp;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;

    elseif strcmp(ext, '.xif')
        image_type = MimImageFileFormat.Xif;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;

    elseif strcmp(ext, '.vtk')
        image_type = MimImageFileFormat.Vtk;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;

    elseif strcmp(ext, '.vff')
        image_type = MimImageFileFormat.MicroCT;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;

    elseif strcmp(ext, '.par') || strcmp(ext, '.rec')
        image_type = MimImageFileFormat.Par;
        par_filename = {[name '.par']};
        rec_filename = {[name '.rec']};
        par_found = CoreDiskUtilities.FileExists(file_path, par_filename);
        rec_found = CoreDiskUtilities.FileExists(file_path, rec_filename);
        if par_found
            if rec_found
                principal_filename = {par_filename};
                secondary_filenames = {fullfile(file_path, rec_filename)};
                return;
            else
                principal_filename = {par_filename};
                secondary_filenames = {};
                return;
            end
        else
            if rec_found
                principal_filename = {rec_filename};
                secondary_filenames = {};
                return;
            end
        end

    % For metaheader files (mhd/mha) we also fetch the filename of the
    % raw image data
    elseif strcmp(ext, '.mhd') || strcmp(ext, '.mha')
        image_type = MimImageFileFormat.Metaheader;
        [is_meta_header, raw_filename] = MimDiskUtilities.IsFileMetaHeader(fullfile(file_path, image_filename), reporting);
        if ~is_meta_header
            reporting.Error('MimGuessFileType:OpenMHDFileFailed', ['Unable to read metaheader file ' image_filename]);
        end
        principal_filename = {image_filename};
        secondary_filenames = {raw_filename}; % ToDo: what if no secondary file?
        return;

    % If a .raw file is selected, look for the corresponding .mha or
    % .mhd file. We thrown an exception if no file is found, it cannot
    % be loaded or the raw filename does not match the raw file we are
    % loading
    elseif strcmp(ext, '.raw')
        [principal_filename, secondary_filenames] = MimDiskUtilities.GetHeaderFileFromRawFile(file_path, name, reporting);
        if isempty(principal_filename)
            reporting.ShowWarning('MimGuessFileType:HeaderFileLoadError', ['Unable to find valid header file for ' fullfile(file_path, image_filename)], []);
        else
            if ~strcmp(secondary_filenames{1}, image_filename)
                reporting.Error('MimGuessFileType:MetaHeaderRawFileMismatch', ['Mismatch between specified image filename and entry in ' principal_filename{1}]);
            end
            image_type = MimImageFileFormat.Metaheader;
            return;
        end
    end

    % Unknown file type. Try looking for a header file
    [principal_filename_mh, secondary_filenames_mh] = MimDiskUtilities.GetHeaderFileFromRawFile(file_path, name, reporting);
    if (~isempty(principal_filename_mh)) && (strcmp(secondary_filenames_mh{1}, image_filename))
        image_type = MimImageFileFormat.Metaheader;
        principal_filename = principal_filename_mh;
        secondary_filenames = secondary_filenames_mh;
        return;
    end

    % Test for a DICOM image
    if DMUtilities.IsDicom(file_path, image_filename)
        image_type = MimImageFileFormat.Dicom;
        principal_filename = {image_filename};
        secondary_filenames = {};
        return;
    end

    % If all else fails, use the guess
    image_type = default_guess;
    principal_filename = {image_filename};
    secondary_filenames = {};
end

author: tomdoel

repo: pulmonarytoolkit


text:
function h = drawPoint3d(varargin)
%DRAWPOINT3D Draw 3D point on the current axis.
%
%   drawPoint3d(X, Y, Z)
%   will draw points defined by coordinates X, Y and Z.
%   X, Y and Z are N*1 array, with N being number of points to be drawn.
%
%   drawPoint3d(COORD) packs coordinates in a single [N*3] array.
%
%   drawPoint3d(..., OPT) will draw each point with given option. OPT is a
%   string compatible with 'plot' model.
%
%   drawPoint3d(AX,...) plots into AX instead of GCA.
%
%   H = drawPoint3d(...) returns a handle H to the line object
%
%   Example
%     % generate points on a 3D circle
%     pts = circleToPolygon([40 30 20], 120);
%     mat = eulerAnglesToRotation3d([30 20 10]);
%     pts3d = transformPoint3d([pts zeros(120,1)],mat);
%     figure; drawPoint3d(pts3d, 'b.');
%     view(3); axis equal;
%
%   See also
%     points3d, clipPoints3d, drawPoint
%

% ---------
% Author : David Legland
% INRA - TPV URPOI - BIA IMASTE
% created the 18/02/2005.
%
%   HISTORY
%   04/01/2007: remove unused variables, and enhance support for plot
%       options
%   12/02/2010 does not clip points anymore
%   12/01/2018 added axes handle input
%

if numel(varargin{1}) == 1 && ishghandle(varargin{1}, 'axes')
    hAx = varargin{1};
    varargin(1)=[];
else
    hAx = gca;
end

if length(varargin) == 1 && size(varargin{1}, 2) == 3
    % points are given as one single array with 3 columns
    px = varargin{1}(:,1);
    py = varargin{1}(:,2);
    pz = varargin{1}(:,3);
    varargin = {};
elseif length(varargin) == 2 && size(varargin{1}, 2) == 3
    % points are given as one single array with 3 columns
    px = varargin{1}(:,1);
    py = varargin{1}(:,2);
    pz = varargin{1}(:,3);
    varargin = varargin(2);
elseif length(varargin) >= 3 && size(varargin{1}, 2) == 3
    % points are given as one single array with 3 columns
    px = varargin{1}(:,1);
    py = varargin{1}(:,2);
    pz = varargin{1}(:,3);
    varargin = varargin(2:end);
elseif length(varargin) == 3 && numel(varargin{1})==numel(varargin{2}) && numel(varargin{1})==numel(varargin{3})
    % points are given as 3 columns with equal lengths
    px = varargin{1};
    py = varargin{2};
    pz = varargin{3};
    varargin = {};
elseif length(varargin) > 3
    % points are given as 3 columns with equal lengths
    px = varargin{1};
    py = varargin{2};
    pz = varargin{3};
    varargin = varargin(4:end);
else
    error('wrong number of arguments in drawPoint3d');
end

% default draw style: no line, marker is 'o'
if length(varargin) ~= 1
    varargin = ['linestyle', 'none', 'marker', 'o', varargin];
end

% plot only points inside the axis.
hh = plot3(hAx, px, py, pz, varargin{:});

if nargout > 0
    h = hh;
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function [ ersp ] = visual_ERSP( data, varargin )
% Description:
%   ERSP (Event-related spectral pergurbation) measures the average dynamic changes
%   in amplitude of the broad band EEG frequency spectrum as a function of time
%   relative to an experimental event
%
% Example code:
%  [dat] = proc_ERSP(data , <OPT>)
%
% Input:
%   data: Data structrue (ex) Epoched data structure
%
% Options:
%   <OPT> :
%      .Channel - Selecting the interested channel in Time-Frequency domain
%                 (e.g. {'Channel', {'C4'}})
%      .Interval - Selecting the interested time intervals
%                 (e.g. {'Interval' , '[-2000 3000]'})
%
% Return:
%    data:  Epoched data structure
%
% See also:
%    opt_cellToStruct , visual_timef
%
% Reference:
%         1. C. Brunner, A. Delorme, and S. Makeig, "EEGLAB?An Open Source
%          Matlab Toolbox for Electrophysiological Research," Biomedical Engineering/Biomedizinische Technik, 2013, pp.1-2.
%         2. EEGLAB tutorial (http://sccn.ucsd.edu/eeglab/)
%          Author: Sigurd Enghoff, Arnaud Delorme & Scott Makeig
%          CNL / Salk Institute 1998- | SCCN/INC, UCSD 2002-
%
%         We used EEGLAB open source toolbox code related in ERSP (timef.m)
%
%
% Ji Hoon, Jeong
% jh_jeong@korea.ac.kr
%
%%
dat = data;
opt = opt_cellToStruct(varargin{:});
chinx = find(strcmp(dat.chan,opt.Channel)==1);
eeg = [];

if ndims(dat.x) == 2
    if ~isfield(opt,'Interval')
        warning('OpenBMI: please input the data intervals');return;
    end
    dat = prep_segmentation(dat, {'interval',opt.Interval});
    dat.x = dat.x(:,:,chinx);
    eeg = reshape(dat.x , [size(dat.x,1)*size(dat.x,2) 1])';
elseif ndims(dat.x) == 3
    dat.x = dat.x(:,:,chinx);
    eeg = reshape(dat.x , [size(dat.x,1)*size(dat.x,2) 1])';
end

figure;
ersp = visual_timef(eeg, size(dat.x,1), [dat.ival(1) dat.ival(end)], dat.fs, [3 0.5], 'plotersp','on','plotitc','off');

author: PatternRecognition

repo: OpenBMI


text:
function determ = jacobi_determinant ( n )

%*****************************************************************************80
%
%% JACOBI_DETERMINANT returns the determinant of the JACOBI matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the order of A.
%
%    Output, real DETERM, the determinant.
%
  if ( mod ( n, 2 ) == 1 )

    determ = 0.0;

  else

    lambda = legendre_zeros ( n );

    determ = prod ( lambda(1:n) );

  end

  return
end


text:
function [ y2, m2, d2, f2 ] = ymdf_next_gregorian ( y1, m1, d1, f1 )

%*****************************************************************************80
%
%% YMDF_NEXT_GREGORIAN returns the Gregorian YMDF date of the next day.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 March 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Output, integer Y2, M2, D2, real F2,
%    tomorrow's YMDF date.
%
  y2 = y1;
  m2 = m1;
  d2 = d1 + 1;
  f2 = f1;

  [ y2, m2, d2 ] = day_carry_gregorian ( y2, m2, d2 );

  return
end


text:
function [E1 E2] = cao_deneme(x,tao,mmax)
%x : time series
%tao : time delay
%mmax : maximum embedding dimension
%reference:Cao, L. (1997), ``Practical method for determining the minimum
%embedding dimension of a scalar time series'', Physcai D, 110, pp. 43-50.
%author:"Merve Kizilkaya"
N=length(x);
for m=1:mmax
    M=N-m*tao;
    Y=psr_deneme(x,m,tao,M);
    for n=1:M
        y0=ones(M,1)*Y(n,:);
        distance=max(abs(Y-y0),[],2);
        [neardis nearpos]=sort(distance);

        newpoint=[Y(n,:) x(n+m*tao)];
        newneig=[Y(nearpos(2),:) x(nearpos(2)+m*tao)];
        R1=max(abs(newpoint-newneig),[],2);

        a(n)=R1/neardis(2);
        d(n)=abs(x(n+m*tao)-x(nearpos(2)+m*tao));
    end
    E(m)=sum(a)/M;
    Ey(m)=sum(d)/M;
end
figure
E1=E(2:end)./E(1:end-1);
E2=Ey(2:end)./Ey(1:end-1);
plot(1:length(E1),E1,'k')
hold on
plot(1:length(E2),E2)
grid on
title('embedding dimension with cao method')
legend('E1','E2',4);

function Y=psr_deneme(x,m,tao,npoint)
%Phase space reconstruction
%x : time series
%m : embedding dimension
%tao : time delay
%npoint : total number of reconstructed vectors
%Y : M x m matrix
% author:"Merve Kizilkaya"
N=length(x);
if nargin == 4
    M=npoint;
else
    M=N-(m-1)*tao;
end

Y=zeros(M,m);

for i=1:m
    Y(:,i)=x((1:M)+(i-1)*tao)';
end

author: Sable

repo: mcbench-benchmarks


text:
U=randn(10,1000);

param.lambda=0.1; % regularization parameter
param.num_threads=-1; % all cores (-1 by default)
param.verbose=true;   % verbosity, false by default
param.pos=false;       % can be used with all the other regularizations
param.intercept=false; % can be used with all the other regularizations

fprintf('First tree example\n');
% Example 1 of tree structure
% tree structured groups:
% g1= {0 1 2 3 4 5 6 7 8 9}
% g2= {2 3 4}
% g3= {5 6 7 8 9}
tree.own_variables=int32([0 2 5]);   % pointer to the first variable of each group
tree.N_own_variables=int32([2 3 5]); % number of "root" variables in each group
                              % (variables that are in a group, but not in its descendants).
                              % for instance root(g1)={0,1}, root(g2)={2 3 4}, root(g3)={5 6 7 8 9}
tree.eta_g=[1 1 1];           % weights for each group, they should be non-zero to use fenchel duality
tree.groups=sparse([0 0 0; ...
                    1 0 0; ...
                    1 0 0]);    % first group should always be the root of the tree
                                % non-zero entriees mean inclusion relation ship, here g2 is a children of g1,
                                % g3 is a children of g1

fprintf('\ntest prox tree-l0\n');
param.regul='tree-l0';
alpha=mexProximalTree(U,tree,param);

fprintf('\ntest prox tree-l2\n');
param.regul='tree-l2';
alpha=mexProximalTree(U,tree,param);

fprintf('\ntest prox tree-linf\n');
param.regul='tree-linf';
alpha=mexProximalTree(U,tree,param);

fprintf('Second tree example\n');
% Example 2 of tree structure
% tree structured groups:
% g1= {0 1 2 3 4 5 6 7 8 9}    root(g1) = { };
% g2= {0 1 2 3 4 5}            root(g2) = {0 1 2};
% g3= {3 4}                    root(g3) = {3 4};
% g4= {5}                      root(g4) = {5};
% g5= {6 7 8 9}                root(g5) = { };
% g6= {6 7}                    root(g6) = {6 7};
% g7= {8 9}                    root(g7) = {8};
% g8 = {9}                     root(g8) = {9};
tree.own_variables=  int32([0 0 3 5 6 6 8 9]);   % pointer to the first variable of each group
tree.N_own_variables=int32([0 3 2 1 0 2 1 1]); % number of "root" variables in each group
tree.eta_g=[1 1 1 2 2 2 2.5 2.5];
tree.groups=sparse([0 0 0 0 0 0 0 0; ...
                    1 0 0 0 0 0 0 0; ...
                    0 1 0 0 0 0 0 0; ...
                    0 1 0 0 0 0 0 0; ...
                    1 0 0 0 0 0 0 0; ...
                    0 0 0 0 1 0 0 0; ...
                    0 0 0 0 1 0 0 0; ...
                    0 0 0 0 0 0 1 0]);  % first group should always be the root of the tree

fprintf('\ntest prox tree-l0\n');
param.regul='tree-l0';
alpha=mexProximalTree(U,tree,param);

fprintf('\ntest prox tree-l2\n');
param.regul='tree-l2';
alpha=mexProximalTree(U,tree,param);

fprintf('\ntest prox tree-linf\n');
param.regul='tree-linf';
alpha=mexProximalTree(U,tree,param);

% mexProximalTree also works with non-tree-structured regularization functions
fprintf('\nprox l1, intercept, positivity constraint\n');
param.regul='l1';
param.pos=true;       % can be used with all the other regularizations
param.intercept=true; % can be used with all the other regularizations
alpha=mexProximalTree([U; ones(1,size(U,2))],tree,param);

% Example of multi-task tree
fprintf('\nprox multi-task tree\n');
param.pos=false;
param.intercept=false;
param.lambda2=param.lambda;
param.regul='multi-task-tree';  % with linf
alpha=mexProximalTree(U,tree,param);

tree.own_variables=int32([0 1 2 3 4 5 6]);   % pointer to the first variable of each group
tree.N_own_variables=int32([1 1 1 1 1 1]); % number of "root" variables in each group
                              % (variables that are in a group, but not in its descendants).
                              % for instance root(g1)={0,1}, root(g2)={2 3 4}, root(g3)={5 6 7 8 9}
tree.eta_g=[1 1 1 1 1 1];           % weights for each group, they should be non-zero to use fenchel duality
tree.groups=sparse([0 0 0; ...
                    1 0 0; ...
                    1 0 0]);    % first group should always be the root of the tree
                                % non-zero entriees mean inclusion relation ship, here g2 is a children of g1,
                                % g3 is a children of g1

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function sp = conventional_bf_1d(R, n, design, wavelength, grid_size, varargin)
%CONVENTIONAL_BF_1D 1D conventional (Bartlett) beamforming.
%Syntax:
%   sp = CONVENTIONAL_BF_1D(R, n, design, wavelength, grid_size, ...);
%   sp = CONVENTIONAL_BF_1D(R, n, f_steering, [], grid_size, ...);
%Inputs:
%   R - Sample covariance matrix.
%   n - Number of sources.
%   design - Array design. Can also be a function handle that generates
%            a steering matrix. This function must take two arguments,
%            wavelength and the doa vector.
%   wavelength - Wavelength.
%   grid_size - Number of grid points used.
%   ... - Options:
%           'Unit' - Can be 'radian', 'degree', or 'sin'. Default value is
%                   'radian'.
%           'RefineEstimates' - If set to true, will refine the estimated
%                               direction of arrivals around the grid.
%Output:
%   sp - Spectrum structure with the following fields:
%           x - An 1 x grid_size vector.
%           y - An 1 x grid_size vector. Calling `plot(x, y)` will plot the
%               spectrum.
%           x_est - An 1 x n vector storing the estimated DOAs. May not
%                   fall on the grid if 'RefineEstimates' is set to true.
%           x_unit - The same as the unit specified by 'Unit'.
%           resolved - True if the number of peaks in the spectrum is
%                      greater or equal to the number of sources.
%           discrete - Constant value false.
%Reference:
%   [1] H. L. Van Trees, Optimum array processing. New York: Wiley, 2002.
unit = 'radian';
refine_estimates = false;
for ii = 1:2:nargin-5
    option_name = varargin{ii};
    option_value = varargin{ii+1};
    switch lower(option_name)
        case 'unit'
            unit = lower(option_value);
        case 'refineestimates'
            refine_estimates = true;
        otherwise
            error('Unknown option ''%s''.', option_name);
    end
end
% discretize and create the corresponding steering matrix
[doa_grid_rad, doa_grid_display, ~] = default_doa_grid(grid_size, unit, 1);
% compute spectrum
sp_intl = compute_spectrum(R, design, wavelength, doa_grid_rad);
[x_est, x_est_idx, resolved] = find_doa_from_spectrum_1d(doa_grid_display, sp_intl, n);
% refine
if resolved && refine_estimates
    switch unit
        case 'radian'
            f_obj = @(x) compute_spectrum(R, design, wavelength, x);
        case 'degree'
            f_obj = @(x) compute_spectrum(R, design, wavelength, deg2rad(x));
        case 'sin'
            f_obj = @(x) compute_spectrum(R, design, wavelength, asin(x));
        otherwise
            error('Invalid unit ''%s''.', unit);
    end
    x_est = refine_grid_estimates(f_obj, doa_grid_rad, x_est_idx);
end
% return
sp = struct();
sp.x = doa_grid_display;
sp.x_est = x_est;
sp.x_unit = unit;
sp.y = sp_intl;
sp.resolved = resolved;
sp.discrete = false;
end

function v = compute_spectrum(R, design, wavelength, theta)
if ishandle(design)
    A = design(wavelength, theta);
else
    A = steering_matrix(design, wavelength, theta);
end
v = real(sum(conj(A) .* (R * A), 1));
end

author: morriswmz

repo: doa-tools

save_path: github-repos/MATLAB/morriswmz-doa-tools


text:
function editCobraToolboxPath(basePath,folderPattern,method)
%recursively finds all the folders within a given basePath that contains a
%given folderPattern and either adds or removes that folder from the matlab
%path

if ~exist('basePath','var') || isempty(basePath)
    aPath = which('initVonBertalanffy');
    basePath = strrep(aPath,'vonBertalanffy/initVonBertalanffy.m','');
end
if ~exist('folderPattern','var') || isempty(folderPattern)
    folderPattern=[filesep 'new'];
end
if ~exist('method','var') || isempty(method)
    method = 'remove';
end

folderlist = dir(fullfile(basePath, ['**' filesep '*.*']));  %get list of files and folders in any subfolder
folderlist = folderlist([folderlist.isdir]);  %remove files from list

folderlistOnly=cell(length(folderlist),1);
for i=1:length(folderlist)
    folderlistOnly{i,1} = folderlist(i).folder;
end
folderlistOnly = unique(folderlistOnly);

for i=1:length(folderlistOnly)
    if contains(folderlistOnly{i,1},folderPattern)
        switch method
            case 'add'
                disp(['adding: ' folderlistOnly{i,1}])
                addpath(folderlistOnly{i,1})
            case 'remove'
                disp(['removing: ' folderlistOnly{i,1}])
                rmpath(folderlistOnly{i,1})
            otherwise
                error([method ' is not a recognised method'])
        end
    end
end

author: opencobra

repo: cobratoolbox


text:
function mv = multiVoxelUI(view, roi , scans, dt);
%   mv = multiVoxelUI(view, [roi , scans, dt])
%
% multiVoxelUI: user interface for displaying time course / amplitude
% data across roi  within an ROI.
%
%
%
% ARGS:
% input:
%       view: mrVista view struct
%
%       roi : coords (by the current view's conventions) of the roi
%       from which to plot data. May also be a string providing the name of
%       the ROI to plot. Default is view's currently-selected ROI.
%
%       scans: scans in the current view/data dt to plot. Default is
%       currently-selected scan. For multiple scans,  will concatenate
%       together after detrending (using the detrend options saved in the
%       dataTYPES struct in the mrSESSION file).
%
%       dt: data type to use. Defaults to currently selected data
%       type.
%
%output:
%       h: handle to time course interface.
%
%       data: array containing time course data from selected voxels /scans.
%
% More notes:
% This is currently implemented as a small set of linked analysis tools,
% all starting with the prefix 'tc_'. This particular piece of code
% decides,  based on the number of input args,  whether to open a new
% mv user interface or just update an existing one,  and then refreshes the
% current mv view according to the time course data and user prefs.
%
% Here's how it decides to open a new UI or not: if a view is passed as the
% first argument,  it opens a new one. If there are no input args,  it
% assumes the current figure has the interface and refreshes it. (This is
% called by all the callbacks on the interface).
%
% The time course data,  condition information,  user prefs,  and handles to
% important controls are all kept in a mv struct stored in the user
% interface's 'UserData' field. I currently treat each UI as a separate
% figure; though in time I may find it useful to group several interfaces
% on multiple axes in the same figure.

% 1/3/04 ras: wrote it.
% 1/04: added signal to noise ratios,  ability to switch off different cond
% displays,  shift onsets,  change colors
% 2/04: added ability to read multiple scans / concat parfiles if it uses
% them
% 2/23/04: broke up into several functions,  under mrLoadRet sub-section
% Plots/TimeCourseUI. Also made an offshoot of this,  fancyplot2,  which
% doesn't use the mrLoadRet view info,  for kalanit. Renamed multiVoxelUI
% (previously ras_tc -- saved that as an older,  stable version).
mrGlobals

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1ST DECISION: extract new time course data (w/ corresponding
% axes figure) or just re-render existing data? This is decided by
% the # args: if 1,  update data; otherwise,  just render.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin >= 1
    openFigFlag = 1;
else
    openFigFlag = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Open the figure if needed %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if openFigFlag
    if notDefined('dt')
        dt = viewGet(view, 'curDataType');
    end

    if notDefined('scans')
        scans = er_getScanGroup(view); %viewGet(view, 'curScan');
    end

    if notDefined('roi')
        r = viewGet(view, 'curRoi');
        if r==0
            myErrorDlg('Sorry, you need to load an ROI first.')
        end
        roi = view.ROIs(r);
    end

    % want roi to be a struct: parse
    % how it's passed in and make a struct
    roi = tc_roiStruct(view, roi);

    % init a voxel data struct from params
    % (this includes loading the data)
    mv = mv_init(view, roi, scans, dt);

    % put up the interface window
    mv = mv_openFig(mv);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RENDER TIME COURSE (taking into account prefs as specified in the
% ui)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get the data struct -- in the current axes' UserData
mv = get(gcf, 'UserData');

% clear previous axes in the figure
old = findobj('Type', 'axes', 'Parent', gcf);
old = [old; findobj('Type', 'uicontrol', 'Parent', gcf)];
delete(old);

% legend (kind of a hack -- legends are tricky)
if mv.params.legend==1
    sel = find(tc_selectedConds(mv));
    co = [];
    for i = 1:length(sel)
        % The legend text should be the specified condition names
        leg{i} = mv.condNames{sel(i)};
        % remove underscores -- the teX interpreter is weird about em
        leg{i}(findstr(leg{i}, '_')) = ' ';
        % also get appropriate colors
        co = [co; mv.condColors{sel(i)}];
    end
    set(gca, 'NextPlot', 'ReplaceChildren');
    if length(sel) > 0
        % Plot some random data,  just to get the legend put up
        set(gca, 'ColorOrder', co);
        plot(rand(2, length(sel)), 'LineWidth', 4);
        legend(gca, leg, -1);
    end
end

% select the current plot type
eval( get(mv.ui.plotHandles(mv.ui.plotType), 'Callback') )

% stash the updated data struct
set(gcf, 'UserData', mv);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function value = pei_condition ( alpha, n )

%*****************************************************************************80
%
%% PEI_CONDITION returns the L1 condition of the PEI matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real ALPHA, the scalar that defines the Pei matrix.  A
%    typical value of ALPHA is 1.0.
%
%    Input, integer N, the order of A.
%
%    Output, real VALUE, the L1 condition.
%
  a_norm = abs ( alpha + 1.0 ) + n - 1;
  b_norm = ( abs ( alpha + n - 1.0 ) + n - 1.0 ) ...
    / abs ( alpha * ( alpha + n ) );
  value = a_norm * b_norm;

  return
end


text:
function ind = gaussian_mutate( ind, prob, domain)
%GAUSSIAN_MUTATE Summary of this function goes here
%   Detailed explanation goes here

if isstruct(ind)
    x = ind.parameter;
else
    x  = ind;
end

   parDim = length(x);
   lowend  = domain(:,1);
   highend =domain(:,2);
   sigma = (highend-lowend)./20;

   newparam = min(max(normrnd(x, sigma), lowend), highend);
   C = rand(parDim, 1)<prob;
   x(C) = newparam(C);

if isstruct(ind)
    ind.parameter = x;
else
    ind = x;
end

author: sfvsfv

repo: Mathematical-modeling


text:
function cvx_optval = norm_largest( x, k )

%NORM_LARGEST Sum of the k largest magnitudes of a vector.
%   NORM_LARGEST( X, k ) computes the 'largest-k' norm; that is, it computes
%   the sum of the magnitudes of the k largest elements in X. X must
%   be a vector, and k must be a real scalar.
%
%   Disciplined convex programming information:
%       NORM_LARGEST is convex and nonmonotonic of X, though it is monotonic for
%       positive values of X. So when used in CVX expressions, X must be affine,
%       monomial, or posynomial. k must be a real scalar constant.

if ~any( size( x ) ~= 1 ),
    cvx_throw( 'The first argument must be a vector.' );
elseif ~isnumeric( k ) || ~isreal( k ) || length( k ) ~= 1,
    cvx_throw( 'Third argument must be a scalar.' );
else
    cvx_optval = sum_largest( abs( x ), k );
end

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

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


text:
function writeProcessedImg(rob,sen,frm,Raw,ExpOpt)

% WRITEPROCESSEDIMG  Write processed image data to file.
%   WRITEPROCESSEDIMG(R,S,F,Raw) writes the processed image in Raw for
%   sensor S in robot R corresponding to frame F. R and S are Rob and Sen
%   structures in SLAMTB. F is an integer corresponding to the frame
%   number.
%
%   The processed image contains simply the coordinates of all the features
%   detected, with an identifier. See below for the file format.
%
%   The file written has the following name format
%       Directory: ./data/ <ExpOpt.site> /
%       Name: ExpOpt.procImgName , e.g. 'procImg-r%02d-s%02d-i%06d.txt'
%   where the first two indices are the robot and sensor numbers (2 digits
%   each) and the third is the frame number (6 digits). The file contains
%   one line per feature, containing identifier and UV coordinates
%   separated by tabs:
%           id1  U1  V1
%           id2  U2  V2
%           ...
%           idn  Un  Vn
%
%   See also READPROCESSEDIMG, WRITECONTROLSIGNAL, READCONTROLSIGNAL.

%   Copyright 2013 Joan Sola

dir = [ExpOpt.root  'data/' ExpOpt.site '/'];
if ~isdir(dir)
    mkdir(dir);
end
filename = sprintf(ExpOpt.procImgName,rob,sen,frm);
filepath = [dir filename];

fid = fopen(filepath,'w');
ids = Raw.data.points.app;
coord = Raw.data.points.coord;
for i = 1:numel(ids)
    count = fprintf(fid,'%02d\t%f\t%f\n',ids(i),coord(1,i),coord(2,i));
end
fclose(fid);

% ========== End of function - Start GPL license ==========

%   # START GPL LICENSE

%---------------------------------------------------------------------
%
%   This file is part of SLAMTB, a SLAM toolbox for Matlab.
%
%   SLAMTB is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   SLAMTB is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with SLAMTB.  If not, see <http://www.gnu.org/licenses/>.
%
%---------------------------------------------------------------------

%   SLAMTB is Copyright:
%   Copyright (c) 2008-2010, Joan Sola @ LAAS-CNRS,
%   Copyright (c) 2010-2013, Joan Sola,
%   Copyright (c) 2014-2015, Joan Sola @ IRI-UPC-CSIC,
%   SLAMTB is Copyright 2009
%   by Joan Sola, Teresa Vidal-Calleja, David Marquez and Jean Marie Codol
%   @ LAAS-CNRS.
%   See on top of this file for its particular copyright.

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
function value = r8mat_norm_fro_affine ( m, n, a1, a2 )

%*****************************************************************************80
%
%% R8MAT_NORM_FRO_AFFINE returns the Frobenius norm of an R8MAT difference.
%
%  Discussion:
%
%    The Frobenius norm is defined as
%
%      value = sqrt ( sum ( 1 <= I <= M ) sum ( 1 <= j <= N ) A(I,J)^2 )
%
%    The matrix Frobenius norm is not derived from a vector norm, but
%    is compatible with the vector L2 norm, so that:
%
%      vec_norm_l2 ( A * x ) <= mat_norm_fro ( A ) * vec_norm_l2 ( x ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the number of rows.
%
%    Input, integer N, the number of columns.
%
%    Input, real A1(M,N), A2(M,N), the matrices for whose difference
%    the Frobenius norm is desired.
%
%    Output, real VALUE, the Frobenius norm of A1 - A2.
%
  value = sqrt ( sum ( sum ( ( a1(1:m,1:n) - a2(1:m,1:n) ).^2 ) ) );

  return
end


text:
classdef tNewInterface < matlab.unittest.TestCase

    properties
        testLoc
        RelTol = 1e-5
        AbsTol = 1e-6
        ToAvoid = [...
            "checkRun", "destinationfile", "estimationinfo", ...
            "BEARpath","datapath","filespath","pref", ...
            "replicationpath","settingspath","sourcefile","settingsm", ...
            "const", "VARtype", "OLS_Bhat", ...
            "IRFt", "IRF", "HD", "Feval", "Fendsmpl", "FEVD", ...
            "F", "CFt", "CF", "ii"]
    end

    methods (TestClassSetup)

        function setup(tc)

            tc.testLoc = fileparts(mfilename('fullpath'));
            % Need to run single threaded to get all the rng defaults
            % correct.
            ps = parallel.Settings;
            if ps.Pool.AutoCreate
                ps.Pool.AutoCreate=false;
                tc.addTeardown(@() set(ps.Pool,'AutoCreate',true))
            end

        end

    end

    methods (TestMethodSetup)

        function prepareTest(tc)
            rng('default');
            s = rng;
            addTeardown(tc, @() rng(s))
        end

    end

    methods (Test, TestTags = {'Git'})

        function tOLSVAR(tc)

            import matlab.unittest.fixtures.TemporaryFolderFixture
            tempFixture = tc.applyFixture(TemporaryFolderFixture);

            resultsFile = "newTest";

            opts= BEARsettings('OLS', 'ExcelFile', fullfile(fullfile(bearroot(),'default_bear_data.xlsx')));
            opts.results_path = tempFixture.Folder;
            opts.results_sub = 'newTest';
            opts.plot = false;

            BEARmain(opts);

            file = exist(fullfile(tempFixture.Folder, resultsFile + ".mat"), 'file');
            tc.verifyEqual(file, 2)
        end

        function tOLSVAR_IRFt2(tc)

            import matlab.unittest.fixtures.TemporaryFolderFixture
            tempFixture = tc.applyFixture(TemporaryFolderFixture);

            resultsFile = "newTest";

            opts= BEARsettings('OLS', 'ExcelFile', fullfile(fullfile(bearroot(),'default_bear_data.xlsx')));
            opts.results_path = tempFixture.Folder;
            opts.results_sub = 'newTest';
            opts.plot = false;
            opts.IRFt = 2;

            BEARmain(opts);

            file = exist(fullfile(tempFixture.Folder, resultsFile + ".mat"), 'file');
            tc.verifyEqual(file, 2)
        end

        function tBVAR_IRFt2(tc)

            import matlab.unittest.fixtures.TemporaryFolderFixture
            tempFixture = tc.applyFixture(TemporaryFolderFixture);

            resultsFile = "newTest";

            opts = BEARsettings('BVAR', 'ExcelFile', fullfile(fullfile(bearroot(),'default_bear_data.xlsx')), ...
                'prior', 'minnesota_univariate', 'IRFt', 4);
            opts.prior=11;
            opts.IRFt=2;
            opts.results_path = tempFixture.Folder;
            opts.results_sub = 'newTest';
            opts.plot = false;
            BEARmain(opts);

            file = exist(fullfile(tempFixture.Folder, resultsFile + ".mat"), 'file');
            tc.verifyEqual(file, 2)
        end

    end

end

author: european-central-bank

repo: BEAR-toolbox


text:
% This script defines amino acids that can be stored by one or more organs
storageAA = {'sink_his_L(c)';
    'sink_ile_L(c)';
    'sink_leu_L(c)';
    'sink_lys_L(c)';
    'sink_met_L(c)';
    'sink_phe_L(c)';
    'sink_thr_L(c)';
    'sink_trp_L(c)';
    'sink_val_L(c)';
    % non-essential AAs
    'sink_ala_L(c)'
'sink_arg_L(c)'
'sink_asn_L(c)'
'sink_asp_L(c)'
'sink_cys_L(c)'
'sink_gln_L(c)'
'sink_glu_L(c)'
'sink_pro_L(c)'
'sink_ser_L(c)'
'sink_tyr_L(c)'
'sink_gly(c)'
    }


text:
%DEMO_SMOOTHPOLYGON  One-line description here, please.
%
%   output = demo_smoothPolygon(input)
%
%   Example
%   demo_smoothPolygon
%
%   See also
%

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% INRAE - BIA Research Unit - BIBS Platform (Nantes)
% Created: 2022-02-18,    using Matlab 9.10.0.1739362 (R2021a) Update 5
% Copyright 2022 INRAE.

% read data
poly = load('leaf_poly.txt');

% display polygon
figure; axis equal; axis([0 600 0 450]); hold on;
drawPolygon(poly, 'lineWidth', 3, 'color', 'k');

% smooth strongly
polyM11 = smoothPolygon(poly, 11);
drawPolygon(polyM11, 'color', 'g', 'linewidth', 2);
polyM31 = smoothPolygon(poly, 31);
drawPolygon(polyM31, 'color', 'r', 'linewidth', 2);
axis([0 200 150 300]);

% print(gcf, 'leafPoly_smoothM31.png', '-dpng');

legend({'Original polygon', 'Smoothed polygon (M=11)', 'Smoothed polygon (M=31)'}, 'Location', 'NorthEast');
print(gcf, 'leafPoly_smoothM31_annot.png', '-dpng');

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
%Function of Circular Convolution for the Overlap Save Method.
function y=mycirconv(x,h)
lx=length(x);
lh=length(h);
l=max(lx,lh);
X=[x zeros(1,l-lx)];
H=zeros(l);
H(1:lh)=h;
for j=1:l-1
for i=1:l-1
H(i+1,j+1)=H(i,j);
end
H(1,j+1)=H(l,j);
end
y=H*X';


text:
% Copyright (C) 1993-2019 Peter I. Corke
%
% This file is part of The Spatial Math Toolbox for MATLAB (SMTB).
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, including without limitation the rights
% to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
% of the Software, and to permit persons to whom the Software is furnished to do
% so, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
% FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
% COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
% IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
% CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
%
% https://github.com/petercorke/spatial-math

function smtbcheck
    % display current versions of MATLAB
    year = version('-release');
    fprintf('You are using:\n - MATLAB release %s\n', year);

    % check how old it is
    today = datevec(now);
    age =  today(1) - str2num(year(1:4));
    if age >= 2
        fprintf('     ** this is at least %d years old, you may have issues\n', age);
    end

    %% display versions of toolboxes (use unique RTB and MVTB functions)

    %SMTB
    smtb = checktoolbox('RTBPose', 'SMTB');

    %% check for shadowed files

    shadows = {};

    shadows = checkshadow(shadows, 'rotx');
    shadows = checkshadow(shadows, 'roty');
    shadows = checkshadow(shadows, 'rotz');
    shadows = checkshadow(shadows, 'angdiff');

    if ~isempty(shadows)
        fprintf('\n*** Some Toolbox files are "shadowed" and will cause problems with the use of this toolbox ***\n');
        for s = shadows'
            fprintf(' - toolbox function "%s" is shadowed by the MathWorks file %s\n', s{1}, s{2} );
        end
        fprintf('Use the pathtool function to:\n 1. move the Toolbox containing the shadowed function to the top of the path\n 2. Delete the problematic MathWorks toolbox (do you really need it?)\n 3. Move the problematic MathWorks toolbox to the end of the path\n')
    end
end

function present = checktoolbox(funcname, tbname)

    present = true;

    p = which(funcname);
    if isempty(p)
        % toolbox not present
        present = false;
        return
    end

    % display versions of toolbox
    fprintf('%s (%s)\n', short2long(tbname), p);
    v = ver( p );

    where = [];
    if strcmp(tbname, 'SMTB')
        % special case for SMTB, is it standalone or part of RTB or MVTB?
        if contains(p, 'vision', 'IgnoreCase', true)
            where = 'included with MVTB';
        elseif contains(p, 'robot', 'IgnoreCase', true)
            where = 'included with RTB';
        end
    end
    if isempty(where)
        if findstr(p, 'Add-Ons')
            where = 'Add-Ons (installed from mltbx file or Add-On Explorer)';
        else
            if exist( fullfile(p, '.git'), 'dir' )
                where = 'local (git clone)';
            else
                where = 'local (zip install)';
            end
        end
    end

    if ~isempty(v)
        fprintf(' - %s%s, %s\n', tbname, v.Version, v.Date);
    end
    fprintf(' - %s\n', where);
end

function out = checkshadow(shadows, funcname)
    % add to the list if function lives below MATLAB root
    funcpath = which(funcname); % find first instance in path
    out = shadows;
    if startsWith(funcpath, matlabroot) || startsWith(funcpath, 'built-in')
        out = [out; {funcname, which(funcname)}];
    end
end

function fullname = short2long(shortname)
    switch shortname
        case 'RTB'
            fullname = 'Robotics Toolbox for MATLAB';
        case 'MVTB'
            fullname = 'Machine Vision Toolbox for MATLAB';
        case 'SMTB'
            fullname = 'Spatial Math Toolbox for MATLAB';
    end
end

function p = getpath(funcname)
    funcpath = which(funcname);
    k = strfind( funcpath, [filesep funcname]);
    p = funcpath(1:k-1);
end

author: petercorke

repo: spatialmath-matlab


text:
function x=fibosearch(fhandle,a,b,npoints)
% fibonacci search for minimum of unknown unimodal function in one variable
%     x = fibosearch(fhandle,a,b,npoints)
% a,b define the search interval with resolution 1/npoints
%create fibonacci sequence of length nfibo
nfibo=22;
fibo=[1,1,zeros(1,nfibo-2)];
for k=1:nfibo-2
    fibo(k+2)=fibo(k+1)+fibo(k);
end
%find number of required iterations
fiboindex=3;
while fibo(fiboindex)<npoints
        fiboindex=fiboindex+1;
end
for k=1:fiboindex-2
    if k==1
        x1 = a+fibo(fiboindex-k-1)/fibo(fiboindex-k+1)*(b-a);
        x2 = b-fibo(fiboindex-k-1)/fibo(fiboindex-k+1)*(b-a);
        fx1 = fhandle(x1);
        fx2 = fhandle(x2);
    end
    if fx1<fx2
        b=x2;
        x2=x1; fx2=fx1;
        x1=a+fibo(fiboindex-k-1)/fibo(fiboindex-k+1)*(b-a);
        fx1=fhandle(x1);
    else
        a=x1;
        x1=x2; fx1=fx2;
        x2=b-fibo(fiboindex-k-1)/fibo(fiboindex-k+1)*(b-a);
        fx2=fhandle(x2);
    end
end
if fx1<fx2
    x=x1;
else
    x=x2;
end
% disp(fiboindex-2)

author: panpanfei


text:
function dx = f_Wrapper(u,x,p,parIdx)
% dx = f(u,x,p)
% parIdx: index of the core (for reentrant purpose)

    if coder.target('MATLAB')
        % Specify your own f(u,x,p) function for normal excution

    else
        % Specify your own f(u,x,p) function for code generation

    end

end

author: deng-haoyang

repo: ParNMPC


text:
% SR_QUANT_EVAL
%
% Script for computing averaged PSNR, SSIM, and IVC.
% The script reproduces the quantitative results on our project page on
% Github: https://github.com/jbhuang0604/SelfExSR
%
% Jia-Bin Huang
% Electrical and Computer Engineering
% University of Illinois, Urbana-Champaign
% www.jiabinhuang.com

clc;
clear;
close all;

% Add pathes for running SSIM and IVC
addpath(genpath('quant_eval'));

% Dataset list
numDataset = 4;
dataset = cell(numDataset,1);
dataset{1}.name = 'Urban100';
dataset{2}.name = 'BSD100';
dataset{3}.name = 'Set5';
dataset{4}.name = 'Set14';

% Super-resolution factor
dataset{1}.SRF  = [2, 4];
dataset{2}.SRF  = [2, 3, 4];
dataset{3}.SRF  = [2, 3, 4];
dataset{4}.SRF  = [2, 3, 4];

% Number of images in the dataset
dataset{1}.numImg = 100;
dataset{2}.numImg = 100;
dataset{3}.numImg = 5;
dataset{4}.numImg = 14;

% Method list
methodList = {'Bicubic', 'ScSR', 'Kim', 'Abhishek', 'Glasner', 'SRCNN', 'A+', 'SelfExSR'};
numMethod  = length(methodList);
indValidMethod = [1, 1, 1, 1, 1, 1, 1, 1]; % Test only

% Initialize result path

resPath = fullfile('quant_eval', 'result');
if(~exist(resPath, 'dir'))
    mkdir(resPath);
end
% =========================================================================
% Start quantitative evaluation
% =========================================================================
for indDataset = 1 : numDataset
    datasetName   = dataset{indDataset}.name;
    numImgDataset = dataset{indDataset}.numImg;
    numSRF        = size(dataset{indDataset}.SRF,2);

    % Run each super-resolution factor
    for indSRF = 1: numSRF
        SRF = dataset{indDataset}.SRF(indSRF);
        imgPath = fullfile('data', datasetName, ['image_SRF_',num2str(SRF)]);

        resName = ['result_', datasetName, '_SRF_', num2str(SRF),'.mat'];
        if(~exist(fullfile(resPath, resName), 'file'))
            % Initialize result table
            SSIM_table = zeros(numImgDataset, numMethod);
            PSNR_table = zeros(numImgDataset, numMethod);
            IFC_table  = zeros(numImgDataset, numMethod);
            reverseStr = '';
            for indImg = 1: numImgDataset
                % Load groundtruth high-resolution image
                imgName = ['img_',num2str(indImg, '%03d'), '_SRF_', num2str(SRF), '_HR.png'];
                imgGT = imread(fullfile(imgPath, imgName));

                % Load super-resolution results
                for indMethod= 1:numMethod
                    if(indValidMethod(indMethod))
                        imgName = ['img_',num2str(indImg, '%03d'), '_SRF_', num2str(SRF), '_', methodList{indMethod}, '.png'];
                        img = imread(fullfile(imgPath, imgName), 'png');

                        % Compute image quality
                        [psnr, ssim, ifc] = compute_difference(img, imgGT, SRF);

                        PSNR_table(indImg,indMethod) = psnr;
                        SSIM_table(indImg,indMethod) = ssim;
                        IFC_table(indImg, indMethod)  = ifc;
                    end
                    % Display progress
                    percentDone = 100*((indImg-1)*numMethod + indMethod)/(numMethod*numImgDataset);
                    msg = sprintf('Evaluating dataset %s on SRF %d, progress: %3.2f', datasetName, SRF, percentDone);
                    fprintf([reverseStr, msg]);
                    reverseStr = repmat(sprintf('\b'), 1, length(msg));
                end
            end

            % Save results
            resName = ['result_', datasetName, '_SRF_', num2str(SRF),'.mat'];
            save(fullfile(resPath, resName), 'PSNR_table', 'SSIM_table', 'IFC_table');
        else
            load(fullfile(resPath, resName));
        end

        % Display results
        avgPSNR = mean(PSNR_table,1);
        avgSSIM = mean(SSIM_table,1);
        avgIFC  = mean(IFC_table, 1);
        fprintf('\n\n=== Quantitative results for dataset %s on SRF %d === \n\n', datasetName, SRF);
        fprintf('Peak signal-to-noise ratio (PSNR) \n')
        fprintf('     %8s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s\t \n', methodList{1}, methodList{2}, ...
            methodList{3}, methodList{4}, methodList{5}, methodList{6}, methodList{7}, methodList{8});
        fprintf('PSNR|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t| \n', ...
            avgPSNR(1), avgPSNR(2), avgPSNR(3), avgPSNR(4), avgPSNR(5), avgPSNR(6), avgPSNR(7), avgPSNR(8));
        fprintf('Structural similarity index (SSIM) \n')
        fprintf('     %8s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s\t \n', methodList{1}, methodList{2}, ...
            methodList{3}, methodList{4}, methodList{5}, methodList{6}, methodList{7}, methodList{8});
        fprintf('SSIM|%8.04f\t|%8.04f\t|%8.04f\t|%8.04f\t|%8.04f\t|%8.04f\t|%8.04f\t|%8.04f\t| \n', ...
            avgSSIM(1), avgSSIM(2), avgSSIM(3), avgSSIM(4), avgSSIM(5), avgSSIM(6), avgSSIM(7), avgSSIM(8));
        fprintf('Information fidelity criterion (IFC) \n')
        fprintf('    %8s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s\t%8s\t \n', methodList{1}, methodList{2}, ...
            methodList{3}, methodList{4}, methodList{5}, methodList{6}, methodList{7}, methodList{8});
        fprintf('IFC|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t|%8.02f\t| \n\n', ...
            avgIFC(1), avgIFC(2), avgIFC(3), avgIFC(4), avgIFC(5), avgIFC(6), avgIFC(7), avgIFC(8));
        fprintf('\n');
    end
end

repo: SelfExSR


text:
function value = i4_sign ( i )

%*****************************************************************************80
%
%% I4_SIGN returns the sign of an I4.
%
%  Discussion:
%
%    The value is +1 if the number is positive or zero, and it is -1 otherwise.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the number whose sign is desired.
%
%    Output, integer VALUE, the sign of I.
%
  if ( i < 0 )
    value = -1;
  else
    value = +1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function ellipsoid_grid_test ( )

%*****************************************************************************80
%
%% ELLIPSOID_GRID_TEST tests the ELLIPSOID_GRID library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 September 2010
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPSOID_GRID_TEST:\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the ELLIPSOID_GRID library.\n' );

  ellipsoid_grid_test01 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'ELLIPSOID_GRID_TEST:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

  return
end


text:
function hout=suptitle(str)
%suptitle puts a title above all subplots.
%
% suptitle('text') adds text to the top of the figure above all subplots.

% Amount of the figure window devoted to subplots
plotregion = .92;

% Y position of title in normalized coordinates
titleypos  = .95;

% Fontsize for supertitle
fs = get(gcf,'defaultaxesfontsize')+4;

% Fudge factor to adjust y spacing between subplots
fudge=1;

haold = gca;
figunits = get(gcf,'units');

% Get the (approximate) difference between full height (plot + title
% + xlabel) and bounding rectangle.

if ~strcmp(figunits,'pixels')
    set(gcf,'units','pixels');
    pos = get(gcf,'position');
    set(gcf,'units',figunits);
else
    pos = get(gcf,'position');
end
ff = (fs-4)*1.27*5/pos(4)*fudge;

% The 5 here reflects about 3 characters of height below
% an axis and 2 above. 1.27 is pixels per point.

% Determine the bounding rectangle for all the plots

h = findobj(gcf,'Type','axes');

oldUnits = get(h, {'Units'});
if ~all(strcmp(oldUnits, 'normalized'))
    % This code is based on normalized units, so we need to temporarily
    % change the axes to normalized units.
    set(h, 'Units', 'normalized');
    cleanup = onCleanup(@()resetUnits(h, oldUnits));
end

max_y=0;
min_y=1;
oldtitle = [];
numAxes = length(h);
thePositions = zeros(numAxes,4);
for i=1:numAxes
    pos=get(h(i),'pos');
    thePositions(i,:) = pos;
    if ~strcmp(get(h(i),'Tag'),'suptitle')
        if pos(2) < min_y
            min_y=pos(2)-ff/5*3;
        end
        if pos(4)+pos(2) > max_y
            max_y=pos(4)+pos(2)+ff/5*2;
        end
    else
        oldtitle = h(i);
    end
end

if max_y > plotregion
    scale = (plotregion-min_y)/(max_y-min_y);
    for i=1:numAxes
        pos = thePositions(i,:);
        pos(2) = (pos(2)-min_y)*scale+min_y;
        pos(4) = pos(4)*scale-(1-scale)*ff/5*3;
        set(h(i),'position',pos);
    end
end

np = get(gcf,'nextplot');
set(gcf,'nextplot','add');
if ~isempty(oldtitle)
    delete(oldtitle);
end
axes('pos',[0 1 1 1],'visible','off','Tag','suptitle');
ht=text(.5,titleypos-1,str);set(ht,'horizontalalignment','center','fontsize',fs);
set(gcf,'nextplot',np);
axes(haold);
if nargout
    hout=ht;
end
end

function resetUnits(h, oldUnits)
    % Reset units on axes object. Note that one of these objects could have
    % been an old supertitle that has since been deleted.
    valid = isgraphics(h);
    set(h(valid), {'Units'}, oldUnits(valid));
end

author: MultiDIC

repo: MultiDIC

save_path: github-repos/MATLAB/MultiDIC-MultiDIC


text:
function V = spm_close_vol(V)
% Close image volume
% See: spm_create_vol and spm_write_plane.
%_______________________________________________________________________
% @(#)spm_close_vol.m	2.4 John Ashburner 02/08/16
for i=1:prod(size(V)),
	if isfield(V,'private') & isfield(V(i).private,'fid') & ~isempty(V(i).private.fid),
		if ~isempty(fopen(V(i).private.fid)),
			fclose(V(i).private.fid);
		end;
		V(i).private = rmfield(V(i).private,'fid');
	end;
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function varargout = process_granger1( varargin )
% PROCESS_GRANGER1: Compute the Granger causality between one signal and all the others, in one file.

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Francois Tadel, 2012-2020

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'Bivariate Granger causality 1xN';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Connectivity';
    sProcess.Index       = 660;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Connectivity';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'data',     'results',  'matrix'};
    sProcess.OutputTypes = {'timefreq', 'timefreq', 'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;

    % === CONNECT INPUT
    sProcess = process_corr1n('DefineConnectOptions', sProcess, 0);
    % === REMOVE EVOKED REPONSE
    sProcess.options.removeevoked.Comment = 'Remove evoked response from each trial';
    sProcess.options.removeevoked.Type    = 'checkbox';
    sProcess.options.removeevoked.Value   = 0;
    sProcess.options.removeevoked.Group   = 'input';
    % === DIRECTION
    sProcess.options.dirlabel.Comment  = 'Direction of the causality:';
    sProcess.options.dirlabel.Type     = 'label';
    sProcess.options.direction.Comment = {'From the selected node (out)', 'To the selected node (in)', 'Both (generates two files)'};
    sProcess.options.direction.Type    = 'radio';
    sProcess.options.direction.Value   = 3;
    % === GRANGER ORDER
    sProcess.options.grangerorder.Comment = 'Maximum Granger model order (default=10):';
    sProcess.options.grangerorder.Type    = 'value';
    sProcess.options.grangerorder.Value   = {10, '', 0};
    % === OUTPUT MODE
    sProcess.options.outputmode.Comment = {'Save individual results (one file per input file)', 'Concatenate input files before processing (one file)', 'Save average connectivity matrix (one file)'};
    sProcess.options.outputmode.Type    = 'radio';
    sProcess.options.outputmode.Value   = 1;
    sProcess.options.outputmode.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFiles = Run(sProcess, sInputA) %#ok<DEFNU>
    % Input options
    OPTIONS = process_corr1n('GetConnectOptions', sProcess, sInputA);
    if isempty(OPTIONS)
        OutputFiles = {};
        return
    end

    % Metric options
    OPTIONS.Method = 'granger';
    OPTIONS.RemoveEvoked = sProcess.options.removeevoked.Value;
    OPTIONS.GrangerOrder = sProcess.options.grangerorder.Value{1};
    OPTIONS.pThresh      = 0.05;
    OPTIONS.Freqs        = 0;

    % Computation depends on the direction
    OutputFiles = {};
    if ismember(sProcess.options.direction.Value, [1 3])
        OPTIONS.GrangerDir = 'out';
        OutputFiles = cat(2, OutputFiles, bst_connectivity({sInputA.FileName}, {sInputA.FileName}, OPTIONS));
    end
    if ismember(sProcess.options.direction.Value, [2 3])
        OPTIONS.GrangerDir = 'in';
        OutputFiles = cat(2, OutputFiles, bst_connectivity({sInputA.FileName}, {sInputA.FileName}, OPTIONS));
    end
end

author: brainstorm-tools


text:
function [p,str] = cluster_surf(varargin)
% Surface plot of clusters on a standard brain
%
%    WARNING: SURFACE_CUTAWAY, TOR_3D, AND CLUSTER_SURF ARE DEPRECATED.
%    THEY USE AN OLDER STYLE OF SURFACE RENDERING THAT IS VERY SLOW. THEY
%    STILL WORK, BUT OBJECT-ORIENTED CANLAB TOOLS HAVE SHIFTED TO USING
%    ADDBRAIN.M COMBINED WITH RENDER_ON_SURFACE() METHOD, WHICH USES
%    MATLAB'S ISOCOLORS FOR DRAMATICALLY FASTER COLOR RENDERING.
%
% :Usage:
% ::
%
%    [suface_handle,colorchangestring] = cluster_surf(varargin)
%
% :Inputs:
%
%   **CLUSTERS:**
%        clusters structures, as created in tor_extract_rois.m
%
%   **COLORS:**
%        cell array of colors for each cluster: {[1 0 0] [0 1 0] [0 0 1]}
%       if number of colors specified is greater than number of clusters
%       structures entered, n+1 and n+2 colors are overlap of 2 and overlap
%       of all clusters, respectively.
%
%   **SURFACE MAT FILE:**
%        file name of mat file containing brain surface vertices and faces
%       as created with isosurface.
%
%   **SPECIAL SURFACE KEYWORDS:**
%        special string: 'bg' 'hipp' (hcmp,thal,amy)
%        number of mm to plot from surface (mmdeep)
%
%          - Special keywords for sets of surfaces are:
%            left, right, bg, limbic, cerebellum, brainstem
%
%          - Other keywords: 'left' 'right' 'amygdala' 'thalamus' 'hippocampus' '
%            'midbrain'  'caudate'   'globus pallidus'  'putamen'  'nucleus accumbens'
%             'hypothalamus' 'cerebellum'
%
%   **EXISTING SURFACE HANDLE(S):**
%        Pass in handles for surface patches, created,
%        e.g., with addbrain.m.  This lets you be very flexible in the
%        surfaces you image onto.
%
%   **'colorscale':**
%        This scales colors by Z-scores of voxels if used
%          - Uses input color, unless also used with 'heatmap'
%          - Z scores should be in ROW vector
%          - use with 'normalize' to scale Z-scores between -1 and 1
%          - will also create transparent effects, mixing
%            blob color with existing surface color in linear
%            proportion to Z-scores
%
%   **'heatmap':**
%        Map Z-scores to surface colors
%          - Used WITH or instead of 'colorscale'
%          - Blobs can have a range of colors
%          - Use with REFERENCE RANGE option below to control scale
%          - solid colors entered as input will be ignored
%          - use with 'colormaps' option below to be flexible
%            in which color maps you apply.
%          - if 'colorscale' is also used, will produce transparent blobs.
%
%   **REFERENCE RANGE:**
%        reference Z-scores range, [zmin_act zmax_act
%        zmax_negact zmin_negact], e.g., [0 5 -5 0], use only
%        with 'heatmap' option
%        to get refZ from clusters, try:
%        ::
%
%            clZ = cat(2,clusters.Z);
%            refZ = [min(clZ(clZ > 0)) max(clZ) min(clZ(clZ < 0)) min(clZ)];
%
%   **'colormaps':**
%        - followed by custom [colors x 3] matrices for positive colors
%          and negative colors.
%        - matlab can create some: e.g., colormap summer, jet, etc.
%          others can be created with colormap_tor.m
%
%        color [0 1 1] (cyan) is reserved for the overlap color btwn cluster sets.
%
% :Examples:
% ::
% % ------------------------------------------------------------
% % ------------------------------------------------------------
%
%    P = 'C:\tor_scripts\3DheadUtility\canonical_brains\surf_single_subj_T1_gray.mat';
%    cluster_surf(tcl,acl,P,10,{[0 1 0] [1 0 0]},'colorscale','heatmap')
%
%    or P = h (surface handle) to use current surface in figure, and refZ
%    cluster_surf(tcl,acl,h,[3 5 -5 -3],10,{[0 1 0] [1 0 0]},'colorscale','heatmap')
%
% % A complete example visualizing a network on several surfaces
% % ------------------------------------------------------------
% % Load atlas
% atl = load_atlas('yeo17networks');
%
% % Select a network
% r = atlas2region(select_atlas_subset(atl, 1));
%
% % Check what it looks like on slices
% figure; montage(r);
%
% % Check what it looks like using surface() method
% % and render_on_surface()
%
% figure; surface(r);
%
% %% Run cluster_surf on several standard left-hemisphere surfaces
% % ------------------------------------------------------------
% create_figure('Left lateral surface', 2, 2);
% cluster_surf(r, 2, 'colors', {[1 1 0]}, 'hires left');
% view(275, 10); lightRestoreSingle;
% title('MNI Colin single-subj, addbrain(...''hires left'')');
%
% subplot(2, 2, 2);
% % Inflated surfaces: These use the pre-calculated mappers from MNI to an
% % inflated freesurfer cortical surface for each hemisphere:
% % which('lh.avgMapping_allSub_RF_ANTs_MNI152_orig_to_fsaverage.mat')
% % which('rh.avgMapping_allSub_RF_ANTs_MNI152_orig_to_fsaverage.mat')
% try
%     cluster_surf(r, 2, 'colors', {[1 1 0]}, 'fsavg_left');
%     view(275, 10); lightRestoreSingle;
%     title('Freesurfer average with Wu 2017 MNI->surf mapping');
% catch
%     disp('cannot find required surface file');
% end
%
% subplot(2, 2, 3);
% % Pial surface from Glasser_et_al_2016 Nature, HCP data
% % 'L.pial_MSMAll_2_d41_WRN_DeDrift.32k.mat';
% try
%     surf_han = addbrain('surface left');
%     cluster_surf(r, 2, 'colors', {[1 1 0]}, surf_han);
%     view(275, 10); lightRestoreSingle;
%     set(surf_han, 'FaceAlpha', 1);
%     title('HCP pial surface, addbrain(...''surface left'')');
% catch
%     title('cannot find required surface file');
% end
%
% subplot(2, 2, 4);
% % BigBrain surface, Amunts et al.
% % 'BigBrainSurfaceLeft.mat' in addbrain
% try
%     surf_han = addbrain('hires surface left');
%     cluster_surf(r, 2, 'colors', {[1 1 0]}, surf_han);
%     view(275, 10); lightRestoreSingle;
%     set(surf_han, 'FaceAlpha', 1);
%     title('Julich BigBrain, addbrain(...''hires surface left'')');
% catch
%     title('cannot find required surface file');
% end
%
% % ------------------------------------------------------------
%
% Example 2:
% Show all 32 L/R Hem Schafer/Yeo networks on the Julich BigBrain
%
% % ------------------------------------------------------------
%
% atl = load_atlas('yeo17networks');
%
% % Parse the atlas into a cell array with one cell per region
% % Each cell will be rendered in a different color
% clear region_cell
%
% for i = 1:num_regions(atl)
%
%     % Select a network
%     region_cell{i} = atlas2region(select_atlas_subset(atl, i));
%
% end
%
% create_figure('surface');
%
% % Seaborn colors....
% colors = seaborn_colors(num_regions(atl) ./ 2)';
% colors = [colors colors];
%
% % Or standard CANlab colors...
% colors = scn_standard_colors(num_regions(atl));
%
% surf_han = addbrain('bigbrain');
%
% cluster_surf(region_cell{:}, 3, 'colors', colors, surf_han);
% view(135, 12);
% lightRestoreSingle;
% set(surf_han, 'FaceAlpha', 1);
% title('HCP pial surface, addbrain(...''surface left'')');
%
% % ------------------------------------------------------------
% % ------------------------------------------------------------
% :More examples:
% ::
% % ------------------------------------------------------------
%
%    cluster_surf(cl,2,'heatmap');     % brain surface.  vertices colored @2 mm
%    cluster_surf(cl,2,'bg','heatmap');    % heatmap on basal ganglia
%    cluster_surf(cl,5,'left','heatmap');  % heatmap on left surface @5 mm
%    cluster_surf(cl,2,'right','heatmap');
%
%    % A multi-color, multi-threshold display on the cerebellum
%    colors = {[1 1 0] [1 .5 0] [1 .3 .3]};
%    tor_fig;
%    sh = cluster_surf(cl{3},colors(3),5,'cerebellum');
%    cluster_surf(cl{2},colors(2),5,sh);
%    cluster_surf(cl{1},colors(1),5,sh);
%
%    % Custom colormaps:
%    create_figure('Brain Surface'); cluster_surf(cl, 2, 'heatmap','left');
%
%    poscm = colormap_tor([.2 .2 .4], [1 1 0], [.9 .6 .1]);  %slate to orange to yellow
%    negcm = colormap_tor([0 0 1], [0 .3 1]);  % light blue to dark blue
%    create_figure('Brain Surface'); cluster_surf(cl, 2, 'heatmap', 'colormaps', poscm, negcm, 'left');
%
%    % Single-color transparent map (green):
%    cluster_surf(cl, 2, {[0 1 0]}, 'colorscale', p3(2), 'normalize');
%
% :See Also: addbrain, img2surf.m, surface() methods for objects,
% cluster_cutaways, render_on_surface()
%
%    WARNING: SURFACE_CUTAWAY, TOR_3D, AND CLUSTER_SURF
%    USE AN OLDER STYLE OF SURFACE RENDERING THAT IS VERY SLOW. THEY
%    STILL WORK, BUT OBJECT-ORIENTED CANLAB TOOLS HAVE SHIFTED TO USING
%    ADDBRAIN.M COMBINED WITH RENDER_ON_SURFACE() METHOD, WHICH USES
%    MATLAB'S ISOCOLORS FOR DRAMATICALLY FASTER COLOR RENDERING.
%    CLUSTER_SURF IS STILL USEFUL FOR INFLATED SURFACES, AND FOR RENDERING
%    MULTIPLE SETS OF REGIONS ON THE SAME BRAIN IN DIFFERENT COLORS, AND
%    WITH FLEXIBLE CONTROL OVER THE DISTANCE TO SURFACE.

% ..
%    Programmers' Notes
%    Created by Tor, a long time ago
%    updated Sept 2015 to keep up with matlab graphics and handle some weird
%    stuff with processing inputs.
%      - Figures are now scalars and we need to check for those first.
%      - Also changed default surface and colormap
%
%    Jan 2020: Updated to add 'noverbose' option, minor cosmetic code
%    cleanup
% ..

disp('Cluster_surf uses a type of surface rendering with cluster_surf that is deprecated')
disp('It still functions, but the new method based on isocolors is dramatically faster')

% -----------------------------------------------------------------------
%    set up input arguments and defaults
% -----------------------------------------------------------------------
mmdeep = 10;
cscale = 0;
heatm = 0;
viewdeg = [135 30];
cl = [];

P = which('surf_spm2_brain.mat');
P = which('surf_spm2_brain_1mm.mat');

% default color maps
% These match fmridisplay:
poscm = colormap_tor([1 .5 0], [1 1 0]);
negcm = colormap_tor([0 0 1], [0 1 1]);

actcolors = [];  % used with heatmap
donormalize = 0; % used with colorscale
adjust_var = [];

doverbose = true;
if any(strcmp(varargin,'noverbose')), doverbose = false; end % do first

% -----------------------------------------------------------------------
%    optional inputs
% -----------------------------------------------------------------------
clind = 1;
for i = 1:length(varargin)

    if isempty(varargin{i})
        % ignore it

    elseif isstruct(varargin{i}) || isa(varargin{i}, 'region')
        cl{clind} = varargin{i}; clind = clind+1;

    elseif iscell(varargin{i}), mycolors = varargin{i};

    elseif isstr(varargin{i})

        if strcmp(varargin{i},'noverbose'), doverbose = false;

        elseif strcmp(varargin{i},'colorscale'), cscale = 1;

        elseif strcmp(varargin{i},'normalize'), donormalize = 1;

        elseif strcmp(varargin{i},'heatmap'), heatm = 1;

        elseif strcmp(varargin{i},'colormaps')
            if doverbose, disp('Using custom color maps.'); end
            poscm = varargin{i + 1}; varargin{i + 1} = [];
            negcm = varargin{i + 2};  varargin{i + 2} = [];

        elseif  strcmp(varargin{i},'left')
            P = which('surf_spm2_left.mat'); %which('surf_single_subj_grayL.mat');
            viewdeg = [90 0];

        elseif  strcmp(varargin{i},'right')
            P = which('surf_spm2_right.mat'); %which('surf_single_subj_grayR.mat');
            viewdeg = [270 0];

        elseif  strcmp(varargin{i},'hires left')
            P = which('surf_spm2_brain_left.mat'); %which('surf_single_subj_grayL.mat');
            viewdeg = [90 0];

        elseif  strcmp(varargin{i},'hires right')
            P = which('surf_spm2_brain_right.mat'); %which('surf_single_subj_grayR.mat');
            viewdeg = [270 0];
        elseif strcmp(varargin{i},'fsavg_right') % uses freesurfer inflated brain
            % with Thomas Yeo group's RF_ANTs mapping
            % from MNI to Freesurfer. (https://doi.org/10.1002/hbm.24213)
            P = which('surf_freesurf_inflated_Right.mat');
            viewdeg = [270 0];
            adjust_var = 'fsavg_right'; % varargin for getVertexColors
        elseif strcmp(varargin{i},'fsavg_left') % uses freesurfer inflated brain
            P = which('surf_freesurf_inflated_Left.mat');
            viewdeg = [90 0];
            adjust_var = 'fsavg_left'; % varargin for getVertexColors
        else P = varargin{i};
        end

    elseif isnumeric(varargin{i}) && isscalar(varargin{i})
        % it's a number, mm deep to render (any face within mmdeep mm of a significant voxel will
        % be colored)
        mmdeep = varargin{i};

    elseif any(ishandle(varargin{i})) && ~all(ishandle(varargin{i}))
        % Note: some weird things are happening as fig handles are class
        % double, and not recognized as figs, but are figure handles.
        % hopefully this will work with strcmp()
        % do nothing, but fig handles may be passed in inadvertently?
        % sometimes passing in integers is interpreted as fig handles...
        % For example, passing in 0 as a double is interpreted as a root
        % graphics object. Not a good design feature of new graphics in
        % Matlab...

        % Do nothing here. Could be all numeric input.

    elseif any(ishandle(varargin{i})) && any(strcmp(get(varargin{i}, 'Type'), 'figure'))
        % all are handles, one is a figure

        disp('You passed in a figure handle.')

    elseif  all(ishandle(varargin{i})) && all(isa(varargin{i}, 'matlab.graphics.primitive.Patch'))  %% all(~isa(varargin{i}, 'matlab.ui.Figure'))
        % all are handles, and patches
        % OLD: Pre-2014:  all(ishandle(varargin{i})) && all(varargin{i} ~= round(varargin{i}))
        if doverbose, disp('Found surface patch handles - plotting on existing surfaces.'); end
        P = varargin{i}; % handle(s) for existing surface

        % get rid of later calls to other surfaces
        wh = find(strcmp(varargin,'left') | strcmp(varargin,'right') | strcmp(varargin,'hires left') | strcmp(varargin,'hires right'));
        for jj = 1:length(wh), varargin{wh(jj)} = []; end

    elseif any(ishandle(varargin{i})) && ~all(isa(varargin{i}, 'matlab.graphics.primitive.Patch'))
        % OLD:  Pre-2014: any(ishandle(varargin{i})) && all(varargin{i} ~= round(varargin{i}))
        % Note: some weird things are happening as fig handles are class
        % double, and not recognized as figs, but are figure handles.
        % hopefully this will work with strcmp()
        if doverbose
            disp('Found existing surface patch handles, but some are invalid.  Check handles.');
            error('Exiting')
        end

    elseif length(varargin{i}) > 1   % it's a vector
        refZ = varargin{i};

    else    % no idea
        warning('cluster_surf: Unknown input');

    end

end

if ~exist('cl', 'var') || isempty(cl)
    if doverbose, disp('cluster_surf.m: No clusters to plot.  Try addbrain for brain surfaces with no activation.'); end
    p = []; str = [];
    return
end

if ~exist('mycolors', 'var')
    mycolors = scn_standard_colors(length(cl));
end

if isempty(P)
    disp(['Cannot find: ' P]);
    P = spm_get(1,'*mat','Choose brain surface file');
end

if doverbose

    disp('cluster_surf')
    disp('___________________________________________')
    fprintf('\t%3.0f cluster structures entered\n',length(cl))
    disp('  Colors are:')

    for i = 1:length(mycolors)
        disp([' ' num2str(mycolors{i})])
    end

end

if length(mycolors) > length(cl)
    if doverbose, disp([' overlap color is ' num2str(mycolors{length(cl)+1})]), end
    ovlc = ['[' num2str(mycolors{length(cl)+1}) ']'];
else
    ovlc = '[0 1 1]';
end

if length(mycolors) > length(cl)+1 && length(cl) > 2
    if doverbose, disp([' all overlap color is ' num2str(mycolors{length(cl)+2})]), end
    aovlc = ['[' num2str(mycolors{length(cl)+2}) ']'];
else
    aovlc = '[1 1 1]';
end

%disp([' Surface stored in: ' P])
if doverbose
    fprintf(' Building XYZ coord list\n');
end

% -------------------------------------------------------------------------
% * build xyz list
%
% also get cscale values for each coordinate, and alphascale values if both
% heatmap and colorscale options are entered
% -------------------------------------------------------------------------
for i = 1:length(cl) % each cell is a whole vector of cl, not a single region

    xyz{i} = cat(2,cl{i}.XYZmm)';

    if cscale || heatm
        for j = 1:length(cl{i})
            if size(cl{i}(j).Z,1) > size(cl{i}(j).Z,2)
                cl{i}(j).Z = cl{i}(j).Z';
            end
        end
        Z{i} = cat(2,cl{i}.Z)';

        % order voxels from lowest to highest, so that peak colors
        % appear because they are plotted last
        tmp = [xyz{i} Z{i}];
        tmp = sortrows(tmp,4);
        xyz{i} = tmp(:,1:3);
        Z{i} = tmp(:,4);1 ./ mad(abs(Z{i}));

        if cscale
            if donormalize
                Z{i} = Z{i} ./ max(Z{i});
            end

            if heatm
                % treat colorscale as alpha scaling to add transparent blobs
                % (preserve existing surface; good for isosurface objects)

                Za = abs(Z{i});
                a = prctile(Za, 85);  % midpoint of Z{i} defines transparency 0.5
                b = 1./mad(Za); % multiplier for Z for sigmoid

                sZ = 1 ./ (1 + exp(-b*(Za-a)));

                % fix, if all constant
                sZ(isnan(sZ)) = 1;

                alphascale{i} = sZ;

                % this is further adusted based on radius
                alphascale{i} = 5 * alphascale{i} ./ mmdeep^3; % should be 3?
                %alphascale

            end

        else
            % if heat map only, set mycolor{1} = [1 1 1]
            mycolors{1} = [1 1 1];
        end
    end
end

% ------------------------------------------------------------
% for heatmap option: get actcolors
% -------------------------------------------------------------
if heatm
    if doverbose, fprintf(' Getting heat-mapped colors\n'); end
    if exist('refZ') == 1
        actcolors = get_actcolors(Z, refZ, poscm, negcm);
    else
        actcolors = get_actcolors(Z, [], poscm, negcm);
    end
end

% rearrange Z to map the negative peak later
if exist('Z', 'var')
    for i = 1:numel(Z)
        neg_idx = Z{i}<0;
        xyz{i}(neg_idx,:)=flipud(xyz{i}(neg_idx,:));
        Z{i}(neg_idx)=flipud(Z{i}(neg_idx));
        actcolors{i}(neg_idx,:)=flipud(actcolors{i}(neg_idx,:));
    end
end

% -------------------------------------------------------------------------
% * build function call
% -------------------------------------------------------------------------
if doverbose, fprintf(' Building color change function call\n'); end

if length(cl) > 2 && exist('aovlc') == 1
    str = ['[c,alld] = getVertexColors(xyz{1},p,mycolors{1},[.5 .5 .5],' num2str(mmdeep) ',''ovlcolor'',' ovlc ',''allcolor'',' aovlc];
else
    str = ['[c,alld] = getVertexColors(xyz{1},p,mycolors{1},[.5 .5 .5],' num2str(mmdeep) ',''ovlcolor'',' ovlc];
end

if heatm
    str = [str ',''colorscale'',actcolors{1}'];
    if cscale
        % treat colorscale as alpha scaling to add transparent blobs
        str = [str ',''alphascale'',alphascale{1}'];
    end
elseif cscale
    % cscale alone - treat cscale as color-mapping index for single
    % colors in actcolors
    str = [str ',''colorscale'',Z{1}'];
end

for i = 2:length(cl)
    str = [str ',''vert'',xyz{' num2str(i) '},mycolors{' num2str(i) '}'];
    if heatm
        str = [str ',''colorscale'',actcolors{' num2str(i) '}'];
        if cscale
            % treat colorscale as alpha scaling to add transparent blobs
            str = [str ',''alphascale'',alphascale{' num2str(i) '}'];
        end
    elseif cscale
        % cscale alone - treat cscale as color-mapping index for single
        % colors in actcolors
        str = [str ',''colorscale'',Z{' num2str(i) '}'];
    end
end

if ~isempty(adjust_var)
    str = [str ',''' adjust_var ''''];
end

if ~doverbose
    str = [str ', ''noverbose'''];
end

str = [str ');'];

%p = P(end);
%co = get(p, 'FaceVertexCData');
if exist('alphascale','var')
    alphascale{1} = alphascale{1} * 12;
end
%eval(str)
%set(p, 'FaceVertexCData', co);

% -------------------------------------------------------------------------
% * run brain surface
% -------------------------------------------------------------------------
if ishandle(P)      % no input file, use existing handle
    if doverbose
        fprintf(' Using existing surface image\n');
        fprintf(' Running color change.\n');
    end
    for i = 1:length(P)
        p = P(i);
        if doverbose, disp([' eval: ' str]), end
        eval(str)
    end
else
    % we have either an input file or a special string ('bg')
    if doverbose, fprintf(' Loading surface image\n'); end
    [dtmp,ftmp,etmp]=fileparts(P);

    if strcmp(etmp,'.mat')

        load(P);

        %%figure
        p = patch('Faces',faces,'Vertices',vertices,'FaceColor',[.5 .5 .5], ...
            'EdgeColor','none','SpecularStrength',.2,'FaceAlpha',1,'SpecularExponent',200);
        lighting gouraud;camlight right
        axis image;
        lightRestoreSingle(gca);
        %myLight = camlight(0,0);set(myLight,'Tag','myLight');
        %set(gcf, 'WindowButtonUpFcn', 'lightFollowView');lightFollowView

        view(viewdeg(1),viewdeg(2));
        drawnow

        % -------------------------------------------------------------------------
        % * run color change
        % -------------------------------------------------------------------------
        if doverbose
            fprintf(' Running color change.\n');
            disp([' eval: ' str])
        end
        eval(str);

        % this for subcortex stuff
    elseif strcmp(P,'bg')
        p = [];
        myp = addbrain('caudate');p = [p myp];
        run_colorchange(myp,str,xyz,mycolors);

        myp = addbrain('globus pallidus');p = [p myp];
        run_colorchange(myp,str,xyz,mycolors);

        myp = addbrain('putamen');p = [p myp];
        run_colorchange(myp,str,xyz,mycolors);

        set(myp,'FaceAlpha',1);

        axis image; axis vis3d; lighting gouraud; lightRestoreSingle(gca)

    elseif strcmp(P,'limbic')
        p = [];

        myp = addbrain('amygdala');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('hypothalamus');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('hippocampus');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('thalamus');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('nucleus accumbens');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);

        myp = addbrain('left');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        set(myp,'FaceAlpha',1);

        axis image; axis vis3d; lighting gouraud; lightRestoreSingle(gca)

    elseif strcmp(P,'brainstem')
        p = [];

        myp = addbrain('thalamus');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('hypothalamus');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('brainstem');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        %myp = addbrain('caudate');p = [p myp];
        %run_colorchange(myp,str,xyz, mycolors, actcolors);

        view(90,10); axis image; axis vis3d; lighting gouraud; lightRestoreSingle(gca)

    elseif strcmp(P,'subcortex')
        p = [];

        myp = addbrain('thalamus');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('hypothalamus');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('brainstem');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('amygdala');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('hippocampus');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('thalamus');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('nucleus accumbens');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('caudate');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        myp = addbrain('putamen');p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);
        %         myp = addbrain('globus pallidus');p = [p myp];  %Seems to be missing from canonical folder
        %         run_colorchange(myp,str,xyz, mycolors, actcolors);

        view(90,10); axis image; axis vis3d; lighting gouraud; lightRestoreSingle(gca)

    elseif strcmp(P,'cerebellum') || strcmp(P,'amygdala') || strcmp(P,'hypothalamus') ...
            || strcmp(P,'thalamus') || strcmp(P,'midbrain') || strcmp(P,'caudate') ...
            || strcmp(P,'globus pallidus') || strcmp(P,'putamen') || strcmp(P,'nucleus accumbens') ...
            || strcmp(P,'hippocampus')
        % this uses addbrain and works with any of its keywords
        p = [];

        myp = addbrain(P);p = [p myp];
        run_colorchange(myp,str,xyz, mycolors, actcolors);

        view(90,10); axis image; axis vis3d; lighting gouraud; lightRestoreSingle(gca)

    else
        error('Must input mat surf file or img file to convert to surf')
    end

end % if ishandle

lighting gouraud
lightRestoreSingle(gca);
material dull
axis off
set(gcf,'Color','w')
%scn_export_papersetup(400);  % this will mess up movies!!

if doverbose
    disp('Finished!')
    disp('___________________________________________')
end

end % main function

function actcolor = get_actcolors(datavaluesets, refZ, poscm, negcm)

actcolor = map_data_to_colormap(datavaluesets, poscm, negcm, refZ);

end

function run_colorchange(myp, str, xyz, mycolors, actcolors)

set(myp,'FaceAlpha',1);

for i = 1:length(myp)
    p = myp(i);

    % get original color
    origcolor = get(p,'FaceColor');

    % color change
    eval(str);

    % change non-active back to original color
    vdat = get(p,'FaceVertexCData');
    wh = find(all(vdat == .5, 2));
    vdat(wh,:) = repmat(origcolor,length(wh),1);
    set(p,'FaceVertexCData',vdat);

end
p = myp;
lighting gouraud; lightRestoreSingle(gca);

end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


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

%%
clear; close all; clc;

%% Syntax
% |L=intesselation(X,TES,XI);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function lambda = circulant2_eigenvalues ( n, x )

%*****************************************************************************80
%
%% CIRCULANT2_EIGENVALUES returns the eigenvalues of the CIRCULANT2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    17 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of rows and columns of A.
%
%    Output, complex LAMBDA(N,1), the eigenvalues.
%
  lambda = zeros ( n, 1 );

  x = zeros ( n, 1 );
  for i = 1 : n
    x(i,1) = i;
  end

  w(1:n,1) = c8vec_unity ( n );

  lambda(1:n,1) =  x(n,1);
  for i = n-1 : -1 : 1
    lambda(1:n,1) = lambda(1:n,1) .* w(1:n,1) + x(i,1);
  end

  return
end


text:
function f3=f3(x)
Bound=[-100 100];

if nargin==0
    f3 = Bound;
else
    f3=sum(cumsum(x).^2);
end