text:
% Matlab Wrapper for SVMlight
% Version 0.92, August 2002
%
% Wrapper functions:
%   svml - Wrapper for SVMlight
%   svmltrain - Wrapper for SVMlight: Training
%   svmlfwd - Wrapper for SVMlight: Prediction
%   demsvml1 - Demo program for Matlab SVMlight wrapper
%
% Auxiliary functions to generate and read data files
%   svmlread - Read a data file generated by SVM light
%   svmlwrite - Write matrix into data file for SVM light
%   svmlopt - Generate/alter options structure for SVM light
%   consist - Check that arguments are consistent.
%
% Calling routines
%   svm_classify - Interface to SVM light, classification module
%   svm_learn - Interface to SVM light, learning module
%
%
% Copyright (c) by Anton Schwaighofer (2002)
% mailto:anton.schwaighofer@gmx.net
%
% This program is released unter the GNU General Public License.
%

% $Revision: 1.3 $ $Date: 2002/08/09 20:25:26 $
%
% 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 2
% 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, write to the
% Free Software Foundation, Inc.,
% 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
%

author: viggin

repo: domain-adaptation-toolbox


text:
function [membership,means,rms] = kmeansML(k,data,varargin)
% [membership,means,rms] = kmeansML(k,data,...)
%
% Multi-level kmeans.
% Tries very hard to always return k clusters.
%
% INPUT
%	k		Number of clusters
% 	data		dxn matrix of data points
%	'maxiter'	Max number of iterations. [30]
%	'dtol'		Min change in center locations. [0]
%	'etol'		Min percent change in RMS error. [0]
%	'ml'		Multi-level? [true]
%	'verbose'	Verbose level. [0]
%			    0 = none
%			    1 = textual
%			    2 = visual
%
% OUTPUT
% 	membership	1xn cluster membership vector
% 	means		dxk matrix of cluster centroids
%	rms		RMS error of model
%
% October 2002
% David R. Martin <dmartin@eecs.berkeley.edu>

% process options
maxIter = 30;
dtol = 0;
etol = 0;
ml = true;
verbose = 0;
for i = 1:2:numel(varargin),
  opt = varargin{i};
  if ~ischar(opt), error('option names not a string'); end
  if i==numel(varargin), error(sprintf('option ''%s'' has no value',opt)); end
  val = varargin{i+1};
  switch opt,
   case 'maxiter', maxIter = max(1,val);
   case 'dtol', dtol = max(0,val);
   case 'etol', etol = max(0,val);
   case 'ml', ml = val;
   case 'verbose', verbose = val;
   otherwise, error(sprintf('invalid option ''%s''',opt));
  end
end

[membership,means,rms] = ...
    kmeansInternal(k,data,maxIter,dtol,etol,ml,verbose,1);

function [membership,means,rms] = kmeansInternal(...
    k,data,maxIter,dtol,etol,ml,verbose,retry)

[d,n] = size(data);
perm = randperm(n);

% compute initial means
rate = 3;
minN = 50;
coarseN = round(n/rate);
if ~ml | coarseN < k | coarseN < minN,
  % pick random points as means
  means = data(:,perm(1:k));
else
  % recurse on random subsample to get means
  coarseData = data(:,perm(1:coarseN));
  [coarseMem,means] = ...
      kmeansInternal(k,coarseData,maxIter,dtol,etol,ml,verbose,0);
end

% Iterate.
iter = 0;
rms = inf;
if verbose>0, fwrite(2,sprintf('kmeansML: n=%d d=%d k=%d [',n,d,k)); end
while iter < maxIter,
  if verbose>0, fwrite(2,'.'); end
  iter = iter + 1;
  % Compute cluster membership and RMS error.
  rmsPrev = rms;
  [membership,rms] = computeMembership(data,means);
  % Compute new means and cluster counts.
  prevMeans = means;
  [means,counts] = computeMeans(k,data,membership);
  % The error should always decrease.
  if rms > rmsPrev, error('bug: rms > rmsPrev'); end
  % Check for convergence.
  rmsPctChange = 2 * (rmsPrev - rms) / (rmsPrev + rms + eps);
  maxMoved = sqrt(max(sum((prevMeans-means).^2)));
  if rmsPctChange <= etol & maxMoved <= dtol, break; end
  % Visualize.
  if verbose>1, kmeansVis(data,membership,means); end
end
[membership,rms] = computeMembership(data,means);
if verbose>0, fwrite(2,sprintf('] rms=%.3g\n',rms)); end

% If there's an empty cluster, then re-run kmeans.
% Retry a fixed number of times.
maxRetries = 3;
if find(counts==0),
  if retry < maxRetries,
    disp('Warning: Re-runing kmeans due to empty cluster.');
    [membership,means] = kmeansInternal( ...
        k,data,maxIter,dtol,etol,ml,verbose,retry+1);
  else
    disp('Warning: There is an empty cluster.');
  end
end

function [membership,rms] = computeMembership(data,means)
z = distSqr(data,means);
[d2,membership] = min(z,[],2);
rms = sqrt(mean(d2));

function [means,counts] = computeMeans(k,data,membership)
[d,n] = size(data);
means = zeros(d,k);
counts = zeros(1,k);
for i = 1:k,
  ind = find(membership==i);
  counts(i) = length(ind);
  means(:,i) = sum(data(:,ind),2) / max(1,counts(i));
end

%  for i = 1:n,
%    j = membership(i);
%    means(:,j) = means(:,j) + data(:,i);
%    counts(j) = counts(j) + 1;
%  end
%  for j = 1:k,
%    means(:,j) = means(:,j) / max(1,counts(j));
%  end

author: Cloud-CV

repo: object-proposals


text:
classdef PTKEmphysemaImage < PTKPlugin
    % PTKEmphysemaImage. Plugin for computing the percentage of emphysema
    %     voxels
    %
    %     This is a plugin for the Pulmonary Toolkit. Plugins can be run using
    %     the gui, or through the interfaces provided by the Pulmonary Toolkit.
    %     See PTKPlugin.m for more information on how to run plugins.
    %
    %     Plugins should not be run directly from your code.
    %
    %
    %
    %     Licence
    %     -------
    %     Part of the TD Pulmonary Toolkit. https://github.com/tomdoel/pulmonarytoolkit
    %     Author: Tom Doel, 2013.  www.tomdoel.com
    %     Distributed under the GNU GPL v3 licence. Please see website for details.
    %

    properties
        ButtonText = 'Emphysema'
        ToolTip = 'Shows emphysema-like voxels'
        Category = 'Analysis'

        AllowResultsToBeCached = true
        AlwaysRunPlugin = false
        Context = PTKContextSet.LungROI
        PluginType = 'ReplaceOverlay'
        HidePluginInDisplay = false
        FlattenPreviewImage = true
        PTKVersion = '2'
        ButtonWidth = 6
        ButtonHeight = 1
        GeneratePreview = false
    end

    methods (Static)
        function emphysema_results_image = RunPlugin(dataset, context, reporting)
            roi = dataset.GetResult('PTKLungROI');
            context_no_airways = dataset.GetResult('PTKGetMaskForContextExcludingAirways', context);
            roi.ResizeToMatch(context_no_airways);
            [~, emphysema_results_image] = PTKComputeEmphysemaFromMask(roi, context_no_airways);
            emphysema_results_image.ImageType = PTKImageType.Colormap;
        end
    end
end

author: tomdoel

repo: pulmonarytoolkit


text:
function varargout = maskToCERRStructure(maskM, isUniform, scanNum, strname, planC)
%"maskToCERRStructure"
%   Adds to planC a new structure whose contours are derived from maskM.
%   MaskM must be registered to either the uniform or non uniform CT scan.
%
%   Uses nearest neighbor interpolation for slices.
%
%   If registered to the uniform, isUniform should be 1, else 0.
%
%   The structure is given the name passed into strname, or 'Imported
%   Structure' if no name is passed.
%
%JRA 08/17/04
%
%Usage:
%GENERAL:
%   function planC = maskToCERRStructure(maskM, isUniform, scanNum, strname, planC)
%OR if CERR is open:
%   function planC = maskToCERRStructure(maskM, isUniform, scanNum, strname)
%
% 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/>.

global stateS

if ~exist('planC','var')
    global planC
end
indexS = planC{end};

if ~exist('strname','var')
    strname = 'Imported Structure';
end

%Check sizes, auto size detection.
siz = size(maskM);
unisiz = getUniformScanSize(planC{indexS.scan}(scanNum));
normsiz = size(getScanArray(planC{indexS.scan}(scanNum)));
if numel(siz) < 3
    siz(3) = 1;
end
if numel(normsiz) < 3
    normsiz(3) = 1;
end
if ~exist('isUniform','var')
    if isequal(siz, unisiz)
        isUniform = 1;
    elseif isequal(siz, normsiz)
        isUniform = 0;
    else
        error('maskM does not match dimension of uniform or nonuniform dataset.');
    end
else
    if (isUniform && ~isequal(siz, unisiz)) || (~isUniform && ~isequal(siz, normsiz))
        error('maskM does not match dimension of uniform or nonuniform dataset.');
    end
end

if ~isUniform
    %If registered to CT, just get contour info.
    [contourS, ~] = maskToPoly(maskM, 1:siz(3), scanNum, planC);
else
    %If registered to uniformized data, use nearest slice neighbor
    %interpolation.
    [~, ~, zUni] = getUniformScanXYZVals(planC{indexS.scan}(scanNum));

%     [xUni, yUni, zUni] = getUniformizedXYZVals(planC);
    [~, ~, zSca] = getScanXYZVals(planC{indexS.scan}(scanNum));

    tmpM = false(normsiz);

    for i=1:normsiz(3)
        zVal = zSca(i);
        uB = find(zUni > zVal, 1 );
        lB = find(zUni <= zVal, 1, 'last' );
        if normsiz(3) > 1 && (isempty(uB) || isempty(lB))
            continue
        end
        if abs(zUni(uB) - zVal) < abs(zUni(lB) - zVal)
            tmpM(:,:,i) = logical(maskM(:,:,uB));
        else
            tmpM(:,:,i) = logical(maskM(:,:,lB));
        end
    end
    %Get contour info.
    [contourS, ~] = maskToPoly(tmpM, 1:normsiz(3), scanNum, planC);
 end

%Make an empty structure, assign name/contour.
newstr = newCERRStructure(scanNum, planC);
newstr.contour = contourS;
newstr.structureName = strname;
newstr.associatedScan = scanNum;
newstr.assocScanUID = planC{indexS.scan}(scanNum).scanUID;
numStructs = length(planC{indexS.structures});

%Append new structure to planC.
if ~isempty(planC{indexS.structures})
    planC{indexS.structures} = dissimilarInsert(planC{indexS.structures}, newstr, numStructs+1, []);
else
    planC{indexS.structures} = newstr;
    planC{indexS.structureArrayMore}(scanNum).indicesArray = [];
    planC{indexS.structureArrayMore}(scanNum).bitsArray = [];
end

%Update uniformized data.
if strcmpi(stateS.optS.createUniformizedDataset,'yes')
    planC = updateStructureMatrices(planC, numStructs+1);
end

%Set varargout if requested.
if nargout > 0
    varargout{1} = planC;
end

if isfield(stateS,'handle') && ishandle(stateS.handle.CERRSliceViewer)
    stateS.structsChanged = 1;
    sliceCallBack('refresh');
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [res] = readResults(separ, offs, iter, avg, peval)
% [dpixc, dpix_ind, blinkmat, p] = readSimul(separ, offs, iter, avg, p)
% Reads simulations with required 'sep', 'offs', 'iter' and makes
% avgeraging of poisson noise 'avg' times

