text:
function el = export_CT_image_module_field(args)
%"export_CT_image_module_field"
%   Given a single scan, return a properly populated CT_image module tag
%   for use with any Composite Image IOD.  See CT_image_module_tags.m.
%
%   For speed, tag must be a decimal representation of the 8 digit
%   hexidecimal DICOM tag desired, ie instead of '00100010', pass
%   hex2dec('00100010');
%
%   Arguments are passed in a structure, arg:
%       arg.tag         = decimal tag of field to fill
%       arg.data        = CERR structure(s) to fill from
%       arg.template    = an empty template of the module created by the
%                         function build_module_template.m
%
%   This function requires arg.data = {scanInfoS, scanS};
%
%JRA 06/19/06
%NAV 07/19/16 updated to dcm4che3
%   replaced ml2dcm_Element to data2dcmElement
%
%Usage:
%   el = export_CT_image_module_field(args)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

%Init output element to empty.
el = [];

%Unpack input parameters.
tag         = args.tag;
scanInfoS   = args.data{1};
scanS       = args.data{2};
template    = args.template;

switch tag
    %Class 1 Tags -- Required, must have data.
    case 524296     %0008,0008 Image Type
        data = {'ORIGINAL', 'PRIMARY', 'AXIAL'};
        el = data2dcmElement(data, tag);

    case 2621442    %0028,0002 Samples per Pixel
        data = 1;               %1 image plane in all CT/MR images.
        el = data2dcmElement(data, tag);

    case 2621444    %0028,0004 Photometric Interpretation
        data = 'MONOCHROME2';   %CT/MR have 0 black, maxVal white.
        el = data2dcmElement(data, tag);

    case 2621456    %0028,0010 Rows
        data = scanInfoS.sizeOfDimension1;
        el = data2dcmElement(data, tag);

    case 2621457    %0028,0011 Columns
        data = scanInfoS.sizeOfDimension2;
        el = data2dcmElement(data, tag);

    case 2621696    %0028,0100 Bits Allocated
        data = 16;              %C.8.2.1.1.4 of PS 3.3 - 2006
        el = data2dcmElement(data, tag);

    case 2621697    %0028,0101 Bits Stored
                                %C.8.2.1.1.4 of PS 3.3 - 2006
                                %Sloppy, consider revising.
        bools = [scanS.scanInfo.zValue] == scanInfoS.zValue;
        vals = scanS.scanArray(:,:,bools);
        maxV = max(vals(:));
        log2s = log2(double(maxV));
        mostSignificantBit = floor(max(log2s)) + 1;
        data = max(mostSignificantBit, 12);
        data = min(data, 16);
        el = data2dcmElement(data, tag);

    case 2621698    %0028,0102 High Bit
                                %C.8.2.1.1.4 of PS 3.3 - 2006
                                %Sloppy, consider revising.
        bools = [scanS.scanInfo.zValue] == scanInfoS.zValue;
        vals = scanS.scanArray(:,:,bools);
        maxV = max(vals(:));
        log2s = log2(double(maxV));
        mostSignificantBit = floor(max(log2s)) + 1;
        data = max(mostSignificantBit, 12);
        data = min(data, 16);
        data = data - 1;
        el = data2dcmElement(data, tag);

    case 2625618    %0028,1052 Rescale Intercept
        ctO = scanInfoS.CTOffset;
        data = -ctO;           %CERR exports stored values as 1*HU + CTOffset.
        el = data2dcmElement(data, tag);

    case 2625619    %0028,1053 Rescale Slope
        %data = 1;
        %data = scanInfoS.rescaleSlope;
        data = args.data{3}; %APA: factor for conversion to uint16 for modalities other than CT
        bools = [scanS.scanInfo.zValue] == scanInfoS.zValue;
        data = data(bools);
        data = data(1); % handle multiple slices at same z-location
        el = data2dcmElement(data, tag);

    %Class 2 Tags -- Must be present, can be NULL.
    case 1572960    %0018,0060 KVP

        %el = org.dcm4che3.data.Attributes;
        %el.setString(tag, template.getVR(tag), template.getString(tag));
        %data = args.data{1};
        data = [];
        el = data2dcmElement(data, tag);

    case 2097170    %0020,0012 Acqusition Number

        %el = org.dcm4che3.data.Attributes;
        %el.setString(tag, template.getVR(tag), template.getString(tag));
        %data = args.data{1};
        data = [];
        el = data2dcmElement(data, tag);

    %Class 3 Tags -- presence is optional, currently undefined.
    case 1572898    %0018,0022 Scan Options
    case 1573008    %0018,0090 Data Collection Diameter
    case 1577216    %0018,1100 Reconstruction Diameter
    case 1577232    %0018,1110 Distance Source to Detector
    case 1577233    %0018,1111 Distance Source to Patient
    case 1577248    %0018,1120 Gantry/Detector Tilt
    case 1577264    %0018,1130 Table Height
    case 1577280    %0018,1140 Rotation Direction
    case 1577296    %0018,1150 Exposure Time
    case 1577297    %0018,1151 X-ray Tube Current
    case 1577298    %0018,1152 Exposure
    case 1577299    %0018,1153 Exposure in microAs
    case 1577312    %0018,1160 Filter Type
    case 1577328    %0018,1170 Generator Power
    case 1577360    %0018,1190 Focal Spot
    case 1577488    %0018,1210 Convolution Kernal
    case 1610501    %0018,9305 Revolution Time
    case 1610502    %0018,9306 Single Collimation Width
    case 1610503    %0018,9307 Total Collimation Width
    case 1610505    %0018,9309 Table Speed
    case 1610512    %0018,9310 Table Feed per Rotation
    case 1610513    %0018,9311 CT Pitch Factor
    case 1610531    %0018,9323 Exposure Modulation Type
    case 1610532    %0018,9324 Estimated Dose Saving
    case 1610565    %0018,9345 CTDIvol

    %Class 1C Tags

    %Class 2C Tags

    otherwise
        warning(['No methods exist to populate DICOM image_pixel module field ' dec2hex(tag,8) '.']);
        return;
end


text:
function eventLog = ma_executeDVManeuver(maneuverEvent, initialState, eventNum, celBodyData)
%ma_executeDVManeuver Summary of this function goes here
%   Detailed explanation goes here

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Do the maneuver according to its sub-type
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    type = maneuverEvent.maneuverType;
    value = maneuverEvent.maneuverValue;
    thruster = maneuverEvent.thruster;
    switch type
        case 'dv_inertial'
            eventLog = ma_executeDVManeuver_dv_inertial(value, thruster, initialState, eventNum);
        case 'dv_orbit'
            eventLog = ma_executeDVManeuver_dv_orbit(value, thruster, initialState, eventNum);
        case 'finite_inertial'
            eventLog = ma_executeDVManeuver_finite_inertial(value, thruster, initialState, eventNum, celBodyData);
        case 'finite_steered'
            eventLog = ma_executeDVManeuver_finite_steered(value, thruster, initialState, eventNum, celBodyData);
        case 'dv_inertial_spherical'
            eventLog = ma_executeDVManeuver_dv_inertial_spherical(value, thruster, initialState, eventNum);
        case 'dv_orbit_spherical'
            eventLog = ma_executeDVManeuver_dv_orbit_spherical(value, thruster, initialState, eventNum);
        case 'finite_inertial_spherical'
            eventLog = ma_executeDVManeuver_finite_inertial_spherical(value, thruster, initialState, eventNum, celBodyData);
        case 'finite_steered_spherical'
            eventLog = ma_executeDVManeuver_finite_steered_spherical(value, thruster, initialState, eventNum, celBodyData);
        case 'circularize'
            eventLog = ma_executeDVManeuver_circularize(thruster, initialState, eventNum, celBodyData);
        otherwise
            error(['Did not recongize maneuver of type ', type]);
    end
end

author: Arrowstar

repo: ksptot

save_path: github-repos/MATLAB/Arrowstar-ksptot


text:
function varargout = plotcoeffs(f, varargin)
%PLOTCOEFFS   Display coefficients of the columns, rows and tubes of a
%   CHEBFUN3 object.
%
% See also CHEBFUN/PLOTCOEFFS and CHEBFUN2/PLOTCOEFFS.

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

% Empty check.
if ( isempty(f) )
    varargout = { [] };
    return
end

% Store the hold state of the current axis:
holdState = ishold;

% Get low rank representation of f:
[~, fCols, fRows, fTubes] = tucker(f);

% PLOTCOEFFS of cols:
ax1 = subplot(1, 3, 1);
plotcoeffs(fCols, varargin{:});
ylim1 = ylim(gca);
% Remove labels from 1D plotcoeff:
xlabel(gca, '')
ylabel(gca, 'Magnitude of coefficient')
title('Cols')

% PLOTCOEFFS of rows:
ax2 = subplot(1, 3, 2);
plotcoeffs(fRows, varargin{:});
ylim2 = ylim(gca);
% Put an xlabel only for the plot in the middle:
if isa(f.cols.funs{1}.onefun,'trigtech')
    xlabel(gca, 'Wave number')
else
    xlabel(gca, 'Degree of Chebyshev polynomial')
end
ylabel(' ')
title('Rows')

% PLOTCOEFFS of tubes:
ax3 = subplot(1, 3, 3);
plotcoeffs(fTubes, varargin{:});
ylim3 = ylim(gca);
xlabel(' ')
ylabel(' ')
title('Tubes')

% Find a proper ylim for all the three subplots:
yLims = [ylim1; ylim2; ylim3];
ylimNew = [min(yLims(:, 1)), max(yLims(:, 2))];

% Set the ylim of the plots again to be similar.
ylim(ax1, ylimNew);
ylim(ax2, ylimNew);
ylim(ax3, ylimNew);

% Return hold state to what it was before:
if ( ~holdState )
    hold off
end

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function order = order_table ( rule )

%*****************************************************************************80
%
%% ORDER_TABLE returns the order of a Lebedev rule.
%
%  Modified:
%
%    13 September 2010
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Vyacheslav Lebedev, Dmitri Laikov,
%    A quadrature formula for the sphere of the 131st
%    algebraic order of accuracy,
%    Russian Academy of Sciences Doklady Mathematics,
%    Volume 59, Number 3, 1999, pages 477-481.
%
%  Parameters:
%
%    Input, integer RULE, the index of the rule, between 1 and 65.
%
%    Output, integer ORDER, the order of the rule.
%
  rule_max = 65;

  table = [ ...
       6,   14,   26,   38,   50,   74,   86,  110,  146,  170, ...
     194,  230,  266,  302,  350,  386,  434,  482,  530,  590, ...
     650,  698,  770,  830,  890,  974, 1046, 1118, 1202, 1274, ...
    1358, 1454, 1538, 1622, 1730, 1814, 1910, 2030, 2126, 2222, ...
    2354, 2450, 2558, 2702, 2810, 2930, 3074, 3182, 3314, 3470, ...
    3590, 3722, 3890, 4010, 4154, 4334, 4466, 4610, 4802, 4934, ...
    5090, 5294, 5438, 5606, 5810 ]';

  if ( rule < 1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ORDER_TABLE - Fatal error!\n' );
    fprintf ( 1, '  RULE < 1.\n' );
    error ( 'ORDER_TABLE - Fatal error!' );
  elseif ( rule_max < rule )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'ORDER_TABLE - Fatal error!\n' );
    fprintf ( 1, '  RULE_MAX < RULE.\n' );
    error ( 'ORDER_TABLE - Fatal error!' );
  end

  order = table(rule);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [X] = spm_dot(X,x,i)
% Multidimensional dot (inner) product
% FORMAT [Y] = spm_dot(X,x,[DIM])
%
% X   - numeric array
% x   - cell array of numeric vectors
% DIM - dimensions to omit (asumes ndims(X) = numel(x))
%
% Y  - inner product obtained by summing the products of X and x along DIM
%
% If DIM is not specified the leading dimensions of X are omitted.
% If x is a vector the inner product is over the leading dimension of X
%
% See also: spm_cross
%__________________________________________________________________________
% Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_dot.m 7314 2018-05-19 10:13:25Z karl $

% initialise dimensions
%--------------------------------------------------------------------------
if iscell(x)
    DIM = (1:numel(x)) + ndims(X) - numel(x);
else
    DIM = 1;
    x   = {x};
end

% omit dimensions specified
%--------------------------------------------------------------------------
if nargin > 2
    DIM(i) = [];
    x(i)   = [];
end

% inner product using recursive summation (and bsxfun)
%--------------------------------------------------------------------------
for d = 1:numel(x)
    s         = ones(1,ndims(X));
    s(DIM(d)) = numel(x{d});
    X         = bsxfun(@times,X,reshape(full(x{d}),s));
    X         = sum(X,DIM(d));
end

% eliminate singleton dimensions
%--------------------------------------------------------------------------
X = squeeze(X);

return

% NB: alternative scheme using outer product
%==========================================================================

% outer product and sum
%--------------------------------------------------------------------------
x      = spm_cross(x);
s      = ones(1,ndims(X));
S      = size(X);
s(DIM) = S(DIM);
x      = reshape(full(x),s);
X      = bsxfun(@times,X,x);
for d  = 1:numel(DIM)
    X  = sum(X,DIM(d));
