text:
%MAIN_doublePendulum.m
%
% This script runs a simulation of a double pendulum

% Use: EoM_Double_Pendulum to write the equations of motion

clear; clc;

%Physical parameters:
P.m1 = 4.0;  % (kg) pendulum mass, link 1
P.m2 = 1.0;  % (kg) pendulum mass, link 2
P.g = 9.81; % (m/s^2) gravity
P.l1 = 1.0; % (m) link 1 length
P.l2 = 2.0; % (m) link 2 length
P.d1 = 0.6*P.l1; % Distance from parent joint to link CoM
P.d2 = 0.4*P.l2; % Distance from parent joint to link CoM
P.I1 = P.m1*P.l1^2/12; %Link 1, moment of inertia
P.I2 = P.m2*P.l2^2/12; %Link 2, moment of inertia

tSpan = [0,10]; %Simulation time interval

z0 = zeros(4,1);
z0(1) = (pi/180)*(-90);  % Link 1, initial abs angle
z0(2) = 0;  % Link 1, initial abs angle rate
z0(3) = (pi/180)*(-90);  % Link 1, initial abs angle
z0(4) = 0;  % Link 1, initial abs angle rate

u = [0;0];

userFunc = @(t,z)doublePendulumDynamics(t,z,u,P);

options = odeset(...
    'AbsTol',1e-6,...
    'RelTol',1e-6,...
    'Vectorized','on');

% Run the simulation!
sol = ode45(userFunc,tSpan,z0,options);

% Break apart solution for plotting
nPlot = 1000;
time = linspace(tSpan(1),tSpan(2),nPlot);
z = deval(sol,time); %Evaluate solution from ode45 at points in time
th1 = z(1,:);
w1 = z(2,:);
th2 = z(3,:);
w2 = z(4,:);
[energy, kinetic, potential] = doublePendulumEnergy(z,P);

% Plotting!

figure(2222);clf;

subplot(3,1,1); hold on;
plot(time,th1,'r-','LineWidth',2);
plot(time,th2,'b-','LineWidth',2);
xlabel('time (s)')
ylabel('angle (rad)');

subplot(3,1,2); hold on;
plot(time,w1,'r-','LineWidth',2);
plot(time,w2,'b-','LineWidth',2);
xlabel('time (s)')
ylabel('angle rate (rad/s)');

subplot(3,1,3);  hold on;
datum = min(potential);
plot(time,energy - datum,'k-','LineWidth',3);
plot(time,kinetic,'LineWidth',2,'color',[0.5,0.1,0.6]);
plot(time,potential - datum,'LineWidth',2,'color',[0.3,0.5,0.2]);
xlabel('time (s)')
ylabel('energy (J)');
legend('total','kinetic','potential');

% Animation
figure(3333); clf;
doublePendulumAnimate(sol,P);


text:
% questdlg2() - questdlg function clone with coloring and help for
%               eeglab().
%
% Usage: same as questdlg()
%
% Warning:
% Case of button text and result might be changed by the function
%
% Author: Arnaud Delorme, CNL / Salk Institute, La Jolla, 11 August 2002
%
% See also: inputdlg2(), errordlg2(), supergui(), inputgui()

% Copyright (C) Arnaud Delorme, CNL / Salk Institute, arno@salk.edu
%
% 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

function [result] = questdlg2(Prompt,Title,varargin);

result = '';
if nargin < 2
   help questdlg2;
   return;
end;
if isempty(varargin)
    varargin = { 'Yes' 'No' 'Cancel' 'Yes' };
end;
result = varargin{end};
if Prompt(end) == 10, Prompt(end) = []; end;
if Prompt(end) == 10, Prompt(end) = []; end;
if Prompt(end) == 10, Prompt(end) = []; end;
if Prompt(end) == 10, Prompt(end) = []; end;

fig = figure('visible', 'off');
set(gcf, 'name', Title);

listui = {};
geometry = {};
if ~isempty(find(Prompt == 10))
    indlines = find(Prompt == 10);
    if indlines(1) ~= 1, indlines = [ 0 indlines ]; end;
    if indlines(end) ~= length(Prompt), indlines = [ indlines length(Prompt)+1 ]; end;
    for index = 1:length(indlines)-1
        geometry{index} = [1];
        listui{index} = { 'Style', 'text', 'string' Prompt(indlines(index)+1:indlines(index+1)-1) };
    end;
else
    for index = 1:size(Prompt,1)
        geometry{index} = [1];
        listui{index} = { 'Style', 'text', 'string' Prompt(index,:) };
    end;
end;
listui{end+1} = {};

