text:
function varargout = Visualization_4d(varargin)
% VISUALIZATION_4D M-file for Visualization_4d.fig
%      VISUALIZATION_4D, by itself, creates a new VISUALIZATION_4D or raises the existing
%      singleton*.
%
%      H = VISUALIZATION_4D returns the handle to a new VISUALIZATION_4D or the handle to
%      the existing singleton*.
%
%      VISUALIZATION_4D('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in VISUALIZATION_4D.M with the given input arguments.
%
%      VISUALIZATION_4D('Property','Value',...) creates a new VISUALIZATION_4D or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before Visualization_4d_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to Visualization_4d_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Copyright 2002-2003 The MathWorks, Inc.

% Edit the above text to modify the response to help Visualization_4d

% Last Modified by GUIDE v2.5 07-Feb-2007 22:44:47

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @Visualization_4d_OpeningFcn, ...
                   'gui_OutputFcn',  @Visualization_4d_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

% --- Executes just before Visualization_4d is made visible.
function Visualization_4d_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to Visualization_4d (see VARARGIN)

% Choose default command line output for Visualization_4d
handles.output = hObject;

handles.folder=[pwd '/Visual_Data/'];
d = dir([handles.folder '*.mat']);
handles.datalist=char(d.name);

if isempty(handles.datalist)==1
    set(handles.select_data,'String','Visual_Data - Folder is Empty')
else
    set(handles.select_data,'String',handles.datalist)
end

handles.cbar=0;
guidata(hObject, handles);

% UIWAIT makes Visualization_4d wait for user response (see UIRESUME)
% uiwait(handles.figure1);

% --- Outputs from this function are returned to the command line.
function varargout = Visualization_4d_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;

% --- Executes on button press in rotate.
function rotate_Callback(hObject, eventdata, handles)
% hObject    handle to rotate (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of rotate

handles=guidata(gcbo);

if get(hObject,'Value')==get(hObject,'Max')
    rotate3d on
else
    rotate3d off
end

guidata(gcbo,handles)

% --- Executes on button press in boundary_slices.
function boundary_slices_Callback(hObject, eventdata, handles)
% hObject    handle to boundary_slices (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of boundary_slices
handles=guidata(gcbo);

if get(hObject,'Value')==get(hObject,'Max')
    set(handles.boundary,'Visible','on')
else
    set(handles.boundary,'Visible','off')
end

guidata(gcbo,handles)

% --- Executes on button press in grid_on_slice.
function grid_on_slice_Callback(hObject, eventdata, handles)
% hObject    handle to grid_on_slice (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of grid_on_slice

handles=guidata(gcbo);

if get(hObject,'Value')==get(hObject,'Max')
    set(handles.boundary,'EdgeColor','k')
    set(handles.sliceplane_x,'EdgeColor','k')
    set(handles.sliceplane_y,'EdgeColor','k')
    set(handles.sliceplane_z,'EdgeColor','k')

else
    set(handles.boundary,'EdgeColor','none')
    set(handles.sliceplane_x,'EdgeColor','none')
    set(handles.sliceplane_y,'EdgeColor','none')
    set(handles.sliceplane_z,'EdgeColor','none')
end

guidata(gcbo,handles)

% --- Executes on button press in grid_on.
function grid_on_Callback(hObject, eventdata, handles)
% hObject    handle to grid_on (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of grid_on

handles=guidata(gcbo);

if get(hObject,'Value')==get(hObject,'Max')
    grid on
else
    grid off
end

guidata(gcbo,handles)

% --- Executes on slider movement.
function x_slider_Callback(hObject, eventdata, handles)
% hObject    handle to x_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider

handles=guidata(gcbo);

i=get(hObject,'Value');

axes(handles.axes1)
delete(handles.sliceplane_x)
handles.sliceplane_x=slice(handles.X,handles.Y,handles.Z,handles.V,i,[],[]);

if get(handles.grid_on_slice,'Value')==0

    delete(handles.sliceplane_x_boundary)
    handles.sliceplane_x_boundary=plot3([i i i i i],[handles.ymin handles.ymax handles.ymax handles.ymin handles.ymin],...
        [handles.zmin handles.zmin handles.zmax handles.zmax handles.zmin],'k');
end

set(handles.sliceplane_x,'FaceColor','interp')

% set(handles.sliceplane_x,'Visible','off')
% set(handles.sliceplane_x(i),'Visible','on')

if get(handles.grid_on_slice,'Value')==get(handles.grid_on_slice,'Min')
    set(handles.sliceplane_x,'EdgeColor','none')
end
set(handles.sliceplane_x,'Visible','on')

set(handles.x_slice,'Value',1,'Enable','on')

guidata(gcbo,handles)

% --- Executes on slider movement.
function y_slider_Callback(hObject, eventdata, handles)
% hObject    handle to y_slider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles=guidata(gcbo);

i=get(hObject,'Value');

axes(handles.axes1)
delete(handles.sliceplane_y)
handles.sliceplane_y=slice(handles.X,handles.Y,handles.Z,handles.V,[],i,[]);

if get(handles.grid_on_slice,'Value')==0

    delete(handles.sliceplane_y_boundary)
    handles.sliceplane_y_boundary=plot3([handles.xmin handles.xmax handles.xmax handles.xmin handles.xmin ],[i i i i i],...
        [handles.zmin handles.zmin handles.zmax handles.zmax handles.zmin],'k');
end

set(handles.sliceplane_y,'FaceColor','interp')

if get(handles.grid_on_slice,'Value')==get(handles.grid_on_slice,'Min')
    set(handles.sliceplane_y,'EdgeColor','none')
end
set(handles.sliceplane_y,'Visible','on')
set(handles.y_slice,'Value',1,'Enable','on')
guidata(gcbo,handles)

% --- Executes on slider movement.
function z_slider_Callback(hObject, eventdata, handles)
% hObject    handle to wf (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles=guidata(gcbo);

i=get(hObject,'Value');
axes(handles.axes1)
delete(handles.sliceplane_z)
handles.sliceplane_z=slice(handles.X,handles.Y,handles.Z,handles.V,[],[],i);

if get(handles.grid_on_slice,'Value')==0
    delete(handles.sliceplane_z_boundary)
    handles.sliceplane_z_boundary=plot3([handles.xmin handles.xmax handles.xmax handles.xmin handles.xmin ],...
        [handles.ymin handles.ymin handles.ymax handles.ymax handles.ymin],[i i i i i],'k');
end

set(handles.sliceplane_z,'FaceColor','interp')

if get(handles.grid_on_slice,'Value')==get(handles.grid_on_slice,'Min')
    set(handles.sliceplane_z,'EdgeColor','none')
end

set(handles.sliceplane_z,'Visible','on')
set(handles.z_slice,'Value',1,'Enable','on')

guidata(gcbo,handles)

% --- Executes on button press in x_slice.
function x_slice_Callback(hObject, eventdata, handles)
% hObject    handle to x_slice (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of x_slice
handles=guidata(gcbo);

if get(hObject,'Value')==get(hObject,'Max')
    set(handles.sliceplane_x,'Visible','on')
    if get(handles.grid_on_slice,'Value')==0
        set(handles.sliceplane_x_boundary,'Visible','on')
    end
else
    set(handles.sliceplane_x,'Visible','off')
    set(handles.sliceplane_x_boundary,'Visible','off')
end

guidata(gcbo,handles)

% --- Executes on button press in y_slice.
function y_slice_Callback(hObject, eventdata, handles)
% hObject    handle to y_slice (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of y_slice
handles=guidata(gcbo);

if get(hObject,'Value')==get(hObject,'Max')
    set(handles.sliceplane_y,'Visible','on')
    if get(handles.grid_on_slice,'Value')==0
        set(handles.sliceplane_y_boundary,'Visible','on')
    end
else
    set(handles.sliceplane_y,'Visible','off')
    set(handles.sliceplane_y_boundary,'Visible','off')
end

guidata(gcbo,handles)

% --- Executes on button press in z_slice.
function z_slice_Callback(hObject, eventdata, handles)
% hObject    handle to z_slice (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of z_slice
handles=guidata(gcbo);

if get(hObject,'Value')==get(hObject,'Max')
    set(handles.sliceplane_z,'Visible','on')
    if get(handles.grid_on_slice,'Value')==0
        set(handles.sliceplane_z_boundary,'Visible','on')
    end

else
    set(handles.sliceplane_z,'Visible','off')
    set(handles.sliceplane_z_boundary,'Visible','off')
end

guidata(gcbo,handles)

% --- Executes on selection change in select_data.
function select_data_Callback(hObject, eventdata, handles)
% hObject    handle to select_data (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: contents = get(hObject,'String') returns select_data contents as cell array
%        contents{get(hObject,'Value')} returns selected item from select_data

handles=guidata(gcbo);

selected=get(hObject,'Value');
aaa=handles.datalist(selected,:);

eval(['load Visual_Data/' aaa])
cla

set([handles.x_slice handles.z_slice],'Value',0)
set(handles.y_slice,'Value',1)

set(handles.boundary_slices,'Value',0)
handles.X=X;
handles.Y=Y;
handles.Z=Z;
handles.V=V;

xmin=(min(min(min(X))));
ymin=(min(min(min(Y))));
zmin=(min(min(min(Z))));

xmax=(max(max(max(X))));
ymax=(max(max(max(Y))));
zmax=(max(max(max(Z))));

xmid=(xmax+xmin)/2;
ymid=(ymax+ymin)/2;
zmid=(zmax+zmin)/2;

handles.xmin=xmin;
handles.ymin=ymin;
handles.zmin=zmin;

handles.xmax=xmax;
handles.ymax=ymax;
handles.zmax=zmax;

alengthx=xmax-xmin;
alengthy=ymax-ymin;
alengthz=zmax-zmin;

xslice=xmin;
yslice=ymin;
zslice=zmin;

cen=[0.5*(xmax+xmin),0.5*(ymax+ymin),0.5*(zmax+zmin)];

[xx,yy,zz]=sphere(10);

xx=xx*(xmax-xmin)/2;
yy=yy*(ymax-ymin)/2;
zz=zz*(zmax-zmin)/2;

xx=xx+cen(1);
yy=yy+cen(2);
zz=zz+cen(3);

axes(handles.axes1)

outer_sphere=mesh(xx,yy,zz);
set(outer_sphere,'Visible','off')
hold on

set(handles.x_slider,'Min',xmin,'Max',xmax,'Value',xmid)%,'SliderStep',[0.01*alengthx 0.1*alengthx])
set(handles.y_slider,'Min',ymin,'Max',ymax,'Value',ymid)%,'SliderStep',[0.01*alengthy 0.1*alengthy])
set(handles.z_slider,'Min',zmin,'Max',zmax,'Value',zmid)%,'SliderStep',[0.01*alengthz 0.1*alengthz])

% qqq
% return

handles.boundary=slice(X,Y,Z,V,xslice,yslice,zslice);
set(handles.boundary,'FaceColor','interp')

if get(handles.boundary_slices,'Value')==0
    set(handles.boundary,'Visible','off')
else
    set(handles.boundary,'Visible','on')
end

hold on

if get(handles.color_bar,'Value')==1

    if handles.cbar==0
        handles.cbar=colorbar('Location','EastOutside');
    else
        delete(handles.cbar)
        handles.cbar=colorbar('Location','EastOutside');
    end
end

box on
axis tight
view(125,20)

% handles.sliceplane_x=slice(X,Y,Z,V,(xmax+xmin)/2,[],[]);
handles.sliceplane_x=slice(X,Y,Z,V,xmid,[],[]);
% handles.sliceplane_x_boundary=plot3([0 0 0 0 0],[ymin ymax ymax ymin ymin],[zmin zmin zmax zmax zmin],'k');
handles.sliceplane_x_boundary=plot3([xmid xmid xmid xmid xmid],[ymin ymax ymax ymin ymin],[zmin zmin zmax zmax zmin],'k');

% handles.sliceplane_y=slice(X,Y,Z,V,[],(ymax+ymin)/2,[]);
handles.sliceplane_y=slice(X,Y,Z,V,[],ymid,[]);
handles.sliceplane_y_boundary=plot3([xmin xmax xmax xmin xmin ],[ymid ymid ymid ymid ymid],[zmin zmin zmax zmax zmin],'k');

% handles.sliceplane_z=slice(X,Y,Z,V,[],[],(zmax+zmin)/2);
handles.sliceplane_z=slice(X,Y,Z,V,[],[],zmid);
handles.sliceplane_z_boundary=plot3([xmin xmax xmax xmin xmin ],[ymin ymin ymax ymax ymin],[zmid zmid zmid zmid zmid],'k');

set([handles.sliceplane_z handles.sliceplane_y handles.sliceplane_z],'FaceColor','interp')
% qqq

if get(handles.grid_on_slice,'Value')==1
    set([ handles.sliceplane_x_boundary handles.sliceplane_y_boundary handles.sliceplane_z_boundary ],'Visible','off')
else
    set(handles.boundary,'EdgeColor','none')
    set([ handles.sliceplane_x handles.sliceplane_y handles.sliceplane_z],'Edgecolor','none')
end

if get(handles.x_slice,'Value')==0
    set(handles.sliceplane_x_boundary ,'Visible','off')
    set(handles.sliceplane_x,'Visible','off')
end

if get(handles.y_slice,'Value')==0
    set(handles.sliceplane_y_boundary ,'Visible','off')
    set(handles.sliceplane_y,'Visible','off')
end

if get(handles.z_slice,'Value')==0
    set(handles.sliceplane_z_boundary ,'Visible','off')
    set(handles.sliceplane_z,'Visible','off')
end

camlookat(outer_sphere)

set([handles.boundary_slices handles.grid_on_slice handles.grid_on handles.x_slider handles.y_slider, ...
    handles.z_slider,handles.capture_image, handles.rotate  handles.x_slice handles.y_slice handles.z_slice, ...
    handles.color_bar, handles.x_label_txt, handles.y_label_txt,handles.z_label_txt,handles.Plot_title],'Enable','on')

set(handles.axes1,'PlotBoxAspectRatio',[1 1 1]);

ttt=title('Variation of \Omega with \alpha , \beta & \gamma');
set(ttt,'FontSize',12)

xlabel('\alpha')
ylabel('\beta')
zlabel('\gamma')

guidata(gcbo,handles)

function x_label_txt_Callback(hObject, eventdata, handles)
% hObject    handle to x_label_txt (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of x_label_txt as text
%        str2double(get(hObject,'String')) returns contents of x_label_txt as a double

handles=guidata(gcbo);

x_label=get(hObject,'String');
xl=xlabel(x_label);

guidata(gcbo,handles)

function z_label_txt_Callback(hObject, eventdata, handles)
% hObject    handle to z_label_txt (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of z_label_txt as text
%        str2double(get(hObject,'String')) returns contents of z_label_txt as a double

handles=guidata(gcbo);

z_label=get(hObject,'String');
% set(handles.z_text,'String',z_label)
zl=zlabel(z_label);

guidata(gcbo,handles)

function y_label_txt_Callback(hObject, eventdata, handles)
% hObject    handle to y_label_txt (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of y_label_txt as text
%        str2double(get(hObject,'String')) returns contents of y_label_txt as a double

handles=guidata(gcbo);

y_label=get(hObject,'String');
% set(handles.y_text,'String',y_label)
yl=ylabel(y_label);

guidata(gcbo,handles)

function Plot_title_Callback(hObject, eventdata, handles)
% hObject    handle to Plot_title (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of Plot_title as text
%        str2double(get(hObject,'String')) returns contents of Plot_title
%        as a double

handles=guidata(gcbo);

plot_title=get(hObject,'String');
% set(handles.Plot_title_static_txt,'String',plot_title)
ttt=title(plot_title)
set(ttt,'FontSize',12)

guidata(gcbo,handles)

% --- Executes on button press in color_bar.
function color_bar_Callback(hObject, eventdata, handles)
% hObject    handle to color_bar (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of color_bar

handles=guidata(gcbo);

if get(handles.color_bar,'Value')==1

    if handles.cbar==0
        handles.cbar=colorbar('Location','EastOutside');
    else
        delete(handles.cbar)
        handles.cbar=colorbar('Location','EastOutside');
    end
elseif handles.cbar~=0
    delete(handles.cbar)
    handles.cbar=0;
end

guidata(gcbo,handles)

% --- Executes on button press in capture_image.
function capture_image_Callback(hObject, eventdata, handles)
% hObject    handle to capture_image (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

handles=guidata(gcbo);
[file,path] = uiputfile('4D_Plot.jpg',['Capture as Image'])
filename=[path,file]

if file~=0
    set([ handles.select_data, handles.color_bar, handles.boundary_slices, handles.grid_on_slice,...
        handles.set_label_txt, handles.grid_on, handles.x_slice, handles.y_slice, handles.z_slice,...
        handles.rotate, handles.capture_image, handles.x_label_txt, handles.y_label_txt, handles.z_label_txt,...
        handles.Plot_title, handles.x_slider, handles.y_slider, handles.z_slider, handles.text1,...
        handles.set_title_txt, handles.text2, handles.text3,handles.uipanel1],'Visible','off')

    saveas(gca,filename)

    set([ handles.select_data, handles.color_bar, handles.boundary_slices, handles.grid_on_slice,...
        handles.set_label_txt,handles.grid_on, handles.x_slice, handles.y_slice, handles.z_slice,...
        handles.rotate, handles.capture_image, handles.x_label_txt, handles.y_label_txt, handles.z_label_txt,...
        handles.Plot_title, handles.x_slider, handles.y_slider, handles.z_slider, handles.text1,...
        handles.set_title_txt, handles.text2, handles.text3,handles.uipanel1],'Visible','on')
end

guidata(gcbo,handles)

function exit_Callback(hObject, eventdata, handles)
% hObject    handle to file_exit (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
handles = guidata(gcbo);
delete(handles.figure1);
% clc
a=0;
save Extra.dat a
delete('Extra.dat')
clear
closereq
load alltools
msg=['                                                                    ';...
     '     j.divahar@yahoo.com                                            '];

button = msgbox(msg,'Thank you For Trying !','custom',Tvanakam,Tcolormap);


text:
function [ y ] = tapas_huge_logit( x )
% Numerical stable calculation of logit function.
%
% INPUTS:
%   x - Array of double.
%
% OUTPUTS:
%   y - logit of x.
%

% Author: Yu Yao (yao@biomed.ee.ethz.ch)
% Copyright (C) 2019 Translational Neuromodeling Unit
%                    Institute for Biomedical Engineering,
%                    University of Zurich and ETH Zurich.
%
% This file is part of TAPAS, which is released under the terms of the GNU
% General Public Licence (GPL), version 3. For further details, see
% <https://www.gnu.org/licenses/>.
%
% This software is provided "as is", without warranty of any kind, express
% or implied, including, but not limited to the warranties of
% merchantability, fitness for a particular purpose and non-infringement.
%
% This software is intended for research only. Do not use for clinical
% purpose. Please note that this toolbox is under active development.
% Considerable changes may occur in future releases. For support please
% refer to:
% https://github.com/translationalneuromodeling/tapas/issues
%

y = log(x./(1-x));
y(y == Inf) = realmax;
y(y == -Inf) = -realmin;

end

author: translationalneuromodeling

repo: tapas


text:
% Read all images and extract point coordinates.
%
% All information needed are stored and retrieved
% from the function read_configuration,
% which in turn uses the --config=FILENAME command-line
% option.

% $Author: svoboda $
% $Revision: 2.0 $
% $Id: im2points.m,v 2.0 2003/06/19 12:07:11 svoboda Exp $
% $State: Exp $

clear all;

% add path to config data (TODO: obsolete?)
addpath ../Cfg
% add path for graphical output if needed
addpath ../OutputFunctions

SHOWFIG	  = 0; % show images during point extraction
STEP4STAT = 5; % step for computing average and std images, if 1 then all images taken

% Read configuration from whatever is specified on command-line (via --config=FILENAME)
config = read_configuration();

im.dir = config.paths.img;
im.ext = config.files.imgext;

NoCams = size(config.files.idxcams,2);	% number of cameras

% load image names
for i=1:NoCams,
  seq(i).camId = config.files.idxcams(i);
  if seq(i).camId > -1
	if strfind(expname,'oscar')
	  seq(i).data = dir([sprintf(im.dir,seq(i).camId),config.files.imnames,'*.',im.ext]);
	else
	  seq(i).data = dir([sprintf(im.dir,seq(i).camId),sprintf(config.files.imnames,seq(i).camId),im.ext]);
	end
  else
	seq(i).data = dir([im.dir,sprintf(config.files.imnames),im.ext]);
  end
  seq(i).size = size(seq(i).data,1);
  if seq(i).size<4
	error('Not enough images found. Wrong image path or name pattern?');
  end
end

% Expected number of 3D points is equal to the number of frames.
% In fact, some frames might be without any calibration point

% Becouse of non-consistent stopping of capturing, the sequences might
% have different number of images, select the minimal value as the right one
NoPoints = min([seq.size]);
pointsIdx = [1:STEP4STAT:NoPoints];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% beginning of the findings

t = cputime;
for i=1:NoCams,
  if ~exist(sprintf(config.files.avIM,seq(i).camId)),
	disp(sprintf('The average image of the camera %d is being computed',seq(i).camId));
	avIM = zeros(size(imread([sprintf(im.dir,seq(i).camId),seq(i).data(1).name])));
	for j=pointsIdx,
	  IM = imread([sprintf(im.dir,seq(i).camId),seq(i).data(j).name]);
	  avIM = avIM + double(IM);
	end
	avIM = uint8(round(avIM./size(pointsIdx,2)));
	imwrite(avIM,sprintf(config.files.avIM,seq(i).camId));
  else	disp('Average files already exist');
  end
end
disp(sprintf('Elapsed time for average images: %d [sec]',cputime-t))
% compute the standard deviations images that will be used for finding LEDs
% if not already computed
t = cputime;
for i=1:NoCams,
  if ~exist(sprintf(config.files.stdIM,seq(i).camId)),
	avIM = double(imread(sprintf(config.files.avIM,seq(i).camId)));
	disp(sprintf('The image of standard deviations of the camera %d is being computed',seq(i).camId));
	stdIM = zeros(size(imread([sprintf(im.dir,seq(i).camId),seq(i).data(1).name])));
	for j=pointsIdx,
	  IM = imread([sprintf(im.dir,seq(i).camId),seq(i).data(j).name]);
	  stdIM = stdIM + (double(IM)-avIM).^2;
	end
	stdIM = uint8(round(sqrt(stdIM./(size(pointsIdx,2)-1))));
	imwrite(stdIM,sprintf(config.files.stdIM,seq(i).camId));
  else
	disp('Image of standard deviations already exist')
  end
end

disp(sprintf('Elapsed time for computation of images [sec]: %d',cputime-t))

% find points in the images
Ws    = [];	  % joint image matrix
Res	  = [];	  % resolution of cameras
IdMat = ones(NoCams,NoPoints);
% IdMat is very important for Martinec&Pajdla filling [ECCV2002]
% it is a NoCams x NoPoints matrix,
% IdMat(i,j) = 0 -> no j-th point in i-th
% IdMat(i,j) = 1 -> point successfully detected

disp('*********************************************')
disp('Finding points (laser projections) in cameras')
disp(sprintf('Totally %d cameras, %d images for each cam', NoCams, NoPoints'))
disp('*********************************************')
for i=1:NoCams,
  t1 = cputime;
  disp(sprintf('Finding points in camera No: %0.2d',config.files.idxcams(i)))
  Points = [];
  avIM  = imread(sprintf(config.files.avIM,seq(i).camId));
  stdIM	= imread(sprintf(config.files.stdIM,seq(i).camId));
  for j=1:NoPoints,
	[pos,err] = getpoint([sprintf(im.dir,seq(i).camId),seq(i).data(j).name], SHOWFIG, config.imgs, avIM, stdIM);
	if err
	  IdMat(i,j) = 0;
	  Points = [Points, [NaN; NaN; NaN]];
	else
	  Points = [Points, [pos; 1]];
	end
  end
  Ws = [Ws; Points];
  Res= [Res; size(avIM,2), size(avIM,1)];
  t2 = cputime;
  disp(sprintf('Elapsed time for finding points in one camera: %d minutes %d seconds',floor((t2-t1)/60), round(mod((t2-t1),60))))
end

%%% End of the findings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if strfind(expname,'oscar')
  % needs special care for handling projector data
  ProjPoints = load(config.files.projdata,'-ASCII');
  Ws = [Ws; ProjPoints(:,2:3)'; ones(size(ProjPoints(:,1)'))];
  IdMat = [IdMat; ones(size(ProjPoints(:,1)'))];
  Res	= [Res; config.imgs.projres];
end

save(config.files.points, 'Ws','-ASCII')
save(config.files.Res, 'Res', '-ASCII')
save(config.files.IdMat, 'IdMat', '-ASCII')

% display the overall statistics
disp('Overall statistics from im2points:  ************************  ')
disp(sprintf('Total number of frames (possible 3D points): %d',NoPoints))
disp(sprintf('Total number of cameras %d', NoCams))
disp('More important statistics: *********************************  ')
disp(sprintf('Detected 3D points:                    %d', sum(sum(IdMat)>0)))
disp(sprintf('Detected 3D points in at least 3 cams: %d', sum(sum(IdMat)>2)))
disp(sprintf('Detected 3D points in ALL cameras:     %d', sum(sum(IdMat)==NoCams)))

author: strawlab

repo: MultiCamSelfCal


text:
function interp_Function_Coeffs(x1,x2)

mat1 = [1 0 0 0 0 0 0 0 0 0 0 0];
mat2 = [0 1 0 0 0 0 0 0 0 0 0 0];
mat3 = [0 0 1 0 0 0 0 0 0 0 0 0];
mat4 = [1 x1 x1^2 x1^3 -1 -x1 -x1^2 -x1^3 0 0 0 0];
mat5 = [0 1 2*x1 3*x1^2 0 -1 -2*x1 -3*x1^2 0 0 0 0];
mat6 = [0 0 2 6*x1 0 0 -2 -6*x1 0 0 0 0];

mat7 = [0 0 0 0 1 x2 x2^2 x2^3 -1 -x2 -x2^2 -x2^3];
mat8 = [0 0 0 0 0 1 2*x2 3*x2^2 0 -1 -2*x2 -3*x2^2];
mat9 = [0 0 0 0 0 0 2 6*x2 0 0 -2 -6*x2];

mat10= [0 0 0 0 0 0 0 0 1 1 1 1];
mat11= [0 0 0 0 0 0 0 0 0 1 2 3];
mat12= [0 0 0 0 0 0 0 0 0 0 2 6];

mat = [mat1; mat2; mat3; mat4; mat5; mat6; mat7; mat8; mat9; mat10; mat11; mat12];

rhs = [0 0 0 0 0 0 0 0 0 1 0 0]';

%mat1 = [-x1^2 0 x1^3 x1^2 x1 1];
%mat2 = [-2*x1 0 3*x1^2 2*x1 1 0];
%mat3 = [-2 0 6*x1 2 0 0];
%mat4 = [0 (x2-1)^2 x2^3 x2^2 x2 1];
%mat5 = [0 2*(x2-1) 3*x2^2 2*x2 1 0];
%mat6 = [0 2 6*x2 2 0 0];
%mat = [mat1; mat2; mat3; mat4; mat5; mat6];
%rhs = [0 0 0 1 0 0]';

coeffs = mat\rhs

a0 = coeffs(1);
a1 = coeffs(2);
a2 = coeffs(3);
a3 = coeffs(4);
b0 = coeffs(5);
b1 = coeffs(6);
b2 = coeffs(7);
b3 = coeffs(8);
c0 = coeffs(9);
c1 = coeffs(10);
c2 = coeffs(11);
c3 = coeffs(12);

ds = 0.01;
x = 0:ds:1;

figure(1)
for i=1:length(x)
    plot(x(i),g(coeffs,x(i),x1,x2),'*'); hold on;
end

figure(2)
for i=1:length(x)
    plot(x(i),gP(coeffs,x(i),x1,x2),'*'); hold on;
end

figure(3)
for i=1:length(x)
    plot(x(i),gPP(coeffs,x(i),x1,x2),'*'); hold on;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: evaluate interpolating polynomial
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = g(coeffs,x,x1,x2)

if x<=x1
    val = coeffs(1) + coeffs(2)*x + coeffs(3)*x^2 + coeffs(4)*x^3;
elseif x<=x2
    val = coeffs(5) + coeffs(6)*x + coeffs(7)*x^2 + coeffs(8)*x^3;
else
    val = coeffs(9) + coeffs(10)*x + coeffs(11)*x^2 + coeffs(12)*x^3;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: evaluate interpolating polynomial 1st derivative
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = gP(coeffs,x,x1,x2)

if x<=x1
    val = coeffs(2) + 2*coeffs(3)*x + 3*coeffs(4)*x^2;
elseif x<=x2
    val = coeffs(6) + 2*coeffs(7)*x + 3*coeffs(8)*x^2;
else
    val = coeffs(10) + 2*coeffs(11)*x + 3*coeffs(12)*x^2;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% FUNCTION: evaluate interpolating polynomial 2nd derivative
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function val = gPP(coeffs,x,x1,x2)

if x<=x1
    val = 2*coeffs(3) + 6*coeffs(4)*x;
elseif x<=x2
    val = 2*coeffs(7) + 6*coeffs(8)*x;
else
    val = 2*coeffs(11) + 6*coeffs(12)*x;
end

%
% figure(1)
% for i=1:length(x1_part)
%     x = x1_part(i);
%     plot(x,a*x^2,'ro'); hold on;
% end
%
% for i=1:length(x2_part)
%     x = x2_part(i);
%     plot(x,c*x^3+d*x^2+g*x+h,'*'); hold on;
% end
%
% for i=1:length(x3_part)
%    x = x3_part(i);
%    plot(x,-b*(x-1)^2+1,'go'); hold on;
% end
%
% xtotal = 0:ds:1;
% for i=1:length(xtotal)
%    x = xtotal(i);
%    plot(x,0.5*(1+tanh(4.3*(x-0.5))),'m*'); hold on;
% end
% title('function vals');
%
%
% figure(2)
% for i=1:length(x1_part)
%     x = x1_part(i);
%     plot(x,2*a*x,'ro'); hold on;
% end
%
% for i=1:length(x2_part)
%     x = x2_part(i);
%     plot(x,3*c*x^2+2*d*x+g,'*'); hold on;
% end
%
% for i=1:length(x3_part)
%    x = x3_part(i);
%    plot(x,-2*b*(x-1),'go'); hold on;
% end
% title('1st deriv');
%
%
% figure(3)
% for i=1:length(x1_part)
%     x = x1_part(i);
%     plot(x,2*a,'ro'); hold on;
% end
%
% for i=1:length(x2_part)
%     x = x2_part(i);
%     plot(x,6*c*x+2*d,'*'); hold on;
% end
%
% for i=1:length(x3_part)
%    x = x3_part(i);
%    plot(x,-2*b,'go'); hold on;
% end
% title('2nd deriv');

author: nickabattista


text:
function m = p01_m ( )

%*****************************************************************************80
%
%% P01_M returns the number of equations M for problem 1.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 April 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, integer M, the number of equations.
%
  m = 5;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [ ax,h ] = m_contfbar( varargin )
% M_CONTFBAR Draws a colour bar for contourf plots
%    M_CONTFBAR([X1,X2],Y,DATA,LEVELS) draws a horizontal colourbar
%    between the normalized coordinates (X1,Y) and (X2,Y) where
%    X/Y are both in the range 0 to 1 across the current axes.
%    Negative values can be used to position outside the axis.
%
%    Differences from COLORBAR are: the bar is divided into solid
%    colour patches exactly corresponding to levels provided by
%    CONTOURF(DATA,LEVELS) instead of showing the whole continuous
%    colourmap, the parent axis is not resized, the axis can be made
%    as large or small as desired, and the presence of values
%    above/below contoured levels is indicated by triangular pieces
%    (MATLAB 2014b or later).
%
%    M_CONTFBAR(X,[Y1,Y2],...) draws a vertical colourbar
%
%    The DATA,LEVELS pair can also be replaced with CS,CH where
%    [CS,CH]=CONTOURF(DATA,LEVELS).
%
%    M_CONTFBAR(...,'parameter','value') lets you specify extra
%    parameter/value pairs for the colourbar in the usual handle-graphics
%    way. Parameters you might set include 'xticks','xticklabels',
%    'xscale','xaxisloc' (or corresponding y-axis parameters for
%    a vertical colourbar) and 'fontsize'.
%
%    Additional parameter/value pairs allow special customization of the
%    colourbar:
%        'axfrac' : width of the colourbar (default 0.03)
%        'endpiece' : 'yes' (default) or 'no' show triangular
%                      endpieces.
%        'levels' : 'set' (default) shows a colourbar with exactly
%                   the levels in the LEVELS argument.
%                   'match' shows only the subsect of levels actually
%                   used in the CONTOURF call. E.g., if your data ranges
%                   from (say) 121 to 192, and LEVELS=[10:10:300],
%                   then only levels in 130:10:190 actually appear in both
%                   the CONTOURF and the colourbar.
%         'edgecolor' : 'none' removes edges between colors.
%
%    [AX,H]=M_CONTFBAR(...) returns the handle to the axis AX and to
%    the contourobject H. This is useful to add titles, xlabels, etc.
%
%    M_CONTFBAR(BAX,...) where BAX is an axis handle draws the colourbar
%    for that axis.
%
%    Note that if you assign a colormap to the specific axes BAX for which
%    you want a colorbar by calling COLORMAP(BAX,...), then in order to
%    import that colormap to the colorbar you must use M_CONTFBAR(BAX,...).
%    However, if you call COLORMAP without specifying an axes (which sets
%    the figure colormap) then you must call M_CONTFBAR without specifying
%    an axis.
%
%    Calling M_CONTFBAR a second time will replace the first colourbar.
%    Calling M_CONTFBAR without arguments will erase a colourbar.
%
%    See also COLORBAR

% R. Pawlowicz Nov/2017
%     Dec/2017 - inherit colormaps from the source axes.
%     Dec/2018 - if no axis specified, inherit colors from gca

% if users enters a '0' or a '1' as the first argument it can
% easily be interpreted as a figure handle - make sure this
% doesn't happen.
if nargin>0 && length(varargin{1})==1 && (varargin{1}==0 || varargin{1}==1 )
    varargin{1}=varargin{1}+eps;
end

% Is first argument an axis handle?
if nargin>0 && length(varargin{1})==1 && ishandle(varargin{1})
   if  strcmp(get(varargin{1},'type'),'axes')
      savax=varargin{1};
      varargin(1)=[];
    %  inheritcolormap=true;
   else
      error(['map: ' mfilename ':invalidAxesHandle'],...
            ' First argument must be an axes handle ');
   end
else
    savax=gca;
   % inheritcolormap=false;
end

% Delete any existing colorbars associated with savax
oldax=get(savax,'userdata');
if ~isempty(oldax) && ishandle(oldax) && strcmp('m_contfbar',get(oldax,'tag'))
    delete(oldax);
end

if isempty(varargin)  % No input arguments? - exit immediately after deleting old colorbar.
    return
elseif length(varargin)>=4
    posx=varargin{1};
    posy=varargin{2};
    Data=varargin{3};
    Levels=varargin{4};
    varargin(1:4)=[];
end

% Check trailing arguments for anything special

axfrac=.03;
endpiece=true;
sublevels=false;
edgeargs={};

k=1;
while k<=length(varargin)
    switch lower(varargin{k}(1:3))
        case 'axf'
            axfrac=varargin{k+1};
            varargin([k k+1])=[];
        case 'end'
            switch lower(varargin{k+1}(1))
                case 'n'
                    endpiece=false;
                otherwise
                    endpiece=true;
            end
            varargin([k k+1])=[];
        case 'lev'
            switch lower(varargin{k+1}(1))
                case 's'
                    sublevels=false;
                otherwise
                    sublevels=true;
            end
            varargin([k k+1])=[];
        case {'edg','lin'}
            edgeargs=varargin([k k+1]);
            varargin([k k+1])=[];
        otherwise
            k=k+2;
    end
end

if isempty(Levels)
    return;
end

if ishandle(Levels)  % CS,CH pair
    Data=get(Levels,'ZData');
    Levels=get(Levels,'LevelList');
end

% Min and max data values is all I need

ii=isfinite(Data(:));
minD=double(min(Data(ii)));
maxD=double(max(Data(ii)));

% I need to get the levels actually contoured
% if its just a set number of levels then we have
% to regenerate what they actually are
% Otherwise use the levels vector

if length(Levels)==1 || sublevels
    CS=contourc([1;1]*[minD maxD],Levels);
    k=1;
    nlevels=0;
    while k<size(CS,2)
        nlevels=nlevels+1;
        Clevel(nlevels)=CS(1,k);
        k=k+CS(2,k)+1;
    end
    Clevel=unique(Clevel);
else
    Clevel=Levels;
end
dC=Clevel(end)-Clevel(1);

% Form the colorbar, with or without the left/right triangles
% as needed depending on data range and levels chosen.

if minD<Clevel(1) && endpiece
    fakedata=[1;1]*[minD Clevel(1)];
    fakex=[1;1]*[Clevel(1)-.1*dC Clevel(1)];
    fakey=[1 2;1 0];
    leftpatch=true;
else
    fakedata=[1;1]*[Clevel(1) ];
    fakex=[1;1]*[Clevel(1)];
    fakey=[2;0];
    leftpatch=false;
end

if maxD>Clevel(end) && endpiece
    fakedata=[ fakedata [1;1]*[Clevel(end) maxD]];
    fakex=[fakex [1;1]*[Clevel(end) Clevel(end)+.1*dC] ];
    fakey=[fakey [2 1;0 1]];

else
    fakedata=[ fakedata [1;1]*[Clevel(end)]];
    fakex=[fakex [1;1]*[Clevel(end)] ];
    fakey=[fakey [2 ;0 ]];

end

axpos=get(savax,'position');
if  ( length(posx)==2 && length(posy)==1)
    horiz=true;
    cpos=[ axpos(1)+posx(1)*axpos(3) ...
         axpos(2)+(posy-1/2*axfrac)*axpos(4) ...
         diff(posx)*axpos(3) ...
         axfrac*axpos(4) ];

elseif  ( length(posx)==1 && length(posy)==2)
    horiz=false;
    cpos=[ axpos(1)+(posx-1/2*axfrac)*axpos(3) ...
         axpos(2)+posy(1)*axpos(4) ...
         axfrac*axpos(3) ...
         diff(posy)*axpos(4) ];
    tmp=fakex;
    fakex=fakey;
    fakey=tmp;
end

ax=axes('position',cpos);

if leftpatch   % If a left triangle is being drawn, colour a "behind"
    if horiz
        patch(Clevel(1)+[0 -.1*dC 0],[0 1 2],get(savax,'color'));
    else
        patch([0 1 2],Clevel(1)+[0 -.1*dC 0],get(savax,'color'));
    end
    hold on;
end

% Finally, draw the colourbar
[~,h]=contourf(fakex,fakey,fakedata,Levels,'clipping','off',edgeargs{:});
set(h,'clipping','off');   % Makes endpieces show in 2014b and later

line(fakex',fakey','color','k'); % long sides
line(fakex,fakey,'color','k');

if horiz
    set(ax,'xlim',Clevel([1 end]),'ytick',[],'clipping','off',...
     'ylim',[0 2],'ycolor','w');
    if posy>0.5, set(ax,'xaxislocation','top'); end
else
    set(ax,'ylim',Clevel([1 end]),'xtick',[],'clipping','off',...
     'xlim',[0 2],'xcolor','w');
    if posx>0.5, set(ax,'yaxislocation','right'); end
end

set(ax,'tickdir','out','box','off','layer','bottom',...
    'ticklength',[.03 .03],'tag','m_contfbar',varargin{:});

% Inherit the colormap. - fix Dec/28/2017
%if inheritcolormap
   drawnow;  % Update all properties - if this is missing
             % then we might be loading a default colormap
             % simply because the correct one is in  a'pending'
             % stack of graphics requests.  For colorbars there
             % appears to be 'peer' property that is perhaps
             % associated with handling this property.
   colormap(ax,colormap(savax));
%end

% Tuck away the info that there is a colorbar
set(savax,'userdata',ax);

% Return to drawing axis
set(gcf,'currentaxes',savax);

% Make the colorlimits match
if strcmp(get(gca,'CLimMode'),'manual')
    caxis(ax,get(gca,'clim'));
else
    caxis(get(ax,'clim'));
end

if nargout==0
  clear ax h
end

end

author: goGPS-Project


text:
function [out] = prep_envelope(dat,varargin)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prep_envelope
%
% Synopsis:
%   [out] = prep_envelope(dat,<var>)
%
% Example :
%    [out] = prep_envelope(dat)
%
% Arguments:
%     dat    - Epoched signal
% Options:
%     Time[ms] - time window. scalar or nx1 vector for weighting (default: 100)
%     Method - 'centered' or 'causal' (default: causal)
%
% Returns:
%     out - Envelope of the signal
%
% Description:
%     This function smoothly outlines the extremes of an oscillating
%     signal, continuous or epoched.
%     continuous data should be [time * channels]
%     epoched data should be [time * channels * trials]
%
% See also 'https://github.com/PatternRecognition/OpenBMI'
%
% Min-ho Lee, 01-2018
% mh_lee@korea.ac.kr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isfield(dat,'x')
    warning('OpenBMI: Data structure must have a field named ''x''')
    return
end

s = size(dat.x);
dat.x= reshape(abs(hilbert(dat.x(:,:))),s);
out= prep_movingAverage(dat,varargin{:});

author: PatternRecognition

repo: OpenBMI


text:
function[]=makefigs_curvemoments
%MAKEFIGS_CURVEMOMENTS  Makes a sample figure for CURVEMOMENTS.

load qgsnapshot

[cv,zeta,N,S,P]=psi2fields(qgsnapshot.psi);
P=frac(P,std(P(:)));

[xc,yc]=closedcurves(qgsnapshot.x,qgsnapshot.y,P,-2);
[xo,yo,kappa,R,L,a,b,theta]=curvemoments(xc,yc);

figure,jpcolor(qgsnapshot.x,qgsnapshot.y,P),axis equal, axis tight,
hold on,colormap gray,flipmap,cellplot(xc,yc,'2b'),
[k,l]=ab2kl(a,b);ellipseplot(k,l,theta,xo+sqrt(-1)*yo,'r')
title('Curves of Okubo-Weiss (blue), moment ellipses (red), and centroids (black)')
caxis([-5 5]),xtick([-5:1:5]*1000),ytick([-5:1:5]*1000),plot(xo,yo,'k.')
noxlabels,noylabels,fontsize 11
h=colorbar;
h.Label.String='Okubo-Weiss Parameter normalized by its own standard deviation';

%To print
if 0
    currentdir=pwd;
    cd([whichdir('jlab_license') '/figures'])
    print -dpng curvemoments
    crop curvemoments.png
    cd(currentdir)
end

author: jonathanlilly

repo: jLab

save_path: github-repos/MATLAB/jonathanlilly-jLab


text:
function varargout = easom(X)
% Easom function
%
%   EASOM([x1, x2]) returns the value of the Easom function at the
%   specified points. [x1] and [x2] may be vectors. The search
%   domain is
%
%               -100 < x_i < 100
%
%   The global minimum is
%
%               f(x1, x2) = f(pi, pi) = -1.

% Author: Rody P.S. Oldenhuis
% Delft University of Technology
% E-mail: oldenhuis@dds.nl
% Last edited 20/Jul/2009

    % if no input is given, return dimensions, bounds and minimum
    if (nargin == 0)
        varargout{1} = 2;  % # dims
        varargout{2} = [-100, -100]; % LB
        varargout{3} = [+100, +100]; % UB
        varargout{4} = [pi,pi]; % solution
        varargout{5} = -1; % function value at solution (too complicated)

    % otherwise, output function value
    else

        % keep values in the serach interval
        X(X < -100) = inf;     X(X > 100) = inf;

        % split input vector X into x1, x2
        if size(X, 1) == 2
            x1 = X(1, :);        x2 = X(2, :);
        else
            x1 = X(:, 1);        x2 = X(:, 2);
        end

        % output function value
        varargout{1} = -cos(x1).*cos(x2).*exp(-((x1-pi).^2 + (x2-pi).^2));

    end

end

author: Sable

repo: mcbench-benchmarks


text:
%% MATLAB code for blur-downsampled (BDx4) degradataion
scale = 4;
sigma = 1.6;

% generate LR image
kernelsize = ceil(sigma * 3) * 2 + 2;
kernel = fspecial('gaussian', kernelsize, sigma);
lq = imfilter(gt, kernel, 'replicate');
lq = lq(scale/2:scale:end-scale/2, scale/2:scale:end-scale/2, :);

author: ckkelvinchan

repo: BasicVSR-IconVSR


text:
function bool = QIBPresent
%"QIBPresent"
%   Returns true only if the matrix files necessary to run QIB are present
%   and in the proper location for this copy of CERR.  Else returns false.
%
%   JRA 10/15/04
%
%Usage:
%   function bool = QIBPresent
%
% 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/>.

bool = 0;
IMRTPdir = fileparts(which('IMRTP'));
if exist(fullfile(IMRTPdir, 'QIBData', 'aahn_6b.dat'), 'file') == 2 && ...
   exist(fullfile(IMRTPdir, 'QIBData', 'aahn_18b.dat'), 'file') == 2 && ...
   exist(fullfile(IMRTPdir, 'QIBData', 'QIB_lin_0pt125.mat'), 'file') == 2 && ...
   exist(fullfile(IMRTPdir, 'QIBData', 'QIBGauss_lin_0pt125.mat'), 'file') == 2
   bool = 1;
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function f = processmulti1D(fun,varargin)

% function f = processmulti1D(fun,varargin)
%
% <fun> is a function that expects a row vector as the first argument
% <varargin> is a set of inputs.  the first input should be one or multiple row vectors
%   concatenated along the first dimension.  the remaining inputs are the
%   same as expected by <fun>.
%
% apply <fun> to each vector and concatenate the results together (along the first dimension).
%
% example:
% isequal(processmulti1D(@mean,[1 2 3; 4 5 6]),[2; 5])

  prev = warning('query');  % make sure stupid warning about obsolete blkproc doesn't come on
  warning('off');
f = [];
for p=1:size(varargin{1},1)
  if p==1
    f = feval(fun,varargin{1}(p,:),varargin{2:end});  % this ensures class is respected
  else
    f = cat(1,f,feval(fun,varargin{1}(p,:),varargin{2:end}));
  end
end
  warning(prev);

author: cvnlab

repo: GLMsingle

save_path: github-repos/MATLAB/cvnlab-GLMsingle


text:
function P3 = reslice_seg_anatomy()

disp(['This script reslices all segmented n*het1*img anatomicals into space of functionals.'])
disp(['Output images are thresholded (1 or 0) to be mask images.'])

mdir = spm_get(-1,'*','Enter Main Directory for Study',pwd);
mypwd = pwd;
eval(['cd ' mdir])

!ls
dwcard = input('Enter wildcard for all individual subjects:','s');

P = get_filename(dwcard,['anatomy' filesep 'n*het1*img']);
disp(['reslicing these images:']), P

d = dir(dwcard);
P2 = get_filename([d(1).name '*'],'scan1/sn*ra*0001.img');
Pt = P2(1,:);

disp(['Reslicing to space of: ' Pt])

reslice_imgs(Pt,P,1);

P3 = get_filename(dwcard,['anatomy' filesep 'rn*het1*img']);

disp('Output in:')
P3

return

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
% runDataExport.m
%   Exports simulation data in JSON format
%
% Author: Ruby Lee
% Author: Jonathan Karr, jkarr@stanford.edu
% Affiliation: Covert Lab, Department of Bioengineering, Stanford University
% Last Updated: 8/22/2012
function runDataExport(simBatch, simIdx)
setWarnings();
setPath();
setPreferences();
warning('off', 'MATLAB:MKDIR:DirectoryExists');

%import classes
import edu.stanford.covert.cell.sim.util.DiskLogger;
import edu.stanford.covert.cell.sim.util.SimulationDiskUtil;
import edu.stanford.covert.cell.sim.util.SimulationEnsemble;
import edu.stanford.covert.util.SparseMat;

%options
verbosity = 1;
exportProps = {
    'Chromosome', 'abasicSites'
    'Chromosome', 'complexBoundSites'
    'Chromosome', 'damagedBases'
    'Chromosome', 'damagedSugarPhosphates'
    'Chromosome', 'doubleStrandedRegions'
    'Chromosome', 'gapSites'
    'Chromosome', 'hollidayJunctions'
    'Chromosome', 'intrastrandCrossLinks'
    'Chromosome', 'linkingNumbers'
    'Chromosome', 'monomerBoundSites'
    'Chromosome', 'ploidy'
    'Chromosome', 'polymerizedRegions'
    'Chromosome', 'segregated'
    'Chromosome', 'singleStrandedRegions'
    'Chromosome', 'strandBreaks'
    'Chromosome', 'superhelicalDensity'
    'FtsZRing', 'numEdges'
    'FtsZRing', 'numEdgesOneStraight'
    'FtsZRing', 'numEdgesTwoBent'
    'FtsZRing', 'numEdgesTwoStraight'
    'FtsZRing', 'numResidualBent'
    'Geometry', 'chamberVolume'
    'Geometry', 'cylindricalLength'
    'Geometry', 'pinchedDiameter'
    'Geometry', 'surfaceArea'
    'Geometry', 'totalLength'
    'Geometry', 'volume'
    'Geometry', 'width'
    'Host', 'isBacteriumAdherent'
    'Host', 'isInflammatoryResponseActivated'
    'Host', 'isNFkBActivated'
    'Host', 'isTLRActivated'
    'Mass', 'cell'
    'Mass', 'cellDry'
    'Mass', 'dnaWt'
    'Mass', 'media'
    'Mass', 'metaboliteWt'
    'Mass', 'proteinWt'
    'Mass', 'rnaWt'
    'Mass', 'total'
    'Mass', 'waterWt'
    'MetabolicReaction', 'doublingTime'
    'MetabolicReaction', 'fluxs'
    'MetabolicReaction', 'growth'
    'Metabolite', 'counts'
    'ProteinComplex', 'counts'
    'ProteinMonomer', 'counts'
    'Ribosome', 'boundMRNAs'
    'Ribosome', 'mRNAPositions'
    'Ribosome', 'nActive'
    'Ribosome', 'nNotExist'
    'Ribosome', 'nStalled'
    'Ribosome', 'stateOccupancies'
    'Ribosome', 'states'
    'Ribosome', 'tmRNAPositions'
    'Rna', 'counts'
    'RNAPolymerase', 'nActive'
    'RNAPolymerase', 'nFree'
    'RNAPolymerase', 'nNonSpecificallyBound'
    'RNAPolymerase', 'nSpecificallyBound'
    'RNAPolymerase', 'positionStrands'
    'RNAPolymerase', 'stateOccupancies'
    'RNAPolymerase', 'states'
    'RNAPolymerase', 'supercoilingBindingProbFoldChange'
    'RNAPolymerase', 'transcriptionFactorBindingProbFoldChange'
    'Stimulus', 'values'
    'Transcript', 'abortedTranscripts'
    'Transcript', 'boundTranscriptChromosome'
    'Transcript', 'boundTranscriptionUnits'
    'Transcript', 'boundTranscriptProgress'
    'Transcript', 'rnaBoundRNAPolymerases'
    };

%sanitize arguments
if ischar(simIdx)
    simIdx = str2double(simIdx);
end

%load simulation object which contains meta data
[~, ~, sim] = SimulationDiskUtil.getSimulation([simBatch filesep num2str(simIdx)]);

%export data to JSON
errorLog = [];
outDir = [SimulationDiskUtil.getBaseDir() filesep simBatch filesep num2str(simIdx) filesep 'json'];

%create output directory
mkdir(outDir);

%summary data
if verbosity > 0, fprintf('Exporting summary ...'), end
summary = load(sprintf('%s%s%s%s%d%ssummary.mat', SimulationDiskUtil.getBaseDir(), filesep, simBatch, filesep, simIdx, filesep));
mkdir([outDir filesep 'Summary']);
fields = fieldnames(summary);
for i = 1:numel(fields)
    data = summary.(fields{i});
    for j = 1:size(data, 1)
        tmp = struct('label', [], 'data', []);
        tmp.label = [fields{i} '_' num2str(j)];
        tmp.data = compressData(data(j, :), 1, 1000);

        writeJSON(tmp, [outDir filesep 'Summary' filesep fields{i} '_' num2str(j) '.json']);
    end

    clear data;
end
clear summary;
if verbosity > 0, fprintf(' done\n'), end

%complete data
for i = 1:numel(sim.states)
    stateID = sim.states{i}.wholeCellModelID(7:end);
    propIDs = [sim.states{i}.stateNames; sim.states{i}.dependentStateNames];
    if verbosity > 0, fprintf('Exporting %s ...\n', stateID), end
    mkdir([outDir filesep stateID]);
    for j = 1:numel(propIDs)
        if ~any(strcmp(exportProps(:, 1), stateID) & strcmp(exportProps(:, 2), propIDs{j}))
            continue;
        end

        if verbosity > 0, fprintf('\t%s ...', propIDs{j}), end
        %try
        %retrieve data
        result = SimulationEnsemble.load(simBatch, {stateID propIDs{j}}, ...
            [], [], 1, 'extract', simIdx);
        result = result.(stateID).(propIDs{j});

        %output data in JSON format
        if isa(sim.states{i}, 'edu.stanford.covert.cell.sim.state.Chromosome') && isa(result, 'SparseMat')
            %chromosome properties
            [subs, vals] = find(result);
            tfs1 = subs(:, 3) < size(result, 3);
            tfs2 = subs(:, 3) > 1;

            result1 = SparseMat(subs(tfs1, :), vals(tfs1, :), size(result) - [0 0 1]);
            result2 = SparseMat([subs(tfs2, 1:2) subs(tfs2, 3)-1], vals(tfs2, :), size(result) - [0 0 1]);

            times = downsampleData([1; 1 + find(full(permute(any(any(result1 ~= result2, 1), 2), [3 1 2]))); size(result, 3)], 1000);
            if numel(times) > 1000
                times = downsampleData((1:size(result, 3))', 1000);
            end

            tmp = struct('time', [], 'strand', [], 'pos', [], 'val', []);
            tmp = tmp(ones(numel(times), 1), 1);
            for k = 1:numel(times)
                [tmpsubs, tmpvals] = find(result(:, :, times(k)));

                tmp(k).time = times(k);
                tmp(k).pos = tmpsubs(:, 1);
                tmp(k).strand = tmpsubs(:, 2);
                tmp(k).val = tmpvals;
            end

            writeJSON(tmp, [outDir filesep stateID filesep propIDs{j} '.json']);
        else
            %generic output

            %sum over compartments
            if ...
                    (isa(sim.states{i}, 'edu.stanford.covert.cell.sim.MoleculeCountState') && strcmp(propIDs{j}, 'counts')) || ...
                    (isa(sim.states{i}, 'edu.stanford.covert.cell.sim.state.Stimulus') && strcmp(propIDs{j}, 'values'))
                result = sum(result, 2);
            end

            %output
            if size(result, 2) == 1
                %save rows separately
                for k = 1:size(result, 1)
                    tmp = struct('label', [], 'data', []);
                    tmp.label = [stateID '_' propIDs{j} '_' num2str(k)];
                    tmp.data = compressData(result(k, :, :), 0, 1000);

                    writeJSON(tmp, [outDir filesep stateID filesep propIDs{j} '_' num2str(k) '.json']);
                end
            else
                %save rows and columns separately
                for k = 1:size(result, 1)
                    for l = 1:size(result, 2)
                        tmp = struct('label', [], 'data', []);
                        tmp.label = [stateID '_' propIDs{j} '_' num2str(k) '_' num2str(l)];
                        tmp.data = compressData(result(k, l, :), 0, 1000);

                        writeJSON(tmp, [outDir filesep stateID filesep propIDs{j} '_' num2str(k) '_' num2str(l) '.json']);
                    end
                end
            end
        end

        clear result tmp;
        %catch exception
        %    errorLog = [errorLog sprintf('Error in exporting %s in simulation %s/%s\n%s', ...
        %        propIDs{j}, simBatch, simIdx, exception.getReport())]; %#ok<AGROW>
        %end
        if verbosity > 0, fprintf(' done\n'), end
    end
    if verbosity > 0, fprintf('done\n'), end
end

%print error log
if ~isempty(errorLog)
    throw(MException('runAnalysis:error', errorLog));
end

function result = compressData(data, timeOffset, maxSamples)
data = full(double(data));
data = data(:);

%get time points where value changes
idxs = unique([1; find(diff(data(1:end-1))); 1 + find(diff(data(1:end-1))); numel(data)]);

if numel(idxs) > maxSamples
    %downsample
    idxs = downsampleData((1:numel(data))', maxSamples);
end

%return as time, value pairs
result = [idxs-timeOffset  data(idxs)];

function idxs = downsampleData(idxs, maxSamples)
if numel(idxs) > maxSamples
    tmp = idxs(1:floor(numel(idxs)/maxSamples):end);
    if tmp(end) ~= idxs(end)
        tmp = [tmp; idxs(end)];
    end
    idxs = tmp;
end

function writeJSON(data, fileName)
fid = fopen(fileName, 'w+');
writeJSON_data(fid, data);
fclose(fid);

function writeJSON_data(fid, data)
if ischar(data)
    writeJSON_char(fid, data)
elseif isnumeric(data)
    writeJSON_numeric(fid, data)
elseif isstruct(data)
    writeJSON_struct(fid, data)
else
    throw(MException('writeJSON_data:error', 'Unsupported data type %s', class(data)))
end

function writeJSON_struct(fid, data)
if numel(data) == 1
    fprintf(fid, '{');
    fields = fieldnames(data);
    for i = 1:numel(fields)
        fprintf(fid, '"%s": ', fields{i});
        writeJSON_data(fid, data.(fields{i}));
        if i < numel(fields)
            fprintf(fid, ', ');
        end
    end
    fprintf(fid, '}');
else
    fprintf(fid, '[');
    for i = 1:numel(data)
        writeJSON_struct(fid, data(i));
        if i < numel(data)
            fprintf(fid, ',');
        end
    end
    fprintf(fid, ']');
end

function writeJSON_char(fid, data)
fprintf(fid, '"%s"', data);

function writeJSON_numeric(fid, data)
fprintf(fid, '[');

if ndims(data) > 2
    throw(MException('writeJSON:error', 'Data with more than 2 dimensions is not supported'))
end

if isempty(data)
elseif size(data, 2) == 1
    if all(abs(data - round(data)) < 1e-10)
        fprintf(fid, '%d', round(data(1)));
        if numel(data) > 1
            fprintf(fid, ',%d', round(data(2:end)));
        end
    else
        fprintf(fid, '%e', data(1));
        if numel(data) > 1
            fprintf(fid, ',%e', data(2:end));
        end
    end
elseif size(data, 2) == 2 && all(abs(data(:, 1) - round(data(:, 1))) < 1e-10)
    for j = 1:size(data, 1)
        fprintf(fid, '[');
        fprintf(fid, '%d', round(data(j, 1)));
        fprintf(fid, ',%e', data(j, 2:end));
        fprintf(fid, ']');
        if j < size(data, 1)
            fprintf(fid, ',');
        end
    end
else
    for j = 1:size(data, 1)
        fprintf(fid, '[');
        fprintf(fid, '%e', data(j, 1));
        fprintf(fid, ',%e', data(j, 2:end));
        fprintf(fid, ']');
        if j < size(data, 1)
            fprintf(fid, ',');
        end
    end
end

fprintf(fid, ']');

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [CV]=faceToVertexMeasure(F,V,CF)
% function [CV]=faceToVertexMeasure(F,V,CF)
% ------------------------------------------------------------------------
%
%
% Change log:
% 2009
% 2019/04/22 Minor changes to improve efficiency through improved allocation
% 2019/04/22 Renamed variables
% 2019/04/23 Updated to handle cell input
% ------------------------------------------------------------------------

%%

%Check if F is a cell of face arrays
cellMode=isa(F,'cell');

if cellMode==1
    CV=repmat({zeros(size(V,1),size(CF{1},2))},1,numel(CF));
    for qc=1:1:numel(F)
        % Get vertex-face connectivity matrix
        conStruct=patchConnectivity(F{qc},V,'vf');
        IND_F=conStruct.vertex.face;
        logicValid=IND_F>0;
        CV_now=ones(size(V,1),size(IND_F,2),size(CF{qc},2));
        cv=nan(size(IND_F));
        for q=1:1:size(CF{qc},2) %Loop over data columns
            cf=CF{qc}(:,q);
            cv(logicValid)=cf(IND_F(logicValid));
            CV_now(:,:,q)=cv;
        end
        CV{qc}=CV_now;
    end
    CV=mean(cell2mat(CV),2,'omitnan'); %Take mean across faces
    if size(CF{1},2)>1
        CV=permute(CV,[1 3 2]); %Get rid of singleton dimension
    end
else
    % Get vertex-face connectivity matrix
    conStruct=patchConnectivity(F,V,'vf');
    IND_F=conStruct.vertex.face;
    logicValid=IND_F>0;
    CV=ones(size(V,1),size(CF,2));
    cv=nan(size(IND_F));
    for q=1:1:size(CF,2) %Loop over data columns
        cf=CF(:,q);
        cv(logicValid)=cf(IND_F(logicValid));
        CV(:,q)=mean(cv,2,'omitnan');
    end
end

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function F = CalHV(points,bounds,k,nSample)
% Calculate the hypervolume-based fitness value of each solution

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

% This function is modified from the code in
% http://www.tik.ee.ethz.ch/sop/download/supplementary/hype/

    [N,M] = size(points);
    if M > 2
        % Use the estimated method for three or more objectives
        alpha = zeros(1,N);
        for i = 1 : k
            alpha(i) = prod((k-[1:i-1])./(N-[1:i-1]))./i;
        end
        Fmin = min(points,[],1);
        S    = unifrnd(repmat(Fmin,nSample,1),repmat(bounds,nSample,1));
        PdS  = false(N,nSample);
        dS   = zeros(1,nSample);
        for i = 1 : N
            x        = sum(repmat(points(i,:),nSample,1)-S<=0,2) == M;
            PdS(i,x) = true;
            dS(x)    = dS(x) + 1;
        end
        F = zeros(1,N);
        for i = 1 : N
            F(i) = sum(alpha(dS(PdS(i,:))));
        end
        F = F.*prod(bounds-Fmin)/nSample;
    else
        % Use the accurate method for two objectives
        pvec  = 1:size(points,1);
        alpha = zeros(1,k);
        for i = 1 : k
            j = 1 : i-1;
            alpha(i) = prod((k-j)./(N-j))./i;
        end
        F = hypesub(N,points,M,bounds,pvec,alpha,k);
    end
end

function h = hypesub(l,A,M,bounds,pvec,alpha,k)
% The recursive function for the accurate method

    h     = zeros(1,l);
    [S,i] = sortrows(A,M);
    pvec  = pvec(i);
    for i = 1 : size(S,1)
        if i < size(S,1)
            extrusion = S(i+1,M) - S(i,M);
        else
            extrusion = bounds(M) - S(i,M);
        end
        if M == 1
            if i > k
                break;
            end
            if alpha >= 0
                h(pvec(1:i)) = h(pvec(1:i)) + extrusion*alpha(i);
            end
        elseif extrusion > 0
            h = h + extrusion*hypesub(l,S(1:i,:),M-1,bounds,pvec(1:i),alpha,k);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function [ a, b ] = p49_lim ( )

%*****************************************************************************80
%
%% P49_LIM returns the integration limits for problem 49.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    04 November 2009
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Output, real A, B, the limits of integration.
%
  a = 0.0;
  b = 3.0;

  return
end


text:
%BILATERALTEXTUREFILTER  Applies the bilateral texture filter to an image
%
%     dst = cv.bilateralTextureFilter(src)
%     dst = cv.bilateralTextureFilter(src, 'OptionName', optionValue, ...)
%
% ## Input
% * __src__ Source image whose depth is 8-bit `uint8` or 32-bit `single`
%
% ## Output
% * __dst__ Destination image of the same size and type as `src`.
%
% ## Options
% * __FR__ Radius of kernel to be used for filtering. It should be positive
%   integer. default 3
% * __NumIter__ Number of iterations of algorithm, It should be positive
%   integer. default 1
% * __SigmaAlpha__ Controls the sharpness of the weight transition from edges
%   to smooth/texture regions, where a bigger value means sharper transition.
%   When the value is negative, it is automatically calculated. default -1
% * __SigmaAvg__ Range blur parameter for texture blurring. Larger value makes
%   result to be more blurred. When the value is negative, it is automatically
%   calculated as described in the paper. default -1
%
% It performs structure-preserving texture filter. For more details about this
% filter see [Cho2014].
%
% ## References
% [Cho2014]:
% > Hojin Cho, Hyunjoon Lee, Henry Kang, and Seungyong Lee.
% > "Bilateral texture filtering". ACM Transactions on Graphics,
% > 33(4):128:1-128:8, July 2014.
%
% See also: cv.rollingGuidanceFilter, cv.bilateralFilter
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
function hb_io_test08 ( )

%*****************************************************************************80
%
%% HB_IO_TEST08 tests HB_RHS_WRITE, HB_GUESS_WRITE, HB_EXACT_WRITE;
%
  neltvl = 0;
  nnzero = 126;
  nrhs = 1;
  nrhsix = 0;
  nrow = 32;

  exact = [ ...
    1.0,   2.0,   3.0,   4.0,   5.0,   6.0,   7.0,   8.0,   9.0,  10.0, ...
   11.0,  12.0,  13.0,  14.0,  15.0,  16.0,  17.0,  18.0,  19.0,  20.0, ...
   21.0,  22.0,  23.0,  24.0,  25.0,  26.0,  27.0,  28.0,  29.0,  30.0, ...
   31.0,  32.0 ]';
  guess = [ ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0,   1.0, ...
    1.0,   1.0 ]';
  indfmt = '(16I5)';
  mxtype = 'RUA';
  output_file = 'rua_32_rhs.txt';
  ptrfmt = '(16I5)';
  rhscrd = 12;
  rhsfmt = '(10F7.1)';
  rhsind = 0;
  rhsptr = 0;
  rhsval = [ ...
    101.0, 102.0, 103.0, 104.0, 107.0, 126.0, 201.0, 202.0, 209.0, 221.0, ...
    228.0, 302.0, 303.0, 306.0, 308.0, 309.0, 329.0, 403.0, 404.0, 405.0, ...
    412.0, 503.0, 505.0, 523.0, 527.0, 601.0, 606.0, 616.0, 703.0, 707.0, ...
    714.0, 721.0 ]';
  rhsvec = 0;
  rhstyp = 'FGX';

  fprintf ( 1, '\n' );
  fprintf ( 1, 'HB_IO_TEST08\n' );
  fprintf ( 1, '  HB_RHS_WRITE writes right hand sides to an HB file.\n' );
  fprintf ( 1, '  HB_GUESS_WRITE writes starting guesses to an HB file.\n' );
  fprintf ( 1, '  HB_EXACT_WRITE writes exact solutions to an HB file.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Writing the file "%s".\n', output_file );

  output_unit = fopen (output_file, 'wt' );

  if ( output_unit < 0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'HB_IO_TEST08 - Fatal error!\n' );
    fprintf ( 1, '  Error opening the file.\n' );
    return;
  end
%
%  Write the right hand sides.
%
  hb_rhs_write ( output_unit, nrow, nnzero, nrhs, nrhsix, ...
    rhscrd, ptrfmt, indfmt, rhsfmt, mxtype, rhstyp, rhsval, ...
    rhsind, rhsptr, rhsvec );
%
%  Write the right hand sides.
%
  hb_guess_write ( output_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp, ...
    guess );
%
%  Write the right hand sides.
%
  hb_exact_write ( output_unit, nrow, nrhs, rhscrd, rhsfmt, rhstyp, ...
    exact );

  fclose ( output_unit );

  return
end


text:
%{
load('dataset/trafficdb/traffic_patches.mat');
[M,m,n,p] = convert_video3d_to_2d(im2double(imgdb{100}));
out = run_algorithm('MC', 'OP-RPCA', M, [])
show_results(M.*out.Omega,out.L,out.S,out.O,p,m,n);
%}

lambda = 0.35; % [0,1] 0.5(+lowrank) 0.4 0.3 0.2(+sparse)
[L,S] = mr_pca_part(M,Omega,lambda);

author: andrewssobral

repo: lrslibrary


text:
function vw = computeResStdMap(vw, scanList, forceSave)
%
% vw = computeResStdMap(vw, [scanList], [forceSave])
%
% Cycles through tSeries, computing the residual std of the
% functional images.  Puts them together into a parameter map and
% calls setParameterMap to set vw.map = stdMap.
% Residual time series have the harmonics of freq removed
%
% scanList:
%   0 - do all scans
%   number or list of numbers - do only those scans
%   default - prompt user via selectScans dialog
%
% forceSave: 1 = true (overwrite without dialog)
%            0 = false (query before overwriting)
%           -1 = do not save
%
%
% If you change this function make parallel changes in:
%    computeCorAnal, computeStdMap, computeMeanMap
%
% rmk, 05/05/99
% djh 2/2001, mrLoadRet-3.0

if notDefined('forceSave'),   forceSave = 0;   end

nScans = viewGet(vw,'numScans');

if strcmp(vw.mapName,'resStdMap')
    % If exists, initialize to existing map
    map=vw.map;
else
    % Otherwise, initialize empty cell array
    map = cell(1,nScans);
end

% (Re-)set scanList
if ~exist('scanList','var')
    scanList = selectScans(vw);
elseif scanList == 0
    scanList = 1:nScans;
end
if isempty(scanList)
  error('Analysis aborted');
end

% Compute it
waitHandle = mrvWaitbar(0,'Computing res std images from the tSeries.  Please wait...');
ncScans = length(scanList);
for iScan = 1:ncScans
    scan    = scanList(iScan);
    dims    = viewGet(vw, 'sliceDims', scanNum);
    nCycles = viewGet(vw, 'numcycles', scanNum);
    datasz  = viewGet(vw, 'dataSize',  scanNum);

    map{scan} = NaN*ones(datasz);
    for slice = sliceList(vw,scan)
        resStd = computeTSResStd(vw,scan,slice,nCycles);
        map{scan}(:,:,slice) = reshape(resStd,dims);
    end
    mrvWaitbar(scan/ncScans)
end
close(waitHandle);

% Set parameter map
vw = setParameterMap(vw,map,'resStdMap');

% Save file
if forceSave >= 0, saveParameterMap(vw, [], forceSave); end

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
%% The Piezoelectricity Tensor
%
%%
% In this chapter we discuss how to compute and visualize piezoellectirc
% properties. At first, let us import some piezoelectric contents for a
% quartz specimen.

CS = crystalSymmetry('32', [4.916 4.916 5.4054], 'X||a*', 'Z||c', 'mineral', 'Quartz');

fname = fullfile(mtexDataPath,'tensor', 'Single_RH_quartz_poly.P');

P = tensor.load(fname,CS,'propertyname','piecoelectricity','unit','C/N','DoubleConvention')

%% Plotting the magnitude surface
% The default plot of the magnitude, which indicates, in which direction we
% have the most polarization. By default, we restrict ourselves to the
% unique region implied by crystal symmetry

% set some colormap well suited for tensor visualisation
setMTEXpref('defaultColorMap',blue2redColorMap);

plot(P)
mtexColorbar

%%
% but also, we can plot the whole crystal behavior

close all
plot(P,'complete','smooth','upper')
mtexColorbar

%%
% Most often, the polarization is illustrated as surface magnitude

close all
surf(P.directionalMagnitude)

%%
% Note, that for directions of negative polarization the surface is mapped
% onto the axis of positive, which then let the surface appear as a double
% coverage

%%
% Quite a famous example in various standard literature is a section through
% the surface because it can easily be described as an analytical
% solution. We just specify the plane normal vector

plotSection(P.directionalMagnitude,vector3d.Z)
xlabel('x')
ylabel('y')
drawNow(gcm)

%%
% so we are plotting the polarization in the xy-plane, or the yz-plane with

plotSection(P.directionalMagnitude,vector3d.X)
ylabel('y')
zlabel('z')
drawNow(gcm)

%% Mean Tensor Calculation
% Let us import some data, which was originally published by Mainprice, D.,
% Lloyd, G.E. and Casey , M. (1993) Individual orientation measurements in
% quartz polycrystals: advantages and limitations for texture and
% petrophysical property determinations. J. of Structural Geology, 15,
% pp.1169-1187
%

fname = fullfile(mtexDataPath,'orientation', 'Tongue_Quartzite_Bunge_Euler');

ori = orientation.load(fname,CS, 'ColumnNames', {'Euler 1' 'Euler 2' 'Euler 3'})

%%
% The figure on p.1184 of the publication

Pm = ori.calcTensor(P)

plot(Pm)
mtexColorbar

%%
%

close all
plot(Pm)
mtexColorbar

setMTEXpref('defaultColorMap',WhiteJetColorMap)

author: mtex-toolbox

repo: mtex

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


text:
classdef ArrayList < handle
    %ARRAYLIST A cell array implementation of java.util.ArrayList
    %   ArrayList(initialCapacity) mimics a fixed sized Java ArrayList
    %   using a cell array to store the data elements. initialCapacity must
    %   be greater than one and less than infinity.
    %
    %   % Example:
    %   %   Append, insert, and remove elements from the ArrayList.
    %
    %   aList = containers.ArrayList(10);
    %   aList.appendElements(1);
    %   aList.appendElements({5,6,7,8,9,10});  % Multiple elements
    %                                          % must be appended as
    %                                          % a cell array.
    %   aList.insertAtIndex(2, {2,3,4});
    %   numberTwoAndThree = aList.removeAtIndex(2:3);
    %   numberNineAndTen = aList.removeLastN(2);
    %   numberSeven = aList.removeElement(7);
    %
    %   ArrayList Methods:
    %       appendElement       - Append an element or elements to the end.
    %       elementAt           - Return the element at a specific index.
    %       insertAtIndex       - Insert an element at a specific index.
    %       isempty             - Returns true if the array is empty.
    %       lastElement         - Returns the last element.
    %       length              - Returns the number of elements.
    %       removeAtIndex       - Remove element at a specific index.
    %       removeElement       - Remove element from the array.
    %       removeFirstN        - Remove the first N elements.
    %       removeLast          - Remove the last element.
    %       removeLastN         - Remove the last N elements.
    %
    %  ArrayList Properties:
    %       Array               - Cell array storing all elements.
    %       Capacity            - The maximum, fixed capacity of the array.
    %       Count               - Number of elements in the array.

    % Copyright (C) 2013, Georgia Tech Research Corporation
    % see the LICENSE file included with this software

    properties (SetAccess = private)
        Count       % Number of elements in the array
        Capacity    % The maximum capacity of the array
    end

    properties (Access = private)
        Array       % Data structure in which elements are stored
    end

    methods

        function obj = ArrayList(initialCapacity)
            %% ARRAYLIST Constructor
            %   obj = ArrayList(initialCapacity) is the default constructor
            %   that creates an array list with a maximum, fixed capacity
            %   of initialCapacity elements.

            assert(all(imag(initialCapacity) == 0) && all(mod(initialCapacity,1) == 0), ...
                   'simiam:NotRealInteger', ...
                   'Initial capacity must be a real integer.')
            assert(initialCapacity > 0 && initialCapacity < Inf, ...
                   'simiam:OutOfBounds', ...
                   'Initial capacity must be greater than zero and less than infinity.');

            obj.Capacity = initialCapacity;
            obj.Array = cell(obj.Capacity, 1);
            obj.Count = 0;
        end

        function appendElement(obj, elementOrElements)
            %% APPENDELEMENT Append an element or elements to the end.
            %   appendElement(elementOrElements) appends an element or
            %   multiple elements to the end of the array. If multiple
            %   elements are appended, then they must be input as a cell
            %   array. An error is thrown if there is not enough space to
            %   append all of the elements.
            if iscell(elementOrElements)
                nElements = length(elementOrElements);
            else
                elementOrElements = {elementOrElements};
                nElements = 1;
            end
            assert(obj.Count+nElements <= obj.Capacity, ...
                   'simiam:NotEnoughSpace', ...
                   'Not enough space to append %d elements.', nElements);
%             assert(iscell(elementOrElements), ...
%                    'simiam:IncorrectDataType', ...
%                    'Input must be a cell array.');
            obj.Array(obj.Count+(1:nElements)) = elementOrElements;
            obj.Count = obj.Count+nElements;
        end

        function insertAtIndex(obj, elementOrElements, index)
            %% INSERTATINDEX Insert an element or elements at a specific index.
            %   insertAtIndex(elemntOrElments, index) inserts an element or
            %   multiple elements at a specific index in the range
            %   [1,min(obj.Count+1,obj.Capacity]. If multiple elements are
            %   inserted, then they must be input as a cell array. An error
            %   is thrown if there is not enough space to insert all of the
            %   elements.

            if iscell(elementOrElements)
                nElements = length(elementOrElements);
            else
                elementOrElements = {elementOrElements};
                nElements = 1;
            end
            assert(obj.Count+nElements <= obj.Capacity, ...
                   'simiam:NotEnoughSpace', ...
                   'Not enough space to insert %d elements.', nElements);
%             assert(iscell(elementOrElements), ...
%                    'simiam:IncorrectDataType', ...
%                    'Input must be a cell array.');
            obj.validateIndex(index);

            obj.Array(obj.Count+(1:nElements)) = obj.Array(index:obj.Count);
            obj.Array(index+{1:nElements}) = elementOrElements;
            obj.Count = obj.Count+nElements;
        end

        function elementOrElements = removeAtIndex(obj, indexOrIndices)
            %% REMOVEATINDEX Remove element(s) at a specific index or indicies.
            %   elementOrElements = removeAtIndex(indexOrIndicies) removes
            %   and returns the elements located at a specific index or
            %   indicies. An error is thrown if there are not enough
            %   elements in the array to remove or the indices are out of
            %   bounds.

            nElements = length(indexOrIndices);
            assert(obj.Count-nElements >= 0, ...
                   'simiam:NotEnoughElementsLeft', ...
                   'Cannot remove more elements than are currently in the array.');
            obj.validateIndex(indexOrIndices);
            elementOrElements = cell(nElements, 1);

            [sortedIndices, sortIndex] = sort(indexOrIndices, 'descend');
            for i = 1:nElements
                index = sortedIndices(i);
                elementOrElements{sortIndex(i)} = obj.Array{index};
                for j = (index+1):obj.Count
                    obj.Array{j-1} = obj.Array{j};
                end
                obj.Array{obj.Count} = [];
                obj.Count = obj.Count-1;
            end
        end

        function removeElement(obj, elementOrElements)
            %% REMOVEELEMENT Remove element(s) from the array.
            %   removeElement(elementOrElements) removes the elements
            %   specified in the input if they are found. An error is
            %   thrown if there are not enough elements in the array
            %   to remove.

            nElements = length(elementOrElements);
            assert(obj.Count-nElements >= 0, ...
                   'simiam:NotEnoughElementsLeft', ...
                   'Cannot remove more elements than are currently in the array.');
            for i = 1:obj.Count
                if(any(obj.Array(i) == elementOrElements))
                    obj.removeAtIndex(i);
                    break;
                end
            end
        end

        function element = removeLast(obj)
            %% REMOVELAST Remove the last element.
            %   element = removeLast() removes the last element from the
            %   array. An error is thrown if the array is empty.

            assert(obj.Count > 0, ...
                   'simiam:NotEnoughElementsLeft', ...
                   'Cannot remove from an empty array.');
            element = obj.removeAtIndex(obj.Count);
        end

        function elements = removeLastN(obj, nElements)
            %% REMOVELASTN Removes the last N elements.
            %   elements = removeLastN(nElements) removes the last N
            %   elements from the array. An error is thrown if there are
            %   not enough elements to remove.

            assert(obj.Count-nElements >=0, ...
                   'simiam:NotEnoughElementsLeft', ...
                   'Cannot remove more elements than are currently in the array.');
            elements = obj.removeAtIndex((obj.Count-nElements):obj.Count);
        end

        function elements = removeFirstN(obj, nElements)
            %% REMOVEFIRSTN Removes the first N elements
            %   elements = removeFirstN(nElements) removes the first N
            %   elements from the array. An error is thrown if there are
            %   not enough elements to remove.

            assert(obj.Count-nElements >= 0, ...
                   'simiam:NotEnoughElementsLeft', ...
                   'Cannot remove more elements than are currently in the array.');
            elements = obj.removeAtIndex(1:nElements);
        end

        function count = length(obj)
            %% LENGTH Returns the number of elements.
            %   count = length() returns the number of elements currently
            %   stored in the array.

            count = obj.Count;
        end

        function bool = isempty(obj)
            %% ISEMPTY Returns true if the array is empty.
            %   bool = isempty() returns true if the array is empty, false
            %   otherwise.

            bool = (obj.Count == 0);
        end

        function element = lastElement(obj)
            %% LASTELEMENT Returns the last element.
            %   element = lastElement() returns the last element in the
            %   array. An error is thrown if there are no elements stored
            %   in the array.

            assert(obj.Count > 0, ...
                   'simiam:NotEnoughElementsLeft', ...
                   'Cannot fetch from an empty array.');
            element = obj.Array{obj.Count};
        end

        function elementOrElements = elementAt(obj, indexOrIndices)
            %% ELEMENTAT Return the element at the specific index.
            %   elementOrElements = elementAt(indexOrIndicies)returns the
            %   element or elements at the specific index or indicies. An
            %   error is thrown if the array is empty.

            assert(obj.Count > 0, ...
                   'simiam:NotEnoughElementsLeft', ...
                   'Cannot fetch from an empty array.');
            obj.validateIndex(indexOrIndices);
            if length(indexOrIndices) > 1
                elementOrElements = obj.Array(indexOrIndices);
            else
                elementOrElements = obj.Array{indexOrIndices};
            end
        end

        function replaceElementAt(obj, element, index)
            %% REPLACEELEMENTAT Replaces the element at the specific index.
            %   replaceElementAt(element, index) replaces the
            %   element or elements at the specific index or indicies. An
            %   error is thrown if the array is empty.

            assert(obj.Count > 0, ...
                   'simiam:NotEnoughElementsLeft', ...
                   'Cannot fetch from an empty array.');
            obj.validateIndex(index);
            obj.Array{index} = element;
        end

        function elementOrElements = allElements(obj)
            assert(obj.Count > 0, ...
                   'simiam:NotEnoughElementsLeft', ...
                   'Cannot fetch from an empty array.');
            elementOrElements = obj.Array(1:obj.Count);
        end

        function disp(obj)
            %% DISP Display the cell array.
            %   disp() displays the contents of the cell array.

            disp(obj.Array);
        end

    end

    methods (Access = private)

        function validateIndex(obj, indexOrIndices)
            %% VALIDATEINDEX Validate the index based on the properties of the array.
            %   validateIndex(indexOrIndices) throws an error if the index
            %   or indices are not positive, real integers than are
            %   not within the bounds of the array.

            assert(all(imag(indexOrIndices) == 0) && all(mod(indexOrIndices,1) == 0), ...
                   'simiam:NotRealInteger', ...
                   'index must be a real integer.')
            assert(all(indexOrIndices > 0) && all(indexOrIndices <= min(obj.Count+1, obj.Capacity)), ...
                   'simiam:OutOfBounds', ...
                   'index must be within bounds of the array.');
        end

    end
end

author: jdelacroix

repo: simiam

save_path: github-repos/MATLAB/jdelacroix-simiam


text:
% just run it.
%
% Prompts you for all necessary info, including image to draw on
% and output image file name.
%
% Writes a mask *img file in analyze format that you can view in spm
%
% Functions called:
% c:\tor_scripts\voistatutility\readim2.m
%   read_hdr.m
%
% c:\tor_scripts\roiutility\imgmanip\getvoxels3.m
%   C:\matlabR12\toolbox\matlab\elmat\ind2sub.m
%
% c:\tor_scripts\roiutility\cluster\tor_extract_rois.m
%
% C:\matlabR12\toolbox\matlab\datatypes\squeeze.m

% c:\tor_scripts\voistatutility\nanmean.m
% center_of_mass.m

% C:\matlabR12\toolbox\matlab\elmat\repmat.m
% (calls other spm functions)
%
% c:\tor_scripts\fixedfxutility\voxel2mm.m
%
% spm_vol.m
% spm_read_vols.m
%
% ..
%    by Tor Wager  last modified 9/22/02
% ..
%
% :Note: readim2 works on little endian systems (Linux).
% If you're using Unix, you should modify this script
% to call the function readim2_b.m instead.
%

P = spm_select(1,'image','Choose whole-brain anatomical mask');
%P = spm_get(1,'*.img','Choose whole-brain anatomical mask');

V = spm_vol(P);
vol = spm_read_vols(V);
vol = double(vol);

dosagg = input('Press 1 to choose saggital slices, or 0 for axial.');
% if saggital, use vol, not rvol
% y is not reversed
% ginput x = brain y, ginput y = brainz, slice = x

if dosagg
    [rvol2] = readim2(vol,'p','sagg');
else
    for i = 1:size(vol,3)
        rvol(:,:,i) = rot90(vol(:,:,i));
    end

    [rvol2] = readim2(rvol,'p');
end

colormap gray
wslices = input('Enter range of slices (e.g, 20:30) ');
close

if dosagg
    [rvol2,hdr,h] = readim2(vol,'p','sagg','noflipy',wslices);
    colormap gray
    [voxels,mask] = getvoxels3(h,vol,wslices,'sagg');
else
    [rvol2,hdr,h] = readim2(rvol,'p','ax','noflipy',wslices);
    colormap gray
    [voxels,mask] = getvoxels3(h,vol,wslices);
end

%for i = 1:size(mask,3)
%    rmask(:,:,i) = rot90(mask(:,:,i),3);
%end

V.fname = input('Enter output filename, no quotes: ','s');
spm_write_vol(V,mask);
[d,fname,e] = fileparts(V.fname);
saveas(gcf,fname,'fig')

% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% cluster structure stuff for ROIs
% +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

CLU.M = V.mat;
CLU.voxSize = diag(V.mat)';
CLU.voxSize = CLU.voxSize(1:3);
CLU.VOX = CLU.voxSize;
CLU.XYZ = voxels';
CLU.XYZmm = voxel2mm(voxels',V.mat);
CLU.Z = ones(1,size(CLU.XYZ,2));
CLU.crit_t = 1;
CLU.cl_size = 0;
CLU.u = CLU.crit_t;
CLU.k = CLU.cl_size;
CLU.title = fname;

clusters = tor_extract_rois([],CLU,CLU);
eval(['save ' fname '_clusters CLU clusters'])

disp(' ')
newim = input('Press 1 to choose a display image, or anything else to use your prior mask.');
if newim, P = spm_get(1,'*.img','Choose overlay image');,end

spm_image('init',P)
for i = 1:length(clusters)
    spm_orthviews('AddColouredBlobs',1,clusters(i).XYZ,clusters(i).Z,CLU.M,rand(1,3))
end


text:
function loadings = computeLoadingPercentPower(V,W,H)
    loadings = [];
    K = size(H,1);
    varv = sum(V(:).^2);
    for fi = 1:K
        WH = helper.reconstruct(W(:,fi,:),H(fi,:));
        loadings(fi) = sum(2*V(:).*WH(:) - WH(:).^2)/varv;
    end
    loadings(loadings<0)=0;
end

author: FeeLab

repo: seqNMF

save_path: github-repos/MATLAB/FeeLab-seqNMF


text:
function [h, compUpAV, compUpVP, compUp] =  lfmComputeH3AV(gamma1_p, gamma1_m, sigma2, t1, ...
    t2, preFactor, mode)

% LFMCOMPUTEH3AV Helper function for computing part of the LFMAV kernel.
% FORMAT
% DESC computes a portion of the LFMAV kernel.
% ARG gamma1 : Gamma value for first system.
% ARG gamma2 : Gamma value for second system.
% ARG sigma2 : length scale of latent process.
% ARG t1 : first time input (number of time points x 1).
% ARG t2 : second time input (number of time points x 1).
% ARG preFactor : precomputed constants.
% ARG mode: indicates the correct derivative.
% RETURN h : result of this subcomponent of the kernel for the given values.
%
% COPYRIGHT : Mauricio Alvarez, 2010

% KERN

% Evaluation of h

if nargout>1
    [compUpAV{1}, compUpVP{1}, compUp{1}] = lfmavComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2, mode);
    [compUpAV{2}, compUpVP{2}, compUp{2}] = lfmavComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
    h = preFactor(1)*compUpAV{1} + preFactor(2)*compUpAV{2};
else
    h = preFactor(1)*lfmavComputeUpsilonMatrix(gamma1_p,sigma2, t1,t2, mode) ...
        + preFactor(2)*lfmavComputeUpsilonMatrix(gamma1_m,sigma2, t1,t2, mode);
end

author: SheffieldML

repo: GPmat

save_path: github-repos/MATLAB/SheffieldML-GPmat


text:
function planC = reviveRS(planC,planCDir,numExt,scanNum)
% Usage: planC = reviveRS(planC,planCDir,numExt,scanNum)
% Extract structure masks in planC and re-import them to eliminate inconsistencies in older versions of saved structures.
% EML 2020-11-26

if ~exist('numExt','var')
    numExt = 0;
end

if ~exist('scanNum','var')
    scanNum = 1;
end

indexS = planC{end};
structureListC = {planC{indexS.structures}.structureName};

for i = 1:numel(structureListC)
    strName = structureListC{i};
    disp(['Processing structure: ' strName]);
    updatedStructureListC = {planC{indexS.structures}.structureName};
    structNum = getMatchingIndex(lower(strName),lower(updatedStructureListC),'exact');
    strMask3M = getStrMask(structNum, planC);

    planC = deleteStructure(planC,structNum);
    planC = maskToCERRStructure(strMask3M,0,scanNum,strName,planC);
end
disp('Generating DICOM UID Relationships');
planC = generate_DICOM_UID_Relationships(planC);
disp('Exporting RTSTRUCT');
export_RS_IOD(planC,planCDir,numExt);


text:
function [ n_data, x, cbi ] = airy_cbi_values ( n_data )

%*****************************************************************************80
%
%% AIRY_CBI_VALUES returns some values of the Airy Bi(x) with complex argument.
%
%  Discussion:
%
%    The Airy functions Ai(X) and Bi(X) are a pair of linearly independent
%    solutions of the differential equation:
%
%      W'' - X * W = 0
%
%    In Mathematica, the function can be evaluated by:
%
%      AiryBi[x]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 April 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz, Irene Stegun,
%    Handbook of Mathematical Functions,
%    National Bureau of Standards, 1964,
%    ISBN: 0-486-61272-4,
%    LC: QA47.A34.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Cambridge University Press, 1999,
%    ISBN: 0-521-64314-7,
%    LC: QA76.95.W65.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, complex X, the argument of the function.
%
%    Output, complex CBI, the value of the Airy BI function.
%
  n_max = 10;

  cbi_vec = [ ...
    1.207423594952871  + 0.0000000000000000 * i, ...
    0.9127160108293936 + 0.3800456133135556 * i, ...
    0.6824453575635721 + 0.3343047153635002 * i, ...
    0.5726265660086474 + 0.3988641086982559 * i, ...
    0.2511841251049547 + 0.3401447690712719 * i, ...
    0.1039973894969446 + 0.0000000000000000 * i, ...
    0.2511841251049547 - 0.3401447690712719 * i, ...
    0.5726265660086474 - 0.3988641086982559 * i, ...
    0.6824453575635721 - 0.3343047153635002 * i, ...
    0.9127160108293936 - 0.3800456133135556 * i ];

  x_vec = [ ...
     1.0000000000000000 + 0.0000000000000000 * i, ...
     0.8090169943749474 + 0.5877852522924731 * i, ...
     0.3090169943749474 + 0.9510565162951536 * i, ...
    -0.3090169943749474 + 0.9510565162951536 * i, ...
    -0.8090169943749474 + 0.5877852522924731 * i, ...
    -1.0000000000000000 + 0.0000000000000000 * i, ...
    -0.8090169943749474 - 0.5877852522924731 * i, ...
    -0.3090169943749474 - 0.9510565162951536 * i, ...
     0.3090169943749474 - 0.9510565162951536 * i, ...
     0.8090169943749474 - 0.5877852522924731 * i ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    x = 0.0;
    cbi = 0.0;
  else
    x = x_vec(n_data);
    cbi = cbi_vec(n_data);
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function showTheseRois = dtiROIShowList(handles)
%
%   showTheseRois = dtiROIShowList(handles);
%
%Author: Dougherty, Wandell
%Purpose
%   Based on the figure settings, including the number of ROIs and the
%   ShowMode build a vector of the ROIs that should be displayed.  Used in
%   various drawing and plotting routines.
%

if (handles.roiShowMode == 1) ...
        || isempty(handles.rois) ...
        || ((length(handles.rois) == 1) && strcmp(handles.rois.name,'Empty'))
    showTheseRois = [];
elseif (handles.roiShowMode == 2),  showTheseRois = handles.curRoi;
elseif (handles.roiShowMode == 3),  showTheseRois = (1:length(handles.rois));
elseif (handles.roiShowMode == 4),  showTheseRois = (1:length(handles.rois));
else error('Bad handles.roiShowMode');
end

return;


text:
function [z_Tji_Ground_m,r_I_IRji_m,r_I_ITji_m] = Rotation_Position_ToTireFixed(z_I_Ground_m,A_K_I,A_F_I,r_I_IAji_m,h_ji_m,re_ji_m)
%% Rotation ground displacement from inertial fixed coordinate system into tire fixed coordinate systems
% Input parameters:
% A_K_I          [---]    Rotation Matrix from Inertial Frame to VehicleFixed Frame
% A_F_I          [---]    Rotation Matrix from Inertial Frame to RearAxleFixed Frame
% h_ji_m         [m]      Suspension Height [SuspHFL SuspHFR SuspHRL SuspHRR]
% re_ji_m        [m]      Effective Wheel Radii [ReFL ReFR ReRL ReRR]
% r_I_IAji       [m]      Absolut Position of Axles in Inertial Reference Frame [r_I_IAFL r_I_IAFR r_I_IARL r_I_IARR] [3x4]
% z_Ground_I_m   [---]    z-Coordinate of Road Displacement in Inertial Reference Frame [1x4]
% Output parameters:
% z_Ground_Tji_m [---]    Transformed Displacement Vector [r_TFL_ITFL r_TFR_ITFR r_TRL_ITRL r_TRR_ITRR] [1x4]

%% Wheel center points in inertial reference frame

% Quarter Vehicle Model along VehicleFixed z-Axis: Rotation from VehicleFixed Axis System into Inertial Frame
% r_I_IRji_m = [(r_I_IAji_m(:,1) + transpose(A_K_I)*[0;0;1]*(h_ji_m(1))) (r_I_IAji_m(:,2) + transpose(A_K_I)*[0;0;1]*(h_ji_m(2))) (r_I_IAji_m(:,3) + transpose(A_K_I)*[0;0;1]*(h_ji_m(3))) (r_I_IAji_m(:,4) + transpose(A_K_I)*[0;0;1]*(h_ji_m(4)))];

% Quarter Vehicle Model along RearAxleFixed z-Axis: Rotation from RearAxleFixed Axis System into Inertial Frame
r_I_IRji_m = [(r_I_IAji_m(:,1) + transpose(A_F_I)*[0;0;1]*(h_ji_m(1))) (r_I_IAji_m(:,2) + transpose(A_F_I)*[0;0;1]*(h_ji_m(2))) (r_I_IAji_m(:,3) + transpose(A_F_I)*[0;0;1]*(h_ji_m(3))) (r_I_IAji_m(:,4) + transpose(A_F_I)*[0;0;1]*(h_ji_m(4)))];

%% Contact points of road-tire-interface in inertial reference frame

% Quarter Vehicle Model along VehicleFixed z-Axis: Rotation from VehicleFixed Axis System into Inertial Frame
% r_I_ITji_m = [(r_I_IAji_m(:,1) + transpose(A_K_I)*[0;0;1]*(re_ji_m(1)+h_ji_m(1))) (r_I_IAji_m(:,2) + transpose(A_K_I)*[0;0;1]*(re_ji_m(2)+h_ji_m(2))) (r_I_IAji_m(:,3) + transpose(A_K_I)*[0;0;1]*(re_ji_m(3)+h_ji_m(3))) (r_I_IAji_m(:,4) + transpose(A_K_I)*[0;0;1]*(re_ji_m(4)+h_ji_m(4)))];

% Quarter Vehicle Model along RearAxleFixed z-Axis: Rotation from RearAxleFixed Axis System into Inertial Frame
r_I_ITji_m = [(r_I_IAji_m(:,1) + transpose(A_F_I)*[0;0;1]*(re_ji_m(1)+h_ji_m(1))) (r_I_IAji_m(:,2) + transpose(A_F_I)*[0;0;1]*(re_ji_m(2)+h_ji_m(2))) (r_I_IAji_m(:,3) + transpose(A_F_I)*[0;0;1]*(re_ji_m(3)+h_ji_m(3))) (r_I_IAji_m(:,4) + transpose(A_F_I)*[0;0;1]*(re_ji_m(4)+h_ji_m(4)))];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% Elaborate Computation of Vertical Displacement
% % The computation of the position vector of the road-tire interface
% % requires information about the rx_I_ITji_m and ry_I_ITji_m coordinates of the
% % road-tire-interface relative to the inertial frame and the vertical
% % displacement z_Ground_I_m.
%
% %% Position Vector of Road-Tire-Interface around Orientation of the Quarter Vehicle Model
%
% % Extracting (x,y)-Coordinates of r_I_ITji_m
% % rx_I_ITji_m = r_I_ITji_m(1,:);
% % ry_I_ITji_m = r_I_ITji_m(2,:);
% %
% % Quarter Vehicle Model along VehicleFixed z-Axis: Rotation from Inertial Frame into VehicleFixed Axis System
% r_Ground_m = A_K_I*vertcat(rx_I_ITji_m,ry_I_ITji_m,+z_Ground_I_m);
%
% % Quarter Vehicle Model along RearAxleFixed z-Axis: Rotation from Inertial Frame into RearAxleFixed Axis System
% r_Ground_m = A_F_I*vertcat(rx_I_ITji_m,ry_I_ITji_m,+z_Ground_I_m);
%
% %% Orientation Change from z-down-Orientation to z-up-Orientation
% z_Ground_Tji_m = -r_Ground_m(3,:);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Robust Computation of Vertical Displacement
% A more robust computation of the Vertical Displacement is possible with
% the simplified z-displacement vector [0;0;1]*z_Ground_I_m under
% neglection of the rx_I_ITji_m and ry_I_ITji_m coordinates of the
% road-tire-interface relative to the inertial frame.

%% Vertical z-Displacement of Road-Tire-Interface around Orientation of the Quarter Vehicle Model

% Quarter Vehicle Model along VehicleFixed z-Axis: Rotation from Inertial Frame into VehicleFixed Axis System
% r_Ground_m = A_K_I*[0;0;1]*z_Ground_I_m;

% Quarter Vehicle Model along RearAxleFixed z-Axis: Rotation from Inertial Frame into RearAxleFixed Axis System
r_Ground_m = A_F_I*[0;0;1]*z_I_Ground_m;

%% Orientation Change from z-down-Orientation to z-up-Orientation
z_Tji_Ground_m = -r_Ground_m(3,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

author: TUMFTM