end
X      = squeeze(X);

author: spm


text:
function [t1,b0] = dtiMakeAnalyze(b0PathName, t1PathName, t1AnalyzePathName, outPathName)
%
% [t1,b0] = dtiMakeAnalyze(b0PathName, t1PathName, t1AnalyzePathName, outPathName)
%
% Takes a set of images from a DTI sequence and saves them out as
% analyze-format files. It will also save out a transform matrix that will
% coregister the DTI images to another scan series (typically a T1 scan
% from the same session).
% Also, if the t1 analyze file is provided and includes a talairach file,
% the transform matrix will also include the affine component of the
% talairach transform.
%
% HISTORY:
% 2003.08.?? RFD (bob@white.stanford.edu) wrote it.

useManual = false;

if ~exist('b0PathName','var') | isempty(b0PathName)
   [f, p] = uigetfile({'*.001','I-files (*.001)';'*.*','All files'}, 'Select one of the B0 I-files...');
   if(isnumeric(f)) error('Need a B0 file to continue...'); end
   b0PathName = fullfile(p, f);
   disp(b0PathName);
end
b0.fname = b0PathName;

if ~exist('t1PathName','var') | isempty(t1PathName)
   [f, p] = uigetfile({'*.001','I-files (*.001)';'*.*','All files'}, 'Select one of the T1 I-files...');
   if(isnumeric(f)) error('Need a T1 file to continue...'); end
   t1PathName = fullfile(p, f);
   disp(t1PathName);
end
t1.fname = t1PathName;

if ~exist('t1AnalyzePathName','var') | isempty(t1AnalyzePathName)
   [f, p] = uigetfile({'*.hdr','Analyze files (*.hdr)';'*.*','All files'},...
       'Select the T1 Analyze file...');
   if(isnumeric(f)) error('Need a T1 file to continue...'); end
   t1AnalyzePathName = fullfile(p, f);
   disp(t1AnalyzePathName);
end
t1.analyzeFname = t1AnalyzePathName;
[p,f,e] = fileparts(t1AnalyzePathName);
t1.analyzeBaseFname = fullfile(p,f);

% Compute the transform that would happen if we made this into an analyze
% format file using our code, which imposes a particular data orientation
% (axial, neurological convention). Note that the transofrm returned here
% should consits of just flips and mirror-reversals. There will be no
% scales, skews or translations.
[b0.cannonical, b0.baseFname, b0.mmPerVox, b0.imDim, b0.notes] = computeCannonicalXformFromIfile(b0.fname);
[t1.cannonical, t1.baseFname, t1.mmPerVox, t1.imDim, t1.notes] = computeCannonicalXformFromIfile(t1.fname);

if(~exist('outPathName','var') | isempty(outPathName))
    outPathName = fullfile(fileparts(fileparts(b0.baseFname)), 'dti_analyses');
end
if(~exist('outPathName','dir'))
    [p,f,e] = fileparts(outPathName);
    mkdir(p, [f e]);
end

swapXY = [0 1 0 0; 1 0 0 0; 0 0 1 0; 0 0 0 1];

t1.talairach = load([t1.analyzeBaseFname '_talairach']);
% We also want a transform that will go into ac-pc space (ie.
% Talairach, but without any scaling). We can compute that by
% adjusting the scale factors in the talairach transform.
t1.talXform = t1.talairach.vol2Tal.transRot'*swapXY;
p = spm_imatrix(t1.talXform);
% We also have to rescale the translations.
scaleDiff = t1.mmPerVox./p(7:9);
p(1:3) = p(1:3).*scaleDiff;
p(7:9) = t1.mmPerVox;
t1.acpcXform = spm_matrix(p);

% Now compute the transform from the standard scanner coordinate space
% to the image coordinate space (eg. scanner2img*[X Y Z 1]' will convert
% the X Y Z scanner coords to indices into the img.)
b0.scanner2img = XformFromIfile(b0.baseFname,'volume');
t1.scanner2img = XformFromIfile(t1.baseFname,'volume');

% create and save the T1 analyze file .mat transform, computed from the
% talairach file.
t1Tal = t1.talairach.vol2Tal.transRot'*swapXY;
M_acpc = t1.acpcXform;
M_tal = t1.talXform;
M = t1.talXform;
save(t1.analyzeBaseFname, 'M', 'M_tal', 'M_acpc');
% Also, make sure the T1 .hdr origin is set to the AC (Talairach 0,0,0)
t1.origin = inv(t1Tal)*[0,0,0,1]';
t1.origin = t1.origin(1:3)';
[dim vox scale type offset origin descrip] = spm_hread(t1.analyzeBaseFname);
spm_hwrite(t1.analyzeBaseFname, dim, vox, scale, type, offset, t1.origin, descrip);

% We use a pre-multiply convention (newCoords = M*oldCoords), so read
% the matrix transforms from right to left:
% 1. undo the b0 cannonical xform (inv(b0.cannonical))
% 2. b0 image space o scanner space (inv(b0.scanner2img))
% 3. scanner space to t1 image space (t1.scanner2img)
% 4. t1 image space to cannonical x-formed t1 image space(t1.cannonical)
% 5. t1 cannonical image space to t1 standard space (t1Hdr.mat)
%
% Note: the SPM convention seems to be to swap X and Y, so we have to undo
% that swap from the t1 transforms when computing ours. Then, we swap our
% transform before saving it. Why? Because it makes everything work. But, I
% really don't understand it fully. Maybe the swap is happening in myCinterp3
% (called from dtiGetSlice)?
M = t1.cannonical * swapXY * t1.scanner2img * inv(b0.scanner2img) * b0.cannonical;
M_tal = t1Tal * M;
M = M*swapXY;
M_tal = M_tal*swapXY;
b0.M = M;
b0.M_tal = M_tal;
% Set the origin to the AC (Talairach 0,0,0)
b0.origin = inv(b0.M_tal)*[0,0,0,1]';
b0.origin = b0.origin(1:3)';

% Build B0, FA and vector images
[b0.img, fa, add, ci, meanDiff, maxEig, medEig, minEig, maxVec] = ...
    dtiLoadTensorcalc(fullfile(fileparts(b0.baseFname), 'Tensor.float.001'));
b0.img = permute(b0.img, [2,1,3]); % x and y are permuted in the Tensor.float file
[b0.imgCannonical, mmPerVoxNew] = applyCannonicalXform(b0.img, b0.cannonical, b0.mmPerVox);
hdr = saveAnalyze(b0.imgCannonical, fullfile(outPathName,'B0'), mmPerVoxNew, b0.notes, b0.origin);
[img, mmPerVoxNew] = applyCannonicalXform(permute(fa*1000, [2,1,3]), b0.cannonical, b0.mmPerVox);
hdr = saveAnalyze(img, fullfile(outPathName,'FA'), mmPerVoxNew, b0.notes, b0.origin);
% The first two dims will get permuted when we write the data out to
% the analyze file (to put them in standard analyze axial form). But,
% we need to permute them here because tensorcalc permuted them
% relative to the original Ifiles.
maxVec = permute(maxVec,[2 1 3 4]);
% Since we just permuted voxel coordinates, we should also swap X and Y,
% but it turns out that the tensorcalc convention makes X
% anterior-posterior and Y left-right, but we want them to be the other
% way to follow Talairach naming convention. So, we leave them as is.
hdr = saveAnalyze(applyCannonicalXform(maxVec(:,:,:,1), b0.cannonical, b0.mmPerVox), fullfile(outPathName,'vecX'), mmPerVoxNew, b0.notes, b0.origin);
hdr = saveAnalyze(applyCannonicalXform(maxVec(:,:,:,2), b0.cannonical, b0.mmPerVox), fullfile(outPathName,'vecY'), mmPerVoxNew, b0.notes, b0.origin);
hdr = saveAnalyze(applyCannonicalXform(maxVec(:,:,:,3), b0.cannonical, b0.mmPerVox), fullfile(outPathName,'vecZ'), mmPerVoxNew, b0.notes, b0.origin);

% Compute and show test images
% t1Hdr should contain the Talairach xform that we just saved.
[t1Img, t1Mm, t1Hdr] = loadAnalyze(t1.analyzeFname);
%t1.img = makeCubeIfiles(t1.baseFname, t1.imDim(1:2), [1:t1.imDim(3)]);

imDims = [-70,70;-120,80;-60,80]';
slice = [0,0,0];

done = false;
newM = b0.M_tal;
while(~done & useManual)
    M = t1Hdr.mat;
    img = t1Img;
    fig = 90;
    [Xsl,x,y,z] = dtiGetSlice(M, img, 3, slice(3), imDims);
    figure(fig); h(1) = subplot(2,3,1); imagesc(imDims(:,1), imDims(:,2), Xsl); colormap(gray); axis equal tight xy;
    [Ysl,x,y,z] = dtiGetSlice(M, img, 2, slice(2), imDims);
    figure(fig); h(2) = subplot(2,3,2); imagesc(imDims(:,3), imDims(:,1), Ysl); colormap(gray); axis equal tight xy;
    [Zsl,x,y,z] = dtiGetSlice(M, img, 1, slice(1), imDims);
    figure(fig); h(3) = subplot(2,3,3); imagesc(imDims(:,3), imDims(:,2), Zsl); colormap(gray); axis equal tight xy;

    t1Xslice = Xsl; t1Yslice = Ysl; t1Zslice = Zsl;

    slice = [0,0,0];
    M = newM;
    %This XY swap business is why we have to do this permute here.
    img = b0.imgCannonical;
    [Xsl,x,y,z] = dtiGetSlice(M, img, 3, slice(3), imDims);
    figure(fig); h(4) = subplot(2,3,4); imagesc(imDims(:,1), imDims(:,2), Xsl); colormap(gray); axis equal tight xy;
    [Ysl,x,y,z] = dtiGetSlice(M, img, 2, slice(2), imDims);
    figure(fig); h(5) = subplot(2,3,5); imagesc(imDims(:,3), imDims(:,1), Ysl); colormap(gray); axis equal tight xy;
    [Zsl,x,y,z] = dtiGetSlice(M, img, 1, slice(1), imDims);
    figure(fig); h(6) = subplot(2,3,6); imagesc(imDims(:,3), imDims(:,2), Zsl); colormap(gray); axis equal tight xy;

    b0Xslice = Xsl; b0Yslice = Ysl; b0Zslice = Zsl;

    disp('select corresponding points on each B0 and T1 image.');
    donePoints = false;
    ptList = [];
    while(~donePoints)
        for(ii=1:6)
            axes(h(ii));
            title('current');
            p = ginput(1);
            title('');
            if(isempty(p))
                donePoints = true;
                break;
            end
            hold on; plot(p(:,1), p(:,2), 'xy'); hold off;
            pts(ii,:) = p;
        end
        % Only add this set if it is complete
        if(size(pts,1)==6)
            if(isempty(ptList)) ptList = pts;
            else ptList = [ptList; pts]; end
        end
    end

    % Adjust newM given points...
    adjustment = eye(4);

    fig = 92;
    t1mv = max([t1Xslice(:); t1Yslice(:); t1Zslice(:)]);
    b0mv = max([b0Xslice(:); b0Yslice(:); b0Zslice(:)]);
    Xsl(:,:,1) = t1Xslice./t1mv.*.5 + b0Xslice./b0mv.*.5;
    Xsl(:,:,2) = t1Xslice./t1mv.*.5; Xsl(:,:,3) = t1Xslice./t1mv.*.5;
    Ysl(:,:,1) = t1Yslice./t1mv.*.5 + b0Yslice./b0mv.*.5;
    Ysl(:,:,2) = t1Yslice./t1mv.*.5; Ysl(:,:,3) = t1Yslice./t1mv.*.5;
    Zsl(:,:,1) = t1Zslice./t1mv.*.5 + b0Zslice./b0mv.*.5;
    Zsl(:,:,2) = t1Zslice./t1mv.*.5; Zsl(:,:,3) = t1Zslice./t1mv.*.5;
    figure(fig); subplot(2,2,1); imagesc(imDims(:,1), imDims(:,2), Xsl); colormap gray; axis equal tight xy;
    figure(fig); subplot(2,2,2); imagesc(imDims(:,3), imDims(:,1), Ysl); colormap(gray); axis equal tight xy;
    figure(fig); subplot(2,2,3); imagesc(imDims(:,3), imDims(:,2), Zsl); colormap(gray); axis equal tight xy;

    bn = questdlg('Is this good?', ...
                  'Continue manual registration...', ...
                  'Yes','No','Abort','No');
    if(strcmpi(bn, 'Yes')) done = true; end
    if(strcmpi(bn, 'Abort')) done = true; useManual = false; end
end

% Save all the transform files
if(useManual)
    M = adjustment*b0.M;
    M_tal = adjustment*b0.M_tal;
