text:
function gui_cb(action)
%
% This file is used by FASTICAG

% This file holds the callbacks to the main window

% @(#)$Id: gui_cb.m,v 1.5 2003/09/10 10:33:41 jarmo Exp $

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variables

% Handle to the main figure
global hf_FastICA_MAIN;

% Handles for needed controls in main figure;
global ht_FastICA_mixedStatus;
global ht_FastICA_dim;
global ht_FastICA_numOfSamp;
global ht_FastICA_newDim;
global ht_FastICA_whiteStatus;
global ht_FastICA_icaStatus;
global hpm_FastICA_approach;
global he_FastICA_numOfIC;
global hpm_FastICA_g;
global hpm_FastICA_stabilization;

% Main values are stored here
global g_FastICA_mixedsig;
global g_FastICA_mixedmean;
global g_FastICA_pca_D;
global g_FastICA_pca_E;
global g_FastICA_white_sig;
global g_FastICA_white_wm;
global g_FastICA_white_dwm;
global g_FastICA_ica_sig;
global g_FastICA_ica_A;
global g_FastICA_ica_W;
global g_FastICA_initGuess;
global g_FastICA_approach;
global g_FastICA_numOfIC;
global g_FastICA_g;
global g_FastICA_finetune;
global g_FastICA_a1;
global g_FastICA_a2;
global g_FastICA_myy;
global g_FastICA_stabilization;
global g_FastICA_epsilon;
global g_FastICA_maxNumIte;
global g_FastICA_maxFinetune;
global g_FastICA_sampleSize;
global g_FastICA_initState;
global g_FastICA_displayMo;
global g_FastICA_displayIn;
global g_FastICA_verbose;

% String values are here
global c_FastICA_appr_strV;
global c_FastICA_g1_strD;
global c_FastICA_g1_strV;
global c_FastICA_g2_strD;
global c_FastICA_g2_strV;
global c_FastICA_finetune_strD;
global c_FastICA_finetune_strV;
global c_FastICA_stabili_strV;
global c_FastICA_iSta_strV;
global c_FastICA_dMod_strV;
global c_FastICA_verb_strV;

% What is the load type of load dialog
global g_FastICA_loadType;

% Global variable for stopping the ICA calculations
global g_FastICA_interrupt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% What ever we do, it will take some time... not much, but some :-)
watchonInFigure = watchon;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch action
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'InitAll'

 % If the data is already loaded, then get the information from data
 % and show to the user (also set g_FastICA_numOfIC)
 if ~isempty(g_FastICA_mixedsig)
   set(ht_FastICA_mixedStatus, 'String', '');
   [dim, numofsamp] = size(g_FastICA_mixedsig);
   set(ht_FastICA_dim, 'String', int2str(dim));
   set(ht_FastICA_numOfSamp, 'String', int2str(numofsamp));
   set(ht_FastICA_newDim, 'String', int2str(dim));
   set(he_FastICA_numOfIC, 'String', int2str(dim));
   g_FastICA_numOfIC = dim;
   g_FastICA_mixedmean = [];
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'LoadData'

 handle = findobj('Tag','f_FastICALoad');  % Check if the window is already
 if isempty(handle)                        % open. If not then open it.
   pos = get(hf_FastICA_MAIN, 'Position');
   % Based on the feedback obtained from some users, it seems
   % that at least in some systems, pos can sometimes be empty. A
   % similar check is done a few lines below.
   if ~isempty (pos),
     gui_l(pos(1), pos(2));
   else
     gui_l (0, 0);
   end
 else
   if strcmp(g_FastICA_loadType, 'data')   % Check if it was the same load
     figure(handle);                       % window. If it wasn't then
   else                                    % close the other window first
     close(handle);                        % and then open the load window
     fprintf('''Load initial guess'' -dialog closed!\n');
     pos = get(hf_FastICA_MAIN, 'Position');
     if ~isempty (pos),
       gui_l(pos(1), pos(2));
     else
       gui_l (0, 0);
     end
   end
 end

 % gui_cb NewData; - is called from the load function if not canceled...

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'NewData'

 % New data is loaded or the old data changed. We need to find out
 % somethings about the new data... and do some other stuff also...
 [dim, numofsamp] = size(g_FastICA_mixedsig);
 set(ht_FastICA_dim, 'String', dim);
 set(ht_FastICA_newDim, 'String', dim);
 set(ht_FastICA_numOfSamp, 'String', numofsamp);
 set(he_FastICA_numOfIC, 'String', int2str(dim));

 g_FastICA_numOfIC = dim;    % Default for numOfIC = the new dimension
			     % PCA needs to be calculated again.
 g_FastICA_pca_E = [];       % We use this to check if PCA is calculated
 gui_cb NullWhite;           % Whitening needs to be done again also

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'NullWhite'

 % Whitening needs to done again next time it's needed
 g_FastICA_white_sig = [];   % We use this to check if whitening is calculated
 gui_cb NullICA;             % The FPICA must be calculated again

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'NullICA'

 % If IC's are needed they have to bee calculated again
 g_FastICA_ica_sig = [];     % We use this to check if FPICA is calculated
 set(ht_FastICA_icaStatus,'String','Not yet done');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'Transpose'

 if isempty(g_FastICA_mixedmean)
   g_FastICA_mixedsig = g_FastICA_mixedsig';
 else
   g_FastICA_mixedsig = (g_FastICA_mixedsig + ...
			 g_FastICA_mixedmean * ...
			 ones(1,size(g_FastICA_mixedsig, 2)))';
   g_FastICA_mixedmean = [];
 end
 gui_cb NewData;             % Data has been changed

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'DoPCA'

 if ~isempty(g_FastICA_mixedsig)
   % We'll remove mean of the data here also just in case...
   if isempty(g_FastICA_mixedmean)
     [g_FastICA_mixedsig, g_FastICA_mixedmean] = remmean(g_FastICA_mixedsig);
   end

   % Do PCA interactively: ask the user for eigenvalues
   [g_FastICA_pca_E, g_FastICA_pca_D] = pcamat(g_FastICA_mixedsig, ...
					       0, 0, 'gui', ...
					       deblank(c_FastICA_verb_strV(g_FastICA_verbose,:)));

   newdim = size(g_FastICA_pca_D, 1);
   set(ht_FastICA_newDim, 'String', int2str(newdim));
   set(he_FastICA_numOfIC, 'String', int2str(newdim));
   g_FastICA_numOfIC = newdim;
   gui_cb NullWhite;           % Whitening needs to be done again also
			       % but we'll do it when it's needed.
 else
   fprintf('Data not loaded yet!\n\n');
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'OrigDim'

 gui_cb NewData;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'ShowMixed'

 if ~isempty(g_FastICA_mixedsig)
   handle = findobj('Tag','f_FastICA_mix');  % Check if the window is already
   if isempty(handle)                        % open. If not then open it.
     figure('Tag', 'f_FastICA_mix', ...
	    'Name', 'FastICA: Plot data', ...
	    'NumberTitle', 'off');
   else
     figure(handle);
     clf;		% clear the figure for next plots
   end
   if isempty(g_FastICA_mixedmean)
     icaplot('dispsig',g_FastICA_mixedsig, 0, 0, 0, 'Mixed signals');
   else
     icaplot('dispsig',g_FastICA_mixedsig + g_FastICA_mixedmean * ...
	     ones(1, size(g_FastICA_mixedsig, 2)), 0, 0, 0, 'Mixed signals');
   end
 else
   fprintf('Data not loaded yet!\n\n');
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'ShowWhite'

 if ~isempty(g_FastICA_mixedsig)
   if isempty(g_FastICA_white_sig)     % if whitening is not done, we need to
     gui_cb Whiten;                    % do it before we can display the
   end                                 % whitened signals

   handle = findobj('Tag','f_FastICA_white');  % Check if the window is already
   if isempty(handle)                          % open. If not then open it.
     figure('Tag', 'f_FastICA_white', ...
	    'Name', 'FastICA: Plot whitened', ...
	    'NumberTitle', 'off');
   else
     figure(handle);
     clf;		% clear the figure for next plots
   end
   icaplot('dispsig',g_FastICA_white_sig,0,0,0,'Whitened signals');
 else
   fprintf('Data not loaded yet!\n\n');
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'Whiten'

 set(ht_FastICA_whiteStatus,'String','Computing...');

 % If PCA is not calculated, we'll have to calculate it now,
 % we'll do it without guestions - we don't reduce the dimension
 % here - but PCAMAT might reduce the dimension automatically.
 if isempty(g_FastICA_pca_E)
   % We'll remove mean of the data here also just in case...
   if isempty(g_FastICA_mixedmean)
     [g_FastICA_mixedsig, g_FastICA_mixedmean] = remmean(g_FastICA_mixedsig);
   end

   [g_FastICA_pca_E, g_FastICA_pca_D] = pcamat(g_FastICA_mixedsig, 1, ...
					       size(g_FastICA_mixedsig, ...
						    1), 'off', ...
					       deblank(c_FastICA_verb_strV(g_FastICA_verbose,:)));

   % Check if the dimension was reduced automatically
   newdim = size(g_FastICA_pca_D, 1);
   set(ht_FastICA_newDim, 'String', int2str(newdim));
   % Check if the numOfIC now has illegal value entered
   % We do that by telling the program that there is new value
   % entered for NumOfIC.
   gui_cb ChangeNumOfIC;
 end

 % And now we can calculate whitening...
 [g_FastICA_white_sig, g_FastICA_white_wm, g_FastICA_white_dwm] = ...
								  whitenv(g_FastICA_mixedsig, g_FastICA_pca_E, g_FastICA_pca_D, deblank(c_FastICA_verb_strV(g_FastICA_verbose,:)));

 set (ht_FastICA_whiteStatus,'String','');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'ChangeApproach'

 % Get the old value for g
 eval(['g_str = c_FastICA_g' int2str(g_FastICA_approach) '_strV;']);
 old_g = deblank(g_str(g_FastICA_g,:));

 % Get and set the new value for approach
 g_FastICA_approach = get(hpm_FastICA_approach, 'Value');

 % The possible values for g depend on the value of approach...
 eval(['g_str = c_FastICA_g' int2str(g_FastICA_approach) '_strD;']);
 set(hpm_FastICA_g, 'String', g_str);

 % Match the old g value from the new g values so that if the
 % old_g can be found from the new values (anywhere), then set new g
 % to that value, and if it's not found then set the new value to 1.
 match = 0;
 eval(['g_str = c_FastICA_g' int2str(g_FastICA_approach) '_strV;']);
 for i=1:size(g_str,1)
   if strcmp(old_g, deblank(g_str(i,:)))
     match = i;
   end
 end
 if match == 0
   match = 1;   % the old g is not availabe anymore, set g = 1.
 end
 g_FastICA_g = match;
 set(hpm_FastICA_g, 'Value', match);

 gui_cb NullICA;      % The options are changed so we must calculate ICA again

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'ChangeNumOfIC'

 % Get the new value... and store it later on after some checks
 numofic = str2num(get(he_FastICA_numOfIC, 'String'));

 % The number of IC can't be less than 1 or more than the reduced dimension.
 numoficmax = str2num(get(ht_FastICA_newDim, 'String'));
 if numofic < 1
   set(he_FastICA_numOfIC, 'String', '1');
   g_FastICA_numOfIC = 1;
 elseif numofic > numoficmax
   set(he_FastICA_numOfIC, 'String', int2str (numoficmax));
   g_FastICA_numOfIC = numoficmax;
 else
   g_FastICA_numOfIC = numofic;
 end

 gui_cb NullICA;      % The options are changed so we must calculate ICA again

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'ChangeG'

 % Get the new value for g.
 g_FastICA_g = get(hpm_FastICA_g, 'Value');

 gui_cb NullICA;      % The options are changed so we must calculate ICA again

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'ChangeStab'

 % Get the new value for g.
 g_FastICA_stabilization = get(hpm_FastICA_stabilization, 'Value');

 gui_cb NullICA;      % The options are changed so we must calculate ICA again

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'AdvOpt'

 handle = findobj('Tag','f_FastICAAdvOpt');
 if isempty(handle)                        % Check to see if the window is
   pos = get(hf_FastICA_MAIN, 'Position');     % already open...
   if ~isempty (pos),
     gui_adv(pos(1), pos(2));
   else
     gui_adv(0, 0);
   end
 else
   figure(handle)
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'ShowICASig'

 if ~isempty(g_FastICA_mixedsig)
   % If the IC's are not already calculated, we'll do it now
   if isempty(g_FastICA_ica_sig)
     gui_cb DoFPICA;
   end

   % The signals may have been already displaued by the FPICA function
   % BUT the FPICA may also have shown either the basis of the filters
   % so the signals still need to be shown - besides the mean was added
   % in only later after FPICA

   % Also notice that in this version if there was something wrong in FPICA
   % Then the results are []. - We don't try to plot them!
   if ~isempty(g_FastICA_ica_sig')
     handle = findobj('Tag','f_FastICA_ica');  % Check if the window is already
     if isempty(handle)                        % open. If not then open it.
       figure('Tag', 'f_FastICA_ica', ...
	      'Name', 'FastICA: Plot ICs', ...
	      'NumberTitle', 'off');
     else
       figure(handle);
       clf;		% clear the figure for next plots
     end

     icaplot('dispsig', g_FastICA_ica_sig, 0, ...
	     0, 0, 'Independent components');
   end
 else
   fprintf('Data not loaded yet!\n\n');
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'DoFPICA'

 gui_cb DisableButtons;
 g_FastICA_interrupt = 0;

 if ~isempty(g_FastICA_mixedsig)
   if isempty(g_FastICA_white_sig)     % We need the whitened signal here
     gui_cb Whiten;
   end

   set(ht_FastICA_icaStatus,'String','Computing...');

   % The possible values for g depend on approach
   eval(['g_str = c_FastICA_g' int2str(g_FastICA_approach) '_strV;']);

   % We'll contruct a command string which we'll later on evaluate...
   % This is where the Fixed point algorithm is used.
   command_str = ['[g_FastICA_ica_A,g_FastICA_ica_W]=' ...
		  'fpica(g_FastICA_white_sig,' ...
		  'g_FastICA_white_wm,' ...
		  'g_FastICA_white_dwm,' ...
		  '''' deblank(c_FastICA_appr_strV(g_FastICA_approach,:)) ...
		  ''',' ...
		  'g_FastICA_numOfIC,' ...
		  '''' deblank(g_str(g_FastICA_g,:)) ''',' ...
		  '''' ...
		  deblank(c_FastICA_finetune_strV(g_FastICA_finetune,:)) ...
		  ''',' ...
		  'g_FastICA_a1,' ...
		  'g_FastICA_a2,' ...
		  'g_FastICA_myy,' ...
		  '''' ...
		  deblank(c_FastICA_stabili_strV(g_FastICA_stabilization,:)) ...
		  ''',' ...
		  'g_FastICA_epsilon,' ...
		  'g_FastICA_maxNumIte,' ...
		  'g_FastICA_maxFinetune,' ...
		  '''' deblank(c_FastICA_iSta_strV(g_FastICA_initState,:)) ...
		  ''',' ...
		  'g_FastICA_initGuess,' ...
		  'g_FastICA_sampleSize,' ...
		  '''' deblank(c_FastICA_dMod_strV(g_FastICA_displayMo,:)) ...
		  ''',' ...
		  'g_FastICA_displayIn,' ...
		  '''' deblank(c_FastICA_verb_strV(g_FastICA_verbose,:)) ...
		  ''');'];

   % If the user wants to plot while computing...
   % let's at least plot it to the right figure then
   if ~strcmp(deblank(c_FastICA_dMod_strV(g_FastICA_displayMo,:)),'off')
     handle = findobj('Tag','f_FastICA_ica');  % Check if the window is already
     if isempty(handle)                        % open. If not then open it.
       figure('Tag', 'f_FastICA_ica', ...
	      'Name', 'FastICA: Plot ICs', ...
	      'NumberTitle', 'off');
     else
       figure(handle);
       clf;		% clear the figure for next plots
     end
   end

   % ... and so let's do it...
   eval(command_str);

   % Also notice that in this version if there was something wrong in FPICA
   % Then the results are [].
   if ~isempty(g_FastICA_ica_W)
     % Add the mean back in.
     g_FastICA_ica_sig = g_FastICA_ica_W * g_FastICA_mixedsig ...
			 + (g_FastICA_ica_W * g_FastICA_mixedmean) ...
			 * ones(1,size(g_FastICA_mixedsig, 2));
     set (ht_FastICA_icaStatus,'String','Done');
   else
     gui_cb NullICA;  % set icasig=[] and do what ever needs to be done then
   end

   if ~(g_FastICA_interrupt)
     gui_cb EnableButtons;
   end
 else
   fprintf('Data not loaded yet!\n\n');
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'Interrupt'
 g_FastICA_interrupt = 1;
 set(ht_FastICA_icaStatus,'String','Interrupted');
 gui_cb EnableButtons;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'DisableButtons'
 set(findobj('Tag','b_Transpose'),'Enable','off');
 set(findobj('Tag','b_ShowMixed'),'Enable','off');
 set(findobj('Tag','b_DoPCA'),'Enable','off');
 set(findobj('Tag','b_OrigDim'),'Enable','off');
 set(findobj('Tag','b_ShowWhite'),'Enable','off');
 set(findobj('Tag','b_advOpt'),'Enable','off');
 set(findobj('Tag','b_ShowICASig'),'Enable','off');
 set(findobj('Tag','b_LoadData'),'Enable','off');
 set(findobj('Tag','b_DoFPICA'),'Enable','off');
 set(findobj('Tag','b_SaveData'),'Enable','off');
 set(findobj('Tag','b_Quit'),'Enable','off');
 set(findobj('Tag','b_Interrupt'),'Visible','on');
 drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'EnableButtons'
 set(findobj('Tag','b_Transpose'),'Enable','on');
 set(findobj('Tag','b_ShowMixed'),'Enable','on');
 set(findobj('Tag','b_DoPCA'),'Enable','on');
 set(findobj('Tag','b_OrigDim'),'Enable','on');
 set(findobj('Tag','b_ShowWhite'),'Enable','on');
 set(findobj('Tag','b_advOpt'),'Enable','on');
 set(findobj('Tag','b_ShowICASig'),'Enable','on');
 set(findobj('Tag','b_LoadData'),'Enable','on');
 set(findobj('Tag','b_DoFPICA'),'Enable','on');
 set(findobj('Tag','b_SaveData'),'Enable','on');
 set(findobj('Tag','b_Quit'),'Enable','on');
 set(findobj('Tag','b_Interrupt'),'Visible','off');
 drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'SaveData'

 handle = findobj('Tag','f_FastICASave');  % Check if the window is already
 if isempty(handle)                        % open. If not then open it.
   pos = get(hf_FastICA_MAIN, 'Position');
   if ~isempty (pos),
     gui_s(pos(1), pos(2));
   else
     gui_s(0, 0);
   end
 else
   figure(handle);                       % window. If it wasn't then
 end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'Quit'

 % We'll close the other dialogs if they are open.
 Tags = ['f_FastICALoad  '
	 'f_FastICAAdvOpt'
	 'f_FastICASave  '
	 'f_FastICA_mix  '
	 'f_FastICA_white'
	 'f_FastICA_ica  '];
 for i=1:size(Tags,1)
   handle = findobj('Tag', deblank(Tags(i,:)));
   if ~isempty(handle)
     close(handle);
   end
 end

 % Close this window
 close(hf_FastICA_MAIN);

 % Clear the used global variables.
 gui_cg;

 % Use return to avoid reaching the watchoff statement at the end
 % (There used to be a 'break' statement here, but it resulted in
 % errors in more recent version of Matlab -- jarmo)
 return;
 % ... and we're done.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'About'

 gui_help('gui_cb_about');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