namedir_data = [peval.prename num2str(100*separ) 'offset_' num2str(offs)];
namefile_data = [namedir_data '-iter_' num2str(iter)];
namefile_res = [namefile_data '_avg' num2str(avg)];
res = load ([peval.path_res namedir_data '/' namefile_res '.mat']);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
function varargout = bbApply( action, varargin )
% Functions for manipulating bounding boxes (bb).
%
% A bounding box (bb) is also known as a position vector or a rectangle
% object. It is a four element vector with the fields: [x y w h]. A set of
% n bbs can be stores as an [nx4] array, most funcitons below can handle
% either a single or multiple bbs. In addtion, typically [nxm] inputs with
% m>4 are ok (with the additional columns ignored/copied to the output).
%
% bbApply contains a number of utility functions for working with bbs. The
% format for accessing the various utility functions is:
%  outputs = bbApply( 'action', inputs );
% The list of functions and help for each is given below. Also, help on
% individual subfunctions can be accessed by: "help bbApply>action".
%
% Compute area of bbs.
%   bb = bbApply( 'area', bb )
% Shift center of bbs.
%   bb = bbApply( 'shift', bb, xdel, ydel )
% Get center of bbs.
%   cen = bbApply( 'getCenter', bb )
% Get bb at intersection of bb1 and bb2 (may be empty).
%   bb = bbApply( 'intersect', bb1, bb2 )
% Get bb that is union of bb1 and bb2 (smallest bb containing both).
%   bb = bbApply( 'union', bb1, bb2 )
% Resize the bbs (without moving their centers).
%   bb = bbApply( 'resize', bb, hr, wr, [ar] )
% Fix bb aspect ratios (without moving the bb centers).
%   bbr = bbApply( 'squarify', bb, flag, [ar] )
% Draw single or multiple bbs to image (calls rectangle()).
%   hs = bbApply( 'draw', bb, [col], [lw], [ls], [prop], [ids] )
% Embed single or multiple bbs directly into image.
%  I = bbApply( 'embed', I, bb, [varargin] )
% Crop image regions from I encompassed by bbs.
%   [patches, bbs] = bbApply('crop',I,bb,[padEl],[dims])
% Convert bb relative to absolute coordinates and vice-versa.
%   bb = bbApply( 'convert', bb, bbRef, isAbs )
% Randomly generate bbs that fall in a specified region.
%   bbs =  bbApply( 'random', pRandom )
% Convert weighted mask to bbs.
%   bbs = bbApply('frMask',M,bbw,bbh,[thr])
% Create weighted mask encoding bb centers (or extent).
%   M = bbApply('toMask',bbs,w,h,[fill],[bgrd])
%
% USAGE
%  varargout = bbApply( action, varargin );
%
% INPUTS
%  action     - string specifying action
%  varargin   - depends on action, see above
%
% OUTPUTS
%  varargout  - depends on action, see above
%
% EXAMPLE
%
% See also bbApply>area bbApply>shift bbApply>getCenter bbApply>intersect
% bbApply>union bbApply>resize bbApply>squarify bbApply>draw bbApply>crop
% bbApply>convert bbApply>random bbApply>frMask bbApply>toMask
%
% Piotr's Computer Vision Matlab Toolbox      Version 3.30
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

%#ok<*DEFNU>
varargout = cell(1,max(1,nargout));
[varargout{:}] = feval(action,varargin{:});
end

function a = area( bb )
% Compute area of bbs.
%
% USAGE
%  bb = bbApply( 'area', bb )
%
% INPUTS
%  bb     - [nx4] original bbs
%
% OUTPUTS
%  a      - [nx1] area of each bb
%
% EXAMPLE
%  a = bbApply('area', [0 0 10 10])
%
% See also bbApply
a=prod(bb(:,3:4),2);
end

function bb = shift( bb, xdel, ydel )
% Shift center of bbs.
%
% USAGE
%  bb = bbApply( 'shift', bb, xdel, ydel )
%
% INPUTS
%  bb     - [nx4] original bbs
%  xdel   - amount to shift x coord of each bb left
%  ydel   - amount to shift y coord of each bb up
%
% OUTPUTS
%  bb     - [nx4] shifted bbs
%
% EXAMPLE
%  bb = bbApply('shift', [0 0 10 10], 1, 2)
%
% See also bbApply
bb(:,1)=bb(:,1)-xdel; bb(:,2)=bb(:,2)-ydel;
end

function cen = getCenter( bb )
% Get center of bbs.
%
% USAGE
%  cen = bbApply( 'getCenter', bb )
%
% INPUTS
%  bb     - [nx4] original bbs
%
% OUTPUTS
%  cen    - [nx1] centers of bbs
%
% EXAMPLE
%  cen = bbApply('getCenter', [0 0 10 10])
%
% See also bbApply
cen=bb(:,1:2)+bb(:,3:4)/2;
end

function bb = intersect( bb1, bb2 )
% Get bb at intersection of bb1 and bb2 (may be empty).
%
% USAGE
%  bb = bbApply( 'intersect', bb1, bb2 )
%
% INPUTS
%  bb1    - [nx4] first set of bbs
%  bb2    - [nx4] second set of bbs
%
% OUTPUTS
%  bb     - [nx4] intersection of bbs
%
% EXAMPLE
%  bb = bbApply('intersect', [0 0 10 10], [5 5 10 10])
%
% See also bbApply bbApply>union
n1=size(bb1,1); n2=size(bb2,1);
if(n1==0 || n2==0), bb=zeros(0,4); return, end
if(n1==1 && n2>1), bb1=repmat(bb1,n2,1); n1=n2; end
if(n2==1 && n1>1), bb2=repmat(bb2,n1,1); n2=n1; end
assert(n1==n2);
lcsE=min(bb1(:,1:2)+bb1(:,3:4),bb2(:,1:2)+bb2(:,3:4));
lcsS=max(bb1(:,1:2),bb2(:,1:2)); empty=any(lcsE<lcsS,2);
bb=[lcsS lcsE-lcsS]; bb(empty,:)=0;
end

function bb = union( bb1, bb2 )
% Get bb that is union of bb1 and bb2 (smallest bb containing both).
%
% USAGE
%  bb = bbApply( 'union', bb1, bb2 )
%
% INPUTS
%  bb1    - [nx4] first set of bbs
%  bb2    - [nx4] second set of bbs
%
% OUTPUTS
%  bb     - [nx4] intersection of bbs
%
% EXAMPLE
%  bb = bbApply('union', [0 0 10 10], [5 5 10 10])
%
% See also bbApply bbApply>intersect
n1=size(bb1,1); n2=size(bb2,1);
if(n1==0 || n2==0), bb=zeros(0,4); return, end
if(n1==1 && n2>1), bb1=repmat(bb1,n2,1); n1=n2; end
if(n2==1 && n1>1), bb2=repmat(bb2,n1,1); n2=n1; end
assert(n1==n2);
lcsE=max(bb1(:,1:2)+bb1(:,3:4),bb2(:,1:2)+bb2(:,3:4));
lcsS=min(bb1(:,1:2),bb2(:,1:2));
bb=[lcsS lcsE-lcsS];
end

function bb = resize( bb, hr, wr, ar )
% Resize the bbs (without moving their centers).
%
% If wr>0 or hr>0, the w/h of each bb is adjusted in the following order:
%  if(hr~=0), h=h*hr; end
%  if(wr~=0), w=w*wr; end
%  if(hr==0), h=w/ar; end
%  if(wr==0), w=h*ar; end
% Only one of hr/wr may be set to 0, and then only if ar>0. If, however,
% hr=wr=0 and ar>0 then resizes bbs such that areas and centers are
% preserved but aspect ratio becomes ar.
%
% USAGE
%  bb = bbApply( 'resize', bb, hr, wr, [ar] )
%
% INPUTS
%  bb     - [nx4] original bbs
%  hr     - ratio by which to multiply height (or 0)
%  wr     - ratio by which to multiply width (or 0)
%  ar     - [0] target aspect ratio (used only if hr=0 or wr=0)
%
% OUTPUT
%  bb    - [nx4] the output resized bbs
%
% EXAMPLE
%  bb = bbApply('resize',[0 0 1 1],1.2,0,.5) % h'=1.2*h; w'=h'/2;
%
% See also bbApply, bbApply>squarify
if(nargin<4), ar=0; end; assert(size(bb,2)>=4);
assert((hr>0&&wr>0)||ar>0);
% preserve area and center, set aspect ratio
if(hr==0 && wr==0), a=sqrt(bb(:,3).*bb(:,4)); ar=sqrt(ar);
  d=a*ar-bb(:,3); bb(:,1)=bb(:,1)-d/2; bb(:,3)=bb(:,3)+d;
  d=a/ar-bb(:,4); bb(:,2)=bb(:,2)-d/2; bb(:,4)=bb(:,4)+d; return;
end
% possibly adjust h/w based on hr/wr
if(hr~=0), d=(hr-1)*bb(:,4); bb(:,2)=bb(:,2)-d/2; bb(:,4)=bb(:,4)+d; end
if(wr~=0), d=(wr-1)*bb(:,3); bb(:,1)=bb(:,1)-d/2; bb(:,3)=bb(:,3)+d; end
% possibly adjust h/w based on ar and NEW h/w
if(~hr), d=bb(:,3)/ar-bb(:,4); bb(:,2)=bb(:,2)-d/2; bb(:,4)=bb(:,4)+d; end
if(~wr), d=bb(:,4)*ar-bb(:,3); bb(:,1)=bb(:,1)-d/2; bb(:,3)=bb(:,3)+d; end
end

function bbr = squarify( bb, flag, ar )
% Fix bb aspect ratios (without moving the bb centers).
%
% The w or h of each bb is adjusted so that w/h=ar.
% The parameter flag controls whether w or h should change:
%  flag==0: expand bb to given ar
%  flag==1: shrink bb to given ar
%  flag==2: use original w, alter h
%  flag==3: use original h, alter w
%  flag==4: preserve area, alter w and h
% If ar==1 (the default), always converts bb to a square, hence the name.
%
% USAGE
%  bbr = bbApply( 'squarify', bb, flag, [ar] )
%
% INPUTS
%  bb     - [nx4] original bbs
%  flag   - controls whether w or h should change
%  ar     - [1] desired aspect ratio
%
% OUTPUT
%  bbr    - the output 'squarified' bbs
%
% EXAMPLE
%  bbr = bbApply('squarify',[0 0 1 2],0)
%
% See also bbApply, bbApply>resize
if(nargin<3 || isempty(ar)), ar=1; end; bbr=bb;
if(flag==4), bbr=resize(bb,0,0,ar); return; end
for i=1:size(bb,1), p=bb(i,1:4);
  usew = (flag==0 && p(3)>p(4)*ar) || (flag==1 && p(3)<p(4)*ar) || flag==2;
  if(usew), p=resize(p,0,1,ar); else p=resize(p,1,0,ar); end; bbr(i,1:4)=p;
end
end