else
    slice = [0,0,0];
    img = t1Img;
    [Xsl,x,y,z] = dtiGetSlice(t1Hdr.mat, img, 3, slice(3), imDims);
    [Ysl,x,y,z] = dtiGetSlice(t1Hdr.mat, img, 2, slice(2), imDims);
    [Zsl,x,y,z] = dtiGetSlice(t1Hdr.mat, img, 1, slice(1), imDims);
    t1Xslice = Xsl; t1Yslice = Ysl; t1Zslice = Zsl;

    img = b0.imgCannonical;
    [Xsl,x,y,z] = dtiGetSlice(b0.M_tal, img, 3, slice(3), imDims);
    [Ysl,x,y,z] = dtiGetSlice(b0.M_tal, img, 2, slice(2), imDims);
    [Zsl,x,y,z] = dtiGetSlice(b0.M_tal, img, 1, slice(1), imDims);
    b0Xslice = Xsl; b0Yslice = Ysl; b0Zslice = Zsl;

    fig = figure;
    t1mv = max([t1Xslice(:); t1Yslice(:); t1Zslice(:)]);
    b0mv = max([b0Xslice(:); b0Yslice(:); b0Zslice(:)]);
    Xsl(:,:,1) = t1Xslice./t1mv.*.5 + b0Xslice./b0mv.*.5;
    Xsl(:,:,2) = t1Xslice./t1mv.*.5; Xsl(:,:,3) = t1Xslice./t1mv.*.5;
    Ysl(:,:,1) = t1Yslice./t1mv.*.5 + b0Yslice./b0mv.*.5;
    Ysl(:,:,2) = t1Yslice./t1mv.*.5; Ysl(:,:,3) = t1Yslice./t1mv.*.5;
    Zsl(:,:,1) = t1Zslice./t1mv.*.5 + b0Zslice./b0mv.*.5;
    Zsl(:,:,2) = t1Zslice./t1mv.*.5; Zsl(:,:,3) = t1Zslice./t1mv.*.5;
    figure(fig); subplot(2,2,1); imagesc(imDims(:,1), imDims(:,2), Xsl); colormap gray; axis equal tight xy;
    figure(fig); subplot(2,2,2); imagesc(imDims(:,3), imDims(:,1), Ysl); colormap(gray); axis equal tight xy;
    figure(fig); subplot(2,2,3); imagesc(imDims(:,3), imDims(:,2), Zsl); colormap(gray); axis equal tight xy;

    M = b0.M;
    M_tal = b0.M_tal;
end
save(fullfile(outPathName,'B0'), 'M', 'M_tal');
save(fullfile(outPathName,'FA'), 'M', 'M_tal');
save(fullfile(outPathName,'vecX'), 'M', 'M_tal');
save(fullfile(outPathName,'vecY'), 'M', 'M_tal');
save(fullfile(outPathName,'vecZ'), 'M', 'M_tal');

% Now build and save the dt6 data
[eigVec, eigVal] = dtiLoadTensor(fullfile(fileparts(b0.baseFname), 'Vectors.float.'));
dt6_tmp = dtiRebuildTensor(eigVec, eigVal);
dt6_tmp = permute(dt6_tmp,[2 1 3 4]);
for(ii=1:6)
    dt6(:,:,:,ii) = applyCannonicalXform(dt6_tmp(:,:,:,ii), b0.cannonical, b0.mmPerVox);
end
clear dt6_tmp;
mmPerVox = mmPerVoxNew;
notes = b0.notes;
xformToAnat = M;
xformToTal = M_tal;
b0 = int16(round(b0.imgCannonical));
anat.img = int16(t1Img);
anat.mmPerVox = t1Mm;
anat.xformToTal = t1Hdr.mat;
save(fullfile(outPathName,'dt6'), 'dt6', 'mmPerVox', 'notes', 'xformToAnat', 'xformToTal', 'b0', 'anat');

return;

[subjects, baseDir] = getSubjects('dc', 'c', 'dtiSubjectsCortexPaper.txt');
%baseDir = pwd;
d = dir(baseDir);
subjects = {};
for(ii=1:length(d))
    if(~strcmp(d(ii).name,'.') & ~strcmp(d(ii).name,'..') ...
        & (strcmp(d(ii).name(3),'0') | strcmp(d(ii).name(4),'0')))
        subjects{end+1} = d(ii).name;
    end
end
for(ii=1:length(subjects))
    sub = subjects{ii};
    base = fullfile(baseDir,sub);
    try
    dtiMakeAnalyze(fullfile(base,'dti','B0.001'),...
        fullfile(base,'t1','spgr1','I.001'), ...
        fullfile(base,'t1',[sub '_t1anat_avg']));
    catch
        disp(['skipped subject ' sub]);
    end
end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [train_users, dev_users] = get_balanced_fold(DISFA_dir, users, au, prop_test, offset)

    % This is for loading the labels
    for i=1:numel(users)
        input_train_label_files{i} = [DISFA_dir, '/ActionUnit_Labels/', users{i}, '/', users{i}];
    end

    % Extracting the labels
    labels_train = extract_au_labels(input_train_label_files, au);

    counts = zeros(numel(users),1);
    for k=1:numel(users)
        counts(k) = sum(labels_train((k-1)*4844+1:k*4844));
    end

    [sorted, inds] = sort(counts);

    dev_users = users(inds(offset:round(1/prop_test):end));
    train_users = setdiff(users, dev_users);

    count_dev = 0;
    count_train = 0;
    for k=1:numel(users)
        if(any(strcmp(dev_users, users{k})))
            count_dev = count_dev + counts(k);
        else
            count_train = count_train + counts(k);
        end

    end
end

author: TadasBaltrusaitis

repo: OpenFace


text:
function [ node_num, face_num, normal_num, order_max ] = obj_size ( ...
  input_file_name )

%*****************************************************************************80
%
%% OBJ_SIZE determines sizes of graphics objects in an Alias OBJ file.
%
%  Discussion:
%
%    The only items of interest to this routine are vertices,
%    faces, and normal vectors.
%
%  Example:
%
%    #  magnolia.obj
%
%    v -3.269770 -39.572201 0.876128
%    v -3.263720 -39.507999 2.160890
%    ...
%    v 0.000000 -9.988540 0.000000
%
%    vn 1.0 0.0 0.0
%    ...
%    vn 0.0 1.0 0.0
%
%    f 8 9 11 10
%    f 12 13 15 14
%    ...
%    f 788 806 774
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 September 2008
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string INPUT_FILE_NAME, the input file name.
%
%    Output, integer NODE_NUM, the number of points.
%
%    Output, integer FACE_NUM, the number of faces.
%
%    Output, integer NORMAL_NUM, the number of normal vectors.
%
%    Output, integer ORDER_MAX, the maximum face order.
%
  ierror = 0;

  face_num = 0;
  node_num = 0;
  normal_num = 0;
  order_max = 0;
  text_num = 0;
%
%  If no file input, try to get one from the user.
%
  if ( nargin < 1 )
    input_file_name = input ( 'Enter the name of the ASCII OBJ file.' );
    if ( isempty ( input_file_name ) )
      return
    end
  end
%
%  Open the file.
%
  input_file_unit = fopen ( input_file_name, 'r' );

  if ( input_file_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'OBJ_SIZE - Fatal error!\n' );
    fprintf ( 1, '  Could not open the file "%s".\n', input_file_name );
    error ( 'OBJ_SIZE - Fatal error!' );
    return
  end
%
%  Read a line of text from the file.
%
  while ( 1 )

    text = fgetl ( input_file_unit );

    if ( text == -1 )
      break
    end

    text_num = text_num + 1;
%
%  Replace any control characters (in particular, TABs) by blanks.
%
    s_control_blank ( text );

    done = 1;
    word_index = 0;
%
%  Read a word from the line.
%
    [ word, done ] = word_next_read ( text, done );
%
%  If no more words in this line, read a new line.
%
    if ( done )
      continue
    end
%
%  If this word begins with '#' or '$', then it is a comment.  Read a new line.
%
    if ( word(1) == '#' || word(1) == '$' )
      continue
    end

    word_index = word_index + 1;

    if ( word_index == 1 )
      word_one = word;
    end
%
%  F V1 V2 V3 ...
%    or
%  F V1/VT1/VN1 V2/VT2/VN2 ...
%    or
%  F V1//VN1 V2//VN2 ...
%
%  Face.
%  A face is defined by the vertices.
%  Optionally, slashes may be used to include the texture vertex
%  and vertex normal indices.
%
    if ( s_eqi ( word_one, 'F' ) )

      face_num = face_num + 1;

      vertex = 0;

      while ( 1 )

        [ word, done ] = word_next_read ( text, done );

        if ( done )
          break
        end

        vertex = vertex + 1;
        order_max = max ( order_max, vertex );
%
%  Locate the slash characters in the word, if any.
%
        i1 = ch_index ( word, '/' );
        if ( 0 < i1 )
          i2 = ch_index ( word(i1+1), '/' ) + i1;
        else
          i2 = 0;
        end
%
%  Read the vertex index.
%
        itemp = s_to_i4 ( word );
%
%  If there are two slashes, then read the data following the second one.
%
        if ( 0 < i2 )
          itemp = s_to_i4 ( word(i2+1) );
        end

      end
%
%  V X Y Z W
%  Geometric vertex.
%
    elseif ( s_eqi ( word_one, 'V' ) )

      node_num = node_num + 1;
      continue
%
%  VN
%  Vertex normals.
%
    elseif ( s_eqi ( word_one, 'VN' ) )

      normal_num = normal_num + 1;
      continue

    end

  end

  fclose ( input_file_unit );

  return
end


text:
% FAST_NONMAX perform non-maximal suppression on FAST features.
%
%     nonmax = FAST_NONMAX(image, threshold, FAST_CORNER_DETECT_9(image, threshold));
%     returns a list of nonmaximally suppressed corners with the X coordinate
%     in nonmax(:,1) and Y in nonmax(:,2).
%
%     If you use this in published work, please cite:
%       Fusing Points and Lines for High Performance Tracking, E. Rosten and T. Drummond, ICCV 2005
%       Machine learning for high-speed corner detection, E. Rosten and T. Drummond, ECCV 2006
%     The Bibtex entries are:
%
%     @inproceedings{rosten_2005_tracking,
%         title       =    "Fusing points and lines for high performance tracking.",
%         author      =    "Edward Rosten and Tom Drummond",
%         year        =    "2005",
%         month       =    "October",
%         pages       =    "1508--1511",
%         volume      =    "2",
%         booktitle   =    "IEEE International Conference on Computer Vision",
%         notes       =    "Oral presentation",
%         url         =    "http://mi.eng.cam.ac.uk/~er258/work/rosten_2005_tracking.pdf"
%     }
%
%     @inproceedings{rosten_2006_machine,
%         title       =    "Machine learning for high-speed corner detection",
%         author      =    "Edward Rosten and Tom Drummond",
%         year        =    "2006",
%         month       =    "May",
%         booktitle   =    "European Conference on Computer Vision (to appear)",
%         notes       =    "Poster presentation",
%         url         =    "http://mi.eng.cam.ac.uk/~er258/work/rosten_2006_machine.pdf"
%     }
%
%
%
%     See also FAST_CORNER_DETECT_9, FAST_CORNER_DETECT_10
%              FAST_CORNER_DETECT_11, FAST_CORNER_DETECT_12
function ret = fast_nonmax(im, barrier, c)
	sz = size(im);
	ysize = sz(1);

	% x, y
	dir=[	0,3
			1,3
			2,2
			3,1
			3,0
			3,-1
			2,-2
			1,-3
			0,-3
			-1,-3
			-2,-2
			-3,-1
			-3,0
			-3,1
			-2,2
			-1,3];

	dir = dir(:,2) + dir(:,1) * ysize;
	centres = c(:,2) + (c(:,1)-1) * ysize;

	%First pass: compute scores for all of the pixels:
	scores = zeros(sz(1) * sz(2), 1);

	for i=1:length(centres)

		p1 = im(centres(i)+dir) -( im(centres(i)) + barrier);
		pos = sum(p1 .* (p1 > 0));
		n1 = im(centres(i)) - barrier - im(centres(i)+dir);
		neg = sum(n1 .* (n1 > 0));

		scores(centres(i)) = max([pos neg]);
	end

	cs = zeros(length(centres),1);

	up = -1;
	down = 1;
	left =-ysize;
	right = ysize;

	%second pass: get local maxima
	for i=1:length(centres)
		p = centres(i);
		square=p +  [ up down left right up+left up+right down+left down+right];

		cs(i) = (sum(scores(p) >= scores(square)) == 8);
	end

	%Get the maxima positions
	maximas = centres(find(cs));

	ret = zeros(length(maximas), 2);

	ret(:,1) = 1 + floor(maximas / ysize);
	ret(:,2) = mod(maximas, ysize);

repo: ekfmonoslam


text:
% Identical to isomapEmbed without the last part which sometimes causes exceptions

function [X, sigma2] = isomap2Embed(Y, dims)

% ISOMAP2EMBED Embed data set with Isomap.

% SHEFFIELDML

% Note: isomap code uses the transpose of a design matrix.
if any(any(isnan(Y)))
  error('Cannot initialise gplvm using isomap when missing data is present.')