case 'Help'

 gui_help('gui_cb_help');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end    % switch

watchoff (watchonInFigure);

author: angeloyeo

repo: gongdols

save_path: github-repos/MATLAB/angeloyeo-gongdols


text:
% Fig. 5.19   Feedback Control of Dynamic Systems, 6e
%             Franklin, Powell, Emami
%
clf
n=1;
d=conv([1 2 0],[1 2 5]);
pzmap(n,d)
axis([-4 4 -3 3])
title('Fig.5.19 Departure angle calculation')
z=0:.1:.9;
 wn=1:1:4;
 sgrid(z, wn)
 hold off

author: Sable

repo: mcbench-benchmarks


text:
function pause(H)
%PAUSE Pause playback of movie in movie player
%  Note that the play and pause methods are identical.
%  Pause is equivalent to calling play.

feval(H.fcns.play,[],[],H.hfig);

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [Problem] = CMATpick_point(...
    Problem,Constraints,Spec,Objectives,Haxes)
satisfied = false;
while ~satisfied
    disp('Choose new point by clicking')
    [x,y] = ginput(1);
    newpth = plot(gca,x,y,'*k','MarkerEdgeColor','k',...
        'MarkerFaceColor','b','MarkerSize',9);

    %%  Change deceision vector to new point
    TestProb = Problem;
    ind = Haxes==gca;
    [ii,jj] = find(ind);
    testsolX = Problem.x0(Problem.activeX);
    testsolX(ii) = y;
    testsolX(jj) = x;
    TestProb.x0(Problem.activeX) = testsolX;
    %%
    commandwindow;
    if get_yes_or_no(...
            'New design point plotted. Evaluate chosen point? [Y]/N: ',true)
        ITERevaluate_X0(TestProb,Constraints,Spec,Objectives,false);
    end
    whattodo = txtmenu('',...
        'Choose new point (zoom/pan now)','Use chosen point',...
        'Cancel - use old design point');
    switch whattodo
        case 1 %satisfied
            Problem.x0 = TestProb.x0;
            satisfied = true;
        case 0 %try again
            delete(newpth);
        case 2 %use old
            delete(newpth);
            satisfied = true;
    end
end

end


