text:
function bnd = prepare_mesh_manual(cfg, mri)

% PREPARE_MESH_MANUAL is called by PREPARE_MESH and opens a GUI to manually
% select points/polygons in an mri dataset.
%
% It allows:
%   Visualization of 3d data in 3 different projections
%   Adjustment of brightness for every slice
%   Storage of the data points in an external .mat file
%   Retrieval of previously saved data points
%   Slice fast scrolling with keyboard arrows
%   Polygons or points selection/deselection
%
% See also PREPARE_MESH_SEGMENTATION, PREPARE_MESH_HEADSHAPE

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

% FIXME: control slice's cmap referred to abs values
% FIXME: clean structure slicedata
% FIXME: check function assign3dpoints

global obj

mri = ft_checkdata(mri, 'datatype', {'volume', 'segmentation'});

hasheadshape = isfield(cfg, 'headshape');
hasbnd       = isfield(cfg, 'bnd');  % FIXME why is this in cfg?
hasmri       = nargin>1;

% check the consistency of the input arguments
if hasheadshape && hasbnd
  ft_error('you should not specify cfg.headshape and cfg.bnd simultaneously');
end

% check the consistency of the input arguments
if ~hasmri
  % FIXME give a warning or so?
  mri.anatomy = [];
  mri.transform = eye(4);
else
  % ensure that it is double precision
  mri.anatomy = double(mri.anatomy);
end

% start with an empty boundary
bnd.pnt = [];
bnd.tri = [];

if hasheadshape
  if ~isempty(cfg.headshape)
    % start with the headshape
    [bnd.pos, bnd.tri] = headsurface([], [], 'headshape', cfg.headshape);
  end
elseif hasbnd
  if ~isempty(cfg.bnd)
    % start with the prespecified boundaries
    bnd = cfg.bnd;
  end
end

% creating the GUI
fig = figure;

% initialize some values
slicedata = []; points = bnd.pnt;
axs = floor(size(mri.anatomy,1)/2);

% initialize properties
set(fig, 'KeyPressFcn',@keypress);
set(fig, 'CloseRequestFcn', @cb_close);
setappdata(fig,'data',mri.anatomy);
setappdata(fig,'prop',{1,axs,0,[]});
setappdata(fig,'slicedata',slicedata);
setappdata(fig,'points',points);
setappdata(fig,'box',[]);

% add GUI elements
cb_creategui(gca);
cb_redraw(gca);
waitfor(fig);

% close sequence
global pnt
try
  tmp = pnt;
  clear global pnt
  pnt = tmp;
  clear tmp
  [tri] = projecttri(pnt);
  bnd.pnt = pnt;
  bnd.tri = tri;
catch
  bnd.pnt = [];
  bnd.tri = [];
end

function cb_redraw(hObject, eventdata, handles)
fig  = get(hObject, 'parent');
prop = getappdata(fig,'prop');
data = getappdata(fig,'data');
slicedata = getappdata(fig,'slicedata');
points = getappdata(fig,'points');

% draw image
try
  cla
  % i,j,k axes
  if (prop{1}==1) %jk
    imh=imagesc(squeeze(data(prop{2},:,:)));
    xlabel('k');
    ylabel('j');
    colormap bone
  elseif (prop{1}==2) %ik
    imh=imagesc(squeeze(data(:,prop{2},:)));
    xlabel('k');
    ylabel('i');
    colormap bone
  elseif (prop{1}==3) %ij
    imh=imagesc(squeeze(data(:,:,prop{2})));
    xlabel('j');
    ylabel('i');
    colormap bone
  end
  axis equal; axis tight
  brighten(prop{3});
  title([ 'slice ' num2str(prop{2}) ])
catch
end
if ~ishold,hold on,end
% draw points and lines
for zz = 1:length(slicedata)
  if(slicedata(zz).proj==prop{1})
    if(slicedata(zz).slice==prop{2})
      polygon = slicedata(zz).polygon;
      for kk=1:length(polygon)
        if ~isempty(polygon{kk})
          [xs,ys]=deal(polygon{kk}(:,1),polygon{kk}(:,2));
          plot(xs, ys, 'g.-');
        end
      end
    end
  end
end % end for

% draw other slices points
points = getappdata(fig,'points');
pnts   = round(points);
if ~isempty(pnts)
  % exclude polygons inside the same slice
  indx  = find(points(:,prop{1})==prop{2});
  tmp   = ones(1,size(points,1));
  tmp(indx)  = 0;
  pnts  = points(find(tmp),:);
  % calculate indexes of other proj falling in current slice
  indx2 = find(round(pnts(:,prop{1}))==prop{2});
  rest  = setdiff([1 2 3],prop{1});
  % plot the points
  for jj=1:length(indx2)
    [x,y] = deal(pnts(indx2(jj),rest(1)),pnts(indx2(jj),rest(2)));
    plot(y,x,'g*')
  end
end

% add blue cross markers in case of box selection
box        = getappdata(fig,'box');
point2mark = getappdata(fig,'point2mark');

if ~isempty(point2mark)
  plot(point2mark.x,point2mark.y,'marker','+')
end

if ~isempty(box)
  for kk=1:size(box,1)
    proj   = box(kk,1);
    slice  = box(kk,2);
    rowmin = box(kk,3);
    rowmax = box(kk,4);
    colmin = box(kk,5);
    colmax = box(kk,6);
    aaa = get(findobj(fig, 'color', 'g'));
    for ii=1:length(aaa)
      if ispolygon(aaa(ii))
        L = length(aaa(ii).YData);
        for jj=1:L
          cond = lt(aaa(ii).YData(jj),rowmax).*gt(aaa(ii).YData(jj),rowmin).* ...
            lt(aaa(ii).XData(jj),colmax).*gt(aaa(ii).XData(jj),colmin);
          if cond
            plot(aaa(ii).XData(jj),aaa(ii).YData(jj),'marker','+')
          end
        end
      elseif ispoint(aaa(ii))
        cond = lt(aaa(ii).YData,rowmax).*gt(aaa(ii).YData,rowmin).* ...
          lt(aaa(ii).XData,colmax).*gt(aaa(ii).XData,colmin);
        % the test after cond takes care the box doesnt propagate in 3D
        if (cond && (proj == prop{1}) && (slice == prop{2}))
          plot(aaa(ii).XData,aaa(ii).YData,'marker','+')
        end
      end
    end
  end
end

if get(findobj(fig, 'tag', 'toggle axes'), 'value')
  axis on
else
  axis off
end

function cb_creategui(hObject, eventdata, handles)

fig = get(hObject, 'parent');
% define the position of each GUI element
% constants
CONTROL_WIDTH  = 0.10;
CONTROL_HEIGHT = 0.075;
CONTROL_HOFFSET = 0.8;
CONTROL_VOFFSET = 0.6;

% control buttons
uicontrol('tag', 'view', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'View', 'value', 1, 'callback', @which_task1);
ft_uilayout(fig, 'tag', 'view',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET+4*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);

uicontrol('tag', 'ins', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'Ins', 'value', 0, 'callback', @which_task2);
ft_uilayout(fig, 'tag', 'ins',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET+2.75*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'del', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'Del', 'value', 0, 'callback', @which_task3);
ft_uilayout(fig, 'tag', 'del',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET+2.75*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'slice', 'parent', fig, 'units', 'normalized', 'style', 'text', 'string', 'slice', 'value', [], 'callback', []);
ft_uilayout(fig, 'tag', 'slice',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET+1.25*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'min', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', '<', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'min',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET+1*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'max', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', '>', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'max',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET+1*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'brightness', 'parent', fig, 'units', 'normalized', 'style', 'text', 'string', 'brightness', 'value', [], 'callback', []);
ft_uilayout(fig, 'tag', 'brightness',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET-0.75*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'plus', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', '+', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'plus',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-1*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'minus', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', '-', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'minus',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-1*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'toggle axes', 'parent', fig, 'units', 'normalized', 'style', 'checkbox', 'string', 'axes', 'value', 0, 'callback', @cb_redraw);
ft_uilayout(fig, 'tag', 'toggle axes',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET-2.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);

uicontrol('tag', 'plane', 'parent', fig, 'units', 'normalized', 'style', 'text', 'string', 'plane', 'value', [], 'callback', []);
ft_uilayout(fig, 'tag', 'plane',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-3.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'one', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'jk', 'value', 0, 'callback', @cb_btnp1);
ft_uilayout(fig, 'tag', 'one',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-3.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);
uicontrol('tag', 'two', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'ik', 'value', 0, 'callback', @cb_btnp2);
ft_uilayout(fig, 'tag', 'two',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-4*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'three', 'parent', fig, 'units', 'normalized', 'style', 'radiobutton', 'string', 'ij', 'value', 0, 'callback', @cb_btnp3);
ft_uilayout(fig, 'tag', 'three',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-4*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'mesh', 'parent', fig, 'units', 'normalized', 'style', 'text', 'string', 'mesh', 'value', [], 'callback', []);
ft_uilayout(fig, 'tag', 'mesh',  'BackgroundColor', [0.8 0.8 0.8], 'width', 2*CONTROL_WIDTH, 'height', CONTROL_HEIGHT, 'vpos', CONTROL_VOFFSET-5.75*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);