geometry = { geometry{:} 1 ones(1,length(varargin)-1) };
for index = 1:length(varargin)-1 % ignoring default val
	listui = {listui{:} { 'width' 80 'align' 'center' 'Style', 'pushbutton', 'string', varargin{index}, 'callback', ['set(gcbf, ''userdata'', ''' varargin{index} ''');'] }  };
	if strcmp(varargin{index}, varargin{end})
		listui{end}{end+1} = 'fontweight';
		listui{end}{end+1} = 'bold';
	end;
end;

%cr = length(find(Prompt == char(10)))+1;
%if cr == 1
%	cr = size(Prompt,1);
%end;
%cr = cr^(7/);
%if cr >= 8, cr = cr-1; end;
%if cr >= 4, cr = cr-1; end;
%[tmp tmp2 allobj] = supergui( 'fig', fig, 'geomhoriz', geometry, 'geomvert', [cr 1 1], 'uilist', listui, ...
[tmp tmp2 allobj] = supergui( 'fig', fig, 'geomhoriz', geometry, 'uilist', listui, ...
    'borders', [0.02 0.015 0.08 0.06], 'spacing', [0 0], 'horizontalalignment', 'left' );

waitfor( fig, 'userdata');
try,
	result = get(fig, 'userdata');
	close(fig);
    drawnow;
end;

author: PatternRecognition

repo: OpenBMI


text:
function [fdr,x] = fdrCurve(fit, MODE, sgn)

% [fdr,x] = fdrCurve(fit, [MODE], [sgn])
%
% Computes the false discovery rate curve.
%
% Input:
%   fit     Poisson regression structure from function empNull.
%   MODE    Type of FDR
%               'tail' for tail FDR (default)
%               'local' for local FDR
%   sgn     1 (default) or -1, for right or left tail in 'FDR' case.
%           Irrelevant in 'fdr' case.
%
% Output:
%   fdr     3 column matrix, containing the fdr curve estimate and
%           the upper and lower limits corresponding to the 95% pointwise
%           confidence bands on the empirical distribution of T.
%   x       values at which fdr is evaluated
%
% See also:
%   fdrEmpNull.m

% HISTORY:
%   2006.05.08 ASH (armins@hsph.harvard.edu) wrote it.
%

if ~exist('fit'),
    error('Not enough arguments')
end
if ~exist('MODE'),
    MODE = 'tail';
end
if ~exist('sgn'),
    sgn = 1;
end

[x, X, W, y, yhat] = deal(fit.x, fit.X, fit.W, fit.y, fit.yhat);
K = length(fit.x);

if strcmp(MODE,'tail'),
    S = triu(ones(K,K),1) + diag(ones(1,K))/2;
    if (sgn == -1), S = S'; end
    fdr(:,1) = (S*yhat)./(S*y);
    V = X * inv(X'*W*diag(y)*W*X) * X' * W;
    V = diag(1./(S*yhat)) * S * diag(yhat) * V - diag(1./(S*y));
else % local
    fdr(:,1) = yhat./y;
    V = X * inv(X'*W*diag(y)*W*X) * X' * W - diag(1./y);
end

% Confidence intervals
cov = V * diag(y) * V';
v = diag(cov);
alpha = 0.05;
fdr(:,2) = fdr(:,1) .* exp(norminv(1-alpha/2) * sqrt(v));
fdr(:,3) = fdr(:,1) .* exp(-norminv(1-alpha/2) * sqrt(v));

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function myResize()
% 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/>.

tables = findobj('tag', 'UITable');
for i=1:length(tables)
    table('redraw', tables(i));
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function savemap()
    % save a current map and its coordinates
    % turned into function by Celso G Reyes 2017

    ZG=ZmapGlobal.Data; % used by get_zmap_globals

    [XX,YY] = meshgrid(gx,gy);
    [nr, nc] = size(valueMap);

    dats = [  reshape(XX,nr*nc,1) reshape(YY,nr*nc,1) reshape(valueMap,nr*nc,1) ];

    [filename, pathname] = uiputfile( ...
        {'*.dat'}, ...
        'Save as');

    fid = fopen([pathname filename],'w') ;
    fprintf(fid,'%9.4f  %9.4f %12.5f \n',dats');
    fclose(fid);
end

author: CelsoReyes


text:
function P=ar1spectrum(ar1,period)
% AR1 power spectrum...
%
% power=ar1spectrum(ar1,period)
%
%
% (c) Aslak Grinsted 2002-2004
%

% -------------------------------------------------------------------------
%   Copyright (C) 2002-2004, Aslak Grinsted
%   This software may be used, copied, or redistributed as long as it is not
%   sold and this copyright notice is reproduced on each copy made.  This
%   routine is provided as is without any express or implied warranties
%   whatsoever.

freq=1./period;
P=(1-ar1.^2)./(abs(1-ar1.*exp(-2*pi*i*freq))).^2; %http://www.madsci.org/posts/archives/may97/864012045.Eg.r.html
%fixed typo in numerical recipes

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
% This program implements the Kinect IMU integration as discussed by Huai
% Formulated in a local frame with camera measurements
% the output trajectory is that of the IMU

% we face the problem that some IMU records are missing, but it is
% generally very small gap, so it is reasonable to use the received
% measurements for Kalman update once the imu epoch exceeds their epoch.

% Test cases (1) start IMU and Kinect in stationary mode, do ZUPT for IMU
% for at least 1 minutes, the initial frame of IMU and the camera are s0
% and c0 frame, then when Kinect is moved smoothly in predefined
% trajectory, the filter is used to estimate gravity in s0 frame, the
% calibration terms of camera and IMU

function imu_cam_calib()
addpath('..\instk'); % imu functions
addpath('..\voicebox\'); % for rotro2qr and rotqr2eu, they are more robuts
% than dcm2quat_v000 and dcm2euler_v000
import java.util.LinkedList
clear variables;
clc; close all; format longg;
fprintf('\n IMU_Kinect_Calib.m to test EKF filtering in calibrating mems IMU and Kinect!\n\n');
rng('default');

experim=4; %test case for IMU CAMera calibration

switch experim

    case 4
        % test on Steval MKI062V2 data integration with camera, Mar 2014
        % output options
        % test observations: block camera measurement in the ZUPT period,
        % or turn off ZUPT, does not change the results much of turn-on ZUPT and camera
        % measurement at the same time.
        % the strange profile of Cs2c has to do with the data,
        % Ts2c wrong becuase of the filter mechanism and data.
        % use of gravity magnitude constraint does not show much improvement
        % it is advised to turn off scale factor estimate, let gravity and
        % RS2C float
        isOutNED=true;
        resdir='C:\Users\huai.3\kinectfusion\KinectFusionExplorer-D2D\temp\';
        filresfile=[resdir, 'filresult.bin']; % navigation states
        imuresfile=[resdir, 'imuresult.bin']; % imu error terms
        % imu options
        options.startTime=2;
        options.endTime=130;
        options.imuErrorModel=4; % how bias and scale factor is modeled
        options.mechanization=2; % 1 for wander azimuth
        options.imutype=6;      %  Steval MKI062V2
        options.dt=1/50;  %sampling interval
        options.maxCovStep=options.dt; %maximum covariance propagation step, if equal to dt, means single speed mode
        options.Cb2imu=eye(3);
        options.Timu2body=zeros(3,1); % h764G is the body frame
        options.imufile=[resdir,'iNemo_02092014_061806.tsv'];

       % Toc2s=[[  0.0970;   -0.2459;    0.9432];-[   -0.0037;   -0.0537;   -0.0680]]; % PoseIMUKinectx302010
       % Toc2s=[[   0.0886;   -0.2499;    0.9400]; -[  -0.0006;   -0.0526;   -0.0681]]; % for PoseIMUKinecty102030
      %  Toc2s=[[0.1194;   -0.2479;    0.9359]; -[  -0.0022;   -0.0530;   -0.0676]];% for PoseIMUKinectz201306
      %  Toc2s=[[   0.177695680320824;-0.00680201661415064;0.962107115439917]; -[  0.0362812169298294;-0.0710746054253426;0.0478464512175324]]; % for PoseIMUCircle
       % Toc2s=[[   -0.219889874055416; -0.230975516372796;0.907607304785909]; -[ -0.0084376;   -0.057877;       -0.068204]]; % for PoseIMUKinect3
       Toc2s=[0.285567755532141, -0.0225589251844046, 1.0874587565859,  0, 3/180*pi, 4/180*pi]'; % for iNEMO2_20140603_160228 anticlock
        % ZUPT start and end Time, determine before the filter
        % format n x 2
        % n is the number of segment
        % first  column is the start time
        % second column is the end time, for example [575908      576243; 576502      576602]
        options.zuptSE=[2, 50]; % for iNEMO2_20140603_160228 anticlock
        imuFileType=3; % iNemo data
       % imuFileType=2; % file provided by Yujia, mix of steval and kinect output
        %Initial PVA of IMU
        options.inillh_ant=[40.00311241687*pi/180; -83.01529294250*pi/180; 227.901];
        % lat, log, height of the antenna, in this case, identical to IMU,
        % guessed from http://www.daftlogic.com/sandbox-google-maps-find-altitude.htm
        options.imuErrors=zeros(12,1);
        options.imuErrors(1:6)=Toc2s; % for the Steval case, comment this line if not steval
        options.Vn=[0;0;0];
        options.Ve=[0;0;0];
        options.qb2n=rotro2qr(R3(0.1477)*[0, -1,0; -1,0, 0; 0,0, -1]);
        % 0.1477 is roughly estimated from earth surface triangulation using data from daftlogic
        options.InvalidateIMUerrors=false;
        options.initAttVar=0.5*pi/180; % 1 deg std for roll and pitch, 2 times 5 deg for yaw std
        % gps options
        useGPS=false;

        % ZUPT options
        options.sigmaZUPTPos = 0.02;% unit m
        options.sigmaZUPTVel = 0.02;% unit m/s
        options.sigmaZUPTAng = 1*180/pi;% unit rad
        rateZUPT=round(sqrt(1/options.dt)); % from tests, we see zupt has adverse effect in this case
        % NHC options
        options.sigmaNHC = 0.1;% unit m/s
        rateNHC=inf; %round(sqrt(1/options.dt));
        % minimum velocity before applying the NHC, this option decouples ZUPT and NHC
        options.minNHCVel=2.0;

        rateGravNorm=inf; %rateZUPT;
        sigmaGravMag=3e-3; % unit m/s^2

        options.useCam=false;
        camFile=options.imufile; % where the camera measurements come from
        options.camPoseFile=[resdir, 'kinectPose.txt']; % output camera position and attitude
        deltaPhi=[0;0;0]; % [-0.487; 1.106; -0.518]/180*pi;
        options.Cimu2cam=rotqr2ro(rvec2quat_v000(deltaPhi))*[0,0,1;-1,0,0;0,-1,0]'; % (Rc2s)'
        options.Tcam2body=[3; -6.2; -4.05]*1e-2; % unit m
        % depth camera in body frame whichStep is assumed to be the imu frame,

        sigmaEuler=[1; 1; 1]*pi/180; % unit rad
        sigmaTrans=[2;2;2]*1e-3; % unit m
    otherwise
        error('Unsupported testing case!');
end
numPrevImuDataToKeep=40; % the number of data put in preimudata
% Initialize the model state and covariance of state, process noise and
% measurment noise
filter =EKF_filter_s0framelet(options);
preimudata=LinkedList();% record the previous imu data

% read in imu data
[fimu, imudata, preimudata]=readimuheader(options.imufile, preimudata, options.startTime, numPrevImuDataToKeep, imuFileType);
lastimu=preimudata.getLast();
preimutime=lastimu(1,end);
imuctr=1;   % to count how many IMU data after the latest GPS observations

%load the camera measurements with timestamp
if(options.useCam)
    camdata=load_KinectDataset(camFile, options.startTime);
    euc2c0imu=zeros(size(camdata,1),4);
    euc2c0imu(:,1)=camdata(:,1);
    imgepoch=camdata(1,1);
    whichStep=1;
else imgepoch=inf;
end
% Start the main INS
initime=preimutime;
imuaccum=zeros(6,1);    % record the accumulated imu measurements
curimutime=imudata(1,end);
covupt_time=preimutime; %the time that we last updated the covariance

fprintf('Estimating trajectory...\n');
ffilres=fopen(filresfile,'Wb'); % navigation states, 'W' use buffer and binary fwrite()
fimures=fopen(imuresfile,'Wb'); % imu errors

while (~feof(fimu)&&curimutime<options.endTime)
    %Write IMU's position and velocity, rpy and accel and gyro bias to the files
    filter.SaveToFile([], preimutime, ffilres);
    if (imudata(1)-initime)>60
        disp(['Process Time:' num2str(imudata(1))]);
        initime=imudata(1);
    end
    % time propagation of the IMU
    imuaccum=filter.ffun_state(imuaccum,[imudata(2:7);preimutime;curimutime]);

    %Update the covariance
    if (((curimutime-covupt_time)>=options.maxCovStep) ||(curimutime>=imgepoch))
        %propagate the covariance
        filter.ffun_covariance(imuaccum, covupt_time, curimutime);
        covupt_time=curimutime;
        imuaccum=zeros(6,1);
        %Record covariance and navigation solution for the smoother
        %Note that I am recording the predicted solutions. That is why the
        %backward part must use filtered solution.
        imuerrors=filter.imuErrors;
        adjcoeff=ones(size(imuerrors));
        % since the scale factor is in ppt rather than in ppm, we need to
        % convert it to ppm
        adjcoeff(filter.imuScaleFactorSIP-filter.imuBiasDriftSIP+(1:6))=1e-3;
        outimuerrors=imuerrors.*adjcoeff;
        fwrite(fimures,[curimutime;outimuerrors;...
            full(sqrt(diag(filter.p_k_k(filter.imuBiasDriftSIP+(0:11),...
            filter.imuBiasDriftSIP+(0:11))))).*adjcoeff],'double');
    end
    %Apply ZUPT. Don't apply zupt for each imu sample. Zupt
    %for each sample must be performed on nominal trajectory, not with
    %the Kalman filter.
    isStatic =~isempty(options.zuptSE) && ~isempty(find(((options.zuptSE(:,1)<=curimutime)&(options.zuptSE(:,2)>=curimutime))==1,1));
    isZUPT =mod(imuctr, rateZUPT)==0;
    if (isStatic&&isZUPT)
        measure=zeros(9,1);
        predict=[filter.rvqs0(1:6); unskew(eye(3)-rotqr2ro(filter.rvqs0(7:10)))];
        H=sparse([eye(9) zeros(9,size(filter.p_k_k,1)-9)]);
        R=zeros(9);
        R(1:3,1:3)=eye(3)*options.sigmaZUPTPos^2;
        R(4:6,4:6)=eye(3)*options.sigmaZUPTVel^2;
        R(7:9,7:9)=eye(3)*options.sigmaZUPTAng^2;
        filter.correctstates(predict,measure, H,R);
    end
    % apply gravity norm measurement
    isGravNorm =mod(imuctr, rateGravNorm)==0;
    if (isStatic&&isGravNorm)
        %Gravity (most time consuming part of ecef implementations)
        xyz_imu=filter.rqs02e(1:3)+quatrot_v000(filter.rqs02e(4:7),filter.rvqs0(1:3),0);
        Llh=ecef2geo_v000(xyz_imu,0);
        [Rn, Re, gn, sL, cL, WIE_E]=geoparam_v000(Llh);
        predict=sqrt(imudata(2:4)'*imudata(2:4));
        gs0=quatrot_v000(filter.rvqs0(7:10), imudata(2:4),1); % predicted gravity in s0 frame
        H=[zeros(1,9) 1/predict*gs0' zeros(1,filter.covDim-12)];
        measure=norm(gn,2);
        R=sigmaGravMag^2;
        filter.correctstates(predict,measure, H,R);
    end

    % image measurement
    if(curimutime>=imgepoch)
        % there are some gap in the imu recordings
        if(abs(curimutime-imgepoch)>options.dt)
            disp(['abs(curimutime-imgepoch)>options.dt at ' num2str(curimutime)...
                ' GTOW sec and image index ' num2str(whichStep) '!']);
        end
        % predict the image cooridnates of points in the states
        RTcnc0=camdata(whichStep,2:7)';
        % 6x1 vector. RTcnc0(1:3) euler angle of Rc2c0, RTcnc0(4:6), Tc0 2c
        [Rc2c0, Tc02c , H]=filter.computeH();
        predict=[unskew(eye(3)-Rc2c0*(roteu2ro('xyz',RTcnc0(1:3)))', 0); Tc02c-RTcnc0(4:6)];
        euc2c0imu(whichStep, :)=[imgepoch; rotro2eu('xyz', Rc2c0)]';
        % predicted error =predicted value-measured value,
        % measured error= measured value- measured value= 0, dummy
        measure=zeros(6,1); % dummy

        R=diag([sigmaEuler.^2;sigmaTrans.^2]);
        filter.correctstates(predict,measure, H,R);
        filter.SaveCamPoseandRqs02e(imgepoch,whichStep);

        % take the next image time
        if(whichStep<size(camdata,1))
            whichStep=whichStep+1;
            imgepoch=camdata(whichStep,1);
        else
            imgepoch=inf; % stop using images
        end
    end
    %     filter.SaveCamPoseandRqs02e(curimutime,whichStep);
    %Read the next imu data
    if(preimudata.size()==numPrevImuDataToKeep)
        preimudata.removeFirst();
    end
    preimudata.addLast(imudata(:,end));    % record previous imudata
    preimutime=curimutime;
    [fimu, imudata]=grabnextimudata(fimu, preimutime, imuFileType);
    if (isempty(imudata))
        break;
    else
        curimutime=imudata(1,end);
        imuctr=imuctr+1;
    end
end
fprintf(' done.\n\n');
% draw Rc2c0 predicted and measured values
% nextFig=1;
% f(nextFig) = figure;
% plot(camdata(:,1)-camdata(1,1),camdata(:, 2:4)*180/pi,'marker','+');
% hold on
% plot(euc2c0imu(:,1)-euc2c0imu(1,1),euc2c0imu(:, 2:4)*180/pi,'marker','v');
% grid
% xlabel('Time [s]')
% ylabel('degree')
% legend('Roll cam','Pitch cam','Yaw/Heading cam','Roll imu','Pitch imu','Yaw/Heading imu');
% title('Cc2c0 by Kinect and predicted by IMU');

fclose all;
% display the navigation results
plotnav;
end

function RoTrKinect=load_KinectDataset(dataFile, startTime)
% output time, euler angle of Rc2c0, T c0 in c
fid = load(dataFile); %'PoseIMUKinect3.tsv');
DataIMU = fid((1:2:size(fid,1)-1),:);
EulerKin = fid((2:2:size(fid,1)),2:4);
TranKin = fid((2:2:size(fid,1)),5:7);
obsindex = zeros(size(EulerKin,1),1);
for i = 2:size(EulerKin,1)
    if EulerKin(i,1) ~= EulerKin(i-1,1)
        obsindex(i) = 1;
    end
end
whom=logical(obsindex)& (DataIMU(:,1)>=startTime);
RoTrKinect=[DataIMU(whom,1), EulerKin(whom, :), TranKin(whom, :)];
end

repo: ekfmonoslam


text:
function fg = dtiXformFiberCoords(fg, xform, coordspace)
%Applies the coordinate xform to all the fiber coordinates
%
%  fg = dtiXformFiberCoords(fg, xform)
%
% Typically fibers are stored in ACPC space.  Sometimes we want the fibers
% in image space or MNI space.  This routine efficiently transforms the
% coordinates from the space they are in the new space defined by the
% homogeneous transform (4x4) in xform.
%
% The code is optimized to transform without an expensive loop, unless it
% seems that memory is tight. In that case, it just does a (slow) loop.
%
% Example:
%  xform = dtiGet(dtiH,'acpc2img xform');
%  fgImg = dtiXformFiberCoords(fg,xform,'img');
%
%  xform = dtiGet(dtiH,'img2acpc xform');
%  fgAcpc = dtiXformFiberCoords(fg,xform,'acpc');
%
% HISTORY:
% 2005.01.25 RFD: wrote it.
% 2006.05.17 RFD: added fallback to slow, memory-efficient loop when it
% seems that space is tight.
% 2006.11.21 RFD: now allow just the fibers cell array as input, so that
% this routine is useful even when you don't have a proper 'fg' struct.
% 2009.01.11 RFD: modified to use a hybrid routine that is a good compromise
% between speed and memory usage.
%
% (c) Stanford VISTA Team

% We take either a fibergroup struct or just the 'fibers' field (a cell
% array of fibers).
if(iscell(fg))
    wasCell = true;
    tmp = fg;
    clear fg;
    fg.fibers = tmp;
else
    wasCell = false;
end

% Operate on 1000 fiber coordinates at a time?
stepSize = 1000;
n = numel(fg.fibers);
fiberLen = cellfun('size',fg.fibers,2);
for jj = 1:stepSize:n
    sInd = jj;
    eInd = min(jj+stepSize-1,n);

    fc = mrAnatXformCoords(xform, horzcat(fg.fibers{sInd:eInd}),0);

    if(~isempty(fc))
        fiberCoord = 1;
        for ii=sInd:eInd
            fg.fibers{ii} = fc(:,(fiberCoord:fiberCoord+fiberLen(ii)-1));
            fiberCoord = fiberCoord + fiberLen(ii);
        end
    end
end

% This was AR's idea, and it was a good one.  But it creates a problem
% because the fg structures have different slots so we can't store them in
% dtiH arrays.  We get an error about unmatched structures.  If the
% dtiH.fiberGroups was a cell array, we would be OK.  But it is an array of
% structs.  Sigh.
if ~notDefined('coordspace')
    fg = fgSet(fg, 'coordspace', coordspace);
else
    fg = fgSet(fg, 'coordspace', []);
end

% Put it back as you found it
if(wasCell)
    fg = fg.fibers;
end

return;


text:
function toks = split_string(line, delim)
ntok = 0;
tmp1 = line;
while length(tmp1)>0
    [tmp2, tmp1] = strtok(tmp1, delim);
    ntok = ntok+1;
end
toks = cell(ntok, 1);
for i = 1:ntok
    [toks{i}, line] = strtok(line, delim);
end

author: zouchuhang

repo: LayoutNet


text:
%ZSAD Sum of absolute differences
%
% M = ZSAD(I1, I2) is the zero-mean sum of absolute differences between the
% two equally sized image patches I1 and I2.  The result M is a scalar that
% indicates image similarity, a value of 0 indicates identical pixel patterns
% and is increasingly positive as image dissimilarity increases.
%
% Notes::
% - The ZSAD similarity measure is invariant to changes in image brightness
%   offset.
%
% See also SAD, SSD, NCC, ISIMILARITY.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function m = zssd(w1, w2)

	w1 = w1 - mean(w1(:));
	w2 = w2 - mean(w2(:));

	m = abs(w1-w2);
    m = sum(m(:));

author: petercorke

repo: machinevision-toolbox-matlab


text:
function output = ...
  SynthesizeByEventBasedMethod(event_index, event_time, ...
                               harmonics_level_at_event, ...
                               harmonics_deviation_at_event, fs, f0i, ...
                               unvoiced_mask)
% Synthesize speech from event based information made from lean
% representation
% output = ...
%  SynthesizeByEventBasedMethod(event_index, event_time, ...
%                               harmonics_level_at_event, ...
%                               harmonics_deviation_at_event, fs, f0i, ...
%                               unvoiced_mask)
%
% output = ...
%  SynthesizeByEventBasedMethod(event_index, event_time, ...
%                               harmonics_level_at_event, ...
%                               harmonics_deviation_at_event, fs, f0i, ...
%                               unvoiced_mask, interp_method)
%
% Input argument
%
% event_index: event location in sample index
% event_time: actual event location (s)
% harmonics_level_at_event: power information (dB)
% harmonics_deviation_at_event: rondom relative power information (dB)
% fs: sampling frequency (Hz)
% f0i: fundamental frequency in sample time (Hz)
% unvoiced_mask: silent indicator (for each event)
%
% Return value
%
% output: synthesized signal

% Copyright 2016 Google Inc. All Rights Reserved
% Author: hidekik@google.com (Hideki Kawahara)

%% convert to FFT-based array
narginchk(7, 7);
[response_for_periodic_trim, response_for_random_trim, ola_fftl] = ...
  GenerateResponses(f0i, fs, harmonics_level_at_event, ...
                    harmonics_deviation_at_event, event_time, event_index);
%%
out_buffer = f0i*0;
noise_out_buffer = out_buffer;
synth_out_buffer = out_buffer;
synth_noise_buffer = out_buffer;
nData = length(out_buffer);
half_width_of_BLIT = 2 * fs / (0.55 * fs);
blit_base = -ceil(half_width_of_BLIT):ceil(half_width_of_BLIT);
initial_index = 1;
noise_spread = max(1,min(1,unvoiced_mask(:))); % 0.2 to 1
base_noise = randn(length(out_buffer),1);
periodic_working_buffer = zeros(ola_fftl, 1);
zero_vector = periodic_working_buffer;
ola_zero = round(fs / min(f0i)) + 1;
ola_base = (1:ola_fftl)' - ola_zero;
sample_index = (1:nData)';
for ii = 1:length(event_time)
  base_index = event_index(ii) - blit_base; % minus (-) is the proper sign
  [synth_out_buffer, out_buffer] = ...
    UpdatePeriodicComponent(response_for_periodic_trim, synth_out_buffer, ...
                            base_index, ii, zero_vector, event_time, ...
                            event_index, fs, ola_zero, ola_base, ola_fftl, ...
                            sample_index, unvoiced_mask, out_buffer);
  [synth_noise_buffer, noise_out_buffer, initial_index] = ...
    UpdateRandomComponent(response_for_random_trim, synth_noise_buffer, ...
                          initial_index, ii, noise_spread, event_time, ...
                          event_index, ola_zero, ola_base, ola_fftl, ...
                          sample_index, base_noise, noise_out_buffer);
end;
raw_blit = GenerateBlitImpulse(blit_base, half_width_of_BLIT);
blit_equalizer_response = DesignBlitEqualizer(raw_blit);
half_width_of_EQBLIT = round((length(blit_equalizer_response) - 1) / 2);
tmp_period = fftfilt(blit_equalizer_response, ...
                     [synth_out_buffer; zeros(2 * half_width_of_EQBLIT, 1)]);
synth_noise_out = synth_noise_buffer ./ sqrt(f0i);
synth_periodic_out = ...
  tmp_period(half_width_of_EQBLIT + (1:nData)) ./ sqrt(f0i);
output = synth_periodic_out + synth_noise_out;
end

function [response_for_periodic_trim, response_for_random_trim, ola_fftl] ...
  = GenerateResponses(f0i, fs, harmonics_level_at_event, ...
                      harmonics_deviation_at_event, event_time, event_index)
kResidualToSNR = 11; % calibration constant from residual to SNR (dB)
spectrum_out.used_f0 = f0i(event_index);
spectrum_out.sampling_frequency = fs;
spectrum_out.temporal_positions = event_time;
spectrum_out.harmonic_power_dB = harmonics_level_at_event;
spectral_envelope = CalculateSpectrumEnvelope(spectrum_out);
noise_out = spectrum_out;
noise_out.refined_f0 = f0i(event_index);
noise_out.frame_time = event_time;
noise_out.aperiodicity_matrix = harmonics_deviation_at_event;
aperiodicity_dB = ...
  CalculateAperiodicitySgram(noise_out) + kResidualToSNR;
aperiodicity_ratio = ConvertDecibelToPower(min(0, aperiodicity_dB));
maximum_response_length = 0.03; % default 30 ms
response_lengh = round(maximum_response_length * fs); % samples
maximum_noise_length = 2 / min(f0i);
noise_length = round(maximum_noise_length * fs); % samples
ola_fftl = 2 ^ ceil(log2(noise_length + response_lengh + 1));
fftl = (size(spectral_envelope ,1) - 1) * 2;

data_frequency_axis = (0:fftl - 1) / fftl * fs;
data_frequency_axis(data_frequency_axis > fs / 2) = ...
  data_frequency_axis(data_frequency_axis > fs / 2) - fs;
% TODO(hidekik) need check: originally 50 maybe 150 better? (hz)
% 1000 is the corner frequenccy (Hz)
low_noise_masker = ...
  GetSigmoidNoiseShaper(data_frequency_axis, 1000, 150);
spectral_envelope_DFTBIN = ...
  [spectral_envelope;spectral_envelope(end-1:-1:2,:)];
random_envelope_DFTBIN = ...
  [aperiodicity_ratio;aperiodicity_ratio(end-1:-1:2,:)];
random_envelope_SHAPED = ...
  diag(low_noise_masker) * random_envelope_DFTBIN;
periodic_part = max(0.0001, (1 - random_envelope_SHAPED));
periodic_spectrum = sqrt(periodic_part .* spectral_envelope_DFTBIN);
random_spectrum = ...
  sqrt(random_envelope_SHAPED .* spectral_envelope_DFTBIN);
%%
time_axis = (0:fftl -1)' / fs;
time_shaper = time_axis * 0 + 1;
time_shaper(time_axis > maximum_response_length) = 0;
time_segment = ...
  (time_axis(time_axis > 0.8 * maximum_response_length & ...
  time_axis <= maximum_response_length) ...
  -0.8 * maximum_response_length) / (0.2 * maximum_response_length);
time_shaper(time_axis > 0.8 * maximum_response_length & ...
  time_axis <= maximum_response_length) ...
  = 0.5 + 0.5 * cos(pi * time_segment);
cepstrum_for_periodic = ifft(log(periodic_spectrum));
complex_cepstrum_for_periodic = cepstrum_for_periodic;
complex_cepstrum_for_periodic(fftl / 2 + 1:end,:) = 0;
complex_cepstrum_for_periodic(2:fftl / 2,:) = ...
  complex_cepstrum_for_periodic(2:fftl / 2,:) * 2;
response_for_periodic = diag(time_shaper) * ...
  real(ifft(exp(fft(complex_cepstrum_for_periodic))));
response_for_periodic_trim = response_for_periodic(1:response_lengh, :);
cepstrum_for_random = ifft(log(random_spectrum));
complex_cepstrum_for_random = cepstrum_for_random;
complex_cepstrum_for_random(fftl / 2 + 1:end,:) = 0;
complex_cepstrum_for_random(2:fftl / 2,:) = ...
  complex_cepstrum_for_random(2:fftl / 2,:) * 2;
response_for_random = diag(time_shaper) * ...
  real(ifft(exp(fft(complex_cepstrum_for_random))));
response_for_random_trim = response_for_random(1:response_lengh, :);
end

function [synth_out_buffer, out_buffer] = ...
  UpdatePeriodicComponent(response_for_periodic_trim, synth_out_buffer, ...
                          base_index, ii, zero_vector, event_time, ...
                          event_index, fs, ola_zero, ola_base, ola_fftl, ...
                          sample_index, unvoiced_mask, out_buffer)
%
nData = length(out_buffer);
half_width_of_BLIT = 2 * fs / (0.55 * fs);
blit_base = -ceil(half_width_of_BLIT):ceil(half_width_of_BLIT);
periodic_working_buffer = zero_vector;
fractional_index = event_time(ii) * fs - event_index(ii) + 1 + blit_base;
tmp_blit = GenerateBlitImpulse(fractional_index, half_width_of_BLIT);
periodic_working_buffer(ola_zero - blit_base) = tmp_blit / sum(tmp_blit);
periodic_working_out_buffer = ...
  real(ifft(fft(periodic_working_buffer) .* ...
            fft(response_for_periodic_trim(:, ii), ola_fftl)));
out_buffer(max(1,min(nData,base_index))) = tmp_blit / sum(tmp_blit);
past_index_condition = event_index(ii) + ola_base > 0 & ...
  event_index(ii) + ola_base < nData;
synth_out_buffer(sample_index(event_index(ii) + ...
  ola_base(past_index_condition))) = ...
  periodic_working_out_buffer(past_index_condition) * ...
  (1 - unvoiced_mask(ii)) + ...
  synth_out_buffer(sample_index(event_index(ii) + ...
  ola_base(past_index_condition)));
end

function [synth_noise_buffer, noise_out_buffer, initial_index] = ...
  UpdateRandomComponent(response_for_random_trim, synth_noise_buffer, ...
                        initial_index, ii, noise_spread, event_time, ...
                        event_index, ola_zero, ola_base, ola_fftl, ...
                        sample_index, base_noise, noise_out_buffer)
%----
nData = length(noise_out_buffer);
prev_index = initial_index:event_index(ii);
if ii == length(event_time)
  post_index = event_index(ii):nData;
else
  post_index = event_index(ii):event_index(ii + 1);
end;
prev_time = (prev_index(:) - event_index(ii)) / length(prev_index);
post_time = (post_index(:) - event_index(ii)) / length(post_index);
noise_shape = 0.5 * [(1 + cos(pi*prev_time / noise_spread(ii))) .* ...
  (abs(prev_time / noise_spread(ii)) <= 1); ...
  (1 + cos(pi * post_time(2:end) / noise_spread(ii))) .* ...
  (abs(post_time(2:end) / noise_spread(ii)) <= 1)];
noise_out_buffer(prev_index(1):post_index(end)) = ...
  noise_out_buffer(prev_index(1):post_index(end)) + noise_shape .* ...
  base_noise(prev_index(1):post_index(end)) / sqrt(sum(noise_shape .^ 2));
initial_index = event_index(ii);
noise_base_index = (prev_index(1):post_index(end))' - event_index(ii);
noise_working_buffer(max(1,ola_zero + noise_base_index)) = ...
  noise_shape .* base_noise(prev_index(1):post_index(end)) / ...
  sqrt(sum(noise_shape .^ 2));
noise_working_out_buffer = ...
  real(ifft(fft(noise_working_buffer(:), ola_fftl) .* ...
            fft(response_for_random_trim(:, ii), ola_fftl)));
past_index_condition = event_index(ii) + ola_base > 0 & ...
  event_index(ii) + ola_base < nData;
synth_noise_buffer(sample_index(event_index(ii) + ...
  ola_base(past_index_condition))) = ...
  noise_working_out_buffer(past_index_condition) + ...
  synth_noise_buffer(sample_index(event_index(ii) + ...
  ola_base(past_index_condition)));
end

author: google


text:
function w=pdeleft(u)
% PDELEFT W=PDELEFT(U) is the nonlinear residual of the left
%         preconditioned pde example with C=20.
%
global rhsf prhsf
%
% Compute the low-order nonlinear term.
%
v=20*u.*(dxmf(u)+dymf(u));
%
% Apply fish2d to the entire pde.
%
w=u+fish2d(v)-prhsf;

repo: Aero-matlab


text:
function [a1, b1, c1, a2, b2, c2] = aaa(A, B, C)
%AAA  gives both solutions to the angle-angle-angle problem, in radians.
%
%   AAA(A, B, C) will result in NaNs if the existence condition
%   |pi - |A|-|B|| <= |C| <= pi - ||A| - |B||
%   is not met.
%
%   See also AAAD.

% Rody P.S. Oldenhuis
% Delft University of Technology
% oldenhuis@gmail.com
%
% Created    : 23/Feb/2009
% Last edited: 30/Nov/2012

    % first solution
    a1 = acos2( (cos(A) + cos(B).*cos(C)) ./ (sin(B).*sin(C)),  A );
    b1 = acos2( (cos(B) + cos(A).*cos(C)) ./ (sin(A).*sin(C)),  B );
    c1 = acos2( (cos(C) + cos(A).*cos(B)) ./ (sin(A).*sin(B)),  C );

    % second solution
    a2 = 2*pi - a1;
    b2 = 2*pi - b1;
    c2 = 2*pi - c1;

    % check constraints
    indices = ( ...
        abs(pi - abs(A)-abs(B)) > abs(C) | ...
        abs(C) > pi - abs(abs(A)-abs(B)) );
    a1(indices) = NaN;    a2(indices) = NaN;
    b1(indices) = NaN;    b2(indices) = NaN;
    c1(indices) = NaN;    c2(indices) = NaN;

end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function varargout = wfilters(wname,o)
%WFILTERS Wavelet filters.
%   [LO_D,HI_D,LO_R,HI_R] = WFILTERS('wname') computes four
%   filters associated with the orthogonal or biorthogonal
%   wavelet named in the string 'wname'.
%   The four output filters are:
%       LO_D, the decomposition low-pass filter
%       HI_D, the decomposition high-pass filter
%       LO_R, the reconstruction low-pass filter
%       HI_R, the reconstruction high-pass filter
%   Available wavelet names 'wname' are:
%   Daubechies: 'db1' or 'haar', 'db2', ... ,'db45'
%   Coiflets  : 'coif1', ... ,  'coif5'
%   Symlets   : 'sym2' , ... ,  'sym8', ... ,'sym45'
%   Discrete Meyer wavelet: 'dmey'
%   Biorthogonal:
%       'bior1.1', 'bior1.3' , 'bior1.5'
%       'bior2.2', 'bior2.4' , 'bior2.6', 'bior2.8'
%       'bior3.1', 'bior3.3' , 'bior3.5', 'bior3.7'
%       'bior3.9', 'bior4.4' , 'bior5.5', 'bior6.8'.
%   Reverse Biorthogonal:
%       'rbio1.1', 'rbio1.3' , 'rbio1.5'
%       'rbio2.2', 'rbio2.4' , 'rbio2.6', 'rbio2.8'
%       'rbio3.1', 'rbio3.3' , 'rbio3.5', 'rbio3.7'
%       'rbio3.9', 'rbio4.4' , 'rbio5.5', 'rbio6.8'.
%
%   [F1,F2] = WFILTERS('wname','type') returns the following
%   filters:
%   LO_D and HI_D if 'type' = 'd' (Decomposition filters)
%   LO_R and HI_R if 'type' = 'r' (Reconstruction filters)
%   LO_D and LO_R if 'type' = 'l' (Low-pass filters)
%   HI_D and HI_R if 'type' = 'h' (High-pass filters)
%
%   See also BIORFILT, ORTHFILT, WAVEINFO.

%   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 12-Mar-96.
%   Last Revision: 05-Jul-1999.
%   Copyright 1995-2001 The MathWorks, Inc.
% $Revision: 1.10 $

% Check arguments.
if errargn(mfilename,nargin,[1 2],nargout,[0 1 2 4 8]), error('*'); end
if errargt(mfilename,wname,'str') , error('*'), end

wname         = deblankl(wname);
[wtype,fname] = wavemngr('fields',wname,'type','file');
mat_f         = findstr('.mat',fname);
if mat_f
   try
     load(fname,'-mat');
   catch
     msg = ['invalid wavelet file : ' fname];
     errargt(mfilename,msg,'msg');
     error('*');
   end
end

if wtype==1                % orth. wavelet
    if ~isempty(mat_f)
        F = eval(wname);
    else
        F = feval(fname,wname);
    end
    [Lo_D,Hi_D,Lo_R,Hi_R] = orthfilt(F);

elseif wtype==2            % biorth. wavelet
    if isempty(mat_f)
        [Rf,Df] = feval(fname,wname);
    else
        if exist('Rf')~=1 | exist('Df')~=1
            msg = ['invalid biorthogonal wavelet file : ' fname];
            errargt(mfilename,msg,'msg');
            error('*');
        end
    end
    [Lo_D,Hi_D1,Lo_R1,Hi_R,Lo_D2,Hi_D,Lo_R,Hi_R2] = biorfilt(Df,Rf,1);
    if (nargout>4) & (nargin<2)
        varargout(5:8) = {Lo_D2,Hi_D1,Lo_R1,Hi_R2};
    end

else
    msg = ['The wavelet ' wname ' is not valid!'];
    errargt(mfilename,msg,'msg');
    error('*');
    return;
end

if nargin==1
    varargout(1:4) = {Lo_D,Hi_D,Lo_R,Hi_R};
else
    o = lower(o(1));
    switch o
        case 'd' , varargout = {Lo_D,Hi_D};
        case 'r' , varargout = {Lo_R,Hi_R};
        case 'l' , varargout = {Lo_D,Lo_R};
        case 'h' , varargout = {Hi_D,Hi_R};
        otherwise
            errargt(mfilename,'invalid argument value','msg');
            error('*');
    end
end

author: xingchenzhang

repo: VIFB

save_path: github-repos/MATLAB/xingchenzhang-VIFB


text:
function [njm] = tapas_zeromat(jm)
%% Zeros all but the first entry in each column of jm.
%
% Input
%   jm      -- Matrix.
%
% Output
%   njm     -- Matrix.
%

% aponteeduardo@gmail.com
% copyright (C) 2015
%

njm = zeros(size(jm));

j = 0;
i = 0;

while j < size(jm, 2)
    j = j + 1;
    i = 0;
    while i < size(jm, 1)
        i = i + 1;
        if jm(i, j);
            njm(i, j) = 1;
            break
        end
    end
end
end

author: translationalneuromodeling

repo: tapas


text:
function varargout = spm_resels_vol(varargin)
% computes the number of resels in a volume - a compiled routine
% FORMAT R = spm_resels_vol(V,W)
% V      -  is a memory mapped image volume.
%           Finite and non-zero values are considered to be part of
%           the search volume.
% W      -  smoothness of the component fields {FWHM in voxels}.
% R      - Resel counts, where:
%          R(1) - Euler Characteristic of the volume (number of connected
%                 components - number of holes).
%          R(2) - Resel Diameter (average over all rotations of the
%                 distance between two parallel planes tangent to the
%                 volume in resel space).
%          R(3) - Resel Surface Area (half the surface area of the
%                 volume in resel space).
%          R(4) - Resel Volume (the volume in resel space).
%_______________________________________________________________________
%
% Reference : Worsley KJ et al 1996, Hum Brain Mapp. 4:58-73
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id: spm_resels_vol.m 1143 2008-02-07 19:33:33Z spm $

%-This is merely the help file for the compiled routine
error('spm_resels_vol.c not compiled - see Makefile')

author: spm


text:
function test_bug1527

% MEM 2gb
% WALLTIME 00:10:00
% DEPENDENCY ft_sourceplot

load(dccnpath('/home/common/matlab/fieldtrip/data/test/bug1527.mat'));

% this should be enough to reproduce the error
ft_sourceplot(cfg, source)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [  ] = vshowSpectralData( vdata, wavelength, opts )
%VSHOWSPECTRALDATA Showing a set of spectral data volume as suplots of a
%figure.
%   See also SPECTRUMRGB, TIGHT_SUBPLOT.

% default parameter setting
npcol    = 6; % number of columns for subplots
% magsize  = 1; % magnification of the image size
width    = 1080; % width of the plot
textpos  = [12 24]; % position of the text on the sub-images
fontsize = 15; % font size of the text (wavelength nm) on images
fontweight = 'normal'; % font weight of the text (wavelength) on images
fontname = 'FixedWidth'; % font name of the text (wavelength nm) on images
% parameters from the input opts
if isfield(opts,'npcol'),           npcol = opts.npcol;      end
if isfield(opts,'magsize'),       magsize = opts.magsize;    end
if isfield(opts,'width'),           width = opts.width;      end
if isfield(opts,'textpos'),       textpos = opts.textpos;    end
if isfield(opts,'fontsize'),     fontsize = opts.fontsize;   end
if isfield(opts,'fontweight'), fontweight = opts.fontweight; end
if isfield(opts,'fontname'),     fontname = opts.fontname;   end

% nw = length(wavelength); % number of wavelength bands
[nr,nc,nw] = size(vdata); % number of rows, columns, and wavelength bands
nprow = ceil(nw/npcol); % number of rows for subplots
% f = figure;
% f = figure('Position',[50 50 1100 900]); % depending on the size of images
% f = figure('Position',[50 50 npcol*nc*magsize nprow*nr*magsize]); % depending on the size of images
f = figure('Position',[50 50 width width/(npcol*nc)*nprow*nr]); % depending on the size of images
h = tight_subplot(nprow,npcol,[.002 .002],[.01 .01],[.01 .01]);
set(f,'Color','white');
% display each wavelength in each black with the same order as the
% wavelength
for iw = 1:nw
    lambda = wavelength(iw); % wavelength (nm)
    % colormap with the corresponding RGB wavelength
    wmap = (gray*kron(ones(3,1),spectrumRGB(lambda)));
    wmap = wmap/max(wmap(:));

    % subplot(nprow,npcol,iw);
    axes(h(iw));
    imshow(vdata(:,:,iw),'Colormap',wmap);
    text(textpos(1),textpos(2),[num2str(lambda,3) ' nm'],'Color','white',...
        'FontWeight',fontweight,'FontName',fontname,'FontSize',fontsize);
end
% display white image for the remaining blacks
for iw = nw+1:nprow*npcol
    axes(h(iw));
    imshow(ones(size(vdata,2)));
end

end

repo: DeSCI


text:
function M = perform_local_dct_transform(M,dir,w)

% perform_local_dct_transform - perform a JPEG-like transfrom
%
%   M = perform_local_dct_transform(M,dir,w, q);
%
%   dir=+1 : forward transform
%   dir=-1 : backward transform
%
%   w is the size of the blocks.
%
%   Copyright (c) 2004 Gabriel Peyre

if nargin<3
    w = 8;
end

n = size(M,1);

for i=1:floor(n/w);
    for j=1:floor(n/w);
        selx = (i-1)*w+1:min(i*w,n);
        sely = (j-1)*w+1:min(j*w,n);
        if dir==1
            M(selx,sely) = perform_dct_transform(M(selx,sely),+1);
        else
            M(selx,sely) = perform_dct_transform(M(selx,sely),-1);
        end
    end
end

author: gpeyre

repo: matlab-toolboxes


text:
function deg = rad2deg(rad)

% RAD2DEG  Radians to degrees conversion.

%   Copyright 2008-2009 Joan Sola @ LAAS-CNRS.

deg = rad/pi*180;

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

%   # START GPL LICENSE

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

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

%   # END GPL LICENSE

author: joansola

repo: slamtb

save_path: github-repos/MATLAB/joansola-slamtb


text:
classdef CC < handle & matlab.mixin.Copyable

    properties (Access = public)
        value
        gradient
    end

    properties (GetAccess = public, SetAccess = private)
        shapeFunctions
        nSF
    end

    properties (Access = private)
        sizeDesigVar
        valueOld
        gradientOld
    end

    methods (Access = protected, Abstract)
        updateFields(obj)
    end

    methods (Access = public)

        function computeFunctionAndGradient(obj)
            obj.initValueAndGradient();
            for iSF = 1:length(obj.shapeFunctions)
                obj.shapeFunctions{iSF}.updateTargetParameters();
                obj.shapeFunctions{iSF}.computeFunctionAndGradient();
                obj.updateFields(iSF);
            end
        end

        function computeFunction(obj)
            obj.initValueAndGradient();
            for iSF = 1:length(obj.shapeFunctions)
                obj.shapeFunctions{iSF}.updateTargetParameters();
                obj.shapeFunctions{iSF}.computeFunction();
                obj.updateFields(iSF);
            end
        end

        function objClone = clone(obj)
            objClone = copy(obj);
        end

        function restart(obj)
            obj.value    = obj.valueOld;
            obj.gradient = obj.gradientOld;
        end

        function updateOld(obj)
           obj.valueOld    = obj.value;
           obj.gradientOld = obj.gradient;
        end

    end

    methods (Access = protected)

        function obj = init(obj,cParams)
            obj.nSF   = 0;
            obj.sizeDesigVar = size(cParams.designVar.value);
            obj.createShapeFunctions(cParams);
        end

    end

    methods (Access = private)

        function initValueAndGradient(obj)
            obj.value = 0;
            obj.gradient = zeros(obj.sizeDesigVar);
        end

        function createShapeFunctions(obj,cParams)
            nS = cParams.nShapeFuncs;
            for iS = 1:nS
                s = cParams.shapeFuncSettings{iS};
                s.designVariable = cParams.designVar;
                s.shNumber = iS;

                if isprop(cParams,'homogenizedVarComputer')
                    s.homogVarComputer = cParams.homogenizedVarComputer;
                end
                if isprop(cParams,'targetParameters')
                    s.targetParameters = cParams.targetParameters;
                end
                shapeFunction = ShapeFunctional.create(s);
                obj.append(shapeFunction);
            end
        end

        function append(obj,shapeFunction)
            obj.shapeFunctions{obj.nSF+1} = shapeFunction;
            obj.nSF = obj.nSF+1;
        end

    end

end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
clear all, close all, clc
A = imread('../../CH01_SVD/DATA/dog.jpg');
B = rgb2gray(A);

%%
n = 2; w = 'db1'; [c,s] = wavedec2(B,n,w);
plotwavelet2(c,s,2,'db1',128,'square');
colormap gray
axis off
axis tight
set(gcf,'Position',[100 100 600 800])
set(gcf,'PaperPositionMode','auto')
print('-dpng', '-loose', '../figures/2DWavelet_Square');

%%
n = 3; w = 'sym2'; [c,s] = wavedec2(B,n,w);
[H1,V1,D1] = detcoef2('all',c,s,1);
A1 = appcoef2(c,s,'haar',1);
V1img = wcodemat(V1,255,'mat',1);
H1img = wcodemat(H1,255,'mat',1);
D1img = wcodemat(D1,255,'mat',1);
A1img = wcodemat(A1,255,'mat',1);
imagesc(A1img)
figure
[H2,V2,D2] = detcoef2('all',c,s,2);
A2 = appcoef2(c,s,'haar',2);
V2img = wcodemat(V2,255,'mat',2);
H2img = wcodemat(H2,255,'mat',2);
D2img = wcodemat(D2,255,'mat',2);
A2img = wcodemat(A2,255,'mat',2);
imagesc(V2img)
figure
colormap(gray); rv = length(gray);
plotwavelet2(c,s,3,'haar',rv,'square')

figure
subplot(2,2,1)
imagesc(A1img)

subplot(2,2,2)
imagesc(V1img)

subplot(2,2,3)
imagesc(H1img)

subplot(2,2,4)
imagesc(D1img)

author: dynamicslab


text:
function [kW] = hp2kW(hp)
% Convert power from mechanical horsepower to kilowatts.
% Chad A. Greene 2012
kW = hp*0.745699872;

author: Sable

repo: mcbench-benchmarks


text:
function s =SNR(x,xrec)
signalPower= mean(x(:).^2);
noisePower = mean((x(:) - xrec(:)).^2);
s =  signalPower/noisePower;

author: thomaskuestner


text:
function m = SchmidTensor(sS)
% Schmid tensors for a list of slip systems
%
% Syntax
%   m = SchmidTensor(sS)
%
%   sigma = stressTensor.uniaxial(xvector)
%   crss = m:sigma
%
% Input
%  sS - list of @slipSystem
%
% Output
%  m - Schmid tensor @velocityGradientTensor
%

% actually we need to multiply with the slipping rates $\dot gamma(t)$ to
% obtain the velocity gradient tensors
m = SchmidTensor(sS.n,sS.b);

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function [h,t,w,w_times,halfh, auc] = fir2htw2(b,varargin)
% Estimates height, time to peak, and width of FIR response
%
% :Usage:
% ::
%
%     [h,t,w,w_times,halfh, auc] = fir2htw2(b,[hconstraint],[doplot],[colors cell])
%
% :Inputs:
%
%   **b:**
%        beta/estimate series for hemodynamic response curve
%
%   hconstraint:**
%        max time in samples that can be considered the peak
%        (default = last sample)
%
%   doplot:**
%        flag for plot, 1/0
%
%   colors:**
%        cell vector of colors for plot
%
% ..
%    tor wager, 2/14/05
% ..
%
% minh = min height
%
% This version uses turning points (zero gradient) to find the largest
% "hump" in the data and the time it occurs.
%
% :Example:
% ::
%
%    hrf = spm_hrf(.5); hrf = hrf ./ max(hrf); hrf = hrf + .1 * randn(length(hrf), 1);
%    create_figure('hrf'); plot(hrf);
%    [h,t,w,w_times,halfh, auc] = fir2htw2(hrf, [], 1);
%

    if size(b,2) < length(b), b = b'; end

    if ~isempty(varargin) && ~isempty(varargin{1})
        hconstraint = varargin{1};
    else
        hconstraint = length(b);
    end

    if hconstraint > length(b)
        warning('fir2htw is trying to use more betas than there are!  limiting.')
        hconstraint = length(b);
    end

    if length(varargin) > 1,
        doplot = varargin{2};
    else
        doplot = 1;
    end

    % exit if empty
    if isempty(b) || all(b==0) || all(isnan(b))
        h = NaN; t = h; w = h; w_times = [NaN NaN]; halfh = NaN;
        warning('fir2htw:  Empty, nan, or all zero values!')
        return
    end

    % exit if all vals are same
    if all(b == mean(b)), h = NaN; t = h; w = h; w_times = [NaN NaN]; halfh = NaN;
        warning('fir2htw:  all values are the same!')
        return
    end

    colors = {'ro-'};
    if length(varargin) > 2, colors = varargin{3}; end

    % find turning points in data
    tmp = diff(b);
    turnpts = find([0 diff(sign(tmp))]);

    if isempty(turnpts)
        % monotonic increasing or decreasing function, can't define
        turnpts = [find(b==min(b)) find(b==max(b))];
    end

    dat = b(1:hconstraint) - mean(b(1:2));   %- b(1); % ; - mean(b(1:2));          % deviations from first point (baseline)

    turnpts(find(turnpts > length(dat))) = [];  % eliminate extra turn points after hconst

    if isempty(turnpts)
        % monotonic increasing or decreasing function, can't define
        turnpts = [find(dat==min(dat)) find(dat==max(dat))];
    end

    tpdat = abs(dat(turnpts));

    t = turnpts(tpdat == max(tpdat));       % time of max absolute turning point
    t = t(1);
    h = dat(t);                             % the value at max

    % calculations for width
    halfh = .5 * h;                         % 1/2 the max to min distance

    %d = distance(halfh,dat);

    if h > 0
        wh = find(dat > halfh);
    else
        wh = find(dat < halfh);
    end

    if isempty(wh), wh = NaN; end

    % first half
    x2 = max(wh(1),1);    % first above halfh
    x1 = max(wh(1)-1,1);  % first above-half - 1
    y2 = dat(x2);
    y1 = dat(x1);
    m = y2 - y1;          % slope, x reduces to 1 (x2 - x1 = 1)
    if m == 0
        w_times(1) = x1;    % exact match
    else
        w_times(1) = x1 + (halfh - y1) ./ m; % solve y = mx + b for x* given m; y1 = b
    end

    % second half
    x1 = min(wh(end),hconstraint);      % last one above halfh
    x2 = min(wh(end)+1,hconstraint);    % + 1
    y2 = dat(x2);  y1 = dat(x1);
    m = y2 - y1;                      % slope, x reduces to 1 (x2 - x1 = 1)
    if m == 0
        w_times(2) = x1;    % exact match
    else
        w_times(2) = x1 + (halfh - y1) ./ m;     % solve y = mx + b for x* given m; y1 = b
    end

    w = w_times(2) - w_times(1);        % width in elements

    % Add area under curve measure
    if nargout > 5, auc = sum(dat); end

    if doplot
        minh = mean(b(1:2));
        try

            hold on
            h1 = arrow([t minh],[t h+minh],'Length',12);        % height arrow
            h2 = arrow([0 h+minh],[t h+minh],'Length',12);      % delay arrow
            h3 = arrow([w_times(1) halfh+minh],[w_times(2) halfh+minh],'Length',12);
            h4 = arrow([w_times(2) halfh+minh],[w_times(1) halfh+minh],'Length',12);

            h5 = text(t + .1 * w, halfh+minh + .4*(halfh+minh),'h','FontSize',16,'FontWeight','bold');
            h6 = text(t - .5 * w, h+minh - .2*(halfh+minh),'t','FontSize',16,'FontWeight','bold');
            h7 = text(t - .5 * w, halfh+minh - .2*(halfh+minh),'w','FontSize',16,'FontWeight','bold');

            set(h1,'Color',colors{1}(1))
            set(h2,'Color',colors{1}(1))
            set(h3,'Color',colors{1}(1))
            set(h4,'Color',colors{1}(1))
            set(h5,'Color',colors{1}(1))
            set(h6,'Color',colors{1}(1))
            set(h7,'Color',colors{1}(1))

        catch
            warning('Error drawing arrows!')
        end
    end

end

% OLD EXTRA STUFF

% % range = max(b(1:hconstraint)) - min(b(1:hconstraint));
% %
% % % find all values within some tolerance of the half-max height
% % % if not enough points found to determine width, then interpolate until we
% % % get values.
% % % it would be better to interpolate!
% % tolval = .02; wh = []; resampleval = 1;
% %
% % tol = tolval * range;   % tolerance of tolval % of range
% % wh = find(d <= tol);    % find all points within tolerance
% %
% % while length(find(wh>t)) < 1 | length(find(wh<t)) < 1
% %     resampleval = resampleval + 1;
% %
% %     b2 = resample(b(1:hconstraint),resampleval,1); % interpolate
% %
% %     d = distance(halfh,b2);
% %     wh = find(d <= tol);    % find all points within tolerance
% %     wh = wh ./ resampleval;  % convert back to original time units
% %
% %     if resampleval > 10, warning('Could not find width!'), break, end
% % end
% %
% % % now find nearest elements to peak that are at max height
% % d_from_t = wh - t;
% % w_times = [max(d_from_t(d_from_t < 0)) min(d_from_t(d_from_t > 0))];
% %
% %
% %
% % if length(w_times) > 1
% %     w = w_times(2) - w_times(1);        % width in elements
% %     w_times = t + w_times;              % time indices of elements before and after peak at 1/2 max height
% % else
% %     w_times = [NaN NaN];
% %     w = NaN;                            % width is undefined; 1st or last tp was max
% % end

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function n = p31_n ( )

%*****************************************************************************80
%
%% P31_N returns the number of variables for problem 31.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    10 January 2001
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer N.  If N is positive, N represents the
%    only legal value for N for this problem.  If N is
%    negative, the absolute value of N is the least legal
%    value of N, but other values are allowable.
%
  n = 4;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function out = FieldMap_applyvdm(job)
% Apply VDM and reslice images
% FORMAT FieldMap_applyvdm(job)
% job.data(sessnum).scans   - images for session/run sessnum
% job.data(sessnum).vdmfile - VDM file for session/run sessnum
% job.roptions.rinterp - interpolation method
% job.roptions.wrap    - perform warp around in specified dimensions
% job.roptions.mask    - perform masking
% job.roptions.which(1) - reslice images in time series only
% job.roptions.which(2) - reslice images in time series and mean
% job.roptions.prefix   - prefix for vdm applied files
% job.roptions.pedir    - phase encode direction (i.e. aplly vdm file along
% this dimension
%__________________________________________________________________________
%
% A VDM (voxel displacement map) created using the FieldMap toolbox
% can be used to resample and reslice realigned images to the original
% subdirectory with the same (prefixed) filename.
%
% Voxels in the images will be shifted according to the values in the VDM
% file along the direction specified by job.roptions.pedir (i.e. this is
% usually the phase encode direction) and resliced to the space of the
% first image in the time series.
%
% Inputs:
% A job structure containing fields for the input data and the processing
% options. The input data contains the series of images conforming to
% SPM data format (see 'Data Format'), the relative displacement of the images
% is stored in their header and a VDM which has (probably) been created
% using the FieldMap toolbox and matched to the first image in the time
% series (this can also be done via the FieldMap toolbox).
%
% Outputs:
% The resampled and resliced images resliced to the same subdirectory with a prefix.
%__________________________________________________________________________
% Copyright (C) 2011-2014 Wellcome Trust Centre for Neuroimaging

% Chloe Hutton
% $Id: FieldMap_applyvdm.m 6258 2014-11-07 18:15:40Z guillaume $

tiny = 5e-2;

% assemble roptions
%--------------------------------------------------------------------------
flags.interp    = job.roptions.rinterp;
flags.wrap      = job.roptions.wrap;
flags.mask      = job.roptions.mask;
flags.which     = job.roptions.which(1);
flags.mean      = job.roptions.which(2);
flags.prefix    = job.roptions.prefix;
flags.pedir     = job.roptions.pedir;
hold = [repmat(flags.interp,1,3) flags.wrap];

% Determine dimension along which to apply vdm
applydim=flags.pedir;
if applydim~=1 &&  applydim~=2 && applydim~=3
    applydim=2;
end

% Gather up data into ds structure which holds images and vdm file
%--------------------------------------------------------------------------
P = {};
for i = 1:numel(job.data)
    P{i} = strvcat(job.data(i).scans{:});
    ds(i).P=spm_vol(P{i});
    if ~isempty(job.data(i).vdmfile)
        sfP{i} = job.data(i).vdmfile{1};
        ds(i).sfP=spm_vol(sfP{i});
    else
        sfP{i} = [];
    end
    ds(i).hold = [1 1 1 0 1 0];
end

ntot = 0;
for i=1:length(ds)
   ntot = ntot + length(ds(i).P);
end

% Set up  x y z for resampling
%--------------------------------------------------------------------------
[x,y,z] = ndgrid(1:ds(1).P(1).dim(1),1:ds(1).P(1).dim(2),1:ds(1).P(1).dim(3));
xyz = [x(:) y(:) z(:) ones(prod(ds(1).P(1).dim(1:3)),1)]; clear x y z;

% Create mask if required (usually default and required to create mean)
%--------------------------------------------------------------------------
if flags.mask || flags.mean,
    spm_progress_bar('Init',ntot,'Computing available voxels',...
        'volumes completed');

    if flags.mean
        Count    = zeros(prod(ds(1).P(1).dim(1:3)),1);
        Integral = zeros(prod(ds(1).P(1).dim(1:3)),1);
    end

    % if flags.mask
    msk = zeros(prod(ds(1).P(1).dim(1:3)),1);
    % end

    % To create mean, read each session specific vdmfile in
    % to the space of the first image of first session
    tv = 1;
    for s=1:length(ds)
        T = ds(s).sfP.mat\ds(1).P(1).mat;
        txyz = xyz * T';
        c = spm_bsplinc(ds(s).sfP,ds(1).hold);
        ds(s).sfield = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),ds(1).hold);
        ds(s).sfield = ds(s).sfield(:);
        clear c txyz;
        sess_msk = zeros(prod(ds(s).P(1).dim(1:3)),1);
        % Read in each images in space of first image of first session
        for i = 1:numel(ds(s).P)
            T = inv(ds(s).P(i).mat) * ds(1).P(1).mat;
            txyz = xyz * T';
            txyz(:,applydim) = txyz(:,applydim) + ds(s).sfield;
            tmp       = false(size(txyz,1),1);
            if ~flags.wrap(1), tmp = tmp | txyz(:,1) < (1-tiny) | txyz(:,1) > (ds(s).P(i).dim(1)+tiny); end
            if ~flags.wrap(2), tmp = tmp | txyz(:,2) < (1-tiny) | txyz(:,2) > (ds(s).P(i).dim(2)+tiny); end
            if ~flags.wrap(3), tmp = tmp | txyz(:,3) < (1-tiny) | txyz(:,3) > (ds(s).P(i).dim(3)+tiny); end
            sess_msk = sess_msk + real(tmp);
            spm_progress_bar('Set',tv);
            tv = tv+1;
        end
        msk = msk + sess_msk;
        if flags.mean, Count = Count + repmat(length(ds(s).P),prod(ds(s).P(1).dim(1:3)),1) - sess_msk; end

        %
        % Include static field in estmation of mask.
        %
        if isfield(ds(s),'sfP') && ~isempty(ds(s).sfP)
            T = inv(ds(1).sfP.mat) * ds(1).P(1).mat;
            txyz = xyz * T';
            tmp  = false(size(txyz,1),1);
            if ~flags.wrap(1), tmp = tmp | txyz(:,1) < (1-tiny) | txyz(:,1) > (ds(1).sfP.dim(1)+tiny); end
            if ~flags.wrap(2), tmp = tmp | txyz(:,2) < (1-tiny) | txyz(:,2) > (ds(1).sfP.dim(2)+tiny); end
            if ~flags.wrap(3), tmp = tmp | txyz(:,3) < (1-tiny) | txyz(:,3) > (ds(1).sfP.dim(3)+tiny); end
            msk = msk + real(tmp);
        end

        if isfield(ds(1),'sfield') && ~isempty(ds(1).sfield)
            ds(1).sfield = [];
        end
    end
    if flags.mask, msk = find(msk ~= 0); end
end

% Apply fieldmap to all files, looping through sessions
%--------------------------------------------------------------------------
spm_progress_bar('Init',ntot,'Reslicing','volumes completed');
tv = 1;
for s = 1:numel(ds)

    % Get transformation between distortion field and first image
    T = ds(s).sfP.mat\ds(1).P(1).mat;
    txyz = xyz * T';
    c = spm_bsplinc(ds(s).sfP,ds(1).hold);
    ds(s).sfield = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),ds(1).hold);
    ds(s).sfield = ds(s).sfield(:);

    % Read in each images in space of first image of first session
    for i = 1:numel(ds(s).P)
        T = inv(ds(s).P(i).mat) * ds(1).P(1).mat;
        txyz = xyz * T';
        txyz(:,applydim) = txyz(:,applydim) + ds(s).sfield;
        c = spm_bsplinc(ds(s).P(i),hold);
        ima = spm_bsplins(c,txyz(:,1),txyz(:,2),txyz(:,3),hold);
        %
        % Write out resliced images
        %
        if flags.which
            PO         = ds(s).P(i);
            [pth,nm,xt,vr] = spm_fileparts(deblank(PO.fname));
            PO.fname       = fullfile(pth,[flags.prefix nm xt vr]);
            PO.mat     = ds(1).P(1).mat;
            PO.descrip = sprintf('spm - applied vdm');
            ivol       = ima;
            if flags.mask
                ivol(msk) = NaN;
            end
            ivol = reshape(ivol,PO.dim(1:3));
            PO   = spm_create_vol(PO);
            spm_write_vol(PO,ivol);
            if nargout > 0
                  out.sess(s).rfiles{i} = PO.fname;
            end
        end
        %
        % Build up mean image if so required.
        %
        if flags.mean
            Integral = Integral + nan2zero(ima);
        end
        spm_progress_bar('Set',tv);
        tv = tv+1;
    end
    if isfield(ds(s),'sfield') && ~isempty(ds(s).sfield)
        ds(s).sfield = [];
    end
end

% Write mean image
%--------------------------------------------------------------------------
if flags.mean
   % Write integral image (16 bit signed)
   %-----------------------------------------------------------------------
   sw = warning('off','MATLAB:divideByZero');
   Integral   = Integral./Count;
   warning(sw);
   PO         = ds(1).P(1);
   [pth,nm,xt,vr] = spm_fileparts(deblank(ds(1).P(1).fname));
   PO.fname       = fullfile(pth,['mean' flags.prefix nm xt vr]);
   PO.pinfo   = [max(max(max(Integral)))/32767 0 0]';
   PO.descrip = 'spm - mean applied vdm image';
   PO.dt      = [spm_type('int16') spm_platform('bigend')];
   ivol = reshape(Integral,PO.dim);
   spm_write_vol(PO,ivol);
end

if nargout > 0
   out.rmean{1} = PO.fname;
end

spm_figure('Clear','Interactive');

%==========================================================================
function vo = nan2zero(vi)
vo = vi;
vo(~isfinite(vo)) = 0;
return;


text:
function [y,fs,bits]=wavread(fn,n)
%WAVREAD  Legacy MATLAB function to read .WAV file [Y,FS,BITS]=(FILENAME,NMAX)
% wavread supports multichannel data, with up to 32 bits per sample, and supports reading 24- and 32-bit .wav files.
%
% Usage:
% y = wavread('filename')               loads a WAVE file specified by the string filename, returning
%                                       the sampled data in y. The .wav extension is appended if no
%                                       extension is given. Amplitude values are in the range [-1,+1].
% [y,Fs,bits] = wavread('filename')     returns the sample rate (Fs) in Hertz and the number of bits
%                                       per sample (bits) used to encode the data in the file.
% [...] = wavread('filename',N)         returns only the first N samples from each channel in the file.
% [...] = wavread('filename',[N1 N2])   returns only samples N1 through N2 from each channel in the file.
% siz = wavread('filename','size')      returns the size of the audio data contained in the file in place
%                                       of the actual audio data, returning the vector siz = [samples channels].

%	   Copyright (C) Mike Brookes 2018
%      Version: $Id: wavread.m 10865 2018-09-21 17:22:45Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin<2
    [y,fs,wm,fx]=v_readwav(fn);
elseif ischar(n)
    if strcmp(n,'size')
        [y,fs,wm,fx]=v_readwav(fn,'',0);
        y=fx(4:5); % number of samples and channels
    else
        error('%s is invalid option',n);
    end
elseif length(n)<2
    [y,fs,wm,fx]=v_readwav(fn,'',n);
else
    [y,fs,wm,fx]=v_readwav(fn,'',n(2)-n(1)+1,n(1)-1);
end
bits=fx(7); % bits precision

author: ImperialCollegeLondon

repo: sap-voicebox


text:
function indx = r8vec_sort_heap_mask_a ( n, a, mask_num, mask )

%*****************************************************************************80
%
%% R8VEC_SORT_HEAP_MASK_A: indexed heap ascending sort of a masked R8VEC.
%
%  Discussion:
%
%    An array A is given.  An array MASK of indices into A is given.
%    The routine produces a vector INDX, which is a permutation of the
%    entries of MASK, so that:
%
%      A(MASK(INDX(I)) <= A(MASK(INDX(J))
%
%    whenever
%
%      I <= J
%
%    In other words, only the elements of A that are indexed by MASK
%    are to be considered, and the only thing that happens is that
%    a rearrangment of the indices in MASK is returned that orders the
%    masked elements.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the array.
%
%    Input, real A(N), an array to be index-sorted.
%
%    Input, integer MASK_NUM, the number of mask elements.
%
%    Input, integer MASK(MASK_NUM), the mask array.  This is
%    simply a list of indices of A.  The entries of MASK should
%    be unique, and each one should be between 1 and N.
%
%    Output, integer INDX(MASK_NUM), the sort index.  There are MASK_NUM
%    elements of A selected by MASK.  If we want to list those elements
%    in order, then the I-th element is A(MASK(INDX(I))).
%
  if ( n < 1 )
    return
  end

  if ( mask_num < 1 )
    return
  end

  if ( mask_num == 1 )
    indx(1) = 1;
    return
  end

  indx = i4vec_indicator1 ( mask_num );

  l = floor ( mask_num / 2 ) + 1;
  ir = mask_num;

  while ( 1 )

    if ( 1 < l )

      l = l - 1;
      indxt = indx(l);
      aval = a(mask(indxt));

    else

      indxt = indx(ir);
      aval = a(mask(indxt));
      indx(ir) = indx(1);
      ir = ir - 1;

      if ( ir == 1 )
        indx(1) = indxt;
        break
      end

    end

    i = l;
    j = l + l;

    while ( j <= ir )

      if ( j < ir )
        if ( a(mask(indx(j))) < a(mask(indx(j+1))) )
          j = j + 1;
        end
      end

      if ( aval < a(mask(indx(j))) )
        indx(i) = indx(j);
        i = j;
        j = j + j;
      else
        j = ir + 1;
      end

    end

    indx(i) = indxt;

  end

  return
end


text:
%  anova2_cell() - compute F-values in cell array using ANOVA.
%
% Usage:
%    >> [FC FR FI dfc dfr dfi] = anova2_cell( data );
%
% Inputs:
%   data       = data consisting of PAIRED arrays to be compared. The last
%                dimension of the data array is used to compute ANOVA.
% Outputs:
%   FC   - F-value for columns.
%   FR   - F-value for rows.
%   FI   - F-value for interaction.
%   dfc  - degree of freedom for columns.
%   dfr  - degree of freedom for rows.
%   dfi  - degree of freedom for interaction.
%
% Note: the advantage over the ANOVA2 function of Matlab statistical
%       toolbox is that this function works on arrays (see examples). Note
%       also that you still need the statistical toolbox to assess
%       significance using the fcdf() function. The other advantage is that
%       this function will work with complex numbers.
%
% Example:
%   a = { rand(1,10) rand(1,10) rand(1,10); rand(1,10) rand(1,10) rand(1,10) }
%   [FC FR FI dfc dfr dfi] = anova2_cell(a)
%   signifC = 1-fcdf(FC, dfc(1), dfc(2))
%   signifR = 1-fcdf(FR, dfr(1), dfr(2))
%   signifI = 1-fcdf(FI, dfi(1), dfi(2))
%
%   % for comparison
%   anova2(  [ a{1,1}' a{1,2}' a{1,3}'; a{2,1}' a{2,2}' a{2,3}' ], 10)
%
%   b = { [ a{1,1}; a{1,1} ] [ a{1,2}; a{1,2} ] [ a{1,3}; a{1,3} ];
%         [ a{2,1}; a{2,1} ] [ a{2,2}; a{2,2} ] [ a{2,3}; a{2,3} ] }
%   [FC FR FI dfc dfr dfi] = anova2_cell(b)
%
%   c{1,1} = reshape(repmat(b{1,1}, [2 1]),2,2,10);
%   c{1,2} = reshape(repmat(b{1,2}, [2 1]),2,2,10);
%   c{1,3} = reshape(repmat(b{1,3}, [2 1]),2,2,10);
%   c{2,3} = reshape(repmat(b{2,3}, [2 1]),2,2,10);
%   c{2,2} = reshape(repmat(b{2,2}, [2 1]),2,2,10);
%   c{2,1} = reshape(repmat(b{2,1}, [2 1]),2,2,10)
%   [FC FR FI dfc dfr dfi] = anova2_cell(c)
%
% Author: Arnaud Delorme, SCCN/INC/UCSD, La Jolla, 2005
%
% Reference:
%   Schaum's outlines in statistics (3rd edition). 1999. Mc Graw-Hill.

% Copyright (C) Arnaud Delorme
%
% 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

function [FC, FR, FI, freeC, freeR, freeI] = anova2_cell(data)

    % compute all means and all std
    % -----------------------------
    a = size(data,1);
    b = size(data,2);
    nd = myndims( data{1} );
    c  = size(data{1}, nd);

    % dataabs if for complex data only
    % --------------------------------
    dataabs = data;
    if ~isreal(data{1})
        for i = 1:a
            for ii = 1:b
                dataabs{i,ii} = abs(data{i,ii});
            end;
        end;
    end;

    if nd == 1

        VE = 0;
        m  = zeros( size(data), 'single' );
        for i = 1:a
            for ii = 1:b
                m(i,ii) = mymean(data{i,ii});
                VE      = VE+sum( (dataabs{i,ii}-m(i,ii)).^2 );
            end;
        end;
        X  = mean(mean(m));
        Xj = mean(m,2);
        Xk = mean(m,1);
        VR = b*c*sum( (Xj-X).^2 );
        VC = a*c*sum( (Xk-X).^2 );

        Xj = repmat(Xj, [1 size(m,2) ]);
        Xk = repmat(Xk, [size(m,1)  1]);
        VI = c*sum( sum( ( m - Xj - Xk + X ).^2 ) );

    elseif nd == 2

        VE = zeros( size(data{1},1),1, 'single');
        m  = zeros( [ size(data{1},1) size(data) ], 'single' );
        for i = 1:a
            for ii = 1:b
                tmpm = mymean(data{i,ii}, 2);
                m(:,i,ii) = tmpm;
                VE        = VE+sum( (dataabs{i,ii}-repmat(tmpm, [1 size(data{i,ii},2)])).^2, 2);
            end;
        end;
        X  = mean(mean(m,3),2);
        Xj = mean(m,3);
        Xk = mean(m,2);
        VR = b*c*sum( (Xj-repmat(X, [1 size(Xj,2)])).^2, 2 );
        VC = a*c*sum( (Xk-repmat(X, [1 1 size(Xk,3)])).^2, 3 );

        Xj = repmat(Xj, [1 1 size(m,3) ]);
        Xk = repmat(Xk, [1 size(m,2)  1]);
        VI = c*sum( sum( ( m - Xj - Xk + repmat(X, [1 size(m,2) size(m,3)]) ).^2, 3), 2 );

    elseif nd == 3

        VE = zeros( size(data{1},1), size(data{1},2), 'single' );
        m  = zeros( [ size(data{1},1) size(data{1},2) size(data) ], 'single' );
        for i = 1:a
            for ii = 1:b
                tmpm = mymean(data{i,ii}, 3);
                m(:,:,i,ii) = tmpm;
                VE          = VE+sum( (dataabs{i,ii}-repmat(tmpm, [1 1 size(data{i,ii},3)])).^2, 3);
            end;
        end;
        X  = mean(mean(m,4),3);
        Xj = mean(m,4);
        Xk = mean(m,3);
        VR = b*c*sum( (Xj-repmat(X, [1 1 size(Xj,3)  ])).^2, 3 );
        VC = a*c*sum( (Xk-repmat(X, [1 1 1 size(Xk,4)])).^2, 4 );

        Xj = repmat(Xj, [1 1 1 size(m,4) ]);
        Xk = repmat(Xk, [1 1 size(m,3)  1]);
        VI = c*sum( sum( ( m - Xj - Xk + repmat(X, [1 1 size(m,3) size(m,4)]) ).^2, 4 ), 3 );

    else % nd == 4

        VE = zeros( size(data{1},1), size(data{1},2), size(data{1},3), 'single' );
        m  = zeros( [ size(data{1},1) size(data{1},2) size(data{1},3) size(data) ], 'single' );
        for i = 1:a
            for ii = 1:b
                tmpm = mymean(data{i,ii}, 4);
                m(:,:,:,i,ii) = tmpm;
                VE            = VE+sum( (dataabs{i,ii}-repmat(tmpm, [1 1 1 size(data{i,ii},4)])).^2, 4);
            end;
        end;
        X  = mean(mean(m,5),4);
        Xj = mean(m,5);
        Xk = mean(m,4);
        VR = b*c*sum( (Xj-repmat(X, [1 1 1 size(Xj,4)  ])).^2, 4 );
        VC = a*c*sum( (Xk-repmat(X, [1 1 1 1 size(Xk,5)])).^2, 5 );

        Xj = repmat(Xj, [1 1 1 1 size(m,5) ]);
        Xk = repmat(Xk, [1 1 1 size(m,4)  1]);
        VI = c*sum( sum( ( m - Xj - Xk + repmat(X, [1 1 1 size(m,4) size(m,5)]) ).^2, 5 ), 4 );

    end;

    SR2 = VR/(a-1);
    SC2 = VC/(b-1);
    SI2 = VI/(a-1)/(b-1);
    SE2 = VE/(a*b*(c-1));

    FR = SR2./SE2; % rows
    FC = SC2./SE2; % columns
    FI = SI2./SE2; % interaction

    freeR = [ a-1 a*b*(c-1) ];
    freeC = [ b-1 a*b*(c-1) ];
    freeI = [ (a-1)*(b-1) a*b*(c-1) ];

function val = myndims(a)
    if ndims(a) > 2
        val = ndims(a);
    else
        if size(a,1) == 1,
            val = 2;
        elseif size(a,2) == 1,
            val = 1;
        else
            val = 2;
        end;
    end;

function res = mymean( data, varargin) % deal with complex numbers
    res = mean( data, varargin{:});
    if ~isreal(data)
        res = abs( res );
    end;

author: PatternRecognition

repo: OpenBMI


text:
function stroud_test30 ( )

%*****************************************************************************80
%
%% TEST30 tests SPHERE_UNIT_**_3D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    06 April 2009
%
%  Author:
%
%    John Burkardt
%
  global FUNC_3D_INDEX;

  fprintf ( 1, '\n' );
  fprintf ( 1, 'TEST30\n' );
  fprintf ( 1, '  For integrals on the unit sphere in 3D:\n' );
  fprintf ( 1, '  SPHERE_UNIT_07_3D uses a formula of degree 7.\n' );
  fprintf ( 1, '  SPHERE_UNIT_11_3D uses a formula of degree 11.\n' );
  fprintf ( 1, '  SPHERE_UNIT_14_3D uses a formula of degree 14.\n' );
  fprintf ( 1, '  SPHERE_UNIT_15_3D uses a formula of degree 15.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Unit sphere area = %f\n', sphere_unit_area_nd ( 3 ) );
  fprintf ( 1, '\n' );
  fprintf ( 1, ...
    '    F(X)    S3S07        S3S11         S3S14         S3S15\n' );
  fprintf ( 1, '\n' );

  num = function_3d_num ( );

  for i = 1 : num

    FUNC_3D_INDEX = i;

    result1 = sphere_unit_07_3d ( 'function_3d' );
    result2 = sphere_unit_11_3d ( 'function_3d' );
    result3 = sphere_unit_14_3d ( 'function_3d' );
    result4 = sphere_unit_15_3d ( 'function_3d' );

    fname = function_3d_name ( i );

    fprintf ( 1, '  %s  %12f  %12f  %12f  %12f\n', ...
      fname, result1, result2, result3, result4 );

  end

  return
end


text:
function tdetx = tdet(x,K)
% tdetx = tdet(x,K)
%
% TDET  Computes twice determinant for Lorentz block
%
% **********  INTERNAL FUNCTION OF SEDUMI **********
%
% See also sedumi

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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.,  51 Franklin Street, Fifth Floor, Boston, MA
% 02110-1301, USA

if isempty(K.q)
    tdetx = zeros(0,1);
else
    ix = K.mainblks;
    tdetx = x(ix(1):ix(2)-1).^2 - ddot(x(ix(2):ix(3)-1),x,K.qblkstart);
end

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function h = mrUtilScattercloud(x,y,n,l,clm,cmap)
%SCATTERCLOUD display density of scatter data
%   SCATTERCLOUD(X,Y) creates a scatterplot of X and Y, displayed over a
%   surface representing the smoothed density of the points.  The density is
%   determined with a 2D histogram, using 25 equally spaced bins in both
%   directions.
%   SCATTERCLOUD(X,Y,N) uses N equally spaced bins.
%   SCATTERCLOUD(X,Y,N,L) uses L as a parameter to the smoothing algorithm.
%    Defaults to 1.  Larger values of L lead to a smoother density, but a
%    worse fit to the original data.
%   SCATTERCLOUD(X,Y,N,L,CLM) uses CLM as the color/linestyle/marker for
%    the scatter plot.  Defaults to 'k+'.
%   SCATTERCLOUD(X,Y,N,L,CLM,CMAP) uses CMAP as the figure's colormap.  The
%    default is 'flipud(gray(256))'.
%   H = SCATTERCLOUD(...) returns the handles for the surface and line
%    objects created.
%
%   Example:
%
%     scattercloud(1:100 + randn(1,100), sin(1:100) + randn(1,100),...
%                  50,.5,'rx',jet(256))

error(nargchk(2,6,nargin),'struct');

x = x(:);
y = y(:);

if length(x) ~= length(y)
    error('SCATTERCLOUDDataVectorSizesDoNotMatch','The number of elements in x and y do not match')
end

if nargin < 6
    cmap = flipud(gray(256));
end

if nargin < 5
    clm = 'k+';
end

if nargin < 4
    l = 1;
end

if nargin < 3
    n = 25;
end

% min/max of x and y
minX = min(x);
maxX = max(x);
minY = min(y);
maxY = max(y);

% edge locations
xEdges = linspace(minX,maxX,n);
yEdges = linspace(minY,maxY,n);

% shift edges
xDiff = xEdges(2) - xEdges(1);
yDiff = yEdges(2) - yEdges(1);
xEdges = [-Inf, xEdges(2:end) - xDiff/2, Inf];
yEdges = [-Inf, yEdges(2:end) - yDiff/2, Inf];

% number of edges
numX = numel(xEdges);
numY = numel(yEdges);

% hold counts
C = zeros(numY,numX);

% do counts
for i = 1:numY-1
    for j = 1:numX-1
        C(i,j) = length(find(x >= xEdges(j) & x < xEdges(j+1) &...
                             y >= yEdges(i) & y < yEdges(i+1)));
    end
end

% get rid of Infs from the edges
xEdges = [xEdges(2) - xDiff,xEdges(2:end-1), xEdges(end-1) + xDiff];
yEdges = [yEdges(2) - yDiff,yEdges(2:end-1), yEdges(end-1) + yDiff];

% smooth the density data, in both directions.
C = localSmooth(localSmooth(C,l)',l)';

% create the graphics
ax = newplot;
s = surf(xEdges,yEdges,zeros(numY,numX),C,...
         'EdgeColor','none',...
         'FaceColor','interp');
view(ax,2);
colormap(ax,cmap);
grid(ax,'off');
holdstate = get(ax,'NextPlot');
set(ax,'NextPlot','add');
p = plot(x,y,clm);
axis(ax,'tight');
set(ax,'NextPlot',holdstate)

% outputs
if nargout
    h = [s;p];
end

function B = localSmooth(A,L)
r = size(A,1);
I = eye(r);
D1 = diff(I);
D2 = diff(I,2);
B = (I + L ^ 2 * D2' * D2 + 2 * L * D1' * D1) \ A;

author: vistalab

repo: vistasoft


text:
%  pageify -- turn databases into one db with one observation and many pages
%
%  ::
%
%
%    db=pageify(pivot_date,varargin)
%
%  Args:
%
%     - **pivot_date** [char|serial date] : reference date, typically end of
%     history.
%
%     - **varargin** [struct|ts] : databases in time series format or in struct
%
%     format
%
%  Returns:
%     :
%
%     - **db** [ts] : time series with one observation and many pages
%
%  Note:
%
%     - This routine is useful for preparing data for conditional forecasting
%
%  Example:
%
%     See also:
%

author: jmaih


text:
function [P] = prolongation(CV,CF,V,varargin)
  % PROLONGATION Build a linear prolongation operator taking solutions on a
  % coarse mesh (CV,CF) and prolongating (upsampling) them onto vertices of a
  % fine mesh (V). Useful for multigrid and embedded mesh deformations.
  %
  % P = prolongation(CV,CF,V)
  %
  % Inputs:
  %   CV  #CV by dim list of coarse mesh vertex positions
  %   CF  #CF by dim+1 list of coarse mesh element indices into CV
  %   V  #V by dim list of fine mesh vertex positions
  %   Optional:
  %     'Extrapolation'  followed by name of method to use for points in V
  %       lying outside (CV,CF). One of the following:
  %         {'linear'}   finds closest element and uses barycentric coordinates
  %           (with some negative weights). This is linearly precise.
  %         'contstant'  finds closest point and uses its barycentric
  %           coordinates (all non-negtive since it lies on a facet). This is
  %           only constantly precise.
  % Outputs:
  %   P  #V by #CV prolongation matrix so that X = P * CX prolongates a
  %     solution CX on the coarse mesh to a solution X on the fine mesh.
  %

  [I] = in_element_aabb(CV,CF,V);
  I0 = find(I==0);
  extrapolation = 'linear';
  % Map of parameter names to variable names
  params_to_variables = containers.Map( ...
    {'Extrapolation'}, ...
    {'extrapolation'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  switch size(CF,2)
  case 4
    [BF,BJ,BK] = boundary_faces(CF);
    BV = V;
    % Extrapolate
    % Snap to boudnary and use snapped points interpolation
    [~,BI,BVI0] = point_mesh_squared_distance(V(I0,:),CV,BF);
    I(I0) = BJ(BI);
    switch extrapolation
    case 'linear'
      % Leave BV
    case 'constant'
      BV(I0,:) = BVI0;
    otherwise
      error('Unknown extrapolation (%s)',extrapolation);
    end
    [IV] = max(I,[],2);
    % recover barycentric coordinates (interpolation & extrapolation)
    B = barycentric_coordinates( ...
      BV,CV(CF(IV,1),:),CV(CF(IV,2),:),CV(CF(IV,3),:),CV(CF(IV,4),:));
    %% Snap to closest vertex
    %B(any(B<0,2),:) = bsxfun(@eq,max(B(any(B<0,2),:),[],2),B(any(B<0,2),:));
  case 3
    % Dirty trick to make V,CV 3D
    BV = V;
    if ~isempty(I0)
      pV = V;
      pCV = CV;
      pBV = BV;
      pV(:,end+1:3) = 0;
      pCV(:,end+1:3) = 0;
      [~,I(I0),pBVI0] = point_mesh_squared_distance(pV(I0,:),pCV,CF);
      switch extrapolation
      case 'linear'
        % Leave BV
      case 'constant'
        BV(I0,:) = pBVI0(:,1:2);
      otherwise
        error('Unknown extrapolation (%s)',extrapolation);
      end
    end
    [IV] = max(I,[],2);
    B = barycentric_coordinates( ...
      BV,CV(CF(IV,1),:),CV(CF(IV,2),:),CV(CF(IV,3),:));
  end

  P = sparse( ...
    repmat(1:size(V,1),size(CF,2),1)', ...
    CF(IV,:), ...
    B, ...
    size(V,1),size(CV,1));
end

author: alecjacobson

repo: gptoolbox


text:
function ui = mrViewMeshPanel(ui, dockFlag, vis);
%
% ui = mrViewMeshPanel(ui, [dockFlag], [vis=0]);
%
% Adds a panel to the mrViewer UI with
% controls for displaying data onto 3D meshes.
%
% dockFlag: if 1,  will attach panel to the
% mrViewer figure. Otherwise,  will make a
% separate figure for the panel.
%
% vis: if 0, will initialize the panel to be hidden (default when
% starting up mrViewer); otherwise, leave it visible (as when
% docking/undocking the panel).
%
% ras,  07/13/06.
if ~exist('ui','var') | isempty(ui),   ui = mrViewGet; end
if ~exist('dockFlag','var') | isempty(dockFlag),   dockFlag = 1; end
if ~exist('vis','var') | isempty(vis),   vis = 0; end

mrGlobals2;

bgColor = [.9 .9 .9];

if dockFlag==1
    ui.panels.mesh = mrvPanel('right',   0.2,  ui.fig,  'normalized');
else
    hfig = figure('Name',   sprintf('Info [%s]', ui.tag),   ...
                  'Units',   'normalized',   ...
                  'Position',  [.88 .5 .12 .5],   ...
                  'MenuBar',   'none',   ...
                  'NumberTitle',   'off',   ...
                  'UserData', ui.tag);
    ui.panels.mesh = uipanel('Parent', hfig, 'Units',   'normalized',   ...
                               'Position',  [0 0 1 1]);
	% close request function for this figure: just toggle visibility
    crf = 'tmp = findobj(''Parent'', gcf, ''Type'', ''uipanel''); ';
    crf = [crf 'mrvPanelToggle(tmp); clear tmp;'];
    set(hfig, 'CloseRequestFcn', crf);
end

set(ui.panels.mesh, 'BackgroundColor', bgColor, 'Title', 'Mesh', ...
    'ShadowColor', bgColor, 'UIContextMenu', meshPanelContextMenu(ui));

% add a popup showing the selected segmentation
uicontrol('Parent', ui.panels.mesh, 'Style', 'text', 'String', 'Selected Segmentation', ...
    'Units', 'normalized', 'Position', [.1 .95 .8 .03], ...
    'HorizontalAlignment', 'left', 'BackgroundColor', [.9 .9 .9]);

ui.controls.segSelect = uicontrol('Parent', ui.panels.mesh, 'Style', 'popup',...
	'String', {'(none loaded)'}, 'Value', 1, ...
    'Units', 'normalized', 'Position', [.1 .91 .8 .04], ...
    'HorizontalAlignment', 'left', 'BackgroundColor', [1 1 1], ...
	'Callback', 'mrViewSet(gcf, ''SegmentationNum'', get(gcbo, ''Value''));');

% add a popup menu showing all meshes loaded (across all segmentations)
uicontrol('Parent', ui.panels.mesh, 'Style', 'text', 'String', 'Selected Mesh', ...
    'Units', 'normalized', 'Position', [.1 .86 .8 .03], ...
    'HorizontalAlignment', 'left', 'BackgroundColor', [.9 .9 .9]);

ui.controls.meshSelect = uicontrol('Parent', ui.panels.mesh, 'Style', 'popup',...
	'String', {'(none loaded)'}, 'Value', 1, ...
    'Units', 'normalized', 'Position', [.1 .82 .8 .04], ...
    'HorizontalAlignment', 'left', 'BackgroundColor', [1 1 1], ...
	'Callback', 'mrViewSelectMesh(gcf, get(gcbo, ''Value''));');

% update all meshes button
ui.controls.meshUpdateAll = uicontrol('Parent', ui.panels.mesh,  ...
    'Style', 'pushbutton',  'String', 'Update All', ...
    'BackgroundColor', [.9 .8 .9], ...
    'Units', 'normalized', 'Position', [.1 .78 .8 .04],   ...
    'Callback', 'mrViewUpdateAllMeshes; ');

% mesh update button
ui.controls.meshUpdate = uicontrol('Parent', ui.panels.mesh,  ...
    'Style', 'pushbutton',  'String', 'Update', ...
    'BackgroundColor', [.9 .8 .8], ...
    'Units', 'normalized', 'Position', [.1 .7 .8 .04],   ...
    'Callback', 'mrViewMesh; ');

% mesh open button
cb = 'mrViewLoad(gcf, [], ''Mesh''); mrViewDisplayMesh; ';
ui.controls.meshUpdate = uicontrol('Parent', ui.panels.mesh,  ...
    'Style', 'pushbutton',  'String', 'Open', ...
    'BackgroundColor', [.9 .8 .8], ...
    'Units', 'normalized', 'Position', [.1 .65 .4 .04],   ...
    'Callback', cb);

% mesh close button
cb = ['segCloseMesh( mrViewGet(gcf, ''Segmentation'') );'];
ui.controls.meshUpdate = uicontrol('Parent', ui.panels.mesh,  ...
    'Style', 'pushbutton',  'String', 'Close', ...
    'BackgroundColor', [.9 .8 .8], ...
    'Units', 'normalized', 'Position', [.5 .65 .4 .04],   ...
    'Callback', cb);

% mesh get ROI button
cb = ['mrViewROI(''mesh''); '];
ui.controls.meshUpdate = uicontrol('Parent', ui.panels.mesh,  ...
    'Style', 'pushbutton',  'String', 'Get ROI', ...
    'BackgroundColor', [.9 .8 .8], ...
    'Units', 'normalized', 'Position', [.1 .6 .4 .04],   ...
    'Callback', cb);

% mesh get Disk button (N.Y.I.; just does mesh ROI layer 1 xform now)
cb = ['mrViewROI(''meshgrow''); '];
ui.controls.meshUpdate = uicontrol('Parent', ui.panels.mesh,  ...
    'Style', 'pushbutton',  'String', 'Grow ROI', ...
    'BackgroundColor', [.9 .8 .8], ...
    'Units', 'normalized', 'Position', [.5 .6 .4 .04],   ...
    'Callback', cb);

% toggle cursor button
ui.controls.meshUpdate = uicontrol('Parent', ui.panels.mesh,  ...
    'Style', 'pushbutton',  'String', 'Cursor', ...
    'BackgroundColor', [.9 .8 .8], ...
    'Units', 'normalized', 'Position', [.1 .55 .4 .04],   ...
    'Callback', 'meshToggleCursor(mrViewGet(gcf, ''mesh''));'); % 'mrViewSetCursorFromMesh;'

% mesh prefs button
ui.controls.meshUpdate = uicontrol('Parent', ui.panels.mesh,  ...
    'Style', 'pushbutton',  'String', 'Prefs', ...
    'BackgroundColor', [.9 .8 .8], ...
    'Units', 'normalized', 'Position', [.5 .55 .4 .04],   ...
    'Callback', 'mrmPreferences; ');

%%%%%%%%%%%%%%%%%%%
%% Mesh Settings %%
%%%%%%%%%%%%%%%%%%%
% listbox w/ settings
mrGlobals2;
names = {''};
uicontrol('Parent', ui.panels.mesh, 'Style', 'text', 'String', 'Settings', ...
    'Units', 'normalized', 'Position', [.1 .45 .8 .05], 'FontSize', 12, ...
    'HorizontalAlignment', 'left', 'BackgroundColor', [.9 .9 .9]);

h.settingsList = uicontrol('Parent', ui.panels.mesh, 'Style', 'listbox', ...
                   'Units', 'normalized', 'Position', [.1 .2 .8 .25], ...
                   'String', names, 'Tag', 'MeshSettingsList');

% store settings button
callback = 'meshStoreSettings(mrViewGet(gcf, ''CurMesh'')); ';
h.store = uicontrol('Parent', ui.panels.mesh, 'Style', 'pushbutton', ...
            'Units', 'normalized', 'Position', [.1 .15 .4 .04], ...
            'String', 'Store', 'Callback', callback);

% retrieve settings button
callback = 'meshRetrieveSettings(mrViewGet(gcf, ''CurMesh'')); ';
h.retrieve = uicontrol('Parent', ui.panels.mesh, 'Style', 'pushbutton', ...
            'Units', 'normalized', 'Position', [.1 .1 .4 .04], ...
            'String', 'Retrieve', 'Callback', callback);

% rename settings button
callback = 'meshRenameSettings(mrViewGet(gcf, ''CurMesh'')); ';
h.rename = uicontrol('Parent', ui.panels.mesh, 'Style', 'pushbutton', ...
            'Units', 'normalized', 'Position', [.5 .15 .4 .04], ...
            'String', 'Rename', 'Callback', callback);

% delete settings button
callback = 'meshDeleteSettings(mrViewGet(gcf, ''CurMesh'')); ';
h.delete = uicontrol('Parent', ui.panels.mesh, 'Style', 'pushbutton', ...
            'Units', 'normalized', 'Position', [.5 .1 .4 .04], ...
            'String', 'Delete', 'Callback', callback);

% % manual zoom slider
% cb = ['MSH = mrViewGet(gcf, ''SelectedMesh''); ' ...
%       'TMP = mrmGet(MSH,''Camera''); TMP.actor = 0; ' ...
%       'TMP.frustum(2) = 600 - get(gcbo, ''Value''); ' ...
%       'mrMesh(MSH.host, MSH.id, ''set'', TMP); clear TMP MSH;'];
% h.zoom = uicontrol('Parent', ui.panels.mesh, 'Style', 'slider', ...
%             'Min', 0, 'Max', 600, ...
%             'Units', 'normalized', 'Position', [.1 .05 .8 .04], ...
%             'String', 'Zoom', 'Callback', cb);
% uicontrol('Parent', ui.panels.mesh, 'Style', 'text', 'Units', 'normalized', ...
%     'Position', [.1 .1 .8 .04], 'String', 'Manual Zoom', ...
%     'HorizontalAlignment', 'left');

ui.meshHandles = h;

%%%%%hide the panel if requested
if vis==0
	mrvPanelToggle(ui.panels.mesh, 'off');
end

return

function h = meshPanelContextMenu(ui);
% creates a UI Context menu for the mesh panel. This initially has a single
% option, allowing you to set the mesh width (in case it gets screwed up).
% In the future, I may add more context-specific operations.
h = uicontextmenu;

cb = ['UI = mrViewGet;  TMP_POS = get(UI.panels.mesh, ''Position''); ' ...
	  'PROMPT = {''Enter New Panel Width: ''}; ' ...
	  'DEF = {num2str(TMP_POS(3))}; ' ...
	  'TMP_RESP = inputdlg(PROMPT, ''Mesh Panel'', 1, DEF); ' ...
	  'TMP_POS(3) = str2num( TMP_RESP{1} ); ' ...
	  'set(ui.panels.mesh, ''Position'', TMP_POS); ' ...
	  'clear UI TMP_POS PROMPT DEF TMP_RESP TMP_POS; '];
uimenu(h, 'Label', 'Set Panel Width');

return