text:
 function job = calcGraph(job, varargin)
 % set up similarity graph for parent grain reconstruction
 %
 % Syntax
 %   job.calcGraph
 %
 % Input
 %  job - @parentGrainReconstructor
 %
 % Output
 %  job.graph - adjacency matrix of the graph
 %
 % Options
 %  threshold - misfit at which the probability is set to 0.5, default is 2 degree
 %  tolerance - range around the threshold where the probability increases from 0 to 1
 %  noC2C     - ignore child to child grain boundaries
 %  noP2C     - ignore parent to child grain boundaries
 %
 % Description
 % The weights of the graph are computed from a cummulative gaussion
 % distribution with mean given by the option |'threshold'| and variance
 % given by the option |'tolerance'|
 %

 threshold = get_option(varargin,'threshold',2*degree);
 tol = get_option(varargin,'tolerance',1.5*degree);

 % OR fit of grain neighbors
 [fit, grainPairs] = calcGBFit(job,varargin{:});

 % turn into probability
 prob = 1 - 0.5 * (1 + erf(2*(fit - threshold)./tol));

 % write into similarity matrix
 job.graph = sparse(grainPairs(:,1),grainPairs(:,2),prob,...
   length(job.grains),length(job.grains));

 % ensure graph is symmetric
 job.graph = max(job.graph, job.graph.');

 end

author: mtex-toolbox

repo: mtex

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


text:
function lines = strflushright(lines, offset)
%STRFLUSHRIGHT Flush right each line of text.
%
%   STRFLUSHRIGHT(LINES, OFFSET) where LINES is a cell array of lines and
%   OFFSET is a non-negative integer, does a right flush of the text in each
%   line so the last character in each line is at the specified offset.  If
%   OFFSET is omitted, 72 is used.

%   Author:      Peter J. Acklam
%   Time-stamp:  2003-10-13 11:13:53 +0200
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   error(nargchk(1, 2, nargin));
   if nargin < 2
      offset = 72;
   end

   spc = ' ';
   tab = sprintf('\t');
   fmt = sprintf('%%%ds', offset);
   lines = lines(:);

   for i = 1:length(lines)
      k = find(lines{i} ~= spc & lines{i} ~= tab);
      if isempty(k)
         lines{i} = '';
      else
         lines{i} = sprintf(fmt, lines{i}(k(1):k(end)));
      end
   end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
%% housekeeping
clc

%% rise the model

m=rise('pdk','rise_flags',{'Sektors',{'C','I'}},...
    'steady_state_file','sstate_model',...
    'saveas',true);

%% get the parameters

[p,priors]=create_parameters(~true);

m=set(m,'parameters',p);

%% get the data

data=create_data();

%% estimate model
clc

ms=estimate(m,'data',data,'estim_priors',priors,'optimizer','fminunc');

author: jmaih


text:
function [Q,c] = compileQuadratic(c,p,onlysquares)
Q = spalloc(length(c),length(c),0);
%c = p.c;
for i = 1:size(p.bilinears,1)
    quadratic = (p.bilinears(i,2)==p.bilinears(i,3));
    if onlysquares == 3 && quadratic
        if c(p.bilinears(i,1))>0
            Q(p.bilinears(i,2),p.bilinears(i,3)) = c(p.bilinears(i,1));
            c(p.bilinears(i,1)) = 0;
        end
    elseif c(p.bilinears(i,1)) && (~onlysquares || quadratic) && ~(onlysquares == 3)
        %if ~quadratic || (quadratic && c(p.bilinears(i,1))>0)
            Q(p.bilinears(i,2),p.bilinears(i,3)) = c(p.bilinears(i,1))/2;
            Q(p.bilinears(i,3),p.bilinears(i,2)) = Q(p.bilinears(i,3),p.bilinears(i,2))+c(p.bilinears(i,1))/2;
            c(p.bilinears(i,1)) = 0;
       % end
    end
end

author: yalmip

repo: YALMIP

save_path: github-repos/MATLAB/yalmip-YALMIP


text:
function z = max( x, y )

% MAX   Maximum, z = max(x,y)

if isa(x,'tfocs_tuple')
    z = x;
    if isa(y,'tfocs_tuple')
        z.value_ = cellfun( @max, x.value_, y.value_, 'UniformOutput', false );
    elseif isscalar(y)
        z.value_ = cellfun( @max, x.value_, {y}, 'UniformOutput', false );
    else
        z.value_ = cellfun( @max, x.value_, {y}, 'UniformOutput', false );
    end
else
    z = y;
    if isscalar(x)
        z.value_ = cellfun( @max, {x}, y.value_, 'UniformOutput', false );
    else
        z.value_ = cellfun( @max, {x}, y.value_, 'UniformOutput', false );
    end
end

% TFOCS v1.3 by Stephen Becker, Emmanuel Candes, and Michael Grant.
% Copyright 2013 California Institute of Technology and CVX Research.
% See the file LICENSE for full license information.

author: cvxr

repo: TFOCS

save_path: github-repos/MATLAB/cvxr-TFOCS


text:
function [ x ] = phi_x_inv( y )

global minus_log_phi_inv_table;
global min_minus_log_phi;
global max_minus_log_phi;
global increment_minus_log_phi;
minus_log_phi = -log(y);
minus_log_phi = max(minus_log_phi, min_minus_log_phi);
minus_log_phi = min(minus_log_phi, max_minus_log_phi);
minus_log_phi_index = round((minus_log_phi - min_minus_log_phi)/increment_minus_log_phi - 0.499) + 1;
x = minus_log_phi_inv_table(minus_log_phi_index);

end

author: tavildar

repo: Polar

save_path: github-repos/MATLAB/tavildar-Polar


text:
function [vacc, hacc, vcm, hcm, pg] = ...
    testSingleSegmentationsCV(imsegs, labdata, maps, vclassifier, hclassifier, ncv)
% [vacc, hacc, vcm, hcm, pg] = testSingleSegmentationsCV(
%                       imsegs, labdata, maps, vclassifier, hclassifier, ncv)

pg = cell(numel(imsegs), 1);

for f = 1:numel(imsegs)

    nsp = imsegs(f).nseg;

    c = ceil(f/numel(imsegs)*ncv);

    pg{f} = zeros(nsp, 7);

    segs = cell(nsp, 1);

    for s = 1:size(labdata{f}, 1)

        ind = find(maps{f}==s);

        if ~isempty(ind)

            vconf = test_boosted_dt_mc(vclassifier(c), labdata{f}(s, :));
            vconf = 1 ./ (1+exp(-vconf));
            vconf = vconf / sum(vconf);

            hconf = test_boosted_dt_mc(hclassifier(c), labdata{f}(s, :));
            hconf = 1 ./ (1+exp(-hconf));
            hconf = hconf / sum(hconf);

            pgs = [vconf(1) vconf(2)*hconf vconf(3)];

            pg{f}(ind, :) = repmat(pgs, numel(ind), 1);

        end

    end

end

[vacc, hacc, vcm, hcm] = mcmcProcessResult(imsegs, pg);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [segs, ind] = checksegs(segs, map, s)
% Checks whether this segment has been seen before

ind = find(map==s);

oldsegs = segs{ind(1)};

for k = 1:numel(oldsegs)
    if (numel(oldsegs{k})==numel(ind)) && all(oldsegs{k}==ind)
        ind = [];
        return;
    end
end

segs{ind(1)}{end+1} = ind;

author: zouchuhang

repo: LayoutNet


text:
function spy_ge ( m, n, a, header )

%*****************************************************************************80
%
%% SPY_GE plots a sparsity pattern for a general storage (GE) matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    16 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, N, the number of rows and columns
%    in the matrix.
%
%    Input, real A(M,N), the matrix.
%
%    Input, string HEADER, the name to be used for the
%    title of the plot, and as part of the names of the data, command
%    and plot files.
%

%
%  Create data file.
%
  data_filename = strcat ( header, '_data.txt' );
  data_unit = fopen ( data_filename, 'wt' );
  nz_num = 0;
  for j = 1 : n
    for i = 1 : m
      if ( a(i,j) ~= 0.0 )
        fprintf ( data_unit, '%d  %d\n', j, i );
        nz_num = nz_num + 1;
      end
    end
  end
  fclose ( data_unit );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Created sparsity data file "%s".\n', data_filename );
%
%  Create command file.
%
  command_filename = strcat ( header, '_commands.txt' );
  command_unit = fopen ( command_filename, 'wt' );

  fprintf ( command_unit, '# %s\n', command_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, '# Usage:\n' );
  fprintf ( command_unit, '#  gnuplot < %s\n', command_filename );
  fprintf ( command_unit, '#\n' );
  fprintf ( command_unit, 'unset key\n' );
  fprintf ( command_unit, 'set term png\n' );

  png_filename = strcat ( header, '.png' );
  fprintf ( command_unit, 'set output "%s"\n', png_filename );
  fprintf ( command_unit, 'set size ratio -1\n' );
  fprintf ( command_unit, 'set xlabel "<--- J --->"\n' );
  fprintf ( command_unit, 'set ylabel "<--- I --->"\n' );
  fprintf ( command_unit, 'set title "%d nonzeros for ''%s''"\n', nz_num, header );
  fprintf ( command_unit, 'set timestamp\n' );
  fprintf ( command_unit, 'plot [y=1:%d] [x=%d:1] "%s" with points pt 5\n', ...
    n, m, data_filename );

  fclose ( command_unit );
  fprintf ( 1, '  Created graphics command file "%s".\n', command_filename );

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function display(t)
% DISPLAY Display object properties

newline = sprintf('\n');
for i = 1:length(t)
    disp('Sqdhandle object properties: ');
    disp('---------------------------');
    get(t(i));
end;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [output,validFrameMask] = F_cmn(input_layer)
input = input_layer.a;
[D,T,N] = size(input);

if N==1
    output = CMN(input')';
    validFrameMask = [];
else
    [validFrameMask, variableLength] = getValidFrameMask(input_layer);
    if variableLength
        input2 = ExtractVariableLengthTrajectory(input, validFrameMask);
        precision = class(gather(input(1)));
        if IsInGPU(input)
            output = gpuArray.zeros(size(input), precision);
        else
            output = zeros(size(input), precision);
        end
        for i=1:N
            output(:,1:size(input2{i},2),i) = CMN(input2{i}')';
        end
    else
        input2 = reshape(permute(input, [1 3 2]), D*N,T);
        output = CMN(input2')';
        output = permute(reshape(output, D, N, T), [1 3 2]);
    end
end
end

author: singaxiong

repo: SignalGraph


text:
%% RATE OF CONVERGENCE OF WEAK GALERKIN FINITE ELEMENT METHOD
%
% This example is to show the rate of convergence of linear finite element
% approximation of the Poisson equation on the unit square with the
% following boundary conditions:
%
% - Non-empty Dirichlet boundary condition.
% - Pure Neumann boundary condition.
% - Robin boundary condition.
%
% The basis, data structure and numerical test is summarized in <a
% href="matlab:ifem Poissonfemrate">Poissonfemrate</a>.
%
% See also PoissonP2femrate, Poissonafemrate
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

clear variables
%% Setting
[node,elem] = squaremesh([0,1,0,1],0.25);
mesh = struct('node',node,'elem',elem);
option.L0 = 2;
option.maxIt = 4;
option.printlevel = 1;
option.elemType = 'WG';

%% Non-empty Dirichlet boundary condition.
option.plotflag = 1;
pde = sincosdata;
mesh.bdFlag = setboundary(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
femPoisson(mesh,pde,option);

%% Pure Neumann boundary condition.
option.plotflag = 0;
pde = sincosNeumanndata;
% pde = sincosdata;
mesh.bdFlag = setboundary(node,elem,'Neumann');
femPoisson(mesh,pde,option);

%% Pure Robin boundary condition.
option.plotflag = 0;
pde = sincosRobindata;
mesh.bdFlag = setboundary(node,elem,'Robin');
femPoisson(mesh,pde,option);

%% Conclusion
%
% The optimal rate of convergence of the H1-norm (1st order) and L2-norm
% (2nd order) is observed. No superconvergence for ||DuI-Duh||.
%
% MGCG converges uniformly in all cases.

repo: ifem


text:
function varargout = imHistogram(img, varargin)
%IMHISTOGRAM  Histogram of 2D/3D grayscale or color images
%
%   Usage
%   H = imHistogram(IMG);
%   imHistogram(IMG);
%
%   Description
%   H = imHistogram(IMG);
%   Compute the histogram of the image IMG. IMG can be a 2D or 3D image.
%   The number of bins is computed automatically depending on image type
%   for integer images, and on image min/max values for floating-point
%   images.
%   If IMG is a color image, the result is a N-by-3 array, containing
%   histograms for the red, green and blue bands in each column.
%
%   H = imHistogram(..., N);
%   Specifies the number of histogram bins. N must be a scalar>0.
%
%   H = imHistogram(..., [GMIN GMAX]);
%   Specifies the gray level extents. This can e especially useful for
%   images stored in float, or for images with more than 256 gray levels.
%
%   H = imHistogram(..., []);
%   Forces the function to compute the histogram limits from values of
%   image gray levels.
%
%   H = imHistogram(..., BINS);
%   Specifies the bin centers.
%
%   H = imHistogram(..., ROI);
%   where ROI is a binary image the same size as IMG, computes the
%   histogram only for pixels/voxels located inside of the specified region
%   of interest.
%
%   [H X] = imHistogram(...);
%   Returns the center of bins used for histogram computation.
%
%   imHistogram(IMG);
%   When called with no output argument, displays the histogram on the
%   current axis.
%
%
%   Examples
%   % Histogram of a grayscale image (similar to imhist)
%     img = imread('cameraman.tif');
%     imHistogram(img);
%
%   % RGB Histogram of a color image
%     img = imread('peppers.png');
%     imHistogram(img);
%
%   % Compute histogram of a 3D image, only for pixel with non null values
%   % (requires image processing toolbox)
%     info = analyze75info('brainMRI.hdr');
%     X = analyze75read(info);
%     imHistogram(X, X>0, 0:88)
%
%   See also
%   imhist, hist
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2010-01-27,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

% HISTORY
% 2010-09-10 code cleanup, update doc

%% Initialise default parameters

% check if image is color
colorImage = size(img, 3)==3;

% compute intensity bounds, based either on type or on image data
if isinteger(img)
    type = class(img);
    minimg = intmin(type);
    maximg = intmax(type);
else
    minimg = min(img(:));
    maximg = max(img(:));
end

% default number of histogram bins
N = 256;

% default roi is empty
roi = [];

%% Process inputs

% process each argument
while ~isempty(varargin)
    var = varargin{1};

    if isempty(var)
        % if an empty variable is given, assumes gray level bounds must be
        % recomputed from image values
        minimg = min(img(:));
        maximg = max(img(:));
    elseif isnumeric(var) && length(var)==1
        % argument is number of bins
        N = var;
    elseif isnumeric(var) && length(var)==2
        % argument is min and max of values to compute
        minimg = var(1);
        maximg = var(2);
    elseif islogical(var)
        % argument is a ROI
        roi = var;
    elseif isnumeric(var)
        % argument is value for histo bins
        x = var;
        minimg = var(1);
        maximg = var(end);
        N = length(x);
    end

    % remove processed argument from the list
    varargin(1) = [];
end

% compute bin centers if they were not specified
if ~exist('x', 'var')
    x = linspace(double(minimg), double(maximg), N);
end

%% Main processing

% compute image histogram
if ~colorImage
    % process 2D or 3D grayscale image
    h = calcHistogram(img, x, roi);
else
    % process color image: compute histogram of each channel
    h = zeros(length(x), 3);
    if ndims(img)==3
        % process 2D color image
        for i=1:3
            h(:, i) = calcHistogram(img(:,:,i), x, roi);
        end
    else
        % process 3D color image
        for i=1:3
            h(:, i) = calcHistogram(img(:,:,i,:), x, roi);
        end
    end
end

%% Process output arguments

% In case of no output argument, display the histogram
if nargout==0
    % display histogram in current axis
    if ~colorImage
        % Display grayscale histogram
        bar(gca, x, h, 'hist');
        % use jet color to avoid gray display
        colormap jet;
    else
        % display each color histogram as stairs, to see the 3 curves
        hh = stairs(gca, x, h);

        % setup curve colors
        set(hh(1), 'color', [1 0 0]); % red
        set(hh(2), 'color', [0 1 0]); % green
        set(hh(3), 'color', [0 0 1]); % blue
    end

    % setup histogram bounds
    xlim([minimg maximg]);

elseif nargout==1
    % return histogram
    varargout = {h};
elseif nargout==2
    % return histogram and x placement
    varargout = {h, x};
elseif nargout==3
    % return red, green and blue histograms as separate outputs
    varargout = {h(:,1), h(:,2), h(:,3)};
elseif nargout==4
    % return red, green and blue histograms as separate outputs as well as
    % the bin centers
    varargout = {h(:,1), h(:,2), h(:,3), x};
end

%% Utilitary functions

function h = calcHistogram(img, x, roi)
% Compute image histogram using specified bins, and eventually a region of
% interest
if isempty(roi)
    % histogram of whole image
    h = hist(double(img(:)), x);
else
    % histogram constrained to ROI
    h = hist(double(img(roi)), x);
end


text:
function [ correct_samples, total_samples, correct_un_samples, total_un_samples ] = rcnn_test_idpr_predictions( expidx, image_set, firstidx, nImgs, bVis )

fprintf('rcnn_test_idpr_predictions()\n');

if (ischar(expidx))
    expidx = str2num(expidx);
end

if (nargin < 3)
    firstidx = 1;
elseif ischar(firstidx)
    firstidx = str2num(firstidx);
end

imdb = exp2imdb(expidx, image_set);

if (nargin < 4)
    nImgs = length(imdb.image_ids);
elseif ischar(nImgs)
    nImgs = str2num(nImgs);
end

if (nargin < 5)
    bVis = false;
end

fprintf('expidx: %d\n',expidx);
fprintf('firstidx: %d\n',firstidx);
fprintf('nImgs: %d\n',nImgs);

p = exp_params(expidx);
exp_dir = [p.expDir '/' p.shortName];
fprintf('exp_dir: %s\n',exp_dir);

if (isfield(p, 'featExpDir'))
    feat_exp_dir = p.featExpDir;
else
    feat_exp_dir = exp_dir;
end

cache_name   = 'v1_finetune_iter_70k';

if (isfield(p,'pwModelExpDir'))
    model_exp_dir = p.pwModelExpDir;
else
    model_exp_dir = exp_dir;
end

conf = rcnn_config('sub_dir', '/cachedir/train', 'exp_dir', model_exp_dir);
net_file = rcnn_get_net_filename(p.netsDir);
assert(~isempty(net_file));

lastidx = firstidx + nImgs - 1;
if (lastidx > length(imdb.image_ids))
    lastidx = length(imdb.image_ids);
end

if (firstidx > lastidx)
    return;
end

conf = rcnn_config('sub_dir', ['/cachedir/' imdb.image_set], 'exp_dir', exp_dir);

joints = zeros(length(firstidx:lastidx), length(p.pidxs), 2);

num_classes = sum(p.idpr_num_clusters)+1;
skip_class = num_classes;

top_k = 2;

total_samples = zeros(num_classes, 1);
correct_samples = zeros(num_classes, 1);
correct_topk_samples = zeros(num_classes, 1);

num_joints = 14;
total_un_samples = zeros(num_joints+1, 1);
correct_un_samples = zeros(num_joints+1, 1);

cumul_num_clusters = zeros(1, length(p.idpr_num_clusters));
for i = 1:length(p.idpr_num_clusters)
    cumul_num_clusters(i) = sum(p.idpr_num_clusters(1:i));
end

for i = firstidx:lastidx

    tic
    fprintf('%s: test (%s) %d/%d\n', procid(), imdb.name, i, length(imdb.image_ids));
    file = sprintf('%s/feat_cache/%s/%s/%s', feat_exp_dir, cache_name, imdb.image_set, imdb.image_ids{i});
    d = rcnn_load_cached_pool5_features(file);

    assert(~isempty(d.feat));
    %if isempty(d.feat)
    %    continue;
    %end

    feat_len = size(d.feat, 2);
    if feat_len == 9216
        featFC8 = rcnn_pool5_to_fcX(d.feat, 8, rcnn_model);
        featFC8 = rcnn_softmax(featFC8);
    else
        featFC8 = d.feat;
    end
    featFC8_score = featFC8(:,2:end);

    overlap = d.overlap;

    boxes = d.boxes;

    num_boxes = size(boxes, 1);
    for j = 1:num_boxes
        [ov, gt_label] = max(overlap(j,:));
        % zero overlap => label = 0 (background)
        if gt_label == skip_class
            continue;
        end
        if ov < 0.5
          gt_label = 0;
        end
        [~, pred_label] = max(featFC8(j,:));
        pred_label = pred_label-1;

        if gt_label == 0
            gt_un_label = 0;
        else
            gt_un_label =  find(cumul_num_clusters >= gt_label, 1);
        end
        if pred_label == 0
            pred_un_label = 0;
        else
            pred_un_label =  find(cumul_num_clusters >= pred_label, 1);
        end

        total_un_samples(gt_un_label+1) = total_un_samples(gt_un_label+1) + 1;
        total_samples(gt_label+1) = total_samples(gt_label+1) + 1;

        if pred_label == gt_label
            correct_samples(gt_label+1) = correct_samples(gt_label+1) + 1;
        end
        if pred_un_label == gt_un_label
            correct_un_samples(gt_un_label+1) = correct_un_samples(gt_un_label+1) + 1;
        end

        [~, max_labels] = sort(featFC8(j,:), 'descend');
        for k = 1:top_k
            pred_k_label = max_labels(k)-1;
            if pred_k_label == gt_label
                correct_topk_samples(gt_label+1) = correct_topk_samples(gt_label+1) + 1;
                break;
            end
        end

    end
end

toc

fprintf('accuracy for positives idpr %f\n', sum(correct_samples(2:end))/sum(total_samples(2:end)));
fprintf('accuracy for positives top k idpr %f\n', sum(correct_topk_samples(2:end))/sum(total_samples(2:end)));
fprintf('accuracy for positives unary %f\n', sum(correct_un_samples(2:end))/sum(total_un_samples(2:end)));

fprintf('done\n');

if (isdeployed)
    close all;
end

end

author: eldar

repo: deepcut

save_path: github-repos/MATLAB/eldar-deepcut


text:
%% STUDY_04_wrist_scaffold_multiple_fiber_evaluation
% Below is a demonstration for:
%
% * Finite element analysis of the performance of additively manufactured scaffolds for scapholunate ligament reconstruction
%___________________________________________________________________________________________________________________________
%  This study presents a patient-specific computational biomechanical evaluation of the effect of scaffold length,
%  and positioning of the bone attachment sites. Through segmentation and image processing of medical image data
%  for natural wrist motion, detailed 3D geometries as well as patient-specific physiological wrist motion could
%  be derived. This data formed the input for detailed finite element analysis, enabling computational of scaffold
%  stress and strain distributions, which are key predictors of scaffold structural integrity.
%___________________________________________________________________________________________________________________________
% * febio_spec version 3.0
% * febio, FEBio
% * hexahedral elements, hex8
% * static, solid
% * hyperelastic, Ogden
% * displacement logfile
% * stress logfile

clear; close all; clc;

% Plot settings
fontSize=40;
numSmoothStepsMain=1;
cParSmoothMain.n=numSmoothStepsMain;
cParSmoothMain.Method='HC';

%%
% Path names
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'2021_Nataliya_Perevoshchikova_wrist_project','data','temp');
loadPathSurface=fullfile(defaultFolder,'2021_Nataliya_Perevoshchikova_wrist_project','data','STL');
loadPathMotion=fullfile(defaultFolder,'2021_Nataliya_Perevoshchikova_wrist_project','data','motion');

%% Building a quadrilateral circular mesh
d=0.35;
h=7;
xSpacing=0.93;
ySpacing=0.525;
nCopies_x=4;
nCopies_y=7;
r=d/2;
ne=2; %Elements in radius
f=0.6; %Fraction (with respect to outer radius) where central square appears

pointSpacingHeight=0.1;
numStepsSweep=ceil(h/pointSpacingHeight);

xc=zeros(numStepsSweep,1);
yc=zeros(numStepsSweep,1);
zc=linspace(-h/2,h/2,numStepsSweep)';
Vc=[xc yc zc];

%SweepLoft
%Create the mesh
[Fq,Vq]=discQuadMesh(ne,r,f);

Vq(:,3)=0;
Vq1=Vq;
Vq1(:,3)=Vq1(:,3)-h/2;
Vq2=Vq;
Vq2(:,3)=Vq2(:,3)+h/2;

% Visualizing mesh
cFigure; hold on;
title('Single fiber','FontSize',fontSize);
plotV(Vq1,'k.','lineWidth',5,'MarkerSize',25);
plotV(Vq2,'r.','lineWidth',5,'MarkerSize',25);
plotV(Vc,'g.-','lineWidth',5,'MarkerSize',25);

axisGeom;
camlight headlight;
drawnow;

[~,~,~,S]=sweepLoft(Vq1,Vq2,[0 0 1],[0 0 1],Vc,numStepsSweep,0,0);

X=S.X'; Y=S.Y'; Z=S.Z'; %Coordinate matrices
V=[X(:) Y(:) Z(:)]; %Create node list

I=size(Vq,1)*((1:1:numStepsSweep)-1);
I=I(ones(size(Fq,1),1),:);
I=I(:);

FQ=repmat(Fq,numStepsSweep,1)+I(:,ones(size(Fq,2),1));
Efib=[FQ(1:end-size(Fq,1),:) FQ(size(Fq,1)+1:end,:)]; %The hexahedral elements

[Efib,V,ind1,ind2]=mergeVertices(Efib,V); %Merge nodes (start and end are not shared yet)
%%
nCopies=nCopies_x*nCopies_y;
VC=repmat({V},nCopies,1);
EC=repmat({Efib},nCopies,1);
numNodesBar=size(V,1);
[Efib,V]=joinElementSets(EC,VC);

c=1;

for qx=1:1:nCopies_x
    for qy=1:1:nCopies_y
        if c>1
            ind1=((c-1)*numNodesBar)+1;
            ind2=(c*numNodesBar);
            V(ind1:ind2,1)=V(ind1:ind2,1)+(qx-1)*xSpacing;
            V(ind1:ind2,2)=V(ind1:ind2,2)+(qy-1)*ySpacing;
        end
        c=c+1;
    end
end

[F]=element2patch(Efib);

%% Get boundary conditions faces
%Get boundary faces
ind=tesBoundary(F,V);
Fb=F(ind,:);
%Get tops/bottoms
Nb=patchNormal(Fb,V);

zVec=[0 0 1];
d=dot(Nb,zVec(ones(size(Nb,1),1),:),2);
Z=V(:,3);
ZF=mean(Z(Fb),2);
logicTop_Fb=(d>0.9) & ZF>=(max(V(:,3))-eps(1));
logicBottom_Fb=(d<-0.9) & ZF<=(min(V(:,3))+eps(1));
F_start=Fb(logicTop_Fb,:);
F_end=Fb(logicBottom_Fb,:);

logic_Lig=(ZF>-1.5) & ZF<=1.5;
F_Lig=Fb(logic_Lig,:);

%%
% Visualizing mesh
cFigure; hold on;
title('Multiple fibers','FontSize',fontSize);
patch('faces',Fb,'vertices',V,'FaceColor','g','FaceAlpha',0.3,'EdgeColor','None');
plotV(V(F_end',:),'k.','lineWidth',5,'MarkerSize',20);
plotV(V(F_start',:),'r.','lineWidth',5,'MarkerSize',20);

camlight headlight;
axisGeom(gca,fontSize);
drawnow;

F_L=F_end;
F_S=F_start;

%% Scaphoid and lunate at first frame from TriPlane, Brown University
loadName=fullfile(loadPathSurface,['modelGeometry.mat']);
dataStruct=load(loadName);

FS_triP=dataStruct.FS;
FL_triP=dataStruct.FL;
FC_triP=dataStruct.FC;
FR_triP=dataStruct.FR;

VL_triP=dataStruct.VL;
VS_triP=dataStruct.VS;
VC_triP=dataStruct.VC;
VR_triP=dataStruct.VR;

%%%
hFig=cFigure;
hold on;
title('Radiography','FontSize',fontSize);

gpatch(FS_triP,VS_triP,'w','none',0.5);
gpatch(FL_triP,VL_triP,'w','none',0.5);
gpatch(FC_triP,VC_triP,'w','none',0.5);
gpatch(FR_triP,VR_triP,'w','none',0.5);

grid on; axis equal;
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Cutting Radius
pointSpacingR=mean(patchEdgeLengths(FR_triP,VR_triP));
%% Cut Capitate bone
% The capitate is cut in the x direction.

%Create a logic for cutting away faces
max_X1=max(VR_triP(:,1))-70*pointSpacingR; %Max x-level used for cutting
logicVertices=VR_triP(:,1)<max_X1; %Logic for the points below this level
logicFaces=all(logicVertices(FR_triP),2); %Logic for the faces
logicFaces=triSurfLogicSharpFix(FR_triP,logicFaces,3); %Altered logic so it is smoother

%%
% Visualize
cFigure; hold on;
gpatch(FS_triP,VS_triP,'w','none',0.5);
gpatch(FL_triP,VL_triP,'w','none',0.5);
gpatch(FR_triP,VR_triP,logicFaces,'none',0.5);
gpatch(FC_triP,VC_triP,'w','none',0.5);

camlight('headlight');
axisGeom(gca,fontSize);
drawnow;

%%
% Cut away faces using logic
FR_triP=FR_triP(logicFaces,:); %The faces to keep
[FR_triP,VR_triP]=patchCleanUnused(FR_triP,VR_triP); %Remove unused points

%Attempt to self triangulate potentially jagged edge
Eb=patchBoundary(FR_triP,VR_triP); %Get boundary edges
indBoundary=edgeListToCurve(Eb); %Convert boundary edges to a curve list
indBoundary=indBoundary(1:end-1); %Trim off last point since it is equal to first on a closed loop
angleThreshold=pi*(120/180); %threshold for self triangulation
[FR_triP,VR_triP,indBoundaryTop]=triSurfSelfTriangulateBoundary(FR_triP,VR_triP,indBoundary,angleThreshold,1);

%Force boundary to have the max X level chosen
VR_triP(indBoundaryTop,1)=max_X1;

%%
% Visualize
cFigure; hold on;
gpatch(FS_triP,VS_triP,'w','none',0.5);
gpatch(FL_triP,VL_triP,'w','none',0.5);
gpatch(FC_triP,VC_triP,'w','none',0.5);
gpatch(FR_triP,VR_triP,'w','k',0.5);
plotV(VR_triP(indBoundaryTop,:),'r-','LineWidth',3);

camlight('headlight');
axisGeom(gca,fontSize);
drawnow;

%% Close over top of radius
% The top boundary curve of the cut surface is filled with triangles. This
% is a 2D method. The x-coordinate is added after.
[F2t,V2t]=regionTriMesh2D({VR_triP(indBoundaryTop,[2 3])},pointSpacingR,0,0);
V2t(:,3)=V2t(:,2);
V2t(:,2)=V2t(:,1);
V2t(:,1)=max_X1; %Add/set X-level

%%
% Visualize
cFigure; hold on;
gpatch(FS_triP,VS_triP,'w','none',0.5);
gpatch(FL_triP,VL_triP,'w','none',0.5);
gpatch(FC_triP,VC_triP,'w','none',0.5);
gpatch(FR_triP,VR_triP,'rw','k',0.5);
gpatch(F2t,V2t,'gw','k',1);
plotV(VR_triP(indBoundaryTop,:),'r-','LineWidth',3);
camlight('headlight');
axisGeom(gca,fontSize);
drawnow;

%% Joining surface features
% Add all surface sets together in joint list of faces, vertices
[FR_triP,VR_triP]=joinElementSets({FR_triP,F2t},{VR_triP,V2t});
%% Merge shared nodes
% The join operation only adds the sets together. Nodes with the same
% coordinates are not seen as the same yet and need to be merged.
[FR_triP,VR_triP]=mergeVertices(FR_triP,VR_triP); %Merge nodes

%Remove tri connected
[FS,VS]=triSurfRemoveThreeConnect(FS_triP,VS_triP);
[FL,VL]=triSurfRemoveThreeConnect(FL_triP,VL_triP);
[FC,VC]=triSurfRemoveThreeConnect(FC_triP,VC_triP);
[FR,VR]=triSurfRemoveThreeConnect(FR_triP,VR_triP);
%Smoothen
[VL]=patchSmooth(FL,VL,[],cParSmoothMain);
[VS]=patchSmooth(FS,VS,[],cParSmoothMain);
[VR]=patchSmooth(FR,VR,[],cParSmoothMain);
[VC]=patchSmooth(FC,VC,[],cParSmoothMain);
%Refine
[FS,VS]=subtri(FS,VS,1);
[FL,VL]=subtri(FL,VL,1);
[FR,VR]=subtri(FR,VR,1);
[FC,VC]=subtri(FC,VC,1);
%Smoothen
[VL]=patchSmooth(FL,VL,[],cParSmoothMain);
[VS]=patchSmooth(FS,VS,[],cParSmoothMain);
[VR]=patchSmooth(FR,VR,[],cParSmoothMain);
[VC]=patchSmooth(FC,VC,[],cParSmoothMain);
%%%%

%Construct placement according to natural SLIL positioning
alphaRot=pi/2;
R=euler2DCM([alphaRot 0 0]); %The rotation tensor for each step
V=V*R; %Rotated 1 part for visualization of stepwise amount

alphaRot=pi/2;
R=euler2DCM([0 alphaRot 0]); %The rotation tensor for each step
V=V*R; %Rotated 1 part for visualization of stepwise amount

alphaRot=-pi/15;
R=euler2DCM([0 0 alphaRot]); %The rotation tensor for each step
V=V*R; %Rotated 1 part for visualization of stepwise amount

V(:,3)=V(:,3)-1.5;
V(:,1)=V(:,1)-15;
V(:,2)=V(:,2)-3.2;

hFig=cFigure;
hold on;

gpatch(FS,VS,'w','none',0.7);
gpatch(FL,VL,'w','none',0.7);

gpatch(FC,VC,'w','none',0.5);
gpatch(FR,VR,'w','none',0.5);

gpatch(Fb,V,'gw','k',1);

grid on; axis equal;
axisGeom(gca,fontSize);
camlight headlight;
drawnow;

%% Joining node sets
[E_F,V,C]=joinElementSets({Efib,FS,FL,FC,FR},{V,VS,VL,VC,VR});
E=E_F{1};
FS1=E_F{2};
FL1=E_F{3};
FC1=E_F{4};
FR1=E_F{5};

% hFig=cFigure;
% hold on;
%
% gpatch(FS1,V,'w','none',0.7);
% gpatch(FL1,V,'w','none',0.7);
%
% gpatch(FC1,V,'w','none',0.5);
% gpatch(FR1,V,'w','none',0.5);
%
% gpatch(E,V,'gw','k',1);
%
% grid on; axis equal;
% axisGeom(gca,fontSize);
% camlight headlight;
% drawnow;

%% FEA control settings
Step=5;
numTimeSteps=round((1989-1000)/Step);%Number of time steps desired
max_refs=300; %Max reforms
max_ups=0; %Set to zero to use full-Newton iterations
opt_iter=150; %Optimum number of iterations
max_retries=50; %Maximum number of retires
dtmin=(1/numTimeSteps)/10; %Minimum time step size
dtmax=1/numTimeSteps; %Maximum time step size
tStep=0.005*Step; %Duration of step
min_residual=1e-20;
Mag_val=1;

%% Control parameters
% Defining file names
febioFebFileNamePart='tempModel';
febioFebFileName=fullfile(savePath,[febioFebFileNamePart,'.feb']); %FEB file name
febioLogFileName=[febioFebFileNamePart,'.txt']; %FEBio log file name
febioLogFileName_disp=[febioFebFileNamePart,'_disp_out.txt']; %Log file name for exporting displacement
febioLogFileName_force=[febioFebFileNamePart,'_force_out.txt']; %Log file name for exporting force
febioLogFileName_stress=[febioFebFileNamePart,'_stress_out.txt']; %Log file name for exporting stress

%% Defining the FEBio input structure
% See also |febioStructTemplate| and |febioStruct2xml| and the FEBio user
% manual.

%Get a template with default settings
[febio_spec]=febioStructTemplate;

%febio_spec version
febio_spec.ATTR.version='3.0';

%Module section
febio_spec.Module.ATTR.type='solid';

%Control section
febio_spec.Control.analysis='STATIC';
febio_spec.Control.time_steps=numTimeSteps;
febio_spec.Control.step_size=tStep;
febio_spec.Control.time_stepper.dtmin=dtmin;
febio_spec.Control.time_stepper.dtmax=dtmax;
febio_spec.Control.time_stepper.max_retries=max_retries;
febio_spec.Control.time_stepper.opt_iter=opt_iter;
febio_spec.Control.solver.max_refs=max_refs;
febio_spec.Control.solver.max_ups=max_ups;
febio_spec.Control.solver.min_residual=min_residual;

%Material section
materialName1='Material1';
febio_spec.Material.material{1}.ATTR.id=1;
febio_spec.Material.material{1}.ATTR.name=materialName1;

%Material section
C1=57.21894783089972;
C2=57.23071414361289;
K=11444.966197451262;
M1=2.00012552045284;
M2=-1.9997730268725347;
t1=25.365;
g1=1.23339;

materialOption='elastic';%'elastic_viscoelastic'
switch materialOption
    case 'elastic'
        febio_spec.Material.material{1}.ATTR.type='Ogden';
        febio_spec.Material.material{1}.c1=C1;
        febio_spec.Material.material{1}.m1=M1;
        febio_spec.Material.material{1}.c2=C2;
        febio_spec.Material.material{1}.m2=M2;
        febio_spec.Material.material{1}.k=K;
    case 'elastic_viscoelastic'
        %Viscoelastic part
        febio_spec.Material.material{1}.ATTR.type='uncoupled viscoelastic';
        febio_spec.Material.material{1}.g1=g1;
        febio_spec.Material.material{1}.t1=t1;

        %Elastic part
        febio_spec.Material.material{1}.elastic{1}.ATTR.type='Ogden';
        febio_spec.Material.material{1}.elastic{1}.c1=C1;
        febio_spec.Material.material{1}.elastic{1}.m1=M1;
        febio_spec.Material.material{1}.elastic{1}.c2=C2;
        febio_spec.Material.material{1}.elastic{1}.m2=M2;
        febio_spec.Material.material{1}.elastic{1}.k=K;

end

materialName2='Lunate_side';
febio_spec.Material.material{2}.ATTR.name=materialName2;
febio_spec.Material.material{2}.ATTR.type='rigid body';
febio_spec.Material.material{2}.ATTR.id=2;
febio_spec.Material.material{2}.density=1;
febio_spec.Material.material{2}.center_of_mass=[0 0 0];

materialName3='Scaphoid_side';
febio_spec.Material.material{3}.ATTR.name=materialName3;
febio_spec.Material.material{3}.ATTR.type='rigid body';
febio_spec.Material.material{3}.ATTR.id=3;
febio_spec.Material.material{3}.density=1;
febio_spec.Material.material{3}.center_of_mass=[0 0 0];

materialName4='Lunate';
febio_spec.Material.material{4}.ATTR.name=materialName4;
febio_spec.Material.material{4}.ATTR.type='rigid body';
febio_spec.Material.material{4}.ATTR.id=4;
febio_spec.Material.material{4}.density=1;
febio_spec.Material.material{4}.center_of_mass=[0 0 0];

materialName5='Scaphoid';
febio_spec.Material.material{5}.ATTR.name=materialName5;
febio_spec.Material.material{5}.ATTR.type='rigid body';
febio_spec.Material.material{5}.ATTR.id=5;
febio_spec.Material.material{5}.density=1;
febio_spec.Material.material{5}.center_of_mass=[0 0 0];

materialName6='Capitate';
febio_spec.Material.material{6}.ATTR.name=materialName6;
febio_spec.Material.material{6}.ATTR.type='rigid body';
febio_spec.Material.material{6}.ATTR.id=6;
febio_spec.Material.material{6}.density=1;
febio_spec.Material.material{6}.center_of_mass=[0 0 0];

materialName7='Radius';
febio_spec.Material.material{7}.ATTR.name=materialName7;
febio_spec.Material.material{7}.ATTR.type='rigid body';
febio_spec.Material.material{7}.ATTR.id=7;
febio_spec.Material.material{7}.density=1;
febio_spec.Material.material{7}.center_of_mass=[0 0 0];

%Geometry section
% -> Nodes
febio_spec.Mesh.Nodes{1}.ATTR.name='nodeSet_all'; %The node set name
febio_spec.Mesh.Nodes{1}.node.ATTR.id=(1:size(V,1))'; %The node id's
febio_spec.Mesh.Nodes{1}.node.VAL=V; %The nodel coordinates

% % -> Elements
partName1='Scaffold';
febio_spec.Mesh.Elements{1}.ATTR.name=partName1; %Name of the element set
febio_spec.Mesh.Elements{1}.ATTR.mat=1; %material index for this set
febio_spec.Mesh.Elements{1}.ATTR.type='hex8'; %Element type of this set
febio_spec.Mesh.Elements{1}.elem.ATTR.id=(1:1:size(E,1))'; %Element id's
febio_spec.Mesh.Elements{1}.elem.VAL=E;

partName2='Lunate side';
febio_spec.Mesh.Elements{2}.ATTR.name=partName2; %Name of the element set
febio_spec.Mesh.Elements{2}.ATTR.mat=2; %material index for this set
febio_spec.Mesh.Elements{2}.ATTR.type='quad4'; %Element type of this set
febio_spec.Mesh.Elements{2}.elem.ATTR.id=size(E,1)+(1:1:size(F_L,1))'; %Element id's
febio_spec.Mesh.Elements{2}.elem.VAL=F_L;

partName3='Scaphoid side';
febio_spec.Mesh.Elements{3}.ATTR.name=partName3; %Name of the element set
febio_spec.Mesh.Elements{3}.ATTR.mat=3; %material index for this set
febio_spec.Mesh.Elements{3}.ATTR.type='quad4'; %Element type of this set
febio_spec.Mesh.Elements{3}.elem.ATTR.id=size(E,1)+(1:1:size(F_S,1))'; %Element id's
febio_spec.Mesh.Elements{3}.elem.VAL=F_S;

partName4='Lunate';
febio_spec.Mesh.Elements{4}.ATTR.name=partName4; %Name of the element set
febio_spec.Mesh.Elements{4}.ATTR.mat=4; %material index for this set
febio_spec.Mesh.Elements{4}.ATTR.type='tri3'; %Element type of this set
febio_spec.Mesh.Elements{4}.elem.ATTR.id=size(E,1)+(1:1:size(FL1,1))'; %Element id's
febio_spec.Mesh.Elements{4}.elem.VAL=FL1;

partName5='Scaphoid';
febio_spec.Mesh.Elements{5}.ATTR.name=partName5; %Name of the element set
febio_spec.Mesh.Elements{5}.ATTR.mat=5; %material index for this set
febio_spec.Mesh.Elements{5}.ATTR.type='tri3'; %Element type of this set
febio_spec.Mesh.Elements{5}.elem.ATTR.id=size(E,1)+(1:1:size(FS1,1))'; %Element id's
febio_spec.Mesh.Elements{5}.elem.VAL=FS1;

%Capitate
partName6='Capitate';
febio_spec.Mesh.Elements{6}.ATTR.name=partName6; %Name of the element set
febio_spec.Mesh.Elements{6}.ATTR.mat=6; %material index for this set
febio_spec.Mesh.Elements{6}.ATTR.type='tri3'; %Element type of this set
febio_spec.Mesh.Elements{6}.elem.ATTR.id=size(E,1)+(1:1:size(FC1,1))'; %Element id's
febio_spec.Mesh.Elements{6}.elem.VAL=FC1;

%Radius
partName7='Radius';
febio_spec.Mesh.Elements{7}.ATTR.name=partName7; %Name of the element set
febio_spec.Mesh.Elements{7}.ATTR.mat=7; %material index for this set
febio_spec.Mesh.Elements{7}.ATTR.type='tri3'; %Element type of this set
febio_spec.Mesh.Elements{7}.elem.ATTR.id=size(E,1)+(1:1:size(FR1,1))'; %Element id's
febio_spec.Mesh.Elements{7}.elem.VAL=FR1;

%MeshDomains section
febio_spec.MeshDomains.SolidDomain.ATTR.name=partName1;
febio_spec.MeshDomains.SolidDomain.ATTR.mat=materialName1;

febio_spec.MeshDomains.ShellDomain{1}.ATTR.name=partName2;
febio_spec.MeshDomains.ShellDomain{1}.ATTR.mat=materialName2;

febio_spec.MeshDomains.ShellDomain{2}.ATTR.name=partName3;
febio_spec.MeshDomains.ShellDomain{2}.ATTR.mat=materialName3;

febio_spec.MeshDomains.ShellDomain{3}.ATTR.name=partName4;
febio_spec.MeshDomains.ShellDomain{3}.ATTR.mat=materialName4;

febio_spec.MeshDomains.ShellDomain{4}.ATTR.name=partName5;
febio_spec.MeshDomains.ShellDomain{4}.ATTR.mat=materialName5;

febio_spec.MeshDomains.ShellDomain{5}.ATTR.name=partName6;
febio_spec.MeshDomains.ShellDomain{5}.ATTR.mat=materialName6;

febio_spec.MeshDomains.ShellDomain{6}.ATTR.name=partName7;
febio_spec.MeshDomains.ShellDomain{6}.ATTR.mat=materialName7;

% Boundary conditions
i=1;
for count=1:numTimeSteps

    FileData = load(fullfile(loadPathMotion, sprintf('Flex_Ext_%d_TriPlane_Lunate.mat',count)));
    MatT_Lunate = FileData.Trans_L;

    FileData = load(fullfile(loadPathMotion, sprintf('Flex_Ext_%d_TriPlane_Scaphoid.mat',count)));
    MatT_Scaphoid = FileData.Trans_S;

    FileData = load(fullfile(loadPathMotion, sprintf('Flex_Ext_%d_TriPlane_Capitate.mat',count)));
    MatT_Capitate = FileData.Trans_C;

    RLunate                       = MatT_Lunate(1:3,1:3);%radians
    RScaphoid                     = MatT_Scaphoid(1:3,1:3); %radians
    RCapitate                     = MatT_Capitate(1:3,1:3); %radians

    Prescribe_Lunate{i}       = MatT_Lunate(1:3,end);
    Prescribe_Scaphoid{i}     = MatT_Scaphoid(1:3,end);
    Prescribe_Capitate{i}     = MatT_Capitate(1:3,end);

    axangL                    = rotm2axang(RLunate);%Axis-Angle {[x, y, z], angle (radians)}
    AngleLunate{i}            = axangL(1:3)*axangL(4);%angle *rot_axis; Axis with angle magnitude (radians) [x, y, z]
    axangS                    = rotm2axang(RScaphoid);
    AngleScaphoid{i}          = axangS(1:3)*axangS(4);%angle*rot_axis;
    axangC                    = rotm2axang(RCapitate);
    AngleCapitate{i}          = axangC(1:3)*axangC(4);%angle*rot_axis;
    tStartNow(i)=count*tStep;
    i=i+1;
end

%Boundary conditions
%Rigid section
% %Start and End
febio_spec.Rigid.rigid_constraint{1}.ATTR.name='Rigid_Lunate_X';
febio_spec.Rigid.rigid_constraint{1}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{1}.rb=2;
febio_spec.Rigid.rigid_constraint{1}.dof='Rx';
febio_spec.Rigid.rigid_constraint{1}.value.ATTR.lc=1;
febio_spec.Rigid.rigid_constraint{1}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{1}.relative=0;

febio_spec.Rigid.rigid_constraint{2}.ATTR.name='Rigid_Lunate_Y';
febio_spec.Rigid.rigid_constraint{2}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{2}.rb=2;
febio_spec.Rigid.rigid_constraint{2}.dof='Ry';
febio_spec.Rigid.rigid_constraint{2}.value.ATTR.lc=2;
febio_spec.Rigid.rigid_constraint{2}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{2}.relative=0;

febio_spec.Rigid.rigid_constraint{3}.ATTR.name='Rigid_Lunate_Z';
febio_spec.Rigid.rigid_constraint{3}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{3}.rb=2;
febio_spec.Rigid.rigid_constraint{3}.dof='Rz';
febio_spec.Rigid.rigid_constraint{3}.value.ATTR.lc=3;
febio_spec.Rigid.rigid_constraint{3}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{3}.relative=0;

febio_spec.Rigid.rigid_constraint{4}.ATTR.name='Rigid_Lunate_Ru';
febio_spec.Rigid.rigid_constraint{4}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{4}.rb=2;
febio_spec.Rigid.rigid_constraint{4}.dof='Ru';
febio_spec.Rigid.rigid_constraint{4}.value.ATTR.lc=4;
febio_spec.Rigid.rigid_constraint{4}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{4}.relative=0;

febio_spec.Rigid.rigid_constraint{5}.ATTR.name='Rigid_Lunate_Rv';
febio_spec.Rigid.rigid_constraint{5}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{5}.rb=2;
febio_spec.Rigid.rigid_constraint{5}.dof='Rv';
febio_spec.Rigid.rigid_constraint{5}.value.ATTR.lc=5;
febio_spec.Rigid.rigid_constraint{5}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{5}.relative=0;

febio_spec.Rigid.rigid_constraint{6}.ATTR.name='Rigid_Lunate_Rw';
febio_spec.Rigid.rigid_constraint{6}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{6}.rb=2;
febio_spec.Rigid.rigid_constraint{6}.dof='Rw';
febio_spec.Rigid.rigid_constraint{6}.value.ATTR.lc=6;
febio_spec.Rigid.rigid_constraint{6}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{6}.relative=0;

febio_spec.Rigid.rigid_constraint{7}.ATTR.name='Rigid_Scaphoid_X';
febio_spec.Rigid.rigid_constraint{7}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{7}.rb=3;
febio_spec.Rigid.rigid_constraint{7}.dof='Rx';
febio_spec.Rigid.rigid_constraint{7}.value.ATTR.lc=7;
febio_spec.Rigid.rigid_constraint{7}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{7}.relative=0;

febio_spec.Rigid.rigid_constraint{8}.ATTR.name='Rigid_Scaphoid_Y';
febio_spec.Rigid.rigid_constraint{8}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{8}.rb=3;
febio_spec.Rigid.rigid_constraint{8}.dof='Ry';
febio_spec.Rigid.rigid_constraint{8}.value.ATTR.lc=8;
febio_spec.Rigid.rigid_constraint{8}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{8}.relative=0;

febio_spec.Rigid.rigid_constraint{9}.ATTR.name='Rigid_Scaphoid_Z';
febio_spec.Rigid.rigid_constraint{9}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{9}.rb=3;
febio_spec.Rigid.rigid_constraint{9}.dof='Rz';
febio_spec.Rigid.rigid_constraint{9}.value.ATTR.lc=9;
febio_spec.Rigid.rigid_constraint{9}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{9}.relative=0;

febio_spec.Rigid.rigid_constraint{10}.ATTR.name='Rigid_Scaphoid_RX';
febio_spec.Rigid.rigid_constraint{10}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{10}.rb=3;
febio_spec.Rigid.rigid_constraint{10}.dof='Ru';
febio_spec.Rigid.rigid_constraint{10}.value.ATTR.lc=10;
febio_spec.Rigid.rigid_constraint{10}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{10}.relative=0;

febio_spec.Rigid.rigid_constraint{11}.ATTR.name='Rigid_Scaphoid_RY';
febio_spec.Rigid.rigid_constraint{11}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{11}.rb=3;
febio_spec.Rigid.rigid_constraint{11}.dof='Rv';
febio_spec.Rigid.rigid_constraint{11}.value.ATTR.lc=11;
febio_spec.Rigid.rigid_constraint{11}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{11}.relative=0;

febio_spec.Rigid.rigid_constraint{12}.ATTR.name='Rigid_Scaphoid_RZ';
febio_spec.Rigid.rigid_constraint{12}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{12}.rb=3;
febio_spec.Rigid.rigid_constraint{12}.dof='Rw';
febio_spec.Rigid.rigid_constraint{12}.value.ATTR.lc=12;
febio_spec.Rigid.rigid_constraint{12}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{12}.relative=0;

% % -> Prescribed boundary conditions on the rigid body bones
febio_spec.Rigid.rigid_constraint{13}.ATTR.name='Rigid_Lunate_side_X';
febio_spec.Rigid.rigid_constraint{13}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{13}.rb=4;
febio_spec.Rigid.rigid_constraint{13}.dof='Rx';
febio_spec.Rigid.rigid_constraint{13}.value.ATTR.lc=1;
febio_spec.Rigid.rigid_constraint{13}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{13}.relative=0;

febio_spec.Rigid.rigid_constraint{14}.ATTR.name='Rigid_Lunate_side_Y';
febio_spec.Rigid.rigid_constraint{14}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{14}.rb=4;
febio_spec.Rigid.rigid_constraint{14}.dof='Ry';
febio_spec.Rigid.rigid_constraint{14}.value.ATTR.lc=2;
febio_spec.Rigid.rigid_constraint{14}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{14}.relative=0;

febio_spec.Rigid.rigid_constraint{15}.ATTR.name='Rigid_Lunate_side_Z';
febio_spec.Rigid.rigid_constraint{15}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{15}.rb=4;
febio_spec.Rigid.rigid_constraint{15}.dof='Rz';
febio_spec.Rigid.rigid_constraint{15}.value.ATTR.lc=3;
febio_spec.Rigid.rigid_constraint{15}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{15}.relative=0;

febio_spec.Rigid.rigid_constraint{16}.ATTR.name='Rigid_Lunate_side_RX';
febio_spec.Rigid.rigid_constraint{16}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{16}.rb=4;
febio_spec.Rigid.rigid_constraint{16}.dof='Ru';
febio_spec.Rigid.rigid_constraint{16}.value.ATTR.lc=4;
febio_spec.Rigid.rigid_constraint{16}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{16}.relative=0;

febio_spec.Rigid.rigid_constraint{17}.ATTR.name='Rigid_Lunate_side_RY';
febio_spec.Rigid.rigid_constraint{17}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{17}.rb=4;
febio_spec.Rigid.rigid_constraint{17}.dof='Rv';
febio_spec.Rigid.rigid_constraint{17}.value.ATTR.lc=5;
febio_spec.Rigid.rigid_constraint{17}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{17}.relative=0;

febio_spec.Rigid.rigid_constraint{18}.ATTR.name='Rigid_Lunate_side_RZ';
febio_spec.Rigid.rigid_constraint{18}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{18}.rb=4;
febio_spec.Rigid.rigid_constraint{18}.dof='Rw';
febio_spec.Rigid.rigid_constraint{18}.value.ATTR.lc=6;
febio_spec.Rigid.rigid_constraint{18}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{18}.relative=0;

% %Scaphoid
febio_spec.Rigid.rigid_constraint{19}.ATTR.name='Rigid_Scaphoid_side_X';
febio_spec.Rigid.rigid_constraint{19}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{19}.rb=5;
febio_spec.Rigid.rigid_constraint{19}.dof='Rx';
febio_spec.Rigid.rigid_constraint{19}.value.ATTR.lc=7;
febio_spec.Rigid.rigid_constraint{19}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{19}.relative=0;

febio_spec.Rigid.rigid_constraint{20}.ATTR.name='Rigid_Scaphoid_side_Y';
febio_spec.Rigid.rigid_constraint{20}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{20}.rb=5;
febio_spec.Rigid.rigid_constraint{20}.dof='Ry';
febio_spec.Rigid.rigid_constraint{20}.value.ATTR.lc=8;
febio_spec.Rigid.rigid_constraint{20}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{20}.relative=0;

febio_spec.Rigid.rigid_constraint{21}.ATTR.name='Rigid_Scaphoid_side_Z';
febio_spec.Rigid.rigid_constraint{21}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{21}.rb=5;
febio_spec.Rigid.rigid_constraint{21}.dof='Rz';
febio_spec.Rigid.rigid_constraint{21}.value.ATTR.lc=9;
febio_spec.Rigid.rigid_constraint{21}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{21}.relative=0;

febio_spec.Rigid.rigid_constraint{22}.ATTR.name='Rigid_Scaphoid_side_RX';
febio_spec.Rigid.rigid_constraint{22}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{22}.rb=5;
febio_spec.Rigid.rigid_constraint{22}.dof='Ru';
febio_spec.Rigid.rigid_constraint{22}.value.ATTR.lc=10;
febio_spec.Rigid.rigid_constraint{22}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{22}.relative=0;

febio_spec.Rigid.rigid_constraint{23}.ATTR.name='Rigid_Scaphoid_side_RY';
febio_spec.Rigid.rigid_constraint{23}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{23}.rb=5;
febio_spec.Rigid.rigid_constraint{23}.dof='Rv';
febio_spec.Rigid.rigid_constraint{23}.value.ATTR.lc=11;
febio_spec.Rigid.rigid_constraint{23}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{23}.relative=0;

febio_spec.Rigid.rigid_constraint{24}.ATTR.name='Rigid_Scaphoid_side_RZ';
febio_spec.Rigid.rigid_constraint{24}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{24}.rb=5;
febio_spec.Rigid.rigid_constraint{24}.dof='Rw';
febio_spec.Rigid.rigid_constraint{24}.value.ATTR.lc=12;
febio_spec.Rigid.rigid_constraint{24}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{24}.relative=0;

%%Capitate
febio_spec.Rigid.rigid_constraint{25}.ATTR.name='Rigid_Capitate_X';
febio_spec.Rigid.rigid_constraint{25}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{25}.rb=6;
febio_spec.Rigid.rigid_constraint{25}.dof='Rx';
febio_spec.Rigid.rigid_constraint{25}.value.ATTR.lc=13;
febio_spec.Rigid.rigid_constraint{25}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{25}.relative=0;

febio_spec.Rigid.rigid_constraint{26}.ATTR.name='Rigid_Capitate_Y';
febio_spec.Rigid.rigid_constraint{26}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{26}.rb=6;
febio_spec.Rigid.rigid_constraint{26}.dof='Ry';
febio_spec.Rigid.rigid_constraint{26}.value.ATTR.lc=14;
febio_spec.Rigid.rigid_constraint{26}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{26}.relative=0;

febio_spec.Rigid.rigid_constraint{27}.ATTR.name='Rigid_Capitate_Z';
febio_spec.Rigid.rigid_constraint{27}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{27}.rb=6;
febio_spec.Rigid.rigid_constraint{27}.dof='Rz';
febio_spec.Rigid.rigid_constraint{27}.value.ATTR.lc=15;
febio_spec.Rigid.rigid_constraint{27}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{27}.relative=0;

febio_spec.Rigid.rigid_constraint{28}.ATTR.name='Rigid_Capitate_RX';
febio_spec.Rigid.rigid_constraint{28}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{28}.rb=6;
febio_spec.Rigid.rigid_constraint{28}.dof='Ru';
febio_spec.Rigid.rigid_constraint{28}.value.ATTR.lc=16;
febio_spec.Rigid.rigid_constraint{28}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{28}.relative=0;

febio_spec.Rigid.rigid_constraint{29}.ATTR.name='Rigid_Capitate_RY';
febio_spec.Rigid.rigid_constraint{29}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{29}.rb=6;
febio_spec.Rigid.rigid_constraint{29}.dof='Rv';
febio_spec.Rigid.rigid_constraint{29}.value.ATTR.lc=17;
febio_spec.Rigid.rigid_constraint{29}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{29}.relative=0;

febio_spec.Rigid.rigid_constraint{30}.ATTR.name='Rigid_Capitate_RZ';
febio_spec.Rigid.rigid_constraint{30}.ATTR.type='prescribe';
febio_spec.Rigid.rigid_constraint{30}.rb=6;
febio_spec.Rigid.rigid_constraint{30}.dof='Rw';
febio_spec.Rigid.rigid_constraint{30}.value.ATTR.lc=18;
febio_spec.Rigid.rigid_constraint{30}.value.VAL=Mag_val;
febio_spec.Rigid.rigid_constraint{30}.relative=0;

% %Radius
febio_spec.Rigid.rigid_constraint{31}.ATTR.name='Rigid_Radiuas';
febio_spec.Rigid.rigid_constraint{31}.ATTR.type='fix';
febio_spec.Rigid.rigid_constraint{31}.rb=7;
febio_spec.Rigid.rigid_constraint{31}.dofs='Rx,Ry,Rz,Ru,Rv,Rw';

k=1;
for i=1:numTimeSteps+1

    if k==1
        Mat_X_Lunate(k,1)=0;
        Mat_X_Lunate(k,2)=0;

        Mat_Y_Lunate(k,1)=0;
        Mat_Y_Lunate(k,2)=0;

        Mat_Z_Lunate(k,1)=0;
        Mat_Z_Lunate(k,2)=0;

        Mat_RX_Lunate(k,1)=0;
        Mat_RX_Lunate(k,2)=0;

        Mat_RY_Lunate(k,1)=0;
        Mat_RY_Lunate(k,2)=0;

        Mat_RZ_Lunate(k,1)=0;
        Mat_RZ_Lunate(k,2)=0;

        Mat_X_Scaphoid(k,1)=0;
        Mat_X_Scaphoid(k,2)=0;

        Mat_Y_Scaphoid(k,1)=0;
        Mat_Y_Scaphoid(k,2)=0;

        Mat_Z_Scaphoid(k,1)=0;
        Mat_Z_Scaphoid(k,2)=0;

        Mat_RX_Scaphoid(k,1)=0;
        Mat_RX_Scaphoid(k,2)=0;

        Mat_RY_Scaphoid(k,1)=0;
        Mat_RY_Scaphoid(k,2)=0;

        Mat_RZ_Scaphoid(k,1)=0;
        Mat_RZ_Scaphoid(k,2)=0;

        Mat_X_Capitate(k,1)=0;
        Mat_X_Capitate(k,2)=0;

        Mat_Y_Capitate(k,1)=0;
        Mat_Y_Capitate(k,2)=0;

        Mat_Z_Capitate(k,1)=0;
        Mat_Z_Capitate(k,2)=0;

        Mat_RX_Capitate(k,1)=0;
        Mat_RX_Capitate(k,2)=0;

        Mat_RY_Capitate(k,1)=0;
        Mat_RY_Capitate(k,2)=0;

        Mat_RZ_Capitate(k,1)=0;
        Mat_RZ_Capitate(k,2)=0;
    else
        Mat_X_Lunate(k,1)=tStartNow(k-1);
        Mat_X_Lunate(k,2)=Prescribe_Lunate{k-1}(1);

        Mat_Y_Lunate(k,1)=tStartNow(k-1);
        Mat_Y_Lunate(k,2)=Prescribe_Lunate{k-1}(2);

        Mat_Z_Lunate(k,1)=tStartNow(k-1);
        Mat_Z_Lunate(k,2)=Prescribe_Lunate{k-1}(3);

        Mat_RX_Lunate(k,1)=tStartNow(k-1);
        Mat_RX_Lunate(k,2)=AngleLunate{k-1}(1);

        Mat_RY_Lunate(k,1)=tStartNow(k-1);
        Mat_RY_Lunate(k,2)=AngleLunate{k-1}(2);

        Mat_RZ_Lunate(k,1)=tStartNow(k-1);
        Mat_RZ_Lunate(k,2)=AngleLunate{k-1}(3);

        Mat_X_Scaphoid(k,1)=tStartNow(k-1);
        Mat_X_Scaphoid(k,2)=Prescribe_Scaphoid{k-1}(1);

        Mat_Y_Scaphoid(k,1)=tStartNow(k-1);
        Mat_Y_Scaphoid(k,2)=Prescribe_Scaphoid{k-1}(2);

        Mat_Z_Scaphoid(k,1)=tStartNow(k-1);
        Mat_Z_Scaphoid(k,2)=Prescribe_Scaphoid{k-1}(3);

        Mat_RX_Scaphoid(k,1)=tStartNow(k-1);
        Mat_RX_Scaphoid(k,2)=AngleScaphoid{k-1}(1);

        Mat_RY_Scaphoid(k,1)=tStartNow(k-1);
        Mat_RY_Scaphoid(k,2)=AngleScaphoid{k-1}(2);

        Mat_RZ_Scaphoid(k,1)=tStartNow(k-1);
        Mat_RZ_Scaphoid(k,2)=AngleScaphoid{k-1}(3);

        Mat_X_Capitate(k,1)=tStartNow(k-1);
        Mat_X_Capitate(k,2)=Prescribe_Capitate{k-1}(1);

        Mat_Y_Capitate(k,1)=tStartNow(k-1);
        Mat_Y_Capitate(k,2)=Prescribe_Capitate{k-1}(2);

        Mat_Z_Capitate(k,1)=tStartNow(k-1);
        Mat_Z_Capitate(k,2)=Prescribe_Capitate{k-1}(3);

        Mat_RX_Capitate(k,1)=tStartNow(k-1);
        Mat_RX_Capitate(k,2)=AngleCapitate{k-1}(1);

        Mat_RY_Capitate(k,1)=tStartNow(k-1);
        Mat_RY_Capitate(k,2)=AngleCapitate{k-1}(2);

        Mat_RZ_Capitate(k,1)=tStartNow(k-1);
        Mat_RZ_Capitate(k,2)=AngleCapitate{k-1}(3);

    end
   k=k+1;
end

%Lunate
febio_spec.LoadData.load_controller{1}.ATTR.id=1;
febio_spec.LoadData.load_controller{1}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{1}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{1}.points.point.VAL=Mat_X_Lunate;

febio_spec.LoadData.load_controller{2}.ATTR.id=2;
febio_spec.LoadData.load_controller{2}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{2}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{2}.points.point.VAL=Mat_Y_Lunate;

febio_spec.LoadData.load_controller{3}.ATTR.id=3;
febio_spec.LoadData.load_controller{3}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{3}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{3}.points.point.VAL=Mat_Z_Lunate;

febio_spec.LoadData.load_controller{4}.ATTR.id=4;
febio_spec.LoadData.load_controller{4}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{4}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{4}.points.point.VAL=Mat_RX_Lunate;

febio_spec.LoadData.load_controller{5}.ATTR.id=5;
febio_spec.LoadData.load_controller{5}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{5}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{5}.points.point.VAL=Mat_RY_Lunate;

febio_spec.LoadData.load_controller{6}.ATTR.id=6;
febio_spec.LoadData.load_controller{6}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{6}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{6}.points.point.VAL=Mat_RZ_Lunate;

% %Scaphoid
febio_spec.LoadData.load_controller{7}.ATTR.id=7;
febio_spec.LoadData.load_controller{7}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{7}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{7}.points.point.VAL=Mat_X_Scaphoid;

febio_spec.LoadData.load_controller{8}.ATTR.id=8;
febio_spec.LoadData.load_controller{8}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{8}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{8}.points.point.VAL=Mat_Y_Scaphoid;

febio_spec.LoadData.load_controller{9}.ATTR.id=9;
febio_spec.LoadData.load_controller{9}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{9}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{9}.points.point.VAL=Mat_Z_Scaphoid;

febio_spec.LoadData.load_controller{10}.ATTR.id=10;
febio_spec.LoadData.load_controller{10}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{10}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{10}.points.point.VAL=Mat_RX_Scaphoid;

febio_spec.LoadData.load_controller{11}.ATTR.id=11;
febio_spec.LoadData.load_controller{11}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{11}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{11}.points.point.VAL=Mat_RY_Scaphoid;

febio_spec.LoadData.load_controller{12}.ATTR.id=12;
febio_spec.LoadData.load_controller{12}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{12}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{12}.points.point.VAL=Mat_RZ_Scaphoid;

% %Capitate
febio_spec.LoadData.load_controller{13}.ATTR.id=13;
febio_spec.LoadData.load_controller{13}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{13}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{13}.points.point.VAL=Mat_X_Capitate;

febio_spec.LoadData.load_controller{14}.ATTR.id=14;
febio_spec.LoadData.load_controller{14}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{14}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{14}.points.point.VAL=Mat_Y_Capitate;

febio_spec.LoadData.load_controller{15}.ATTR.id=15;
febio_spec.LoadData.load_controller{15}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{15}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{15}.points.point.VAL=Mat_Z_Capitate;

febio_spec.LoadData.load_controller{16}.ATTR.id=16;
febio_spec.LoadData.load_controller{16}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{16}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{16}.points.point.VAL=Mat_RX_Capitate;

febio_spec.LoadData.load_controller{17}.ATTR.id=17;
febio_spec.LoadData.load_controller{17}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{17}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{17}.points.point.VAL=Mat_RY_Capitate;

febio_spec.LoadData.load_controller{18}.ATTR.id=18;
febio_spec.LoadData.load_controller{18}.ATTR.type='loadcurve';
febio_spec.LoadData.load_controller{18}.interpolate='LINEAR';
febio_spec.LoadData.load_controller{18}.points.point.VAL=Mat_RZ_Capitate;

febio_spec.Output.logfile.ATTR.file=febioLogFileName;
febio_spec.Output.logfile.node_data{1}.ATTR.file=febioLogFileName_disp;
febio_spec.Output.logfile.node_data{1}.ATTR.data='ux;uy;uz';
febio_spec.Output.logfile.node_data{1}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{1}.VAL=1:size(V,1);

febio_spec.Output.logfile.node_data{2}.ATTR.file=febioLogFileName_force;
febio_spec.Output.logfile.node_data{2}.ATTR.data='Rx;Ry;Rz';
febio_spec.Output.logfile.node_data{2}.ATTR.delim=',';
febio_spec.Output.logfile.node_data{2}.VAL=1:size(V,1);

febio_spec.Output.logfile.element_data{1}.ATTR.file=febioLogFileName_stress;
febio_spec.Output.logfile.element_data{1}.ATTR.data='sx;sy;sz';
febio_spec.Output.logfile.element_data{1}.ATTR.delim=',';
%% Quick viewing of the FEBio input file structure
%The febView function can be used to view the xml structure in a MATLAB figure window.
%
%disp('Viewing the febio file');
%febView(febio_spec); %Viewing the febio file

%% Exporting the FEBio input file
% Exporting the febio_spec structure to an FEBio input file is done using
% the |febioStruct2xml| function.
febioStruct2xml(febio_spec,febioFebFileName); %Exporting to file and domNode

% %
% %
% % %% Running the FEBio analysis
% % % To run the analysis defined by the created FEBio input file the
% % % |runMonitorFEBio| function is used. The input for this function is a
% % % structure defining job settings e.g. the FEBio input file name. The
% % % optional output runFlag informs the user if the analysis was run
% % % succesfully.

febioAnalysis.run_filename=febioFebFileName; %The input file name
febioAnalysis.run_logname=febioLogFileName; %The name for the log file
febioAnalysis.disp_on=1; %Display information on the command window
febioAnalysis.disp_log_on=1; %Display convergence information in the command window
febioAnalysis.runMode='internal';%'internal';
febioAnalysis.t_check=0.25; %Time for checking log file (dont set too small)
febioAnalysis.maxtpi=1e99; %Max analysis time
febioAnalysis.maxLogCheckTime=3; %Max log file checking time

[runFlag]=runMonitorFEBio(febioAnalysis);%START FEBio NOW!!!!!!!!

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
clear all

disp('Data loading and generation...');
ima = double(imread('lena.png'));
L = 3;
s = zeros(size(ima));
for k = 1:L
    s = s + abs(randn(size(ima)) + i * randn(size(ima))).^2 / 2;
end
ima_nse = ima .* sqrt(s / L);

disp('Non-it PPB filter...');
hW = 10;
hD = 3;
alpha = 0.88;
T = inf;
nbit = 1;
ima_fil_1 = ppb_nakagami(ima_nse, L, ...
                         hW, hD, ...
                         alpha, T, ...
                         nbit);

disp('It PPB filter...');
hW = [1 3 5 10];
hD = [0 1 2 3];
alpha = 0.92;
T = 0.2;
nbit = [1 2 3 4];
ima_fil_2 = ppb_nakagami(ima_nse, L, ...
                         hW, hD, ...
                         alpha, T, ...
                         nbit);

figure
subplot(2,2,1), plotimage_sar(ima);
title('Noise-free image');
subplot(2,2,2), plotimage_sar(ima_nse, ima);
title(sprintf('Noisy image (PSNR %.3f)', snr(ima_nse, ima)));
subplot(2,2,3), plotimage_sar(ima_fil_1, ima);
title(sprintf('Non-it PPB filter (PSNR %.3f)', snr(ima_fil_1, ima)));
subplot(2,2,4), plotimage_sar(ima_fil_2, ima);
title(sprintf('It PPB filter (PSNR %.3f)', snr(ima_fil_2, ima)));

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
% Synchronoise the streams such that they will work with the network.
% Things did in this function:
%   1. upsample streams to match other streams of higher frame rate. This
%   is only done when streams are specified by DataSyncSet to have the same
%   frame rate. Also discarded last frames longer streams to make the
%   streams have the same number of frames.
%   2. apply VAD to some streams. Here VAD is stored in one stream
%   (normally the last stream). We first upsample/downsample VAD framerate
%   to match the stream to be VADed. Then remove nonspeech segments. VAD
%   should be stored as 1-D feature stream.
%
% Author: Xiong Xiao, Temasek labs, NTU, Singapore.
% Date Created: 10 Oct 2015
% Last Modified: 21 Apr 2016
%
function [data_sync, nFr_actual, isTensor, isVAD, isFramewise] = SynchronizeDataStreams2(alldata, para)
[nStream, nUtt] = size(alldata);

isTensor = para.IO.isTensor;
isVAD = zeros(nStream,1);
isFramewise = zeros(nStream,1);
data_sync = alldata;
nFr_actual = zeros(nUtt, nStream);

ApplyVADSet = para.IO.ApplyVADSet;
DataSyncSet = para.IO.DataSyncSet;
if isempty(DataSyncSet) && isempty(ApplyVADSet)
    for utt_i = 1:nUtt
        for i=1:nStream
            [~,n2,n3] = size(data_sync{i, utt_i});
            if isTensor(i)
                if n2==0
                    nFr_actual(utt_i,i) = 0;
                else
                    nFr_actual(utt_i,i) = n3;
                end
            else
                nFr_actual(utt_i,i) = n2;
            end
            if n2>1
                isFramewise(i) = 1;     % if any sentence has more than one frame, it is framewise feature or label
            end
        end
    end
    return;
end

for utt_i = 1:nUtt
    % Apply VAD and segmentation first, if required
    for vi=1:length(ApplyVADSet)
        isVAD(ApplyVADSet{vi}(end))=1;
        curr_vad_set = ApplyVADSet{vi};
        if isempty(curr_vad_set); continue; end
        data_sync{curr_vad_set(1), utt_i} = SynchronizeStreamsVAD( data_sync(curr_vad_set, utt_i), para.IO.vadAction{vi} );
    end

    % Apply Frame rate synchronization to all Sync Sets
    for di=1:length(DataSyncSet)
        curr_sync_set = DataSyncSet{di};
        if isempty(curr_sync_set); continue; end
        data_sync(curr_sync_set, utt_i) = SynchronizeStreamsFrameRate(  data_sync(curr_sync_set, utt_i) , para.IO.frame_rate(curr_sync_set) , isTensor(curr_sync_set)  );
    end

%     for i=1:nStream
%         [~,~,n3] = size(data_sync{i, utt_i});
%         if n3>1; isTensor(i) = 1;   end
%     end
end
for utt_i = 1:nUtt
    for i=1:nStream
        [~,n2,n3] = size(data_sync{i, utt_i});
        if isTensor(i)
            if n2==0
                nFr_actual(utt_i,i) = 0;
            else
                nFr_actual(utt_i,i) = n3;
            end
        else
            nFr_actual(utt_i,i) = n2;
        end
        if n2>1
            isFramewise(i) = 1;     % if any sentence has more than one frame, it is framewise feature or label
        end
    end
end
end


text:
function obstacle_pts = random_obstacles(dim, n_obs, lb, ub, scale)
% Generate a randomly distributed field of random obstacles. Each
% obstacle consists of 2^dim points uniformly randomly distributed
% within a 2*scale-length cube around a randomly selected center point
if nargin < 5
  scale = 1;
end
offsets = scale * (rand(dim, 2^dim * n_obs) * 2 - 1);
centers = bsxfun(@plus, bsxfun(@times, rand(dim, n_obs), ub - lb), lb);
centers = reshape(centers, dim, 1, []);
obstacle_pts = bsxfun(@plus, centers, reshape(offsets ./ (n_obs^(1/dim)), [dim, 2^dim, n_obs]));

end

author: rdeits

repo: iris-distro

save_path: github-repos/MATLAB/rdeits-iris-distro


text:
% LFUnpackRawBuffer - Unpack a buffer of packed raw binary data into an image
%
% Usage:
%
%     ImgOut = LFUnpackRawBuffer( Buff, BitPacking, ImgSize )
%
% Used by LFReadRaw and LFReadLFP, this helper function unpacks a raw binary data buffer in one of several formats.
%
% Inputs :
%
%     Buff : Buffer of chars to be unpacked
%     BitPacking : one of '12bit', '10bit' or '16bit'; default is '12bit'
%     ImgSize : size of the output image
%
% Outputs:
%
%      Img : an array of uint16 gray levels. No demosaicing (decoding Bayer pattern) is performed.
%
% User guide: LFToolbox.pdf
% See also: LFDecodeLensletImageDirect, demosaic

% Copyright (c) 2013-2020 Donald G. Dansereau

function ImgOut = LFUnpackRawBuffer( Buff, BitPacking, ImgSize )

switch( BitPacking )
    case '8bit'
        ImgOut = reshape(Buff, ImgSize);

    case '10bit'
        t0 = uint16(Buff(1:5:end));
        t1 = uint16(Buff(2:5:end));
        t2 = uint16(Buff(3:5:end));
        t3 = uint16(Buff(4:5:end));
        lsb = uint16(Buff(5:5:end));

        t0 = bitshift(t0,2);
        t1 = bitshift(t1,2);
        t2 = bitshift(t2,2);
        t3 = bitshift(t3,2);

        t0 = t0 + bitand(lsb,bin2dec('00000011'));
        t1 = t1 + bitshift(bitand(lsb,bin2dec('00001100')),-2);
        t2 = t2 + bitshift(bitand(lsb,bin2dec('00110000')),-4);
        t3 = t3 + bitshift(bitand(lsb,bin2dec('11000000')),-6);

        ImgOut = zeros(ImgSize, 'uint16');
        ImgOut(1:4:end) = t0;
        ImgOut(2:4:end) = t1;
        ImgOut(3:4:end) = t2;
        ImgOut(4:4:end) = t3;

    case '12bit'
        t0 = uint16(Buff(1:3:end));
        t1 = uint16(Buff(2:3:end));
        t2 = uint16(Buff(3:3:end));

        a0 = bitshift(t0,4) + bitshift(bitand(t1,bin2dec('11110000')),-4);
        a1 = bitshift(bitand(t1,bin2dec('00001111')),8) + t2;

        ImgOut = zeros(ImgSize, 'uint16');
        ImgOut(1:2:end) = a0;
        ImgOut(2:2:end) = a1;

    case '16bit'
        t0 = uint16(Buff(1:2:end));
        t1 = uint16(Buff(2:2:end));
        a0 = bitshift(t1, 8) + t0;
        ImgOut = zeros(ImgSize, 'uint16');
        ImgOut(:) = a0;

    otherwise
        error('Unrecognized bit packing');
end

repo: LFToolbox


text:
function stress=outputResult(displacements,numberElements,...
        elementNodes,nodeCoordinates,D)
stress=zeros(3,1);
for e=1:numberElements
  numNodePerElement = length(elementNodes(e,:));
  numEDOF = 2*numNodePerElement;
  elementDof=zeros(1,numEDOF);
  for i = 1:numNodePerElement
      elementDof(2*i-1)=2*elementNodes(e,i)-1;
      elementDof(2*i)=2*elementNodes(e,i);
  end

  %  B matrix
  x1 = nodeCoordinates(elementNodes(e,1),1);
  y1 = nodeCoordinates(elementNodes(e,1),2);
  x2 = nodeCoordinates(elementNodes(e,2),1);
  y2 = nodeCoordinates(elementNodes(e,2),2);
  x3 = nodeCoordinates(elementNodes(e,3),1);
  y3 = nodeCoordinates(elementNodes(e,3),2);
  A = 1/2*det([1 x1 y1; 1 x2 y2; 1 x3 y3]);
  B = 1/(2*A).*[y2-y3 0 y3-y1 0 y1-y2 0;
                        0 x3-x2 0 x1-x3 0 x2-x1;
                        x3-x2 y2-y3 x1-x3 y3-y1 x2-x1 y1-y2];

  stress(e,:)=stress+D*B*displacements(elementDof)/2;
end

repo: Aero-matlab


text:
%mdd is a class.
%    obj = mdd(theta_draws, log_post_kern, lb, ub, maximization)
%
%    Reference page in Doc Center
%       doc mdd
%
%


text:
function poly2 = densifyPolygon(poly, N)
%DENSIFYPOLYGON Add several points on each edge of the polygon.
%
%   POLY2 = densifyPolygon(POLY, N)
%   POLY is a NV-by-2 array containing polygon coordinates. The function
%   iterates on polygon edges, divides it into N subedges (by inserting N-1
%   new vertices on each edges), and return the resulting polygon.
%   The new polygon POLY has therefore N*NV vertices.
%
%   Example
%     % Densifies a simple polygon
%     poly = [0 0 ; 10 0;5 10;15 15;5 20;-5 10];
%     poly2 = densifyPolygon(poly, 10);
%     figure; drawPolygon(poly); axis equal
%     hold on; drawPoint(poly2);
%
%   See also
%     drawPolygon, edgeToPolyline
%

% ------
% Author: David Legland
% E-mail: david.legland@grignon.inra.fr
% Created: 2011-11-25, using Matlab 7.9.0.529 (R2009b)
% Copyright 2011-2022 INRA - Cepia Software Platform

% number of vertices, and of edges
Nv = size(poly, 1);

% number of vertices in new polygon
N2 = N * Nv;
poly2 = zeros(N2, 2);

% iterate on polygon edges
for i = 1:Nv
    % extract current edge
    v1 = poly(i, :);
    v2 = poly(mod(i, Nv) + 1, :);

    % convert current edge to polyline
    newVertices = edgeToPolyline([v1 v2], N);

    % indices of current polyline to resulting polygon
    i1 = (i-1)*N + 1;
    i2 = i * N;

    % fill up polygon
    poly2(i1:i2, :) = newVertices(1:end-1, :);
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function meshI = genIVmesh(mesh)

Np = size(mesh.p,1);
node = [mesh.p;mesh.eIntP];
allCutElem = mesh.t(mesh.tLoc<0,:);
isCutElem = (mesh.tLoc<0);
vSign = [mesh.pLoc;zeros(size(mesh.eIntP,1),1)];
isInterfaceNode = false(size(node,1), 1);
isInterfaceNode(allCutElem(:)) = true; % include the vertices of allCutElem
numCutElemV1 = sum(isInterfaceNode);
isInterfaceNode(Np+1:end) = true;    % and the cut points and the aux points
allCutElemReoderTmp = zeros(size(node,1),1);
numCutElemV2 = sum(isInterfaceNode);
allCutElemReoderTmp(isInterfaceNode) = 1:numCutElemV2;
allCutElemReoder = allCutElemReoderTmp(allCutElem);
interfaceNode = node(isInterfaceNode,:);
ntI = -min(mesh.tLoc);
intID = find(mesh.tLoc<0);

tetElem = zeros(12*ntI,4);
tetElemLoc = zeros(12*ntI,1);
idx2cube = zeros(12*ntI,1);
tetcount = 0;

for i = 1:ntI
    if i == 3
        stp = 1;
    end
    tID = intID(i);
    t_e = mesh.t_e(tID,:);
    pLocK = mesh.pLoc(mesh.t(tID,:));
    vert0 = mesh.p(mesh.t(tID,:),:);
    nodeid = [allCutElemReoder(i,:),numCutElemV1-mesh.eLoc(t_e(mesh.eLoc(t_e)<0))'];
    vert2 = vert0(pLocK>0,:); % plus domain
    vert1 = vert0(pLocK<0,:); % minus domain
    intpt0 = mesh.eIntP(-mesh.eLoc(t_e(mesh.eLoc(t_e)<0)),:);
    id1 = [find(pLocK<0)', 5:4+size(intpt0,1)];
    id2 = [find(pLocK>0)', 5:4+size(intpt0,1)];
    p1 = [vert1;intpt0]; DT = delaunayTriangulation(p1); t1 = DT.ConnectivityList;
    p2 = [vert2;intpt0]; DT = delaunayTriangulation(p2); t2 = DT.ConnectivityList;
    tetElem(tetcount+1:tetcount+size(t1,1),:) = nodeid(id1(t1));
    idx2cube(tetcount+1:tetcount+size(t1,1)) = tID;
    tetElemLoc(tetcount+1:tetcount+size(t1,1)) = 1; % inside subdomain
    tetcount = tetcount+size(t1,1);
    tetElem(tetcount+1:tetcount+size(t2,1),:) = nodeid(id2(t2));
    idx2cube(tetcount+1:tetcount+size(t2,1)) = tID;
    tetElemLoc(tetcount+1:tetcount+size(t2,1)) = 2; % outside subdomain
    tetcount = tetcount+size(t2,1);
end
tetElem(tetcount+1:end,:) = [];
tetElemLoc(tetcount+1:end,:) = [];
idx2cube(tetcount+1:end,:) = [];
[tetElem, ortidx, volume] = fixorder3(interfaceNode, tetElem);

% there are some tets which are coplane, need to get rid of them
d12 = interfaceNode(tetElem(:,2),:) - interfaceNode(tetElem(:,1),:);
d13 = interfaceNode(tetElem(:,3),:) - interfaceNode(tetElem(:,1),:);
d14 = interfaceNode(tetElem(:,4),:) - interfaceNode(tetElem(:,1),:);
d12 = sum(d12.^2,2).^(1/2); d13 = sum(d13.^2,2).^(1/2); d14 = sum(d14.^2,2).^(1/2);
ld = max([d12,d13,d14],[],2);
idPlane = (volume./ld<=10^(-16));

tetElem(idPlane,:) = [];
volume(idPlane,:) = [];
tetElemLoc(idPlane,:) = [];
idx2cube(idPlane,:) = [];
localidx2globalidx = find(isInterfaceNode); % tetElem points to interfaceNode
tetElem = localidx2globalidx(tetElem);  % map to the global index of node
% there are some tets which are on the interface, but contained in the surrounding tets
% need to get rid of them (but the following algorithm may not be robust)
tetSign = vSign(tetElem);
idSliver = (sum(abs(tetSign),2)==0);
tetElem(idSliver,:) = [];
volume(idSliver,:) = [];
tetElemLoc(idSliver,:) = [];
idx2cube(idSliver,:) = [];
PolyVolume = accumarray([idx2cube,tetElemLoc],volume);
PolyVolume = PolyVolume(mesh.tLoc<0,:);

%% Get triangular faces for interrior elements and interface
localFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2];
NT = size(tetElem,1);
tface = zeros(4*NT, 3);
tface2elem = zeros(4*NT, 1);
iface = zeros(4*NT, 3);
iface2elem = zeros(4*NT, 1);
% find the interior tet elements
isIntTet1 = min(vSign(tetElem),[], 2) == -1; % can not be == -1 as there is sliver
isIntTet2 = sum(abs(vSign(tetElem)),2) == 0;
isIntTet = isIntTet1 | isIntTet2;
intTet = tetElem(isIntTet,:);
% find the corresponding cube indices
intIdx2cube = idx2cube(isIntTet);
% find triangular interface
T = auxstructure3(intTet);
neighbor = T.neighbor; % if a face is on the boundary, then the neighbor element index is itself
clear T;
tmp = (1:size(intTet, 1))';
ct = 0;
ci = 0;
for i = 1:4
    face = intTet(:, localFace(i,:));
    % find the triangle faces of polyhedron
    % 1. face and its neighbor associated to different cubes, and
    % 2. face is not on the boundary of all cut elements (which are squares)
    isPolyTriFace = ((neighbor(:, i) == tmp) | (intIdx2cube ~= intIdx2cube(neighbor(:,i)))) &...
        (sum(abs(vSign(face)), 2) >= 10^(-12));% & (sum(abs(vSign(face)), 2) > 0);
    c2 = ct + sum(isPolyTriFace);
    tface((ct+1):c2,:) = face(isPolyTriFace,:);
    tface2elem((ct+1):c2,:) = intIdx2cube(isPolyTriFace); % the indices of the polyhedron
    ct = c2;
    % note that only interior elements are treated

    % find the triangle faces on interface with normal points to exterior
    % 1. face is on the boundary and
    % 2. all vertices are on the interface
    isInterface = (sum(abs(vSign(face)), 2) <= 10^(-12));

    % add to interface
    c4 = ci + sum(isInterface);
    iface((ci+1):c4,:) = face(isInterface,:); % the interface tri faces.
    iface2elem((ci+1):c4,:) = intIdx2cube(isInterface); % the indices of the polyhedron
    ci = c4;
end
iface((ci+1):end,:) = [];
iface2elem((ci+1):end,:) = [];
c2old = c2;
% plot to check
%trisurf(tface(1:c2,:),node(:,1),node(:,2),node(:,3))
%trisurf(iface(1:ci,:),node(:,1),node(:,2),node(:,3))

% Find the triangular faces for exterior elements
extTet = tetElem(~isIntTet, :);
extIdx2cube = idx2cube(~isIntTet);

T = auxstructure3(extTet);
neighbor = T.neighbor;
clear T;
tmp = (1:size(extTet,1))';
for i = 1:4
    face = extTet(:, localFace(i,:));
    % find the triangle faces of polyhedron
    % 1. face and its neighbor associated to different cubes, and
    % 2. face is not on the boundary of all cut elements (which are squares)
    isPolyTriFace = ((neighbor(:, i) == tmp) | (extIdx2cube ~= extIdx2cube(neighbor(:,i)))) &...
        (sum(abs(vSign(face)), 2) > 0);
    c2 = ct + sum(isPolyTriFace);
    tface((ct+1):c2,:) = face(isPolyTriFace,:);
    tface2elem((ct+1):c2,:) = extIdx2cube(isPolyTriFace);
    % index of exterior polyhedron is append to the end of elem
    ct = c2;

end

% plot to check
%trisurf(tface(c2old+1:c2,:),node(:,1),node(:,2),node(:,3))

tface((ct+1):end,:) = [];    % remove empty meomory
tface2elem((ct+1):end) = [];
face2elemLoc = [ones(c2old,1);2*ones(c2-c2old,1)]; % face2elemLoc is the same as face location

meshI.tface = tface;
meshI.tface2elem = tface2elem;
meshI.iface = iface;
meshI.iface2elem = iface2elem;
meshI.face2elemLoc = face2elemLoc;
meshI.PolyVolume = PolyVolume;
% meshI.vSign = vSign;
meshI.node = node;
meshI.tetElem = tetElem;
meshI.tetElemLoc = tetElemLoc;
meshI.idx2cube = idx2cube;
meshI.tetVolume = volume;
meshI.vSign = vSign;

end