function hs = draw( bb, col, lw, ls, prop, ids )
% Draw single or multiple bbs to image (calls rectangle()).
%
% To draw bbs aligned with pixel boundaries, subtract .5 from the x and y
% coordinates (since pixel centers are located at integer locations).
%
% USAGE
%  hs = bbApply( 'draw', bb, [col], [lw], [ls], [prop], [ids] )
%
% INPUTS
%  bb     - [nx4] standard bbs or [nx5] weighted bbs
%  col    - ['g'] color or [kx1] array of colors
%  lw     - [2] LineWidth for rectangle
%  ls     - ['-'] LineStyle for rectangle
%  prop   - [] other properties for rectangle
%  ids    - [ones(1,n)] id in [1,k] for each bb into colors array
%
% OUTPUT
%  hs     - [nx1] handles to drawn rectangles (and labels)
%
% EXAMPLE
%  im(rand(3)); bbApply('draw',[1.5 1.5 1 1 .5],'g');
%
% See also bbApply, bbApply>embed, rectangle
[n,m]=size(bb); if(n==0), hs=[]; return; end
if(nargin<2 || isempty(col)), col=[]; end
if(nargin<3 || isempty(lw)), lw=2; end
if(nargin<4 || isempty(ls)), ls='-'; end
if(nargin<5 || isempty(prop)), prop={}; end
if(nargin<6 || isempty(ids)), ids=ones(1,n); end
% prepare display properties
prop=['LineWidth' lw 'LineStyle' ls prop 'EdgeColor'];
tProp={'FontSize',10,'color','w','FontWeight','bold',...
  'VerticalAlignment','bottom'}; k=max(ids);
if(isempty(col)), if(k==1), col='g'; else col=hsv(k); end; end
if(size(col,1)<k), ids=ones(1,n); end; hs=zeros(1,n);
% draw rectangles and optionally labels
for b=1:n, hs(b)=rectangle('Position',bb(b,1:4),prop{:},col(ids(b),:)); end
if(m==4), return; end; hs=[hs zeros(1,n)]; bb=double(bb);
for b=1:n, hs(b+n)=text(bb(b,1),bb(b,2),num2str(bb(b,5),4),tProp{:}); end
end

function I = embed( I, bb, varargin )
% Embed single or multiple bbs directly into image.
%
% USAGE
%  I = bbApply( 'embed', I, bb, varargin )
%
% INPUTS
%  I      - input image
%  bb     - [nx4] or [nx5] input bbs
%  varargin   - additional params (struct or name/value pairs)
%    .col       - [0 255 0] color for rectangle or nx3 array of colors
%    .lw        - [3] width for rectangle in pixels
%    .fh        - [35] font height (if displaying weight), may be 0
%    .fcol      - [255 0 0] font color or nx3 array of colors
%
% OUTPUT
%  I      - output image
%
% EXAMPLE
%  I=imResample(imread('cameraman.tif'),2); bb=[200 70 70 90 0.25];
%  J=bbApply('embed',I,bb,'col',[0 0 255],'lw',8,'fh',30); figure(1); im(J)
%  K=bbApply('embed',J,bb,'col',[0 255 0],'lw',2,'fh',30); figure(2); im(K)
%
% See also bbApply, bbApply>draw, char2img

% get additional parameters
dfs={'col',[0 255 0],'lw',3,'fh',35,'fcol',[255 0 0]};
[col,lw,fh,fcol]=getPrmDflt(varargin,dfs,1);
n=size(bb,1); bb(:,1:4)=round(bb(:,1:4));
if(size(col,1)==1), col=col(ones(1,n),:); end
if(size(fcol,1)==1), fcol=fcol(ones(1,n),:); end
if( ismatrix(I) ), I=I(:,:,[1 1 1]); end
% embed each bb
x0=bb(:,1); x1=x0+bb(:,3)-1; y0=bb(:,2); y1=y0+bb(:,4)-1;
j0=floor((lw-1)/2); j1=ceil((lw-1)/2); h=size(I,1); w=size(I,2);
x00=max(1,x0-j0); x01=min(x0+j1,w); x10=max(1,x1-j0); x11=min(x1+j1,w);
y00=max(1,y0-j0); y01=min(y0+j1,h); y10=max(1,y1-j0); y11=min(y1+j1,h);
for b=1:n
  for c=1:3, I([y00(b):y01(b) y10(b):y11(b)],x00(b):x11(b),c)=col(b,c); end
  for c=1:3, I(y00(b):y11(b),[x00(b):x01(b) x10(b):x11(b)],c)=col(b,c); end
end
% embed text displaying bb score (inside upper-left bb corner)
if(size(bb,2)<5 || fh==0), return; end
bb(:,1:4)=intersect(bb(:,1:4),[1 1 w h]);
for b=1:n
  M=char2img(sprintf('%.4g',bb(b,5)),fh); M=M{1}==0; [h,w]=size(M);
  y0=bb(b,2); y1=y0+h-1; x0=bb(b,1); x1=x0+w-1;
  if( x0>=1 && y0>=1 && x1<=size(I,2) && y1<=size(I,1))
    Ir=I(y0:y1,x0:x1,1); Ig=I(y0:y1,x0:x1,2); Ib=I(y0:y1,x0:x1,3);
    Ir(M)=fcol(b,1); Ig(M)=fcol(b,2); Ib(M)=fcol(b,3);
    I(y0:y1,x0:x1,:)=cat(3,Ir,Ig,Ib);
  end
end
end

function [patches, bbs] = crop( I, bbs, padEl, dims )
% Crop image regions from I encompassed by bbs.
%
% The only subtlety is that a pixel centered at location (i,j) would have a
% bb of [j-1/2,i-1/2,1,1].  The -1/2 is because pixels are located at
% integer locations. This is a Matlab convention, to confirm use:
%  im(rand(3)); bbApply('draw',[1.5 1.5 1 1],'g')
% If bb contains all integer entries cropping is straightforward. If
% entries are not integers, x=round(x+.499) is used, eg 1.2 actually goes
% to 2 (since it is closer to 1.5 then .5), and likewise for y.
%
% If ~isempty(padEl), image is padded so can extract full bb region (no
% actual padding is done, this is fast). Otherwise bb is intersected with
% the image bb prior to cropping. If padEl is a string ('circular',
% 'replicate', or 'symmetric'), uses padarray to do actual padding (slow).
%
% USAGE
%  [patches, bbs] = bbApply('crop',I,bb,[padEl],[dims])
%
% INPUTS
%  I        - image from which to crop patches
%  bbs      - bbs that indicate regions to crop
%  padEl    - [0] value to pad I or [] to indicate no padding (see above)
%  dims     - [] if specified resize each cropped patch to [w h]
%
% OUTPUTS
%  patches  - [1xn] cell of cropped image regions
%  bbs      - actual integer-valued bbs used to crop
%
% EXAMPLE
%  I=imread('cameraman.tif'); bb=[-10 -10 100 100];
%  p1=bbApply('crop',I,bb); p2=bbApply('crop',I,bb,'replicate');
%  figure(1); im(I); figure(2); im(p1{1}); figure(3); im(p2{1});
%
% See also bbApply, ARRAYCROP, PADARRAY, IMRESAMPLE

% get padEl, bound bb to visible region if empty
if( nargin<3 ), padEl=0; end; h=size(I,1); w=size(I,2);
if( nargin<4 ), dims=[]; end;
if(isempty(padEl)), bbs=intersect([.5 .5 w h],bbs); end
% crop each patch in turn
n=size(bbs,1); patches=cell(1,n);
for i=1:n, [patches{i},bbs(i,1:4)]=crop1(bbs(i,1:4)); end

  function [patch, bb] = crop1( bb )
    % crop single patch (use arrayCrop only if necessary)
    lcsS=round(bb([2 1])+.5-.001); lcsE=lcsS+round(bb([4 3]))-1;
    if( any(lcsS<1) || lcsE(1)>h || lcsE(2)>w )
      if( ischar(padEl) )
        pt=max(0,1-lcsS(1)); pb=max(0,lcsE(1)-h);
        pl=max(0,1-lcsS(2)); pr=max(0,lcsE(2)-w);
        lcsS1=max(1,lcsS); lcsE1=min(lcsE,[h w]);
        patch = I(lcsS1(1):lcsE1(1),lcsS1(2):lcsE1(2),:);
        patch = padarray(patch,[pt pl],padEl,'pre');
        patch = padarray(patch,[pb pr],padEl,'post');
      else
        if(ndims(I)==3); lcsS=[lcsS 1]; lcsE=[lcsE 3]; end
        patch = arrayCrop(I,lcsS,lcsE,padEl);
      end
    else
      patch = I(lcsS(1):lcsE(1),lcsS(2):lcsE(2),:);
    end
    bb = [lcsS([2 1]) lcsE([2 1])-lcsS([2 1])+1];
    if(~isempty(dims)), patch=imResample(patch,[dims(2),dims(1)]); end
  end
end

function bb = convert( bb, bbRef, isAbs )
% Convert bb relative to absolute coordinates and vice-versa.
%
% If isAbs==1, bb is assumed to be given in absolute coords, and the output
% is given in coords relative to bbRef. Otherwise, if isAbs==0, bb is
% assumed to be given in coords relative to bbRef and the output is given
% in absolute coords.
%
% USAGE
%  bb = bbApply( 'convert', bb, bbRef, isAbs )
%
% INPUTS
%  bb     - original bb, either in abs or rel coords
%  bbRef  - reference bb
%  isAbs  - 1: bb is in abs coords, 0: bb is in rel coords
%
% OUTPUTS
%  bb     - converted bb
%
% EXAMPLE
%  bbRef=[5 5 15 15]; bba=[10 10 5 5];
%  bbr = bbApply( 'convert', bba, bbRef, 1 )
%  bba2 = bbApply( 'convert', bbr, bbRef, 0 )
%
% See also bbApply
if( isAbs )
  bb(1:2)=bb(1:2)-bbRef(1:2);
  bb=bb./bbRef([3 4 3 4]);
else
  bb=bb.*bbRef([3 4 3 4]);
  bb(1:2)=bb(1:2)+bbRef(1:2);
end
end

function bbs = random( varargin )
% Randomly generate bbs that fall in a specified region.
%
% The vector dims defines the region in which bbs are generated. Specify
% dims=[height width] to generate bbs=[x y w h] such that: 1<=x<=width,
% 1<=y<=height, x+w-1<=width, y+h-1<=height. The biggest bb generated can
% be bb=[1 1 width height]. If dims is a three element vector the third
% coordinate is the depth, in this case bbs=[x y w h d] where 1<=d<=depth.
%
% A number of constraints can be specified that control the size and other
% characteristics of the generated bbs. Note that if incompatible
% constraints are specified (e.g. if the maximum width and height are both
% 5 while the minimum area is 100) no bbs will be generated. More
% generally, if fewer than n bbs are generated a warning is displayed.
%
% USAGE
%  bbs =  bbApply( 'random', pRandom )
%
% INPUTS
%  pRandom    - parameters (struct or name/value pairs)
%   .n          - ['REQ'] number of bbs to generate
%   .dims       - ['REQ'] region in which to generate bbs [height,width]
%   .wRng       - [1 inf] range for width of bbs (or scalar value)
%   .hRng       - [1 inf] range for height of bbs (or scalar value)
%   .aRng       - [1 inf] range for area of bbs
%   .arRng      - [0 inf] range for aspect ratio (width/height) of bbs
%   .unique     - [1] if true generate unique bbs
%   .maxOverlap - [1] max overlap (intersection/union) between bbs
%   .maxIter    - [100] max iterations to go w/o changes before giving up
%   .show       - [0] if true show sample generated bbs
%
% OUTPUTS
%  bbs        - [nx4] array of randomly generated integer bbs
%
% EXAMPLE
%  bbs=bbApply('random','n',50,'dims',[20 20],'arRng',[.5 .5],'show',1);
%
% See also bbApply