uicontrol('tag', 'smoothm', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'smooth', 'value', [], 'callback', @smooth_mesh);
ft_uilayout(fig, 'tag', 'smoothm',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-6*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'viewm', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'view', 'value', [], 'callback', @view_mesh);
ft_uilayout(fig, 'tag', 'viewm',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-6*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

uicontrol('tag', 'cancel', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'cancel', 'value', [], 'callback', @cancel_mesh);
ft_uilayout(fig, 'tag', 'cancel',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-6.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET);
uicontrol('tag', 'ok', 'parent', fig, 'units', 'normalized', 'style', 'pushbutton', 'string', 'ok', 'value', [], 'callback', @cb_btn);
ft_uilayout(fig, 'tag', 'ok',  'BackgroundColor', [0.8 0.8 0.8], 'width', CONTROL_WIDTH, 'height', CONTROL_HEIGHT/2, 'vpos', CONTROL_VOFFSET-6.5*CONTROL_HEIGHT, 'hpos', CONTROL_HOFFSET+CONTROL_WIDTH);

function cb_btn(hObject, eventdata, handles)
fig  = get(gca, 'parent');
prop = getappdata(fig,'prop');
slice = prop{2};
br    = prop{3};
datdim = size(getappdata(fig,'data'));
% p1 = [];p2 = [];p3 = [];
p1  = get(findobj('string','jk'),'value');
p2  = get(findobj('string','ik'),'value');
p3  = get(findobj('string','ij'),'value');
set(findobj('string','-'),'Value',prop{3});
set(findobj('string','+'),'Value',prop{3});
beta = get(findobj('string','-'),'Value');

if (p1)     %jk
  setappdata(fig,'prop',{1,round(datdim(1)/2),br,prop{4}});
  cb_redraw(gca);
elseif (p2) %ik
  setappdata(fig,'prop',{2,round(datdim(2)/2),br,prop{4}});
  cb_redraw(gca);
elseif (p3) %ij
  setappdata(fig,'prop',{3,round(datdim(3)/2),br,prop{4}});
  cb_redraw(gca);
end
if strcmp(get(hObject,'string'),'-')
  beta=beta-0.05;
  set(findobj('string','-'),'Value',beta);
  set(findobj('string','+'),'Value',beta);
  setappdata(fig,'prop',{prop{1},prop{2},beta,prop{4}});
  cb_redraw(gca);
elseif strcmp(get(hObject,'string'),'+')
  beta=beta+0.05;
  set(findobj('string','+'),'Value',beta);
  setappdata(fig,'prop',{prop{1},prop{2},beta,prop{4}});
  cb_redraw(gca);
end

if strcmp(get(hObject,'string'),'<')
  setappdata(fig,'prop',{prop{1},slice-1,0,prop{4}});
  cb_redraw(gca);
end
if strcmp(get(hObject,'string'),'>')
  setappdata(fig,'prop',{prop{1},slice+1,0,prop{4}});
  cb_redraw(gca);
end
if strcmp(get(hObject,'string'),'OK')
  close(fig)
end

function keypress(h, eventdata, handles, varargin)
fig   = get(gca, 'parent');
prop  = getappdata(fig,'prop');
dat   = guidata(gcbf);
key   = get(gcbf, 'CurrentCharacter');
slice = prop{2};
if  key
  switch key
    case 28
      setappdata(fig,'prop',{prop{1},slice-1,0,prop{4}});
      cb_redraw(gca);
    case 29
      setappdata(fig,'prop',{prop{1},slice+1,0,prop{4}});
      cb_redraw(gca);
    otherwise
  end
end
uiresume(h)

function build_polygon
fig   = get(gca, 'parent');
prop  = getappdata(fig,'prop');
data  = getappdata(fig,'data');
ss = size(data);
proj  = prop{1};
slice = prop{2};
thispolygon =1;
polygon{thispolygon} = zeros(0,2);
maskhelp = [ ...
  '------------------------------------------------------------------------\n' ...
  'specify polygons for masking the topographic interpolation\n' ...
  'press the right mouse button to add another point to the current polygon\n' ...
  'press backspace on the keyboard to remove the last point\n' ...
  'press "c" on the keyboard to close this polygon and start with another\n' ...
  'press "q" or ESC on the keyboard to continue\n' ...
  ];
again = 1;
if ~ishold,hold,end
fprintf(maskhelp);
fprintf('\n');

while again
  for i=1:length(polygon)
    fprintf('polygon %d has %d points\n', i, size(polygon{i},1));
  end

  [x, y, k] = ginput(1);
  okflag = 0;

  % check points do not fall out of image boundaries
  if get(findobj('string','Ins'),'value')
    if (proj == 1 && y<ss(2) && x<ss(3) && x>1 && y>1)
      okflag = 1;
    elseif (proj == 2 && y<ss(1) && x<ss(3) && x>1 && y>1)
      okflag = 1;
    elseif (proj == 3 && y<ss(1) && x<ss(2) && x>1 && y>1)
      okflag = 1;
    else
      okflag = 0;
    end
  end

  if okflag
    switch lower(k)
      case 1
        polygon{thispolygon} = cat(1, polygon{thispolygon}, [x y]);
        % add the last line segment to the figure
        if size(polygon{thispolygon},1)>1
          x = polygon{i}([end-1 end],1);
          y = polygon{i}([end-1 end],2);
        end
        plot(x, y, 'g.-');

      case 8 % backspace
        if size(polygon{thispolygon},1)>0
          % redraw existing polygons
          cb_redraw(gca);
          % remove the last point of current polygon
          polygon{thispolygon} = polygon{thispolygon}(1:end-1,:);
          for i=1:length(polygon)
            x = polygon{i}(:,1);
            y = polygon{i}(:,2);
            if i~=thispolygon
              % close the polygon in the figure
              x(end) = x(1);
              y(end) = y(1);
            end
            set(gca,'nextplot','new')
            plot(x, y, 'g.-');
          end
        end

      case 'c'
        if size(polygon{thispolygon},1)>0
          % close the polygon
          polygon{thispolygon}(end+1,:) = polygon{thispolygon}(1,:);
          % close the polygon in the figure
          x = polygon{i}([end-1 end],1);
          y = polygon{i}([end-1 end],2);
          plot(x, y, 'g.-');
          % switch to the next polygon
          thispolygon = thispolygon + 1;
          polygon{thispolygon} = zeros(0,2);
          slicedata = getappdata(fig,'slicedata');
          m = length(slicedata);
          slicedata(m+1).proj    = prop{1};
          slicedata(m+1).slice   = prop{2};
          slicedata(m+1).polygon = polygon;
          setappdata(fig,'slicedata',slicedata);
        end

      case {'q', 27}
        if size(polygon{thispolygon},1)>0
          % close the polygon
          polygon{thispolygon}(end+1,:) = polygon{thispolygon}(1,:);
          % close the polygon in the figure
          x = polygon{i}([end-1 end],1);
          y = polygon{i}([end-1 end],2);
          plot(x, y, 'g.-');
        end
        again = 0;
        %set(gcf,'WindowButtonDownFcn','');
        slicedata = getappdata(fig,'slicedata');
        m = length(slicedata);
        if ~isempty(polygon{thispolygon})
          slicedata(m+1).proj    = prop{1};
          slicedata(m+1).slice   = prop{2};
          slicedata(m+1).polygon = polygon;
          setappdata(fig,'slicedata',slicedata);
        end
      otherwise
        ft_warning('invalid button (%d)', k);
    end

  end
  assign_points3d
end

function cb_btn_dwn(hObject, eventdata, handles)
build_polygon
uiresume

function cb_btn_dwn2(hObject, eventdata, handles)
global obj
erase_points2d(obj)
uiresume

function erase_points2d (obj)
fig       = get(gca, 'parent');
prop      = getappdata(fig,'prop');
slicedata = getappdata(fig,'slicedata');
pntsslice = [];

% 2d slice selected box boundaries
if strcmp(get(obj,'string'),'box')
  [x, y] = ft_select_box;
  rowmin = min(y);    rowmax = max(y);
  colmin = min(x);    colmax = max(x);
  box = getappdata(fig,'box');
  box = [box; prop{1} prop{2} rowmin rowmax colmin colmax];
  setappdata(fig,'box',box);
else
  % converts lines to points
  pos = slicedata2pnts(prop{1},prop{2});
  tmp = ft_select_point(pos);
  point2mark.x = tmp(1); point2mark.y = tmp(2);
  setappdata(fig,'point2mark',point2mark);
end

maskhelp = [ ...
  '------------------------------------------------------------------------\n' ...
  'Delete points:\n' ...
  'press "c" on the keyboard to undo all selections\n' ...
  'press "d" or DEL on the keyboard to delete all selections\n' ...
  ];
fprintf(maskhelp);
fprintf('\n');
again = 1;
if ~ishold,hold,end
cb_redraw(gca);

% waits for a key press
while again
  [junk1, junk2, k] = ginput(1);

  switch lower(k)
    case 'c'
      if strcmp(get(obj,'string'),'box')
        % undoes all the selections, but only in the current slice
        ind = [];
        for ii=1:size(box,1)
          if box(ii,1)==prop{1} && box(ii,2)==prop{2}
            ind = [ind,ii];
          end
        end
        box(ind,:) = [];
        setappdata(fig,'box',box);
        cb_redraw(gca);
        again = 0;
        set(gcf,'WindowButtonDownFcn','');
      else
        setappdata(fig,'point2mark',[]);
        cb_redraw(gca);
        again = 0;
        set(gcf,'WindowButtonDownFcn','');
      end
    case {'d', 127}
      if strcmp(get(obj,'string'),'box')
        update_slicedata
        % deletes only the points selected in the current slice
        ind = [];
        for ii=1:size(box,1)
          if box(ii,1)==prop{1} && box(ii,2)==prop{2}
            ind = [ind,ii];
          end
        end
        box(ind,:) = [];
        setappdata(fig,'box',box);
        cb_redraw(gca);
        again = 0;
        set(gcf,'WindowButtonDownFcn','');
      else
        update_slicedata
        setappdata(fig,'point2mark',[]);
        cb_redraw(gca);
        again = 0;
        set(gcf,'WindowButtonDownFcn','');
      end
    otherwise
  end
end
assign_points3d

function assign_points3d
fig       = get(gca, 'parent');
slicedata = getappdata(fig,'slicedata');

points = [];
for zz = 1:length(slicedata)
  slice   = slicedata(zz).slice;
  proj    = slicedata(zz).proj;
  polygon = slicedata(zz).polygon;
  for kk=1:length(polygon)
    if ~isempty(polygon{kk})
      [xs,ys]=deal(polygon{kk}(:,1),polygon{kk}(:,2));
      if (proj==1)     %jk
        tmp    = [ slice*ones(length(xs),1),ys(:),xs(:) ];
        points = [points; unique(tmp,'rows')];
      elseif (proj==2) %ik
        tmp    = [ ys(:),slice*ones(length(xs),1),xs(:) ];
        points = [points; unique(tmp,'rows')];
      elseif (proj==3) %ij
        tmp    = [ ys(:),xs(:),slice*ones(length(xs),1) ];
        points = [points; unique(tmp,'rows')];
      end
    end
  end
end
setappdata(fig,'points',points);

function update_slicedata
fig       = get(gca, 'parent');
prop      = getappdata(fig,'prop');
slicedata = getappdata(fig,'slicedata');
box       = getappdata(fig,'box');
point2mark = getappdata(fig,'point2mark');
% case of box selection
if ~isempty(box) && ~isempty(slicedata)
  for zz = 1:length(slicedata)
    slice   = slicedata(zz).slice;
    proj    = slicedata(zz).proj;
    polygon = slicedata(zz).polygon;
    for uu=1:size(box,1)
      if (box(uu,1)==proj) && (box(uu,2)==slice)
        for kk=1:length(polygon)
          if ~isempty(polygon{kk})
            [xs,ys] = deal(polygon{kk}(:,1),polygon{kk}(:,2));
            tmpind = lt(ys,ones(length(xs),1)*box(uu,4)).*gt(ys,ones(length(xs),1)*box(uu,3)).* ...
              lt(xs,ones(length(xs),1)*box(uu,6)).*gt(xs,ones(length(xs),1)*box(uu,5));
            slicedata(zz).polygon{kk} = [ xs(find(~tmpind)),ys(find(~tmpind)) ];
          end
        end
      end
    end
  end
  % case of point selection
else
  if ~isempty(slicedata)
    for zz = 1:length(slicedata)
      slice   = slicedata(zz).slice;
      proj    = slicedata(zz).proj;
      polygon = slicedata(zz).polygon;
      for kk=1:length(polygon)
        if ~isempty(polygon{kk})
          [xs,ys] = deal(polygon{kk}(:,1),polygon{kk}(:,2));
          tmpind = eq(xs,point2mark.x).*eq(ys,point2mark.y);
          slicedata(zz).polygon{kk} = [ xs(find(~tmpind)),ys(find(~tmpind)) ];
        end
      end
    end
  end
end
setappdata(fig,'slicedata',slicedata)

% FIXME: obsolete: replaced by headshape at the beginning
function smooth_mesh(hObject, eventdata, handles)
fig = get(gca,'parent');
disp('not yet implemented')

% FIXME: make it work for pre-loaded meshes
function view_mesh(hObject, eventdata, handles)
fig  = get(gca, 'parent');
assign_points3d
pnt  = getappdata(fig,'points');
pnt_ = pnt;
if ~isempty(pnt)
  tri = projecttri(pnt);
  bnd.pnt = pnt_;
  bnd.tri = tri;
  slicedata = getappdata(fig,'slicedata');
  figure
  ft_plot_mesh(bnd,'vertexcolor','k');
end

function cancel_mesh(hObject, eventdata, handles)
fig  = get(gca, 'parent');
close(fig)

function cb_close(hObject, eventdata, handles)
% get the points from the figure
fig    = hObject;
global pnt
pnt   = getappdata(fig, 'points');
set(fig, 'CloseRequestFcn', @delete);
delete(fig);

function cb_btnp1(hObject, eventdata, handles)
set(findobj('string','jk'),'value',1);
set(findobj('string','ik'),'value',0);
set(findobj('string','ij'),'value',0);
cb_btn(hObject);

function cb_btnp2(hObject, eventdata, handles)
set(findobj('string','jk'),'value',0);
set(findobj('string','ik'),'value',1);
set(findobj('string','ij'),'value',0);
cb_btn(hObject);

function cb_btnp3(hObject, eventdata, handles)
set(findobj('string','jk'),'value',0);
set(findobj('string','ik'),'value',0);
set(findobj('string','ij'),'value',1);
cb_btn(hObject);

function which_task1(hObject, eventdata, handles)
set(gcf,'WindowButtonDownFcn','');
set(findobj('string','View'),'value',1);
set(findobj('string','Ins'), 'value',0);
set(findobj('string','Del'), 'value',0);
set(gcf,'WindowButtonDownFcn','');
change_panel

function which_task2(hObject, eventdata, handles)
set(findobj('string','View'),'value',0);
set(findobj('string','Ins'), 'value',1);
set(findobj('string','Del'), 'value',0);
set(gcf,'WindowButtonDownFcn',@cb_btn_dwn);
change_panel

function which_task3(hObject, eventdata, handles)
set(gcf,'WindowButtonDownFcn','');
set(findobj('string','View'),'value',0);
set(findobj('string','Ins'), 'value',0);
set(findobj('string','Del'), 'value',1);
change_panel

function change_panel
% affect panel visualization
w1 = get(findobj('string','View'),'value');
w2 = get(findobj('string','Ins'), 'value');
w3 = get(findobj('string','Del'), 'value');

if w1
  set(findobj('tag','slice'),'string','slice');
  set(findobj('tag','min'),'string','<');
  set(findobj('tag','max'),'string','>');
  set(findobj('tag','min'),'style','pushbutton');
  set(findobj('tag','max'),'style','pushbutton');
  set(findobj('tag','min'),'callback',@cb_btn);
  set(findobj('tag','max'),'callback',@cb_btn);
  set(findobj('tag','mesh'),'visible','on');
  set(findobj('tag','openm'),'visible','on');
  set(findobj('tag','viewm'),'visible','on');
  set(findobj('tag','brightness'),'visible','on');
  set(findobj('tag','plus'),'visible','on');
  set(findobj('tag','minus'),'visible','on');
  set(findobj('tag','toggle axes'),'visible','on');
  set(findobj('tag','plane'),'visible','on');
  set(findobj('tag','one'),'visible','on');
  set(findobj('tag','two'),'visible','on');
  set(findobj('tag','three'),'visible','on');
elseif w2
  set(findobj('tag','slice'),'string','select points');
  set(findobj('tag','min'),'string','close');
  set(findobj('tag','max'),'string','next');
  set(findobj('tag','min'),'style','pushbutton');
  set(findobj('tag','max'),'style','pushbutton');
  set(findobj('tag','min'),'callback',@tins_close);
  set(findobj('tag','max'),'callback',@tins_next);
  set(findobj('tag','mesh'),'visible','off');
  set(findobj('tag','openm'),'visible','off');
  set(findobj('tag','viewm'),'visible','off');
  set(findobj('tag','brightness'),'visible','off');
  set(findobj('tag','plus'),'visible','off');
  set(findobj('tag','minus'),'visible','off');
  set(findobj('tag','toggle axes'),'visible','off');
  set(findobj('tag','plane'),'visible','off');
  set(findobj('tag','one'),'visible','off');
  set(findobj('tag','two'),'visible','off');
  set(findobj('tag','three'),'visible','off');
elseif w3
  set(findobj('tag','slice'),'string','select points');
  set(findobj('tag','min'),'string','box');
  set(findobj('tag','max'),'string','point');
  set(findobj('tag','min'),'style','pushbutton');
  set(findobj('tag','max'),'style','pushbutton');
  set(findobj('tag','min'),'callback',@tdel_box);
  set(findobj('tag','max'),'callback',@tdel_single);
  set(findobj('tag','mesh'),'visible','off');
  set(findobj('tag','openm'),'visible','off');
  set(findobj('tag','viewm'),'visible','off');
  set(findobj('tag','brightness'),'visible','off');
  set(findobj('tag','plus'),'visible','off');
  set(findobj('tag','minus'),'visible','off');
  set(findobj('tag','toggle axes'),'visible','off');
  set(findobj('tag','plane'),'visible','off');
  set(findobj('tag','one'),'visible','off');
  set(findobj('tag','two'),'visible','off');
  set(findobj('tag','three'),'visible','off');
end

function tins_close(hObject, eventdata, handles)
set(gcf,'WindowButtonDownFcn','');

function tins_next(hObject, eventdata, handles)
set(gcf,'WindowButtonDownFcn','');

function tdel_box(hObject, eventdata, handles)
global obj
obj = gco;
set(gcf,'WindowButtonDownFcn',@cb_btn_dwn2);

function tdel_single(hObject, eventdata, handles)
global obj
obj = gco;
set(gcf,'WindowButtonDownFcn',@cb_btn_dwn2);

% accessory functions:
function [pnts] = slicedata2pnts(proj,slice)
fig       = get(gca, 'parent');
slicedata = getappdata(fig,'slicedata');
pnts = [];
for i=1:length(slicedata)
  if slicedata(i).slice==slice && slicedata(i).proj == proj
    for j=1:length(slicedata(i).polygon)
      if ~isempty(slicedata(i).polygon{j})
        tmp = slicedata(i).polygon{j};
        xs = tmp(:,1); ys = tmp(:,2);
        pnts = [pnts;[unique([xs,ys],'rows')]];
      end
    end
  end
end

function h = ispolygon(x)
h = length(x.XData)>1;

function h = ispoint(x)
h = length(x.XData)==1;

function [T,P] = points2param(pnt)
x = pnt(:,1); y = pnt(:,2); z = pnt(:,3);
[phi,theta,R] = cart2sph(x,y,z);
theta = theta + pi/2;
phi   = phi + pi;
T     = theta(:);
P     = phi(:);

function [bnd2,a,Or] = spherical_harmonic_mesh(bnd, nl)

% SPHERICAL_HARMONIC_MESH realizes the smoothed version of a mesh contained in the first
% argument bnd. The boundary argument (bnd) contains typically 2 fields
% called .pnt and .tri referring to vertices and triangulation of a mesh.
% The degree of smoothing is given by the order in the second input: nl
%
% Use as
%   bnd2 = spherical_harmonic_mesh(bnd, nl);

% nl    = 12; % from van 't Ent paper

bnd2 = [];

% calculate midpoint
Or = mean(bnd.pnt);
% rescale all points
x = bnd.pnt(:,1) - Or(1);
y = bnd.pnt(:,2) - Or(2);
z = bnd.pnt(:,3) - Or(3);
X = [x(:), y(:), z(:)];

% convert points to parameters
[T,P] = points2param(X);

% basis function
B     = shlib_B(nl, T, P);
Y     = B'*X;

% solve the linear system
a     = pinv(B'*B)*Y;

% build the surface
Xs = zeros(size(X));
for l = 0:nl-1
  for m = -l:l
    if m<0
      Yml = shlib_Yml(l, abs(m), T, P);
      Yml = (-1)^m*conj(Yml);
    else
      Yml = shlib_Yml(l, m, T, P);
    end
    indx = l^2 + l + m+1;
    Xs = Xs + Yml*a(indx,:);
  end
  fprintf('%d of %d\n', l, nl);
end
% Just take the real part
Xs = real(Xs);
% reconstruct the matrices for plotting.
xs = reshape(Xs(:,1), size(x,1), size(x,2));
ys = reshape(Xs(:,2), size(x,1), size(x,2));
zs = reshape(Xs(:,3), size(x,1), size(x,2));

bnd2.pnt = [xs(:)+Or(1) ys(:)+Or(2) zs(:)+Or(3)];
[bnd2.tri] = projecttri(bnd.pnt);

function B = shlib_B(nl, theta, phi)
% function B = shlib_B(nl, theta, phi)
%
% Constructs the matrix of basis functions
% where b(i,l^2 + l + m) = Yml(theta, phi)
%
% See also: shlib_Yml.m, shlib_decomp.m, shlib_gen_shfnc.m
%
% Dr. A. I. Hanna (2006).
B = zeros(length(theta), nl^2+2*nl+1);
for l = 0:nl-1
  for m = -l:l
    if m<0
      Yml = shlib_Yml(l, abs(m), theta, phi);
      Yml = (-1)^m*conj(Yml);
    else
      Yml = shlib_Yml(l, m, theta, phi);
    end
    indx = l^2 + l + m+1;
    B(:, indx) = Yml;
  end
end
return;

function Yml = shlib_Yml(l, m, theta, phi)
% function Yml = shlib_Yml(l, m, theta, phi)
%
% MATLAB function that takes a given order and degree, and the matrix of
% theta and phi and constructs a spherical harmonic from these. The
% analogue in the 1D case would be to give a particular frequency.
%
% Inputs:
%  m - order of the spherical harmonic
%  l - degree of the spherical harmonic
%  theta - matrix of polar coordinates \theta \in [0, \pi]
%  phi - matrix of azimuthal cooridinates \phi \in [0, 2\pi)
%
% Example:
%
% [x, y] = meshgrid(-1:.1:1);
% z = x.^2 + y.^2;
% [phi,theta,R] = cart2sph(x,y,z);
% Yml = spharm_aih(2,2, theta(:), phi(:));
%
% See also: shlib_B.m, shlib_decomp.m, shlib_gen_shfnc.m
%
% Dr. A. I. Hanna (2006)
Pml=legendre(l,cos(theta));
if l~=0
  Pml=squeeze(Pml(m+1,:,:));
end
Pml = Pml(:);
% Yml = sqrt(((2*l+1)/4).*(factorial(l-m)/factorial(l+m))).*Pml.*exp(sqrt(-1).*m.*phi);
% new:
Yml = sqrt(((2*l+1)/(4*pi)).*(factorial(l-m)/factorial(l+m))).*Pml.*exp(sqrt(-1).*m.*phi);
return;


text:
function pdf = burr_pdf ( x, a, b, c, d )

%*****************************************************************************80
%
%% BURR_PDF evaluates the Burr PDF.
%
%  Discussion:
%
%    PDF(X)(A,B,C,D) = ( C * D / B ) * ( ( X - A ) / B )**( - C - 1 )
%      * ( 1 + ( ( X - A ) / B )**( - C ) )**( - D - 1 )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    M E Johnson,
%    Multivariate Statistical Simulation,
%    Wiley, New York, 1987.
%
%  Parameters:
%
%    Input, real X, the argument of the PDF.
%    A <= X
%
%    Input, real A, B, C, D, the parameters of the PDF.
%    0 < B,
%    0 < C.
%
%    Output, real PDF, the value of the PDF.
%
  if ( x <= a )
    pdf = 0.0;
  else

    y = ( x - a ) / b;

    pdf = ( c * d / b ) * y^( - c - 1.0 ) * ( 1.0 + y^( -c ) )^( - d - 1.0 );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   demo script to insert new nodes to a tetrahedral mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load rat_head.mat

[node,elem,face]=vol2mesh(volimage>0.05,1:size(volimage,1),1:size(volimage,2),...
                           1:size(volimage,3),2,2,1);

plotmesh(node,face);

dt=pi/10;
r=3;
x0=30;
y0=28;
theta=dt:dt:2*pi;
x=x0+r*cos(theta);
y=y0+r*sin(theta);

p0=[x;y;ones(1,length(theta))*50]';
v0=[0,0,-1];

[t,u,v,idx,xnode]=raysurf(p0,v0,node,face(:,1:3));

hold on
plotmesh(p0,'r+');

goodpt=find(~isnan(xnode(:,1)) & ~isnan(xnode(:,2)) & ~isnan(xnode(:,3)));
[newnode,newelem,newface]=meshrefine(node,elem,face,xnode(goodpt,:));

figure
plotmesh(newnode,newface);
hold on
plotmesh(p0,'r+');
plotmesh(xnode,'r.')

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function [R,eff] = randmio_und_signed(W, ITER)
% RANDMIO_UND_SIGNED	Random graph with preserved signed degree distribution
%
%   R       = randmio_und_signed(W,ITER);
%   [R,eff] = randmio_und_signed(W,ITER);
%
%   This function randomizes an undirected network with positively and
%   negatively signed connections, while preserving the positively and
%   negatively signed degree distribution. The function does not preserve
%   the strength distribution in weighted networks.
%
%   Input:      W,      undirected (binary/weighted) connection matrix
%               ITER,   rewiring parameter
%                       (each edge is rewired approximately ITER times)
%
%   Output:     R,      randomized network
%               eff,    number of actual rewirings carried out
%
%   Reference:  Maslov and Sneppen (2002) Science 296:910
%
%
%   2011-2015
%   Dani Bassett, UCSB
%   Olaf Sporns,  Indiana U
%   Mika Rubinov, U Cambridge

%   Modification History:
%   Mar 2011: Original (Dani Bassett, based on randmio_und.m)
%   Mar 2012: Limit number of rewiring attempts,
%             count number of successful rewirings (Olaf Sporns)
%   Dec 2015: Rewritten the core of the rewiring algorithm to allow
%             unbiased exploration of all network configurations. The new
%             algorithm allows positive-positive/negative-negative
%             rewirings, in addition to the previous positive-positive/0-0
%             and negative-negative/0-0 rewirings (Mika Rubinov).

if nargin('randperm')==1
    warning('This function requires a recent (>2011) version of MATLAB.')
end

R     = double(W);              % sign function requires double input
n     = size(R,1);
ITER  = ITER*n*(n-1)/2;

% maximal number of rewiring attempts per 'iter'
maxAttempts = round(n/2);
% actual number of successful rewirings
eff = 0;

for iter=1:ITER
    att=0;
    while (att<=maxAttempts)    %while not rewired
        %select four distinct vertices
        nodes = randperm(n,4);
        a = nodes(1);
        b = nodes(2);
        c = nodes(3);
        d = nodes(4);

        r0_ab = R(a,b);
        r0_cd = R(c,d);
        r0_ad = R(a,d);
        r0_cb = R(c,b);

        %rewiring condition
        if      (sign(r0_ab)==sign(r0_cd)) && ...
                (sign(r0_ad)==sign(r0_cb)) && ...
                (sign(r0_ab)~=sign(r0_ad))

            R(a,d)=r0_ab; R(a,b)=r0_ad;
            R(d,a)=r0_ab; R(b,a)=r0_ad;
            R(c,b)=r0_cd; R(c,d)=r0_cb;
            R(b,c)=r0_cd; R(d,c)=r0_cb;

            eff = eff+1;
            break;
        end %rewiring condition
        att=att+1;
    end %while not rewired
end %iterations

author: canlab

repo: CanlabCore

save_path: github-repos/MATLAB/canlab-CanlabCore


text:
function sphere_cubed_grid_points_display ( ns, xyz, filename )

%*****************************************************************************80
%
%% SPHERE_CUBED_GRID_POINTS_DISPLAY displays the points on a cubed sphere grid.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 May 2015
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer NS, the number of points.
%
%    Output, real XYZ(NS,3), distinct points on the unit sphere
%    generated by a cubed sphere grid.
%
  figure ( )
  clf
  hold on
  [ x, y, z ] = sphere ( 20 );
  c = ones ( size ( z ) );
  surf ( x, y, z, c );
  plot3 ( xyz(:,1), xyz(:,2), xyz(:,3), 'b.', 'Markersize', 20 );
  axis equal
  grid on
  view ( 3 )
  xlabel ( '<--X-->' )
  ylabel ( '<--Y-->' )
  zlabel ( '<--Z-->' )
  title_string = s_escape_tex ( filename )
  title ( title_string, 'FontSize', 24 );
  hold off

  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Plot file saved to "%s".\n', filename );

  return
end


text:
function logistic
% Logistic growth
%    using MATLAB for analytical solution
%
%   $Ekkehard Holzbecher  $Date: 2006/04/20 $
%--------------------------------------------------------------------------
T = 10;                  % maximum time
r = 1;                   % rate
kappa = 1;               % capacity
c0 = 0.01;               % initial value

%----------------------execution-------------------------------------------

t = linspace (0,T,100);
e = exp(r*t);
c = c0*kappa*e./(kappa+c0*(e-1));

%---------------------- graphical output ----------------------------------

plot (t,c); grid;
xlabel ('time'); legend ('population');
title ('logistic growth');

author: Sable

repo: mcbench-benchmarks


text:
classdef DOC4 < PROBLEM
% <multi> <real> <constrained>
% Benchmark MOP with constraints in both decision and objective spaces

%------------------------------- Reference --------------------------------
% Z. Liu and Y. Wang, Handling constrained multiobjective optimization
% problems with constraints in both the decision and objective spaces. IEEE
% Transactions on Evolutionary Computation, 2019, 23(5): 870-884.
%------------------------------- Copyright --------------------------------
% Copyright (c) 2023 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

    methods
        %% Default settings of the problem
        function Setting(obj)
            obj.M = 2;
            obj.D = 8;
            obj.lower    = [0 -10 -10 -10 -10 -10 -10 -10];
            obj.upper    = [ 1 10 10 10 10 10 10 10];
            obj.encoding = ones(1,obj.D);
        end
        %% Calculate objective values and constraint violations
        function Population = Evaluation(obj,varargin)
            X = varargin{1};
            X = max(min(X,repmat(obj.upper,size(X,1),1)),repmat(obj.lower,size(X,1),1));
            g_temp = (X(:, 2) - 10).^2 + 5 * (X(:, 3) - 12).^2 + X(:, 4).^4 + 3 * (X(:, 5) - 11).^2 + 10 * X(:, 6).^6 + ...
            7 * X(:, 7).^2 + X(:, 8).^4 - 4 * X(:, 7).* X(:, 8) - 10 * X(:, 7) - 8 * X(:, 8);
            g = g_temp-680.6300573745 +1;
            PopObj(:,1) = X(:,1);
            PopObj(:,2) = g.*(1-sqrt(PopObj(:,1))./g);
            % Constraints in objective space
            c(:,1) = max( -(PopObj(:,1) + PopObj(:,2)-1), 0);
            c(:,2) = max(- ( PopObj(:,1)+ PopObj(:,2) - 1 - abs(sin(10*pi*(PopObj(:,1) - PopObj(:,2) + 1) ))), 0);
            % Constraints in decision space
            c(:,3) = -127 + 2 * X(:, 2).^2 + 3 * X(:, 3).^4 + X(:, 4) + 4 * X(:, 5).^2 + 5 * X(:, 6);
            c(:,4) = -282 + 7 * X(:, 2) + 3 * X(:, 3) + 10 * X(:, 4).^2 + X(:, 5) - X(:, 6);
            c(:,5) = -196 + 23 * X(:, 2) + X(:, 3).^2 + 6 * X(:, 7).^2 - 8 * X(:, 8);
            c(:,6) = 4 * X(:, 2).^2 + X(:, 3).^2 - 3 * X(:, 2).* X(:, 3) + 2 * X(:, 4).^2 + 5 * X(:, 7) - 11 * X(:, 8);
            Population = SOLUTION(X,PopObj,c,varargin{2:end});
            obj.FE     = obj.FE + length(Population);
        end
        %% Generate points on the Pareto front
        function R = GetOptimum(obj,N)
            R(:,1) = (0:20)/20;
            R(:,2) = 1 - R(:,1);
        end
    end
end

author: BIMK

repo: PlatEMO

save_path: github-repos/MATLAB/BIMK-PlatEMO


text:
function voxelThicknessV = deduceVoxelThicknesses(scanNum, planC)
%"deduceVoxelThicknesses"
%   Deduces voxel thicknesses when the treatment planning system provides
%   only zValues.
%
%   Note:  these voxel thicknesses may well be different than CT slice
%   thicknesses as read out by the CT scanners!  This is because slice
%   thickness on a CT scanner is not determined by slice spacing necessarily,
%   whereas we are assigning voxels absolute volumes which they are to
%   represent in 3-D.  Therefore, our slice thicknesses need to be consistent
%   with z values, and every point in the CT scan needs to be assigned to
%   a unique voxel.
%
%   The planC argument is optional.  If not specified, the global planC is
%   used.
%
%Created            :  30 Apr 03, JOD
%Latest modifications: 08 May 03, JOD, check for thicknesses < 0.
%                      28 May 03, JOD, raised tolerance for consistency
%                                      check to 10000 * eps due to failure
%                                      at 100 * eps.
%                      05 Sep 03, JOD, changed assignment of thickness for
%                                       first and last slices: previously
%                                       not consistent and produced an error.
%                      02 Feb 16, JRA, Added multiscan support.
%Usage:
%   function voxelThicknessV = deduceVoxelThicknesses(scanNum, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

if ~exist('planC')
    global planC
end

indexS = planC{end};

zValuesV = [planC{indexS.scan}(scanNum).scanInfo(:).zValue];

voxelThicknessV = ones(size(zValuesV)) * NaN;

for i = 2 : length(zValuesV) - 1

  nextDelta = abs(zValuesV(i+1) - zValuesV(i));

  lastDelta = abs(zValuesV(i) - zValuesV(i-1));

  if nextDelta == lastDelta

    voxelThicknessV(i) = nextDelta;

  else

    %split thicknesses:
    voxelThicknessV(i) = 0.5 * lastDelta + 0.5 * nextDelta;

  end

end

if length(zValuesV) > 1
    voxelThicknessV(1) = abs(zValuesV(2) - zValuesV(1));  %JOD, 5 Sept 03

    voxelThicknessV(end) = abs(zValuesV(end) - zValuesV(end - 1)); %JOD, 5 Sept 03
else
    %voxelThicknessV = optS.alternateLimitUniformCTSliceSpacing;
    voxelThicknessV = 1; % dummy value for single slice
end

%Check

if any(isnan(voxelThicknessV))
  error('Error in determining slice thicknesses: not all slices were assigned thicknesses')
end

if length(zValuesV)>1 && abs(sum(voxelThicknessV) - (abs(zValuesV(end)-zValuesV(1)) + 0.5 * voxelThicknessV(1) + 0.5 * voxelThicknessV(end))) > 10000 * eps
  error('Voxel thicknesses inconsistent with z values.')
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To run:  Edit user options in npscalc.m and execute the code.
%
% This code is intended to accompany the manuscript entitled:
%
% A simple approach to measure computed tomography (CT) modulation transfer
% function (MTF) and noise-power spectrum (NPS) using the American College
% of Radiology (ACR) accreditation phantom
%
% Please do not distribute.
%
% This program requires the following 4 files:
%   npscalc.m
%   SliceBrowser.m
%   SliceBrowser.fig
%   license.txt
%
% Purpose:  To calculate the 3D NPS using CT data of the American College
%           of Radiology (ACR) accreditation phantom.
%
% Input:    Two consecutive scans of the phantom are needed.  The program
%           requires a data directory to be selected in which there are
%           only two subdirectories containing only CT slices corresponding
%           to the third module of the phantom.  Be careful of partial
%           volume effects with surrounding modules.
%
%           i.e.,  datadir
%                      |-> scan 1 dir
%                      |           | -> only module 3 slices
%                      |
%                      |-> scan 2 dir
%                                  |-> only module 3 slices
%
% Copyright 2012 Saul N. Friedman
%   Distributed under the terms of the "New BSD License."  Please see
%   license.txt.
%
% N.B.:  SliceBrowser.m is an altered version of code by Marian Uhercik.
%        Please see license.txt.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ======================================================================
%> SLICEBROWSER M-file for SliceBrowser.fig
%>       SliceBrowser is an interactive viewer of 3D volumes,
%>       it shows 3 perpendicular slices (XY, YZ, ZX) with 3D pointer.
%>   Input:  a) VOLUME - a 3D matrix with volume data
%>           b) VOLUME - a 4D matrix with volume data over time
%>   Control:
%>       - Clicking into the window changes the location of 3D pointer.
%>       - 3D pointer can be moved also by keyboard arrows.
%>       - Pressing +/- will switch to next/previous volume.
%>       - Pressing 1,2,3 will change the focus of current axis.
%>       - Pressing 'e' will print the location of 3D pointer.
%>       - Pressing 'c' switches between color-mode and grayscale.
%>       - Pressing 'q' switches scaling of axis (equal/normal).
%>   Example of usage:
%>       load mri.dat
%>       volume = squeeze(D);
%>       SliceBrowser(volume);
%>
%> Last modified by Saul N. Friedman
% ======================================================================
function varargout = SliceBrowser(varargin)

% Documentation generated GUIDE:
%
%SLICEBROWSER M-file for SliceBrowser.fig
%      SLICEBROWSER, by itself, creates a new SLICEBROWSER or raises the existing
%      singleton*.
%
%      H = SLICEBROWSER returns the handle to a new SLICEBROWSER or the handle to
%      the existing singleton*.
%
%      SLICEBROWSER('Property','Value',...) creates a new SLICEBROWSER using the
%      given property value pairs. Unrecognized properties are passed via
%      varargin to SliceBrowser_OpeningFcn.  This calling syntax produces a
%      warning when there is an existing singleton*.
%
%      SLICEBROWSER('CALLBACK') and SLICEBROWSER('CALLBACK',hObject,...) call the
%      local function named CALLBACK in SLICEBROWSER.M with the given input
%      arguments.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

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

% Last Modified by GUIDE v2.5 01-Sep-2008 13:14:20

% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @SliceBrowser_OpeningFcn, ...
                   'gui_OutputFcn',  @SliceBrowser_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 SliceBrowser is made visible.
function SliceBrowser_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   unrecognized PropertyName/PropertyValue pairs from the
%            command line (see VARARGIN)

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

% Update handles structure
guidata(hObject, handles);

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

if (length(varargin) <=0)
    error('Input volume has not been specified.');
end;
volume = varargin{1};
if (ndims(volume) ~= 3 && ndims(volume) ~= 4)
    error('Input volume must have 3 or 4 dimensions.');
end;
handles.volume = volume;

handles.axis_equal = 0;
handles.color_mode = 1;
if (size(volume,4) ~= 3)
    handles.color_mode = 0;
end;

% set main wnd title
set(gcf, 'Name', 'Slice Viewer')

% init 3D pointer
vol_sz = size(volume);
if (ndims(volume) == 3)
    vol_sz(4) = 1;
end;
pointer3dt = floor(vol_sz/2)+1;
handles.pointer3dt = pointer3dt;
handles.vol_sz = vol_sz;

plot3slices(hObject, handles);

% stores ID of last axis window
% (0 means that no axis was clicked yet)
handles.last_axis_id = 0;

% Update handles structure
guidata(hObject, handles);

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

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

% --- Executes on mouse press over axes background.
function Subplot1_ButtonDownFcn(hObject, eventdata, handles)
% hObject    handle to Subplot1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% This object contains the XY slice

xaxis = getappdata(0,'uaxis');
yaxis = getappdata(0,'vaxis');
zaxis = getappdata(0,'waxis');
axlabel = getappdata(0,'axlabel');
deltayaxis = yaxis(2) - yaxis(1);
deltaxaxis = xaxis(2) - xaxis(1);
deltazaxis = zaxis(2) - zaxis(1);

%disp('Subplot1:BtnDown');
pt=get(gca,'currentpoint');
if length(xaxis) > 0 && length(yaxis)
    xpos=round((pt(1,2)+xaxis(end))/deltaxaxis); ypos=round((pt(1,1)+yaxis(end))/deltayaxis);
else
xpos=round(pt(1,2)); ypos=round(pt(1,1));
end
zpos = handles.pointer3dt(3);
tpos = handles.pointer3dt(4);
handles.pointer3dt = [xpos ypos zpos tpos];
handles.pointer3dt = clipointer3d(handles.pointer3dt,handles.vol_sz);
plot3slices(hObject, handles);
% store this axis as last clicked region
handles.last_axis_id = 1;
% Update handles structure
guidata(hObject, handles);

% --- Executes on mouse press over axes background.
function Subplot2_ButtonDownFcn(hObject, eventdata, handles)
% hObject    handle to Subplot2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% This object contains the YZ slice

xaxis = getappdata(0,'uaxis');
yaxis = getappdata(0,'vaxis');
zaxis = getappdata(0,'waxis');
axlabel = getappdata(0,'axlabel');
deltayaxis = yaxis(2) - yaxis(1);
deltaxaxis = xaxis(2) - xaxis(1);
deltazaxis = zaxis(2) - zaxis(1);

%disp('Subplot2:BtnDown');
pt=get(gca,'currentpoint');
if length(xaxis) > 0 && length(zaxis)
    xpos=round((pt(1,2)+xaxis(end))/deltaxaxis); zpos=round((pt(1,1)+zaxis(end))/deltazaxis);
else
xpos=round(pt(1,2)); zpos=round(pt(1,1));
end
ypos = handles.pointer3dt(2);
tpos = handles.pointer3dt(4);
handles.pointer3dt = [xpos ypos zpos tpos];
handles.pointer3dt = clipointer3d(handles.pointer3dt,handles.vol_sz);
plot3slices(hObject, handles);
% store this axis as last clicked region
handles.last_axis_id = 2;
% Update handles structure
guidata(hObject, handles);

% --- Executes on mouse press over axes background.
function Subplot3_ButtonDownFcn(hObject, eventdata, handles)
% hObject    handle to Subplot3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% This object contains the XZ slice

xaxis = getappdata(0,'uaxis');
yaxis = getappdata(0,'vaxis');
zaxis = getappdata(0,'waxis');
axlabel = getappdata(0,'axlabel');
deltayaxis = yaxis(2) - yaxis(1);
deltaxaxis = xaxis(2) - xaxis(1);
deltazaxis = zaxis(2) - zaxis(1);

%disp('Subplot3:BtnDown');
pt=get(gca,'currentpoint');
if length(xaxis) > 0 && length(zaxis)
    ypos=round((pt(1,1)+yaxis(end))/deltayaxis); zpos=round((pt(1,2)+zaxis(end))/deltazaxis);
else
zpos=round(pt(1,2)); ypos=round(pt(1,1));
end
xpos = handles.pointer3dt(1);
tpos = handles.pointer3dt(4);
handles.pointer3dt = [xpos ypos zpos tpos];
handles.pointer3dt = clipointer3d(handles.pointer3dt,handles.vol_sz);
plot3slices(hObject, handles);
% store this axis as last clicked region
handles.last_axis_id = 3;
% Update handles structure
guidata(hObject, handles);

% --- Executes on key press with focus on SliceBrowserFigure and no controls selected.
function SliceBrowserFigure_KeyPressFcn(hObject, eventdata, handles)
% hObject    handle to SliceBrowserFigure (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

%disp('SliceBrowserFigure_KeyPressFcn');
curr_char = int8(get(gcf,'CurrentCharacter'));
if isempty(curr_char)
    return;
end;

xpos = handles.pointer3dt(1);
ypos = handles.pointer3dt(2);
zpos = handles.pointer3dt(3);
tpos = handles.pointer3dt(4);
% Keys:
% - up:   30
% - down:   31
% - left:   28
% - right:   29
% - '1': 49
% - '2': 50
% - '3': 51
% - 'e': 101
% - plus:  43
% - minus:  45
switch curr_char
    case 99 % 'c'
        handles.color_mode = 1 - handles.color_mode;
        if (handles.color_mode ==1 && size(handles.volume,4) ~= 3)
            handles.color_mode = 0;
        end;

    case 113 % 'q'
        handles.axis_equal = 1 - handles.axis_equal;

    case 30
        switch handles.last_axis_id
            case 1
                xpos = xpos -1;
            case 2
                xpos = xpos -1;
            case 3
                zpos = zpos -1;
            case 0
        end;
    case 31
        switch handles.last_axis_id
            case 1
                xpos = xpos +1;
            case 2
                xpos = xpos +1;
            case 3
                zpos = zpos +1;
            case 0
        end;
    case 28
        switch handles.last_axis_id
            case 1
                ypos = ypos -1;
            case 2
                zpos = zpos -1;
            case 3
                ypos = ypos -1;
            case 0
        end;
    case 29
        switch handles.last_axis_id
            case 1
                ypos = ypos +1;
            case 2
                zpos = zpos +1;
            case 3
                ypos = ypos +1;
            case 0
        end;
    case 43
        % plus key
        tpos = tpos+1;
    case 45
        % minus key
        tpos = tpos-1;
    case 49
        % key 1
        handles.last_axis_id = 1;
    case 50
        % key 2
        handles.last_axis_id = 2;
    case 51
        % key 3
        handles.last_axis_id = 3;
    case 101
        disp(['[' num2str(xpos) ' ' num2str(ypos) ' ' num2str(zpos) ' ' num2str(tpos) ']']);
    otherwise
        return
end;
handles.pointer3dt = [xpos ypos zpos tpos];
handles.pointer3dt = clipointer3d(handles.pointer3dt,handles.vol_sz);
plot3slices(hObject, handles);
% Update handles structure
guidata(hObject, handles);

% --- Plots all 3 slices XY, YZ, XZ into 3 subplots
function [sp1,sp2,sp3] = plot3slices(hObject, handles)
% pointer3d     3D coordinates in volume matrix (integers)

handles.pointer3dt;
size(handles.volume);
value3dt = handles.volume(handles.pointer3dt(1), handles.pointer3dt(2), handles.pointer3dt(3), handles.pointer3dt(4));

text_str = ['[X:' int2str(handles.pointer3dt(1)) ...
           ', Y:' int2str(handles.pointer3dt(2)) ...
           ', Z:' int2str(handles.pointer3dt(3)) ...
           ', Time:' int2str(handles.pointer3dt(4)) '/' int2str(handles.vol_sz(4)) ...
           '], value:' num2str(value3dt)];
set(handles.pointer3d_info, 'String', text_str);
guidata(hObject, handles);

if (handles.color_mode ==1)
    sliceXY = squeeze(handles.volume(:,:,handles.pointer3dt(3),:));
    sliceYZ = squeeze(handles.volume(handles.pointer3dt(1),:,:,:));
    sliceXZ = squeeze(handles.volume(:,handles.pointer3dt(2),:,:));

    max_xyz = max([ max(sliceXY(:)) max(sliceYZ(:)) max(sliceXZ(:)) ]);
    min_xyz = min([ min(sliceXY(:)) min(sliceYZ(:)) min(sliceXZ(:)) ]);
    clims = [ min_xyz max_xyz ];
else
    sliceXY = squeeze(handles.volume(:,:,handles.pointer3dt(3),handles.pointer3dt(4)));
    sliceYZ = squeeze(handles.volume(handles.pointer3dt(1),:,:,handles.pointer3dt(4)));
    sliceXZ = squeeze(handles.volume(:,handles.pointer3dt(2),:,handles.pointer3dt(4)));

    max_xyz = max([ max(sliceXY(:)) max(sliceYZ(:)) max(sliceXZ(:)) ]);
    min_xyz = min([ min(sliceXY(:)) min(sliceYZ(:)) min(sliceXZ(:)) ]);
    clims = [ min_xyz max_xyz ];
end;
sliceZY = squeeze(permute(sliceYZ, [2 1 3]));

xaxis = getappdata(0,'uaxis');
yaxis = getappdata(0,'vaxis');
zaxis = getappdata(0,'waxis');
axlabel = getappdata(0,'axlabel');
if length(xaxis) > 0 && length(yaxis) && length(zaxis)
deltayaxis = yaxis(2) - yaxis(1);
deltaxaxis = xaxis(2) - xaxis(1);
deltazaxis = zaxis(2) - zaxis(1);
end

sp1 = subplot(2,2,1);
%colorbar;
if length(xaxis) > 0 && length(yaxis)
    imagesc([yaxis(1) yaxis(end)],[xaxis(1) xaxis(end)],sliceXY,clims)
else
imagesc(sliceXY, clims);
end
if (handles.axis_equal == 1)
    axis image;
else
    axis normal;
end;
title('Slice XY');
if exist('axlabel')
    xlabel(['Y',axlabel]);ylabel(['X',axlabel])
else
ylabel('X');xlabel('Y');
end
if length(xaxis) > 0 && length(yaxis)
    line([handles.pointer3dt(2).*deltayaxis - yaxis(end) handles.pointer3dt(2).*deltayaxis - yaxis(end)], [-size(handles.volume,1) size(handles.volume,1)]);
    line([-size(handles.volume,2) size(handles.volume,2)], [handles.pointer3dt(1).*deltaxaxis - xaxis(end) handles.pointer3dt(1).*deltaxaxis - xaxis(end)]);
else
line([handles.pointer3dt(2) handles.pointer3dt(2)], [0 size(handles.volume,1)]);
line([0 size(handles.volume,2)], [handles.pointer3dt(1) handles.pointer3dt(1)]);
end
%set(allchild(gca),'ButtonDownFcn',@Subplot1_ButtonDownFcn);
set(allchild(gca),'ButtonDownFcn','SliceBrowser(''Subplot1_ButtonDownFcn'',gca,[],guidata(gcbo))');

sp2 = subplot(2,2,2);
if length(xaxis) > 0 &&length(zaxis) > 0
   imagesc([zaxis(1) zaxis(end)],[xaxis(1) xaxis(end)],sliceXZ,clims)
else
imagesc(sliceXZ, clims);
end
if (handles.axis_equal == 1)
    axis image;
else
    axis normal;
end;
title('Slice XZ');
if exist('axlabel')
    ylabel(['X',axlabel]);xlabel(['Z',axlabel])
else
ylabel('X');xlabel('Z');
end
if length(xaxis) > 0 && length(zaxis)
    line([handles.pointer3dt(3).*deltazaxis - zaxis(end) handles.pointer3dt(3).*deltazaxis - zaxis(end)], [-size(handles.volume,1).*deltaxaxis size(handles.volume,1).*deltaxaxis]);
line([-size(handles.volume,3).*deltazaxis size(handles.volume,3).*deltazaxis], [handles.pointer3dt(1).*deltaxaxis - xaxis(end) handles.pointer3dt(1).*deltaxaxis - xaxis(end)]);
else
line([handles.pointer3dt(3) handles.pointer3dt(3)], [0 size(handles.volume,1)]);
line([0 size(handles.volume,3)], [handles.pointer3dt(1) handles.pointer3dt(1)]);
end
%set(allchild(gca),'ButtonDownFcn',@Subplot2_ButtonDownFcn);
set(allchild(gca),'ButtonDownFcn','SliceBrowser(''Subplot2_ButtonDownFcn'',gca,[],guidata(gcbo))');

sp3 = subplot(2,2,3);
if length(yaxis) > 0 && length(zaxis) > 0
   imagesc([yaxis(1) yaxis(end)],[zaxis(1) zaxis(end)],sliceZY,clims)
else
imagesc(sliceZY, clims);
end
if (handles.axis_equal == 1)
    axis image;
else
    axis normal;
end;
title('Slice ZY');
if exist('axlabel')
    ylabel(['Z',axlabel]);xlabel(['Y',axlabel])
else
ylabel('Z');xlabel('Y');
end

if length(yaxis) > 0 && length(zaxis)
    line([-size(handles.volume,2) size(handles.volume,2)], [handles.pointer3dt(3).*deltazaxis - zaxis(end) handles.pointer3dt(3).*deltazaxis - zaxis(end)]);
    line([handles.pointer3dt(2).*deltayaxis - yaxis(end) handles.pointer3dt(2).*deltayaxis - yaxis(end)], [-size(handles.volume,3).*deltazaxis size(handles.volume,3).*deltazaxis]);
else
line([0 size(handles.volume,2)], [handles.pointer3dt(3) handles.pointer3dt(3)]);
line([handles.pointer3dt(2) handles.pointer3dt(2)], [0 size(handles.volume,3)]);
end
%set(allchild(gca),'ButtonDownFcn',@Subplot3_ButtonDownFcn);
set(allchild(gca),'ButtonDownFcn','SliceBrowser(''Subplot3_ButtonDownFcn'',gca,[],guidata(gcbo))');

function pointer3d_out = clipointer3d(pointer3d_in,vol_size)
pointer3d_out = pointer3d_in;
for p_id=1:4
    if (pointer3d_in(p_id) > vol_size(p_id))
        pointer3d_out(p_id) = vol_size(p_id);
    end;
    if (pointer3d_in(p_id) < 1)
        pointer3d_out(p_id) = 1;
    end;
end;


text:
function [t,tm,b] = get_max_t(Zpop,sterr,tt)
% :Usage:
% ::
%
%     [t,tm,b] = get_max_t(Zpop,sterr,tt)
%
% :Outputs:
%
%   **t:**
%        t-value timeseries
%
%   **tm:**
%        max t-value (abs)
%
%   **b:**
%        time (index) of max t-value
%

mu = mean(Zpop(1:tt));                              % population mean
t = (Zpop - mu) ./ sterr;
[tm, b] = max(abs(t(tt+1:end)));              % Calculate maximum absolute t-value
b = b+tt;                                     % max t time
tm = t(b);                                      % put sign back in

return


text:
function y=mvaveragec(x,N)

%MVAVERAGEC   signal smoothing through the moving average method.
%   Y = MVAVERAGEC(X,N) Quickly smooths the signal X via averaging each
%   sampling with the previous and afterwards N samples.
%
%   Example:
%      t = 2*pi*linspace(-1,1);
%      xn = cos(x) + 0.25 - 0.5*rand(size(x));
%      xs = mvaveragec(xn,4);
%      plot(t,xn,t,xs), legend('noisy','smooth'), axis tight
%
%   See also MVAVERAGE2 and MVAVERAGEC

%   by Yi Cao
%   Cranfield University
%   11-09-2007
%   y.cao@cranfield.ac.uk
%

y0=filter(ones(1,2*N+1)/(2*N+1),1,x);
y=x;
y(N+1:end-N)=y0(2*N+1:end);


text:
 function [kb, alpha, kb_m] = kaiser_bessel(x, J, alpha, kb_m, K_N)
%function [kb, alpha, kb_m] = kaiser_bessel(x, J, alpha, kb_m)
%function [kb, alpha, kb_m] = kaiser_bessel(x, J, 'best', 0, K_N)
%
% generalized Kaiser-Bessel function for x in support [-J/2,J/2]
% shape parameter "alpha" (default 2.34 J)
% order parameter "kb_m" (default 0)
% see (A1) in lewitt:90:mdi, JOSA-A, Oct. 1990
% in
%	x	[M,1]	arguments
% out
%	kb	[M,1]	KB function values, if x is numbers
%			or string for kernel(k,J), if x is 'string'
%			or inline function, if x is 'inline'
%	alpha
%	kb_m
%
% Copyright 2001-3-30, Jeff Fessler, The University of Michigan

% Modification 2002-10-29 by Samuel Matej
% - for Negative & NonInteger kb_m the besseli() function has
%	singular behavior at the boundaries - KB values shooting-up/down
%	(worse for small alpha) leading to unacceptable interpolators
% - for real arguments and higher/reasonable values of alpha the
%	besseli() gives similar values for positive and negative kb_m
%	except close to boundaries - tested for kb_m=-2.35:0.05:2.35
%	(besseli() gives exactly same values for integer +- kb_m)
%	=> besseli(kb_m,...) approximated by besseli(abs(kb_m),...), which
%	behaves well at the boundaries
% WARNING: it is not clear how correct the FT formula (JOSA) is
%	for this approximation (for NonInteger Negative kb_m)
% NOTE: Even for the original KB formula, the JOSA FT formula
%	is derived only for m > -1 !

% if no arguments, make example plots
if nargin < 2
	help(mfilename)
	J = 8; alpha = 2.34 * J;
	x = linspace(-(J+1)/2, (J+1)/2, 1001)';
%	x = linspace(J/2-1/4, J/2+1/4, 1001)';

	mlist = [-4 0 2 7];
	leg = {};
	for ii=1:length(mlist)
		kb_m = mlist(ii);
		yy(:,ii) = kaiser_bessel(x, J, alpha, kb_m);
		func = kaiser_bessel('inline', 0, alpha, kb_m);
		yf = func(x, J);
		if any(yf ~= yy(:,ii)),
		[yf yy(:,ii)]
		error 'bug', end
		leg{ii} = sprintf('m=%d', kb_m);
	end
	yb = kaiser_bessel(x, J, 'best', [], 2);
	plot(	x, yy(:,1), 'c-', x, yy(:,2), 'y-', ...
		x, yy(:,3), 'm-', x, yy(:,4), 'g-', x, yb, 'r--')
	leg{end+1} = 'best';
	axis tight, legend(leg)
%	axisy(0, 0.01)	% to see endpoints
	xlabel \kappa, ylabel F(\kappa)
	title(sprintf('KB functions, J=%g \\alpha=%g', J, alpha) )
return
end

if ~isvar('J'), J = 6; end
if ~isvar('alpha') | isempty('alpha'), alpha = 2.34 * J; end
if ~isvar('kb_m') | isempty('kb_m'), kb_m = 0; end

if ischar(alpha)
	[alpha kb_m] = kaiser_bessel_params(alpha, J, K_N);
end

if ischar(x)
	if ischar(alpha)
		if ~isvar('K_N'), error 'K_N required', end
		kb = 'kaiser_bessel(k, J, ''%s'', [], %g)';
		kb = sprintf(kb, alpha, K_N);
	else
		kernel_string = 'kaiser_bessel(k, J, %g, %g)';
		kb = sprintf(kernel_string, alpha, kb_m);
	end
	if streq(x, 'inline')
		kb = inline(kb, 'k', 'J');
	elseif ~streq(x, 'string')
		error '1st argument must be "inline" or "string"'
	end
return
end

%
% Warn about use of modified formula for negative kb_m
%
if (kb_m < 0) & ((abs(round(kb_m)-kb_m)) > eps)
	persistent warned
	if isempty(warned)	% only print this reminder the first time
		printf('\nWarning: Negative NonInt kb_m=%g in kaiser_bessel()', kb_m)
		printf('	- using modified definition of KB function\n')
		warned = 1;
	end
end

kb_m_bi = abs(kb_m);		% modified "kb_m" as described above
ii = abs(x) < J/2;
f = sqrt(1 - (x(ii)/(J/2)).^2);
denom = besseli(kb_m_bi,alpha);
if ~denom
	printf('m=%g alpha=%g', kb_m, alpha)
end
kb = zeros(size(x));
kb(ii) = f.^kb_m .* besseli(kb_m_bi, alpha*f) / denom;
kb = reale(kb);

%
% optimized shape and order parameters
%
function [alpha, kb_m] = kaiser_bessel_params(alpha, J, K_N)
if streq(alpha, 'best')
	if K_N ~= 2
		warning 'kaiser_bessel optimized only for K/N=2'
		printf('using good defaults: m=0 and alpha = 2.34*J')
		kb_m = 0;
		alpha = 2.34 * J;
	else
		kb_m = 0;	% hardwired, because it was nearly the best!
		try
			s = 'private/kaiser,m=0';
			s = load(s);
			ii = find(J == s.Jlist);
			if isempty(ii)
				ii = imin(abs(J - s.Jlist));
				warning(sprintf('J=%d not found, using %d', ...
					J, s.Jlist(ii)))
			end
			alpha = J * s.abest.zn(ii);
		catch
			warning(['could not open file "' s '" so using default alpha = 2.34 J which should be fine.'])
			alpha = 2.34 * J;
		end
	end
else
	error 'unknown alpha mode'
end

author: thomaskuestner


text:
% drawcloud ... draw point cloud
%
% [fig] = drawcloud(X,fig,{color})
%
% X ... 3xn matrix containing the points
%       if X is 4xn only the first 3 rows are used
% fig . figure handle
% color color of plotting; defaults to blue
%
% fig . return the figure handle
%
% $Id: drawcloud.m,v 2.0 2003/06/19 12:07:02 svoboda Exp $

function [fig] = drawcloud(X,fig,color)
if nargin < 3
  color = 'b';	% default color
end

figure(fig), hold on
plot3(X(1,:),X(2,:),X(3,:),[color,'o'])

view([1,1,1]);
axis('equal');
grid on
rotate3d on
fig=fig;
return

author: strawlab

repo: MultiCamSelfCal