else
  D = L2_distance(Y', Y', 1);
  options.dims = 1:dims;
  neighbours = 7;
   options.display = 0;

  [Xstruct, sigma2, E] = Isomap(D, 'k', neighbours, options);
  X = zeros(size(Y, 1), dims);
  if length(Xstruct.index) ~= size(Y, 1)
    % We don't really deal with this problem correctly here ...
    warning('Isomap graph is not fully connected');
  end
  X(Xstruct.index, :) = Xstruct.coords{dims}';
  % Rescale X so that variance is 1 and mean is zero.
  meanX = mean(X);
  X = X-ones(size(Y, 1), 1)*meanX;
  varX = var(X);
  X = X*diag(sqrt(1./varX));
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function corner_4_coordinates = Get4Points( corner_coordinates )
%GET4POINTS Summary of this function goes here
%
variance = 10;   % delta_x + delta_y: The distance of the 4 corners should not <= the designated variance
numberOfPoints = 0; % should be 4 corners eventually for a quadrilateral

%corner_coordinates = sort(corner_coordinates,2);
sizeOfCornerPoints = size(corner_coordinates); % 1 element of the vector is the number of row

delta_x = variance; delta_y = 1;    % ensure (delta_x + delta_y) > variance to start with

% for the 4 corners
numberOfPoint1 = 0; numberOfPoint2 = 0; numberOfPoint3 = 0; numberOfPoint4 = 0;

% get the 4 designated points
        for i=1:sizeOfCornerPoints(1)
            if (i+1) < sizeOfCornerPoints(1)
                % appoint it as the candidate point
                if (numberOfPoints < 4) & ((delta_x + delta_y) > variance )
                    numberOfPoints = numberOfPoints+1;
                    designatedPoint(numberOfPoints,:) = corner_coordinates(i,:);
                end

                delta_x = abs(designatedPoint(numberOfPoints,1) - corner_coordinates(i+1,1));    % difference of the 2 x values
                delta_y = abs(designatedPoint(numberOfPoints,2) - corner_coordinates(i+1,2));    % difference of the 2 y values

            end
        end

        corner_4_coordinates = designatedPoint;
end

% Notes
% =====
%  Recommendation: The points may be averaged as well
%

author: Sable

repo: mcbench-benchmarks


text:
function y = filtdn(x, f, dim, extmod, shift)
% FILTDN   Filter and downsample (by 2) along a dimension
%
%       y = filtdn(x, f, dim, extmod, shift)
%
% Input:
%   x:      input signal
%   f:      1-D filter
%   dim:    the processing dimension
%   extmod: extension mode (e.g. 'per' or 'sym')
%   shift:  specifies the window over which filtering occurs
%
% Output:
%   y:      filtered and dowsampled signal
%
% Note:
%   The origin of the filter f is assumed to be floor(size(f)/2) + 1.
%   Amount of shift should be no more than floor((size(f)-1)/2).

% Skip singleton dimension
if size(x, dim) == 1
    y = x;
    return
end

% Cell array of indexes for each dimension
nd = ndims(x);
I = cell(1, nd);
for d = 1:nd
    I{d} = 1:size(x,d);
end

% Border extend
n = size(x, dim);
hlf = (length(f) - 1) / 2;
% Amount of extension at two ends
e1 = floor(hlf) + shift;
e2 = ceil(hlf) - shift;

switch extmod
    case 'per'
        I{dim} = [ly-e1+1:n , 1:n , 1:e2];

    case 'sym'
        I{dim} = [e1+1:-1:2 , 1:n , n-1:-1:e2];

    otherwise
        error('Invalid input for EXTMOD')

end
y = x(I{:});

% Filter, downsample, and return only the 'valid' part
y = filter(f, 1, y, [], dim);

I{dim} = (1:2:n) + length(f) - 1;
y = y(I{:});


text:
function dImgLow = fDownSample( dImgHigh, out_size )

out_x_sz = out_size(1); out_y_sz = out_size(2); out_z_sz = out_size(3);
[in_x_sz,in_y_sz, in_z_sz, nCha, nTime] = size( dImgHigh );