% get parameters
rng=[1 inf]; dfs={ 'n','REQ', 'dims','REQ', 'wRng',rng, 'hRng',rng, ...
  'aRng',rng, 'arRng',[0 inf], 'unique',1, 'maxOverlap',1, ...
  'maxIter',100, 'show',0 };
[n,dims,wRng,hRng,aRng,arRng,uniqueOnly,maxOverlap,maxIter,show] ...
  = getPrmDflt(varargin,dfs,1);
if(length(hRng)==1), hRng=[hRng hRng]; end
if(length(wRng)==1), wRng=[wRng wRng]; end
if(length(dims)==3), d=5; else d=4; end

% generate random bbs satisfying constraints
bbs=zeros(0,d); ids=zeros(0,1); n1=min(n*10,1000);
M=max(dims)+1; M=M.^(0:d-1); iter=0; k=0;
tid=ticStatus('generating random bbs',1,2);
while( k<n && iter<maxIter )
  ys=1+floor(rand(2,n1)*dims(1)); ys0=min(ys); ys1=max(ys); hs=ys1-ys0+1;
  xs=1+floor(rand(2,n1)*dims(2)); xs0=min(xs); xs1=max(xs); ws=xs1-xs0+1;
  if(d==5), ds=1+floor(rand(1,n1)*dims(3)); else ds=zeros(0,n1); end
  if(arRng(1)==arRng(2)), ws=hs.*arRng(1); end
  ars=ws./hs; ws=round(ws); xs1=xs0+ws-1; as=ws.*hs;
  kp = ys0>0 & xs0>0 & ys1<=dims(1) & xs1<=dims(2) & ...
    hs>=hRng(1) & hs<=hRng(2) & ws>=wRng(1) & ws<=wRng(2) & ...
    as>=aRng(1) & as<=aRng(2) & ars>=arRng(1) & ars<=arRng(2);
  bbs1=[xs0' ys0' ws' hs' ds']; bbs1=bbs1(kp,:);
  k0=k; bbs=[bbs; bbs1]; k=size(bbs,1); %#ok<AGROW>
  if( maxOverlap<1 && k ), bbs=bbs(1:k0,:);
    for j=1:size(bbs1,1), bbs0=bbs; bb=bbs1(j,:);
      if(d==5), bbs=bbs(bbs(:,5)==bb(5),:); end
      if(isempty(bbs)), bbs=[bbs0; bb]; continue; end
      ws1=min(bbs(:,1)+bbs(:,3),bb(1)+bb(3))-max(bbs(:,1),bb(1));
      hs1=min(bbs(:,2)+bbs(:,4),bb(2)+bb(4))-max(bbs(:,2),bb(2));
      o=max(0,ws1).*max(0,hs1); o=o./(bbs(:,3).*bbs(:,4)+bb(3).*bb(4)-o);
      if(max(o)<=maxOverlap), bbs=[bbs0; bb]; else bbs=bbs0; end
    end
  elseif( uniqueOnly && k )
    ids=[ids; sum(bbs1.*M(ones(1,size(bbs1,1)),:),2)]; %#ok<AGROW>
    [ids,o]=sort(ids); bbs=bbs(o,:); kp=[ids(1:end-1)~=ids(2:end); true];
    bbs=bbs(kp,:); ids=ids(kp,:);
  end
  k=size(bbs,1); if(k0==k), iter=iter+1; else iter=0; end
  if(k>n), bbs=bbs(randSample(k,n),:); k=n; end;
  tocStatus(tid,max(k/n,iter/maxIter));
end
if( k<n ), warning('only generated %i of %i bbs',k,n); n=k; end %#ok<WNTAG>

% optionally display a few bbs
if( show )
  k=8; figure(show); im(zeros(dims)); cs=uniqueColors(1,k,0,0);
  if(n>k), bbs1=bbs(randsample(n,k),:); else bbs1=bbs; end
  bbs1(:,1:2)=bbs1(:,1:2)-.5;
  for i=1:min(k,n), rectangle('Position',bbs1(i,:),...
      'EdgeColor',cs(i,:),'LineStyle','--'); end
end

end

function bbs = frMask( M, bbw, bbh, thr )
% Convert weighted mask to bbs.
%
% Pixels in mask above given threshold (thr) indicate bb centers.
%
% USAGE
%  bbs = bbApply('frMask',M,bbw,bbh,[thr])
%
% INPUTS
%  M      - mask
%  bbw    - bb target width
%  bbh    - bb target height
%  thr    - [0] mask threshold
%
% OUTPUTS
%  bbs    - bounding boxes
%
% EXAMPLE
%  w=20; h=10; bbw=5; bbh=8; M=double(rand(h,w)); M(M<.95)=0;
%  bbs=bbApply('frMask',M,bbw,bbh); M2=bbApply('toMask',bbs,w,h);
%  sum(abs(M(:)-M2(:)))
%
% See also bbApply, bbApply>toMask
if(nargin<4), thr=0; end
ids=find(M>thr); ids=ids(:); h=size(M,1);
if(isempty(ids)), bbs=zeros(0,5); return; end
xs=floor((ids-1)/h); ys=ids-xs*h; xs=xs+1;
bbs=[xs-floor(bbw/2) ys-floor(bbh/2)];
bbs(:,3)=bbw; bbs(:,4)=bbh; bbs(:,5)=M(ids);
end

function M = toMask( bbs, w, h, fill, bgrd )
% Create weighted mask encoding bb centers (or extent).
%
% USAGE
%  M = bbApply('toMask',bbs,w,h,[fill],[bgrd])
%
% INPUTS
%  bbs    - bounding boxes
%  w      - mask target width
%  h      - mask target height
%  fill   - [0] if 1 encodes extent of bbs
%  bgrd   - [0] default value for background pixels
%
% OUTPUTS
%  M      - hxw mask
%
% EXAMPLE
%
% See also bbApply, bbApply>frMask
if(nargin<4||isempty(fill)), fill=0; end
if(nargin<5||isempty(bgrd)), bgrd=0; end
if(size(bbs,2)==4), bbs(:,5)=1; end
M=zeros(h,w); B=true(h,w); n=size(bbs,1);
if( fill==0 )
  p=floor(getCenter(bbs)); p=sub2ind([h w],p(:,2),p(:,1));
  for i=1:n, M(p(i))=M(p(i))+bbs(i,5); end
  if(bgrd~=0), B(p)=0; end
else
  bbs=[intersect(round(bbs),[1 1 w h]) bbs(:,5)]; n=size(bbs,1);
  x0=bbs(:,1); x1=x0+bbs(:,3)-1; y0=bbs(:,2); y1=y0+bbs(:,4)-1;
  for i=1:n, y=y0(i):y1(i); x=x0(i):x1(i);
    M(y,x)=M(y,x)+bbs(i,5); B(y,x)=0; end
end
if(bgrd~=0), M(B)=bgrd; end
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function varargout = implies(varargin)
%IMPLIES Logical implication
%
% IMPLIES(X,Y) creates a mixed integer representation of
% the constraint X --> Y, i.e. Y is true if X is true.
%
% Syntax
%   F = implies(X,Y,tol)
%
% Input
%   X : binary SDPVAR variable or a constraint
%   Y : binary SDPVAR variable or a constraint
%  tol: Optional threshhold for defining zero (see NOTE)
%
% Output
%   F : SET object
%
% Examples
%
%  binvar X Y; F = implies(X,Y);
%  binvar X;sdpvar Y; F = [implies(X,Y>=5), -10 <= Y <= 100];
%  binvar X;Y=sdpvar(3,1); F = [implies(X,[sum(Y);Y(2)]>=[5;0]), -1<= Y <= 10];
%
% Note
%
%  All variables in the expressions have to be explicitly bounded somewhere
%  in the model (implicit constraints are not sufficients such as [Y <= Z,...,Z<= 10]
%
%  Using implies with X non-binary is highly sensitive numerically.
%  The problem comes from the definition of 0 in a floating-point
%  environment, and precision in the solver. To account for this,
%  the user can supply a third argument to define a dead-zone around
%  zero, i.e Implies(X<=0,Y) will be replaced with IMPLIES(X<=-tol,Y)
%  Note, you typically need to tweak this number for your
%  application/solver. By default, YALMIP uses tol = 0, which means you
%  easily can get garbage... A positive number means YALMIP is cautious in
%  terms of activating the condition, while a negative number means YALMIP
%  will be aggressive  in activating the condition.
%
%   See also @SDPVAR/AND, @SDPVAR/OR, IFF

% There are some cases to take care of...
%
% X --> Y     binary/binary :                     Implemented
% X --> Y     binary/(LP,equality,sdp)            Implemented
% X --> Y     (LP,equality,sdp)/binary            Not implemented
% X --> Y     (LP,equality,sdp)/(LP,equality,sdp) Not implemented

X = varargin{1};
Y = varargin{2};

if isempty(X)
    varargout{1} = [];
end

switch class(X)

    case {'sdpvar','constraint','lmi'}

        if nargin > 2
            if ~isa(varargin{3},'double')
                error('Third argument in IMPLIES should be a numerical value (tolerance)');
            end
        end
        varargout{1} = setupMeta(lmi([]), mfilename,varargin{:});

    case 'char'
        varargout{1} = implies_internal(varargin{3:end});

    case 'logical'
        if length(X)==1
            if X
                varargout{1} = Y;
            else
                varargout{1} = [];
            end
        else
            if length(X) == length(Y)
                i = find(X);
                if isempty(i)
                    varargout{1} = [];
                else
                    varargout{1} = Y(i);
                end
            else
                error('Size mismatch in input arguments');
            end
        end
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function tv = Fval(u, img, alpha, huber)
% total variation part of the criterion

[H W] = size(img);
N = W * H;

nabla = make_derivatives_mine(H, W);

if huber
    tv = nabla * u(:);
    tv_size = sqrt(tv(1:N).^2 + tv(N+1:end).^2);
    idx1 = tv_size <= alpha;
    idx2 = tv_size > alpha;
    tv = sum(tv_size(idx1).^2/(2*alpha)) + sum(tv_size(idx2) - alpha/2);
else
    tv = nabla * u(:);
    tv = sum(sqrt(tv(1:N).^2 + tv(N+1:end).^2));
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function wathen_test12 ( )

%*****************************************************************************80
%
%% WATHEN_TEST12 assemble, factor and solve using WATHEN_DAVIS + CG_SPARSE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 June 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'WATHEN_TEST12\n' );
  fprintf ( 1, '  Assemble, factor and solve a Wathen system\n' );
  fprintf ( 1, '  defined by WATHEN_DAVIS and CG_SPARSE.\n' );
  fprintf ( 1, '\n' );

  nx = 1;
  ny = 1;
  fprintf ( 1, '  Elements in X direction NX = %d\n', nx );
  fprintf ( 1, '  Elements in Y direction NY = %d\n', ny );
  fprintf ( 1, '  Number of elements = %d\n', nx * ny );
%
%  Compute the number of unknowns.
%
  n = wathen_order ( nx, ny );
  fprintf ( 1, '  Number of nodes N = %d\n', n );
%
%  Set up a random solution X1.
%
  seed = 123456789;
  [ x1, seed ] = r8vec_uniform_01 ( n, seed );
%
%  Compute the matrix.
%
  seed = 123456789;
  [ a, seed ] = wathen_davis ( nx, ny, n, seed );
%
%  Compute the corresponding right hand side B.
%
  b = a * x1;
%
%  Solve the linear system.
%
  x2 = ones ( n, 1 );
  x2 = cg_sparse ( n, a, b, x2 );
%
%  Compute the maximum solution error.
%
  e = max ( abs ( x1 - x2 ) );
  fprintf ( 1, '  Maximum solution error is %g\n', e );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function pass = test_subsref( pref )

if ( nargin < 1 )
    pref = chebfunpref;
end

tol = 1000*pref.cheb2Prefs.chebfun2eps;

% Evaluation at x,y,z
ff = @(x,y,z) z.*sin(x.*(y-.1));
f = spherefun(ff);
x = 1/sqrt(3); y = 1/sqrt(3); z = -1/sqrt(3);
pass(1) = ( abs( f(x,y,z) - ff(x,y,z) ) < tol );

% Evaluation at (lambda,theta)
ff = @(lam,th) exp(cos(lam).*sin(th).*cos(th));
f = spherefun(ff);
lam = pi/3; th = pi/4;
pass(2) = ( abs( f(lam,th ) - ff(lam,th) ) < tol );

% Slices in theta
ff = @(lam,th) exp(cos(lam).*sin(th).*cos(th));
f = spherefun(ff);
th1 = 0.2;
slice1 = chebfun(@(lam) feval(ff,lam,th1),[-pi pi], 'trig');
th2 = 0.7;
slice2 = chebfun(@(lam) feval(ff,lam,th2),[-pi pi], 'trig');
pass(3) = ( norm( f(:,[th1 th2]).' - [slice1 slice2] ) < tol );

% Slices in lambda
ff = @(lam,th) exp(cos(lam).*sin(th).*cos(th));
f = spherefun(ff);
lam1 = 0.2;
slice1 = chebfun(@(th) feval(ff,lam1,th),[-pi pi], 'trig');
lam2 = 0.7;
slice2 = chebfun(@(th) feval(ff,lam2,th),[-pi pi], 'trig');
pass(4) = ( norm( f([th1 th2],:) - [slice1 slice2] ) < tol );

% Slice in z
ff = @(x,y,z) z.*sin(x.*(y-.1));
f = spherefun(ff);
z = 0.1; th = acos(z);
slice = chebfun(@(lam) feval(ff,cos(lam).*sin(th),sin(lam).*sin(th),z),[-pi pi], 'trig');
pass(5) = ( norm( f(:,:,z) - slice ) < tol );

% Slice in x
ff = @(x,y,z) z.*sin(x.*(y-.1));
f = spherefun(ff);
x = 0.1;
slice = chebfun(@(t) feval(ff,x,sqrt(1-x.^2).*cos(t),sqrt(1-x.^2).*sin(t)),[-pi pi], 'trig');
pass(6) = ( norm( f(x,:,:) - slice ) < tol );

% Slice in y
ff = @(x,y,z) z.*sin(x.*(y-.1));
f = spherefun(ff);
y = -0.4;
slice = chebfun(@(t) feval(ff,sqrt(1-y.^2).*cos(t),y,sqrt(1-y.^2).*sin(t)),[-pi pi], 'trig');
pass(7) = ( norm( f(:,y,:) - slice ) < tol );

% GET properties
f = spherefun(@(lam,th) cos(lam).*sin(th));
pass(8) = ( norm(f.rows - chebfun(@(lam) cos(lam),[-pi pi],'trig')) < tol || ...
    norm(f.rows + chebfun(@(lam) cos(lam),[-pi pi],'trig')) < tol );
pass(9) = ( norm(f.cols - chebfun(@(th) sin(th),[-pi pi],'trig')) < tol || ...
    norm(f.cols + chebfun(@(th) sin(th),[-pi pi],'trig')) < tol );
pass(10) = ( norm( f.domain - [-pi pi 0 pi] ) < tol );

% Composition of a spherefun with a chebfun (1 and 3 columns):
f = spherefun(@(x,y,z) z + sin(pi*x.*y));
g = chebfun(@(t) t.^2, [ -1.5, 1.5 ]);
h_true = spherefun(@(x,y,z) (z + sin(pi*x.*y)).^2);
h = g(f);
pass(11) = ( norm(h - h_true) < tol );

G = chebfun(@(t) [ t.^2, t, -t.^2 ], [ -1.5, 1.5 ]);
H_true = spherefunv(@(x,y,z) (z + sin(pi*x.*y)).^2, ...
    @(x,y,z) z + sin(pi*x.*y), @(x,y,z) -(z + sin(pi*x.*y)).^2);
H = G(f);
pass(12) = ( norm(H - H_true) < tol );

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pass = test_isPeriodicTech(pref)
% Test isPeriodicTech().

if ( nargin == 0 )
    pref = chebfunpref;
end
tol = 1000*pref.cheb3Prefs.chebfun3eps;

f = chebfun3v(@(x,y,z) x, @(x,y,z) y);
pass(1) = ~isPeriodicTech(f);

f1 = chebfun3(@(x,y,z) cos(pi*x), 'trig');
f2 = chebfun3(@(x,y,z) sin(pi*y), 'trig');
pass(2) = isPeriodicTech([f1; f2]);
pass(3) = isPeriodicTech([f1; f2; f2]);
F = [f1; f2];
pass(4) = isPeriodicTech([F; f2]);

end


text:
function traffic ( cycle_num )

%*****************************************************************************80
%
%% TRAFFIC simulates the cars waiting at one traffic light.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2009
%
%  Author:
%
%    Original MATLAB version by Brian Hahn, Dan Valentine.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Brian Hahn, Dan Valentine,
%    Essential MATLAB for Engineers and Scientists,
%    Academic Press, 2009,
%    ISBN13: 978-0123748836,
%    LC: TA345.V34.
%
%  Parameters:
%
%    Input, integer CYCLE_NUM, the number of 10-second time cycles to model.
%
%  Local Parameters:
%
%    Local, integer CARS, the number of cars waiting at the light.
%
%    Local, integer CARS_IN, the total number of cars that have come.
%
%    Local, integer CARS_OUT, the total number of cars that have left.
%
%    Local, integer CYCLE, the number of time cycles that have elapsed.
%
%    Local, integer CYCLE_LENGTH, the number of seconds in one time cycle.
%
%    Local, integer GREEN_CYCLES, the number of 10-second time cycles that
%    a green light lasts.
%
%    Local, integer GREEN_TIMER, keeps track of the number of time cycles the
%    green light has been on.
%
%    Local, integer LIGHT, the state of the light.
%    'r', the light is now red.
%    'g', the light is now green.
%
%    Local, real P, the probability that a new car will come to the light
%    in the next second.
%
%    Local, integer RED_CYCLES, the number of 10-second time cycles that
%    a red light lasts.
%
%    Local, integer RED_TIMER, keeps track of the number of time cycles the
%    red light has been on.
%

%
%  Initialize.
%
  cars = 0;
  cars_in = 0;
  cars_out = 0;
  car_wait_cycles = 0;
  cycle = 0;
  cycle_length = 10;
  green_cycles = 2;
  green_timer = 0;
  light = 'r';
  p = 0.3;
  red_cycles = 4;
  red_timer = 0;
%
%  Set up the plot data.
%
  plot_data = zeros(2,cycle_num+1);
%
%  Handle the "0"-th cycle.
%
  plot_data(1,cycle+1) = cycle;
  plot_data(2,cycle+1) = cars;

  prq ( cars, light, cycle );
%
%  Handle cycles 1 through CYCLE_NUM.
%
  for cycle = 1 : cycle_num
%
%  Each second of the cycle, choose a random number.
%  If it is less than P, then a new car appeared at the light at that second.
%
    r = rand ( cycle_length, 1 );
    cars_new = sum ( r < p );
    cars = cars + cars_new;
    cars_in = cars_in + cars_new;
%
%  Handle this time cycle depending on whether the light is green or red.
%
    if ( light == 'g' )
      [ cars, cars_out, light, green_timer ] = go ( green_cycles, cars, ...
        cars_out, light, green_timer );
    else
      [ cars, light, red_timer ] = stop ( red_cycles, cars, light, red_timer );
    end
%
%  At the end of this cycle, how many cars are waiting?
%
    car_wait_cycles = car_wait_cycles + cars;
%
%  Print the current status.
%
    prq ( cars, light, cycle );

    plot_data(1,cycle+1) = cycle;
    plot_data(2,cycle+1) = cars;

  end

  plot ( plot_data(1,1:cycle_num+1), plot_data(2,1:cycle_num+1) )
  xlabel ( 'Time Cycles' )
  ylabel ( 'Cars Waiting' )
  title ( 'Traffic waiting at a Light' )

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Number of cycles =       %d\n', cycle_num );
  fprintf ( 1, '  Simulated time =         %d seconds\n', cycle_num * cycle_length );
  fprintf ( 1, '  Number of cars in =      %d\n', cars_in );
  fprintf ( 1, '  Number of cars waiting = %d\n', cars );
  fprintf ( 1, '  Number of cars out =     %d\n', cars_out );
  fprintf ( 1, '  Percentage Out/In = %7.1f%%\n', 100 * cars_out / cars_in );
  wait_average_seconds = car_wait_cycles * cycle_length / cars_in;
  fprintf ( 1, '  Average wait = %7.2f seconds\n', wait_average_seconds );
  wait_average_lights = car_wait_cycles / cars_in / ( red_cycles + green_cycles );
  fprintf ( 1, '  Average wait = %7.2f light cycles\n', wait_average_lights );

  return
end
function [ cars, cars_out, light, green_timer ] = go ( green_cycles, cars, ...
  cars_out, light, green_timer )

%*****************************************************************************80
%
%% GO simulates traffic when the light is green.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 November 2009
%
%  Author:
%
%    Original MATLAB version by Brian Hahn, Dan Valentine.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Brian Hahn, Dan Valentine,
%    Essential MATLAB for Engineers and Scientists,
%    Academic Press, 2009,
%    ISBN13: 978-0123748836,
%    LC: TA345.V34.
%
%  Parameters:
%
%    Input, integer GREEN_CYCLES, the number of 10-second time cycles that
%    a green light lasts.
%
%    Input, integer CARS, the number of cars stopped at the light.
%
%    Input, integer CARS_OUT, the total number of cars that have gone
%    through the light.
%
%    Input, integer LIGHT, the state of the light.
%    'r', the light is now red.
%    'g', the light is now green.
%
%    Input, integer GREEN_TIMER, keeps track of the number of time cycles the
%    green light has been on.
%
%    Output, integer CARS, the number of cars stopped at the light.
%
%    Output, integer CARS_OUT, the total number of cars that have gone
%    through the light.
%
%    Output, integer LIGHT, the state of the light.
%    'r', the light is now red.
%    'g', the light is now green.
%
%    Output, integer GREEN_TIMER, keeps track of the number of time cycles the
%    green light has been on.
%

%
%  In one 10-second time cycle, we estimate 8 cars can move out.
%
  cars_through = min ( 8, cars );

  cars = cars - cars_through;
  cars_out = cars_out + cars_through;
%
%  Advance the timer.  If the green light has timed out, reset the timer
%  and switch to red.
%
  green_timer = green_timer + 1;

  if ( green_cycles <= green_timer )
    light = 'r';
    green_timer = 0;
  end

  return
end
function [ cars, light, red_timer ] = stop ( red_cycles, cars, light, ...
  red_timer )

%*****************************************************************************80
%
%% STOP simulates the traffic when the light is red.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2009
%
%  Author:
%
%    Original MATLAB version by Brian Hahn, Dan Valentine.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Brian Hahn, Dan Valentine,
%    Essential MATLAB for Engineers and Scientists,
%    Academic Press, 2009,
%    ISBN13: 978-0123748836,
%    LC: TA345.V34.
%
%  Parameters:
%
%    Input, integer RED_CYCLES, the number of 10-second time cycles that
%    a red light lasts.
%
%    Input, integer CARS, the number of cars stopped at the light.
%
%    Input, integer LIGHT, the state of the light.
%    'r', the light is now red.
%    'g', the light is now green.
%
%    Input, integer RED_TIMER, keeps track of the number of time cycles the
%    red light has been on.
%
%    Output, integer CARS, the number of cars stopped at the light.
%
%    Output, integer LIGHT, the state of the light.
%    'r', the light is now red.
%    'g', the light is now green.
%
%    Output, integer RED_TIMER, keeps track of the number of time cycles the
%    red light has been on.
%

%
%  Advance the timer.
%  If the red light has timed out, reset the timer and switch to green.
%
  red_timer = red_timer + 1;

  if ( red_cycles <= red_timer )
    light = 'g';
    red_timer = 0;
  end

  return
end
function prq ( cars, light, cycle )

%*****************************************************************************80
%
%% PRQ prints the current traffic waiting at the light.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 November 2009
%
%  Author:
%
%    Original MATLAB version by Brian Hahn, Dan Valentine.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Brian Hahn, Dan Valentine,
%    Essential MATLAB for Engineers and Scientists,
%    Academic Press, 2009,
%    ISBN13: 978-0123748836,
%    LC: TA345.V34.
%
%  Parameters:
%
%    Input, integer CARS, the number of cars stopped at the light.
%
%    Input, integer LIGHT, the state of the light.
%    'r', the light is now red.
%    'g', the light is now green.
%
%    Input, integer CYCLE, the current 10-second time cycle.
%
  fprintf ( 1, '%4d ', cycle );
  if ( light == 'r' )
    fprintf ( 'R  ' );
  else
    fprintf ( 'G  ' );
  end
  i = cars;
  c = floor ( i / 100 );
  i = i - 100 * c;
  for j = 1 : c
    fprintf ( 'C' );
  end
  x = floor ( i / 10 );
  i = i - 10 * x;
  for j = 1 : x
    fprintf ( 'X' );
  end
  for j = 1 : i
    fprintf ( 'I' );
  end
  fprintf ( 1, '\n' );

  return
end


text:
% ------------------------------------------------------
% SwarmOps - Heuristic optimization for Matlab
% Copyright (C) 2003-2010 Magnus Erik Hvass Pedersen.
% Please see the file license.txt for license details.
% SwarmOps on the internet: http://www.Hvass-Labs.org/
% ------------------------------------------------------

% Pick a random point from a bounded range.
% Parameters:
%     n; dimensionality of the search-space.
%     x; center position of the sampling-range.
%     range; desired sampling-range.
%     lowerBound; search-space lower-boundary.
%     upperBound; search-space upper-boundary.
% Returns:
%     fitness; the measure to be minimized.
function y = samplebounded(n, x, range, lowerBound, upperBound)

    % Adjust sampling range so it does not exceed bounds.
    l = max(x - range, lowerBound);
    u = min(x + range, upperBound);

    % Pick a sample from within the bounded range.
    y = initagent(n, l, u);
end

% ------------------------------------------------------

author: Sable

repo: mcbench-benchmarks


text:
function euler = qua2euler(qin)
% qua2euler: transforms quaternion to Euler angles.
%
% INPUT
%   qin: 4x1 quaternion.
%
% OUTPUT
%   euler: 3x1 Euler angles [roll pitch yaw] (rad, rad, rad).
%
%   Copyright (C) 2014, Rodrigo Gonzalez, all rights reserved.
%
%   This file is part of NaveGo, an open-source MATLAB toolbox for
%   simulation of integrated navigation systems.
%
%   NaveGo is free software: you can redistribute it and/or modify
%   it under the terms of the GNU Lesser General Public License (LGPL)
%   version 3 as published by the Free Software Foundation.000000000
%
%   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 Lesser General Public License for more details.
%
%   You should have received a copy of the GNU Lesser General Public
%   License along with this program. If not, see
%   <http://www.gnu.org/licenses/>.
%
% References:
%
%		Dr. Paolo Zoccarato's comments at
% https://github.com/rodralez/NaveGo/pull/9
%
%		Crassidis, J.L. and Junkins, J.L. (2011). Optimal Esti-
% mation of Dynamic Systems, 2nd Ed. Chapman and Hall/CRC, USA.
% Eq. 7.39, p. 458.
%
% Version: 005
% Date:    2017/12/05
% Author:  Rodrigo Gonzalez <rodralez@frm.utn.edu.ar>
% URL:     https://github.com/rodralez/navego

% Quaternion format from Crassidis' book.

DCMbn = qua2dcm(qin);

phi   = atan2( DCMbn(3,2), DCMbn(3,3) );    % roll
theta = asin (-DCMbn(3,1) );                % pitch
psi   = atan2( DCMbn(2,1), DCMbn(1,1) );    % yaw

euler = [phi theta psi];

end

author: rodralez

repo: NaveGo

save_path: github-repos/MATLAB/rodralez-NaveGo


text:
function findingVolumeMatch

vS = findVolume('SmoothRectangle');
vR = findVolume('Rectangle');

end

function volV = findVolume(micro)

d = load(['/home/alex/git-repos/SwanLab/Swan/Output/',micro,'/',micro,'.mat']);

var = d.d.variables;

mxV = d.d.domVariables.mxV;
myV = d.d.domVariables.myV;

for imx = 1:length(mxV)
    for imy = 1:length(myV)
        volV(imx,imy) = var{imx,imy}.volume;
    end
end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function  y = At_xy_nonorm(z, Phi)
%z = z./Phi_sum;
y = bsxfun(@times, z, Phi);

end

repo: DeSCI


text:
function [pvec, pstruct] = tapas_cdfgaussian_obs_transp(r, ptrans)
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2012-2015 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% This is a dummy function

pvec    = [];
pstruct = struct;

return;

author: translationalneuromodeling

repo: tapas


text:
% Copyright (C) 2018  Symeon Symeonidis, Stefanos Tsantilas, Stelios Mitilineos
% simos421@gmail.com, steftsantilas@gmail.com, smitil@gmail.com
% 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.

function CstDiscretePort(mws,PortNumber,SetP1,SetP2)

DiscretePort = invoke(mws,'DiscretePort');
invoke(DiscretePort,'Reset');
invoke(DiscretePort,'PortNumber',int2str(PortNumber));
invoke(DiscretePort,'Type','SParameter');
invoke(DiscretePort,'Label','');
invoke(DiscretePort,'Impedance','50');
invoke(DiscretePort,'VoltagePortImpedance','0.0');
invoke(DiscretePort,'Voltage','1.0');
invoke(DiscretePort,'Current','1.0');
invoke(DiscretePort,'SetP1','False',int2str(SetP1(1)),int2str(SetP1(2)),int2str(SetP1(3)));
invoke(DiscretePort,'SetP2','False',int2str(SetP2(1)),int2str(SetP2(2)),int2str(SetP2(3)));
invoke(DiscretePort,'InvertDirection','False');
invoke(DiscretePort,'LocalCoordinates','False');
invoke(DiscretePort,'Monitor','True');
invoke(DiscretePort,'Radius','0.0');
invoke(DiscretePort,'Wire','');
invoke(DiscretePort,'Position','end1');
invoke(DiscretePort,'Create');
release(DiscretePort);
end

repo: CST-MATLAB-API


text:
%  INTERNAL FUNCTION: Numerical evaluation of the hessian of the policy objective
%
%  ::
%
%    H=evaluate_policy_objective_hessian_numerically(funcs,y,x,ss,param,sparam,def,s0,s1)
%
%  Args:
%
%     - **funcs** [fhandle|cell array]: function or functions to be
%       differentiated
%     - **y** [vector]: values of endogenous variables
%     - **x** [vector]: values of exogenous variables
%     - **ss** [vector]: steady state
%     - **param** [vector]: parameter vector
%     - **sparam** [vector]: vector of parameters appearing with a lead
%     - **def** [vector]: values of definitions
%     - **s0** [scalar]: state today
%     - **s1** [scalar]: state tomorrow
%
%  Returns:
%     :
%
%     - **H** [matrix]: Numerical Hessian of **funcs** at [y]
%
%  Note:
%     It is assumed that the inputs are y,x,ss,param,sparam,def,s0,s1 as
%     ordered in parser.input_list()
%
%

author: jmaih


text:
classdef RMMEDA_F5 < PROBLEM
% <multi> <real> <large/none>
% Benchmark MOP for testing RM-MEDA

%------------------------------- Reference --------------------------------
% Q. Zhang, A. Zhou, and Y. Jin, RM-MEDA: A regularity model-based
% multiobjective estimation of distribution algorithm, IEEE Transactions on
% Evolutionary Computation, 2008, 12(1): 41-63.
%------------------------------- 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".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            if isempty(obj.D); obj.D = 30; end
            obj.lower    = zeros(1,obj.D);
            obj.upper    = ones(1,obj.D);
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values
        function PopObj = CalObj(obj,X)
            g = 1 + 9*mean((X(:,2:end).^2-repmat(X(:,1),1,size(X,2)-1)).^2,2);
            PopObj(:,1) = X(:,1);
            PopObj(:,2) = g.*(1-sqrt(PopObj(:,1)./g));
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = linspace(0,1,N)';
            R(:,2) = 1 - sqrt(R(:,1));
        end
        %% Generate the image of Pareto front
        function R = GetPF(obj)
            R = obj.GetOptimum(100);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function model = gpSubspaceOptimise(model,varargin)

% GPSUBSPACEOPTIMISE
%
% COPYRIGHT : Carl Henrik Ek, 2008

% GP

model = gpOptimise(model,varargin{:});

return;

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function varargout=lpcpp2cw(varargin)
%
% For calling details please see v_lpcpp2cw.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcpp2cw.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcpp2cw.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcpp2cw(varargin{:});
else
    v_lpcpp2cw(varargin{:});
end

author: ImperialCollegeLondon

repo: sap-voicebox


text:
% COMP_UGRID_UV
% Compare depth-averaged current from 2 unstructured grid
% models with UGRID conventions (http://bit.ly/cf_ugrid), allowing
% comparison with no model specific code.  This code makes a movie
% if there is more than one time step specified

comp_titl='Comparison of 2D No-Wave runs for Hurricane Ike';

titl{1}='ADCIRC';
uris{1}='http://testbedapps.sura.org/thredds/dodsC/in/und/adcirc/ike/ultralite/lr/vardrag/nowave/3d'
%uris{1}='http://testbedapps.sura.org/threddsdev/dodsC/inundation/ADCIRC/rita/2Dvrwoww'
uname{1}='u-vel';
vname{1}='v-vel';
tname{1}='time';

titl{2}='FVCOM';
uris{2}='http://testbedapps.sura.org/thredds/dodsC/inundation/FVCOM/ike/3Dvrwww';
%uris{2}='http://testbedapps.sura.org/threddsdev/dodsC/inundation/FVCOM/rita/2Dvrwoww';
uname{2}='ua';
vname{2}='va';
tname{2}='time';

% for just one plot, define one time:
jdi=datenum([2008 9 13 03 0 0]); %Ike
%jdi=datenum([2005 9 24 06 0 0]); %Rita

% for a movie define more than one time
%jdi=datenum([2008 9 12 21 0 0]):3/24:datenum([2008 9 13 15 0 0]);
if length(jdi)>1
  movie_name='ike';
end

% scale factor for arrows
vfac=0.03;
% bounding box for figures
ax=[-95.2  -94.4   28.8   29.9]
% color range for figures
cax=[0 3];

if length(jdi)>1
  vidobj=VideoWriter(movie_name,'Uncompressed AVI');
  vidobj.FrameRate = 1;
  open(vidobj);
end

% There is nothing model specific below this point

%  loop through each model to read the grid information
for i=1:length(uris)
  tic
  %initialize dataset object
  nc{i}=ncgeodataset(uris{i});
  %get geovariable object
  uvar{i}=nc{i}.geovariable(uname{i});
  vvar{i}=nc{i}.geovariable(vname{i});
  jd{i}=nc{i}.time(tname{i});
  % get mesh variable name (inference array)
  gvar_name{i}=uvar{i}.attribute('mesh');
  if isempty(gvar_name{i}),
    gtype{i}='structured';
  else
    gtype{i}='unstructured';
  end
  switch gtype{i}
    case 'structured'
      g=uvar{i}.grid(itime,:,:);
      lon{i}=g.lon;
      lat{i}=g.lat;
    case 'unstructured'
      gridvar=nc{i}.variable(gvar_name{i});
      a=textscan(gridvar.attribute('node_coordinates'),'%s %s');
      lon{i}=nc{i}{a{1}}(:);
      lat{i}=nc{i}{a{2}}(:);
      location{i}=uvar{i}.attribute('location');
      % find connectivity array variable from mesh variable
      trivar_name=gridvar.attribute('face_node_connectivity');
      % get connnectivity array data
      tri{i}=nc{i}{trivar_name}(:);
      [m,n]=size(tri{i});
      % check/fix orientation of connectivity array
      if m==3,
        tri{i}=tri{i}.';
      elseif n~=3
        disp('Error:Currently handling triangles only');return
      end
  end
end
% loop through time, reading data at each step from all models
fig=figure(1);
clf;set(gcf,'renderer','zbuffer');set(gcf,'color','white');
for k=1:length(jdi);
  for i=1:length(uris)
    itime=date_index(jd{i},jdi(k));
    subplot(1,2,i);cla
    switch gtype{i}
      case 'structured'
        u=double(squeeze(uvar{i}.data(itime,:,:)));
        v=double(squeeze(vvar{i}.data(itime,:,:)));
        U=complex(u,v);
        pcolorjw(lon{i},lat{i},abs(U));shading interp;colorbar;...
          axis(ax);caxis(cax);
        arrows(lon{i},lat{i},U,0.03,'black');
      case 'unstructured'
        % read data at specified time step for all nodes
        u=double(squeeze(uvar{i}.data(itime,:)));
        v=double(squeeze(vvar{i}.data(itime,:)));
        U=complex(u,v);
        U(U==0)=nan;
        igood=isfinite(U);
        % handle case of 'face' location (e.g. FVCOM)
        if strmatch('face',location{i}),
          % calculate u,v lon,lat positions (faster than reading lonc,latc)
          lon_vel=1/3*(lon{i}(tri{i}(:,1))+lon{i}(tri{i}(:,2))+lon{i}(tri{i}(:,3)));
          lat_vel=1/3*(lat{i}(tri{i}(:,1))+lat{i}(tri{i}(:,2))+lat{i}(tri{i}(:,3)));
          lonp=lon{i}(tri{i});
          latp=lat{i}(tri{i});
          patch(lonp.',latp.',abs(U));view(2);shading flat % shading 'interp' doesn't work with patch
          arrows(lon_vel(igood),lat_vel(igood),U(igood),vfac,'black');
        else
          trisurf(tri{i},lon{i},lat{i},0*lon{i}-5,abs(U));view(2);shading interp;
          arrows(lon{i}(igood),lat{i}(igood),U(igood),vfac,'black');
        end
        colorbar;axis(ax);caxis(cax);...
          fac=cos(mean(lat{1}(:))*pi/180);...
          title(sprintf('%s %s (%s): %sZ',titl{i},uname{i},...
          uvar{i}.attribute('units'),datestr(jd{i}(itime),'yyyy-mm-dd HH:MM')));...
          set (gca, 'DataAspectRatio', [1 fac 1] );...
          set (gca, 'tickdir','out');
        toc

    end
  end
  if(length(jdi)>1)
    F=getframe(fig);
    writeVideo(vidobj,F);
  end
end
if(length(jdi)>1)
  close(vidobj)
end

author: nctoolbox

repo: nctoolbox

save_path: github-repos/MATLAB/nctoolbox-nctoolbox


text:
function [T, Eps] = estimateRigidTransform(x, y)
% ESTIMATERIGIDTRANSFORM
%   [T, EPS] = ESTIMATERIGIDTRANSFORM(X, Y) estimates the rigid transformation
%   that best aligns x with y (in the least-squares sense).
%
%   Reference: "Estimating Rigid Transformations" in
%   "Computer Vision, a modern approach" by Forsyth and Ponce (1993), page 480
%   (page 717(?) of the newer edition)
%
%   Input:
%       X: 3xN, N 3-D points (N>=3)
%       Y: 3xN, N 3-D points (N>=3)
%
%   Output
%       T: the rigid transformation that aligns x and y as:  xh = T * yh
%          (h denotes homogenous coordinates)
%          (corrspondence between points x(:,i) and y(:,i) is assumed)
%
%       EPS: the smallest singular value. The closer this value it is
%          to 0, the better the estimate is. (large values mean that the
%          transform between the two data sets cannot be approximated
%          well with a rigid transform.
%
%   Babak Taati, 2003
%   (revised 2009)

if nargin ~= 2
    error('Requires two input arguments.')
end

if size(x,1)~=3 || size(y,1)~=3
    error('Input point clouds must be a 3xN matrix.');
end

if size(x, 2) ~= size(y,2)
    error('Input point clouds must be of the same size');
end

if size(x,2)<3 || size(y,2)<3
    error('At least 3 point matches are needed');
end

pointCount = length(x); % since x has N=3+ points, length shows the number of points

x_centroid = sum(x,2) / pointCount;
y_centroid = sum(y,2) / pointCount;

x_centrized = [x(1,:)-x_centroid(1) ; x(2,:)-x_centroid(2); x(3,:)-x_centroid(3)];
y_centrized = [y(1,:)-y_centroid(1) ; y(2,:)-y_centroid(2); y(3,:)-y_centroid(3)];

R12 = y_centrized' - x_centrized';
R21 = x_centrized - y_centrized;
R22_1 = y_centrized  + x_centrized;
R22 = crossTimesMatrix(R22_1(1:3,:));

B = zeros(4, 4);
A = zeros(4, 4, pointCount);
for ii=1:pointCount
    A(1:4,1:4,ii) = [0, R12(ii,1:3); R21(1:3,ii), R22(1:3,1:3,ii)];
    B = B + A(:,:,ii)' * A(:,:,ii);
end

[~, S, V] = svd(B);
quat = V(:,4);
rot = quat2rot(quat);

T1 = [eye(3,3), -y_centroid ; 0 0 0 1];
T2 = [rot, [0; 0; 0]; 0 0 0 1];
T3 = [eye(3,3), x_centroid ;  0 0 0 1];

T = T3 * T2 * T1;
Eps = S(4,4);


text:
function y = log_sum_exp( x, dim )

%LOG_SUM_EXP   CVX internal version.

error( nargchk( 1, 2, nargin ) );
cvx_expert_check( 'log_sum_exp', x );

sx = size( x );
if nargin < 2 || isempty( dim ),
    dim = cvx_default_dimension( sx );
elseif ~cvx_check_dimension( dim, true ),
    error( 'Second argument must be a valid dimension.' );
end

%
% Quick exits
%

sx( end + 1 : dim ) = 1;
nx = sx( dim );
sy = sx;
sy( dim ) = 1;
if nx == 0,
    sx( dim ) = 1; %#ok
    y = -Inf * ones( sy );
    return
elseif nx == 1,
    y = x;
    return;
elseif any( sx == 0 ),
    y = zeros( sy );
    return
end

%
% Determine the expression types
%

persistent remap
if isempty( remap ),
    remap_2 = cvx_remap( 'real' );
    remap_1 = cvx_remap( 'convex' ) & ~remap_2;
    remap = remap_1 + 2 * remap_2;
end
v = reshape( remap( cvx_classify( x ) ), sx );
v = min( v, [], dim );

%
% Process each type of expression one piece at a time
%

vu = sort( v(:) );
vu = vu([true;diff(vu)~=0]);
nv = length( vu );
if nv > 1,
    y = cvx( sy, [] );
    if prod(sx(1:dim+1))>1 && prod(sx(dim+1:end))>1,
        perm = [ dim, 1:dim-1, dim+1:length(sx) ];
        x  = permute( x, perm );
        v  = permute( v, perm );
        y  = permute( y, perm );
        dim = 1;
    end
end
for k = 1 : nv,

    %
    % Select the category of expression to compute
    %

    vk = vu( k );
    if nv == 1,
        xt = x;
        sz = sy; %#ok
    else
        t = v == vk;
        xt = cvx_subsref( x, cvx_expand_dim( t, dim, nx ) );
        sx = size( xt );
        sz = sx;
        sz( dim ) = 1; %#ok
    end

    %
    % Perform the computations
    %

    switch vk,
        case 0,
            % Invalid
            error( 'Disciplined convex programming error:\n    Illegal operation: log_sum_exp( {%s} ).', cvx_class( xt ) );
        case 1,
            % Affine, convex
            cvx_begin
                variable w( sx )
                epigraph variable z( sz )
                { cvx_accept_convex( x ) - cvx_expand_dim( z, dim, nx ), 1, w } == exponential( sx );
                sum( w, dim ) == 1;
            cvx_end
        case 2,
            % Constant
            cvx_optval = cvx( log_sum_exp( cvx_constant( xt ) ) );
        otherwise,
            error( 'Shouldn''t be here.' );
    end

    %
    % Store the results
    %

    if nv == 1,
        y = cvx_optval;
    else
        y = cvx_subsasgn( y, t, cvx_optval );
    end

end

% Reshape again, just in case
y = reshape( y, sy );

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
%kkprdata 'Print Data in Formatted ASCII'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kprdata.pane file
%
% Parameters:
% InputFile: i 'Input ', required: 'Input Data Object'
% Toggle: attr 'Attributes', default: 0: 'Print Object Information & Attributes'
% Toggle: loc ' ', default: 0: 'Print Location Data'
% Toggle: time ' ', default: 0: 'Print Time Data'
% Toggle: val ' ', default: 0: 'Print Value Data'
% Toggle: mask ' ', default: 0: 'Print Mask Data'
% Toggle: ind ' ', default: 0: 'Print Implicit Indexing'
% String: replace 'Replacement for Masked Values', default: '*': 'Replacement String for Masked Values'
% OutputFile: f 'Output', optional: 'Resulting Output File'
% Toggle: mind ' ', default: 0: 'Print Implicit Map Indexing'
% Toggle: map ' ', default: 0: 'Print Map Data'
% Toggle: mval 'Map Value Data Through Map', default: 0: 'Map Value Data Through Map'
%
% Example: f = kkprdata(i, {'i','';'attr',0;'loc',0;'time',0;'val',0;'mask',0;'ind',0;'replace','*';'f','';'mind',0;'map',0;'mval',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kprdata - Print Data in Formatted ASCII
%
%  DESCRIPTION
% .I kprdata
% prints the selected segments of data from the input object in
% ASCII format to a file or to the screen.
%
% If the Print Attributes option (-att) is selected, all header and data
% information is printed before any data is printed.  (These non-data
% lines begin with a #, which indicates that the line is a "comment".)
%
% After the attribute lines, data fields selected by the user are formatted in
% the following order.
%
%
% Implicit Indexing: Location Data: Time Data: Value Data: Mask Data.
%
% This data can be followed by map data (if selected):
%
%
% Map Indexing: Map Data.
%
% This data can be followed by connectivity data.
%
% When map data exists, selecting the "Map Value Data Through Map" option
% will cause the value data to be mapped through the map data before being
% presented.
%
% The implicit indexing options, Print Implicit Indexing (-ind), and Print Map
% Indexing (-mind), if selected, will print the indices of the data.
% The index values index both the value data and explicit location and time
% data.  For instance, in the following example, the index [0][0][0][1]
% points to the 3D location data [4][27][38], and the value of 104.5.
%
%
% t  d  h  w    Location  Value
% 0  0  0  1    4 27 38   104.5
%
% The \fBReplacement for Masked Values" option (-replace) allows value data
% which has mask data associated with it, to be printed with the masked values
% replaced by the given string.
%
% Complex data is printed as real/imaginary pairs (real, imaginary).
%
% If data segments that do not exist in the input data object are selected,
% kprdata will return a warning message indicating that the segments do not
% exist.  Information about other selected segments will be printed.
%
% If the location grid is uniform or rectilinear and only the location
% data is being printed, then the location data will be printed in a
% manner consistent with the grid type.  That is, a uniform grid implies
% that the location begin and end will be printed, and a rectilinear
% grid implies that the data for the width, height, and depth axes will
% be printed.  When location data is printed with any other parts of the
% data model, it is presented as if it were curvilinear location data.
%
%
%
%  EXAMPLES
%
% kprdata -i in_object -attr -val -ind
% This command sequence will cause the attributes, followed by
% the value data, with indices, to be printed to stdout.
%
% kprdata -i in_obj -map -mind -val -f out_file
% This command sequence will print all of the value data,
% followed by the map data, with indices, to the file out_file.
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkprdata(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkprdata(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'attr', 0;'loc', 0;'time', 0;'val', 0;'mask', 0;'ind', 0;'replace', '*';'f', '__output';'mind', 0;'map', 0;'mval', 0};
maxval={0,0,0,0,0,0,0,0,1,0,0,0};
minval={0,0,0,0,0,0,0,0,1,0,0,0};
istoggle=[0,1,1,1,1,1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','String','OutputFile','Toggle','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=0; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kprdata"  '],Inputs,narglist);

repo: MATLAB


text:
%% Compare the two implementations
% Copyright 2013 The MathWorks, Inc.

%% Run the data through the two implementations
adjustedImageCPU = whitebalance(imageData);
adjustedImageGPU = whitebalance_gpu(imageData);

%% Visualize the results and their difference
subplot(1,2,1); imshow(adjustedImageCPU); title('CPU Output');
subplot(1,2,2); imshow(adjustedImageGPU); title('GPU Output');

difference = double(adjustedImageCPU(:) - adjustedImageGPU(:));
msgbox(['The norm of the difference is: ' num2str(norm(difference))]);


text:
function varargout=lpcaa2ao(varargin)
%
% For calling details please see v_lpcaa2ao.m
%
% This dummy routine is included for backward compatibility only
% and will be removed in a future release of voicebox. Please use
% v_lpcaa2ao.m in future and/or update with v_voicebox_update.m
%
%      Copyright (C) Mike Brookes 2018
%      Version: $Id: lpcaa2ao.m 10863 2018-09-21 15:39:23Z dmb $
%
if nargout
    varargout=cell(1,nargout);
    [varargout{:}]=v_lpcaa2ao(varargin{:});
else
    v_lpcaa2ao(varargin{:});
end


text:
function IM = generateQIBInfluence(IM, structROIV, sampleRateV)
%"generateQIBInfluence"
%   Uses the QIB engine to populate the beamlet fields of an IM structure.
%   StructROIV is the list of structures, sampleRateV is a vector of sample
%   rates, one for each structROIV.
%   The doses are stored in sparse format.
%   The stored index is with respect to the structure mask registered to the
%   uniformized CT scan.
%
%   This code was broken out of a function by JOD and CZ.
%
%JRA  26 Aug 2004
%LM:  14 Sept 05, JOD, fixed bug in call to mtoxyz
%     17 Dec  05, JOD, fixed bug in un-downsampled ROI
%                      definition; added comments to clarify calculation, and a small speedup change.
%     05 Jul  06, JJW, added missing inverse square law; added sigma_100
%
% 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/>.

%-----------Generate influence matrices---------------------%
%Loop over structures (i.e. goal terms)

global planC
global stateS
indexS = planC{end};

%obtain associated scanNum for structures. It is assumed that all the
%structures are associated to same scan (which is checked in IMRTP.m)
scanNum = getStructureAssociatedScan(structROIV(1));

for i = 1 : length(structROIV)

    if sampleRateV(i) ~= 1

        %Ensures that interpolative downsampling won't miss edge points:
        maskSingle3D = getSurfaceExpand(structROIV(i),0.5,1);
        if rem(log2(sampleRateV(i)),1) ~= 0
            error('Sample factor must (currently) be a power of 2.')
        end
        maskSample3D = getDown3Mask(maskSingle3D, sampleRateV(i), 1);
        %Get a mask of where to sample points
        tmp = logical(maskSample3D) & maskSingle3D;
        clear maskSample3D, maskSingle3D;

    else

        tmp = getUniformStr(structROIV(i));

    end

    scanIndV = find(tmp(:));  %Indices with respect to the uniformized scan.
    [rowV, colV, sliceV] = find3d(tmp); clear tmp;     %returns locations of mask voxels

    %[xV,yV,zV] = mtoxyz(rowV,colV,sliceV,1,planC,'uniform'); change by
    %JOD.
    [xV,yV,zV] = mtoxyz(rowV,colV,sliceV,scanNum,planC,'uniform');

    doseV = zeros(size(xV));

    numPts = length(xV);

    pM = [xV(:),yV(:),zV(:)];   %Location of dose calc points, downsampled.

    PBCounter = 1;

    QIBDataS = loadPBData;

    for j = 1 : length(IM.beams)

        disp(['Compute doses to structure number ' num2str(structROIV(i)) ' for beam ' num2str(j) '.'])

        %Get beamlet data for this beam and this structure

        sourceM = repmat([IM.beams(j).x, IM.beams(j).y, IM.beams(j).z],numPts,1);

        pRelM = pM - sourceM;    %get relative vector direction from source to dose calc points.

        RTOGPBVectorsM = IM.beams(j).RTOGPBVectorsM;

        str = int2str(size(RTOGPBVectorsM,1));

        disp(['numPBs = ' str])

        PBM = [];
        for PBNum = 1 : size(RTOGPBVectorsM,1)

            if mod(PBNum, 25) == 0 | (PBNum == size(RTOGPBVectorsM,1))
                disp(['Computed ' int2str(PBNum) ' out of ' str]); pause(0.003);
                try
                    IMRTPGui('status', j, length(IM.beams), structROIV(i), PBNum, size(RTOGPBVectorsM,1));
                end
            end

            PBV = RTOGPBVectorsM(PBNum,:);

            distSamplePts = IM.beams(j).CTTraceS(PBNum).distSamplePts;
            cumDensity    = IM.beams(j).CTTraceS(PBNum).cumDensityRay;

            distV = pRelM * PBV';    %Each row is the dot product of (the un-normalized) source-to-calc point
                                     %direction and the PB unit vector.
                                     %Hence, this is the distance along the PB ray line for each dose calc point.

            qM = sourceM + distV * PBV;  %The last term gives the vector i, j, and k components of
                                         %distance to the depth of closest approach.
                                         %Each row of sourceM is the i,j,k position of the source.
                                         %So vector sum locates the positions of closest approach
                                         %to the dose calc points.

            rM  = pM - qM;    %Each row of pM is the x, y, z location of a dose calc point.
                              %rM is then the vector pointing from the dose calc point to the point
                              %at which the PB vector makes closest approach.

            %Trap points which are too far away to need dose calcs
            %Uses sepsq if possible, much faster.
            try
                rDistVSquared = sepsq([0 0 0]', rM');  %mex file speedup
                goV = rDistVSquared < (IM.params.cutoffDistance^2);
            catch
                rDistV = (rM(:,1).^2 + rM(:,2).^2 + rM(:,3).^2);    %if mex not there
                goV = rDistV < IM.params.cutoffDistance^2;         %Mod by JOD to speedup, Dec 05.
            end

            rGoM = rM(goV,:);
            %       rGoDistV = rDistV(goV);
            distGoV = distV(goV);

            %We ignore the slight tilt of PB cross-sections with respect to the longitudinal axis.

            [gantryVectorsM] = RTOGVectors2Gantry(rGoM, IM.beams(j).gantryAngle);
            Xb = gantryVectorsM(:,1);
            Yb = gantryVectorsM(:,2);

            tmpV=clip(distGoV,0.001,max(distSamplePts) - 0.001,'limits');

            %get cum density at that point
            radDepthV = interp1([0, distSamplePts],[0, cumDensity],tmpV);

            %Get the widths of the PBs at that distance
            PBWidth_YbV = IM.beams(j).beamletDelta_y * distGoV/IM.beams(j).isodistance; %This is the width 'out-of-plane'
            PBWidth_XbV = IM.beams(j).beamletDelta_x * distGoV/IM.beams(j).isodistance; %This is the width 'in-plane'

            %compute dose using QIB
            [A_zV, a_zV, B_zV, b_zV] = GetPBConsts(radDepthV, IM.beams(j).beamEnergy, QIBDataS, 'nearest');

            doseFlag = IM.params.DoseTerm;

            sigmaVal_100 = IM.beams(j).sigma_100;

            doseV = getQIBDose([Xb,Yb], radDepthV, PBWidth_XbV, PBWidth_YbV, ...
                QIBDataS, A_zV, a_zV, B_zV, b_zV, IM.beams(j).beamEnergy, doseFlag, sigmaVal_100, distGoV);

            % apply inverse square law:
            % (divide by depth dependent area of beamlet scaled to isocenter distance)
            doseV = IM.beams(j).beamletDelta_x * IM.beams(j).beamletDelta_y * doseV ./ (PBWidth_XbV .* PBWidth_YbV);

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    %This function selectively eliminates some scatter components, if invoked.
            doseV = applyIMRTCompression(IM.params, doseV);
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

            %-----------Construct the dose matrix---------------------%

            indNZV = [doseV == 0];
            whereV = goV;     %First trap on distance
            ind0V = find(whereV);
            whereV(ind0V(indNZV)) = 0;    %zero out indices where dose was not computed from within QIB routine.
            [indV] = find(whereV);        %indV returns locations where dose is non-zero.
            doseV(indNZV) = [];           %eliminate dose entries which were zero.

            beamlet = createIMBeamlet(doseV, scanIndV(indV), j, length(goV));

            IM.beamlets(structROIV(i),PBCounter) = beamlet;
            IM.beamlets(structROIV(i),PBCounter).structureName = planC{planC{end}.structures}(structROIV(i)).structureName;
            IM.beamlets(structROIV(i),PBCounter).sampleRate = sampleRateV(i);

            PBCounter = PBCounter + 1;
        end

    end

    clear xV yV zV sourceM pM pRelM qM rM rDistV rGoM rowV sliceV;

end

repo: CERR