text:
%% inRange Thresholding Operations
%
% In this demo, we show how to:
%
% * Perform basic thresholding operations using OpenCV function |cv.inRange|
% * Detect an object based on the range of pixel values it has
%
% For improving detection by color, it is common to perform thresholding in
% other colorspaces like HSV or LAB, where luma and chroma are represented
% separately.
%
% Sources:
%
% * <https://docs.opencv.org/3.2.0/da/d97/tutorial_threshold_inRange.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/ImgProc/Threshold_inRange.cpp>
%

function varargout = threshold_inrange_demo_gui(im)
    % load source image
    if nargin < 1
        src = cv.imread(fullfile(mexopencv.root(),'test','monster.jpg'));
    elseif ischar(im)
        src = cv.imread(im);
    else
        src = im;
    end
    validateattributes(src, {'uint8'}, {'size',[NaN NaN 3]});

    % create the UI
    h = buildGUI(src);
    if nargout > 0, varargout{1} = h; end
end

function onType(~,e,h)
    %ONTYPE  Event handler for key press on figure

    % handle keys
    switch e.Key
        case 'h'
            helpdlg({
                'Hot keys:'
                'h - this help dialog'
                'q - quit the program'
                's - save thresholded image'
            });

        case {'q', 'escape'}
            close(h.fig);

        case {'s', 'space'}
            % save image
            img = get(h.img(2), 'CData');
            fname = fullfile(tempdir(), ...
                sprintf('out_%s.png', datestr(now(),'yyyymmddTHHMMSS')));
            cv.imwrite(fname, img);
            disp(['Saved ' fname]);
    end
end

function onChange(~,~,h)
    %ONCHANGE  Event handler for UI controls

    % retrieve current values from UI controls
    rlo = round(get(h.slid(1), 'Value'));
    rhi = round(get(h.slid(2), 'Value'));
    glo = round(get(h.slid(3), 'Value'));
    ghi = round(get(h.slid(4), 'Value'));
    blo = round(get(h.slid(5), 'Value'));
    bhi = round(get(h.slid(6), 'Value'));

    % permform thresholding
    lowerb = min([rlo glo blo], [rhi ghi bhi]);
    upperb = max([rlo glo blo], [rhi ghi bhi]);
    mask = cv.inRange(h.src, lowerb, upperb);

    % apply mask
    out = cv.copyTo(h.src, 'Mask',mask);

    % update UI and show result
    set(h.txt(1), 'String',sprintf(' Low R: %3d',rlo));
    set(h.txt(2), 'String',sprintf('High R: %3d',rhi));
    set(h.txt(3), 'String',sprintf(' Low G: %3d',glo));
    set(h.txt(4), 'String',sprintf('High G: %3d',ghi));
    set(h.txt(5), 'String',sprintf(' Low B: %3d',blo));
    set(h.txt(6), 'String',sprintf('High B: %3d',bhi));
    set(h.img(2), 'CData',out);
    drawnow;
end

function h = buildGUI(img)
    %BUILDGUI  Creates the UI

    % parameters
    sz = size(img);
    rlo = 100; rhi = 255;
    glo = 0; ghi = 100;
    blo = 0; bhi = 100;

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.src = img;
    h.fig = figure('Name','Threshold Demo', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2)*2 sz(1)+155-1]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    imgs = {img, false(sz(1:2))};
    for i=1:numel(imgs)
        h.ax(i) = axes('Parent',h.fig, 'Units','pixels', ...
            'Position',[1+sz(2)*(i-1) 155 sz(2) sz(1)]);
        if ~mexopencv.isOctave()
            h.img(i) = imshow(imgs{i}, 'Parent',h.ax(i));
        else
            %HACK: https://savannah.gnu.org/bugs/index.php?45473
            axes(h.ax(i));
            h.img(i) = imshow(imgs{i});
        end
    end
    h.txt(1) = uicontrol('Parent',h.fig, 'Style','text', ...
        'Position',[5 5 130 20], 'FontSize',11, 'ForegroundColor','r', ...
        'String',sprintf(' Low R: %3d',rlo));
    h.txt(2) = uicontrol('Parent',h.fig, 'Style','text', ...
        'Position',[5 30 130 20], 'FontSize',11, 'ForegroundColor','r', ...
        'String',sprintf('High R: %3d',rhi));
    h.txt(3) = uicontrol('Parent',h.fig, 'Style','text', ...
        'Position',[5 55 130 20], 'FontSize',11, 'ForegroundColor','g', ...
        'String',sprintf(' Low G: %3d',glo));
    h.txt(4) = uicontrol('Parent',h.fig, 'Style','text', ...
        'Position',[5 80 130 20], 'FontSize',11, 'ForegroundColor','g', ...
        'String',sprintf('High G: %3d',ghi));
    h.txt(5) = uicontrol('Parent',h.fig, 'Style','text', ...
        'Position',[5 105 130 20], 'FontSize',11, 'ForegroundColor','b', ...
        'String',sprintf(' Low B: %3d',blo));
    h.txt(6) = uicontrol('Parent',h.fig, 'Style','text', ...
        'Position',[5 130 130 20], 'FontSize',11, 'ForegroundColor','b', ...
        'String',sprintf('High B: %3d',bhi));
    h.slid(1) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Position',[135 5 sz(2)-135-5 20], 'Value',rlo, ...
        'Min',0, 'Max',255, 'SliderStep',[1 20]./(255-0));
    h.slid(2) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Position',[135 30 sz(2)-135-5 20], 'Value',rhi, ...
        'Min',0, 'Max',255, 'SliderStep',[1 20]./(255-0));
    h.slid(3) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Position',[135 55 sz(2)-135-5 20], 'Value',glo, ...
        'Min',0, 'Max',255, 'SliderStep',[1 20]./(255-0));
    h.slid(4) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Position',[135 80 sz(2)-135-5 20], 'Value',ghi, ...
        'Min',0, 'Max',255, 'SliderStep',[1 20]./(255-0));
    h.slid(5) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Position',[135 105 sz(2)-135-5 20], 'Value',blo, ...
        'Min',0, 'Max',255, 'SliderStep',[1 20]./(255-0));
    h.slid(6) = uicontrol('Parent',h.fig, 'Style','slider', ...
        'Position',[135 130 sz(2)-135-5 20], 'Value',bhi, ...
        'Min',0, 'Max',255, 'SliderStep',[1 20]./(255-0));

    % hook event handlers, and trigger default start
    opts = {'Interruptible','off', 'BusyAction','cancel'};
    set(h.slid, 'Callback',{@onChange,h}, opts{:});
    set(h.fig, 'WindowKeyPressFcn',{@onType,h}, opts{:});
    onChange([],[],h);
end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function [ y2, j2, f2 ] = ymdf_to_yjf_islamic ( y1, m1, d1, f1 )

%*****************************************************************************80
%
%% YMDF_TO_YJF_ISLAMIC converts from YMDF to YJF form in the Islamic calendar.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer Y1, M1, D1, real F1,
%    the YMDF date.
%
%    Output, integer Y2, J2, real F2, the YJF date.
%

%
%  Check the date.
%
  [ y1, m1, d1, ierror ] = ymd_check_islamic ( y1, m1, d1 );

  if ( ierror ~= 0 )
    y2 = 0;
    j2 = 0;
    f2 = 0.0;
    return
  end

  y2 = y1;
  j2 = d1;
  f2 = f1;
%
%  Add in the days of the elapsed months.
%
  for m = 1 : m1 - 1
    j2 = j2 + month_length_islamic ( y2, m );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [s,lon,lat] = m_geodesic(lon1,lat1,lon2,lat2,Npoints,varargin)
% M_GEODESIC - Compute points along geodesics of an ellipsoidal earth.
%
%       [S,LON,LAT] = m_geodesic(lon1,lat1,lon2,lat2,Npoints,spheroid)
%
% lat1 = GEODETIC latitude of first point (degrees)
% lon1 = longitude of first point (degrees)
% lat2, lon2 = second point (degrees)
% Npoints = number of points along geodesic
% spheroid = (Optional) spheroid, defaults to 'wgs84'
% S = distance in meters
% LON,LAT = points on geodesics.
%
% Note that inputs can be the same size, or a mixture of scalars and matrices.
% However, output curves will be on columns of a matrix (i.e. form of
% inputs will not be preserved).
%
% See M_FDIST, M_IDIST, M_LLDIST

% R. pawlowicz (rich@ocgy.ubc.ca) 10/Jan/2005
%
%
% This software is provided "as is" without warranty of any kind. But
% it's mine, so you can't sell it.

% Do the equivalent of meshgrid-type calls to make all
% matrix sizes match. To do this we can have no more than
% two different numbers in any particular dimension, and one
% of these has to be a '1'.
allsize=[size(lon1);size(lat1);size(lon2);size(lat2)];
i1=ones(1,size(allsize,2));
for k=1:size(allsize,2)
 rs=unique(allsize(:,k));
 if length(rs)==2 && rs(1)==1
   j1=i1;j1(k)=rs(2);
   if allsize(1,k)==1,lon1=repmat(lon1,j1); end
   if allsize(2,k)==1,lat1=repmat(lat1,j1); end
   if allsize(3,k)==1,lon2=repmat(lon2,j1); end
   if allsize(4,k)==1,lat2=repmat(lat2,j1); end
 elseif length(rs)>2
  error('incompatible array sizes!');
 end
end

% Get the distances and bearings

[S,A12,A21]=m_idist(lon1,lat1,lon2,lat2,varargin{:});

s=[0:Npoints-1]'/(Npoints-1)*S(:)';

