text:
function varargout = metchgui(varargin)
%    alldata = metchgui(node,elem,points,pface) or metchgui(volume,points,pface)
%
%    A GUI to register a point cloud to a mesh or volumetric image
%
%    author: Qianqian Fang <q.fang at neu.edu>
%    date: 12/16/2008
%
%   parameters:
%        node: node coordinate of the surface mesh (nn x 3)
%        elem: element list of the surface mesh (3 columns for
%              triangular mesh, 4 columns for cubic surface mesh)
%        points: the coordinates (3 columns for x/y/z) of the
%              point cloud which you want to register
%        pface:trianglular surface defined on the point cloud.
%              pface is optional; if presents, metch will display
%              a surface object instead of a point cloud.
%
%   the input can also be two parameters in form of metchgui(volume,points),
%    where volume is a 3D image (array).
%
%   outputs:
%        alldata: a structrure containing all processing outputs
%        the fields include:
%         .node: the input node
%         .elem: the input surface mesh elements
%         .volume: if the input volumetric image
%         .A0: the affine rotation for selected point pairs (after Initialize)
%         .b0: the affine translation for selected point pairs (after Initialize)
%         .A: the affine rotation for the point cloud (after Optimize)
%         .b: the affine translation for the point cloud (after Optimize)
%         .points: the input point cloud
%         .pointsinit: the point cloud after initialization
%         .pointsopt: the point cloud after optimization
%         .pointsproj: the point cloud after projecting to the surface
%         .initplot: the handle to the point cloud plot after init
%         .optplot: the handle to the point cloud plot after optimization
%         .projplot: the handle to the point cloud plot after projection
%
%   If user supplys an output variable, the GUI will not return until the
%   user hits the "close" button or close the window; if user does not
%   supply any output, the call will return immediately; any data user
%   intends to save, he has to click on "Save Session" button and provides
%   a mat-file file name. A single structure named "metchsession" will be
%   stored in this file.
%
%   example: (meshasphere/meshunitsphere are defined in iso2mesh http://iso2mesh.sf.net)
%
%       [noderef,faceref,elemref]=meshunitsphere(0.08,10);
%       [no,fc]=removeisolatednode(noderef(:,1:3),faceref(:,1:3));
%       [node,face,elem]=meshasphere([10 20 15],3,0.5,10);
%       [no2,fc2]=removeisolatednode(node(:,1:3),face(:,1:3));
%       alldata = metchgui(no,fc,no2);
%       % or alldata = metchgui(no,fc,no2,fc2);
%
%   Please find more information at http://iso2mesh.sf.net/cgi-bin/index.cgi?metch
%
%   this function is part of "metch" toobox, see COPYING for license

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

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

% --- Executes just before metchgui is made visible.
function metchgui_OpeningFcn(hObject, eventdata, handles, varargin)
handles.hasoutput=0;
if(isempty(varargin))
    fprintf(1,'Metch GUI must be called with parameters:\nFormat: alldata = metchgui(node,elem,points,pface);\n');
    close(handles.MetchGUI);
    return;
end
if(ischar(varargin{end}) && strcmp(varargin{end},'hasoutput'))
        handles.hasoutput=1;
        varargin(end)=[];
end
handles.output = hObject;

set(handles.btAddMeshPt,'userdata',[handles.axMesh,handles.btAddMeshPt,handles.btAddCloudPt]);
set(handles.btAddCloudPt,'userdata',[handles.axMesh,handles.btAddCloudPt,handles.btAddMeshPt]);

% if uses supplied 2 input variables, assume a volume image and a point cloud
if(isnumeric(varargin{1}) && length(size(varargin{1}))==3)
       vol=varargin{1};
       pt=varargin{2};

       dat.volume=vol;
       dat.points=pt;

       slice=round(size(vol,3)/2);
       hs=imagesc(vol(:,:,slice),'parent',handles.axMesh);
       set(handles.slPos,'max',size(vol,3),'min',1,'value',slice);

       if(length(varargin)>=3)
       	 pface=varargin{3};
       	 dat.pface=pface;
	 trisurf(pface(:,1:3),pt(:,1),pt(:,2),pt(:,3),'parent',handles.axPoints);
       else
         %plot3(pt(:,1),pt(:,2),pt(:,3),'.','parent',handles.axPoints);
         ptcolor=(pt-repmat(min(pt),size(pt,1),1))./repmat(max(pt)-min(pt),size(pt,1),1);
         drawnow;
         scatter3(pt(:,1),pt(:,2),pt(:,3),3,ptcolor,'filled');
       end
       set(handles.MetchGUI,'userdata',dat);

       axis(handles.axMesh,'equal');
       axis(handles.axPoints,'equal');
       axis(handles.axMesh,'off');
       grid(handles.axPoints,'on');
       %axis(handles.axPoints,'off');

       set(handles.axMesh,'tag','axMesh');
       set(handles.axPoints,'tag','axPoints');
       set(handles.slPos,'visible','on');
       set(handles.lbZPos,'visible','on');
       rotate3d(handles.axPoints,'on');
       rotate3d(gcf,'on');
end

% if uses supplied 3 input variables, assume a surface mesh and a point cloud/surface
if(length(varargin)>=3 && length(size(varargin{1}))==2)
       node=varargin{1};
       elem=varargin{2};
       pt=varargin{3};

       dat.node=node;
       dat.elem=elem;
       dat.points=pt;
       if(length(varargin)>=4)
       	 pface=varargin{4};
       	 dat.pface=pface;
       end
       set(handles.MetchGUI,'userdata',dat);

       drawinit(handles,dat)
end
guidata(hObject, handles);

if(handles.hasoutput)
        uiwait(handles.MetchGUI);
end

%---------------------------------------------------------------------------
function drawinit(handles,dat)

node=dat.node;
elem=dat.elem;
pt=dat.points;

hs=trisurf(elem,node(:,1),node(:,2),node(:,3),'parent',handles.axMesh);
%set(hs,'linestyle','none');
%set(hs,'facecolor','b','facealpha',0.8);

%plot3(pt(:,1),pt(:,2),pt(:,3),'.','parent',handles.axPoints);
if(~isfield(dat,'pface'))
	ptcolor=(pt-repmat(min(pt),size(pt,1),1))./repmat(max(pt)-min(pt),size(pt,1),1);
	drawnow;
	scatter3(pt(:,1),pt(:,2),pt(:,3),3,ptcolor,'filled');
else
	trisurf(dat.pface,pt(:,1),pt(:,2),pt(:,3),'parent',handles.axPoints);
end
%hold(handles.axPoints,'on');
%plot3(pt(5:7,1),pt(5:7,2),pt(5:7,3),'ro','parent',handles.axPoints);

axis(handles.axMesh,'equal');
axis(handles.axPoints,'equal');
axis(handles.axMesh,'off');
%axis(handles.axPoints,'off');

set(handles.axMesh,'tag','axMesh');
set(handles.axPoints,'tag','axPoints');

rotate3d(handles.axPoints,'on');
rotate3d(handles.axMesh,'on');
rotate3d(gcf,'on');

%---------------------------------------------------------------------------
function varargout = metchgui_OutputFcn(hObject, eventdata, handles)
if(length(handles) && isfield(handles,'output') && isfield(handles,'hasoutput') && handles.hasoutput)
        handles.output=get(handles.MetchGUI,'userdata');
        varargout{1} =handles.output;
        close(handles.MetchGUI);
end
%---------------------------------------------------------------------------
function lbMesh_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%---------------------------------------------------------------------------
function lbPoints_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%---------------------------------------------------------------------------
function edit1_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%---------------------------------------------------------------------------
function edit2_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%---------------------------------------------------------------------------
function edit3_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

%---------------------------------------------------------------------------
function isSelect_Callback(hObject, eventdata, handles)
if(get(hObject,'Value'))
        datacursormode(gcf,'on');
        set(datacursormode(gcf),'UpdateFcn',@myupdatefcn)
else
        datacursormode(gcf,'off');
        rotate3d(gcf,'on');
end

%---------------------------------------------------------------------------
% the respond function when there is a data-tip to popup
%---------------------------------------------------------------------------
function txt=myupdatefcn(empt,event_obj)
pos = get(event_obj,'Position');
idx=  get(event_obj,'DataIndex');

if(length(pos)==3)
  txt = {['x: ',num2str(pos(1))],...
         ['y: ',num2str(pos(2))],['z: ',num2str(pos(3))],['index:', num2str(idx)]};
elseif(length(pos)==2)
  txt = {['x: ',num2str(pos(1))],['y: ',num2str(pos(2))]};
end
targetup=get(get(event_obj,'Target'),'parent');
set(targetup,'userdata',struct('pos',pos,'idx',idx));
if(targetup==findobj('tag','axMesh'))
         set(findobj('tag','btAddMeshPt'),'enable','on');
         set(findobj('tag','btAddCloudPt'),'enable','off');
elseif(targetup==findobj('tag','axPoints'))
         set(findobj('tag','btAddMeshPt'),'enable','off');
         set(findobj('tag','btAddCloudPt'),'enable','on');
end

%---------------------------------------------------------------------------
function bInit_Callback(hObject, eventdata, handles)
mapto=get(handles.lbMesh,'userdata');
mapfrom=get(handles.lbPoints,'userdata');

maptoidx=get(handles.txMapTo,'userdata');
mapfromidx=get(handles.txMapFrom,'userdata');

if(length(mapto)<4 | length(mapfrom)<4)
        msgbox('You have to select >3 points from the point cloud plot and corresponding points from the mesh','Error','error');
        return;
end
[A,b]=affinemap(mapfrom,mapto)
dat=get(handles.MetchGUI,'userdata');
dat.A0=A;
dat.b0=b;

newpt=(A*dat.points'+repmat(b(:),1,size(dat.points,1)))';
dat.pointsinit=newpt;

hold(handles.axMesh,'on');
if(isfield(dat,'initplot'))
        delete dat.initplot;
        dat.initplot=0;
end
dat.initplot=plot3(newpt(:,1),newpt(:,2),newpt(:,3),'r.','parent',handles.axMesh);

dat.fromidx=mapfromidx;
dat.toidx=maptoidx;
set(handles.MetchGUI,'userdata',dat);

%---------------------------------------------------------------------------
function axPoints_ButtonDownFcn(hObject, eventdata, handles)
if(get(handles.btSelect,'value')==1)
        pp=getCursorInfo(datacursormode(gcf));
        str=[get(handles.lbPoints,'string');mat2str(pp.Position)];
        set(handles.lbPoints,'string',str);
end

%---------------------------------------------------------------------------
function addselectedpt(pos,idx,lb)
if(isempty(get(lb,'value')))
    set(lb,'value',1);
end
listpt=get(lb,'string');
listpt{end+1}=[num2str(idx) ':' mat2str(pos)];
set(lb,'string',listpt);
set(lb,'userdata',[get(lb,'userdata');pos]);

%---------------------------------------------------------------------------
function btAddMeshPt_Callback(hObject, eventdata, handles)
dat=get(handles.axMesh,'userdata');

dat0=get(handles.MetchGUI,'userdata');
if(isfield(dat0,'volume'))
    dat.pos(:,3)=get(handles.slPos,'value');
end

if(isfield(dat,'pos'))
    addselectedpt(dat.pos,dat.idx,handles.lbMesh);
    set(handles.txMapTo,'userdata',[get(handles.txMapTo,'userdata');dat.idx]);
	dat0.mapto=get(handles.lbMesh,'userdata');
    dat0.maptoidx=get(handles.txMapTo,'userdata');
    set(handles.MetchGUI,'userdata',dat0);
else
        msgbox('No point was selected. Please click on "Select" and select a point on the mesh or point cloud','Error','error');
        return;
end

%---------------------------------------------------------------------------
function btAddCloudPt_Callback(hObject, eventdata, handles)
dat=get(handles.axPoints,'userdata');
dat0=get(handles.MetchGUI,'userdata');
if(isfield(dat,'pos'))
    addselectedpt(dat.pos,dat.idx,handles.lbPoints);
    set(handles.txMapFrom,'userdata',[get(handles.txMapFrom,'userdata');dat.idx]);
	dat0.mapfrom=get(handles.lbPoints,'userdata');
    dat0.mapfromidx=get(handles.txMapFrom,'userdata');
    set(handles.MetchGUI,'userdata',dat0);
else
        msgbox('No point was selected. Please click on "Select" and select a point on the mesh or point cloud','Error','error');
        return;
end

%---------------------------------------------------------------------------
function btOptimize_Callback(hObject, eventdata, handles)
dat=get(handles.MetchGUI,'userdata');
if(isfield(dat,'A0') && isfield(dat,'b0')& isfield(dat,'node')& isfield(dat,'elem')& ...
   isfield(dat,'pointsinit')& isfield(dat,'toidx')& isfield(dat,'fromidx'))
        pmask=-1*ones(size(dat.pointsinit,1),1);
        pmask(dat.fromidx)=dat.toidx;
	if(isfield(dat,'A') && isfield(dat,'b'))
	        [Anew,bnew,posnew]=regpt2surf(dat.node,dat.elem,dat.points,pmask,dat.A,dat.b,ones(12,1),20);
	else
                [Anew,bnew,posnew]=regpt2surf(dat.node,dat.elem,dat.points,pmask,dat.A0,dat.b0,ones(12,1),20);
	end
        dat.A=Anew;
        dat.b=bnew;
        dat.pointsopt=posnew;
        if(isfield(dat,'optplot') && dat.optplot)
                delete dat.optplot;
                dat.optplot=0;
        end
        dat.optplot=plot3(posnew(:,1),posnew(:,2),posnew(:,3),'g+','parent',handles.axMesh);
else
        msgbox('You have to select 4 points and click "Initialize" button first','Error','error');
        return;
end
set(handles.MetchGUI,'userdata',dat);

%---------------------------------------------------------------------------
function slPos_Callback(hObject, eventdata, handles)
dat=get(handles.MetchGUI,'userdata');
if(isfield(dat,'volume'))
        hold(handles.axMesh,'off');
        imagesc(dat.volume(:,:,round(get(hObject,'value'))),'parent',handles.axMesh);
        set(handles.axMesh,'tag','axMesh');
        set(handles.lbZPos,'string',num2str(round(get(hObject,'value'))));
end
set(handles.MetchGUI,'userdata',dat);

%---------------------------------------------------------------------------
function slPos_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor',[.9 .9 .9]);
end

%--------------------------------------------------------------------------
function btProj_Callback(hObject, eventdata, handles)
dat=get(handles.MetchGUI,'userdata');
if(isfield(dat,'pointsopt'))
        if(isfield(dat,'projplot') && dat.projplot)
                delete dat.projplot;
                dat.projplot=0;
        end
        nv=nodesurfnorm(dat.node,dat.elem);
        [d2surf,cn]=dist2surf(dat.node,nv,dat.pointsopt);
        [dat.pointsproj dat.elemid dat.weight]=proj2mesh(dat.node,dat.elem,dat.pointsopt,nv,cn);
        hold(handles.axMesh,'on');
        dat.projplot=plot3(dat.pointsproj(:,1),dat.pointsproj(:,2),dat.pointsproj(:,3),'c.','parent',handles.axMesh);
else
        msgbox('You have to first select 4 points, then click "Initialize" and "Optimize" button','Error','error');
        return;
end
set(handles.MetchGUI,'userdata',dat);

%--------------------------------------------------------------------------
function btSaveRes_Callback(hObject, eventdata, handles)
[filename, pathname] = uiputfile('*.mat', 'Save Metch Workspace as');
mapto=get(handles.lbMesh,'userdata');
mapfrom=get(handles.lbPoints,'userdata');

maptoidx=get(handles.txMapTo,'userdata');
mapfromidx=get(handles.txMapFrom,'userdata');

metchsession=get(handles.MetchGUI,'userdata');
if(~isempty(mapto))    metchsession.mapto=mapto; end
if(~isempty(mapfrom))    metchsession.mapfrom=mapfrom; end
if(~isempty(maptoidx))    metchsession.maptoidx=maptoidx; end
if(~isempty(mapfromidx))    metchsession.mapfromidx=mapfromidx; end
fname=[pathname filename];
save(fname,'metchsession');

%--------------------------------------------------------------------------
function btLoadSession_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile('*.mat', 'Load Metch Workspace from');
fname=[pathname filename];
load(fname);
handle.output=metchsession;
cla(handles.axMesh);
cla(handles.axPoints);
drawinit(handles,metchsession.node,metchsession.elem,metchsession.points);
%--------------------------------------------------------------------------
function btPlotResults_Callback(hObject, eventdata, handles)

function btClose_Callback(hObject, eventdata, handles)
guidata(hObject,handles);
uiresume;

function btHelp_Callback(hObject, eventdata, handles)
helpmsg={
'Metch GUI: A mesh/volume registration toolbox',
'',
'Author: Qianqian Fang <q.fang at neu.edu>',
'        Martinos Center for Biomedical Imaging',
'        Charlestown, MA 02129, USA',
'',
'== Description of the workflow ==',
'',
' 1. when the GUI pops up, it will display the mesh and the points,',
'    you can rotate both plots so that you can identify the matching ',
'    features',
' 2. switch on "Select" mode, then, click on a land-mark point on the point',
'    plot, when a data-tip shows up, click "Add Selected" button',
' 3. click on the corresponding position on the mesh, and click',
'    "Add Selected"      ',
' 4. repeat the above for at least 4 point pairs (you can select more);',
'    if you want to change views, switch off "Select" box and rotate;',
'    after rotation, switch on "Select" box again',
' 5. click "Initialize": this will create the initial mapping using the',
'    selected point pairs',
' 6. click "Optimize": this will fit the surface with the whole point cloud',
' 7. click "Proj2Mesh": this will project the fitted point clouds onto the',
'    mesh',
' 8. you can quit the GUI by hit "Close", your results will be saved to reg',
' 9. close the window '};

helpdlg(helpmsg);

author: fangq


text:
function [b,S] = bin_packing(a,V,varargin)
  % BIN_PACKING Given a list of n items of varying sizes (a), pack them into the
  % smallest number (b) of bins of equal size (V).
  %
  % Inputs:
  %   a  #a list of item sizes
  %   V  scalar size of bins
  % Outputs:
  %   b  minimal number of bins
  %   S  #a list of indices into 1:b, assigning each item to a bin
  %
  % Example:
  %   a = ceil(rand(100,1)*10);
  %   V = 20;
  %   [b,S] = bin_packing(a,V);
  %   barh(full(sparse(S,1:numel(a),a)),'stacked');
  %

  max_iters = 5;
  params_to_variables = containers.Map( ...
    {'MaxIters'},{'max_iters'});
  v = 1;
  while v <= numel(varargin)
    param_name = varargin{v};
    if isKey(params_to_variables,param_name)
      assert(v+1<=numel(varargin));
      v = v+1;
      % Trick: use feval on anonymous function to use assignin to this workspace
      feval(@()assignin('caller',params_to_variables(param_name),varargin{v}));
    else
      error('Unsupported parameter: %s',varargin{v});
    end
    v=v+1;
  end

  assert(all(a<V),'No item can be larger than single bin');
  n = numel(a);
  b = inf;
  S = zeros(n,1);
  % Try five times (first time is sort)
  for iter = 1:max_iters
    switch iter
    case 1
      P = (1:numel(a))';
    case 2
      [~,P] = sort(a,'descend');
    otherwise
      P = randperm(n);
    end
    ar = a(P);
    % start with no assignments
    Sr = zeros(n,1);
    % running remaining bin volumes
    Br = repmat(V,n,1);

    % naive O(n²) first fit:
    for i = 1:n
      % O(n) search to find first available bin
      j = find(Br>ar(i),1,'first');
      if j>=b
        break;
      end
      Sr(i) = j;
      Br(j) = Br(j) - ar(i);
    end
    if j >= b
      continue;
    end

    br = find(Br==V,1,'first')-1;
    if br < b
      b = br;
      B = Br;
      S(P) = Sr;
    end
  end

end

author: alecjacobson

repo: gptoolbox


text:
clear

load('results/results_clnf_cross-data_general.mat');
load('results/menpo_labels.mat');
[clnf_error, frontal_ids] = compute_error_menpo_small( labels,  experiment.shapes);
clnf_error_frontal = clnf_error(frontal_ids);
clnf_error_profile = clnf_error(~frontal_ids);

load('results/menpo_train_chehra.mat');
[drmf_error, frontal_ids] = compute_error_menpo_small( labels,  shapes);
drmf_error_frontal = drmf_error(frontal_ids);
drmf_error_profile = drmf_error(~frontal_ids);

load('results/menpo_train_sdm.mat');
for i = 1:numel(experiment.shapes)
    shapes{i} = shapes{i}+0.5;
end
[sdm_error, frontal_ids] = compute_error_menpo_small( labels,  shapes);
sdm_error_frontal = sdm_error(frontal_ids);
sdm_error_profile = sdm_error(~frontal_ids);

load('results/Menpo_train_pocr.mat');
[pocr_error, frontal_ids] = compute_error_menpo_small( labels,  experiments.shapes);
pocr_error_frontal = pocr_error(frontal_ids);
pocr_error_profile = pocr_error(~frontal_ids);

load('results/results_ceclm_cross-data.mat');
[ceclm_error, frontal_ids] = compute_error_menpo_small( labels,  experiment.shapes);
ceclm_error_frontal = ceclm_error(frontal_ids);
ceclm_error_profile = ceclm_error(~frontal_ids);

load('results/tcdcn_menpo.mat');
for i = 1:numel(shapes)
    shapes{i} = shapes{i};
end

[tcdcn_error, frontal_ids] = compute_error_menpo_small(labels, shapes);
tcdcn_error_frontal = tcdcn_error(frontal_ids);
tcdcn_error_profile = tcdcn_error(~frontal_ids);

load('results/CFAN_menpo_train.mat');
for i = 1:numel(shapes)
    shapes{i} = shapes{i}-0.5;
end

[cfan_error, frontal_ids] = compute_error_menpo_small(labels, shapes);
cfan_error_frontal = cfan_error(frontal_ids);
cfan_error_profile = cfan_error(~frontal_ids);

load('results/menpo_train_3DDFA.mat');
for i = 1:numel(shapes)
    shapes{i} = shapes{i}-0.5;
end

[error_3ddfa, frontal_ids] = compute_error_menpo_small(labels, shapes);
error_3ddfa_frontal = error_3ddfa(frontal_ids);
error_3ddfa_profile = error_3ddfa(~frontal_ids);

load('results/Menpo-CFSS_train.mat');
shapes = cell(size(estimatedPoseFull,1),1);

for i = 1:numel(shapes)
    shape = cat(2, estimatedPoseFull(i,1:68)', estimatedPoseFull(i,69:end)');
    shapes{i} = shape-0.5;
end

[cfss_error, frontal_ids] = compute_error_menpo_small(labels, shapes);
cfss_error_frontal = cfss_error(frontal_ids);
cfss_error_profile = cfss_error(~frontal_ids);

author: TadasBaltrusaitis

repo: OpenFace


text:
function mrmMakeMovie(id,rotBegin,rotEnd)
%
%
% Function to make a movie of brain rotating left to ventral view.  This
% should really be a more general function, but this is a start.
%
%   written by amr Jun 2010
%

if ~exist('id','var'), id = 500; end
if ~exist('rotBegin','var')
    mrmRotateCamera(id,'left')
    [rotBegin,zoom1] = mrmGetRotation(id);
end

if ~exist('rotEnd','var')
    mrmRotateCamera(id,'bottom')
    [rotEnd,zoom2] = mrmGetRotation(id);
end

%keyboard
% Set the background color
p.color=[0,0,0,1];
%p.color=[.3,.3,.3,1];

host = 'localhost';
%id = 174;
mrMesh(host, id, 'background', p);

% % rotate left to bottom
rotY = [rotBegin(2):-pi/128:rotEnd(2)];
rotZ = [rotBegin(3):pi/128:rotEnd(3)];
rotX = ones(1,length(rotY))*pi;

% rotate bottom to left
% rotY = [rotEnd(2):pi/64:rotBegin(2)];
% rotZ = [rotEnd(3):-pi/64:rotBegin(3)];
% rotX = ones(1,length(rotY))*pi;

n = length(rotX);
%pitch = linspace(.5*pi, 1.5*pi, n);
%pitch = -pi/2.5;
%zoom = 1;
movDir = '/biac3/wandell7/data/Words/Meshes/';
movFile = 'mrmRotateVentralToLeft.avi';

%mkdir('/tmp', 'mrmMovie');
clear M;
%for(ii=1:length(pitch))
f.filename = 'nosave';
for(ii=1:length(rotX))
    mrmRotateCamera(id, [rotX(ii) rotY(ii) rotZ(ii)], zoom1);

    [id,stat,res] = mrMesh(host, id, 'screenshot', f);
    M((1-1)*length(rotX)+ii) = im2frame(permute(res.rgb, [2,1,3])./255);

    %fname = sprintf('%c%0.2d.png', ltr(ii), jj);
    %fname = fullfile(movDir, fname);
    %imwrite(permute(res.rgb,[2,1,3])./255, fname);
end
%end
%figure; movie(M,-3)
movie2avi(M,fullfile(movDir,movFile)); %'/tmp/mrmMovieAllFrames.avi');

return

function [rot,zoom] = mrmGetRotation(id)
% Get rotation and zoom
p.actor=0; p.get_all=1;
[id,stat,r] = mrMesh('localhost', id, 'get', p);
zoom = diag(chol(r.rotation'*r.rotation))';
rotMat = r.rotation/diag(zoom);
% Note- there may be slight rounding errors allowing the inputs to
% asin/atan go outside of the range (-1,1). May want to clip those.
rot(2) = asin(rotMat(1,3));
if (abs(rot(2))-pi/2).^2 < 1e-9,
	rot(1) = 0;
	rot(3) = atan2(-rotMat(2,1), -rotMat(3,1)/rotMat(1,3));
else
	c      = cos(rot(2));
	rot(1) = atan2(rotMat(2,3)/c, rotMat(3,3)/c);
	rot(3) = atan2(rotMat(1,2)/c, rotMat(1,1)/c);
end
rot(1) = -rot(1); % flipped OpenGL Y-axis.
rot(3) = -rot(3);
fprintf('rot=[%0.6f %0.6f %0.6f];\nzoom=[%0.3f %0.3f %0.3f];\nfrustum=[%0.6f %0.6f %0.6f %0.6f];\n',rot,zoom,r.frustum);

return

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function varargout = img2mesh(varargin)
%
%  Format:
%      newworkspace = img2mesh or imgmesh(workspace)
%
%  A GUI for Iso2Mesh for streamlined mesh data processing
%
%  Author: Qianqian Fang <q.fang at neu.edu>
%
%  Input:
%        workspace (optional): a struct containing the below fields
%           .graph: a digraph object containing the i2m workspace data
%  Output:
%        newworkspace (optional): the updated workspace, with the same
%        subfields as the input.
%
%   If a user supplys an output variable, the GUI will not return until
%   the user closes the window; if a user does not provide any output,
%   the call will return immediately.
%
%   Please find more information at http://iso2mesh.sf.net/
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

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

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

% --- Executes just before i2m is made visible.
function i2m_OpeningFcn(hObject, eventdata, handles, varargin)

cm = uicontextmenu;
uimenu(cm,'Label','Plot','CallBack',{@processdata,handles});
uimenu(cm,'Label','Rename','CallBack',{@processdata,handles});
uimenu(cm,'Label','Delete','CallBack',{@processdata,handles});

mimeshing = uimenu(cm,'Label','Meshing');
miv2s = uimenu(mimeshing,'Label','Volume to surface','CallBack',{@processdata,handles});
uimenu(mimeshing,'Label','Volume to mesh','CallBack',{@processdata,handles});
uimenu(mimeshing,'Label','Surface to mesh','CallBack',{@processdata,handles});
uimenu(mimeshing,'Label','Surface to volume','CallBack',{@processdata,handles});
uimenu(mimeshing,'Label','Close and fill volume','CallBack',{@processdata,handles});
uimenu(mimeshing,'Label','Extract surface','CallBack',{@processdata,handles});
uimenu(mimeshing,'Label','Tessellate surface','CallBack',{@processdata,handles});

mirepair = uimenu(cm,'Label','Surface repair');
uimenu(mirepair,'Label','Clean surface','CallBack',{@processdata,handles});
uimenu(mirepair,'Label','Repair surface','CallBack',{@processdata,handles});
uimenu(mirepair,'Label','Smooth surface','CallBack',{@processdata,handles});
uimenu(mirepair,'Label','Simplify surface','CallBack',{@processdata,handles});
uimenu(mirepair,'Label','Remesh surface','CallBack',{@processdata,handles});
uimenu(mirepair,'Label','Simplify surface','CallBack',{@processdata,handles});
uimenu(mirepair,'Label','Reorient mesh elements','CallBack',{@processdata,handles});

mibool = uimenu(cm,'Label','Surface boolean');
uimenu(mibool,'Label','Or','CallBack',{@processdata,handles});
uimenu(mibool,'Label','And','CallBack',{@processdata,handles});
uimenu(mibool,'Label','All','CallBack',{@processdata,handles});
uimenu(mibool,'Label','Diff','CallBack',{@processdata,handles});
uimenu(mibool,'Label','First','CallBack',{@processdata,handles});
uimenu(mibool,'Label','Second','CallBack',{@processdata,handles});

mireport = uimenu(cm,'Label','Report');
uimenu(mireport,'Label','Containing data','CallBack',{@processdata,handles});
uimenu(mireport,'Label','Mesh quality histogram','CallBack',{@processdata,handles});
uimenu(mireport,'Label','Element volume histogram','CallBack',{@processdata,handles});
uimenu(mireport,'Label','Total volume','CallBack',{@processdata,handles});

mirefresh=uimenu(cm,'Label','Refresh','CallBack',{@miRefresh_Callback,handles});
uimenu(cm,'Label','Save as','CallBack',{@processdata,handles});

miv2s.Separator='on';
mimeshing.Separator='on';
mirefresh.Separator='on';

root=get(handles.fgI2M,'userdata');

if(isempty(root))
    root=struct('graph',digraph,'menu',cm);
end
set(handles.fgI2M,'userdata',root);
set(handles.axFlow,'position',[0 0 1 1]);
set(handles.axPreview,'position',[0 0 1 1]);

set(handles.fgI2M,'UIContextMenu',handles.meCreate);

axis(handles.axFlow,'off');
axis(handles.axPreview,'off');

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

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes i2m wait for user response (see UIRESUME)
% uiwait(handles.fgI2M);

function processdata(source,callbackdata,handles)

try
obj=get(handles.fgI2M,'currentobject');
if(strcmp(class(obj),'matlab.graphics.chart.primitive.GraphPlot')==0)
    if(~isempty(obj.UserData) && strcmp(class(obj.UserData),'matlab.graphics.chart.primitive.GraphPlot'))
        obj=obj.UserData;
    else
        return
    end
end
root=get(handles.fgI2M,'userdata');
pos=get(handles.axFlow,'currentpoint');
[nodedata,nodetype,nodeid]=getnodeat(root,obj,pos);

newtype=dummytype;
prefix='x';
switch source.Label
    case 'Volume to surface'
        if(isstruct(nodetype) && isfield(nodetype,'hasvol') && nodetype.hasvol)
            [newdata,newtype]=v2sgui(nodedata);
            prefix='Vol2Surf';
        else
            error('no volume data found');
        end
    case 'Volume to mesh'
        if(nodetype.hasvol)
            [newdata,newtype]=v2mgui(nodedata);
            prefix='Vol2Mesh';
        else
            error('no volume data found');
        end
    case 'Surface to mesh'
        if(nodetype.hasnode && nodetype.hasface)
            [newdata,newtype]=s2mgui(nodedata);
            prefix='Surf2Mesh';
        else
            error('no surface data found');
        end
    case 'Surface to volume'
        if(nodetype.hasnode && nodetype.hasface)
            ndiv = inputdlg('Division number along the shortest dimension:',...
                'surf2vol - rasterizing a surface mesh',1,{'50'});
            if(isempty(ndiv))
                return;
            end
            newdata.vol=s2v(nodedata.node,nodedata.face,str2num(ndiv{1}));
            newtype.hasvol=1;
            prefix='Surf2Vol';
        else
            error('no surface data found');
        end
    case 'Close and fill volume'
        if(nodetype.hasvol)
            rad=inputdlg('maximum gap length in voxel (scalar)','Close and fill a volume',1,{'1'});
            newdata.vol=fillholes3d(nodedata.vol,str2num(rad{1}));
            newtype=nodetype;
            prefix='FillVol';
        else
            error('no volume data found');
        end
    case 'Extract surface'
        if(isstruct(nodetype) && isfield(nodetype,'hasvol') && nodetype.hasvol)
            [newdata.node,newdata.face]=binsurface(nodedata.vol);
            prefix='BinSurf';
        elseif(nodetype.hasnode && nodetype.haselem)
            newdata.node=nodedata.node;
            newdata.face=volface(nodedata.elem(:,1:min(4,size(nodedata.elem,2))));
            [newdata.node,newdata.face]=removeisolatednode(newdata.node,newdata.face);
            prefix='VolSurf';
        elseif(nodetype.hasnode && nodetype.hasface)
            [newdata.node,newdata.face]=deal(nodedata.node,nodedata.face);
            prefix='CopySurf';
        end
        newtype.hasnode=1;
        newtype.hasface=1;
    case 'Clean surface'
        if(nodetype.hasnode && nodetype.hasface)
            [newdata.node,newdata.face]=meshcheckrepair(nodedata.node,nodedata.face,'deep');
            newtype=nodetype;
            prefix='CleanSurf';
        else
            error('no surface data found');
        end
    case 'Repair surface'
        if(nodetype.hasnode && nodetype.hasface)
            [newdata.node,newdata.face]=meshcheckrepair(nodedata.node,nodedata.face,'meshfix');
            newtype=nodetype;
            prefix='RepairSurf';
        else
            error('no surface data found');
        end
    case 'Smooth surface'
        if(nodetype.hasnode && nodetype.hasface)
            res = inputdlg({'Method (laplacian,laplacianhc,lowpass):','Iteration (integer):','Alpha (scalar):'},...
                'sms - smoothing a surface mesh',[1,1,1],{'lowpass','20','0.5'});
            if(isempty(res))
                return;
            end
            newdata=nodedata;
            newtype=nodetype;
            newdata.node=sms(nodedata.node,nodedata.face,str2num(res{2}),str2num(res{3}),res{1});
            prefix='SmoothSurf';
        else
            error('no surface data found');
        end
    case 'Simplify surface'
        if(nodetype.hasnode && nodetype.hasface)
            res = inputdlg('Percentage of edges to keep (0-1):',...
                'Simplify mesh',[1],{'1'});
            if(isempty(res))
                return;
            end
            [newdata.node,newdata.face]=meshresample(nodedata.node,nodedata.face,str2num(res{1}));
            newtype=nodetype;
            prefix='SimplifySurf';
        else
            error('no surface data found');
        end
    case 'Remesh surface'
        if(nodetype.hasnode && nodetype.hasface)
            res = inputdlg({'Rasterization voxel size (scalar):','Max gap to close (scalar):','Max surface element radis (scalar):'},...
                'Remesh surface',[1,1,1],{'1','20','3'});
            if(isempty(res))
                return;
            end
            opt.gridsize=str2num(res{1});
            opt.closesize=str2num(res{2});
            opt.elemsize=str2num(res{3});
            [newdata.node, newdata.face]=remeshsurf(nodedata.node,nodedata.face,opt);
            newtype.hasnode=1;
            newtype.hasface=1;
            prefix='RemeshSurf';
        else
            error('no surface data found');
        end
    case 'Tessellate surface'
        if(nodetype.hasnode && nodetype.hasface)
            [newdata.node,newdata.elem]=fillsurf(nodedata.node,nodedata.face);
            newdata.face=volface(newdata.elem(:,1:min(4,size(newdata.elem,2))));
            newtype=nodetype;
            newtype.haselem=1;
            prefix='TessMesh';
        else
            error('no surface data found');
        end
    case 'Reorient mesh elements'
        if(nodetype.hasnode && nodetype.haselem)
            [newdata.node,newdata.elem]=meshreorient(nodedata.node,nodedata.elem);
            newtype=nodetype;
            prefix='ReorientMesh';
        elseif(nodetype.hasnode && nodetype.hasface)
            [newdata.node,newdata.face]=surfreorient(nodedata.node,nodedata.face);
            newtype=nodetype;
            prefix='ReorientSurf';
        else
            error('no surface or tetrahedral mesh data found');
        end
    case 'Mesh quality histogram'
        if(nodetype.hasnode && nodetype.haselem)
            quality=meshquality(nodedata.node,nodedata.elem);
        elseif(nodetype.hasnode && nodetype.hasface)
            quality=meshquality(nodedata.node,nodedata.face);
        end
        if(exist('quality','var'))
            figure;
            hist(quality,50);
        end
        return;
    case 'Element volume histogram'
        if(nodetype.hasnode && nodetype.haselem)
            evol=elemvolume(nodedata.node,nodedata.elem);
        elseif(nodetype.hasnode && nodetype.hasface)
            evol=elemvolume(nodedata.node,nodedata.face);
        end
        if(exist('evol','var'))
            figure;
            hist(evol,50);
        end
        return;
    case 'Total volume'
        if(nodetype.hasnode && nodetype.haselem)
            evol=elemvolume(nodedata.node,nodedata.elem);
        elseif(nodetype.hasnode && nodetype.hasface)
            [no,el]=fillsurf(nodedata.node,nodedata.face);
            evol=elemvolume(no,el);
        end
        if(exist('evol','var'))
            msgbox(sprintf('Total volume is %f cubic voxel',sum(evol)),'Total volume');
        end
        return;
    case 'Containing data'
        msg='';
        if(nodetype.hasnode)
            msg=[msg sprintf('\nContaining %d nodes (%d columns)',size(nodedata.node,1),size(nodedata.node,2))];
        end
        if(nodetype.hasface)
            msg=[msg sprintf('\nContaining %d triangles (%d columns)',size(nodedata.face,1),size(nodedata.face,2))];
        end
        if(nodetype.haselem)
            msg=[msg sprintf('\nContaining %d tetrehedra (%d columns)',size(nodedata.elem,1),size(nodedata.elem,2))];
        end
        if(nodetype.hasvol)
            msg=[msg sprintf('\nContaining [%d x %d x %d ] volume',size(nodedata.vol,1),size(nodedata.vol,2),size(nodedata.vol,3))];
        end
        msgbox(msg,'Mesh data report');
        return;
    case 'Plot'
        if(isstruct(nodetype) && isfield(nodetype,'hasnode') && nodetype.hasnode)
            if(isfield(nodetype,'haselem') && nodetype.haselem)
                figure('keypressfcn',@plotfigevent,'userdata',struct('node',nodedata.node,'face',[],'elem',nodedata.elem));
                plotmesh(nodedata.node,[],nodedata.elem);
            else
                figure('keypressfcn',@plotfigevent,'userdata',struct('node',nodedata.node,'elem',[],'face',nodedata.face));
                plotmesh(nodedata.node,nodedata.face);
            end
        else
            figure('keypressfcn',@plotfigevent,'userdata',struct('vol',nodedata.vol));
            hs=slice(double(nodedata.vol),[],[ceil(size(nodedata.vol,2)*0.5)],ceil(size(nodedata.vol,3)*0.5));
            set(hs,'linestyle','none');
        end
    case {'Or','And','Diff','All','First','Second'}
        if(isstruct(nodetype) && isfield(nodetype,'hasnode'))
            if((nodetype.hasface || nodetype.haselem) && nodetype.hasnode)
                pt=ginput(1);
                [nodedata2,nodetype2,nodeid2]=getnodeat(root,obj,pt);
                if(~nodetype2.hasnode || ~nodetype2.hasface)
                    if(nodetype2.hasnode && nodetype2.haselem)
                        nodedata2.face=volface(nodedata2.elem);
                    else
                        error('Second operand does not contain a surface');
                    end
                end
                op=source.Label;
                if(strcmp(op,'Intersect'))
                    op='inter';
                end
                if(~nodetype.hasface)
                    nodedata.face=volface(nodedata.elem);
                end
                [newdata.node,newdata.face]=surfboolean(nodedata.node,nodedata.face,lower(op),nodedata2.node,nodedata2.face(:,[1 3 2]));
                newtype.hasnode=1;
                newtype.hasface=1;
                prefix=source.Label;
            else
                warndlg('Selected node does not contain a surface mesh');
            end
        end
    case 'Delete'
        button = questdlg('Are you sure to delete the selected node?', 'Confirm','No');
        if strcmpi(button, 'No') ||  strcmpi(button, 'Cancel')
	        return;
        end
        root.graph=rmnode(root.graph,root.graph.Nodes.Name{nodeid});
        updategraph(root,handles);
    case 'Rename'
        newname = inputdlg('Define a new name:',...
                'Rename',1,{root.graph.Nodes.Name{nodeid}});
        if(isempty(newname))
            return;
        end
        if(isempty(newname{1}) || ~isempty(cell2mat(regexp( root.graph.Nodes.Name, ['^' newname{1} '$']))))
            error('empty or duplicated node name');
        end
        root.graph.Nodes.Name{nodeid}=newname{1};
        updategraph(root,handles);
    case 'Save as'
        if(~nodetype.haselem && ~nodetype.hasnode)
            error('selected data does not have a mesh');
            return;
        end
        filter = {'*.jmesh';'*.*'};
        [file, path] = uiputfile(filter,'Export mesh');
        if ~isequal(file,0) && ~isequal(path,0)
            if(nodetype.haselem)
                savejmesh(nodedata.node,nodedata.face,nodedata.elem,fullfile(path,file));
            else
                savejmesh(nodedata.node,nodedata.face,fullfile(path,file));
            end
        end
end

if(exist('newdata','var') && exist('newtype','var'))
    cla(handles.axPreview);
    cla(handles.axFlow);
    newdata.preview=getpreview(newdata,newtype,[400,400]);
    [newkey,root.graph]=addnodewithdata(handles,newdata,newtype,prefix);
    root.graph=addedge(root.graph,{root.graph.Nodes.Name{nodeid}},{newkey});
    if(strcmp(source.Parent.Type, 'uimenu') && strcmp(source.Parent.Label,'Surface boolean'))
        root.graph=addedge(root.graph,{root.graph.Nodes.Name{nodeid2}},{newkey});
    end
    updategraph(root,handles);
end

catch ME
    msg=sprintf('Error: \n%s\n',ME.message);
    for e=1:length(ME.stack)
       msg=sprintf('%s\nFile: %s\nFunction: %s\nLine: %d\n\n',msg,ME.stack(e).file,ME.stack(e).name,ME.stack(e).line);
    end
	uiwait(warndlg(msg,'I2M ERROR'));
    updategraph(root,handles);
end

function plotfigevent(hobject,event)
data=get(hobject,'userdata');
plotpos=[];
if(isfield(data,'plotpos'))
    plotpos=data.plotpos;
end
if(isempty(plotpos))
    switch event.Key
        case {'rightarrow','uparrow'}
            plotpos=0;
        case {'leftarrow','downarrow'}
            plotpos=9;
        otherwise
    end
end
if(isfield(data,'node'))
    pmax=max(data.node);
    pmin=min(data.node);
    cla;
    switch event.Key
        case 'rightarrow'
            plotpos=min(plotpos+1,9);
            plotmesh(data.node,data.face,data.elem,sprintf('x>%f',(pmax(1)-pmin(1))*plotpos*0.1+pmin(1)))
        case 'leftarrow'
            plotpos=max(plotpos-1,0);
            plotmesh(data.node,data.face,data.elem,sprintf('x>%f',(pmax(1)-pmin(1))*plotpos*0.1+pmin(1)))
        case 'uparrow'
            plotpos=min(plotpos+1,9);
            plotmesh(data.node,data.face,data.elem,sprintf('y>%f',(pmax(2)-pmin(2))*plotpos*0.1+pmin(2)))
        case 'downarrow'
            plotpos=max(plotpos-1,0);
            plotmesh(data.node,data.face,data.elem,sprintf('y>%f',(pmax(2)-pmin(2))*plotpos*0.1+pmin(2)))
        otherwise
    end
end
data.plotpos=plotpos;
set(hobject,'userdata',data);
%----------------------------------------------------------------
function [nodedata,nodetype,nodeid]=getnodeat(root,obj,pos)
nodedist=[obj.XData(:)-pos(1,1) obj.YData(:)-pos(1,2)];
nodedist=sum(nodedist.*nodedist,2);
[mindist, nodeid]=min(nodedist);
nodedata=root.graph.Nodes.Data{nodeid};
nodetype=root.graph.Nodes.Type{nodeid};

%----------------------------------------------------------------
function mytype=dummytype
mytype.hasnode=0;
mytype.hasface=0;
mytype.haselem=0;
mytype.hasvol=0;

%----------------------------------------------------------------
function [newdata, newtype]=v2sgui(data)
prompt = {'Threshold (scalar or array):',...
    'Surface element radius bound (scalar):',...
    'Surface element distance bound (scalar)','Method: (cgalsurf,cgalmesh,simplify)'};
title = 'vol2surf - extracting surface mesh from volume';
dims = [1 1 1 1];
definput = {'0.5','5','1','cgalsurf'};
res = inputdlg(prompt,title,dims,definput);
newdata=[];
newtype=dummytype;
if(isempty(res))
    return;
end

opt=struct('radbound',str2num(res{2}),'distbound',str2num(res{3}));
[newdata.node,newdata.face]=v2s(data.vol,eval(res{1}),opt, res{4});
newtype.hasnode=1;
newtype.hasface=1;
%----------------------------------------------------------------
function [newdata, newtype]=v2mgui(data)
prompt = {'Threshold (scalar or []):',...
    'Surface element radius bound (scalar):',...
    'Surface element distance bound (scalar)',...
    'Max element volume (scalar):',...
    'Method (cgalsurf,cgalmesh,simplify):'};
title = 'vol2mesh - extracting tet mesh from volume';
dims = [1 1 1 1 1];
definput = {'[]','5','1','30','cgalmesh'};
res = inputdlg(prompt,title,dims,definput);
newdata=[];
newtype=dummytype;
if(isempty(res))
    return;
end

opt=struct('radbound',str2num(res{2}),'distbound',str2num(res{3}));
[newdata.node,newdata.elem,newdata.face]=v2m(data.vol,eval(res{1}),...
   opt, res{4});
newtype.hasnode=1;
newtype.hasface=1;
newtype.haselem=1;

%----------------------------------------------------------------
function img=getpreview(nodedata,nodetype,imsize)
hfpreview=figure('visible', 'off');
ax=axes('parent',hfpreview,'Units','pixels','position',[1, 1, imsize(1), imsize(2)]);
if(isfield(nodetype,'haselem') && nodetype.haselem)
    plotmesh(nodedata.node,[],nodedata.elem,'linestyle',':','edgealpha',0.3,'parent',ax);
elseif(isfield(nodetype,'hasface') && nodetype.hasface)
    plotmesh(nodedata.node,nodedata.face,'linestyle','-','parent',ax);
elseif(isfield(nodetype,'hasvol') && nodetype.hasvol)
    hs=slice(double(nodedata.vol),[],[ceil(size(nodedata.vol,2)*0.5)],ceil(size(nodedata.vol,3)*0.5),'parent',ax);
    set(hs,'linestyle','none');
elseif(isfield(nodetype,'hasnode') && nodetype.hasnode)
    plotmesh(nodedata.node,'.','parent',ax);
end
set(ax,'color','none')
axis(ax,'equal');
axis(ax,'off');
img=getframe(gca);
img=flipud(img.cdata);
delete(ax);
close(hfpreview);

%----------------------------------------------------------------
function [newdata, newtype]=s2mgui(data)
prompt = {'Simplification ratio (%edges to keep, 0-1):',...
    'Max element volume (scalar):',...
    'Method (tetgen,tetgen1.5,cgalpoly):',...
    'Region seeds (N x 3 array):',...
    'Hole seeds (N x 3 array):'};
title = 'surf2mesh - creating tet mesh from surfaces';
dims = [1 1 1 1 1];
definput = {'1','30','tetgen','[]','[]'};
res = inputdlg(prompt,title,dims,definput);
newdata=[];
newtype=dummytype;
if(isempty(res))
    return;
end

[newdata.node,newdata.elem,newdata.face]=...
   s2m(data.node,data.face,str2num(res{1}),...
     str2num(res{2}),res{3},eval(res{4}),eval(res{5}));
newtype.hasnode=1;
newtype.hasface=1;
newtype.haselem=1;

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

% Get default command line output from handles structure
handles.output=get(handles.fgI2M,'userdata');
varargout{1} = handles.output;

% --------------------------------------------------------------------
function miWeb_Callback(hObject, eventdata, handles)
web('http://iso2mesh.sourceforge.net');

% --------------------------------------------------------------------
function miDoc_Callback(hObject, eventdata, handles)
web('http://iso2mesh.sourceforge.net/cgi-bin/index.cgi?Doc');

% --------------------------------------------------------------------
function miAbout_Callback(hObject, eventdata, handles)
helpmsg={
'\bf\fontsize{12}I2M: An Integrated GUI for Iso2Mesh Meshing Toolbox\rm\fontsize{10}',
'',
'Copyright (c) 2018 Qianqian Fang <q.fang at neu.edu>',
''
'Computational Optics&Translational Imaging Lab (http://fanglab.org)',
'Department of Bioengineering',
'Northeastern University',
'360 Huntington Ave, Boston, MA 02115, USA',
'',
'URL:    http://iso2mesh.sourceforge.net',
''};

opt.Interpreter = 'tex';
opt.WindowStyle = 'modal';

msgbox(helpmsg,'About','help',opt);

% --------------------------------------------------------------------
function miSphere_Callback(hObject, eventdata, handles)
prompt = {'Center:','Radius (scalar):',...
    'Surface element radius bound (scalar)',...
    'Max element volume (scalar, 0 - only create surface):'};
title = 'Create Mesh';
dims = [1 1 1 1];
definput = {'[0 0 0]','50','6','20'};
res = inputdlg(prompt,title,dims,definput);
if(isempty(res))
    return;
end
newtype=dummytype;
opt=str2num(res{3});
if(str2num(res{4})==0)
   [newdata.node,newdata.face]=meshasphere(eval(res{1}),...
     str2num(res{2}),opt);
else
   [newdata.node,newdata.face,newdata.elem]=meshasphere(eval(res{1}),...
     str2num(res{2}),opt, str2num(res{4}));
   newtype.haselem=1;
end
newtype.hasnode=1;
newtype.hasface=1;

if(exist('newdata','var') && exist('newtype','var'))
    newkey=addnodewithdata(handles,newdata,newtype,'Sphere');
end

% --------------------------------------------------------------------
function miBox_Callback(hObject, eventdata, handles)
prompt = {'Diagonal end point 1 (1x3 vector):','Diagonal end point 2 (1x3 vector):',...
    'Surface element radius bound (scalar)',...
    'Max element volume (scalar, 0 - only create surface):'};
title = 'Create Mesh';
dims = [1 1 1 1];
definput = {'[0 0 0]','[100 60 30]','6','30'};
res = inputdlg(prompt,title,dims,definput);
if(isempty(res))
    return;
end
newtype=dummytype;
opt=str2num(res{3});
if(str2num(res{4})==0)
   [newdata.node,newdata.face]=meshabox(eval(res{1}),eval(res{2}),opt);
else
   [newdata.node,newdata.face,newdata.elem]=meshabox(eval(res{1}),...
     eval(res{2}),opt, str2num(res{4}));
   newtype.haselem=1;
end
newtype.hasnode=1;
newtype.hasface=1;

if(exist('newdata','var') && exist('newtype','var'))
    newkey=addnodewithdata(handles,newdata,newtype,'Box');
end

% --------------------------------------------------------------------
function miCylinder_Callback(hObject, eventdata, handles)
prompt = {'Axis end-point 1','Axis end-point 2','Radius (scalar):',...
    'Surface element radius bound (scalar)',...
    'Max element volume (scalar, 0 - only create surface):',...
    'Circle division:'};
title = 'Create Mesh';
dims = [1 1 1 1 1 1];
definput = {'[0 0 0]','[0 0 50]','10','3','20','20'};
res = inputdlg(prompt,title,dims,definput);
if(isempty(res))
    return;
end
newtype=dummytype;
opt=str2num(res{4});
maxvol=str2num(res{5});

if(maxvol==0)
   [newdata.node,newdata.face]=meshacylinder(eval(res{1}),eval(res{2}),...
     str2num(res{3}),opt);
else
   [newdata.node,newdata.face,newdata.elem]=meshacylinder(eval(res{1}),eval(res{2}),...
     str2num(res{3}),opt,maxvol);
   newtype.haselem=1;
end
newtype.hasnode=1;
newtype.hasface=1;

if(exist('newdata','var') && exist('newtype','var'))
    newkey=addnodewithdata(handles,newdata,newtype,'Cyl');
end

% --------------------------------------------------------------------
function miLoadVol_Callback(hObject, eventdata, handles)

nodedata=struct;
nodetype=dummytype;
filters={'*.nii;*.hdr;*.img;*.tif;*.tiff;*.inr;*.bin;*.ubj','3D volume file (*.nii;*.hdr;*.img;*.tif;*.tiff;*.inr;*.bin;*.ubj)';...
    '*.nii','Nifti file (*.nii)';...
    '*.hdr;*.img','Analyze 7.5 file (*.hdr;*.img)';...
    '*.tif;*.tiff','Multipage TIFF file (*.tif)';...
    '*.inr','INR image (*.inr)';...
    '*.bin','Binary file (*.bin)';...
    '*.ubj','Universal JSON (*.ubj)';...
    '*.*','All (*.*)'};
[file,path,idx] = uigetfile(filters);
if isequal(file,0)
   return;
else
   if(regexp(file,'\.[Nn][Ii][Ii]$'))
       im=readnifti(fullfile(path,file));
       nodedata.vol=im.img;
       nodetype.hasvol=1;
   elseif(regexp(file,'(\.[Hh][Dd][Rr]$|\.[Ii][Mm][Gg]$)'))
       im=readnifti(fullfile(path,file));
       nodedata.vol=im.img;
       nodetype.hasvol=1;
   elseif(regexp(file,'\.[Tt][Ii][Ff][Ff]*$'))
       nodedata.vol=readmptiff(fullfile(path,file));
       nodetype.hasvol=1;
   elseif(regexp(file,'\.[Ii][Nn][Rr]$'))
       nodedata.vol=readinr(fullfile(path,file));
       nodetype.hasvol=1;
   elseif(regexp(file,'\.[Bb][Ii][Nn]$'))
       prompt = {'Dimension (1x3 vector):',...
            'Datatype (short,float,double,integer,...):'};
       title = 'Load generic binary file';
       dims = [1 1];
       definput = {'[]','short'};
       [res,isok] = inputdlg(prompt,title,dims,definput);
       if(isok==0)
           return;
       end
       nodedata.vol=loadmc2(fullfile(path,file),eval(res{1}),res{2});
       nodetype.hasvol=1;
   elseif(regexp(file,'\.[Uu][Bb][Jj]$'))
       nodedata=loadbj(fullfile(path,file));
       if(isstruct(nodedata) && isfield(nodedata,'vol'))
           nodetype.hasvol=1;
       end
   end
end

if(exist('nodedata','var'))
    nodetype=getnodetype(nodedata);
    if(nodetype.hasvol)
        addnodewithdata(handles,nodedata,nodetype,'Vol');
    end
else
    warndlg('no valid mesh data found','Warning');
end

% --------------------------------------------------------------------
function miLoadMesh_Callback(hObject, eventdata, handles)

nodedata=struct;
nodetype=dummytype;
filters={'*.jmesh;*.off;*.medit;*.smf;*.json','3D Mesh files (*.jmesh;*.off;*.medit;*.smf;*.json)';...
    '*.jmesh','JSON mesh (*.jmesh)';...
    '*.off','OFF file (*.off)';...
    '*.medit','Medit file (*.medit)';...
    '*.ele','Tetgen element mesh file (*.ele)';...
    '*.json','JSON file (*.json)';'*.*','All (*.*)'};
[file,path,idx] = uigetfile(filters);
if isequal(file,0)
   return;
else
   if(regexp(file,'\.[Oo][Ff][Ff]$'))
       [nodedata.node, nodedata.face]=readoff(fullfile(path,file));
   elseif(regexp(file,'\.[Mm][Ee][Dd][Ii][Tt]$'))
       [nodedata.node, nodedata.elem]=readmedit(fullfile(path,file));
   elseif(regexp(file,'\.[Ee][Ll][Ee]$'))
       [pathstr,name,ext] = fileparts(fullfile(path,file));
       [nodedata.node, nodedata.elem]=readtetgen(fullfile(pathstr,name));
   elseif(regexp(file,'\.[Jj][Mm][Ee][Ss][Hh]$'))
       nodedata=importjmesh(fullfile(path,file));
   elseif(regexp(file,'\.[Jj][Ss][Oo][Nn]$'))
       nodedata=loadjson(fullfile(path,file));
   end
end
if(exist('nodedata','var'))
    adddatatograph(handles,nodedata);
else
    warndlg('no valid mesh data found','Warning');
end

% --------------------------------------------------------------------
function miLoadSurf_Callback(hObject, eventdata, handles)
nodedata=struct;
nodetype=dummytype;
filters={'*.jmesh;*.off;*.asc;*.smf;*.smf;*.json','3D Mesh files (*.jmesh;*.off;*.asc;*.smf;*.smf;*.json)';...
    '*.jmesh','JSON mesh (*.jmesh)';...
    '*.off','OFF file (*.off)';...
    '*.asc','ASC file (*.asc)';...
    '*.gts','GNU Trangulated Surface file (*.gts)';...
    '*.smf','Simple Model Format (*.smf)';...
    '*.json','JSON file (*.json)';'*.*','All (*.*)'};
[file,path,idx] = uigetfile(filters);
if isequal(file,0)
   return;
else
   if(regexp(file,'\.[Oo][Ff][Ff]$'))
       [nodedata.node, nodedata.face]=readoff(fullfile(path,file));
   elseif(regexp(file,'\.[Aa][Ss][Cc]$'))
       [nodedata.node, nodedata.face]=readasc(fullfile(path,file));
   elseif(regexp(file,'\.[Gg][Tt][Ss]$'))
       [nodedata.node, nodedata.face]=readgts(fullfile(path,file));
   elseif(regexp(file,'\.[Ss][Mm][Ff]$'))
       [nodedata.node, nodedata.face]=readsmf(fullfile(path,file));
   elseif(regexp(file,'\.[Jj][Mm][Ee][Ss][Hh]$'))
       nodedata=importjmesh(fullfile(path,file));
   elseif(regexp(file,'\.[Jj][Ss][Oo][Nn]$'))
       nodedata=loadjson(fullfile(path,file));
   end
end
if(exist('nodedata','var'))
    adddatatograph(handles,nodedata);
else
    warndlg('no valid mesh data found','Warning');
end

% --------------------------------------------------------------------
function nodedata=importjmesh(filename)
data=loadjson(filename);
nodedata=struct;
if(isfield(data,'MeshNode'))
   nodedata.node=data.MeshNode;
end
if(isfield(data,'MeshElem'))
   nodedata.elem=data.MeshElem;
end
if(isfield(data,'MeshSurf'))
   nodedata.face=data.MeshSurf;
end
if(isfield(data,'MeshNodeVal'))
   nodedata.node(:,end+1:end+size(data.MeshNodeVal,2))=data.MeshNodeVal;
end
if(isfield(data,'MeshTetraVal'))
   nodedata.elem(:,end+1:end+size(data.MeshTetraVal,2))=data.MeshTetraVal;
end

% --------------------------------------------------------------------
function adddatatograph(handles,nodedata)
nodetype=getnodetype(nodedata);
if(nodetype.haselem)
    addnodewithdata(handles,nodedata,nodetype,'Tet');
elseif(nodetype.hasface)
    addnodewithdata(handles,nodedata,nodetype,'Surf');
elseif(nodetype.hasnode)
    addnodewithdata(handles,nodedata,nodetype,'Point');
elseif(nodetype.hasvol)
    addnodewithdata(handles,nodedata,nodetype,'Vol');
end

% --- Executes during object creation, after setting all properties.
function axFlow_CreateFcn(hObject, eventdata, handles)
% hObject    handle to axFlow (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axFlow

% --- Executes during object creation, after setting all properties.
function fgI2M_CreateFcn(hObject, eventdata, handles)
% hObject    handle to fgI2M (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

%--------------------------------------------------------------------
function nodetype=getnodetype(nodedata)
nodetype=dummytype;
if(~isstruct(nodedata))
    return;
end
names=fieldnames(nodedata);
for i=1:length(names)
      switch names{i}
           case 'node'
                nodetype.hasnode=1;
           case 'face'
                nodetype.hasface=1;
           case 'elem'
                nodetype.haselem=1;
           case 'vol'
                nodetype.hasvol=1;
      end
      dat=nodedata.(names{i});
      if((isnumeric(dat) || islogical(dat)) && ndims(dat)==3)
           nodetype.hasvol=1;
      end
end
% --------------------------------------------------------------------

function [key,newgraph]=addnodewithdata(handles,nodedata,nodetype,name)

root=get(handles.fgI2M,'userdata');

if(isempty(root))
    root=struct('graph',digraph,'menu',uicontextmenu);
end
if(nargin<4)
    name='x';
end

id=1;
if(~isempty(root.graph.Nodes))
  while(find(strcmp(root.graph.Nodes.Name,sprintf('%s%d',name,id))))
    id=id+1;
  end
end
key=sprintf('%s%d',name,id);

cla(handles.axPreview);
cla(handles.axFlow);
nodedata.preview=getpreview(nodedata,nodetype,[400,400]);

nodeprop=table({key},{nodedata},{nodetype},'VariableNames',{'Name','Data','Type'});
root.graph=addnode(root.graph,nodeprop);
if(nargout>1)
    newgraph=root.graph;
end
updategraph(root, handles);

% --------------------------------------------------------------------
function hobj=updatepreview(root,obj,handles)
cla(handles.axPreview);
view(handles.axPreview,2)
nx=obj.XData(:);
ny=obj.YData(:);
nn=length(nx);
hold(handles.axPreview,'on');
dx=get(handles.axFlow,'xlim');
dy=get(handles.axFlow,'ylim');
wd=min([diff(dx) diff(dy)])/15;
hobj=zeros(1,nn);
set(handles.axPreview,'xlim',dx);
set(handles.axPreview,'ylim',dy);
dim=get(handles.axPreview,'dataaspectratio');
[wfig, hfig]=getwindowsize(handles.fgI2M);
wfig=wfig*dim(2);
hfig=hfig*dim(1);

for i=1:nn
    if(isfield(root.graph.Nodes.Data{i},'preview'))
        hobj(i)=imagesc(nx(i)+[-2*wd 0], ny(i)+[-wd wd]*(wfig/hfig), imresize(root.graph.Nodes.Data{i}.preview,0.25), ...
            'parent',handles.axPreview);
    end
end
set(hobj,'userdata',obj);
set(hobj,'UIContextMenu',root.menu);
set(handles.axPreview,'xlim',dx);
set(handles.axPreview,'ylim',dy);
axis(handles.axPreview,'off');
hold(handles.axPreview,'off');
% --------------------------------------------------------------------

function [width, height]=getwindowsize(fig)
oldunits = get(fig, 'Units');
set(fig, 'Units', 'pixels');
figpos = get(fig, 'Position');
set(fig, 'Units', oldunits);
width=figpos(3);
height=figpos(4);

function [hg,hobj]=updategraph(root, handles)
set(handles.fgI2M,'userdata',root);
hg=plot(root.graph,'parent',handles.axFlow,'ArrowSize',15);
hobj=updatepreview(root,hg,handles);
% set(hg,'Selected','on');
axis(handles.axFlow,'off');
set(handles.axFlow,'XColor','none')
set(hg,'UIContextMenu',root.menu);

% --------------------------------------------------------------------
function miEllipsoid_Callback(hObject, eventdata, handles)
prompt = {'Center (1x3 vector):',...
    'Radii (a scalar, or 1x3 or 1x5 vector):',...
    'Max element volume (scalar, 0 - only create surface):',...
    'Circle division:'};
title = 'Create an Ellipsoid Mesh';
dims = [1 1 1 1];
definput = {'[0 0 0]','[50 30 20]','3','30'};
res = inputdlg(prompt,title,dims,definput);
if(isempty(res))
    return;
end
newtype=dummytype;
opt=str2num(res{3});
maxvol=str2num(res{4});

if(maxvol==0)
   [newdata.node,newdata.face]=meshanellip(eval(res{1}),eval(res{2}),opt);
else
   [newdata.node,newdata.face,newdata.elem]=meshanellip(eval(res{1}),...
     eval(res{2}),opt, maxvol);
   newtype.haselem=1;
end
newtype.hasnode=1;
newtype.hasface=1;

if(exist('newdata','var') && exist('newtype','var'))
    newkey=addnodewithdata(handles,newdata,newtype,'Cyl');
end

% --------------------------------------------------------------------
function miLattice_Callback(hObject, eventdata, handles)
prompt = {'X-lattice range (a vector):',...
    'Y-lattice range (a vector):',...
    'Z-lattice range (a vector):',...
    'Max element volume (scalar, 0 - only create surface):'};
title = 'Create Lattice Grid Mesh';
dims = [1 1 1 1];
definput = {'[1 100]','[1 50]','[1 10 30]','30'};
res = inputdlg(prompt,title,dims,definput);
if(isempty(res))
    return;
end
newtype=dummytype;
maxvol=str2num(res{4});
if(maxvol==0)
   [newdata.node,newdata.face]=latticegrid(eval(res{1}),eval(res{2}),eval(res{3}));
else
   [no,fc,c0]=latticegrid(eval(res{1}),eval(res{2}),eval(res{3}));
   [newdata.node,newdata.elem,newdata.face]=surf2mesh(no,fc,[],[],1,maxvol,c0);
   newtype.haselem=1;
end
newtype.hasnode=1;
newtype.hasface=1;

if(exist('newdata','var') && exist('newtype','var'))
    newkey=addnodewithdata(handles,newdata,newtype,'Lattice');
end

% --------------------------------------------------------------------
function miMeshgrid5_Callback(hObject, eventdata, handles)
prompt = {'X-lattice range (a vector):',...
    'Y-lattice range (a vector):',...
    'Z-lattice range (a vector):'};
title = 'Create Meshgrid (5 tet/cell) Mesh';
dims = [1 1 1];
definput = {'1:10','1:8','1:5'};
res = inputdlg(prompt,title,dims,definput);
if(isempty(res))
    return;
end
newtype=dummytype;

[newdata.node,newdata.elem]=meshgrid5(eval(res{1}),eval(res{2}),eval(res{3}));
newdata.face=volface(newdata.elem);
newtype.hasnode=1;
newtype.hasface=1;
newtype.haselem=1;

if(exist('newdata','var') && exist('newtype','var'))
    newkey=addnodewithdata(handles,newdata,newtype,'Meshgrid5_');
end

% --------------------------------------------------------------------
function miMeshgrid6_Callback(hObject, eventdata, handles)
prompt = {'X-lattice range (a vector):',...
    'Y-lattice range (a vector):',...
    'Z-lattice range (a vector):'};
title = 'Create Meshgrid (6 tet/cell) Mesh';
dims = [1 1 1];
definput = {'1:10','1:8','1:5'};
res = inputdlg(prompt,title,dims,definput);
if(isempty(res))
    return;
end
newtype=dummytype;

[newdata.node,newdata.elem]=meshgrid5(eval(res{1}),eval(res{2}),eval(res{3}));
newdata.face=volface(newdata.elem);
newtype.hasnode=1;
newtype.hasface=1;
newtype.haselem=1;

if(exist('newdata','var') && exist('newtype','var'))
    newkey=addnodewithdata(handles,newdata,newtype,'Meshgrid6_');
end

% --------------------------------------------------------------------
function miOpen_Callback(hObject, eventdata, handles)
filter = {'*.mat';'*.*'};
root=get(handles.fgI2M,'userdata');

[file, path] = uigetfile(filter,'Load workspace');
if isequal(file,0)
   return;
else
   data=load(fullfile(path,file));
end
if(isfield(data,'i2mworkspace'))
    root.graph=data.i2mworkspace;
    updategraph(root,handles);
else
    warndlg('no saved workspace found','Warning');
end

% --------------------------------------------------------------------
function miSaveAll_Callback(hObject, eventdata, handles)
root=get(handles.fgI2M,'userdata');
i2mworkspace=root.graph;
filter = {'*.mat';'*.*'};
[file, path] = uiputfile(filter,'Save workspace');
if ~isequal(file,0) && ~isequal(path,0)
   save(fullfile(path,file),'i2mworkspace');
end

% --------------------------------------------------------------------
function miLoadVar_Callback(hObject, eventdata, handles)

nodedata=struct;
[file,path] = uigetfile({'*.mat','MATLAB data (*.mat)'});
if isequal(file,0)
   return;
else
   data=load(fullfile(path,file));
   vars=fieldnames(data);
   [idx, isok]=listdlg('ListString',vars,...
       'PromptString','Select a 3D array:');
   if(~isok)
       return;
   end
   for i=1:length(idx)
       dat=data.(vars{idx(i)});
       if(ndims(dat)==3 && ~isfield(nodedata,'vol'))
           nodedata.vol=dat;
       else
           nodedata.(vars{idx(i)})=dat;
       end
   end
end
adddatatograph(handles,nodedata);

% --------------------------------------------------------------------
function miExit_Callback(hObject, eventdata, handles)
close(handles.fgI2M);

% --------------------------------------------------------------------
function miRefresh_Callback(hObject, eventdata, handles)
root=get(handles.fgI2M,'userdata');
updategraph(root, handles);


text:
function [Q fcnEvals iter] = adaptiveSimpson(fcn, a, b, varargin)
% adaptiveSimpson - Numerically evaluate integral, adaptive Simpson quadrature.
%
% function [Q fcnEvals iter] = adaptiveSimpson(fcn, a, b, varargin)
%
% (c) Matthias Conrad and Nils Papenberg (2007-08-03)
%
% Authors:
%   Matthias Conrad (e-mail: conrad@tiaco.de)
%   Nils Papenberg  (e-mail: papenber@math.uni-luebeck.de)
%
% Version:
%		Release date: 2008-08-12   Version: 1.2
%   MATLAB Version 7.5.0.338 (R2007b)
%
% Description:
%   The adaptive Simpson algorithm programmed in an iterative not recursive
%   manner
%
% Input arguments:
%   fcn             - function to be integrated
%   a               - first point of interval
%   b               - final point of interval
%   #varargin       - further options of algorithm
%     tol           - tolerance accuracy of quadrature [ 1e-6 ]
%     parts         - initial number of partitions [ 2 ]
%     maxFcnEvals   - maximal number of function evaluations allowed [ 20000 ]
%     maxParts      - maximal number of partitions allowed [ 8000 ]
%
% Output arguments:
%   Q               - numerical integral of function fcn on [a,b]
%   fcnEvals        - number of function evaluations
%   iter            - number of iterations
%
% Details:
%   This function behavior is similar to of Matlab integrated function "quadv".
%
%   Example:
%     Q = adaptiveSimpson(@(x) [-cos(50*x); sin(x)], 0, pi, 'tol', 1e-6)
%
% References:
%   [1] Gander, W. & Gautschi, W. Adaptive Quadrature - Revisited
%       Eidgenoessische technische Hochschule Zuerich, 2000.

% check scalar limits of interval
if ~isscalar(a) || ~isscalar(b)
  error('Matlab:adaptiveSimpson:Limits',...
    'The limits of integration must be scalars.');
end

% default values
tol = 1e-6; parts = 2; maxFcnEvals = 20000; maxParts = 8000;

% rewrite default options if needed
for j = 1 : length(varargin) / 2
  eval([varargin{2 * j - 1},'=varargin{',int2str(2 * j),'};']);
end

% initial values, termination constant (incl. quadr factor 15), parts of interval and integral value
tol = 15 * tol; m = parts;  parts = 4 * parts + 1; Q = 0;
iter = 0; fcnEvals = 0; maxResolution = 0; minH = eps(b - a) / 1024;
poleWarning = 0;

% check if interval has infinite boundaries, in case substitute function
if ~isfinite(a) || ~isfinite(b)
  warning('Matlab:adaptiveSimpson:infiniteInterval',...
    'The integral has an infinite interval; proceed with a substitution of function on finite interval.')
  if ~isfinite(a) && isfinite(b)
    [Q fcnEvals iter] = adaptiveSimpson(fcn, 0, b, varargin);
    fcn = @(t) infiniteLeft(t, fcn);
    a = 0; b = 1;
  elseif isfinite(a) && ~isfinite(b)
    [Q fcnEvals iter] = adaptiveSimpson(fcn, a, 0, varargin);
    fcn = @(t) infiniteRight(t, fcn);
    a = 0; b = 1;
  else
    fcn = @(t) infiniteBoth(t, fcn);
    a = - pi / 2; b = pi / 2;
  end
end

% choice of initial evaluation points
t = [a + (b - a) / m * (kron(ones(1,m), [0, 0.27158, 0.72842]) + kron(0: m - 1, [1 1 1])), b];
H = diff(t);
t = [t(1:end-1); t(1:end-1) + H/4; t(1:end-1) + H/2; t(1:end-1) + 3 * H/4;];
t = [t(:);b]';

% initialize equidistant mesh and dimension of fcn
y = fcn(t); fcnEvals = fcnEvals + length(y); n = size(y,1);

% avoid infinities at start point of interval
if any(~isfinite(y(:,1)))
  y(:,1) = fcn(a + eps(superiorfloat(a,b)) * (b - a));
  fcnEvals = fcnEvals + 1;
  poleWarning = 1;
end

% avoid infinities at end point of interval
if any(~isfinite(y(:, end)))
  y(:, end) = fcn(b - eps(superiorfloat(a,b)) * (b - a));
  fcnEvals = fcnEvals + 1;
  poleWarning = 1;
end

% poles at initial points
if ~isempty(find(~isfinite(max(abs(y))))), poleWarning = 1; end

% initialize interval boundaries and values
A = t(1:4:end-1); yA = y(:, 1:4:end-1);
B = t(5:4:end);   yB = y(:, 5:4:end);
C = t(3:4:end-1); yC = y(:, 3:4:end-1);
D = t(2:4:end-1); yD = y(:, 2:4:end-1);
E = t(4:4:end-1); yE = y(:, 4:4:end-1);

% adaptive Simpson iteration
while 1

  % number of iteration
  iter = iter + 1;

  % Simpson formulas (on rough and fine grid)
  Q1 = kron(H, ones(n,1)) / 6 .* (yA + 4 * yC + yB);
  Q2 = kron(H, ones(n,1)) / 12 .* (yA + 4 * yD + 2 * yC + 4 * yE + yB);

  % difference of Simpson formulas
  diffQ = Q2 - Q1; diffQ(find(isnan(diffQ))) = 0;

  % intervals which do not fulfill termination criterion
  idx = find(max(abs(diffQ), [], 1) > tol);

  % intervals fulfill termination criterion
  idxQ = setdiff(1:length(A), idx);

  % check stop criterions
  STOP1 = isempty(idx); % check regular termination
  STOP2 = fcnEvals > maxFcnEvals;  % check maximal function evaluations
  STOP3 = 2 * length(idx) > maxParts; % check maximal partition

  % regular termination
  if STOP1
    Q = Q + sum(Q2 + diffQ / 15, 2);
    break
  end

  % check if maximal resolution reached
  idxH = find(abs(H) < minH);
  if ~isempty(idxH)
    Q = Q + sum(Q2(idxH) + diffQ(idxH) / 15, 2);
    idx = setdiff(idx, idxH);
    idxQ = setdiff(idxQ, idxH);
    maxResolution = 1;
    % termination criterion
    if isempty(idx), break, end
  end

  % maximal function evaluations reached
  if STOP2
    warning('Matlab:adaptiveSimpson:MaxEvaluations',...
      'The maximal number of function evaluations reached; singularity likely.')
    Q = Q + sum(Q2 + diffQ / 15, 2);
    break
  end

  % maximal partition reached
  if STOP3
    warning('Matlab:adaptiveSimpson:parts',...
      'The maximal number of parts reached.')
    Q = Q + sum(Q2 + diffQ / 15, 2);
    break
  end

  % update quadrature value
  Q = Q + sum(Q2(:,idxQ) + diffQ(:,idxQ) / 15, 2);

  % update A, B, and C
  t = zeros(1, 2*length(idx));
  t(1:2:end) = A(idx); t(2:2:end) = C(idx); A = t;
  t(1:2:end) = C(idx); t(2:2:end) = B(idx); B = t;
  t(1:2:end) = D(idx); t(2:2:end) = E(idx); C = t;

  % update yA, yB, and yC
  y = zeros(n, 2*length(idx));
  y(:,1:2:end) = yA(:,idx); y(:,2:2:end) = yC(:,idx); yA = y;
  y(:,1:2:end) = yC(:,idx); y(:,2:2:end) = yB(:,idx); yB = y;
  y(:,1:2:end) = yD(:,idx); y(:,2:2:end) = yE(:,idx); yC = y;

  % update interval length
  H = B - A;

  % update D and E by interval bisection
  D = (A + H / 4); E = (A + 3 * H / 4);

  % calculate yD and yE
  y = fcn([D,E]); fcnEvals = fcnEvals + 4 * length(idx);

  % poles at new points
  if ~isempty(find(~isfinite(max(abs(y))))), poleWarning = 1; end

  % assign new values ob yD and yE
  yD = y(:,1:end/2); yE = y(:,end/2+1: end);

end

% display warnings
if any(~isfinite(Q))
  warning('Matlab:adaptiveSimpson:Infinite',...
    'The Quadrature of the function reached infinity or is Not-a-Number.')
end
if maxResolution
  warning('Matlab:adaptiveSimpson:MaxResolution',...
    'The maximal resolution of partial interval reached; singularity likely.')
end
if poleWarning
  warning('Matlab:adaptiveSimpson:PoleDetection',...
    'A detection of a pole; singularity likely.')
end

return

% substitute function interval [-inf, 0] on [0, 1]
function f = infiniteLeft(t, fcn)
f = fcn(log(t));
f = f ./ kron(ones(size(f,1),1), t);
return

% substitute function interval [0, inf] on [0, 1]
function f = infiniteRight(t, fcn)
f = fcn(-log(t));
f = f ./ kron(ones(size(f,1),1), t);
return

% substitute function interval [-inf, inf] on [-pi / 2, pi / 2]
function f = infiniteBoth(t, fcn)
f = fcn(tan(t));
f = f ./ kron(ones(size(f,1),1), cos(t).^2);
return

author: Sable

repo: mcbench-benchmarks


text:
function out = iembed(ob, x)
% recall that imask/omask and idim/odim are swapped on ctranspose
out = fatrix2_embed(ob.imask, ob.idim, x);

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function c = unifytracelengths(c)
   %unifytracelengths   make all traces the same length (mode length)
   %   C = unifytracelengths(C) if traces are of unequal lengths, then they
   %   are either trimmed or zero-padded to the same size.  The final length
   %   will be the mode length of the original set of traces.
   %
   %   See also mode

   if isempty(c.traces)
      return
   end

   lengths = [c.traces.nsamples];
   ideal = mode(lengths);

   tooShort = lengths < ideal;
   tooLong = lengths > ideal;

   for i = 1:numel(lengths)
      if tooShort(i)
         c.traces(i).data(ideal) = 0;
      elseif tooLong(i)
         c.traces(i).data = c.traces(i).data(1:ideal);
      end
   end
end

author: geoscience-community-codes

repo: GISMO


text:
% -------------------------------------------------------------------------
% function hfssSweepAnalysis(fid, Name, Analysis, Type, Variables, Data,
%                            Units, [Sync])
%
% Description :
% -------------
% Generates a simple linear sweep analysis for a given variable.
%
% Parameters :
% ------------
% fid       - file identifier of the HFSS script file.
% Name      - name of the sweep analysis to be created.
% Analysis  - name of the associated analysis setup.
% Type      - cell array with sweep types. Supported:
%                SINGLE: array of single values.
%                LIN:    follows the structure "LIN Start Stop Step", i.e.
%                        Start:Step:Stop
%                LINC:   follows the structure "LINC Start Stop Count",
%                        i.e. linspace(Start, Stop, Count)
% Variables - cell array with the name of the variable of which the sweep
%             will be performed.
% Data      - cell array containing numerical values of [Start, Stop,
%             Step] for each variable.
% Units     - for Data, can be 'mm', 'cm', 'meters', etc.
% Sync      - (optional, cell array) set to 1, 2, ... to sync variable
%             sweeps in different groups. Default: 0.
%
% Note :
% ------
%
% Examples :
% ---------
% hfssSweepAnalysis(fid, 'ParSetup1', 'MySetup', 'LIN', 'var', ...
%    [1, 9, 1], 'mm');
% hfssSweepAnalysis(fid, 'ParSetup2', 'MySetup', {'LIN', 'LIN'}, ...
%     {'pB', 'pA'}, {[0.1, 3.5, 0.2], [0.1, 3.5, 0.2]}, 'mm', {1, 1});
%
% hfssSweepAnalysis(fid, 'ParSetup3', 'MySetup', {'Single', 'LIN'}, ...
%     {'pC', 'pD'}, {[pi, exp(1), -sqrt(2)], [0.1, 3.5, 0.2]}, 'rad', ...
%     {1, 1});
%
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% CHANGELOG
%
% 03-Sep-2020: *Initial release (only LIN sweep type supported).
% 09-Sep-2020: *LINC sweep type supported.
%              *Sync is now a cell array, allowing to syncronize variables
%               in different groups.
% 22-Dec-2020: *Array of single value support through SINGLE keyword.
% -------------------------------------------------------------------------

% -------------------------------------------------------------------------
% Written by Daniel R. Prado
% danysan@gmail.com / drprado@tsc.uniovi.es
% 03 September 2020
% -------------------------------------------------------------------------
function hfssSweepAnalysis(fid, Name, Analysis, Type, Variables, Data, ...
                           Units, Sync)

% Arguments processor.
if (nargin < 7)
	error('Insufficient # of arguments !');
end
if (nargin == 7)
    Sync = zeros(numel(Variables), 1);
end

if (~iscell(Type))
    Type = {Type};
end
if (~iscell(Variables))
    Variables = {Variables};
end
if (~iscell(Data))
    Data = {Data};
end
if (~iscell(Sync))
    Sync = num2cell(Sync);
end

% Preamble
fprintf(fid, '\n');
fprintf(fid, 'Set oModule = oDesign.GetModule("Optimetrics")\n');

% Command
fprintf(fid, 'oModule.InsertSetup "OptiParametric", _\n');
fprintf(fid, '\tArray("NAME:%s", _\n', Name);
fprintf(fid, '\t"IsEnabled:=", true, _\n');
fprintf(fid, '\tArray("NAME:ProdOptiSetupDataV2", _\n');
fprintf(fid, '\t\t"SaveFields:=", false, _\n');
fprintf(fid, '\t\t"CopyMesh:=", false, _\n');
fprintf(fid, '\t\t"SolveWithCopiedMeshOnly:=", true), _\n');
fprintf(fid, '\tArray("NAME:StartingPoint"), _\n');
fprintf(fid, '\t"Sim. Setups:=", Array("%s"), _\n', Analysis);
fprintf(fid, '\tArray("NAME:Sweeps", _\n');

for n = 1:numel(Variables)
    fprintf(fid, '\t\tArray("NAME:SweepDefinition", _\n');
    fprintf(fid, '\t\t"Variable:=", "%s", _\n', Variables{n});
    switch Type{n}
        case 'SINGLE'
            strUnits  = repmat('deg ', numel(Data{n}), 1);
            strData   = num2str(Data{n});
            strAux    = [strData, strUnits];
            str       = strtrim(reshape(strAux.', 1, numel(strAux)));
            fprintf(fid, ['\t\t"Data:=", "', str, '", _\n']);
        case 'LIN'
            fprintf(fid, '\t\t"Data:=", "LIN %f%s %f%s %f%s", _\n', ...
                Data{n}(1), Units, Data{n}(2), Units, Data{n}(3), Units);
        case 'LINC'
            fprintf(fid, '\t\t"Data:=", "LINC %f%s %f%s %i", _\n', ...
                Data{n}(1), Units, Data{n}(2), Units, Data{n}(3));
        otherwise
            error('Wrong type!');
    end
    fprintf(fid, '\t\t"OffsetF1:=", false, _\n');
    fprintf(fid, '\t\t"Synchronize:=", %i)', Sync{n});

    if (n < numel(Variables))
        fprintf(fid, ', _\n'); % Don't include comma in last item.
    end
end
fprintf(fid, '\t\t), _\n');
fprintf(fid, '\tArray("NAME:Sweep Operations"), _\n');
fprintf(fid, '\tArray("NAME:Goals"))\n');

author: yuip

repo: hfss-api

save_path: github-repos/MATLAB/yuip-hfss-api


text:
function ebsd = subGrid(ebsd,ind)
% indexing of EBSD data
%
% Syntax
%   subGrid(ebsd,ind)
%

id = ebsd.id(ind);

% determine bounding rectangle
if islogical(ind)
  yMin = find(any(ind,2),1,"first");
  yMax = find(any(ind,2),1,"last");
  xMin = find(any(ind,1),1,"first");
  xMax = find(any(ind,1),1,"last");
else
  [x,y] = ind2sub(ebsd,ind);
  yMin = min(y(:));
  yMax = max(y(:));
  xMin = min(x(:));
  xMax = max(x(:));
end

mask = false(size(ebsd));
mask(yMin:yMax,xMin:xMax) = true;

ebsd = subSet(ebsd,mask,'keepGrid');
ebsd = reshape(ebsd,yMax-yMin+1,xMax-xMin+1);

% fill with NaN
nanId = ~ismember(ebsd.id,id);
ebsd.phaseId(nanId) = NaN;
ebsd.rotations(nanId) = NaN;

author: mtex-toolbox

repo: mtex

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


text:
function [magTS, phaseTS] = cReconBW(fName, slice)

% [magTS, phaseTS] = cReconBW(fName, slice);
%
% Perform a complex recon using Alex Wade's tstrecon code.
% Returns pair of single-slice time-series arrays of size:
% ny x nx x nFrames -- magTS is magnitude and phaseTS is
% phase. Time-duration and spatial cropping are performed,
% but no trend removal. Inputs are pFile name [fName], and
% slice index [slice]. Assumes pwd points to valid session
% directory, and uses global mrSESSION.
%
% DBR 8/00

global mrSESSION

% Do recon, extract complex time series:
cTS = tstrecon3Tb([], fName, slice, slice);
cTS = cTS{1}';

% Remove junk frames:
f0 = mrSESSION.junkFirstFrames+1;
nFrames = mrSESSION.nFrames;
fEnd = f0 + nFrames - 1;
cTS = cTS(f0:fEnd, :);

% Crop:
cTS = reshape(cTS, [nFrames, mrSESSION.fullSize]); % reshape the array to full size
x0 = mrSESSION.tseriesCrop(1, 1);
xN = mrSESSION.tseriesCrop(2, 1);
y0 = mrSESSION.tseriesCrop(1, 2);
yN = mrSESSION.tseriesCrop(2, 2);
cTS = cTS(:, y0:yN, x0:xN);

%Shuffle to standard t-series shape:
cTS = reshape(cTS, nFrames, (yN-y0+1)*(xN-x0+1));

% Get magnitude and phase:
magTS = abs(cTS);
phaseTS = angle(cTS);


text:
function [ value, t ] = r4_random ( t, n )

%*****************************************************************************80
%
%% R4_RANDOM is a portable pseudorandom number generator.
%
%  Discussion:
%
%    This random number generator is portable amoung a wide variety of
%    computers.  It generates a random number between 0.0 and 1.0
%    according to the algorithm presented by Bays and Durham.
%
%    The motivation for using this scheme, which resembles the
%    Maclaren-Marsaglia method, is to greatly increase the period of the
%    random sequence.  If the period of the basic generator is P,
%    then the expected mean period of the sequence generated by this
%    generator is given by
%
%      new mean P = sqrt ( pi * factorial ( N ) / ( 8 * P ) ),
%
%    where factorial ( N ) must be much greater than P in this
%    asymptotic formula.  Generally, N should be 16 to maybe 32.
%
%  Modified:
%
%    25 September 2011
%
%  Author:
%
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Carter Bays, Stephen Durham,
%    Improving a Poor Random Number Generator,
%    ACM Transactions on Mathematical Software,
%    Volume 2, Number 1, March 1976, pages 59-64.
%
%  Parameters:
%
%    Input, integer N.  The absolute value of N is the number
%    of random numbers in an auxiliary table.  Note though that abs(N)+1 is
%    the number of items in array T.  If N is positive and differs from its
%    value in the previous invocation, then the table is initialized for
%    the new value of N.  If N is negative, abs(N) is the number of items
%    in an auxiliary table, but the tables are now assumed already to
%    be initialized.  This option enables the user to save the table T at
%    the end of a long computer run and to restart with the same sequence.
%    Normally, this function would be called at most once with negative N.
%    Subsequent invocations would have N positive and of the correct magnitude.
%
%    Input/output, real T(abs(N)+1), an array of random numbers
%    from a previous invocation of this function.  Whenever N is positive
%    and differs from the old N, the table is initialized.  The first
%    abs(N) numbers are the table discussed in the reference, and the
%    last value is Y.  This array may be saved in order to restart a sequence.
%
%    Output, real VALUE, a random number between 0.0 and 1.0.
%
  persistent floatn
  persistent nold

  if ( isempty ( floatn ) )
    floatn = -1.0;
    nold = -1;
  end

  if ( n ~= nold )

    nold = abs ( n );
    floatn = nold;

    if ( n < 0 )

      dummy = r4_rand ( t(nold+1) );

    else

      for i = 1 : nold
        t(i) = r4_rand ( 0.0 );
      end
      t(nold+1) = r4_rand ( 0.0 );

    end

  end

  j = r4_aint ( t(nold+1) * floatn + 1.0 );
  t(nold+1) = t(j);
  r4_random = t(j);
  t(j) = r4_rand ( 0.0 );

  return
end

author: johannesgerer

repo: jburkardt-m