% build grid vectors for the up/down sampling
% ============================================
% if the input is even & output is odd-> use floor for all
% if the output is even & input is odd -> use ceil for all
% other cases - don't care
% for downsampling -> the opposite
if (~mod( in_x_sz,2 ) & (out_x_sz>in_x_sz)) | (mod( in_x_sz,2 ) & (out_x_sz<in_x_sz))
    x_output_space  = max(floor((out_x_sz-in_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
    x_input_space   = max(floor((in_x_sz-out_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
else
    x_output_space  = max(ceil((out_x_sz-in_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
    x_input_space   = max(ceil((in_x_sz-out_x_sz)/2),0) + [1:min(in_x_sz,out_x_sz)];
end
if (~mod( in_y_sz,2 ) & (out_y_sz>in_y_sz)) | (mod( in_y_sz,2 ) & (out_y_sz<in_y_sz))
   y_output_space  = max(floor((out_y_sz-in_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
   y_input_space   = max(floor((in_y_sz-out_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
else
   y_output_space  = max(ceil((out_y_sz-in_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
   y_input_space   = max(ceil((in_y_sz-out_y_sz)/2),0) + [1:min(in_y_sz,out_y_sz)];
end
if (~mod( in_z_sz,2 ) & (out_z_sz>in_z_sz)) | (mod( in_z_sz,2 ) & (out_z_sz<in_z_sz))
   z_output_space  = max(floor((out_z_sz-in_z_sz)/2),0) + [1:min(in_z_sz,out_z_sz)];
   z_input_space   = max(floor((in_z_sz-out_z_sz)/2),0) + [1:min(in_z_sz,out_z_sz)];
else
   z_output_space  = max(ceil((out_z_sz-in_z_sz)/2),0) + [1:min(in_z_sz,out_z_sz)];
   z_input_space   = max(ceil((in_z_sz-out_z_sz)/2),0) + [1:min(in_z_sz,out_z_sz)];
end

dImgLow = zeros( out_x_sz, out_y_sz, out_z_sz, size(dImgHigh,4), size(dImgHigh,5) );
dImgLow(x_output_space,y_output_space, z_output_space, :, : ) = dImgHigh(x_input_space,y_input_space, z_input_space, :, :);

end

author: thomaskuestner


text:
%%% void show_results(matrix,matrix,matrix,matrix,int,int,int)
% I - Input sequence
% L - Low-rank sequence
% S - Sparse sequence
% O - Outliers sequence
% nFrame - Number of frames
% vidHeight - Video height
% vidWidth - Video width
%
function show_results(I,L,S,O,nFrame,vidHeight,vidWidth)
  warning('off','all');
  clf;
  for i = 1:nFrame
    Input = reshape(I(:,i),vidHeight,vidWidth);
    %Input = uint8(Input);

    if(size(L) == size(I))
      LowRank = reshape(L(:,i),vidHeight,vidWidth);
      %LowRank1 = reshape(L(:,i),vidHeight,vidWidth);
      %LowRank = mat2gray(LowRank1);
      %LowRank = im2uint8(LowRank);
    else
      LowRank = uint8(zeros(size(Input)));
    end

    if(size(S) == size(I))
      Sparse = reshape(S(:,i),vidHeight,vidWidth);
      %Sparse1 = reshape(S(:,i),vidHeight,vidWidth);
      %Sparse = mat2gray(Sparse1);
      %Sparse = im2uint8(Sparse);
    else
      Sparse = uint8(zeros(size(Input)));
    end

    if(~isempty(O))
      Outlier = reshape(O(:,i),vidHeight,vidWidth);
      %Outlier1 = reshape(O(:,i),vidHeight,vidWidth);
      %Outlier = mat2gray(Outlier1);
      %Outlier = im2uint8(Outlier);
      %Outlier = medfilt2(Outlier, [5 5]);
    else
      Outlier = uint8(zeros(size(Input)));
    end

    subplot(2,3,1), imshow(Input,[]), title('Input (I)');
    subplot(2,3,2), imshow(LowRank,[]), title('Low Rank (L)');
    subplot(2,3,3), imshow(Sparse,[]), title('Sparse (S)');
    subplot(2,3,4), imshow(Outlier,[]), title('Outliers (O)');
    subplot(2,3,5), imshow(medfilt2(Outlier, [5 5]),[]), title('Filtered Outliers');
    %disp(i);
    pause(0.01);
  end
  warning('on','all');
end

author: andrewssobral

repo: lrslibrary


text:
function [dop_buf] = rayleigh_dop(in_val)
% Function to generate Doppler-filtered Rayleigh fading simulator

% Doppler parameters
vo = (1.2*1000/3600);    % vo = 1.2 km/h  ->  0.333 m/s
lambda = 3e8 / 5.4e9;    % lambda = c / fc;
fd = vo/lambda;          % Doppler Spread (around 6 Hz)

% Generate Doppler Sf (256 points)
f_rng = 10*fd;
f = -f_rng/2 : f_rng/255 : f_rng/2;
dop_Sf = 1./(1+9*(f/fd).^2);

% Generate 256 samples of Doppler-filtered, Rayleigh fading sim.
dop_buf1(129:256) = randn(1,128)+j*randn(1,128);
dop_buf1(  1:128) = conj(dop_buf1(256:-1:129));
dop_buf2(129:256) = randn(1,128)+j*randn(1,128);
dop_buf2(  1:128) = conj(dop_buf2(256:-1:129));
dop_buf1 = ifft(dop_buf1 .* dop_Sf);
dop_buf2 = ifft(dop_buf2 .* dop_Sf);
dop_buf = sqrt(dop_buf1.^2 + dop_buf2.^2);


text:
% Lists files only (not directories).
system('ls -l | egrep -v ''^d''');

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function handle = vwhist(x,y,dx,dy,direction,FaceColor,EdgeColor)
%VWHIST   Variable width histogram
%   Draws a variable width histogram.
%   Does not calculate data, it is just a drawing utility.
%   Draws histogram with defined x, so data length(data)=length(x)-1,
%   or, if lengths are equal, bars are centred at x.
%
%   Syntax:
%      HANDLES = VWHIST(X,Y,DX,DY,DIREC,FaceColor,EdgeColor)
%
%   Inputs:
%      X           X locations
%      Y           Data
%      DX, DY      Offset added to bars [ 0 ]
%      DIREC       Direction, horizontal or vertical [ {0} | 1 ]
%      FaceColor   Bars FaceColor [ 'y' ]
%      EdgeColor   Bars EdgeColor [ 'r' ]
%
%   Output:
%      HANDLES   Handles to fill objects (bars)
%
%   Example:
%      x= [1.5 3 4 6 7 7.5];
%      y =[  3  7  2  0 3  ];
%      figure,
%
%      subplot(2,2,1)
%      handle = vwhist(x,y);
%      axis([0 10 0 10]);
%
%      subplot(2,2,2)
%      handle = vwhist(x,y,0,2,1); set(handle,'facecolor','g');
%      axis([0 10 0 10]);
%
%      y = [4 y];
%      subplot(2,2,3)
%      handle = vwhist(x,y,0,2); set(handle,'facecolor','w')
%      axis([0 10 0 10]);
%
%
%   MMA 8-10-2004, martinho@fis.ua.pt

%   Department of physics
%   University of Aveiro, Portugal

handle=[];

if nargin < 7
  EdgeColor = 'r';
end
if nargin < 6
  FaceColor = 'y';
end
if nargin < 5
  direction = 0;
end
if nargin < 4
  dy = 0;
end
if nargin < 3
  dx = 0;
end
if nargin < 2
  y = x;
  x=1:length(y);
end
if nargin < 1
  disp('# nothing to do...');
end

% check lengths:
nx = length(x);
ny = length(y);

if nx == ny
  x_ = (x(2:end)+x(1:end-1))/2;
  dx1 = x_(1) - x(1);
  dxe = x(end) - x_(end);
  x = [x(1)-dx1 x_ x(end)+dxe];
elseif nx ~= ny+1
  disp('# wrong lengths of x and y, nx=ny or nx=ny+1, only')
end

for i=1:length(x)-1
  xx = [x(i) x(i+1) x(i+1) x(i) x(i) ] + dx;
  yy = [0    0      y(i)   y(i) 0    ] + dy;

  if ~direction
    handle(i) = fill(xx,yy,FaceColor);
  else
    handle(i) = fill(yy,xx,FaceColor);
  end
  set(handle(i),'EdgeColor',EdgeColor);
  if i==1
    h = ishold;
    hold on
  end
end

if ~h
  hold off
end


text:
% creates cell masks for fluorescence computation
% also creates exclusion mask for neuropil (cellPix)
function [stat, centerMasks, surroundMasks, neuropMasks] = ...
    createCenterSurroundMasks2(ops, stat, Ny, Nx, allow_overlap, radius0)

if nargin > 4
    aov = allow_overlap;
else
    aov = 0;
end

if nargin > 4
    radius = radius0;
else
    radius = 3;
end

Nk = length(stat);
cellPix = zeros(Ny, Nx);
centerMasks = zeros(Nk, Ny, Nx, 'single');
surroundMasks = zeros(Nk, Ny, Nx, 'single');

for k = 1:Nk
    % only use non-overlapping pixels of cell
    temp = zeros(Ny, Nx);
    if aov
        ipix = stat(k).ipix(stat(k).isoverlap==0);
        ypix = stat(k).ypix(stat(k).isoverlap==0);
        xpix = stat(k).xpix(stat(k).isoverlap==0);
        lam  = stat(k).lam(stat(k).isoverlap==0);
    else
        ipix = stat(k).ipix;
        ypix = stat(k).ypix;
        xpix = stat(k).xpix;
        lam  = stat(k).lam;
    end

    % fit MV gaussian to cell mask
    % define cell as all pixels within 2 std's of lambda
    if ~isempty(ypix)
        params      = FitMVGaus(ypix, xpix, lam, 2);
        % cell radius
        stat(k).radius = sqrt(mean(params.eval));

        % add masks for cells with radius greater than
        temp        = zeros(Ny, Nx);
        tempS       = zeros(Ny, Nx);
        if stat(k).radius >= radius
            pdist = ((stat(k).ypix - stat(k).med(1)).^2 + ...
                (stat(k).xpix - stat(k).med(2)).^2);
            temp(stat(k).ipix( pdist   < (stat(k).radius*.7)^2)) = 1;
            tempS(stat(k).ipix(pdist >= (stat(k).radius*.7)^2)) = 1;
        end

        centerMasks(k, :, :)      = temp / sum(temp(:));
        surroundMasks(k, :, :)    = tempS / sum(tempS(:));
    else
        stat(k).radius = 0;
    end

    % use all pixels for neuropil mask computation
    temp        = zeros(Ny, Nx);
    ipix = stat(k).ipix;
    lam  = stat(k).lam;
    temp(ipix)  = lam;

    % input thresholded pixels to cellPix to exclude in neuropil computation
    cellPix = cellPix + (temp > 0);
end

% create surround neuropil masks
[~, neuropMasks] = createNeuropilMasks(ops, stat, cellPix);

author: cortex-lab


text:
%--- help for set_inputs ---
%
%  INTERNAL FUNCTION: Sets inputs for the VAR
%
%  ::
%
%     var_obj = set_inputs(var_obj,varargin);
%
%  Args:
%     var_obj (var object): VAR object
%     varargin: Options to set. Must come in pairs
%
%        - **'data'**: Data for the VAR. Can be either a struct or ts object
%        - **'prior'**: Prior to use if BVAR. Available priors are
%
%           - minnesota
%           - jeffrey
%           - nw
%           - normal-wishart
%           - inw
%           - sz
%
%        - **'linear_restrictions'**: Linear restrictions. For the format, see XXXXXXXXX
%
%  Returns:
%     :
%
%        - var_obj (var object): New VAR with the given inputs updated.
%
%

author: jmaih


text:
function v=mai(x,lag)
%Syntax: v=mai(x,lag)
%____________________
%
% Calculates the mutual average information of a time series x for
% some time lag.
%
% v is the the value of the mutual average information.
% x is the time series.
% lag is the time lag.
%
% Alexandros Leontitsis
% Institute of Mathematics and Statistics
% University of Kent at Canterbury
% Canterbury
% Kent, CT2 7NF
% U.K.
% University e-mail: al10@ukc.ac.uk
% Lifetime e-mail: leoaleq@yahoo.com
% Homepage: http://www.geocities.com/CapeCanaveral/Lab/1421
%
% May 25, 2001.

if nargin<1 | isempty(x)==1
   error('You should provide a time series.');
else
   % x must be a vector
   if min(size(x))>1
      error('Invalid time series.');
   end
   x=x(:);
   % n is the time series length
   n=length(x);
end

if nargin<2 | isempty(lag)==1
   lag=0:min(n/2-1,20);
else
   % lag must be a vector
   if min(size(lag))>1
      error('The time lag must be a scalar or a vector.');
   end
   % lag must contain integers
   lag=round(lag);
   % lag values must be between 0 and n/2-1
   lag=lag(find(lag>=0 & lag<n/2));
   % lag must not be empty
   if isempty(lag)==1
      error('You must give another set of values for lag.')
   end
end

% The mutual average information
x=x-min(x);
x=x/max(x);
for i=1:length(lag)

   % Define the number of bins
   k=floor(1+log2(n-lag(i))+0.5);

   % If the time series has no variance then the MAI is 0
   if var(x,1)==0
      v(i)=0;
   else
      v(i)=0;
      for k1=1:k
         for k2=1:k
            ppp=find((k1-1)/k<x(1:n-lag(i)) & x(1:n-lag(i))<=k1/k ...
               & (k2-1)/k<x(1+lag(i):n) & x(1+lag(i):n)<=k2/k);
            ppp=length(ppp);
            px1=find((k1-1)/k<x(1:n-lag(i)) & x(1:n-lag(i))<=k1/k);
            px2=find((k2-1)/k<x(1+lag(i):n) & x(1+lag(i):n)<=k2/k);
            if ppp>0
               ppp=ppp/(n-lag(i));
               px1=length(px1)/(n-lag(i));
               px2=length(px2)/(n-lag(i));
               v(i)=v(i)+ppp*log2(ppp/px1/px2);
            end
         end
      end
   end
end


text:
function compareFigures(fig1, fig2)
%COMPAREFIGURES - places the second figure next to the first one
%
% Syntax:
%   compareFigures(fig1, fig2) places the contents of fig2 next to the
%   contents of fig1, inside fig2

% Copyright 2008 The MathWorks, Inc.

pos1 = get(fig1,'Position');
pos2 = get(fig2,'Position');

%set absolute units to prevent resizing
setChildrenToPixel(fig1);
setChildrenToPixel(fig2);

%move the item from the second figure to the first
moveChildren(fig2,fig1,pos1(2));
%get rid of the second figure
close(fig2);

%resize the first figure to accodate the second
pos1(3) = pos1(3) + pos2(3);
pos1(4) = max(pos1(4),pos2(4));
set(fig1, 'Position', pos1);

%center on screen
set(fig1,'units','normalized');
pos1 = get(fig1,'Position');
pos1(1) = (1 - pos1(3)) / 2;
pos1(2) = (1 - pos1(4)) / 2;
set(fig1,'Position',pos1);
%bring to the front
figure(fig1);

function setChildrenToPixel(fig)
children = get(fig,'Children');
for i=1:length(children)
    set(children(i),'Units','pixels')
end

function moveChildren(from,to,x)
children = get(from,'Children');
for i=1:length(children)
    pos = get(children(i),'Position');
    pos(1) = x;
    set(children(i),'Parent',to,'Position', pos);
end


text:
function obj_cell = split(obj, varargin)
% Split fmri_data object into separate objects with different images in each
% default: Use the obj.images_per_session field
%
% :Usage:
% ::
%
%     obj_cell = split(obj, [optional inputs])
%
% For objects: Type methods(object_name) for a list of special commands
%              Type help object_name.method_name for help on specific
%              methods.
%
% ..
%     Author and copyright information:
%
%     Copyright (C)2016 Tor Wager
%
%     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/>.
% ..
%
% :Inputs:
%
%   **obj:**
%        An fmri_display object with multiple images, to split
%
% :Optional Inputs:
%   **'imgs_per_obj'**
%        Followed by a vector of number of images in each object
%        e.g., for 100 images, [25 25 25 25] would split into 4 objects of
%        25 images each
%        Default: Use obj.images_per_session
%
% :Outputs:
%
%   **objcell:**
%        A cell array with each output fmri_data object in a cell.
%
%
% :Examples:
% ::
%
%
% :References:
%   None.
%
% :See also:
%   - list other functions related to this one, and alternatives*
%

% ..
%    Programmers' notes:
%    List dates and changes here, and author of changes
% ..
% Tor Wager, Aug 2016
% Added  >= en(i) to some lines - bug fix. Seemed to have little effect on most cases. Tor, June 2017

% ..
%    DEFAULTS AND INPUTS
% ..

imgs_per_obj = obj.images_per_session;

% optional inputs with default values
for i = 1:length(varargin)
    if ischar(varargin{i})
        switch varargin{i}

            case 'imgs_per_obj', imgs_per_obj = varargin{i+1}; varargin{i+1} = [];

            otherwise, warning(['Unknown input string option:' varargin{i}]);
        end
    end
end

if isempty(imgs_per_obj)
    obj_cell = {obj};
    return
end

obj = replace_empty(obj);  % to replace removed images

if ~isrow(imgs_per_obj), imgs_per_obj = imgs_per_obj'; end

en = cumsum(imgs_per_obj);
st = [0 en(1:end-1)] + 1;

k = length(imgs_per_obj);
obj_cell = cell(1, k);

% warn user about fields that aren't handled but appear like they should be
% split. Some of these (e.g. additional_info) are handled by
% fmri_data.cat(), so it may surprise the user that split doesn't invert
% the process.
%
% criteria for warning: if length of field is equal to number of splits
fnames = {'image_names', 'fullpath', 'files_exist', 'X', 'Y'};
supportedFields = {fnames{:}, 'dat', 'images_per_session', 'removed_images'};
allFields = fieldnames(obj);

for i = 1:length(allFields)
    if sum(strcmp(allFields{i},supportedFields)) == 0
        if isa(obj.(allFields{i}),'table')
            if height(obj.(allFields{i})) == k
                warning(['Object''s ''', allFields{i}, ''' field length matches number of splits, but ''', allFields{i}, ''' field will not be split'])
            end
        elseif length(obj.(allFields{i})) == k
            warning(['Object''s ''', allFields{i}, ''' field length matches number of splits, but ''', allFields{i}, ''' field will not be split'])
        end
    end

end

for i = 1:k

    obj_cell{i} = obj;

    wh = st(i):en(i);

    obj_cell{i}.dat = obj_cell{i}.dat(:, wh);

    obj_cell{i}.images_per_session = en(i) - st(i) + 1;

    for j=1:length(fnames)

        if ~isempty(obj_cell{i}.(fnames{j})) && size(obj_cell{i}.(fnames{j}), 1) >= en(i)

            obj_cell{i}.(fnames{j}) = obj_cell{i}.(fnames{j})(wh, :);

        end
    end

    obj_cell{i} = check_image_filenames(obj_cell{i}, 'noverbose');

    if ~isempty(obj_cell{i}.removed_images) && size(obj_cell{i}.removed_images, 1) >= en(i)

        obj_cell{i}.removed_images = obj_cell{i}.removed_images(wh, :);

    else

        obj_cell{i}.removed_images = 0; % add in case

    end

end

end % function

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function [x_pca, m] = PCAMyself(x)
% Feature extraction by PCA

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

    [COEFF, SCORE, latent]=pca(x);
    answer=0;s1=0;s2=sum(latent);count=1;
    while answer < 0.95
        s1=s1+latent(count);
        answer=s1/s2;
        count=count+1;
    end
    m=COEFF(:,1:count-1);
    x_pca=x*m;
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function varargout = gradient(f)
%GRADIENT   Gradient of a CHEBFUN3.
%   [FX, FY, FZ] = GRAD(F) returns the gradient of the CHEBFUN3 object F,
%   where
%   FX is the partial derivative of F in the first variable,
%   FY is the partial derivative of F in the second varaiable, and
%   FZ is the partial derivative of F in the third variable.
%
%   G = GRAD(F) returns a CHEBFUN3V object which represents
%            G = [FX; FY; FZ].
%
% See also CHEBFUN3/GRAD.

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

F = chebfun3v(diffx(f), diffy(f), diffz(f));

if ( nargout <= 1 )
    varargout = {F};
else
    varargout = F.components;
end

end


text:
classdef MimEditCyan < MimGuiPlugin
    % MimEditCyan. Gui Plugin for setting paint colour
    %
    %     You should not use this class within your own code. It is intended to
    %     be used by the gui of the TD MIM Toolkit.
    %
    %     MimEditCyan is a Gui Plugin for the MIM Toolkit.
    %
    %
    %     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.
    %

    properties
        ButtonText = 'Cyan'
        SelectedText = 'Cyan'
        ToolTip = 'Edit with cyan label'
        Category = 'Segmentation label'
        Visibility = 'Dataset'
        Mode = 'Edit'

        HidePluginInDisplay = false
        PTKVersion = '1'
        ButtonWidth = 5
        ButtonHeight = 1
        Icon = 'paint.png'
        IconColour = GemMarkerPoint.DefaultColours{4}
        Location = 34
    end

    methods (Static)
        function RunGuiPlugin(gui_app)
            gui_app.ImagePanel.PaintBrushColour = 4;
        end

        function enabled = IsEnabled(gui_app)
             enabled = gui_app.IsDatasetLoaded && gui_app.ImagePanel.OverlayImage.ImageExists && ...
                isequal(gui_app.ImagePanel.SelectedControl, 'Paint');
        end

        function is_selected = IsSelected(gui_app)
            is_selected = gui_app.ImagePanel.PaintBrushColour == 4;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function varargout = min(varargin)
% Converts first or second inputs to double if they are char
if nargin>=1 && ischar(varargin{1})
    varargin{1} = double(varargin{1});
end
if nargin>=2 && ischar(varargin{2})
    varargin{2} = double(varargin{2});
end
if nargout==0
    nargout = 1; % if called without outputs: produce 1 output
end
varargout = cell(1,nargout);
if ~any(cellfun(@(x) isa(x,'sym'), varargin))
    [varargout{:}] = builtin('min', varargin{:});
else
    [varargout{:}] = builtin('@sym/min', varargin{:});
end
end

author: lmendo

repo: MATL

save_path: github-repos/MATLAB/lmendo-MATL


text:
function att = Cbn2att(Cbn);
%-------------------------------------------------------
% Yudan Yi, May 26, 2005
%-------------------------------------------------------
% attitude: roll, pitch and heading
att(1) = atan2(Cbn(3,2),Cbn(3,3));
att(2) = asin(-Cbn(3,1));
att(3) = atan2(Cbn(2,1),Cbn(1,1));
att = att(:);
return;


text:
function [cal] = kJ2cal(kJ)
% Convert energy or work from kilojoules to calories.
% Chad A. Greene 2012
cal = kJ*238.84589663;


text:
function varargout = vl_nnmax(M, varargin)
% VL_NNMAX Element-wise maximum
%   Y = VL_NNMAX(M, X1, ..., XM) applies the elementwise max operator to
%   the given M input tensors, each of which should have the same input
%   dimension.
%
%   [DZDX1,..,DZDXM] = VL_NNMAX(M,X1, ...,XM, DZDY) computes the
%   derivatives of the block projected onto DZDY.
%
% Copyright (C) 2017 Samuel Albanie
% Licensed under The MIT License [see LICENSE.md for details]

  assert(M >= 2, 'at least two input tensors expected') ;
  ins = varargin(1:M) ; varargin(1:M) = [] ;
  [~, dzdy] = vl_argparsepos(struct(), varargin) ;

  shared = cat(5, ins{:}) ;
  if isempty(dzdy)
     varargout{1} = max(shared, [], 5) ;
  else
    % recompute forward max and insert derivatives at max locations
    [~,I] = max(shared, [], 5) ;
    dzdx = zeros(size(shared), 'like', shared) ;
    offsets = (I - 1) .* numel(dzdy{1}) ;
    idx = offsets(:) + (1:numel(dzdy{1}))' ;
    dzdx(idx) = dzdy{1} ;
    varargout = arrayfun(@(x) {dzdx(:,:,:,:,x)}, 1:M) ;
  end

repo: MFT


text:
function C = dot_product(A,B)

% Computes the dot product of the corresponding rows of the matrices A and B

C = sum(A.*B,2);

author: InverseTampere

repo: TreeQSM


text:
%% Filename and Format string
% This code loads in the contents of the data file with textscan a block at a
% time and converts the cellarray from textscan into a double array

filename='waferdata.csv';
FormatString=['%*f%*f%*f%*f' repmat('%u8',1,9) repmat('%*f',1,7)];

%% Parameters
Headers=3;
NumLines=1e6; % Total number of lines to read
BlockLines=10000; % Size of block
NumBlocks=NumLines/BlockLines; % Number of blocks

%% Open file
fid = fopen(filename); % Open file

%% Preallocate data
data=zeros(NumLines,9,'uint8'); % Pre-allocate space for data

%% Get Headers
cellchunk=textscan(fid,'%s',3,'delimiter','\n'); % Get first 3 lines

%% Read in blocks
for Block=1:NumBlocks
    data((Block-1)*BlockLines+(1:BlockLines),:) = ...
	cell2mat(textscan(fid,FormatString,BlockLines,'delimiter',','));
    disp(['Block ' num2str(Block) ]); % Display current block to show progress
end

%% Close file
fclose(fid);


text:
%Copyright 2013 The MathWorks, Inc
function [sSAnt, sSAntPlat] = splitSubArrays(sAnt)
faceSlope = 45; % Change face slope
sAnt.SubarrayNormal = [0 90 180 -90;90-faceSlope 90-faceSlope 90-faceSlope 90-faceSlope];

%% Build 4 separate arrays
epos = getElementPosition(sAnt);
aznorm = [0 90 180 -90];
for i=1:4,
    sSAnt{i} = phased.ConformalArray; %#ok<*AGROW>
    set(sSAnt{i},'ElementPosition',epos(:,(i-1)*16+(1:16)),...
                 'ElementNormal',repmat([aznorm(i);faceSlope],1,16));
    % Use isotropic antenna elements
    sSAnt{i}.Element.FrequencyRange = sAnt.Subarray.Element.FrequencyRange;
    sSAnt{i}.Element.BackBaffled = true;
end

%% Antenna Platform
sSAntPlat = phased.Platform([2.5e4; 0; 0]);


text:
function sgmga_write_tests ( )

%****************************************************************************80
%
%% SGMGA_WRITE_TESTS calls SGMGA_WRITE_TEST with various arguments.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 June 2010
%
%  Author:
%
%    John Burkardt
%
%  Local Parameters:
%
%    Local, real TOL, a tolerance for point equality.
%    A value of sqrt ( eps ) is reasonable, and will allow the code to
%    consolidate points which are equal, or very nearly so.  A value of
%    -1.0, on the other hand, will force the code to use every point, regardless
%    of duplication.
%
  addpath ( '../sandia_rules' );

  timestamp ( );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_WRITE_TESTS\n' );
  fprintf ( 1, '  Call SGMGA_WRITE_TEST with various arguments.\n' );
%
%  Set the point equality tolerance.
%
  tol = sqrt ( eps );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  All tests will use a point equality tolerance of %e\n', tol );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = 1.0;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 1, 1 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l2_ccxcc_iso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );

  dim_num = 2;
  for dim = 1 : dim_num
    importance(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 1, 1 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l2_ccxcc_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );

  dim_num = 3;
  for dim = 1 : dim_num
    importance(dim,1) = 1.0;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 1, 1, 1 ]';
  growth = [ 6, 6, 6 ]';
  np = [ 0, 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d3_l2_ccxccxcc_iso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );

  dim_num = 3;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 1, 1, 1 ]';
  growth = [ 6, 6, 6 ]';
  np = [ 0, 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d3_l2_ccxccxcc_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );

  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 3;
  rule = [ 1, 3 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l3_ccxgp_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );

  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 1, 4 ]';
  growth = [ 6, 3 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l2_ccxgl_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );

  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 1, 7 ]';
  growth = [ 6, 3 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l2_ccxlg_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );

  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 1, 8 ]';
  growth = [ 6, 3 ]';
  np = [ 0, 1 ]';
  np_sum = sum ( np(1:dim_num) );
  p(1) = 1.5;
  file_name = 'sgmga_d2_l2_ccxglg_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );

  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 2, 9 ]';
  growth = [ 6, 3 ]';
  np = [ 0, 2 ]';
  np_sum = sum ( np(1:dim_num) );
  p(1,1) = 0.5;
  p(2,1) = 1.5;
  file_name = 'sgmga_d2_l2_f2xgj_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );

  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 6, 10 ]';
  growth = [ 3, 4 ]';
  np = [ 1, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p(1) = 1.0;
  file_name = 'sgmga_d2_l2_gghxhgk_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );
%
%  LEVEL_MAX 1
%
  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 1;
  rule = [ 1, 1 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l1_ccxcc_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );
%
%  LEVEL_MAX 2 already done
%
%  LEVEL_MAX 3
%
  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 3;
  rule = [ 1, 1 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l3_ccxcc_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );
%
%  LEVEL_MAX 4
%
  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 4;
  rule = [ 1, 1 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l4_ccxcc_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );
%
%  LEVEL_MAX 5
%
  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 5;
  rule = [ 1, 1 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l5_ccxcc_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );
%
%  Dimension 3
%
  dim_num = 3;
  for dim = 1 : dim_num
    important(dim,1) = dim;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 2;
  rule = [ 1, 4, 5 ]';
  growth = [ 6, 3, 3 ]';
  np = [ 0, 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d3_l2_ccxglxgh_aniso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );
%
%  Rule 3, LEVEL_MAX 4
%
  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = 1.0;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 4;
  rule = [ 3, 3 ]';
  growth = [ 6, 6 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l4_gpxgp_iso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );
%
%  Rule 3, Growth 4, LEVEL_MAX 4
%
  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = 1.0;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 4;
  rule = [ 3, 3 ]';
  growth = [ 4, 4 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l4_gpsexgpse_iso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );
%
%  Rule 3, Growth 5, LEVEL_MAX 4
%
  dim_num = 2;
  for dim = 1 : dim_num
    important(dim,1) = 1.0;
  end
  level_weight = sgmga_importance_to_aniso ( dim_num, importance );
  level_max = 4;
  rule = [ 3, 3 ]';
  growth = [ 5, 5 ]';
  np = [ 0, 0 ]';
  np_sum = sum ( np(1:dim_num) );
  p = [];
  file_name = 'sgmga_d2_l4_gpmexgpme_iso';
  sgmga_write_test ( dim_num, level_weight, level_max, rule, growth, np, p, tol, ...
    file_name );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SGMGA_WRITE_TESTS:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  fprintf ( 1, '\n' );
  timestamp ( );

  rmpath ( '../sandia_rules' );

  return
end


text:
function candoit = canGetPrecon(problem)
% Checks whether a preconditioner was specified in the problem description.
%
% function candoit = canGetPrecon(problem)
%
% Returns true if a preconditioner was specified, false otherwise. Notice
% that even if this function returns false, it is still possible to call
% getPrecon, as the default preconditioner is simply the identity operator.
% This check function is mostly useful to tell whether that default
% preconditioner will be in use or not.
%
% See also: getPrecon getSqrtPrecon canGetSqrtPrecon getHessian

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, July 3, 2013.
% Contributors:
% Change log:

    candoit = isfield(problem, 'precon') || canGetSqrtPrecon(problem);

end

author: NicolasBoumal

repo: manopt


text:
function [ vS ] = RecSignalFreqSamples( vF, mF, vFIdx )
% ----------------------------------------------------------------------------------------------- %
%[ estFreq ] = EstimateSineFreqKay( vX, samplingFreq, estType )
% Estimates the frequency of a single Real Harmonic signal with arbitrary
% phase.
% Input:
%   - vX                -   Input Samples.
%                           The vector to be optimized. Initialization of
%                           the iterative process.
%                           Structure: Vector (numSamples X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - paramLambda       -   Parameter Lambda.
%                           The L1 Regularization parameter.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: (0, inf).
%   - numIterations     -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% Output:
%   - estFreq           -   Number of Iterations.
%                           Number of iterations of the algorithm.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range {1, 2, ...}.
% References
%   1.  Steven Kay - A Fast and Accurate Single Frequency Estimator.
% Remarks:
%   1.  It would work with complex numbers (Harmonic Signla) by changing:
%       vX(ii) * vX(ii + 1) -> vX(ii)' * vX(ii + 1).
%   2.  fds
% Known Issues:
%   1.  C
% TODO:
%   1.  D
% Release Notes:
%   -   1.0.000     09/08/2021  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

% arguments
%     vX (:, 1) {mustBeNumeric, mustBeReal}
%     samplingFreq (1, 1) {mustBeNumeric, mustBeReal, mustBePositive}
%     estType (1, 1) {mustBeNumeric, mustBeReal, mustBePositive, mustBeInteger, mustBeMember(estType, [1, 2])} = 2
% end

numSamples = size(mF, 2);

mFR = real(mF);
mFI = imag(mF);

vS = real(SolveBasisPursuitLp001([mFR(vFIdx, :); mFI(vFIdx, :)], [real(vF); imag(vF)]));
% vSRec = real(SolveBasisPursuitLp002(mF(vFIdx, :), vF));

% vS = complex(vSRec(1:numSamples), vSRec((numSamples + 1):end));

end

author: RoyiAvital

repo: StackExchangeCodes


text:
%
% LibQPEP: A Library for Globally Optimal Solving Quadratic Pose Estimation Problems (QPEPs),
%          It also gives highly accurate uncertainty description of the solutions.
%
%
% Article:
%      Wu,    J.,    Zheng,    Y.,    Gao,    Z.,    Jiang,    Y.,    Hu,    X.,    Zhu,    Y.,    Jiao,    J.,    Liu,    M. (2020)
%           Quadratic Pose Estimation Problems: Unified Solutions,
%           Solvability/Observability Analysis and Uncertainty Description
%           in A Globally Optimal Framework.
%
%
% Authors:      Jin Wu and Ming Liu
% Affiliation:  Hong Kong University of Science and Technology (HKUST)
% Emails:       jin_wu_uestc@hotmail.com; eelium@ust.hk
% Websites:     https://zarathustr.github.io
%               https://ram-lab.com

function sols = solver_WQ_approx(data)
[C0, C1_] = setup_elimination_template_solver_WQ_approx(data);
C1 = C0 \ C1_;
RR = [ - C1(end-12:end, :); eye(27)];
AM_ind = [38, 16, 1, 19, 2, 3, 21, 22, 4, 25, 5, 6, 28, 7, 8, 30, 31, 9, 34, 10, 11, 36, 37, 12, 39, 40, 13];
AM = RR(AM_ind, :);
[V, D] = eig(AM);
V = V ./ (ones(size(V, 1), 1) * V(1, :));
sols = complex(zeros(4, 27));
sols(2, :) = sqrt(V(10, :));
sols(1, :) = V(2, :) ./ (sols(2, :));
sols(3, :) = V(4, :) ./ (sols(1, :));
sols(4, :) = V(3, :) ./ (sols(1, :) .* V(16, :));
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
% Fit a piece-wise linear regression model.
% Here is the model
%
%  X \
%  | |
%  Q |
%  | /
%  Y
%
% where all arcs point down.
% We condition everything on X, so X is a root node. Q is a softmax, and Y is a linear Gaussian.
% Q is hidden, X and Y are observed.

X = 1;
Q = 2;
Y = 3;
dag = zeros(3,3);
dag(X,[Q Y]) = 1;
dag(Q,Y) = 1;
ns = [1 2 1]; % make X and Y scalars, and have 2 experts
dnodes = [2];
onodes = [1 3];
bnet = mk_bnet(dag, ns, 'discrete', dnodes, 'observed', onodes);

w = [-5 5];  % w(:,i) is the normal vector to the i'th decisions boundary
b = [0 0];  % b(i) is the offset (bias) to the i'th decisions boundary

mu = [0 0];
sigma = 1;
Sigma = repmat(sigma*eye(ns(Y)), [ns(Y) ns(Y) ns(Q)]);
W = [-1 1];
W2 = reshape(W, [ns(Y) ns(X) ns(Q)]);

bnet.CPD{1} = root_CPD(bnet, 1);
bnet.CPD{2} = softmax_CPD(bnet, 2, w, b);
bnet.CPD{3} = gaussian_CPD(bnet, 3, 'mean', mu, 'cov', Sigma, 'weights', W2);

% Check inference

x = 0.1;
ystar = 1;

engine = jtree_inf_engine(bnet);
[engine, loglik] = enter_evidence(engine, {x, [], ystar});
Qpost = marginal_nodes(engine, 2);

% eta(i,:) = softmax (gating) params for expert i
eta = [b' w'];

% theta(i,:) = regression vector for expert i
theta = [mu' W'];

% yhat(i) = E[y | Q=i, x] = prediction of i'th expert
x1 = [1 x]';
yhat = theta * x1;

% gate_prior(i,:) = Pr(Q=i | x)
gate_prior = normalise(exp(eta * x1));

% cond_lik(i) = Pr(y | Q=i, x)
cond_lik = (1/(sqrt(2*pi)*sigma)) * exp(-(0.5/sigma^2) * ((ystar - yhat) .* (ystar - yhat)));

% gate_posterior(i,:) = Pr(Q=i | x, y)
[gate_posterior, lik] = normalise(gate_prior .* cond_lik);

assert(approxeq(gate_posterior(:), Qpost.T(:)));
assert(approxeq(log(lik), loglik));

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [X,labels,class_names] = ml_load_data(filename,format,label_j)
%ML_LOAD_DATA Loads a data set in which the rows are samples and the
%columns are features/dimensions. The last column should contain the class
%label. If not specified
%
%
%
%   input -----------------------------------------------------------------
%
%       o filename : string,   path to file
%
%       o formate  : string,  [.txt,.csv,...]
%
%       o label_j  : (1 x 1), column index which containts the class labels
%
%   output ----------------------------------------------------------------
%
%       o X         : (N x D), dataset of N samples of dimension D
%
%       o labels    : (N x 1), N labels
%

data = [];
if ~exist('label_j','var'), label_j = []; end

if strcmp(format,'csv') == true

   data = readtable(filename);

elseif strcmp(format,'txt') == true

    data = load(filename);
else
   warning(['No such format: ' froamt ' implemented!']);
end

if isempty(label_j)
   X = data;
else

   X            = table2array(data(:,1:end-1));
   labels       = table2array(data(:,end));
   class_names  = unique(labels);

   % if data is string convert to numbers
   if iscellstr(X) == true
      x_unique   = unique(X);
      cellfind   = @(string)(@(cell_contents)(strcmp(string,cell_contents)));
      X_tmp      = zeros(size(X));
   %   bUnknown   = false;

   %   if sum(cellfun(cellfind('?'),x_unique)) == true, bUnknown=true;end

      v = 0;
      for i=1:length(x_unique)

           idx = cellfun(cellfind(x_unique{i}),X);

           if strcmp(x_unique{i},'?') == true
               X_tmp(idx) = -1;
           else
               X_tmp(idx) = v;
               v          = v + 1;
           end
      end

        X = X_tmp;
   end

   % if labels are strings convert them to numbers
   if iscellstr(class_names) == true
       labels_tmp = zeros(size(labels,1),1);
       cellfind   = @(string)(@(cell_contents)(strcmp(string,cell_contents)));

       for i=1:size(class_names,1)
           idx              = cellfun(cellfind(class_names{i}),labels);
           labels_tmp(idx)  = i;
       end

       labels = labels_tmp;

   end

end

end

author: epfl-lasa


text:
function W = randInitializeWeights(L_in, L_out)
%RANDINITIALIZEWEIGHTS Randomly initialize the weights of a layer with L_in
%incoming connections and L_out outgoing connections
%   W = RANDINITIALIZEWEIGHTS(L_in, L_out) randomly initializes the weights
%   of a layer with L_in incoming connections and L_out outgoing
%   connections.
%
%   Note that W should be set to a matrix of size(L_out, 1 + L_in) as
%   the column row of W handles the "bias" terms
%

% You need to return the following variables correctly
W = zeros(L_out, 1 + L_in);
% ====================== YOUR CODE HERE ======================
% Instructions: Initialize W randomly so that we break the symmetry while
%               training the neural network.
%
% Note: The first row of W corresponds to the parameters for the bias units
%

epsilon_init = sqrt(6) / sqrt(L_in + L_out);
W = rand(L_out, 1 + L_in) * 2 * epsilon_init - epsilon_init;

% =========================================================================

end

author: scruel

repo: Notes-ML-AndrewNg


text:
  function img = feldkamp(cg, ig, proj, varargin)
%|function img = feldkamp(cg, ig, proj, varargin)
%|
%| version modified by james balter and others to try accelereys "jacket"
%| UNDER DEVELOPMENT!
%|
%| FBP reconstruction of cone-beam tomography data collected with
%| a circular source trajectory.
%| See feldkamp_example.m for example.
%|
%| in:
%|	cg			ct_geom()
%|	ig			image_geom()
%|	proj	[ns,nt,na]	cone-beam projection views (line integrals)
%|
%| options:
%|	'window' [npad]		'ramp' (default), or 'hann', or array.
%|				if array, then use samples [-K/2, K/2).
%|	'offset_source'		distance from isocenter to perpendicular ray
%|				[the same units (e.g., mm) as pixel_size etc.]
%|				caution: probably should not be used
%|	'ia_skip' [int]		downsample in angle to save time for tests
%|	'use_mex' 0|1		backprojector: 0 for matlab, 1 for mex (default)
%|	'nthread'		default: jf('ncore')
%|
%| out:
%|	img	[nx,ny,nz]	reconstructed image
%|
%| References: Feldkamp, Davis, Kress, JOSA-A, 1(6):612-9, June 1984.
%| Notation here follows Fessler tomography book chapter (ask if interested).
%|
%| Copyright 2004-8-28 Nicole Caparanis, Patty Laskowsky, Taka Masuda,
%| and Jeff Fessler, University of Michigan
%| arc detector case contributed by Yingying Zhang 2005-6-13

if nargin == 1 && streq(cg, 'test')
	run_mfile_local('feldkamp_example')
return
end
if nargin < 3, ir_usage, end

if isnumeric(cg)
	img = feldkamp_old(cg, ig, proj, varargin{:});
return
end

% defaults
arg.use_mex = 1;
arg.window = 'ramp';
arg.offset_source = 0; % r_off, distance between rotation iso-center
			% and ray from source that is orthogonal to detector.
arg.ia_skip = 1;
arg.nthread = jf('ncore');
arg = vararg_pair(arg, varargin);

if cg.pitch ~= 0 || any(cg.zshifts ~= 0)
	fail('sorry, helical CT unsupported')
end

img = feldkamp_do(proj, ...
	cg.ds, cg.dt, cg.offset_s, cg.offset_t, arg.offset_source, ...
	cg.dsd, cg.dso, cg.dfs, cg.orbit, cg.orbit_start, ...
	ig.mask_or, ig.nz, ig.dx, ig.dy, ig.dz, ...
	[ig.offset_x ig.offset_y ig.offset_z], ...
	arg.window, arg.ia_skip, arg.use_mex, arg.nthread);
end % feldkamp()

 function img = feldkamp_old(proj, mask, varargin)
%function img = feldkamp_old(proj, mask, [options])
%
% FBP reconstruction of cone-beam tomography data collected with
% a circular source trajectory.
% See feldkamp_example.m for example.
%
% in:
%	proj	[ns,nt,na]	cone-beam projection views (line integrals)
%	mask	[nx,ny]		logical array of desired (2d) fov support
%
% options:
%	window	[npad]		'ramp' (default), or 'hann', or array.
%				if array, then use samples [-K/2, K/2).
%	'dx' 'dy' 'dz'		voxel dimensions
%	'nz'			# of slices
%	'center_xyz'		image center [pixels] (default: 0)
%	'ds' 'dt'		detector sample spacing (horiz, vert)
%	'offset_st'		detector offset_s,t [integer number of samples
%				or fraction thereof, for 1/4-detector offset.
%	'orbit'			default 360 degrees
%	'orbit_start'		default 0 degrees
%	'dis_src_det'
%	'dis_iso_det'
%	'dis_foc_src'		default inf (flat), or 0 for arc
%	'offset_source'		distance from isocenter to perpendicular ray
%				[the same units (e.g., mm) as pixel_size etc.]
%	'ia_skip' [int]		downsample in angle to save time for tests
%	'use_mex' 0|1		backprojector: 0 for matlab, 1 for mex (default)
%
% out:
%	img	[nx,ny,nz]		reconstructed image
%
% offset_s is an integer or fraction thereof, e.g., 0.25,
% relative to centerline between two central channels.
%
% References: Feldkamp, Davis, Kress, JOSA-A, 1(6):612-9, June 1984.
% Fessler tomography chapter.
%
% Copyright 2004-8-28 Nicole Caparanis, Patty Laskowsky, Taka Masuda,
% and Jeff Fessler, University of Michigan
% arc detector case contributed by Yingying Zhang 2005-6-13

warn 'this is obsolete use of feldkamp.m - use new style!'

% defaults
arg.use_mex = 1; % force MATLAB
arg.window = 'ramp';
arg.dx = nan;
arg.dy = [];
arg.dz = [];
arg.nz = size(proj,2);
arg.center_xyz = [0 0 0]; % center_x,y,z (pixels)
arg.flip_y = 1;

arg.dsd = []; % 'dis_src_det'
arg.dod	= []; %	'dis_iso_det'
arg.dfs = inf; % 'dis_foc_src' default flat
arg.offset_source = 0; % r_off, distance between rotation iso-center
			% and ray from source that is orthogonal to detector.

arg.ds = nan;
arg.dt = [];
arg.offset_st = [0 0]; % offset_s,t
arg.ia_skip = 1;
arg.orbit = 360;
arg.orbit_start = 0;
arg = vararg_pair(arg, varargin, 'subs', ...
	{'dis_src_det', 'dsd';
	'dis_iso_det', 'dod';
	'dis_foc_src', 'dfs'});
if isnan(arg.ds) || isnan(arg.dx), fail('ds and dx required'), end

if isempty(arg.dy), arg.dy = -arg.dx; end
if isempty(arg.dz), arg.dz = arg.dx; end
if isempty(arg.dt), arg.dt = arg.ds; end

arg.dso = arg.dsd - arg.dod; % src to "isocenter" distance

if ~isempty(arg.zshifts), error 'todo: helical feldkamp not done', end

img = feldkamp_do(proj, ...
	arg.ds, arg.dt, arg.offset_st(1), arg.offset_st(2), ...
	arg.offset_source, ...
	arg.dsd, arg.dso, arg.dfs, ...
	arg.orbit, arg.orbit_start, ...
	mask, arg.nz, arg.dx, arg.dy, arg.dz, arg.center_xyz, ...
	arg.window, arg.ia_skip, arg.use_mex, 1);

end % feldkamp_old()

%
% feldkamp_do()
%
function img = feldkamp_do(proj, ...
	ds, dt, offset_s, offset_t, offset_source, ...
	dsd, dso, dfs, orbit, orbit_start, ...
	mask, nz, dx, dy, dz, offset_xyz, ...
	window, ia_skip, use_mex, nthread)

% step 1: weight the projections as in fan-beam case
proj = feldkamp_weight1(proj, ds, dt, offset_s, offset_t, dsd, dso, dfs);

% step 2: filter the (zero padded) projections
[ns nt na] = size(proj);
if use_mex
	pad1 = 0;
else
	pad1 = 2;
end
proj = feldkamp_filter(proj, window, dsd, dfs, ds, pad1);

if pad1 % trick: zero at end saves indexing in loop
%	proj = [proj; zeros(2, nt, na)];
	proj(:,end+1,:) = 0;
end

% step 3: cone-beam backprojection of the filtered views
cpu etic
img = feldkamp_back(proj, ns, nt, na, ...
	ds, dt, offset_s, offset_t, offset_source, ...
	dsd, dso, dfs, orbit, orbit_start, ...
	mask, nz, dx, dy, dz, offset_xyz, ...
	ia_skip, use_mex, nthread);
cpu etoc 'fdk backprojection cpu time:'

end % feldkamp_do()

%
% feldkamp_weight1()
% step 1: weight the projections as in fan-beam case
%
function proj = feldkamp_weight1(proj, ds, dt, offset_s, offset_t, ...
	dsd, dso, dfs);
[ns nt na] = size(proj);
ss = ([-(ns-1)/2:(ns-1)/2]' - offset_s) * ds;
tt = ([-(nt-1)/2:(nt-1)/2]' - offset_t) * dt;

[ss tt] = ndgrid(ss, tt);
if isinf(dfs) % flat
	ww1 = dso * sqrt(1 + (tt/dsd).^2) ./ sqrt(dsd^2 + ss.^2 + tt.^2);
elseif dfs == 0 % arc
	ww1 = (dso/dsd) * cos(ss ./ (dsd * sqrt(1 + (tt/dsd).^2)));
else
	error 'other configurations not implemented'
end

for ia=1:na % same weighting for each view angle
	proj(:,:,ia) = proj(:,:,ia) .* ww1;
end
end % feldkamp_weight1()

%
% feldkamp_filter()
% step 2: filter the (zero padded) projections
%
function proj = feldkamp_filter(proj, window, dsd, dfs, ds, pad1)
[ns nt na] = size(proj);
npadh = 2^ceil(log2(2*ns-1));
printf('ns=%d npadh=%d', ns, npadh)

if isinf(dfs)
	H = fan_filter('flat', npadh, ds, [], window);	% [nb,1]
elseif dfs == 0
	H = fan_filter('arc', npadh, ds, dsd, window);
end
H = ds * H; % differential for discrete-space convolution vs integral

%% START JACKET
GPU = 1;
if GPU,
    proj = gsingle( proj );
    H = gsingle( H );
    proj_large = gzeros( npadh, size(proj,2), size(proj,3) );
else,
    proj_large = zeros( npadh, size(proj,2), size(proj,3) );
end
%% END JACKET

tic;
proj_large(1:size(proj,1),:,:) = proj;
proj = ifft( fft(proj_large) .* repmat(H, [1 nt na]) );
proj = proj(1:(ns+pad1),:,:); % trick: extra zero at end saves indexing in loop
toc

%% START JACKET
if GPU,
    proj = double( proj );
end
%% END JACKET

end % feldkamp_filter()

%
% fan_filter()
% apodized filter frequency response
%
function H = fan_filter(type, n, ds, dsd, window)

if streq(type, 'flat')
	h = fbp_ramp('flat', n, ds);
else
	h = fbp_ramp('arc', n, ds, dsd);
end
H = reale(fft(fftshift(h)));

if ischar(window)
	if streq(window, 'ramp')
		window = ones(n,1);
	elseif streq(window, 'hann')
		window = ir_hann_periodic(n);
	else
		error 'unknown window'
	end
elseif length(window) ~= n
	error 'bad window length'
end

H = H .* fftshift(window);
end % fan_filter()

%
% feldkamp_back()
% step 3: cone-beam backprojection of the filtered views
%
function img = feldkamp_back(proj, ns, nt, na, ...
	ds, dt, offset_s, offset_t, offset_source, ...
	dsd, dso, dfs, orbit, orbit_start, ...
	mask, nz, dx, dy, dz, offset_xyz, ...
	ia_skip, use_mex, nthread)
[nx ny] = size(mask);

betas = deg2rad(orbit_start + orbit * [0:na-1] / na); % [na] source angles

% mex backprojector
if use_mex
	proj = single(proj);
	nthread = int32(nthread);

	if use_mex == 1 % loop in mex
		proj = permute(proj, [2 1 3]); % ts
		tmp = jf_mex('fdk,ts,back', int32([nx ny nz]), [dx dy dz], ...
				offset_xyz, uint8(mask), ...
				dso, dsd, dfs, [ds dt], [offset_s offset_t], ...
				proj, betas, nthread);
		img = double6(tmp);
		img = permute(img, [2 3 1]); % zxy -> xyz

	elseif use_mex == 2 % loop in matlab (for testing)
		proj = permute(proj, [2 1 3]); % ts
		img = 0;
		for ia=1:na
			ticker(mfilename, ia, na)
			% note: 2006-5-30: replaced -dy with dy
			tmp = jf_mex('fdk,ts,back', ...
				int32([nx ny nz]), [dx dy dz], ...
				offset_xyz, uint8(mask), ...
				dso, dsd, dfs, [ds dt], [offset_s offset_t], ...
				proj(:,:,ia), betas(ia), nthread);
			tmp = double6(tmp);
			img = img + tmp;
		end
		img = permute(img, [2 3 1]); % zxy -> xyz

	elseif use_mex == 3 % fdk,st (for testing only - slower!)
		tmp = jf_mex('fdk,st,back', int32([nx ny nz]), [dx dy dz], ...
				offset_xyz, uint8(mask), ...
				dso, dsd, dfs, [ds dt], [offset_s offset_t], ...
				proj, betas, nthread);
		img = double6(tmp);

	else
		fail 'bug'
	end

	% final "\der angle" scale:
	img = (0.5 * deg2rad(abs(orbit)) / (na/ia_skip)) * img;
return
end

% matlab backprojector (slower)

fprintf('We are using slower version!\n');

% precompute as much as possible
wx = (nx-1)/2 + offset_xyz(1);
wy = (ny-1)/2 + offset_xyz(2);
wz = (nz-1)/2 + offset_xyz(3);
[xc yc] = ndgrid(([0:nx-1] - wx) * dx, ([0:ny-1] - wy) * dy);
zc = ([0:nz-1] - wz) * dz;

if 0 % limit back-projection to FOV?  removed 2008-10-9
	rr = sqrt(xc.^2 + yc.^2); % [nx,ny]
	smax = ((ns-1)/2-abs(offset_s)) * ds; % maximum detector s coordinate

	if isinf(dfs)
		gamma_max = atan(smax/dsd);
	elseif dfs == 0
		gamma_max = smax / dsd;
	end

	rmax = dso * sin(gamma_max);
	mask = mask & (rr < rmax);
end
clear wx wy wz rr smax rmax

xc = xc(mask); % [np] pixels within mask
yc = yc(mask);

% loop over slices
GPU = 0;
if GPU,
    img = gzeros([size(mask) nz]) + i;
    img = img - i;
else,
    img = zeros([size(mask) nz]);
end;
ticker reset
tic;
for iz=1:nz

	% loop over each projection angle
	img2 = 0;
	for ia=1:na
		ticker(mfilename, [iz ia], [nz na])
		beta = betas(ia);

		x_beta = +xc * cos(beta) + yc * sin(beta);
		y_beta = -xc * sin(beta) + yc * cos(beta);

		% detector indices
		mag = dsd ./ (dso - y_beta);

		if isinf(dfs)
			sprime = mag .* x_beta;
		elseif (dfs == 0)
			d_loop = dso - y_beta;
			r_loop = x_beta - offset_source;
			sprime = dsd * atan2(r_loop, d_loop);
		end

		tprime = mag * zc(iz);
		bh = sprime / ds + (ns+1)/2 + offset_s;
		bv = tprime / dt + (nt+1)/2 + offset_t;

		% bi-linear interpolation:
		ih = floor(bh); % left bin
		iv = floor(bv);
		igood = 1<=ih & ih<ns & 1<= iv & iv<nt;

		ih(~igood) = ns+1; % trick! point at harmless zeros
		iv(~igood) = 1;

		wr = bh - ih;	% left weight
		wl = 1 - wr;	% right weight
		wu = bv - iv;	% upper weight
		wd = 1 - wu;	% lower weight

		sdim = size(proj);
		ia1 = ia * ones(size(ih));

		p1 =	wl .* proj(sub2ind(sdim, ih,iv,ia1)) + ...
			wr .* proj(sub2ind(sdim, ih+1,iv,ia1));
		p2 =	wl .* proj(sub2ind(sdim, ih,iv+1,ia1)) + ...
			wr .* proj(sub2ind(sdim, ih+1,iv+1,ia1));

		p0 = wu .* p1 + wd .* p2; % vertical interpolation

		if isinf(dfs)
			p0 = p0 .* mag.^2; % back-projection weighting for flat
		elseif dfs == 0
			p0 = p0 .* (dsd.^2) ./ (r_loop.^2 + d_loop.^2);
		end

		img2 = img2 + p0;
	end % ia

	img(:,:,iz) = embed(img2, mask);
end % iz
toc

img = (0.5 * deg2rad(orbit) / (na/ia_skip)) * img; % final "\der angle" scale
img = double( img );

end % feldkamp_back()

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function [knn_dist, nn_dist] = smooth_knn_dist(distances, k, varargin)
%SMOOTH_KNN_DIST Compute a continuous version of the distance to the kth
% nearest neighbor. That is, this is similar to knn-distance but allows
% continuous k values rather than requiring an integral k. In esscence we
% are simply computing the distance such that the cardinality of fuzzy set
% we generate is k.
%
% [knn_dist, nn_dist] = SMOOTH_KNN_DIST(distances, k, local_connectivity,
% n_iter, bandwidth)
%
% Parameters
% ----------
% distances: array of size (n_samples, n_neighbors)
%     Distances to nearest neighbors for each samples. Each row should be a
%     sorted list of distances to a given sample's nearest neighbors.
%
% k: double
%     The number of nearest neighbors to approximate for.
%
% n_iter: double (optional, default 64)
%     We need to binary search for the correct distance value. This is the
%     max number of iterations to use in such a search.
%
% local_connectivity: double (optional, default 1)
%     The local connectivity required -- i.e. the number of nearest
%     neighbors that should be assumed to be connected at a local level.
%     The higher this value the more connected the manifold becomes
%     locally. In practice this should be not more than the local intrinsic
%     dimension of the manifold.
%
% bandwidth: double (optional, default 1)
%     The target bandwidth of the kernel, larger values will produce
%     larger return values.
%
% Returns
% -------
% knn_dist: array of size (n_samples, 1)
%     The distance to kth nearest neighbor, as suitably approximated.
%
% nn_dist: array of size (n_samples, 1)
%     The distance to the 1st nearest neighbor for each point.
%
%   AUTHORSHIP
%   Math Lead & Primary Developer:  Connor Meehan <connor.gw.meehan@gmail.com>
%   Secondary Developer: Stephen Meehan <swmeehan@stanford.edu>
%   Bioinformatics Lead:  Wayne Moore <wmoore@stanford.edu>
%   Provided by the Herzenberg Lab at Stanford University
%   License: BSD 3 clause
%
    p=parseArguments();
    parse(p,varargin{:});
    args=p.Results;
    local_connectivity = args.local_connectivity;
    n_iter = args.n_iter;
    bandwidth = args.bandwidth;
    same_set = args.same_set;

    SMOOTH_K_TOLERANCE = 1e-5;
    MIN_K_DIST_SCALE = 1e-3;

    height = size(distances,1);

    target = log2(k)*bandwidth;

    lo = zeros(height,1);
    hi = Inf(height,1);
    mid = ones(height,1);

    zero_dists = sum(distances == 0, 2);

    if any(zero_dists==size(distances,2))
        warning('There are at least n_neighbors identical data points in the raw data. Results may be inaccurate.');
    end

    index = floor(local_connectivity);
    interpolation = local_connectivity - index;
    aug_dists = [lo distances repmat(max(distances, [], 2), [1 index+1])];

    idx = sub2ind(size(aug_dists), (1:height)', zero_dists + index + 1);
    rho = aug_dists(idx) + interpolation*(aug_dists(idx) - aug_dists(idx+height));

    if same_set
        d = distances(:,2:end) - rho;
    else
        d = distances(:,1:end) - rho;
    end

    for n = 1:n_iter

        summands = exp(-max(0, d./repmat(mid, [1 size(d, 2)])));

        psum = sum(summands, 2);

        if all(abs(psum - target) < SMOOTH_K_TOLERANCE)
            break
        end

        b = ~(psum > target).*hi;
        b(isnan(b)) = 0;

        hi = (psum > target).*mid + b;
        lo = (psum > target).*lo + ~(psum > target).*mid;

        c = (psum > target).*((lo + hi)/2);
        c(isnan(c)) = 0;

        mid = c + ~(psum > target).*min(2*lo, (lo + hi)/2);

    end

    result = mid;

    result = max(result, MIN_K_DIST_SCALE * ((rho > 0).*mean(distances,2) + (rho == 0).*mean(mean(distances))));

    knn_dist = result;
    nn_dist = rho;
end

    function p=parseArguments(varargin)
        p = inputParser;
        addParameter(p,'local_connectivity', 1);
        addParameter(p,'n_iter', 64);
        addParameter(p,'bandwidth', 1);
        addParameter(p,'same_set','true');
    end


text:
function str = convertCell2Str(cellStrArray)
%% Convert a cell array of strings to a single char array
%
%  Parameters:
%      cellStrArray   A char array or cell array containing char
%      str            (output) A single string representation of input

%% Test cellStrArray and convert
if isempty(cellStrArray)
    str = '';
elseif ischar(cellStrArray)
    str = ['[' cellStrArray ']'];
else
    str = ['[' cellStrArray{1} ']'];
    for k = 2:length(cellStrArray)
        str = [str ' [' cellStrArray{k} ']']; %#ok<*AGROW>
    end
end

author: VisLab

repo: EEG-Clean-Tools