[lon,lat]=m_fdist(lon1(:)',lat1(:)',A12(:)',s,varargin{:});

author: goGPS-Project


text:
function opa=decimate_optimums(z,op,opv)
% z  zeros
% op optimums from d/dt=0
% opv optimums values
Lz=length(z);
Lop=length(op);
zop=[z; op];
zopv=[zeros(size(z)); opv]; % values
[zops ind]=sort(zop);
zopsv=zopv(ind);
isz=(ind<=Lz); % if zero in zops list
iszind=find(isz);
L1=length(iszind);
opa=[];
opva=[];
wasused=false(size(zops)); % what optimums was already used
for c=1:L1-1
    ind1=iszind(c);
    z1=zops(ind1);
%     if z1>1.47e4
%         'here'
%     end
    ind2=iszind(c+1);
    z2=zops(ind2);
    %ii=ind1+1:ind2-1;

    if (ind2+1)<=length(isz)
        inOnePlace=(~isz(ind2+1))&&(zops(ind2)==zops(ind2+1));
    else
        inOnePlace=false;
    end

    if inOnePlace
        % optimum and zero in one place
        ii=ind1+1:ind2+1;
    else
        ii=ind1+1:ind2;
    end
    wasused1=wasused(ii);

    ii=ii(~wasused1);

%     isz1=isz(ii);
%     isz2=~isz1; % what is optimums
%     isz3=isz2&(~wasused1); % optimums that was not used
%     %if all(isz(ii))
%     if any(isz3) % if any optimums that was not used
%
%     else
%         % case when zeros and maximum are same
%         ii=ind1+1:ind2+1;
%     end
    if ~isempty(ii)
        opt=zops(ii); % optimums between z1 and z2
        optv=zopsv(ii); % values
        [mx mxi]=max(abs(optv));
        opa=[opa; opt(mxi)];
        opva=[opva; optv(mxi)];
        wasused(ii(mxi))=true;
    end
end

mdf=0.005;

% detete optimums that to smaller in compare to bigest optimum:
bo=max(abs(opva));
indo=find(abs(opva)>=mdf*bo);
opa=opa(indo);

% delete repeated points:
opa=unique(opa);

author: Sable

repo: mcbench-benchmarks


text:
function [maskOut]=kGaussian_color_EM(imageFile,k)
%% this function uses EM algorithm and gaussian distribution function to do
% image segmentation on the color space
% try:
%   img=imread('onion.png');
%   [maskOut]=kGaussian_color_EM(img,9); figure;  imshow(maskOut)
% written by Rongwen Lu on 12/11/2011; email: rongwen@uab.edu
img=double(imageFile);

[M,N,P]=size(img);
n=M*N;
imgR=img(:,:,1);
imgG=img(:,:,2);
imgB=img(:,:,3);
[cy,cx]=ind2sub([M,N],1:n);

%normalize each vector; delete it if normalization is not needed; weights
%are also assigned here;
% imgR=mat2gray(imgR);
% imgG=mat2gray(imgG);
% imgB=mat2gray(imgB);
% cy=mat2gray(cy);
% cx=mat2gray(cx);
imgR=imgR/255;
imgG=imgG/255;
imgB=imgB/255;
cy=cy/M;
cx=cx/N;

% %% Gaussian filter
% w=fspecial('gaussian',[5,5],4);
% imgR=imfilter(imgR,w);
% imgG=imfilter(imgG,w);
% imgB=imfilter(imgB,w);

%% assign vectors into the matrix raw
raw=zeros(n,3);
raw(:,1)=imgR(:);
raw(:,2)=imgG(:);
raw(:,3)=imgB(:);

% raw=zeros(n,5);
% raw(:,1)=cy.';
% raw(:,2)=cx.';
% raw(:,3)=imgR(:);
% raw(:,4)=imgG(:);
% raw(:,5)=imgB(:);

%% get assignment matrix p, which is also memebership probability here; u
%% is vector of the estimated means of Gaussian function, v is the vector ot the estimated SD
[p,u,v]=em(raw,k);

imgRe=zeros(n,3);
kColor=jet(k);
kColor=u(:,1:3);
imgRe=p*kColor;

maskOut=zeros(M,N,3);
for ii=1:3
    maskOut(:,:,ii)=reshape(imgRe(:,ii),[M,N]);
end

% figure; imshow((maskOut))
% title('based on k Gaussian by EM algorithm on color space')

 function [p,u,v]=em(raw,k)
 % this function uses Expectation Maximization method to estimate k
 % Gaussian distribution functions. There are 2 input arguments. Raw's size
 % is [n,dim], where n is the NO of data and dim is the dimention of data.
 % k means using k Gaussin functions to do image segmentation.The output is
 % p, which is assignment matrix.its [ii,jj]th element means the
 % probability that Xii is generated by jjth Gaussian function.

%% the input of following codes are raw,and K

%% initialize centroid matrix u; u has k raws. NO of column is the same as
%% raw
[n,dim]=size(raw);
u=raw(randint(k,1,[1,n]),:);
%u=[0,0;40,40;-40,-40];

%% initialize standard diviation v, v has k raws, one column.
v=zeros(k,1);
for ii=1:k
    raw_tmp=raw(ii:k:end,1);
    v(ii,:)=std(raw_tmp);
end
%v=.5*[4,4,4]';
%% initialize weight w
w=ones(k,1)/k;

%% initialize membership probability matrix (assignment matrix) p, which is p(k|x)
p=zeros(n,k);

%% do interation to get best r
u0=u*0;
v0=0*v;
w0=w*0;
energy=sum(sum((u-u0).^2))+sum(sum((v-v0).^2))+(sum((w-w0).^2));
iteration=1;
x_u=zeros(size(raw));
while energy>10^(-6)

    %% calculate membership probability, which is also assignment matrix
    for jj=1:k
        for ss=1:dim
            x_u(:,ss)=raw(:,ss)-u(jj,ss)*ones(n,1);
        end
        x_u=x_u.*x_u;
        p(:,jj)=power(sqrt(2*pi)*v(jj),-1*dim)*exp((-1/2)*sum(x_u,2)./(v(jj).^2));
        p(:,jj)=p(:,jj)*w(jj);

    end
    %% normalize p on the x dimention
    pSum=sum(p,2);
    for jj=1:k
        p(:,jj)=p(:,jj)./pSum;
    end

    %% normlaize p on the y dimention, yielding pNorm
    pSum2=sum(p,1);
    pNorm=p*0;
    for jj=1:k
        pNorm(:,jj)=p(:,jj)/pSum2(jj);
    end

    %% save current u, v, and w as u0, v0 and w0
    u0=u;v0=v;w0=w;

    %%update u
    u=(pNorm.')*raw;

    %% update v
    for jj=1:k
        for ss=1:dim
            x_u(:,ss)=raw(:,ss)-u(jj,ss)*ones(n,1);
        end
        x_u=x_u.*x_u;
        x_uSum=sum(x_u,2);
        v(jj)=sqrt(1/dim*(pNorm(:,jj).')*x_uSum);
    end

    %% update w
    w=(sum(p)/n).';

    %% update and display iteration and energy

    %disp(sprintf(['iteration=',num2str(iteration),'; energy=',num2str(energy,'%g')]))
    iteration=iteration+1;
    energy=sum(sum((u-u0).^2))+sum(sum((v-v0).^2))+(sum((w-w0).^2));

end


text:
function [xi] = log_multiSE3(chi)
C = chi(1:3,1:3);
r = chi(1:3,4:end);

[vecs,eigs] = eig(C);
[~,idx] = min(abs(diag(eigs)-1));
a = vecs(:,idx);

phi = acos(1/2*(trace(C)-1));
if phi == 0
    a = zeros(3,1);
    iJ = eye(3);
else
    A = [0 -a(3) a(2);a(3) 0 -a(1);-a(2) a(1) 0];
    a_ta = a*transpose(a);
    C1 = cos(phi)*eye(3) + (1-cos(phi))*a_ta + sin(phi)*A;
    C2 = cos(-phi)*eye(3) + (1-cos(-phi))*a_ta + sin(-phi)*A;
    if norm(C2-C)<norm(C1-C)
        phi = -phi;
    end
    iJ = phi/2*cot(phi/2)*eye(3) + (1-phi/2*cot(phi/2))*a_ta - phi/2*A;
end
rho = iJ*r;
xi = real([phi*a;rho(:)]);
end

author: mbrossar


text:
function L = laplacian(f)
%LAPLACIAN   Laplacian of a SEPARABLEAPPROX.
%   L = LAPLACIAN(F) returns a SEPARABLEAPPROX representing the Laplacian of F.
%
% See also SEPARABLEAPPROX/LAP.

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

% laplacian(f) = f_xx + f_yy:
L = diff(f, 2, 2) + diff(f, 2, 1);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function pde = sincosRobindata
%% SINCOSDATA trigonometric  data for Poisson equation
%
%     f = 8*pi^2*sin(2*pi*x)*cos(2*pi*y);
%     u = sin(2*pi*x)*cos(2*pi*y);
%     Du = (2*pi*cos(2*pi*x)*cos(2*pi*y), -2*pi*sin(2*pi*x)*sin(2*pi*y));
%
%  Impose Robin boundary condition g_R*u+ du/dn = g_N on [0,1]^2.
%
% Copyright (C)  Long Chen. See COPYRIGHT.txt for details.

pde = struct('f',@f,'exactu',@exactu,'Du',@Du,'g_D',@g_D,'g_R',@g_R,'g_N',@g_N);

% load data (right hand side function)
function rhs =  f(p)
    x = p(:,1); y = p(:,2);
    rhs =  8*pi^2*sin(2*pi*x).*cos(2*pi*y);
end
% exact solution
function u =  exactu(p)
    x = p(:,1); y = p(:,2);
    u =  sin(2*pi*x).*cos(2*pi*y);
end
% derivative of the exact solution
function uprime =  Du(p)
    x = p(:,1); y = p(:,2);
    uprime(:,1) = 2*pi*cos(2*pi*x).*cos(2*pi*y);
    uprime(:,2) = -2*pi*sin(2*pi*x).*sin(2*pi*y);
end
% Dirichlet boundary condition
function u =  g_D(p)
    u =  exactu(p);
end
% coefficients of Robin boundary condition g_R*u+ du/dn = g_N
function f = g_R(p)
    f = ones(size(p,1),1);
end
%
function f = g_N(p)
    f = zeros(size(p,1),1);
    x = p(:,1); y = p(:,2);
    uprime = Du(p);
    leftbd = (abs(x)<eps);  % n = (-1,0);
    f(leftbd) = g_R(p(leftbd,:)).*exactu(p(leftbd,:)) - uprime(leftbd,1);
    rightbd = (abs(x-1)<eps); % n = (1,0);
    f(rightbd) = g_R(p(rightbd,:)).*exactu(p(rightbd,:)) + uprime(rightbd,1);
    topbd = (abs(y-1)<eps);   % n = (0,1)
    f(topbd) = g_R(p(topbd,:)).*exactu(p(topbd,:)) + uprime(topbd,2);
    bottombd = (abs(y)<eps);% n = (0,-1)
    f(bottombd) = g_R(p(bottombd,:)).*exactu(p(bottombd,:)) - uprime(bottombd,2);
end
end

repo: ifem


text:
% clc; clear all; close all;

%% hydro data
hydro = struct();
hydro = readAQWA(hydro, 'cubes.AH1', 'cubes.LIS');
hydro = radiationIRF(hydro,40,[],[],[],[]);
hydro = radiationIRFSS(hydro,[],[]);
hydro = excitationIRF(hydro,50,[],[],[],[]);
writeBEMIOH5(hydro)

%% Plot hydro data
% plotBEMIO(hydro)

author: WEC-Sim

repo: WEC-Sim

save_path: github-repos/MATLAB/WEC-Sim-WEC-Sim


text:
% This subroutine assigns creates a 3D grid with
% spacing dx,dy, dz (in degreees). The size will
% be selected interactiVELY. The pvalue in each
% volume around a grid point containing ni earthquakes
% will be calculated as well as the magnitude
% of completness
%   Stefan Wiemer 1/98

report_this_filefun(mfilename('fullpath'));
global no1 bo1 inb1 inb2

if sel == 'in'
    % get the grid parameter
    % initial values
    %
    dx = 0.1;
    dy = 0.1 ;
    dz = 5.00 ;

    def = {'0.1','0.1',num2str(dz),num2str(max(a.Depth)), num2str(min(a.Depth))};

    tit ='Three dimesional b-value analysis';
    prompt={ 'Spacing in Lat/Lon (dx in [deg])',...
        'Sample Raduis [km])',...
        'Spacing in Depth    (dz in [km ])',...
        'Depth Range: deep limit [km] ',...
        'Depth Range: shallow limit',...
        };

    ni2 = inputdlg(prompt,tit,1,def);

    l = ni2{1}; dx= str2num(l);dy = dx;
    l = ni2{2}; R= str2double(l);
    l = ni2{3}; dz= str2double(l);
    l = ni2{4}; z1= str2double(l);
    l = ni2{5}; z2= str2double(l);

    sel = 'ca'; density_3D

end   % if sel == 'in'

% get the grid-size interactively and
% calculate the b-value in the grid by sorting
% thge seimicity and selectiong the ni neighbors
% to each grid point

if sel == 'ca'
    selgp3dB

    gz = zvect;
    itotal = length(t5);
    zmap_message_center.set_info(' ','Running... ');think
    %  make grid, calculate start- endtime etc.  ...
    %
    bvg = ones(length(gx),length(gy),length(gz))*nan;

    t0b = min(a.Date)  ;
    n = a.Count;
    teb = a(n,3) ;
    tdiff = round((teb - t0b)*365/par1);
    loc = zeros(3, length(gx)*length(gy));
    Rconst = R;
    % loop over  all points
    %
    i2 = 0.;
    i1 = 0.;
    allcount = 0.;
    wai = waitbar(0,' Please Wait ...  ');
    set(wai,'NumberTitle','off','Name',' 3D gridding - percent done');;
    drawnow
    %
    %

    z0 = 0; x0 = 0; y0 = 0; dt = 1;
    % loop over all points
    for il =1:length(t5)

        x = t5(il,1);
        y = t5(il,2);
        z = t5(il,3);

        allcount = allcount + 1.;

        % calculate distance from center point and sort wrt distance
        l = sqrt(((a.Longitude-x)*cos(pi/180*y)*111).^2 + ((a.Latitude-y)*111).^2 + ((a.Depth - z)).^2 ) ;
        [s,is] = sort(l);
        b = a(is(:,1),:) ;       % re-orders matrix to agree row-wise

        l3 = l <= R;
        b = a.subset(l3);      % new data per grid point (b) is sorted in distanc
        rd = length(b(:,1));

        bvg(t5(il,5),t5(il,6),t5(il,7)) = rd;

        waitbar(allcount/itotal)
    end  % for t5

    % save data
    %
    gz = -gz;
    zv2 = bvg;
    zvg = bvg;

    close(wai)
    watchoff

    sel = 'no';

    ButtonName=questdlg('Which viwer would you like to use?', ...
        'Question', ...
        'Slicer - map view','Slicer - 3D ','Help','none');

    switch ButtonName
        case 'Slicer - map view'
            slm = 'new'; slicemap;
        case 'Slicer - 3D '
            ac2 = 'new'; myslicer;
        case 'Help'
            showweb('3dbgrids')
    end % switch

    uicontrol('Units','normal',...
        'Position',[.90 .95 .04 .04],'String','Slicer',...
         'Callback','')

    %ac2 = 'new'; myslicer;

end  % if cal

author: CelsoReyes


text:
function outputImage = perform_binary_operation(this, otherImage, ...
    functionHandle)
% Performs binary operation (i.e. 2 inputs) on this and a second image
% for given dimension and returns a new output image
%
% NOTE: If voxel dimensions of 2nd image do not match - or a scalar is
% given as 2nd argument - data in the 2nd argument is automatically
% replicated to match this image geometry.
%
%   Y = MrDataNd()
%   outputImage = perform_binary_operation(this, otherImage, ...
%   functionHandle)
%
% This is a method of class MrDataNd.
%
% IN
%   otherImage              2nd operand for binary operation
%   functionHandle          handle of function to be applied to images (e.g.
%                           @plus, @minus, @times and @rdivide for the 4
%                           arithmetic operations )
%
% OUT
%   outputImage             new MrDataNd with possibly new image dimensions,
%                           output of binary operation performed on this
%                           and otherImage
% EXAMPLE
%
%   % Compute difference of 2 images
%		Y = MrDataNd();
%		Z = MrDataNd();
%		X = Y.perform_binary_operation(Z, @minus);
%
%	% Scale image (multiply) by a factor of 3
%		Y = MrDataNd();
%		Y = Y.perform_binary_operation(3, @mult)
%
%	% Compute ratio of 2 images
%		Y 			= MrDataNd();
%		Z 			= MrDataNd();
%		ratioYZ 	= Y.perform_binary_operation(Z, @rdivide);
%
%
%
%   See also MrDataNd perform_unary_operation

% Author:   Saskia Bollmann & Lars Kasper
% Created:  2014-11-13
% Copyright (C) 2014 Institute for Biomedical Engineering
%                    University of Zurich and ETH Zurich
%
% This file is part of the TAPAS UniQC Toolbox, which is released
% under the terms of the GNU General Public Licence (GPL), version 3.
% You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version).
% For further details, see the file COPYING or
%  <http://www.gnu.org/licenses/>.
%

% binary operation with scalar etc. => make MrDataNd first
if ~isa(otherImage, 'MrDataNd')
    otherData = otherImage;
    otherName = 'dataMatrix';
    % Create image from data matrix with similar dimInfo
    % TODO: include proper dimInfo-adaptation (right now, res is assumed
    % the same, not FOV!)
    nSamplesOther = size(otherData);
    nSamplesOther(end+1:this.dimInfo.nDims) = 1; % avoid singleton dim error!
    dimInfoOther = this.dimInfo.copyobj;
    dimInfoOther.nSamples = nSamplesOther;

    otherImage = this.copyobj();
    otherImage.dimInfo = dimInfoOther;
    otherImage.data = otherData;
    otherImage.name = otherName;
end

%% TODO: FOV first, then adapt matrix sizes?
% TODO: Check FOVs first, if they don't match crop or zero-fill otherImage

outputImage = otherImage.resize(this.dimInfo);
outputImage.name = this.name;

%% Perform computation and store!
outputImage.data 	= functionHandle(this.data, outputImage.data);

outputImage.info{end+1,1} = sprintf('%s( %s, %s )', func2str(functionHandle), ...
    outputImage.name, otherImage.name);

author: translationalneuromodeling

repo: tapas


text:
function dataS = populate_planC_field(cellName, dcmdir_patient, optS, varargin)
%"populate_planC_field"
%   Given the name of a child cell to planC, such as 'scan', 'dose',
%   'comment', etc. return a copy of that cell with all fields properly
%   populated with data from the files contained in a dcmdir.PATIENT
%   structure.
%
%JRA 06/15/06
%YWU Modified 03/01/08
%NAV 07/19/16 updated to dcm4che3
%   replaced dcm2ml_element with getTagValue
%   and used getValue instead of get
%
%Usage:
%   dataS = populate_planC_field(cellName, dcmdir);
%
% 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/>.

%Get template for the requested cell.
persistent rtPlans scanOriS
structS = initializeCERR(cellName);
names   = fieldnames(structS);

dataS = [];
switch cellName
    case 'header'
        rtPlans = []; % Clear persistent object
        for i = 1:length(names)
            dataS.(names{i}) = populate_planC_header_field(names{i}, dcmdir_patient);
        end

    case 'scan'

        disp('Loading Scan ...')

        % supportedModalities = {'CT'};
        scansAdded = 0;

        %Extract all series contained in this patient.
        [seriesC, typeC] = extract_all_series(dcmdir_patient);
        % Extract CT, MR, PET, NM... in that order, since frame of
        % reference is reused
        ctIndV = strcmpi(typeC,'CT');
        mrIndV = strcmpi(typeC,'MR');
        newIndV = [];
        newIndV = [newIndV, find(ctIndV), find(mrIndV), find(~(ctIndV | mrIndV))];
        seriesC = seriesC(newIndV);
        typeC = typeC(newIndV);

        %ctSeries = length(find(seriesC(strcmpi(typeC, 'CT'))==1));

        %Initialize structure to store image orientation per scan.
        scanOriS = struct();
        %Place each series (CT, MR, etc.) into its own array element.
        for seriesNum = 1:length(seriesC)
            % Test IOP here to find if it is "nominal" or "non-nominal"
            % % %             outIOP = getTest_Scan_IOP(seriesC{seriesNum}.Data(1).file);
            if ismember(typeC{seriesNum},{'CT','OT','NM','MR','PT','ST','MG','SM'})

                %Populate each field in the scan structure.
                for i = 1:length(names)
                    dataS(scansAdded+1).(names{i}) = ...
                        populate_planC_scan_field(names{i}, ...
                        seriesC{seriesNum}, typeC{seriesNum}, seriesNum, optS);
                end

                if isempty(dataS(scansAdded+1).scanInfo(1).rescaleIntercept)
                    dataS(scansAdded+1).scanInfo(1).rescaleIntercept = 0;
                end

                %Apply ReScale Intercept and Slope
                scanArray3M = zeros(size(dataS(scansAdded+1).scanArray),'single');
                numSlcs = size(dataS(scansAdded+1).scanArray,3);
                rescaleSlopeV = ones(numSlcs,1);
                realWorldImageFlag = false;
                for slcNum = 1:numSlcs
                    rescaleSlope = dataS(scansAdded+1).scanInfo(slcNum).rescaleSlope;
                    rescaleIntrcpt = dataS(scansAdded+1).scanInfo(slcNum).rescaleIntercept;
                    realWorldValueSlope = dataS(scansAdded+1).scanInfo(slcNum).realWorldValueSlope;
                    realWorldValueIntercept = dataS(scansAdded+1).scanInfo(slcNum).realWorldValueIntercept;
                    realWorldMeasurCodeMeaning = dataS(scansAdded+1).scanInfo(slcNum).realWorldMeasurCodeMeaning;
                    if ~isempty(realWorldValueSlope) && strcmpi(typeC{seriesNum}, 'MR') && ...
                            ~strcmpi(realWorldMeasurCodeMeaning,'no units')
                        realWorldImageFlag = true;
                        scanArray3M(:,:,slcNum) = ...
                            single(dataS(scansAdded+1).scanArray(:,:,slcNum)) * single(realWorldValueSlope) + single(realWorldValueIntercept);
                    else
                        scanArray3M(:,:,slcNum) = ...
                            single(dataS(scansAdded+1).scanArray(:,:,slcNum)) * single(rescaleSlope) + single(rescaleIntrcpt);
                    end
                    rescaleSlopeV(slcNum) = rescaleSlope;
                end
                minScanVal = min(scanArray3M(:));
                ctOffset = max(0,-minScanVal);
                scanArray3M = scanArray3M + ctOffset;
                minScanVal = min(scanArray3M(:));
                maxScanVal = max(scanArray3M(:));
                if ~realWorldImageFlag && ~any(abs(rescaleSlopeV-1) > eps*1e5) % convert to uint if rescale slope is not 1.
                   if minScanVal >= -32768 && maxScanVal <= 32767
                       scanArray3M = uint16(scanArray3M);
                   else
                       scanArray3M = uint32(scanArray3M);
                   end
                end
                for slcNum = 1:numSlcs
                    dataS(scansAdded+1).scanInfo(slcNum).CTOffset = ctOffset;
                end
                dataS(scansAdded+1).scanArray = scanArray3M;

                %Apply scale slope & intercept for Philips data if not
                %realWorldValue
                if strcmpi(typeC{seriesNum}, 'MR') && ...
                        strcmpi(optS.store_scan_as_mr_philips_precise_value, 'yes')
                    % Ref: Chenevert, Thomas L., et al. "Errors in quantitative image analysis due to platform-dependent image scaling."
                    manufacturer = dataS(scansAdded+1).scanInfo(1).manufacturer;
                    if ~isempty(strfind(lower(manufacturer),'philips')) && ...
                            ~isempty(dataS(scansAdded+1).scanInfo(1).scaleSlope) && ...
                            ~realWorldImageFlag
                        scaleSlope = dataS(scansAdded+1).scanInfo(1).scaleSlope;
                        dataS(scansAdded+1).scanArray = single(dataS(scansAdded+1).scanArray)./(rescaleSlope*scaleSlope);
                    end
                end

                scansAdded = scansAdded + 1;

            elseif strcmpi(typeC{seriesNum}, 'US')
                %Populate each field in the scan structure.
                for i = 1:length(names)
                    dataS(scansAdded+1).(names{i}) = populate_planC_USscan_field(names{i}, seriesC{seriesNum}, typeC{seriesNum});
                end
                scansAdded = scansAdded + 1;
            end
            scanOriS(scansAdded).scanUID = dataS(scansAdded).scanUID;
            scanOriS(scansAdded).imageOrientationPatient = dataS(scansAdded).scanInfo(1).imageOrientationPatient;
        end

    case 'structures'
        [seriesC, typeC]    = extract_all_series(dcmdir_patient);
        structsAdded          = 0;

        % Get scan object from varargin if it is not empty.
        % for inserting structures to existing planC
        if ~isempty(varargin)
            scanOriS = varargin{1};
        end

        %hWaitbar = waitbar(0,'Loading Structures Please wait...');

        %strSeries = length(find(strcmpi(typeC, 'RTSTRUCT')==1));

        disp('Loading RTSTRUCT ...')

        %Place each structure into its own array element.
        for seriesNum = 1:length(seriesC)
            %if ismember(typeC{seriesNum}, supportedTypes)
            if strcmpi(typeC{seriesNum}, 'RTSTRUCT')

                RTSTRUCT = seriesC{seriesNum}.Data;
                for k = 1:length(RTSTRUCT)
                    strobj  = scanfile_mldcm(RTSTRUCT(k).file);

                    % StructureSetROISequence
                    SSRS = strobj.getValue(805699616); %org.dcm4che3.data.Tag.StructureSetROISequence; %hex2dec('30060020')
                    numSSRS = SSRS.size();
                    roiNummV = zeros(1,numSSRS,'int32');
                    for iStr = 1:numSSRS
                        ssObj = SSRS.get(iStr - 1);
                        roiNummV(iStr) = ssObj.getInts(805699618); %org.dcm4che3.data.Tag.ROINumber; % vr=IS
                    end

                    %ROI Contour Sequence.
                    RCS = strobj.getValue(805699641); %org.dcm4che3.data.Tag.ROIContourSequence; % hex2dec('30060039')

                    %If non-empty sequence, get item count, else set to zero
                    if ~isempty(RCS)
                        nStructures = RCS.size();
                    else
                        nStructures = 0;
                    end

                    % ROI = strobj.getInt(org.dcm4che2.data.Tag.ROIContourSequence);

                    structuresToImportC = optS.structuresToImport;

                    structuresImportMatchCriteria = optS.structuresImportMatchCriteria;

                    %curStructNum = 1; %wy modified for suppport multiple RS files
                    for structNum = 1:nStructures

                        cObj = RCS.get(structNum - 1);

                        %Referenced ROI Number
                        %RRN = getTagValue(cObj, '30060084');
                        RRN = cObj.getInts(805699716); %org.dcm4che3.data.Tag.ReferencedROINumber; %IS

                        indSSRS = find(ismember(roiNummV,RRN));

                        % Get Structure name
                        ssObj = SSRS.get(indSSRS - 1);
                        %structureName = getTagValue(ssObj, '30060026');
                        structureName = char(ssObj.getStrings(805699622)); %org.dcm4che3.data.Tag.ROIName;

                        if ~isempty(structuresToImportC)
                            matchIndex = getMatchingIndex(structureName,...
                                structuresToImportC,...
                                structuresImportMatchCriteria);
                        else
                            matchIndex = 1;
                        end

                        if isempty(matchIndex)
                            continue;
                        end

                        %Populate each field in the structure field set
                        for i = 1:length(names)
                            dataS(structsAdded+1).(names{i}) = ...
                                populate_planC_structures_field(names{i}, ...
                                RTSTRUCT, scanOriS, strobj, ssObj, cObj, optS);
                        end
                        %curStructNum = curStructNum + 1;
                        structsAdded = structsAdded + 1;

                        %waitbar(structsAdded/(nStructures*length(RTSTRUCT)*strSeries), hWaitbar, 'Loading structures, Please wait...');

                        %a temporary limit of 52 structs
                        %if (structsAdded>=52)
                        %    return;
                        %end

                    end
                end

            end

        end
        %close(hWaitbar);
        %pause(0.1);

        %     case 'structureArray'
        %         populate_planC_structureArray_field(fieldName, dcmdir);

    case 'dose'
        [seriesC, typeC]    = extract_all_series(dcmdir_patient);
        dosesAdded          = 0;
        frameOfRefUIDC       = {};
        %Place each RTDOSE into its own array element.
        for seriesNum = 1:length(seriesC)
            %             if ismember(typeC{seriesNum}, supportedTypes)
            if strcmpi(typeC{seriesNum}, 'RTDOSE')

                RTDOSE = seriesC{seriesNum}.Data; %wy RTDOSE{1} for import more than one dose files;
                for doseNum = 1:length(RTDOSE)
                    doseobj  = scanfile_mldcm(RTDOSE(doseNum).file);

                    % Frame of Reference UID
                    %frameOfRefUID = getTagValue(doseobj, '00200052');
                    frameOfRefUID = char(doseobj.getStrings(2097234)); %org.dcm4che3.data.Tag.FrameOfReferenceUID;

                    %check if it is a DVH
                    dvhsequence = populate_planC_dose_field('dvhsequence', ...
                        RTDOSE(doseNum), doseobj, rtPlans, optS);

                    %Check if doesArray is present
                    dose3M = populate_planC_dose_field('doseArray', ...
                        RTDOSE(doseNum), doseobj, rtPlans, optS);

                    if isempty(dvhsequence) || ~isempty(dose3M)
                        %Populate each field in the dose structure.
                        for i = 1:length(names)
                            dataSeriesS(doseNum).(names{i}) = ...
                                populate_planC_dose_field(names{i}, ...
                                RTDOSE(doseNum), doseobj, rtPlans, optS);
                        end

                        if isempty(frameOfRefUIDC)
                            frameOfRefUIDC{1} = frameOfRefUID;
                        else
                            frameOfRefUIDC{end+1} = frameOfRefUID;
                        end

                        % Check if frame of reference UID matches any
                        % existing doses and add slices.

                        % If frame of reference UID does not match, then
                        % create a new dose
                        %dosesAdded = dosesAdded + 1;
                    end
                end

                multiFrameC = cellfun(@(x) x > 1,{dataSeriesS.numberMultiFrameImages},'UniformOutput',false);
                multiFrameIndV = [];
                for iM = 1:length(multiFrameC)
                    if isempty(multiFrameC{iM}) || (~isempty(multiFrameC{iM}) && multiFrameC{iM}==0)
                        multiFrameIndV(iM) = 0;
                    else
                        multiFrameIndV(iM) = 1;
                    end
                end
                multiFrameDoseNumsV = find(multiFrameIndV);
                for doseNum = multiFrameDoseNumsV
                    dosesAdded = dosesAdded + 1;
                    if isempty(dataS)
                        dataS = dataSeriesS(doseNum);
                    else
                        dataS(dosesAdded) = dataSeriesS(doseNum);
                    end
                end

                singleFrameIndV = ~multiFrameIndV;
                dataSeriesS = dataSeriesS(singleFrameIndV);
                doseSummationTypes = unique({dataSeriesS.doseSummationType});
                for doseTypeNum = 1:length(doseSummationTypes)
                    doseSumType = doseSummationTypes{doseTypeNum};
                    indV = strcmpi(doseSumType,{dataSeriesS.doseSummationType});
                    doseSeries = dataSeriesS(indV);
                    switch doseSumType
                        case 'BEAM'
                            beamNumberV = unique([doseSeries.refBeamNumber]);
                            for iB = 1:length(beamNumberV)
                                beamNum = beamNumberV(iB);
                                indV = [doseSeries.refBeamNumber] == beamNum;
                                doseS = doseSeries(indV);
                                [dataS,dosesAdded] = addDoseToSeries(doseS,dataS,dosesAdded);
                            end
                        case 'FRACTION'
                            fractionNumberV = unique(doseSeries.refFractionGroupNumber);
                            for iF = 1:length(fractionNumberV)
                                fractionNumber = fractionNumberV(iF);
                                indV = [doseSeries.refFractionGroupNumber] == fractionNumber;
                                doseS = doseSeries(indV);
                                [dataS,dosesAdded] = addDoseToSeries(doseS,dataS,dosesAdded);
                            end
                        otherwise
                            doseS = doseSeries;
                            [dataS,dosesAdded] = addDoseToSeries(doseS,dataS,dosesAdded);

                    end

                end

            end

        end

        % Get doses that have same frame of reference UIDs
        %[uniqFrameOfRefUIDC,jnk,uniqIndV] = unique(frameOfRefUIDC);
        %dataS(find(uniqIndV==1))

    case 'DVH'
        [seriesC, typeC]    = extract_all_series(dcmdir_patient);
        supportedTypes      = {'DVH'};
        dvhsAdded           = 0;

        %Place each RTDOSE into its own array element.
        for seriesNum = 1:length(seriesC)

            %             if ismember(typeC{seriesNum}, supportedTypes)
            if strcmpi(typeC{seriesNum}, 'RTDOSE')

                RTDOSE = seriesC{seriesNum}.Data; %wy RTDOSE{1} for import more than one dose files;
                for doseNum = 1:length(RTDOSE)
                    doseobj  = scanfile_mldcm(RTDOSE(doseNum).file);

                    %check if it is a DVH
                    dvhsequence = populate_planC_dose_field('dvhsequence', ...
                        RTDOSE(doseNum), doseobj, rtPlans, optS);

                    if ~isempty(dvhsequence) && ~dvhsequence.isEmpty

                        structureNameC = {};
                        structureNumberV = [];

                        % get a list of Structure Names
                        for seriesNumStr = 1:length(seriesC)

                            if strcmpi(typeC{seriesNumStr}, 'RTSTRUCT')

                                RTSTRUCT = seriesC{seriesNumStr}.Data;
                                for k = 1:length(RTSTRUCT)
                                    strobj  = scanfile_mldcm(RTSTRUCT(k).file);

                                    % StructureSetROISequence
                                    SSRS = strobj.getValue(805699616); %org.dcm4che3.data.Tag.StructureSetROISequence; %hex2dec('30060020');
                                    numSSRS = SSRS.size();
                                    for js = 1:numSSRS
                                        % Get Structure name
                                        ssObj = SSRS.get(js - 1);

                                        %structureName = getTagValue(ssObj, '30060026');
                                        structureNameC{end+1} = char(ssObj.getStrings(805699622)); %org.dcm4che3.data.Tag.ROIName;
                                        structureNumberV(end+1) = ssObj.getInts(805699618); %org.dcm4che3.data.Tag.ROINumber;

                                    end
                                end

                            end

                        end

                        %DVH_items = fieldnames(dvhsequence);
                        DVH_items = dvhsequence.size();
                        for i = 1:DVH_items
                            dvhsAdded = dvhsAdded + 1;
                            for j = 1:length(names)
                                dataS(dvhsAdded).(names{j}) = populate_planC_DVH_field(names{j}, RTDOSE(doseNum), doseobj, rtPlans);
                            end
                            dvhObj = dvhsequence.get(i-1);
                            dataS(dvhsAdded).volumeType = char(dvhObj.getStrings(805568513)); %org.dcm4che3.data.Tag.DVHType;
                            dataS(dvhsAdded).doseType = char(dvhObj.getStrings(805568516)); %org.dcm4che3.data.Tag.DoseType;
                            dataS(dvhsAdded).doseUnits = char(dvhObj.getStrings(805568514)); %org.dcm4che3.data.Tag.DoseUnits;

                            %dataS(dvhsAdded+1).volumeType = dvhsequence.(['Item_',num2str(i)]).DVHType;
                            %dataS(dvhsAdded+1).doseType = dvhsequence.(['Item_',num2str(i)]).DoseType;
                            %dataS(dvhsAdded+1).doseUnits = dvhsequence.(['Item_',num2str(i)]).DoseUnits;

                            dvhReferencedROISequence = dvhObj.getValue(805568608); %org.dcm4che3.data.Tag.DVHReferencedROISequence;
                            if ~isempty(dvhReferencedROISequence) %isfield(dvhsequence.(['Item_',num2str(i)]),'DVHReferencedROISequence')
                                %structureNumber = dvhsequence.(['Item_',num2str(i)]).DVHReferencedROISequence.Item_1.ReferencedROINumber;
                                structureNumber = dvhReferencedROISequence.get(0).getInts(805699716); %org.dcm4che3.data.Tag.ReferencedROINumber); % .(['Item_',num2str(i)]).DVHReferencedROISequence.Item_1.ReferencedROINumber;
                                indROINumber = find(structureNumberV==structureNumber);
                                if ~isempty(indROINumber)
                                    %dataS(dvhsAdded+1).structureName = structureNameC{indROINumber};
                                    dataS(dvhsAdded).structureName = structureNameC{indROINumber};
                                end
                            end
                            %dataS(dvhsAdded+1).doseScale = dvhsequence.(['Item_',num2str(i)]).DVHDoseScaling;
                            dataS(dvhsAdded).doseScale = dvhObj.getDoubles(805568594); %org.dcm4che3.data.Tag.DVHDoseScaling;
                            %binWidthsV = dvhsequence.(['Item_',num2str(i)]).DVHData(1:2:end);
                            binWidthsV = dvhObj.getDoubles(805568600); %org.dcm4che3.data.Tag.DVHData;
                            binWidthsV = binWidthsV(1:2:end);
                            volumeBinsV = dvhObj.getDoubles(805568600); %org.dcm4che3.data.Tag.DVHData;
                            volumeBinsV = volumeBinsV(2:2:end);

                            %maxDVHDose = dvhsequence.(['Item_',num2str(i)]).DVHMaximumDose;
                            %minDVHDose = dvhsequence.(['Item_',num2str(i)]).DVHMinimumDose;
                            doseBinsV = [];
                            doseBinsV(1) = 0;
                            for iBin = 2:length(binWidthsV)
                                doseBinsV(iBin) = doseBinsV(iBin-1) + binWidthsV(iBin);
                            end
                            dataS(dvhsAdded).DVHMatrix(:,1) = doseBinsV(:);
                            if strcmpi(dataS(dvhsAdded).volumeType,'cumulative')
                                if length(doseBinsV) > 1
                                    %volumeBinsV = dvhsequence.(['Item_',num2str(i)]).DVHData(2:2:end);
                                    volumeBinsV = diff(volumeBinsV(1)-volumeBinsV);
                                    volumeBinsV = [volumeBinsV(1); volumeBinsV(:)];
                                else
                                    volumeBinsV = dataS(dvhsAdded).DVHMatrix(:,1)*NaN;
                                end
                            else
                                %volumeBinsV = dvhsequence.(['Item_',num2str(i)]).DVHData(2:2:end);
                                volumeBinsV = volumeBinsV(:);
                            end
                            dataS(dvhsAdded).DVHMatrix(:,2) = volumeBinsV;
                            %dvhsAdded = dvhsAdded + 1;
                        end

                    end
                end
            end

        end

        %     case 'digitalFilm'
        %         populate_planC_digitalFilm_field(fieldName, dcmdir);
        %     case 'RTTreatment'
        %         populate_planC_RTTreatment_field(fieldName, dcmdir);
        %     case 'IM'
        %         populate_planC_IM_field(fieldName, dcmdir);

    case 'beams'
        [seriesC, typeC] = extract_all_series(dcmdir_patient);
        %%%DK Re-writing RTPLAN file to not use MATLAB image processing
        %%%tool.

        seriesNum = find(strncmpi('RTPLAN', typeC, 6));

        planCount = 0;

        for j = 1:length(seriesNum)
            RTPLAN = seriesC{seriesNum(j)}.Data;

            for planNum = 1:length(RTPLAN)

                planCount = planCount + 1;

                planobj  = scanfile_mldcm(RTPLAN(planNum).file);

                %Populate each field in the dose structure.
                for i = 1:length(names)
                    dataS(planCount).(names{i}) = populate_planC_beams_field(names{i}, RTPLAN(planNum), planobj);
                end

            end

        end
        rtPlans= dataS;

    case 'beamGeometry'

        dataS = initializeCERR('beamGeometry');

        for i = 1:length(rtPlans)
            beamGeometryS = populate_planC_beamGeometry_field(rtPlans(i), dataS);

            for j = 1:length(beamGeometryS)
                dataS = dissimilarInsert(dataS, beamGeometryS(j), length(dataS)+1);
            end
        end

        %%%OLD RTPLAN import code. commented DK?
        % %         %Place RTPLAN into planC{indexS.beams}
        % %         plansAdded = 0;
        % %         for seriesNum = 1:length(seriesC)
        % %             RTPLAN = seriesC{seriesNum}.Data;
        % %             for planNum =  1:length(RTPLAN)
        % %
        % %                 if strcmpi(typeC{seriesNum}, 'RTPLAN')
        % %
        % %                     try
        % %                         if plansAdded == 0m
        % %                             dataS =
        % dicominfo(seriesC{seriesNum}.Data(planNum).file);
        % %                         else
        % %                             dataS(plansAdded + 1) = dicominfo(seriesC{seriesNum}.Data(planNum).file);
        % %                         end
        % %                         plansAdded = plansAdded + 1;
        % %                         warning('Matlab''s Image Processing Toolbox was used to read RTPLAN...')
        % %                     catch
        % %                         warning('Matlab''s Image Processing Toolbox is required to read RTPLAN. Ignoring...')
        % %                     end
        % %
        % %                 end
        % %
        % %             end
        % %         end

    case 'GSPS'
        [seriesC, typeC]    = extract_all_series(dcmdir_patient);
        supportedTypes      = {'PR'};
        gspsAdded          = 0;

        %hWaitbar = waitbar(0,'Loading GSPS. Please wait...');

        numGspsSeries = length(find(strcmpi(typeC, 'PR')==1));

        %Place each structure into its own array element.
        for seriesNum = 1:length(seriesC)

            %if ismember(typeC{seriesNum}, supportedTypes)
            if strcmpi(typeC{seriesNum}, 'PR')
                GSPS = seriesC{seriesNum}.Data;
                for k = 1:length(GSPS)
                    gspsobj  = scanfile_mldcm(GSPS(k).file);

                    %Graphic Annotation Sequence.
                    el = gspsobj.getValue(7340033); %org.dcm4che3.data.Tag.GraphicAnnotationSequence; %hex2dec('00700001')

                    % ROI = strobj.getInt(org.dcm4che2.data.Tag.ROIContourSequence);

                    if ~isempty(el)
                        nGsps = el.size();
                    else
                        nGsps = 0;
                    end
                    curGspsNum = 1;
                    for j = 1:nGsps

                        %Populate each field in the structure field set
                        for i = 1:length(names)
                            dataS(gspsAdded+1).(names{i}) = populate_planC_gsps_field(names{i}, GSPS, curGspsNum, gspsobj);
                        end
                        curGspsNum = curGspsNum + 1;
                        gspsAdded = gspsAdded + 1;

                        %waitbar(gspsAdded/(nGsps*length(GSPS)*numGspsSeries), hWaitbar, 'Loading Annotations, Please wait...');

                    end
                end
            end

        end
        %close(hWaitbar);
        pause(0.1);

    case 'registration'
        [seriesC, typeC]    = extract_all_series(dcmdir_patient);
        supportedTypes      = {'REG'};
        registrationsAdded  = 0;
        frameOfRefUIDC      = {};
        %Place each REG into its own array element.
        for seriesNum = 1:length(seriesC)

            %             if ismember(typeC{seriesNum}, supportedTypes)
            if strcmpi(typeC{seriesNum}, 'REG')

                REG = seriesC{seriesNum}.Data; %wy RTDOSE{1} for import more than one dose files;
                for regNum = 1:length(REG)
                    regobj  = scanfile_mldcm(REG(regNum).file);

                    % Frame of Reference UID
                    %frameOfRefUID = getTagValue(regobj, '00200052');
                    frameOfRefUID = char(regobj.getStrings(2097234)); %org.dcm4che3.data.Tag.FrameOfReferenceUID;

                    %Populate each field in the dose structure.
                    for i = 1:length(names)
                        dataS(registrationsAdded+1).(names{i}) = populate_planC_registration_field(names{i}, REG(regNum), regobj);
                    end

                    if isempty(frameOfRefUIDC)
                        frameOfRefUIDC{1} = frameOfRefUID;
                    else
                        frameOfRefUIDC{end+1} = frameOfRefUID;
                    end

                    % Check if frame of reference UID matches any
                    % existing doses and add slices.

                    % If frame of reference UID does not match, then
                    % create a new dose
                    registrationsAdded = registrationsAdded + 1;
                end
            end

        end

    case 'importLog'
        %Implementation is unnecessary.

    case 'CERROptions'
        pathStr = getCERRPath;
        optName = [pathStr 'CERROptions.json'];
        dataS = opts4Exe(optName);

    case 'indexS'
        %Implementation is unnecessary.

    otherwise
        % disp(['DICOM Import has no methods defined for import into the planC{indexS.' cellName '} structure, leaving empty.']);
end

end

function [dataS,dosesAdded] = addDoseToSeries(doseS,dataS,dosesAdded)
if length(doseS) > 1
    doseArray = [];
    zValues = [];
    for slcNum = 1:length(doseS)
        doseArray(:,:,slcNum) = doseS(slcNum).doseArray;
        zValues(slcNum) = doseS(slcNum).zValues;
    end
    % Sort doseArray as per zValues
    [zValues,indSort] = sort(zValues);
    doseArray(:,:,1:end) = doseArray(:,:,indSort);
    doseS(1).doseArray = doseArray;
    doseS(1).zValues = zValues;
end
dosesAdded = dosesAdded + 1;

if isempty(dataS)
    dataS = doseS(1);
else
    dataS(dosesAdded) = doseS(1);
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function generateDoseMeshes()
%function clearDoseMeshes()
%This function generatessurface meshes for the dose whose meshRep flag is
%set to 1. Note that the The contourLevels are selected from current options.
%
%APA, 06/21/07
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

global planC
indexS = planC{end};

%Set Matlab path to directory containing the library
currDir = cd;
meshDir = fileparts(which('libMeshContour.dll'));
cd(meshDir)

%Generate new dose meshes
contourLevels = getIsoDoseLevels;
waitbarH = waitbar(0,'Generating surface meshes for dose...');
for doseNum = 1:length(planC{indexS.dose})
    if isfield(planC{indexS.dose}(doseNum),'meshRep') && isnumeric(planC{indexS.dose}(doseNum).meshRep) && planC{indexS.dose}(doseNum).meshRep == 1
        doseUID = planC{indexS.dose}(doseNum).doseUID;
        doseVolume = planC{indexS.dose}(doseNum).doseArray;
        doseVolume = permute(doseVolume,[2 1 3]);
        [xVals, yVals, zVals] = getDoseXYZVals(planC{indexS.dose}(doseNum));
        calllib('libMeshContour','loadVolumeData',doseUID,xVals, yVals, zVals, double(doseVolume))
        for level = 1:length(contourLevels)
            calllib('libMeshContour','generateSurface', doseUID, [doseUID,'_',num2str(contourLevels(level))], double(contourLevels(level)), uint16(10));
        end
        calllib('libMeshContour','clear',doseUID)
    end
    waitbar(doseNum/length(planC{indexS.dose}),waitbarH)
end
close(waitbarH)

%switch back the current irectory
cd(currDir)


text:
function [t,EA,EAd] = fig_animate(tout,Xout,Uout,Xdout,Udout,Uext,p)

flag_movie = p.flag_movie;

if flag_movie

    try
        name = ['test.mp4'];
        vidfile = VideoWriter(name,'MPEG-4');
    catch ME
        name = ['test'];
        vidfile = VideoWriter(name,'Motion JPEG AVI');
    end
    open(vidfile);
end

%% smoothen for animation
t = (tout(1):p.simTimeStep:tout(end));
X = interp1(tout,Xout,t);
U = interp1(tout,Uout,t);
Xd = interp1(tout,Xdout,t);
Ud = interp1(tout,Udout,t);
Ue = interp1(tout,Uext,t);

%% loop through frames
figure('Position',[200 100 1000 600]);
set(0, 'DefaultFigureRenderer', 'opengl');
set(gcf, 'Color', 'white')
N = 3;M = 3;    % subplot size

nt = length(t);
EA = zeros(nt,3);
EAd = zeros(nt,3);
for ii = 1:nt
    EA(ii,:) = fcn_X2EA(X(ii,:));
    EAd(ii,:) = fcn_X2EA(Xd(ii,:));
end

% for ZOH force
t2 = repelem(t,2);
t2(1) = []; t2(end+1) = t2(end);
U2 = repelem(U,2,1);

%%%%%% subfigure handles %%%%%%%
h_x = subplot(N,M,3);
h_dx = subplot(N,M,6);
h_w = subplot(N,M,9);
h_u = subplot(N,M,[7,8]);

for ii = 1:p.playSpeed:nt
    %% The main animation
    % plot setting
    pcom = X(ii,1:3)';
    h_main = subplot(N,M,[1,2,4,5]);
    hold on; grid on;axis square;axis equal;
    h_main.XLim = [pcom(1)-0.5 pcom(1)+0.5];
    h_main.YLim = [pcom(2)-0.5 pcom(2)+0.5];
    h_main.ZLim = [-0.2 0.6];

    viewPt = [0.2,0.5,0.2];
    view(viewPt);

    % plot robot & GRF
    % real
    fig_plot_robot(X(ii,:)',U(ii,:)',Ue(ii,:)',p)
    % desired
    fig_plot_robot_d(Xd(ii,:)',0*Ud(ii,:)',p)

    % text
    txt_time = ['t = ',num2str(t(ii),2),'s'];
    text(pcom(1),pcom(2),0.4,txt_time)
    txt_vd = ['vd = ',num2str(Xd(ii,4),2),'m/s'];
    text(pcom(1),pcom(2),0.5,txt_vd)
    txt_v = ['v = ',num2str(X(ii,4),2),'m/s'];
    text(pcom(1),pcom(2),0.45,txt_v)

    %% states
    %%%%%%%%% position %%%%%%%%%
    plot(h_x,t(1:ii),X(1:ii,1),'r',...
         t(1:ii),X(1:ii,2),'g',...
         t(1:ii),X(1:ii,3),'b',...
         t(1:ii),Xd(1:ii,1),'r--',...
         t(1:ii),Xd(1:ii,2),'g--',...
         t(1:ii),Xd(1:ii,3),'b--','linewidth',1)
    h_x.XLim = [t(1) t(end)];
    set( get(h_x,'Title'), 'String', 'Position [m]');

    %%%%%%%%% velocity %%%%%%%%%
    plot(h_dx,t(1:ii),X(1:ii,4),'r',...
         t(1:ii),X(1:ii,5),'g',...
         t(1:ii),X(1:ii,6),'b',...
         t(1:ii),Xd(1:ii,4),'r--',...
         t(1:ii),Xd(1:ii,5),'g--',...
         t(1:ii),Xd(1:ii,6),'b--','linewidth',1)
    h_dx.XLim = [t(1) t(end)];
    set( get(h_dx,'Title'), 'String', 'Velocity [m/s]');

    %%%%%%%%% Angular velocity %%%%%%%%%
    plot(h_w,t(1:ii),X(1:ii,16),'r',...
         t(1:ii),X(1:ii,17),'g',...
         t(1:ii),X(1:ii,18),'b',...
         t(1:ii),Xd(1:ii,16),'r--',...
         t(1:ii),Xd(1:ii,17),'g--',...
         t(1:ii),Xd(1:ii,18),'b--','linewidth',1)
    h_w.XLim = [t(1) t(end)];
    set( get(h_w,'Title'), 'String', 'Angular velocity [rad/s]' );

    %% control
    plot(h_u,t2(1:2*ii),U2(1:2*ii,3),'r',...
         t2(1:2*ii),U2(1:2*ii,6),'g',...
         t2(1:2*ii),U2(1:2*ii,9),'b',...
         t2(1:2*ii),U2(1:2*ii,12),'k',...
         t(1:ii),Ud(1:ii,3),'r--',...
         t(1:ii),Ud(1:ii,6),'g--',...
         t(1:ii),Ud(1:ii,9),'b--',...
         t(1:ii),Ud(1:ii,12),'k--','linewidth',1)
    h_u.XLim = [t(1) t(end)];
    set( get(h_u,'Title'), 'String', 'Fz [N]' );

    %% make movie
    if flag_movie
        writeVideo(vidfile, getframe(gcf));
    end

    drawnow
    if ii < nt
        cla(h_main);
    end

end

if flag_movie
    close(vidfile);
end

author: YanranDing

repo: RF-MPC

save_path: github-repos/MATLAB/YanranDing-RF-MPC


text:
function [ Q, sumQ] = OffpolicyQlearning150816( qldata3 , gamma, alpha, numtraces)
% OFF POLICY Q LEARNING

%initialisation of variables
sumQ=zeros(numtraces,1);  %record sum of Q after each iteration
nact=numel(unique(qldata3(:,3)))-1;   %nr of actions
ncl=numel(unique(qldata3(:,2)));
Q=zeros (ncl, nact);
maxavgQ=1;
modu=100;
listi=find(qldata3(:,1)==1);   %position of 1st step of each episodes in dataset
nrepi=numel(listi);  %nr of episodes in the dataset
jj=1;

 for j=1:numtraces

    i=listi(floor(rand()*(nrepi-2))+1);  %pick one episode randomly (not the last one!)
    trace = [];

    while qldata3(i+1,1)~=1
    S1=qldata3(i+1,2);
    a1=qldata3(i+1,3);
    r1=qldata3(i+1,4);
     step = [ r1, S1, a1 ];
     trace = [trace ; step];
    i=i+1;
    end

    tracelength = length(trace(:,1));
    return_t = trace(tracelength,1); % get last reward as return for penultimate state and action.

    for t=tracelength-1:-1:1       %Step through time-steps in reverse order
        s = trace(t,2); % get state index from trace at time t
        a = trace(t,3); % get action index
        Q(s,a) = (1-alpha)*Q(s,a) + alpha*return_t; % update Q.
        return_t = return_t*gamma + trace(t,1); % return for time t-1 in terms of return and reward at t
    end

     sumQ(jj,1)=sum(sum(Q));
     jj=jj+1;

 if mod(j,500*modu)==0  %check if can stop iterating (when no more improvement is seen)
%      sumQ(jj,1)=sum(sum(Q));
%      jj=jj+1;
     s=mean(sumQ(j-49999:j));
     d=(s-maxavgQ)/maxavgQ;
     if abs(d)<0.001
         break   %exit routine
     end
     maxavgQ=s;
 end

 end

 sumQ(jj:end)=[];

end

author: matthieukomorowski


text:
function mel = frq2mel(frq)
%FRQ2ERB  Convert Hertz to Mel frequency scale MEL=(FRQ)
%	mel = frq2mel(frq) converts a vector of frequencies (in Hz)
%	to the corresponding values on the Mel scale which corresponds
%	to the perceived pitch of a tone

%	The relationship between mel and frq is given by:
%
%	m = ln(1 + f/700) * 1000 / ln(1+1000/700)
%
%  	This means that m(1000) = 1000
%

mel = sign(frq).*log(1+abs(frq)/700)*1127.01048;
if ~nargout
    plot(frq,mel,'-x');
    xlabel(['Frequency (' xticksi 'Hz)']);
    ylabel(['Frequency (' yticksi 'Mel)']);
end

author: bastamon


text:
%% Tilt and Twist Boundaries
%
%%
% If a material deforms through the movement of dislocations, rearrangement
% of dislocations to a low-energy configuration may happen during
% deformation (i.e. in slow, geologic deformation) or or afterwards (in
% many metals). In any case, the arrangement of dislocation walls can lead
% to so-called subgrains boundaries. If such a boundary is composed of edge
% dislocations, it is called a tilt boundary and the rotation axis relating
% both parts of the grain at each side can be expected to be within the
% boundary plane (ideally parallel to the edge dislocation line). If the
% boundary is composed of screw dislocations, the rotation axis should be
% normal to the boundary. Between those end-members, there are general
% boundaries where the rotation axis is not easily related to the type of
% dislocations unless further information is available.
%
% In this chapter we discuss the computation of the misorientation axes at
% subgrain boundaries and discuss whether they vote for twist or tilt
% boundaries. We start by importing an sample EBSD data set and computing
% all subgrain boundaries as it is described in more detail in the chapter
% <subGrainBoundaries.html Subgrain Boundaries>.

% load some test data
mtexdata forsterite silent

% remove one pixel grains
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'));
ebsd(grains(grains.grainSize<5)) = [];

% compute subgrain boundaries with 1.5 degree threshold angle
[grains,ebsd.grainId] = calcGrains(ebsd('indexed'),'threshold',[1*degree, 15*degree]);

% lets smooth the grain boundaries a bit
grains = smooth(grains,5);

% set up the ipf coloring
cKey = ipfColorKey(ebsd('fo').CS.properGroup);
cKey.inversePoleFigureDirection = xvector;
color = cKey.orientation2color(ebsd('fo').orientations);

% plot the forsterite phase
plot(ebsd('fo'),color,'faceAlpha',0.5,'figSize','large')

% init override mode
hold on

% plot grain boundares
plot(grains.boundary,'linewidth',2)

% compute transparency from misorientation angle
alpha = grains('fo').innerBoundary.misorientation.angle / (5*degree);

% plot the subgrain boundaries
plot(grains('fo').innerBoundary,'linewidth',1.5,'edgeAlpha',alpha,'linecolor','b');

% stop override mode
hold off

%%
% In the above plot we have marked all subgrain boundaries in blue and
% adjusted the transparency value according to the misorientation angle.
%
%% Misorientation Axes
%
% When analysing the misorientation axes of the subgrain boundary
% misorientations we need to distinguish whether we look at the
% misorientation axes in crystal coordinates or in specimen coordinates.
% Lets start with the misorientation axes in crystal coordinates which can
% directly be computed by the command <orientation.axis.html |axis|>.

% extract the Forsterite subgrain boundaries
subGB = grains('fo').innerBoundary;

% plot the misorientation axes in the fundamental sector
plot(subGB.misorientation.axis,'fundamentalRegion','figSize','small')

%%
% Obviously from the above plot it is not easy to judge about prefered
% misorientation axes. We get more insight if we <DensityEstimation.html
% compute the density distribution> of the misorientation axes and look for
% <S2FunOperations.html#4 local extrema>.

% compute the density distribution of misorientation axes
density = calcDensity(subGB.misorientation.axis,'halfwidth',3*degree);

% plot them
plot(density,'figSize','small')
mtexColorbar

% find the two prefered misorientation axes
[~,hkl] = max(density,'numLocal',2); round(hkl)

%%
% We find two preferred misorientation axes - (001) and (071). *TODO*: can
% this be interpreted?
%
%% The misorientation axis in specimen coordinates
%
% The computation of the misorientation axis in specimen coordinates is a
% little bit more complicated as it is impossible using only the
% misoriention. In fact we require the adjacent orientations on both sides
% of the subgrain boundaries. We can find those by making use of the
% |ebsdId| stored in the grain boundaries. The command

oriGB = ebsd('id',subGB.ebsdId).orientations

%%
% results in a $N \times 2$ matrix of orientations with rows corresponding
% to the boundary segments and two columns for both sides of the boundary.
% The misorientation axis in specimen coordinates is again computed by the
% command <orientation.axis.html |axis|>

axS = axis(oriGB(:,1),oriGB(:,2),'antipodal')

% plot the misorientation axes
plot(axS,'MarkerAlpha',0.2,'MarkerSize',2,'figSize','small')

%%
% We have used here the option |antipodal| as we have no fixed ordering of
% the grains at the two sides of the grain boundaries. For a more
% quantitative analysis we again compute the corresponding density
% distribution and find the preferred misorientation axes in specimen
% coordinates

density = calcDensity(axS,'halfwidth',5*degree);
plot(density,'figSize','small')
mtexColorbar

[~,pos] = max(density)
annotate(pos)

%% Tilt and Twist Boundaries
%
% Subgrain boundaries are often assumed to form during deformation by the
% accumulation of edge or screw dislocations. In the first extremal case of
% exclusive edge dislocations the misorientation axis is parallel to the
% deformation line and within the boundary plane. Such boundaries are
% called *tilt boundaries*. In the second extremal case of exclusive screw
% dislocations the misorientation axis is the screw axis and is parallel to
% the boundary normal. Such boundaries are called *twist boundaries*.
%
% In the case of 2d EBSD data one usually has not the full boundary
% information, but only the trace of the boundary with the measurement
% surface. Hence, it is impossible to distinguish tilt and twist
% boundaries. However, for twist boundaries misorientation axis must be normal
% to the boundary trace. This means, if the misorientation axis lays in the
% measurement plane and normal to the boundary trace, the boundary is quite
% likely to be a twist boundary. At the other hand, if the misorientation axis
% is parallel to the trace of a boundary, the boundary is quite likely to be a
% tilt boundary.
% We can be easily check the latter situation from our EBSD data, which allows us to
% exclude certain boundaries to be twist boundaries and to be most likely tilt
% boundaries To do so, we colorize in the following plot all subgrain boundaries
% according to the angle between the boundary trace and the misorientation axis.
% Blue subgrain boundaries are very likely tilt boundaries, while red subgrain
% boundaries are can be either tilt or twist boundaries.

plot(ebsd('fo'),color,'faceAlpha',0.5,'figSize','large')

% init override mode
hold on

% plot grain boundares
plot(grains.boundary,'linewidth',2)

% colorize the subgrain boundaries according the angle between boundary
% trace and misorientation axis
plot(subGB,angle(subGB.direction,axS)./degree,'linewidth',2)
mtexColorMap blue2red
mtexColorbar

hold off

author: mtex-toolbox

repo: mtex

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


text:
function DataSet = prtDataGenManual
%prtDataGenManual Manually specify 2-D data via mouse clicking
%
%   dataSet = prtDataGenManual outputs data X and labels Y defined
%   interactively by clicking on a figure.  Left clickes place points on
%   the figure, and right clicks change the labels for the points to be
%   placed.  To finish, close the figure, or click the small "Done" button
%   on the bottom left of the figure.
%
%   Example:
%
%   dataSet = prtDataGenManual;
%
%   See also: prtDataSetClass, prtDataGenBiModal, prtDataGenIris,
%   prtDataGenManual, prtDataGenMary, prtDataGenNoisySinc,
%   prtDataGenOldFaithful,prtDataGenProtate, prtDataGenSprial,
%   prtDataGenSpiral3Regress, prtDataGenUnimodal, prtDataGenSwissRoll,
%   prtDataGenUnimodal, prtDataGenXor

hF = figure;
X = zeros(0,2);
Y = zeros(0,1);
currentLabel = 0;
maxLabel = 0;

hA = axes;
set(hA,'buttondownfcn',@prtDataGenManualAxesButtonDownFcn);
hB = uicontrol('style','pushbutton','string','Done','callback',@(o,e) close(hF),'units','normalized','Position',[0.05 0.05 0.05 0.05]);
dataManualTitle;

DataSet = prtDataSetClass(X,Y,'name','Manualy Defined DataSet');

uiwait(hF);

    function dataManualTitle
        title(sprintf('Left-click to add an H_{%d} observation. Right-click to change hypothesis.',currentLabel));
    end

    function prtDataGenManualAxesButtonDownFcn(o,e)
        persistent previous;
        sType = get(hF,'selectionType');
        previous = sType;
        point = get(gca,'CurrentPoint');

        switch sType
            case 'normal'
                DataSet = catObservations(DataSet,prtDataSetClass(point(2,1:2),currentLabel));
                V = axis;
                hP = plot(DataSet);
                set(hP,'hittest','off');
                set(gca,'buttondownfcn',@prtDataGenManualAxesButtonDownFcn);
                axis(V);
                if currentLabel == maxLabel
                    maxLabel = maxLabel + 1;
                end
                dataManualTitle;
            case 'alt'
                currentLabel = currentLabel + 1;
                if currentLabel > maxLabel
                    currentLabel = 0;
                end
                dataManualTitle;
            %this happens when double-clicking *either* mouse button, so
            %it's impossible to say which mouse button was pushed.  Which
            %is not really cool.  In either casel clicking too fast results
            %in no changes to the GUI.
            % case 'open'
            %      disp('No double clicks allowed');
            %      % keyboard
        end
    end

end

author: covartech

repo: PRT

save_path: github-repos/MATLAB/covartech-PRT


text:
function [w] = mne_read_w_file(filename)
%
% [w] = mne_read_w_file(filename)
%
% Reads a binary w file into the structure w with the following fields
%
% vertices - vector of vertex indices (0-based)
% data     - vector of data values
%

%
%   Author : Matti Hamalainen, MGH Martinos Center
%   License : BSD 3-clause
%
%
%     $Id$
%
%     Revision 1.6  2006/04/23 15:29:41  msh
%     Added MGH to the copyright
%
%     Revision 1.5  2006/04/10 23:26:54  msh
%     Added fiff reading routines
%
%     Revision 1.4  2005/12/05 20:23:21  msh
%     Added fiff_save_evoked. Improved error handling.
%
%     Revision 1.3  2005/11/21 03:19:12  msh
%     Improved error handling
%
%     Revision 1.2  2005/11/21 02:15:51  msh
%     Added more routines
%
%     Revision 1.1  2005/11/21 01:41:57  msh
%     Introduced structures and start all function names with mne_
%
%
me='MNE:mne_read_w_file';
if(nargin ~= 1)
   error(me,'usage: [w] = mne_read_w_file(filename)');
end

% open it as a big-endian file
[fid,message] = fopen(filename, 'rb', 'b') ;
if (fid < 0)
   error(me,message);
end

fread(fid, 1, 'int16') ;
vnum       = mne_fread3(fid) ;
w.data     = zeros(vnum,1) ;
w.vertices = zeros(vnum,1) ;
for i=1:vnum
  w.vertices(i) = mne_fread3(fid) ; % vertex number (0-based)
  w.data(i)     = fread(fid, 1, 'float') ; % vertex value
end
w.vertices = uint32(w.vertices);

fclose(fid) ;
return;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% SUMMARY:  Calculate probability of gauss mix model
% AUTHOR:   QIUQIANG KONG, Queen Mary University of London
% Created:  19-09-2015
% Modified: 15-11-2015 Modify output size
%           20-11-2015 debug the order of [p,M]
% -----------------------------------------------------------
% input
%   X               input data; size: N*p; dim 1: num of data, dim 2: feature dim
%   pi              prior of mix; dim 1: mix num
%   mu              size: p*M; dim 1: feature dim, dim 2: mix num
%   Sigma           size: p*p*M; dim 1,2: feature dim, dim 3: mix num
% output
%   probs           probability of input data, size: N*1
% ===========================================================
function probs = Gmmpdf(X, prior, mu, Sigma)
N = size(X,1);          % num of data
[p,M] = size(mu);       % mix num & feature dim
probs = zeros(N,1);     % init output array
for m = 1:M
    probs = probs + prior(m) * mvnpdf(X, mu(:,m)', Sigma(:,:,m));
end
end

author: qiuqiangkong

repo: matlab-hmm


text:
function PlotPacketTimeStamps(aedat, newFigure, diffs)

%{
Takes 'aedat' - a data structure containing an imported .aedat file,
as created by ImportAedat.
When imported from aedat3, there is an info field packetTimeStamps. Plot
this.
%}

if ~isfield(aedat, 'info') || ~isfield(aedat.info, 'packetTimeStamps')
    disp('No packet timestamps found')
    return
end

% newFigure - default is to create a new figure in which to plot - supress
% this with the newFigure flag
if ~exist('newFigure', 'var') || newFigure
    figure
end
if exist('diffs', 'var') && diffs
    % In the diffs mode, plots the difference in timestamp from one packet
    % to the next. Exclude packets before last timestanmp reset
    aedat = FindFirstPacketAfterLastTimeStampReset(aedat);
    firstPacket = aedat.info.firstPacketAfterLastTimeStampReset;
    ts = double(aedat.info.packetTimeStamps(firstPacket : end)) / 1e6;
    diffs = [0; ts(2 : end) - ts(1 : end - 1)] * 1e3;
    plot(ts, diffs, '.-')
    xlabel('Packet timestamp (s)')
    ylabel('time difference from previous timestamp (ms)')
else
    plot(aedat.info.packetTimeStamps, '.-')
    xlabel('Packet number')
    ylabel('Time (us)')
end

author: panpanfei


text:
function [p,d]=berk2prob(b)
%BERK2PROB convert Berksons to probability
%
%  Inputs:  B(M,N)       matrix containing Berkson values
%
% Outputs:  P(M,N)       Corresponding probability values
%           D(M,N)       Corresponding derivatives dP/dB
%
% Berksons, or log-odds, are a nonlinear scale for measuring
% probability defined by B = log2(P./(1-P)).
% When Berksons are used to measure probability, a logistic
% psychometric function becomes linear.
%
% The inverse function is berk2prob()

%      Copyright (C) Mike Brookes 2014
%      Version: $Id: berk2prob.m 4501 2014-04-24 06:28:21Z 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
p=1-1./(1+pow2(b));
if nargout>1
    d=log(2)*p.*(1-p);
end

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD


text:
function title = p14_title ( )

%*****************************************************************************80
%
%% P14_TITLE returns the title for problem 14.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, string TITLE, the title of the problem.
%
  title = 'sin ( exp(-x) + exp(-4x) )';

  return
end

author: johannesgerer


text:
function [IDX,sep] = otsu(I,n)

%OTSU Global image thresholding/segmentation using Otsu's method.
%   IDX = OTSU(I,N) segments the image I into N classes by means of Otsu's
%   N-thresholding method. OTSU returns an array IDX containing the cluster
%   indices (from 1 to N) of each point. Zero values are assigned to
%   non-finite (NaN or Inf) pixels.
%
%   IDX = OTSU(I) uses two classes (N=2, default value).
%
%   [IDX,sep] = OTSU(...) also returns the value (sep) of the separability
%   criterion within the range [0 1]. Zero is obtained only with data
%   having less than N values, whereas one (optimal value) is obtained only
%   with N-valued arrays.
%
%   Notes:
%   -----
%   It should be noticed that the thresholds generally become less credible
%   as the number of classes (N) to be separated increases (see Otsu's
%   paper for more details).
%
%   If I is an RGB image, a Karhunen-Loeve transform is first performed on
%   the three R,G,B channels. The segmentation is then carried out on the
%   image component that contains most of the energy.
%
%   Example:
%   -------
%   load clown
%   subplot(221)
%   X = ind2rgb(X,map);
%   imshow(X)
%   title('Original','FontWeight','bold')
%   for n = 2:4
%     IDX = otsu(X,n);
%     subplot(2,2,n)
%     imagesc(IDX), axis image off
%     title(['n = ' int2str(n)],'FontWeight','bold')
%   end
%   colormap(gray)
%
%   Reference:
%   ---------
%   Otsu N, A Threshold Selection Method from Gray-Level Histograms,
%   IEEE Trans. Syst. Man Cybern. 9:62-66;1979
%
%   See also GRAYTHRESH, IM2BW
%
%   -- Damien Garcia -- 2007/08, revised 2010/03
%   Visit my website for more details about OTSU

error(nargchk(1,2,nargin))

% Check if is the input is an RGB image
isRGB = isrgb(I);

assert(isRGB | ndims(I)==2,...
    'The input must be a 2-D array or an RGB image.')

%% Checking n (number of classes)
if nargin==1
    n = 2;
elseif n==1;
    IDX = NaN(size(I));
    sep = 0;
    return
elseif n~=abs(round(n)) || n==0
    error('MATLAB:otsu:WrongNValue',...
        'n must be a strictly positive integer!')
elseif n>255
    n = 255;
    warning('MATLAB:otsu:TooHighN',...
        'n is too high. n value has been changed to 255.')
end

I = single(I);

%% Perform a KLT if isRGB, and keep the component of highest energy
if isRGB
    sizI = size(I);
    I = reshape(I,[],3);
    [V,D] = eig(cov(I));
    [tmp,c] = max(diag(D));
    I = reshape(I*V(:,c),sizI(1:2)); % component with the highest energy
end

%% Convert to 256 levels
I = I-min(I(:));
I = round(I/max(I(:))*255);

%% Probability distribution
unI = sort(unique(I));
nbins = min(length(unI),256);
if nbins==n
    IDX = ones(size(I));
    for i = 1:n, IDX(I==unI(i)) = i; end
    sep = 1;
    return
elseif nbins<n
    IDX = NaN(size(I));
    sep = 0;
    return
elseif nbins<256
    [histo,pixval] = hist(I(:),unI);
else
    [histo,pixval] = hist(I(:),256);
end
P = histo/sum(histo);
clear unI

%% Zeroth- and first-order cumulative moments
w = cumsum(P);
mu = cumsum((1:nbins).*P);

%% Maximal sigmaB^2 and Segmented image
if n==2
    sigma2B =...
        (mu(end)*w(2:end-1)-mu(2:end-1)).^2./w(2:end-1)./(1-w(2:end-1));
    [maxsig,k] = max(sigma2B);

    % segmented image
    IDX = ones(size(I));
    IDX(I>pixval(k+1)) = 2;

    % separability criterion
    sep = maxsig/sum(((1:nbins)-mu(end)).^2.*P);

elseif n==3
    w0 = w;
    w2 = fliplr(cumsum(fliplr(P)));
    [w0,w2] = ndgrid(w0,w2);

    mu0 = mu./w;
    mu2 = fliplr(cumsum(fliplr((1:nbins).*P))./cumsum(fliplr(P)));
    [mu0,mu2] = ndgrid(mu0,mu2);

    w1 = 1-w0-w2;
    w1(w1<=0) = NaN;

    sigma2B =...
        w0.*(mu0-mu(end)).^2 + w2.*(mu2-mu(end)).^2 +...
        (w0.*(mu0-mu(end)) + w2.*(mu2-mu(end))).^2./w1;
    sigma2B(isnan(sigma2B)) = 0; % zeroing if k1 >= k2

    [maxsig,k] = max(sigma2B(:));
    [k1,k2] = ind2sub([nbins nbins],k);

    % segmented image
    IDX = ones(size(I))*3;
    IDX(I<=pixval(k1)) = 1;
    IDX(I>pixval(k1) & I<=pixval(k2)) = 2;

    % separability criterion
    sep = maxsig/sum(((1:nbins)-mu(end)).^2.*P);

else
    k0 = linspace(0,1,n+1); k0 = k0(2:n);
    [k,y] = fminsearch(@sig_func,k0,optimset('TolX',1));
    k = round(k*(nbins-1)+1);

    % segmented image
    IDX = ones(size(I))*n;
    IDX(I<=pixval(k(1))) = 1;
    for i = 1:n-2
        IDX(I>pixval(k(i)) & I<=pixval(k(i+1))) = i+1;
    end

    % separability criterion
    sep = 1-y;

end

IDX(~isfinite(I)) = 0;

%% Function to be minimized if n>=4
    function y = sig_func(k)

        muT = sum((1:nbins).*P);
        sigma2T = sum(((1:nbins)-muT).^2.*P);

        k = round(k*(nbins-1)+1);
        k = sort(k);
        if any(k<1 | k>nbins), y = 1; return, end

        k = [0 k nbins];
        sigma2B = 0;
        for j = 1:n
            wj = sum(P(k(j)+1:k(j+1)));
            if wj==0, y = 1; return, end
            muj = sum((k(j)+1:k(j+1)).*P(k(j)+1:k(j+1)))/wj;
            sigma2B = sigma2B + wj*(muj-muT)^2;
        end
        y = 1-sigma2B/sigma2T; % within the range [0 1]

    end

end

function isRGB = isrgb(A)
% --- Do we have an RGB image?
% RGB images can be only uint8, uint16, single, or double
isRGB = ndims(A)==3 && (isfloat(A) || isa(A,'uint8') || isa(A,'uint16'));
% ---- Adapted from the obsolete function ISRGB ----
if isRGB && isfloat(A)
    % At first, just test a small chunk to get a possible quick negative
    mm = size(A,1);
    nn = size(A,2);
    chunk = A(1:min(mm,10),1:min(nn,10),:);
    isRGB = (min(chunk(:))>=0 && max(chunk(:))<=1);
    % If the chunk is an RGB image, test the whole image
    if isRGB, isRGB = (min(A(:))>=0 && max(A(:))<=1); end
end
end

author: Sable


text:
function dist_alias(H,G,nband)

% dist_alias        Calculate and Plot Distortion and Aliasing of Filter Banks
%
% Arguments:
% H                 Impulse responses of analysis filter
% G                 Impulse responses of synthesis filter
% nband             Number of subbands
%
% Reference:        N.J. Fliege, Multirate Digital Signal Processing, Wiley 1994
%
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

N = length(H);                  % Length of prototype filter

A = zeros(nband,2*N-1);
for m=0:nband-1
     for k=0:nband-1
        A(m+1,:)=A(m+1,:)+conv(H(k+1,:).*exp(-i*[0:N-1]*2*pi*m/nband),G(k+1,:));
     end
end

A0 = zeros(1,1024);
for m=2:nband
    A0=A0+(abs(fft(A(m,:),1024))/nband).^2;
end
A0 = sqrt(A0);

% Plot the distortion function

figure;
subplot(211);
% N1=256;
DFTpoint = 4096;
[h,w] = FreqResp(A(1,:),DFTpoint);
% [h,w]=freqz(A(1,:),1,N1*nband/2);
plot(w/pi,10*log10(abs(h)));
xlabel('Frequency,\omega (\pi)'); ylabel('Amplitude distortion (dB)');
title('Distortion plot');
xlim([0 1]);

% Plot the total aliasing distortion

subplot(212);
plot([0:1/512:1],A0(1:513));
xlabel('Frequency,\omega (\pi)'); ylabel('Aliasing error');
title('Total aliasing distortion as shown in (7.22) in Fliege book');

author: CharlesThaCat

repo: acoustic-interference-cancellation


text:
function [mm, order] = vargplvmReduceModel(model, P, dims)

% VARGPLVMREDUCEMODEL prunes out dimensions of the model.
% FORMAT
% DESC order the latent dimensions acorrding to the inputScales and
% reduces the model to have smaller number of latent dimensions.
% ARG model : the model to be reduced.
% ARG P : the number of dimensions to move to (setting to model.q will
% just reorder the dimensions in the model).
% ARG dims : (optional) explicit set of dimensions to use
% RETURN model : the model with the reduced number of dimensions.
%
% COPYRIGHT : Michalis K. Titsias, 2009
%
% MODIFICATIONS : Neil D. Lawrence, 2009, Patrick Sauer, 2011
%
% SEEALSO : vargplvmCreate

% VARGPLVM

if nargin == 3
    P = length(dims);
end

% create temporary model
options = vargplvmOptions('dtcvar');
options.kern =[];
options.numActive = model.k;
if ~strcmp(model.kern.type,'cmpnd')
  options.kern = model.kern.type;
elseif strcmp(model.kern.type,'cmpnd')
  options.kern{1} = model.kern.comp{1}.type;
  for i = 2:length(model.kern.comp)
    options.kern{i} = model.kern.comp{i}.type;
  end
end

if isstruct( model.betaTransform )
    options.betaTransform = model.betaTransform;
end

mm = vargplvmCreate(P, model.d, model.y, options);
N = size(model.vardist.means,1);

if ~strcmp(model.kern.type,'cmpnd')
  %
  if strcmp(model.kern.type,'rbfardjit') | strcmp(model.kern.type,'linard2') | strcmp(model.kern.type,'rbfard2')
    %
    if nargin == 2
        [vals, order] = sort(-model.kern.inputScales);
        mm.kern.inputScales = model.kern.inputScales(order(1:P));
    else
        order = [dims setdiff(1:length(model.kern.inputScales), dims)];
        mm.kern.inputScales = model.kern.inputScales(order);
    end
    %
  end
  %
else
  %
  for i = 1:length(model.kern.comp)
    %
    if strcmp(model.kern.comp{i}.type,'rbfardjit') | strcmp(model.kern.comp{i}.type,'linard2') | strcmp(model.kern.comp{i}.type,'rbfard2')
      %
      if nargin == 2
          [vals, order] = sort(-model.kern.comp{i}.inputScales);
          mm.kern.comp{i}.inputScales = model.kern.comp{i}.inputScales(order(1:P));
      else
          order = [dims setdiff(1:length(model.kern.comp{i}.inputScales), dims)];
          mm.kern.comp{i}.inputScales = model.kern.comp{i}.inputScales(order);
      end
      % you order only wrt the first ARD kernel you find
      break;
      %
    end
    %
  end
  %
end

mm.vardist.means  = model.vardist.means(:,order(1:P));
mm.vardist.covars = model.vardist.covars(:,order(1:P));

mm.X_u = model.X_u(:,order(1:P));
mm.X   = model.vardist.means(:,order(1:P));
mm.inputSclsOrder = order;

initParams   = vargplvmExtractParam(mm);
mm.numParams = length(initParams);

% This forces kernel computation.
mm = vargplvmExpandParam(mm, initParams);

author: SheffieldML

repo: vargplvm


text:
%function ind = MS_nearest(x,tau,v);
%
% returns the row vector containing the indicies of the nearest
% neighbours to each of the columns of x. Each point and its tau
% temporal neighbours are excluded from the search.
% v is an array (not necessarily logical) indicating which columns of x to
% use or, the relative importance of these columns, in the computation
% (i.e. use v(i)*x(i,:), not x(i,:)).
%
% default : tau=0
%           v=ones(1,length(x(:,1)))
%
% this is a mex version of nearneigh, and provides a speed-up of
% atleast 1000%.
%
% Michael Small
% michael.small@uwa.edu.au, http://school.maths.uwa.edu.au/~small/
% 15/7/04
% For further details, please see M. Small. Applied Nonlinear Time Series
% Analysis: Applications in Physics, Physiology and Finance. Nonlinear Science
% Series A, vol. 52. World Scientific, 2005. (ISBN 981-256-117-X) and the
% references therein.

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%     NeuroSLAM System Copyright (C) 2018-2019
%     NeuroSLAM: A Brain inspired SLAM System for 3D Environments
%
%     Fangwen Yu (www.yufangwen.com), Jianga Shang, Youjian Hu, Michael Milford(www.michaelmilford.com)
%
%     The NeuroSLAM V1.0 (MATLAB) was developed based on the OpenRatSLAM (David et al. 2013).
%     The RatSLAM V0.3 (MATLAB) developed by David Ball, Michael Milford and Gordon Wyeth in 2008.
%
%     Reference:
%     Ball, David, Scott Heath, Janet Wiles, Gordon Wyeth, Peter Corke, and Michael Milford.
%     "OpenRatSLAM: an open source brain-based SLAM system." Autonomous Robots 34, no. 3 (2013): 149-176.
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

% Jan 2, 2019
groundTruthFile = 'C:\NeuroSLAM_Datasets\02_NeuroSLAM_Groudtruth\02_SynPanData_GT.txt';
expMapFile = 'C:\NeuroSLAM_Datasets\03_NeuroSLAM_Experiments_Results\SynPanData\01_exp_map_ml.txt';
% plot_3d_multilayer_experience_map(groundTruthFile, expMapFile, xExpMapScaling, yExpMapScaling, zExpMapScaling, xExpMapTrans, yExpMapTrans, zExpMapTrans, xGtScaling, yGtScaling, zGtScaling)
plot_3d_multilayer_experience_map(groundTruthFile, expMapFile, 0.3, -0.35, 0.213, 0.5, 0, 0, 20,20,20);

author: cognav

repo: NeuroSLAM

save_path: github-repos/MATLAB/cognav-NeuroSLAM


%==============================================================================
% This code is part of the Matlab-based toolbox
% FAIR - Flexible Algorithms for Image Registration.
% For details see
% - https://github.com/C4IR and
% - http://www.siam.org/books/fa06/
%==============================================================================
%
% function varargout = regularizer(varargin)
%
% Main function for the regularization model, uses persistent parameter
%
% The regularization functional is either based on a linear differential
% operator B, i.e.
%
%        S(Y) = 0.5*alpha*hd*|B*Y|^2
%
%  or a general (non-linear) functional S(Y)
%
% Typical call
%    [Sc,dS,d2S] = regularizer(Y,omega,m)
%
%  1.     initialize (pick regularizer and set parameters)
%
%         regularizer('reset','regularizer','mfElastic','alpha',1e3,...
%               'mu',1,'lambda',0','rigid2D','center',center);
%
%  2.     [Sc,dS,d2S] = regularizer(Y,omega,m)
%        evaluates the regularization model, e.g. Sc = 0.5*yc'*B'*B*yc and its derivatives
%
%  A = alpha*hd*B'*B is made persistent for eficiency
%  hd = prod(omega./m)
%  an option grid is introduce to indicate the appropriate discretization:
%  staggered for elastic and diffusive,
%  cell-centered for curvature
%  nodal for TV, hyper elastic or diffusive EPI regularization
%==============================================================================

function varargout = regularizer(varargin)

persistent OPTN A

if nargin == 0 && nargout == 0 && isempty(OPTN),
    help(mfilename);
    return;
end;

% -----------------------------------------------------------------------------
% handle options
[method,OPTN,task,stop] = dealOptions(mfilename,OPTN,varargin{:});

% check the grid according to the regularizer to be used
% setup default solver for Gauss-Newton systems

if strcmp(task,'set') || strcmp(task,'reset'),

    switch method,
        case 'mfElastic',
            scheme     = 'elastic';
            matrixFree = 1;
            grid       = 'staggered';
            solver     = 'MG-elastic';

        case 'mbElastic',
            scheme     = 'elastic';
            matrixFree = 0;
            grid       = 'staggered';
            solver     = 'backslash';

        case 'mbElasticNodal',
            scheme     = 'elasticNodal';
            matrixFree = 0;
            grid       = 'nodal';
            solver     = 'backslash';

        case 'mfElasticNodal',
            scheme     = 'elasticNodal';
            matrixFree = 1;
            grid       = 'nodal';
            solver     = 'PCG';

        case 'mbCurvature',
            scheme     = 'curvature';
            matrixFree = 0;
            grid       = 'cell-centered';
            solver     = 'backslash';

        case 'mfCurvature',
            scheme     = 'curvature';
            matrixFree = 1;
            grid       = 'cell-centered';
            solver     = 'PCG-curvature';

        case 'mbHyperElastic',
            scheme     = 'hyperElastic';
            matrixFree = 0;
            grid       = 'nodal';
            solver     = 'backslash';

        case 'mfHyperElastic',
            scheme     = 'hyperElastic';
            matrixFree = 1;
            grid       = 'nodal';
            solver     = 'PCG-hyperElastic';

        case 'mbHyperElasticFEM',
            scheme     = 'hyperElasticFEM';
            matrixFree = 0;
            grid       = 'FEM';
            solver     = 'backslash';

        case 'mbElasticFEM',
            scheme     = 'elasticFEM';
            matrixFree = 0;
            grid       = 'FEM';
            solver     = 'backslash';

      case 'mfElasticFEM',
            scheme     = 'elasticFEM';
            matrixFree = 1;
            grid       = 'FEM';
            solver     = 'MG';

        otherwise
            scheme = method;
            [grid,matrixFree,solver] = feval(scheme,'para',[],[],OPTN{:});
    end;

    [dummy,OPTN] = dealOptions(mfilename,OPTN,'set','scheme',scheme,...
        'grid',grid,'matrixFree',matrixFree,'solver',solver);
end;

% return, if no further work has to be handled
if stop,
    varargout{1} = method;
    if nargout > 1, varargout{2} = OPTN;  varargout{3} = A;  end;
    return;
end
% -----------------------------------------------------------------------------
% do the work

% extract regularization parameters
scheme      = dealOptions(mfilename,OPTN,'get','scheme');
alpha       = dealOptions(mfilename,OPTN,'get','alpha');
matrixFree  = dealOptions(mfilename,OPTN,'get','matrixFree');

% extract variables
yc     = varargin{1};
omega  = varargin{2};
m      = varargin{3};
doDerivative = (nargout>1);
varargin = varargin(4:end);

for k=1:2:length(varargin), % overwrites defaults
    eval([varargin{k},'=varargin{',int2str(k+1),'};']);
end;
[Sc,dS,d2S] = feval(scheme,yc,omega,m,...
    'alpha',alpha,'matrixFree',matrixFree,...
    'doDerivative',doDerivative,OPTN{:},varargin{:});
varargout = {Sc,dS,d2S};

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

function runMinimalExample
fprintf('run minimal example for %s\n',mfilename)
keyboard
%==============================================================================


text:
function caxisChange(inputLimits,varargin)
	% caxis wrapper to also change colorbar range to be correct automatically.
	% Biafra Ahanonu
	% started: 2019.10.08 [19:15:59]
	% inputs
		% inputLimits - [lowerLimit upperLimit] vector
	% outputs
		%

	ciapkg.view.caxisChange(inputLimits,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
function channels = demCmu35VargplvmLoadChannels(Ytest, skel)
% DEMCMU35VARGPLVMLOADCHANNELS Given a dataset load the skeleton channels
%
% SEEALSO : demCmu35VargplvmAnimate, demCmu35gplvmVargplvm3.m
% VARGPLVM

% playSkel(Ytest, startInd)
    if nargin < 2 || isempty(skel)
        skel = acclaimReadSkel('35.asf');
        [tmpchan, skel] = acclaimLoadChannels('35_01.amc', skel);
    end

    %left indices
    xyzInd = [2];
    xyzDiffInd = [1 3];
    rotInd = [4 6];
    rotDiffInd = [5];
    generalInd = [7:38 41:47 49:50 53:59 61:62];
    startInd = 1;
    endInd = length(generalInd);
    channels(:, generalInd) = 180*Ytest(:, startInd:endInd)/pi;
    startInd = endInd + 1;
    endInd = endInd + length(xyzDiffInd);
    channels(:, xyzDiffInd) = cumsum(Ytest(:, startInd:endInd), 1);
    startInd = endInd + 1;
    endInd = endInd + length(xyzInd);
    channels(:, xyzInd) = Ytest(:, startInd:endInd);
    startInd = endInd + 1;
    endInd = endInd + length(xyzDiffInd);
    channels(:, xyzDiffInd) = cumsum(Ytest(:, startInd:endInd), 1);
    startInd = endInd + 1;
    endInd = endInd + length(rotInd);
    channels(:, rotInd) = asin(Ytest(:, startInd:endInd))*180/pi;
    channels(:, rotInd(end)) = channels(:, rotInd(end))+270;
    startInd = endInd + 1;
    endInd = endInd + length(rotDiffInd);
    channels(:, rotDiffInd) = 0;%cumsum(asin(Ytest(:, startInd:endInd)), 1))*180/pi;
    % skelPlayData(skel, channels, 1/25);


text:
function q = p04_q ( m, c, w )

%*****************************************************************************80
%
%% P04_Q evaluates the integral for problem p04.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Alan Genz,
%    A Package for Testing Multiple Integration Subroutines,
%    in Numerical Integration: Recent Developments, Software
%    and Applications,
%    edited by Patrick Keast and Graeme Fairweather,
%    Reidel, 1987, pages 337-340,
%    ISBN: 9027725144,
%    LC: QA299.3.N38.
%
%  Parameters:
%
%    Input, integer M, the dimension of the argument.
%
%    Input, real C(M), W(M), the problem parameters.
%
%    Output, real Q, the integral.
%
  q = 1.0;

  for i = 1 : m

    q = q * sqrt ( pi ) ...
      * ( erf ( c(i) * ( 1.0 - w(i) ) ) ...
        + erf ( c(i) *         w(i) ) ) ...
      / ( 2.0 * c(i) );

  end

  return
end

repo: jburkardt-m