text:
function [ t_stop, y_stop ] = p00_stop ( test, neqn )

%*****************************************************************************80
%
%% P00_STOP returns the end data for any test problem.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    23 February 2013
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer TEST, the test problem index.
%
%    Input, integer NEQN, the number of equations.
%
%    Output, real T_STOP, the end time.
%
%    Output, real Y_STOP(NEQN), the end value of the dependent
%    variables.
%
  y_stop = zeros ( neqn, 1 );

  if ( test == 1 )
    [ t_stop, y_stop ] = p01_stop ( neqn );
  elseif ( test == 2 )
    [ t_stop, y_stop ] = p02_stop ( neqn );
  elseif ( test == 3 )
    [ t_stop, y_stop ] = p03_stop ( neqn );
  elseif ( test == 4 )
    [ t_stop, y_stop ] = p04_stop ( neqn );
  elseif ( test == 5 )
    [ t_stop, y_stop ] = p05_stop ( neqn );
  elseif ( test == 6 )
    [ t_stop, y_stop ] = p06_stop ( neqn );
  elseif ( test == 7 )
    [ t_stop, y_stop ] = p07_stop ( neqn );
  elseif ( test == 8 )
    [ t_stop, y_stop ] = p08_stop ( neqn );
  elseif ( test == 9 )
    [ t_stop, y_stop ] = p09_stop ( neqn );
  elseif ( test == 10 )
    [ t_stop, y_stop ] = p10_stop ( neqn );
  elseif ( test == 11 )
    [ t_stop, y_stop ] = p11_stop ( neqn );
  elseif ( test == 12 )
    [ t_stop, y_stop ] = p12_stop ( neqn );
  elseif ( test == 13 )
    [ t_stop, y_stop ] = p13_stop ( neqn );
  elseif ( test == 14 )
    [ t_stop, y_stop ] = p14_stop ( neqn );
  elseif ( test == 15 )
    [ t_stop, y_stop ] = p15_stop ( neqn );
  elseif ( test == 16 )
    [ t_stop, y_stop ] = p16_stop ( neqn );
  elseif ( test == 17 )
    [ t_stop, y_stop ] = p17_stop ( neqn );
  elseif ( test == 18 )
    [ t_stop, y_stop ] = p18_stop ( neqn );
  elseif ( test == 19 )
    [ t_stop, y_stop ] = p19_stop ( neqn );
  elseif ( test == 20 )
    [ t_stop, y_stop ] = p20_stop ( neqn );
  elseif ( test == 21 )
    [ t_stop, y_stop ] = p21_stop ( neqn );
  elseif ( test == 22 )
    [ t_stop, y_stop ] = p22_stop ( neqn );
  elseif ( test == 23 )
    [ t_stop, y_stop ] = p23_stop ( neqn );
  elseif ( test == 24 )
    [ t_stop, y_stop ] = p24_stop ( neqn );
  elseif ( test == 25 )
    [ t_stop, y_stop ] = p25_stop ( neqn );
  elseif ( test == 26 )
    [ t_stop, y_stop ] = p26_stop ( neqn );
  elseif ( test == 27 )
    [ t_stop, y_stop ] = p27_stop ( neqn );
  elseif ( test == 28 )
    [ t_stop, y_stop ] = p28_stop ( neqn );
  elseif ( test == 29 )
    [ t_stop, y_stop ] = p29_stop ( neqn );
  elseif ( test == 30 )
    [ t_stop, y_stop ] = p30_stop ( neqn );
  elseif ( test == 31 )
    [ t_stop, y_stop ] = p31_stop ( neqn );
  elseif ( test == 32 )
    [ t_stop, y_stop ] = p32_stop ( neqn );
  elseif ( test == 33 )
    [ t_stop, y_stop ] = p33_stop ( neqn );
  elseif ( test == 34 )
    [ t_stop, y_stop ] = p34_stop ( neqn );
  elseif ( test == 35 )
    [ t_stop, y_stop ] = p35_stop ( neqn );
  elseif ( test == 36 )
    [ t_stop, y_stop ] = p36_stop ( neqn );
  elseif ( test == 37 )
    [ t_stop, y_stop ] = p37_stop ( neqn );
  elseif ( test == 38 )
    [ t_stop, y_stop ] = p38_stop ( neqn );
  elseif ( test == 39 )
    [ t_stop, y_stop ] = p39_stop ( neqn );
  elseif ( test == 40 )
    [ t_stop, y_stop ] = p40_stop ( neqn );
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P00_STOP - Fatal error!\n' );
    fprintf ( 1, '  Illegal problem index TEST = %d\n', test );
    error ( 'P00_STOP - Fatal error!' );
  end

  return
end


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

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

% Last Modified by GUIDE v2.5 24-Nov-2008 12:22:09

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

% Choose default command line output for Blister_GUI
handles.output = hObject;
axes(handles.axes1);
imshow('sample 01.jpg');  % Show a default sample image

% Create video input objects
handles.vid_obj = videoinput('winvideo');
handles.obj_src= get(handles.vid_obj,'Source');

% Set device parameters
set(handles.obj_src,'Brightness',160);
set(handles.obj_src,'Saturation',160);
set(handles.vid_obj,'FramesPerTrigger',1);
triggerconfig(handles.vid_obj,'manual');
set(handles.vid_obj,'TriggerRepeat',inf);
start(handles.vid_obj);

handles.vidRes = get(handles.vid_obj, 'VideoResolution');
handles.nBands = get(handles.vid_obj, 'NumberOfBands');

axes(handles.axes2);
himage = image( zeros(handles.vidRes(2), handles.vidRes(1), handles.nBands) );
preview(handles.vid_obj,himage);
% Update handles structure
guidata(hObject, handles);

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

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

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

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

% Update the status messages
set(handles.text9,'string','Status:  Taking a snapshot...');

% Send a manual trigger to the camera
trigger(handles.vid_obj);

% Get data from the buffer and display
imdata = getdata(handles.vid_obj);
axes(handles.axes1);
imshow(imdata);
axes(handles.axes1);
set(handles.text9,'string','Status:  Searching for broken pills...');
drawnow

% Call our code in Blister Read
[numPills, numBroken] = BlisterRead(imdata);

% Display information regarding broken and full pills
set(handles.text7,'string',num2str(numPills));
set(handles.text8,'string',num2str(numBroken));
set(handles.text9,'string','Status:  Estimated broken pills.');
hold off;

% --- Executes when user attempts to close figure1.
function figure1_CloseRequestFcn(hObject, eventdata, handles)
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
stop(handles.vid_obj);
% Hint: delete(hObject) closes the figure
delete(hObject);

author: Sable

repo: mcbench-benchmarks


text:
function db_insert(tablename, s)

% DB_INSERT inserts a structure into a database table. Each field of
% the structure should correspond with one of the fields in the table.
%
% Use as
%   db_insert(tablename, s)

% Copyright (C) 2007, Robert Oostenveld
%
% $Log: db_insert.m,v $
% Revision 1.2  2007/12/03 12:35:24  roboos
% fixed bug when input is double([])
%
% Revision 1.1  2007/11/07 10:50:46  roboos
% created helper functions for easy access to a MySQL database table using a structure for representing the data
%

if isempty(s)
  % the input consists of an empty structure. so nothing to do
  return
end

f  = fieldnames(s);
v = {};
t = {};
l = length(f);
for i=1:l
  v{i} = getfield(s, f{i});
  t{i} = class(v{i});
end

for j=1:numel(s)
  % loop over all elements of a structure array

  cmd = ['INSERT INTO ' tablename ' ( '];

  % these are the fieldnames
  for i=1:l
    cmd = [cmd f{i} ', '];
  end
  cmd = cmd(1:(end-2));  % remove the last ', '
  cmd = [cmd ' ) VALUES ( '];

  % these are the corresponding values
  for i=1:l
    switch t{i}
      case 'char'
        if isempty(v{i})
          cmd = [cmd '""' ];
        else
          cmd = [cmd '"' v{i} '"' ];
        end

      case {'int8', 'int16', 'int32', 'uint8', 'uint16', 'uint32', 'double'}
        if isempty(v{i})
          cmd = [cmd '""' ];
        else
          cmd = [cmd num2str(v{i})];
        end

      otherwise
        ft_error('unsuported data type in structure');
    end
    cmd = [cmd ', '];
  end

  cmd = cmd(1:(end-2));  % remove the last ', '
  cmd = [cmd ' ) ;'];

  % execute the query
  mysql(cmd);

end %for numel(s)

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [M,linearInd]=subImage(M,nSub)

% function M=subImage(M,nSub)
% ------------------------------------------------------------------------
%
%This function refines the input image M by splitting its voxels in half during nSub times
%which contains all folders and sub-folders within the folder specified by
%the input pathName.
%
% Kevin Mattheus Moerman
% kevinmoerman@hotmail.com
% 18/04/2013
%------------------------------------------------------------------------

%%
%TO DO: ADD WARNINGS AND INPUT PARSING HERE
% nSub=round(nSub);

%%

if numel(nSub)==1
    nSub=nSub.*ones(1,3);
end

stepSize_I=1/nSub(1);
stepSize_J=1/nSub(2);
stepSize_K=1/nSub(3);

%Define image coordinates of new voxel centres within image coordinate
%system of original coordinate mesh
I_range=linspace((0.5+(stepSize_I/2)),(size(M,1)-(stepSize_I/2))+0.5,size(M,1)*nSub(1));
J_range=linspace((0.5+(stepSize_J/2)),(size(M,2)-(stepSize_J/2))+0.5,size(M,2)*nSub(2));
K_range=linspace((0.5+(stepSize_K/2)),(size(M,3)-(stepSize_K/2))+0.5,size(M,3)*nSub(3));

[In,Jn,Kn]=ndgrid(I_range,J_range,K_range); %Coordinate matrices

%Rounding the coordinates so they "snap" to the mother voxel indices they
%are found in.
In=round(In); Jn=round(Jn); Kn=round(Kn);

%Convert these subscript indiced to linear indices in the image
linearInd =sub2ind(size(M),In,Jn,Kn);

%Get intensity values
Mc=M(:); %Image as column
M=Mc(linearInd); %The refined image

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
% This function computes the entropy of a label set.

% Copyright (C) 2012 Quan Wang <wangq10@rpi.edu>,
% Signal Analysis and Machine Perception Laboratory,
% Department of Electrical, Computer, and Systems Engineering,
% Rensselaer Polytechnic Institute, Troy, NY 12180, USA
%
% You are free to use this software for academic purposes if you cite our paper:
% Q. Wang, Y. Ou, A.A. Julius, K.L. Boyer, M.J. Kim,
% Tracking tetrahymena pyriformis cells using decision trees,
% in: 2012 International Conference on Pattern Recognition, Tsukuba Science City, Japan.
%
% For commercial use, please contact the authors.

function label_ent=getEntropy(label)

% label: one label vector
% label_ent: resulting entropy

n1=sum(1-label);
n2=sum(label);
if n1+n2==0
    label_ent=0;
    return;
end
p1=n1/(n1+n2);
p2=n2/(n1+n2);
if p1==0
    label_ent=-p2*log(p2);
    return;
elseif p2==0
    label_ent=-p1*log(p1);
    return;
end
label_ent=-p1*log(p1)-p2*log(p2);


text:
function [newnode,newelem,newelem0]=surfboolean(node,elem,varargin)
%
% [newnode,newelem,newelem0]=surfboolean(node1,elem1,op2,node2,elem2,op3,node3,elem3,...)
%
% merge two or more triangular meshes and resolve intersecting elements
%
% author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
%
% input:
%      node: node coordinates, dimension (nn,3)
%      elem: tetrahedral element or triangle surface (ne,3)
%      op:  a string of a boolean operator, possible op values include
%           'union' or 'or': the outter surface of the union of the enclosed space
%           'inter' or 'and': the surface of the domain contained by both meshes
%           'diff' or '-': the surface of the domain in mesh 1 excluding that of
%                   mesh 2
%           'all' or 'xor' or '+': the output contains 4 subsurfaces, identified by the 4th
%                  column of newelem:
%                    1: mesh 1 outside of mesh 2
%                    2: mesh 2 outside of mesh 1
%                    3: mesh 1 inside of mesh 2
%                    4: mesh 2 inside of mesh 1
%                  you can use newelem(find(mod(newelem(:,4),2)==1),:) to
%                  get mesh 1 cut by mesh 2, or newelem(find(mod(newelem(:,4),2)==0),:)
%                  to get mesh 2 cut by mesh 1;
%           'first': combine 1 and 3 from the output of 'all'
%           'second': combine 2 and 4 from the output of 'all'
%           'self': test for self-intersections; only the first mesh is
%                   tested; other inputs are ignored.
%           'decouple': separate two shells and make sure there is no intersection;
%                   the input surfaces must be closed and ordered from outer to inner
%
% output:
%      newnode: the node coordinates after boolean operations, dimension (nn,3)
%      newelem: tetrahedral element or surfaces after boolean operations (nn,4) or (nhn,5)
%      newelem0: when the operator is 'self', return the intersecting
%               element list in terms of the input node list (experimental)
%
% example:
%
%   [node1,face1,elem1]=meshabox([0 0 0],[10 10 10],1,1);
%   [node2,face2,elem2]=meshabox([0 0 0]+5,[10 10 10]+5,1,1);
%   [newnode,newface]=surfboolean(node1,face1,'union',node2,face2);
%   plotmesh(newnode,newface);
%   figure;
%   [newnode,newface]=surfboolean(node1,face1,'diff',node2,face2);
%   plotmesh(newnode,newface,'x>5');
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

allinputs=varargin;
opt=struct;
if(length(allinputs)>0 && isstruct(allinputs{end}))
    opt=allinputs{end};
    allinputs{end}=[];
end
len=length(varargin);
newnode=node;
newelem=elem;
if(len>0 && mod(len,3)~=0)
   error('you must give operator, node and element in trilet forms');
end

exesuff=fallbackexeext(getexeext,'gtsset');

for i=1:3:len
   op=varargin{i};
   no=varargin{i+1};
   el=varargin{i+2};
   opstr=op;
   if(strcmp(op,'or'))   opstr='union'; end
   if(strcmp(op,'xor'))  opstr='all';   end
   if(strcmp(op,'and'))  opstr='inter'; end
   if(strcmp(op,'-'))    opstr='diff';  end
   if(strcmp(op,'self')) opstr='inter -s';  end
   if(strcmp(op,'first') || strcmp(op,'second') || strcmp(op,'+'))
       opstr='all';
   end

   deletemeshfile(mwpath('pre_surfbool*.gts'));
   deletemeshfile(mwpath('post_surfbool.off'));
   if(strcmp(opstr,'all'))
      deletemeshfile(mwpath('s1out2.off'));
      deletemeshfile(mwpath('s1in2.off'));
      deletemeshfile(mwpath('s2out1.off'));
      deletemeshfile(mwpath('s2in1.off'));
   end
   if(strcmp(op,'decouple'))
       if(exist('node1','var')==0)
          node1=node;
          elem1=elem;
          newnode(:,4)=1;
          newelem(:,4)=1;
       end
       opstr=['-q --shells 2'];
       saveoff(node1(:,1:3),elem1(:,1:3),mwpath('pre_decouple1.off'));
       if(isstruct(el))
           if(isfield(el,'MoreOptions'))
               opstr=[opstr el.MoreOptions];
           end
       else
           opstr=[opstr ' --decouple-inin 1'];
       end
       if(size(no,2)~=3)
           opstr=['-q --shells ' num2str(no)];
           cmd=sprintf('cd "%s" && "%s%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
               mwpath('pre_decouple1.off'),opstr);
       else
           saveoff(no(:,1:3),el(:,1:3),mwpath('pre_decouple2.off'));
           cmd=sprintf('cd "%s" && "%s%s" "%s" "%s" %s',mwpath,mcpath('meshfix'),exesuff,...
               mwpath('pre_decouple1.off'),mwpath('pre_decouple2.off'),opstr);
       end
   else
       savegts(newnode(:,1:3),newelem(:,1:3),mwpath('pre_surfbool1.gts'));
       savegts(no(:,1:3),el(:,1:3),mwpath('pre_surfbool2.gts'));
       cmd=sprintf('cd "%s" && "%s%s" %s "%s" "%s" -v > "%s"',mwpath,mcpath('gtsset'),exesuff,...
           opstr,mwpath('pre_surfbool1.gts'),mwpath('pre_surfbool2.gts'),mwpath('post_surfbool.off'));
   end
   [status outstr]=system(cmd);
   if(status~=0 && strcmp(op,'self')==0)
       error(sprintf('surface boolean command failed:\n%s\nERROR: %s\n',cmd,outstr));
   end
   if(status~=0 && strcmp(op,'self') && ~isempty(strfind(outstr,'(new_ear): assertion failed')))
       fprintf(1,'no self-intersection was found! (ignore the above error)\n');
       newnode=[];
       newelem=[];
       newelem0=[];
       return;
   end
   if(strcmp(opstr,'all'))
      % tag the 4 piceses of meshes, this tag do not propagate to the next boolean operation
      [nnode nelem]=readoff(mwpath('s1out2.off'));
      newelem=[nelem ones(size(nelem,1),1)];
      newnode=[nnode ones(size(nnode,1),1)];

      [nnode nelem]=readoff(mwpath('s1in2.off'));
      newelem=[newelem; nelem+size(newnode,1) 3*ones(size(nelem,1),1)];
      newnode=[newnode; nnode 3*ones(size(nnode,1),1)];

      [nnode nelem]=readoff(mwpath('s2out1.off'));
      newelem=[newelem; nelem+size(newnode,1) 2*ones(size(nelem,1),1)];
      newnode=[newnode; nnode 2*ones(size(nnode,1),1)];

      [nnode nelem]=readoff(mwpath('s2in1.off'));
      newelem=[newelem; nelem+size(newnode,1) 4*ones(size(nelem,1),1)];
      newnode=[newnode; nnode 4*ones(size(nnode,1),1)];

      if(strcmp(op,'first'))
          newelem=newelem(find(mod(newelem(:,4),2)==1),:);
          [newnode,nelem]=removeisolatednode(newnode,newelem(:,1:3));
          newelem=[nelem newelem(:,4)];
      elseif(strcmp(op,'second'))
          newelem=newelem(find(mod(newelem(:,4),2)==0),:);
          [newnode,nelem]=removeisolatednode(newnode,newelem(:,1:3));
          newelem=[nelem,newelem(:,4)];
      end
   elseif(strcmp(op,'decouple'))
      [node1,elem1]=readoff(mwpath('pre_decouple1_fixed.off'));
      newelem=[newelem;elem1+size(newnode,1) (i+1)*ones(size(elem1,1),1)];
      newnode=[newnode;node1 (i+1)*ones(size(node1,1),1)];
   else
      [newnode,newelem]=readoff(mwpath('post_surfbool.off'));
      if(strcmp(op,'self'))
          fprintf(1,'a total of %d self-intersecting elements were found\n',size(newelem,1));
          if(nargout>=3)
              [found,newelem0]=ismember(newnode,node,'rows');
              if(~all(found))
                  error('self intersecting elements contain new nodes');
              end
              newelem0=newelem0(newelem);
          end
          return;
      end
   end
end

author: vigente

repo: gerardus

save_path: github-repos/MATLAB/vigente-gerardus


text:
function ber0_values_test ( )

%*****************************************************************************80
%
%% BER0_VALUES_TEST demonstrates the use of BER0_VALUES.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 February 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BER0_VALUES_TEST:\n' );
  fprintf ( 1, '  BER0_VALUES stores values of \n' );
  fprintf ( 1, '  the Kelvin function BER of order 0.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '      X           FX\n' );
  fprintf ( 1, '\n' );

  n_data = 0;

  while ( 1 )

    [ n_data, x, fx ] = ber0_values ( n_data );

    if ( n_data == 0 )
      break
    end

    fprintf ( 1, '  %12f  %24.16f\n', x, fx );

  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
%   [all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
%   logisitc regression classifiers and returns each of these classifiers
%   in a matrix all_theta, where the i-th row of all_theta corresponds
%   to the classifier for label i

% Some useful variables
m = size(X, 1);
n = size(X, 2);

% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);

% Add ones to the X data matrix
X = [ones(m, 1) X];

% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
%               logistic regression classifiers with regularization
%               parameter lambda.
%
% Hint: theta(:) will return a column vector.
%
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
%       whether the ground truth is true/false for this class.
%
% Note: For this assignment, we recommend using fmincg to optimize the cost
%       function. It is okay to use a for-loop (for c = 1:num_labels) to
%       loop over the different classes.
%
%       fmincg works similarly to fminunc, but is more efficient when we
%       are dealing with large number of parameters.
%
% Example Code for fmincg:
%
%     % Set Initial theta
%     initial_theta = zeros(n + 1, 1);
%
%     % Set options for fminunc
%     options = optimset('GradObj', 'on', 'MaxIter', 50);
%
%     % Run fmincg to obtain the optimal theta
%     % This function will return theta and the cost
%     [theta] = ...
%         fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
%                 initial_theta, options);
%

for c = 1:num_labels
    options = optimset('GradObj', 'on', 'MaxIter', 50);

    all_theta(c, :) = fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
            zeros(n + 1, 1), options);

% =========================================================================

end

author: SaveTheRbtz

repo: ml-class


text:
%ARROWEDLINE  Draws an arrow segment pointing from the first point to the second one
%
%     img = cv.arrowedLine(img, pt1, pt2)
%     [...] = cv.arrowedLine(..., 'OptionName', optionValue, ...)
%
% ## Input
% * __img__ Image where the arrow is drawn.
% * __pt1__ The point `[x,y]` the arrow starts from.
% * __pt2__ The point `[x,y]` the arrow points to.
%
% ## Output
% * __img__ Output image.
%
% ## Options
% * __Color__ Line color. default is a black color
% * __Thickness__ Line thickness. default 1.
% * __LineType__ Type of the line. One of:
%   * __4__ 4-connected line
%   * __8__ 8-connected line (default)
%   * __AA__ anti-aliased line
% * __Shift__ Number of fractional bits in the point coordinates. default 0.
% * __TipLength__ The length of the arrow tip in relation to the arrow length.
%   default 0.1
%
% The function cv.arrowedLine draws an arrow between `pt1` and `pt2` points
% in the image.
%
% See also: cv.line
%

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%% Brightness and contrast adjustments
%
% In this demo we show how to perform the operation
% $g(i,j) = \alpha \cdot f(i,j) + \beta$.
%
% Sources:
%
% * <https://docs.opencv.org/3.3.0/d3/dc1/tutorial_basic_linear_transform.html>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/ImgProc/BasicLinearTransforms.cpp>
% * <https://github.com/opencv/opencv/blob/3.2.0/samples/cpp/tutorial_code/HighGUI/BasicLinearTransformsTrackbar.cpp>
% * <https://github.com/opencv/opencv/blob/3.3.0/samples/cpp/tutorial_code/ImgProc/changing_contrast_brightness_image/changing_contrast_brightness_image.cpp>
%

%% Theory
%
% A general image processing operator is a function that takes one or more
% input images and produces an output image. Image transforms can be seen as:
%
% * Point operators (pixel transforms)
% * Neighborhood (area-based) operators
%
% In pixel transforms, each output pixel's value depends on only the
% corresponding input pixel value (plus, potentially, some globally collected
% information or parameters). Examples of such operators include
% *brightness and contrast adjustments* as well as color correction and
% transformations.
%
% For brightness and contrast adjustments, two commonly used point processes
% are *multiplication* and *addition* with a constant:
%
% $$g(x) = \alpha f(x) + \beta$$
%
% The parameters $\alpha > 0$ and $\beta$ are often called the *gain* and
% *bias* parameters; sometimes these parameters are said to control *contrast*
% and *brightness* respectively. You can think of $f(x)$ as the source image
% pixels and $g(x)$ as the output image pixels. Then, more conveniently we can
% write the expression as:
%
% $$g(i,j) = \alpha \cdot f(i,j) + \beta$$
%
% where $i$ and $j$ indicates that the pixel is located in the *i-th* row and
% *j-th* column.
%
% In the implementation below, instead of using a for-loop to access each
% pixel, we simply use the function |cv.convertTo| which effectively performs
% *new_image = saturate(a*image + beta)*. This is more optimized than
% accessing each pixel and works a lot faster. Also notice that the operation
% can give values out of range or not integers (if $\alpha$ is float), in
% which case |cv.convertTo| makes sure the values are valid.
%

%% Example
%
% We will put into practice what we have learned to correct an underexposed
% image by adjusting the brightness and the contrast of the image. We will
% also see another technique to correct the brightness of an image called
% gamma correction.
%
%% Brightness and contrast adjustments
%
% Increasing/decreasing the $\beta$ value will add/subtract a constant value
% to every pixel. Pixel values outside of the |[0 ; 255]| range will be
% saturated (i.e. a pixel value higher/lesser than 255/0 will be clamp to
% 255/0).
%
% <<https://docs.opencv.org/3.3.0/Basic_Linear_Transform_Tutorial_hist_beta.png>>
%
% _In light gray, histogram of the original image, in dark gray when
% |brightness = 80| in Gimp_
%
% The histogram represents for each color level the number of pixels with that
% color level. A dark image will have many pixels with low color value and
% thus the histogram will present a peak in his left part. When adding a
% constant bias, the histogram is shifted to the right as we have added a
% constant bias to all the pixels.
%
% The $\alpha$ parameter will modify how the levels spread. If $\alpha < 1$,
% the color levels will be compressed and the result will be an image with
% less contrast.
%
% <<https://docs.opencv.org/3.3.0/Basic_Linear_Transform_Tutorial_hist_alpha.png>>
%
% _In light gray, histogram of the original image, in dark gray when
% |contrast < 0| in Gimp_
%
% Note that these histograms have been obtained using the Brightness-Contrast
% tool in the Gimp software. The brightness tool should be identical to the
% $\beta$ bias parameters but the contrast tool seems to differ to the
% $\alpha$ gain where the output range seems to be centered with Gimp (as you
% can notice in the previous histogram).
%
% It can occur that playing with the $\beta$ bias will improve the brightness
% but in the same time the image will appear with a slight veil as the
% contrast will be reduced. The $\alpha$ gain can be used to diminue this
% effect but due to the saturation, we will lose some details in the original
% bright regions.
%
%% Gamma correction
%
% <https://en.wikipedia.org/wiki/Gamma_correction Gamma correction> can be
% used to correct the brightness of an image by using a non linear
% transformation between the input values and the mapped output values:
%
% $$O = \left( \frac{I}{255} \right)^{\gamma} \times 255$$
%
% As this relation is non linear, the effect will not be the same for all the
% pixels and will depend to their original value.
%
% <<https://docs.opencv.org/3.3.0/Basic_Linear_Transform_Tutorial_gamma.png>>
%
% When $\gamma < 1$, the original dark regions will be brighter and the
% histogram will be shifted to the right whereas it will be the opposite with
% $ gamma > 1$.
%
% For the gamma correction, a look-up table can be used to improve the
% performance of the computation as only 256 values needs to be calculated
% once.
%
% Let's an example of how to correct an underexposed image.
%
% The following image has been corrected with: $\alpha = 1.3$ and $\beta = 40$.
%
% <<https://docs.opencv.org/3.3.0/Basic_Linear_Transform_Tutorial_linear_transform_correction.jpg>>
%
% The overall brightness has been improved but you can notice that the clouds
% are now greatly saturated due to the numerical saturation of the
% implementation used
% (<https://en.wikipedia.org/wiki/Clipping_(photography) highlight clipping>
% in photography).
%
% The following image has been corrected with: $\gamma = 0.4$.
%
% <<https://docs.opencv.org/3.3.0/Basic_Linear_Transform_Tutorial_gamma_correction.jpg>>
%
% The gamma correction should tend to add less saturation effect as the
% mapping is non linear and there is no numerical saturation possible as in
% the previous method.
%
% <<https://docs.opencv.org/3.3.0/Basic_Linear_Transform_Tutorial_histogram_compare.png>>
%
% * Left: histogram after alpha, beta correction
% * Center: histogram of the original image
% * Right: histogram after the gamma correction
%
% The previous figure compares the histograms for the three images (the
% y-ranges are not the same between the three histograms). You can notice that
% most of the pixel values are in the lower part of the histogram for the
% original image. After $\alpha$, $\beta$ correction, we can observe a big
% peak at 255 due to the saturation as well as a shift in the right. After
% gamma correction, the histogram is shifted to the right but the pixels in
% the dark regions are more shifted (see the gamma curves figure) than those
% in the bright regions.
%
% In this tutorial, you have seen two simple methods to adjust the contrast
% and the brightness of an image. They are basic techniques and are not
% intended to be used as a replacement of a raster graphics editor!
%
%% Additional resources
%
% * <https://learnopengl.com/#!Advanced-Lighting/Gamma-Correction
%   Gamma correction in graphics rendering>
% * <http://web.archive.org/web/20170106081601/http://www.graphics.cornell.edu/~westin/gamma/gamma.html
%   Gamma correction and images displayed on CRT monitors>
% * <http://www.cambridgeincolour.com/tutorials/digital-exposure-techniques.htm
%   Digital exposure techniques>
%

%% Code

function varargout = linear_transform_demo_gui(im)
    % load source image
    if nargin < 1
        im = fullfile(mexopencv.root(),'test','lena.jpg');
        img = cv.imread(im);
    elseif ischar(im)
        img = cv.imread(im, 'Color',true);
    else
        img = im;
    end

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

function onLinearTransform(~,~,h)
    %ONLINEARTRANSFORM  Event handler for UI controls

    % retrieve current values from UI controls
    a = get(h.slid(1), 'Value');
    b = round(get(h.slid(2), 'Value'));
    set(h.txt(1), 'String',sprintf('Contrast: %.2f',a));
    set(h.txt(2), 'String',sprintf('Brightness: %2d',b));

    % linear transformation
    out = cv.convertTo(h.src, 'Alpha',a, 'Beta',b);

    % show result
    out = cv.putText(out, 'Brightness/Contrast Adjustment', [10 20], ...
        'FontScale',0.5, 'Color',[0 255 0], 'LineType','AA');
    set(h.img(1), 'CData',out);
    drawnow;
end

function onGammaCorrection(~,~,h)
    %ONGAMMACORRECTION  Event handler for UI controls

    % retrieve current values from UI controls
    g = get(h.slid(3), 'Value');
    set(h.txt(3), 'String',sprintf('Gamma: %.2f',g));

    % gamma correction
    lookUpTable = uint8((((0:255)/255) .^ g) * 255);
    if true
        out = cv.LUT(h.src, lookUpTable);
    elseif mexopencv.require('images')
        out = intlut(h.src, lookUpTable);
    else
        out = lookUpTable(double(h.src) + 1);
    end

    % show result
    out = cv.putText(out, 'Gamma Correction', [10 20], ...
        'FontScale',0.5, 'Color',[0 255 0], 'LineType','AA');
    set(h.img(2), 'CData',out);
    drawnow;
end

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

    % parameters
    a = 1.0;  % alpha gain (contrast)
    b = 0;    % beta bias (brightness)
    g = 1.0;  % gamma correction
    sz = size(img);
    sz(2) = max(sz(2), 300);  % minimum figure width

    % build the user interface (no resizing to keep it simple)
    h = struct();
    h.src = img;
    h.fig = figure('Name','Linear Transform Demo', ...
        'NumberTitle','off', 'Menubar','none', 'Resize','off', ...
        'Position',[200 200 sz(2)*2 sz(1)+80-1]);
    if ~mexopencv.isOctave()
        %HACK: not implemented in Octave
        movegui(h.fig, 'center');
    end
    h.ax(1) = axes('Parent',h.fig, 'Units','pixels', 'Position',[1 80 sz(2) sz(1)]);
    h.ax(2) = axes('Parent',h.fig, 'Units','pixels', 'Position',[sz(2)+1 80 sz(2) sz(1)]);
    if ~mexopencv.isOctave()
        h.img(1) = imshow(img, 'Parent',h.ax(1));
        h.img(2) = imshow(img, 'Parent',h.ax(2));
    else
        %HACK: https://savannah.gnu.org/bugs/index.php?45473
        axes(h.ax(1)); h.img(1) = imshow(img);
        axes(h.ax(2)); h.img(2) = imshow(img);
    end
    h.txt(1) = uicontrol('Parent',h.fig, 'Style','text', 'FontSize',11, ...
        'Position',[5 5 130 20], 'String',sprintf('Contrast: %.2f',a));
    h.txt(2) = uicontrol('Parent',h.fig, 'Style','text', 'FontSize',11, ...
        'Position',[5 30 130 20], 'String',sprintf('Brightness: %2d',b));
    h.txt(3) = uicontrol('Parent',h.fig, 'Style','text', 'FontSize',11, ...
        'Position',[5 55 130 20], 'String',sprintf('Gamma: %.2f',g));
    h.slid(1) = uicontrol('Parent',h.fig, 'Style','slider', 'Value',a, ...
        'Min',0.1, 'Max',3, 'SliderStep',[0.01 0.2]./(3-0.1), ...
        'Position',[135 5 sz(2)-135-5 20]);
    h.slid(2) = uicontrol('Parent',h.fig, 'Style','slider', 'Value',b, ...
        'Min',-100, 'Max',100, 'SliderStep',[1 20]./(100+100), ...
        'Position',[135 30 sz(2)-135-5 20]);
    h.slid(3) = uicontrol('Parent',h.fig, 'Style','slider', 'Value',g, ...
        'Min',0.1, 'Max',3, 'SliderStep',[0.01 0.2]./(3-0.1), ...
        'Position',[135 55 sz(2)-135-5 20]);

    % hook event handlers, and trigger default start
    opts = {'Interruptible','off', 'BusyAction','cancel'};
    set(h.slid(1:2), 'Callback',{@onLinearTransform,h}, opts{:});
    set(h.slid(3), 'Callback',{@onGammaCorrection,h}, opts{:});
    onLinearTransform([],[],h);
    onGammaCorrection([],[],h);
end


text:
function y = squish(x)
% y = squish(x);    SQUISH "x" to remove all singleton dimensions.
%   Since singleton dimensions can confuse many operations, this function removes
%   absolutely ALL singleton dimensions in "x". SQUEEZE is similar, however
%   it will not operate on 2D arrays, of which row vectors are included. Thus the
%   expected result may not always occur with SQUEEZE.
%   NOTE! This function will convert all row vectors to column vectors!
%   Example:        [1;2;3;4;5] = squish(shiftdim([1:5]',-5))

%   To see the differences between SQUISH and SQUEEZE compare the results
%   of the following for any positive or negative n:
%       size(squeeze(shiftdim([1:3]',n)))
%       size(squish(shiftdim([1:3]',n)))
%   created 08/16/2006 by Mirko Hrovat with Matlab ver.7.2
%   modified 02/25/2010 by MIH as per Jan Simon's suggestion (slight speed improvement)
%   contact: mhrovat@email.com

dims = size(x);
y = reshape(x,[dims(dims~=1),1,1]);    % the extra 1's help in case x is a scalar or vector.
%**************** END ***************


text:
function sfc_pres_temp_rd ( )

%*****************************************************************************80
%
%% SFC_PRES_TEMP_RD reads pressure and temperature arrays from a NETCDF file.
%
%  Discussion:
%
%    We are reading 2D data, 6 x 12 lat-lon grids for pressure and temperature.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    03 December 2009
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Russ Rew, Glenn Davis, Steve Emmerson, Harvey Davies, Ed Hartne,
%    The NETCDF User's Guide,
%    Unidata Program Center, March 2009.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SFC_PRES_TEMP_RD:\n' );
  fprintf ( 1, '  MATLAB version.\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Read 2D arrays of pressure and temperature\n' );
  fprintf ( 1, '  from a NETCDF file.\n' );
%
%  Open the file.
%
  nc_filename = 'sfc_pres_temp.nc';
  ncid = netcdf.open ( nc_filename, 'NOWRITE' );
%
%  There are a number of inquiry functions in NETCDF which can be
%  used to learn about an unknown NETCDF file.  INQ tells how many
%  NETCDF variables, dimensions, and global attributes are in the
%  file; also the dimension id of the unlimited dimension, if there
%  is one.
%
  [ ndims_in, nvars_in, ngatts_in, unlimdimid_in ] = netcdf.inq ( ncid );
%
%  In this case, we know that there are:
%    2 NETCDF dimensions,
%    4 NETCDF variables,
%    0 global attributes,
%    the unlimited dimension id is -1 (because there is no unlimited dimension).
%
  if ( ndims_in ~= 2 | ...
       nvars_in ~= 4 | ...
       ngatts_in ~= 0 | ...
       unlimdimid_in ~= -1 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SFC_PRES_TEMP_RD - Fatal error!\n' );
    fprintf ( 1, '  File has wrong properties.\n' );
    return
  end

  fprintf ( 1, '\n' );
  fprintf ( 1, '  The file properties are correct.\n' );
%
%  Get the varids of the latitude and longitude coordinate variables.
%
  lat_name = 'latitude';
  lat_varid = netcdf.inqVarid ( ncid, lat_name );

  lon_name = 'longitude';
  lon_varid = netcdf.inqVarid ( ncid, lon_name );
%
%  Read the latitude and longitude data.
%
  lats = netcdf.getVar ( ncid, lat_varid );
  lons = netcdf.getVar ( ncid, lon_varid );
%
%  Check to make sure we got what we expected.
%
  lat_num = 6;
  lat_start = 25.0;

  for lat = 1 : lat_num
    if (lats(lat) ~= lat_start + ( lat - 1 ) * 5.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SFC_PRES_TEMP_RD - Fatal error!\n' );
      fprintf ( 1, '  LAT data is wrong.\n' );
    end
  end
  fprintf ( 1, '  The LAT data is correct.\n' );

  lon_num = 12;
  lon_start = -125.0;

  for lon = 1 : lon_num
    if (lons(lon) ~= lon_start + ( lon - 1 ) * 5.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'SFC_PRES_TEMP_RD - Fatal error!\n' );
      fprintf ( 1, '  LON data is wrong.\n' );
    end
  end
  fprintf ( 1, '  The LON data is correct.\n' );
%
%  Get the varids of the pressure and temperature variables.
%
  press_name = 'pressure';
  press_varid = netcdf.inqVarid ( ncid, press_name );
  temp_name = 'temperature';
  temp_varid = netcdf.inqVarid ( ncid, temp_name );
%
%  Read the pressure and temperature data.
%
  press = netcdf.getVar ( ncid, press_varid );
  temp = netcdf.getVar ( ncid, temp_varid );
%
%  Check the data. It should be the same as the data we wrote.
%
  sample_pressure = 900.0;
  sample_temp = 9.0;

  for lon = 1 : lon_num
    for lat = 1 : lat_num

      if ( press(lon,lat) ~= sample_pressure + ( lon - 1 ) * lat_num + ( lat - 1 ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SFC_PRES_TEMP_RD - Fatal error!\n' );
        fprintf ( 1, '  PRESS data is wrong.\n' );
      end

      if ( temp(lon,lat) ~= sample_temp + 0.25 * ( ( lon - 1 ) * lat_num + ( lat - 1 ) ) )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'SFC_PRES_TEMP_RD - Fatal error!\n' );
        fprintf ( 1, '  TEMP data is wrong.\n' );
      end

    end
  end

  fprintf ( 1, '  The PRESS data is correct.\n' );
  fprintf ( 1, '  The TEMP data is correct.\n' );
%
%  Each of the NETCDF variables has a "units" attribute.
%  Let's read them and check them.
%
  units = 'units';

  lat_units = netcdf.getAtt ( ncid, lat_varid, units );
  if ( lat_units ~= 'degrees north' )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SFC_PRES_TEMP_RD - Fatal error!\n' );
    fprintf ( 1, '  LAT units is wrong.\n' );
  end

  lon_units = netcdf.getAtt ( ncid, lon_varid, units );
  if ( lon_units ~= 'degrees east' )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SFC_PRES_TEMP_RD - Fatal error!\n' );
    fprintf ( 1, '  LON units is wrong.\n' );
  end

  press_units = netcdf.getAtt ( ncid, press_varid, units );
  if ( press_units ~= 'hPa' )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SFC_PRES_TEMP_RD - Fatal error!\n' );
    fprintf ( 1, '  PRESS units is wrong.\n' );
  end

  temp_units = netcdf.getAtt ( ncid, temp_varid, units );
  if ( temp_units ~= 'celsius' )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'SFC_PRES_TEMP_RD - Fatal error!\n' );
    fprintf ( 1, '  TEMP units is wrong.\n' );
  end

  fprintf ( 1, '  The units data is correct.\n' );
%
%  Close the file.
%
  netcdf.close ( ncid );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The data was read from the file "%s".\n', nc_filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SFC_PRES_TEMP_RD:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end


text:
function eps_trial = eps_RCUs_Alamouti_stable(R, Mr, L, np, nc, rho, spa,...
                RELSPREAD_MAX, N_MC_MIN, N_MC_MAX, N_REP_TRIALS, d_verbose)
%EPS_RCUS_ALAMOUTI_STABLE returns the error probability computed via RCUs
%in a system using Alamouti coding (Mt=2 tx antennas) with parameters:
%R (rate), Mr (number of rx antennas), L (number of blocks),
%np (number of pilots), nc (size of coherence block), rho (SNR linear).

	% Default values for parameters that are not passed
    if nargin < 12
        d_verbose = 2;      % set verbosity
    end

    if nargin < 11
        ncores = feature('numcores');   % number of cores on current machine
        nworkers = ncores;              % use ncores-1 in parallel
        poolobj = gcp('nocreate');      % if no pool do not create new one
        if isempty(poolobj)
            parpool(nworkers);          % initialise the pool
        end
        % choose N_REP_TRIALS a multiple of no. of workers for high
        % efficiency
        N_REP_TRIALS = 3 * poolobj.NumWorkers;
    end

    if nargin < 10
        N_MC_MAX = 2^19;    % maximum number of Monte Carlo iterations
    end

    if nargin < 9
        N_MC_MIN = 2^10;    % initial number of Monte Carlo iterations
    end

    if nargin < 8
        RELSPREAD_MAX = 0.50;   % max relative spread of outcomes
    end

    if nargin < 7
        spa = 0;
    end
    % Algorithm parameters

    N_MC = min([N_MC_MIN N_MC_MAX]);	% initial number of Monte Carlo trials

    % Algorithm
    n = nc * L;
    b_iter = 1;

    if d_verbose >= 1
        fprintf([' Estimating Pe with ' ...
            '%d trials...\n'],N_MC);
    end

    while b_iter == 1

        eps_trial = zeros(N_REP_TRIALS,1);

        parfor i_trial = 1:N_REP_TRIALS
            if spa == 0
                % Generate information density samples
                i_s = idsamples_Alamouti(Mr, L, np, nc, rho, N_MC);

                % Compute error probability
                eps_trial(i_trial) = eps_RCUs(i_s, n, R);
            else %saddlepoint approximation
                % Generate information density samples with L = 1
                i_s = idsamples_Alamouti(Mr, 1, np, nc, rho, N_MC);

                % Compute error probability using saddlepoint approximation
                eps_trial(i_trial) = eps_RCUs_SPA(i_s, n, L, R);
            end

        end

        % Check spread of data
        relspread_eps = relspread(eps_trial); % relative spread

        if d_verbose >= 2
            fprintf('  %.3e', eps_trial);
            fprintf('\n  Relative spread: %.3e\n', relspread_eps);
        end

        % Stopping conditions
        if relspread_eps < RELSPREAD_MAX % accuracy reached
            b_iter = 0;
            if d_verbose >= 1
                fprintf(' Estimation accuracy reached with %d trials.\n',N_MC);
            end
        else
            if 2 * N_MC <= N_MC_MAX
                N_MC = 2 * N_MC;
                if d_verbose >= 1
                    fprintf(['  Accuracy not reached: '...
                        'trying estimation of Pe with ' ...
                        '%d trials...\n'], N_MC);
                end
            else
                b_iter = 0;
                if d_verbose >= 1
                    warning('::Relative accuracy target not reached');
                end
            end
        end
    end

end % of function

author: yp-mit

repo: spectre

save_path: github-repos/MATLAB/yp-mit-spectre


text:
function test_failed = test_fwt2

disp('========= TEST FWT2 ============');
global LTFAT_TEST_TYPE;
tolerance = 1e-8;
if strcmpi(LTFAT_TEST_TYPE,'single')
   tolerance = 1e-4;
end

test_failed = 0;

dims = { [20,30], [150,151],[226,253], };
flags = {'standard','tensor'};
filt = {{'mband1',2},{'db10',4},{'sym8',4},{'spline4:4',4},};

for ii=1:numel(dims)
   f = tester_rand(dims{ii});
   for jj=1:numel(flags)
      for ff=1:numel(filt)
         c = fwt2(f,filt{ff}{1},filt{ff}{2},flags{jj});
         fhat = ifwt2(c,filt{ff}{1},filt{ff}{2},dims{ii},flags{jj});
         err = norm(f-fhat,'fro');
         [test_failed,fail]=ltfatdiditfail(err,test_failed,tolerance);
         fprintf('J=%d, %5.5s, dim=[%3.d,%3.d], flag=%8.8s, err=%.4e %s\n',filt{ff}{2},filt{ff}{1},size(f,1),size(f,2),flags{jj},err,fail);
      end
   end
end

author: ltfat

repo: ltfat

save_path: github-repos/MATLAB/ltfat-ltfat


text:
% AXCOPY - Copy a Matlab figure axis and its graphic objects to a new pop-up window
%            using the left mouse button.
%
% Usage:  >> axcopy
%         >> axcopy(fig)
%         >> axcopy('noticks')
%         >> axcopy(fig, command)
%
% Notes:
%   1) Clicking the left mouse button on a Matlab figure axis copies the graphic objects
%      in the axis to a new (pop-up) figure window.
%   2) Option 'noticks' does not make x and y tickloabelmodes 'auto' in the pop-up.
%   2) The command option is an optional string that is evaluated in the new window
%      when it pops up. This allows the user to customize the pop-up display.
%   3) Deleting the pop-up window containing the copied axis leaves the selected axis
%      as the current graphic axis (gca).

% Authors: Tzyy-Ping Jung & Scott Makeig, SCCN/INC/UCSD, La Jolla, 2000

% Copyright (C) 2000 Tzyy-Ping Jung & Scott Makeig, SCCN/INC/UCSD, 3-16-00
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% requires copyaxes.m
% 4-2-00 added 'noticks' -sm
% 01-25-02 reformated help & license -ad
% 02-16-02 debuged exist('fig') -> exist('fig') == 1 -ad
% 02-16-02 added the ignore parent size option -ad
% 03-11-02 remove size option and added command option -ad

function axcopy(fig, command)

if (exist('fig') == 1) && strcmp(fig,'noticks')
   noticks = 1;
   if nargin> 1
     shift
   else
     fig = [];
   end
end
if ~(exist('fig') ==1) || isempty(fig) || fig == 0
   fig = gcf;
end

if ~strcmpi(get(fig, 'type'), 'axes')
    %hndl= findobj('parent',fig,'type','axes');
    %--
    childs = get(fig,'Children');
    hndl   = childs(find(strcmpi(get(childs,'Type'),'axes')));
else
    hndl=fig;
end
offidx=[];
if exist('command') ~= 1
    comstr = 'copyaxis';
else
   command_dbl = double(command);
   comstr = double(['copyaxis(''' char(command_dbl) ''')']);
end
for a=1:length(hndl)                    % make all axes visible
    %allobjs = findobj('parent',hndl(a));
    %--
    allobjs = get(hndl(a),'Children');
    for index = 1:length(allobjs)
        if isempty(get(allobjs(index), 'ButtonDownFcn'))
            set(allobjs(index), 'ButtonDownFcn', char(comstr));
        end
    end
end
if ~strcmpi(get(fig, 'type'), 'axes')
    figure(fig);
else
    figure(get(fig, 'parent'));
end
if exist('command') ~= 1
    set(hndl(a),'ButtonDownFcn','copyaxis');
else
    % set(hndl,'ButtonDownFcn',['copyaxis(''' command ''')']);
    comstr = double(['copyaxis(''' char(command_dbl) ''')']);
    set(hndl,'ButtonDownFcn',char(comstr));
end;
%set(hndl,'ButtonDownFcn','copyaxis');
%if ~exist('noticks')
  %axis on
  %set(gca,'xticklabelmode','auto')
  %set(gca,'yticklabelmode','auto')
%end

author: sccn

repo: eeglab

save_path: github-repos/MATLAB/sccn-eeglab


text:
function TCP = jeongLungTCPmodel(paramS,doseBinsV,volHistV)
% Usage: TCP = jeongLungTCPmodel(paramS);
% Lung TCP model
% Based on code by Jeho Jeong jeongj@mskcc.org
%------------------------------------------------------------------------------------
% INPUTS:
% paramS : Parameter dictionary with fields
%          -frxSize (Fraction size (Gy)),
%          -treatmentSchedule(Vector indicating treatment days)
% Example:
% paramS.frxSize.val = 2;
% paramS.treatmentSchedule.val = [ 1  2  3  4  5 ...
%                           8  9 10 11 12 ...
%                           15 16 17 18 19 ...
%                           22 23 24 25 26];
%
%----------------------------------------------------------------------------
% AI 8/30/18

%% Get fractionation
fx_in = paramS.frxSize.val;
nFrx = paramS.numFractions.val;

%% Get treatment days
schedule_in = paramS.treatmentSchedule.val;
if  isfield(paramS,'scheduleType')
    scheduleType = paramS.scheduleType;
else
    scheduleType = 'weekday';
end
%Check for numeric input
if ~isnumeric(schedule_in)
    scheduleV = str2num(schedule_in);
else
    scheduleV = schedule_in;
end
%Otherwise, assume function specified
if isempty(scheduleV)
    scheduleV = eval([schedule_in,'(',num2str(nFrx),',scheduleType)']);
    %scheduleV = getTreatmentSchedule(nFrx); %every weekday with weekend
                                             %breaks
end

%% Input variables for the analysis
alpha_p_ori=0.305;
a_over_b=2.8;
oer_i=1.7;
rho_t=10^6;
v_t_ref=3e4;
f_s=0.01;
t_c=2;
f_p_pro_in=0.5;
ht_loss=2;
k_m=0.3;
ht_lys=3;
oer_h=1.37;
F_p_cyc=[0.56;0.24;0.2];
Alpha_ratio_p_cyc=[2;3];

d_t=15;
clf_in=0.92;
gf_in=0.25;

beta_p_ori=alpha_p_ori/a_over_b;

%% EQD2 estimation for each cohort

EQD2=[];
n_pt=[];
v_t=3e4;
alpha_p=alpha_p_ori;
beta_p=beta_p_ori;
n_t=rho_t*v_t;                %No. cells
n_t_ref=rho_t*v_t_ref;
total_clono_cell=n_t*f_s;
delta_t=d_t/(60*24);          %dt in days
t_start=0;

IC=[];
GF=[];
TCP=[];
TD50=[];
BED=[];
Reox_time=[];
Reox_time2=[];
Treat_duration=[];
vec_leng=[];
comp_size(1)=0;       %Compartment size (P)
comp_size(2)=0;       %Compartment size (H)
comp_size(3)=0;       %Compartment size (I)
comp_size_ref(1)=0;
comp_size_ref(2)=0;
comp_size_ref(3)=0;
p_pre=[];
i_pre=[];
h_pre=[];
T_end=[];

clf=clf_in;
gf=gf_in;

%% Run sub-routine for specific CLF and GF
%---Variables for the initial st-st distribution--%
f_p_pro=f_p_pro_in;

comp_size(1)=gf/f_p_pro*n_t;
comp_size(2)=(1-gf*(1/f_p_pro_in+clf*ht_loss/t_c))*n_t;
comp_size(3)=clf*gf*ht_loss/t_c*n_t;

comp_size_ref(1)=gf/f_p_pro*n_t_ref;
comp_size_ref(2)=(1-gf*(1/f_p_pro_in+clf*ht_loss/t_c))*n_t_ref;
comp_size_ref(3)=clf*gf*ht_loss/t_c*n_t_ref;
%--- end ----%

%Record  number of cells
f_p=comp_size(1)/sum(comp_size);
f_i=comp_size(2)/sum(comp_size);
f_h=comp_size(3)/sum(comp_size);

for d=fx_in
    Treat_day=scheduleV;

    %Cell cycle and dose-dependent radiosensitivity
    f = @(alpha_s)F_p_cyc(1)*exp(-Alpha_ratio_p_cyc(1)*alpha_s*2-...
        Alpha_ratio_p_cyc(1)*(alpha_s/a_over_b)*4)+F_p_cyc(2)*...
        exp(-alpha_s*2-(alpha_s/a_over_b)*4)+F_p_cyc(3)*...
        exp(-Alpha_ratio_p_cyc(2)*alpha_s*2-...
        Alpha_ratio_p_cyc(2)*(alpha_s/a_over_b)*4)-exp(-alpha_p*2-...
        (alpha_p/a_over_b)*4);

    alpha_s=0.3; grid=0.1; pre_f=f(alpha_s);
    while abs(f(alpha_s)) >= eps
        if pre_f*f(alpha_s)<0
            grid=grid*0.1;
        end
        pre_f = f(alpha_s);
        if f(alpha_s)>0
            alpha_s=alpha_s+grid;
        else
            alpha_s=alpha_s-grid;
        end
    end
    Alpha_p_cyc(2)= alpha_s;

    Alpha_p_cyc(1)=Alpha_p_cyc(2)*Alpha_ratio_p_cyc(1);
    Alpha_p_cyc(3)=Alpha_p_cyc(2)*Alpha_ratio_p_cyc(2);

    %Effective alpha, beta from survival fractions
    Su_p=F_p_cyc(1)*exp(-Alpha_p_cyc(1)*d-(Alpha_p_cyc(1)/a_over_b)*d^2)...
        +F_p_cyc(2)*exp(-Alpha_p_cyc(2)*d-(Alpha_p_cyc(2)/a_over_b)*d^2)...
        +F_p_cyc(3)*exp(-Alpha_p_cyc(3)*d-(Alpha_p_cyc(3)/a_over_b)*d^2);
    alpha_p_eff=-log(Su_p)/(d*(1+(d/a_over_b)));
    beta_p_eff=(alpha_p_eff/a_over_b);

    Su_i_2gy=exp(-alpha_p/oer_i*2-(alpha_p/a_over_b)/(oer_i^2)*2^2);
    oer_i_g1=(-(Alpha_p_cyc(1)*2)-sqrt((Alpha_p_cyc(1)*2)^2-...
        4*log(Su_i_2gy)*(Alpha_p_cyc(1)/a_over_b)*2^2))/(2*log(Su_i_2gy));
    Su_h_2gy=exp(-alpha_p/oer_h*2-(alpha_p/a_over_b)/(oer_h^2)*2^2);
    oer_h_g1=(-(Alpha_p_cyc(1)*2)-sqrt((Alpha_p_cyc(1)*2)^2-...
        4*log(Su_h_2gy)*(Alpha_p_cyc(1)/a_over_b)*2^2))/(2*log(Su_h_2gy));
    alpha_i=Alpha_p_cyc(1)/oer_i_g1;
    beta_i=(Alpha_p_cyc(1)/a_over_b)/(oer_i_g1^2);
    alpha_h=Alpha_p_cyc(1)/oer_h_g1;
    beta_h=(Alpha_p_cyc(1)/a_over_b)/(oer_h_g1^2);

    alpha_p=alpha_p_eff;
    beta_p=beta_p_eff;

    % Run sub-routine for a specific CLF and GF
    %--- RT fractional dose for SBRT schedule ---%

    % Assign proliferating fraction to the initial value
    f_p_pro=f_p_pro_in;

    % Cell distribution in each compartment
    % (1:Pv, 2:Pd, 3:Iv, 4:Id, 5:Hv, 6:Hd, 7:lysis)
    % Initially all compartments are fully filled with viable cells
    % "comp_size" is the size of each compartment (1:P, 2:I, 3:H)

    cell_dist=[];
    cell_dist(1)=comp_size(1);
    cell_dist(2)=0;
    cell_dist(3)=comp_size(2);
    cell_dist(4)=0;
    cell_dist(5)=comp_size(3);
    cell_dist(6)=0;
    cell_dist(7)=0;

    % variables (t:time(day), j:# of fraction, add_time:additional time for
    %           weekend break, cum_cell_dist: cumulative cell distribution for
    %           each time increment)
    t=0;
    j=0;
    cum_cell_dist_sbrt=[];

    % Treat for specific SBRT schedule
    while t<t_start+(max(Treat_day)-1)+delta_t/2

        % Change in f_p_pro (k_p) as blood supply improves
        f_p_pro=1-0.5*(cell_dist(1)+cell_dist(2))/comp_size(1);

        % RT fraction
        if  t>(t_start+(Treat_day(j+1)-1)-delta_t/2) &&...
                t<(t_start+(Treat_day(j+1)-1)+delta_t/2)

            cell_dist(2)=cell_dist(2)+cell_dist(1)*(1-exp(-alpha_p*d-beta_p*d^2));
            cell_dist(1)=cell_dist(1)*exp(-alpha_p*d-beta_p*d^2);
            cell_dist(4)=cell_dist(4)+cell_dist(3)*(1-exp(-alpha_i*d-beta_i*d^2));
            cell_dist(3)=cell_dist(3)*exp(-alpha_i*d-beta_i*d^2);
            cell_dist(6)=cell_dist(6)+cell_dist(5)*(1-exp(-alpha_h*d-beta_h*d^2));
            cell_dist(5)=cell_dist(5)*exp(-alpha_h*d-beta_h*d^2);

            j=j+1;
        end

        % Cell Proliferation & Death
        cell_dist(1)=cell_dist(1)*(2)^(f_p_pro*delta_t/t_c);
        h_pre=cell_dist(5)+cell_dist(6);
        cell_dist(5)=cell_dist(5)*(0.5)^(delta_t/ht_loss);
        cell_dist(6)=cell_dist(6)*(0.5)^(delta_t/ht_loss);
        p_d_pre=cell_dist(2);
        cell_dist(2)=cell_dist(2)*(2)^(f_p_pro*(2*k_m-1)*delta_t/t_c);

        % Mitotically dead cell in 1 time step
        md=p_d_pre-cell_dist(2)+(h_pre-cell_dist(5)-cell_dist(6));
        cell_dist(7)=cell_dist(7)+md;
        cell_dist(7)=cell_dist(7)*(0.5)^(delta_t/ht_lys);

        % Recompartmentalization of the cell
        if cell_dist(1)+cell_dist(2)>=comp_size(1)
            p_ex=(cell_dist(1)+cell_dist(2))-comp_size(1);
            p_ratio=cell_dist(1)/(cell_dist(1)+cell_dist(2));
            cell_dist(1)=comp_size(1)*p_ratio;
            cell_dist(2)=comp_size(1)*(1-p_ratio);
            cell_dist(3)=cell_dist(3)+p_ex*p_ratio;
            cell_dist(4)=cell_dist(4)+p_ex*(1-p_ratio);
        else
            if cell_dist(3)+cell_dist(4)>0
                if cell_dist(3)+cell_dist(4)>comp_size(1)-...
                        (cell_dist(1)+cell_dist(2))
                    p_def=comp_size(1)-(cell_dist(1)+cell_dist(2));
                    i_ratio=cell_dist(3)/(cell_dist(3)+cell_dist(4));
                    cell_dist(1)=cell_dist(1)+p_def*i_ratio;
                    cell_dist(2)=cell_dist(2)+p_def*(1-i_ratio);
                    cell_dist(3)=cell_dist(3)-p_def*i_ratio;
                    cell_dist(4)=cell_dist(4)-p_def*(1-i_ratio);
                else
                    cell_dist(1)=cell_dist(1)+cell_dist(3);
                    cell_dist(2)=cell_dist(2)+cell_dist(4);
                    cell_dist(3)=0; cell_dist(4)=0;
                    if cell_dist(5)+cell_dist(6)>0
                        if cell_dist(5)+cell_dist(6)>comp_size(1)-...
                                (cell_dist(1)+cell_dist(2))
                            p_def=comp_size(1)-(cell_dist(1)+cell_dist(2));
                            h_ratio=cell_dist(5)/(cell_dist(5)+cell_dist(6));
                            cell_dist(1)=cell_dist(1)+p_def*h_ratio;
                            cell_dist(2)=cell_dist(2)+p_def*(1-h_ratio);
                            cell_dist(5)=cell_dist(5)-p_def*h_ratio;
                            cell_dist(6)=cell_dist(6)-p_def*(1-h_ratio);
                        else
                            cell_dist(1)=cell_dist(1)+cell_dist(5);
                            cell_dist(2)=cell_dist(2)+cell_dist(6);
                            cell_dist(5)=0; cell_dist(6)=0;
                        end
                    end
                end
            end
        end
        if cell_dist(3)+cell_dist(4)>=comp_size(2)
            i_ex=(cell_dist(3)+cell_dist(4))-comp_size(2);
            i_ratio=cell_dist(3)/(cell_dist(3)+cell_dist(4));
            cell_dist(3)=comp_size(2)*i_ratio;
            cell_dist(4)=comp_size(2)*(1-i_ratio);
            cell_dist(5)=cell_dist(5)+i_ex*i_ratio;
            cell_dist(6)=cell_dist(6)+i_ex*(1-i_ratio);
        else
            if cell_dist(5)+cell_dist(6)>0
                if cell_dist(5)+cell_dist(6)>comp_size(2)-...
                        (cell_dist(3)+cell_dist(4))
                    i_def=comp_size(2)-(cell_dist(3)+cell_dist(4));
                    h_ratio=cell_dist(5)/(cell_dist(5)+cell_dist(6));
                    cell_dist(3)=cell_dist(3)+i_def*h_ratio;
                    cell_dist(4)=cell_dist(4)+i_def*(1-h_ratio);
                    cell_dist(5)=cell_dist(5)-i_def*h_ratio;
                    cell_dist(6)=cell_dist(6)-i_def*(1-h_ratio);
                else
                    cell_dist(3)=cell_dist(3)+cell_dist(5);
                    cell_dist(4)=cell_dist(4)+cell_dist(6);
                    cell_dist(5)=0; cell_dist(6)=0;
                end
            end
        end

        % time step increase and store the number of cells in each compartment
        t=t+delta_t;
        cum_cell_dist_sbrt=[cum_cell_dist_sbrt cell_dist'];

    end

    s_sbrt=cell_dist(1)+cell_dist(3)+cell_dist(5);
    sf_sbrt=s_sbrt/sum(comp_size);
    ntd2=length(Treat_day)*d*(1+(d/a_over_b))/(1+(2/a_over_b));
    d_sbrt=d;
    n_frac_sbrt=length(Treat_day);
    duration_sbrt=max(Treat_day);
    t_sbrt=t;
    %-------------------------------------------------------------------------------%

    %% EQD2 calculation
    d=2;
    alpha_p=alpha_p_ori;            beta_p=beta_p_ori;
    alpha_i=alpha_p_ori/oer_i;      beta_i=beta_p_ori/(oer_i^2);
    alpha_h=alpha_p_ori/oer_h;      beta_h=beta_p_ori/(oer_h^2);
    s_eqd2=0;       sf_eqd2=0;      eqd2=0;

    %----- RT fractional dose for EQD2 estimation ----%

    % Assign proliferating fraction to the initial value

    f_p_pro=f_p_pro_in;

    % Cell distribution in each compartment
    % (1:Pv, 2:Pd, 3:Iv, 4:Id, 5:Hv, 6:Hd, 7:lysis)
    % Initially all compartments are fully filled with viable cells
    % "comp_size" is the size of each compartment (1:P, 2:I, 3:H)

    cell_dist=[];
    cell_dist(1)=comp_size_ref(1);
    cell_dist(2)=0;
    cell_dist(3)=comp_size_ref(2);
    cell_dist(4)=0;
    cell_dist(5)=comp_size_ref(3);
    cell_dist(6)=0;
    cell_dist(7)=0;

    % variables (t:time(day), j:# of fraction, add_time:additional time for
    %           weekend break, cum_cell_dist: cumulative cell distribution for
    %           each time increment)
    t=0;
    j=0;
    add_time=0;
    cum_cell_dist=[];

    % Treat until the SF becomes equivalent to SBRT regime
    while (cell_dist(1)+cell_dist(3)+cell_dist(5))>s_sbrt

        % Change in f_p_pro (k_p) as blood supply improves
        f_p_pro=1-0.5*(cell_dist(1)+cell_dist(2))/comp_size(1);

        % RT fraction
        if  t>(t_start+j+add_time-delta_t/2) && t<(t_start+j+add_time+delta_t/2)

            cell_dist(2)=cell_dist(2)+cell_dist(1)*(1-exp(-alpha_p*d-beta_p*d^2));
            cell_dist(1)=cell_dist(1)*exp(-alpha_p*d-beta_p*d^2);
            cell_dist(4)=cell_dist(4)+cell_dist(3)*(1-exp(-alpha_i*d-beta_i*d^2));
            cell_dist(3)=cell_dist(3)*exp(-alpha_i*d-beta_i*d^2);
            cell_dist(6)=cell_dist(6)+cell_dist(5)*(1-exp(-alpha_h*d-beta_h*d^2));
            cell_dist(5)=cell_dist(5)*exp(-alpha_h*d-beta_h*d^2);

            j=j+1;

            % Week-end break
            if rem(j,5)==0
                add_time=add_time+2;
            end

        end

        % Cell Proliferation & Death
        cell_dist(1)=cell_dist(1)*(2)^(f_p_pro*delta_t/t_c);
        h_pre=cell_dist(5)+cell_dist(6);
        cell_dist(5)=cell_dist(5)*(0.5)^(delta_t/ht_loss);
        cell_dist(6)=cell_dist(6)*(0.5)^(delta_t/ht_loss);
        p_d_pre=cell_dist(2);
        cell_dist(2)=cell_dist(2)*(2)^(f_p_pro*(2*k_m-1)*delta_t/t_c);

        % Mitotically dead cell in 1 time step
        md=p_d_pre-cell_dist(2)+(h_pre-cell_dist(5)-cell_dist(6));
        cell_dist(7)=cell_dist(7)+md;
        cell_dist(7)=cell_dist(7)*(0.5)^(delta_t/ht_lys);

        % Recompartmentalization of the cell
        if cell_dist(1)+cell_dist(2)>=comp_size(1)
            p_ex=(cell_dist(1)+cell_dist(2))-comp_size(1);
            p_ratio=cell_dist(1)/(cell_dist(1)+cell_dist(2));
            cell_dist(1)=comp_size(1)*p_ratio;
            cell_dist(2)=comp_size(1)*(1-p_ratio);
            cell_dist(3)=cell_dist(3)+p_ex*p_ratio;
            cell_dist(4)=cell_dist(4)+p_ex*(1-p_ratio);
        else
            if cell_dist(3)+cell_dist(4)>0
                if cell_dist(3)+cell_dist(4)>comp_size(1)-...
                        (cell_dist(1)+cell_dist(2))
                    p_def=comp_size(1)-(cell_dist(1)+cell_dist(2));
                    i_ratio=cell_dist(3)/(cell_dist(3)+cell_dist(4));
                    cell_dist(1)=cell_dist(1)+p_def*i_ratio;
                    cell_dist(2)=cell_dist(2)+p_def*(1-i_ratio);
                    cell_dist(3)=cell_dist(3)-p_def*i_ratio;
                    cell_dist(4)=cell_dist(4)-p_def*(1-i_ratio);
                else
                    cell_dist(1)=cell_dist(1)+cell_dist(3);
                    cell_dist(2)=cell_dist(2)+cell_dist(4);
                    cell_dist(3)=0; cell_dist(4)=0;
                    if cell_dist(5)+cell_dist(6)>0
                        if cell_dist(5)+cell_dist(6)>comp_size(1)-...
                                (cell_dist(1)+cell_dist(2))
                            p_def=comp_size(1)-(cell_dist(1)+cell_dist(2));
                            h_ratio=cell_dist(5)/(cell_dist(5)+cell_dist(6));
                            cell_dist(1)=cell_dist(1)+p_def*h_ratio;
                            cell_dist(2)=cell_dist(2)+p_def*(1-h_ratio);
                            cell_dist(5)=cell_dist(5)-p_def*h_ratio;
                            cell_dist(6)=cell_dist(6)-p_def*(1-h_ratio);
                        else
                            cell_dist(1)=cell_dist(1)+cell_dist(5);
                            cell_dist(2)=cell_dist(2)+cell_dist(6);
                            cell_dist(5)=0; cell_dist(6)=0;
                        end
                    end
                end
            end
        end
        if cell_dist(3)+cell_dist(4)>=comp_size(2)
            i_ex=(cell_dist(3)+cell_dist(4))-comp_size(2);
            i_ratio=cell_dist(3)/(cell_dist(3)+cell_dist(4));
            cell_dist(3)=comp_size(2)*i_ratio;
            cell_dist(4)=comp_size(2)*(1-i_ratio);
            cell_dist(5)=cell_dist(5)+i_ex*i_ratio;
            cell_dist(6)=cell_dist(6)+i_ex*(1-i_ratio);
        else
            if cell_dist(5)+cell_dist(6)>0
                if cell_dist(5)+cell_dist(6)>comp_size(2)-...
                        (cell_dist(3)+cell_dist(4))
                    i_def=comp_size(2)-(cell_dist(3)+cell_dist(4));
                    h_ratio=cell_dist(5)/(cell_dist(5)+cell_dist(6));
                    cell_dist(3)=cell_dist(3)+i_def*h_ratio;
                    cell_dist(4)=cell_dist(4)+i_def*(1-h_ratio);
                    cell_dist(5)=cell_dist(5)-i_def*h_ratio;
                    cell_dist(6)=cell_dist(6)-i_def*(1-h_ratio);
                else
                    cell_dist(3)=cell_dist(3)+cell_dist(5);
                    cell_dist(4)=cell_dist(4)+cell_dist(6);
                    cell_dist(5)=0; cell_dist(6)=0;
                end
            end
        end

        % time step increase and store the number of cells in each compartment
        t=t+delta_t;
        cum_cell_dist=[cum_cell_dist cell_dist'];

        s_eqd2_pre=s_eqd2;
        sf_eqd2_pre=sf_eqd2;
        eqd2_pre=eqd2;

        s_eqd2=cell_dist(1)+cell_dist(3)+cell_dist(5);
        sf_eqd2=s_eqd2/sum(comp_size);
        tcp=exp(-s_eqd2*f_s);
        eqd2=j*d;

    end

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

    eqd2=eqd2_pre+((eqd2-eqd2_pre)/(s_eqd2_pre-s_eqd2))*(s_eqd2_pre-s_sbrt);

end

%% Compute TCP
TD_50 = 62.1;
gamma_50 = 1.5;
TCP_upper_bound = 0.95;

TCP=TCP_upper_bound/(1+(TD_50/eqd2)^(4*gamma_50));

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
%Autocalibration.
%
%This add on allows to compute the hand to eye calibration based ont he
%calibration toolbox from Jean-Yves Bouguet
%see http://www.vision.caltech.edu/bouguetj/calib_doc/
%more information can be found here:
%http://www.vision.ee.ethz.ch/~cwengert/calibration_toolbox.php
%
%It performs the calibration 100% automatically, but uses a different
%pattern for calibration which can be found here:
%http://www.vision.ee.ethz.ch/~cwengert/calibration_toolbox.php
%
%See also "Fully Automatic Endoscope Calibration for Intraoperative Use,
%Wengert C., Reeff M., Cattin P., Szekely G."
%
%You can set up some variables to change the behavior of the code:
%doShow	[default=0]
%Set doShow = 1 if you want to see what the calibration is doing
%doRefineGrid [default = 0]
%Set doRefineGrid = 1 if the grid should be refined after detection
%(does not really change much)
%isStronglyDistorted [default = 0]
%Set isStronglyDistorted = 1, if your images show severe distortion, such
%as from an endoscope or extremely wide optics
%doIterateCalibration [default = 1]
%Set doIterateCalibration = 1, if you want the system to give the best
%results. This will try to detect wrong correspondences and calibrate
%again, set to 0 if you dont want to use this feature.
%
%Christian Wengert
%Computer Vision Laboratory
%ETH Zurich
%Sternwartstrasse 7
%CH-8092 Zurich
%www.vision.ee.ethz.ch/cwengert
%wengert@vision.ee.ethz.ch

%Needed for the gridextractor, depending on your images, change this
minArea =15;maxArea = 14000;
%maximum backprojkection error
MAX_ERROR = 4000;

%With the doShow flag you can show the results
if(~exist('doShow'))
    doShow = 0;
end

if(~exist('dX_default'))
	dX_default = 2;
end

if(~exist('dY_default'))
	dY_default = 2;
end

if(~exist('dX'))
    dX = dX_default;
end
if(~exist('dY'))
    dY = dY_default;
end

%Give best accuracy
if(~exist('doIterateCalibration'))
    doIterateCalibration =0;
end
%Check the flag for refining the grid
if(~exist('doRefineGrid'))
    doRefineGrid = 0;
end
%Set it to zero anyway, it is not yet ready!
doRefineGrid = 0;

%Check the distortion flag
if(~exist('isStronglyDistorted'))
    isStronglyDistorted = 0;
end

%Check whether we have images
if(n_ima<=0)
    disp(['autocalibrate:: No images available']);
    return
end
    %CORNER EXTRACTION PART
    active_images = zeros(1,n_ima);
    for i=1:n_ima
        disp(['autocalibrate:: Processing image ' num2str(i)]);
        %Needs to be done for toolbox
        eval(['x_' num2str(i) ' = NaN*ones(2,1);']);
        eval(['X_' num2str(i) ' = NaN*ones(3,1);']);
        eval(['im = I_' num2str(i) ';'])
        if(doShow)  %Show if necessary
            figure;hold on;title(['Image ' num2str(i)]);
        end
        %Extract the grid points
        [success, stats, cnt, id0, id1, searchPoints] = gridextractor(im, doShow, minArea, maxArea);

        %Sometimes if we can find the two main marks, we dont process this
        %image
        if(success & cnt==2 & length(stats)>0)
            %Fit the grid
            if(isStronglyDistorted)
                [success, x_, X_, searchPoints] = fitgridDistorted(stats, searchPoints, id0, id1, dX, dY, doShow);
            else
                [success, x_, X_, searchPoints, err] = fitgrid(stats, searchPoints, id0, id1, dX, dY, doShow);
            end
            %Grid was fitted succesfully
            if(success)
                %Refine grid if wanted
                if(doRefineGrid)
                    refineGrid(im, success, stats, cnt, id0, id1, searchPoints);
                end
                %Update calib toollbox
                active_images(i) = i;
                eval(['X_' num2str(i) '=X_;']);
                eval(['x_' num2str(i) '=x_;']);

                %Check for error
                if(exist('err')==1)
                    if(err>MAX_ERROR)
                        active_images(i) = 0;
                        disp(['autocalibrate:: Image ' num2str(i) ' has rather huge gridprojection error. Setting it to inactive.']);
                    end
                end
            end
        else
            disp(['autocalibrate:: Not enough image information in image ' num2str(i) '. Setting it to inactive.']);
        end
    end
    %Make a backup
    old_active_images = active_images;
    %NOW WE START THE CAMERA CALIBRATION
    %Do this for calibtoolbox compatibility
    [ny,nx,bpp] = size(im);
    center_optim=1;
    est_aspect_ratio = 1;
    est_alpha=0;
    est_dist = [1 1 1 1 0]';
    init_intrinsic_param;
    %Calibrate camera
    try
        while(~exist('err_std'))
            go_calib_optim_iter;
        end
    catch
        disp(['autocalibrate:: Could not calibrate']);
        s = lasterror;
        disp(['autocalibrate:: ' s.message ' in function ' s.stack.name ' in ' s.stack.file ' @line ' num2str(s.stack.line)])
    end

    if(doIterateCalibration)
        %Improve
        maxIter=10;
        iter=0;

        %Max Error to eliminate points which have a backprojectionerror bigger than
        %that
        BP_ERROR = 1.5;
        while(norm(err_std)>0.5 & iter<maxIter)
            for i=active_images
                if(i>0)
                    eval(['error = ex_' num2str(i) ';']);
                    ix=find(error(1,:)>BP_ERROR);
                    iy=find(error(2,:)>BP_ERROR);
                    ixx = union(ix,iy);
                    eval(['x_' num2str(i) '(:,ixx) = [];']);
                    eval(['X_' num2str(i) '(:,ixx) = [];']);
                end

            end
            go_calib_optim_iter;
           iter=iter+1;
        end
    end

author: christianwengert


text:
function poly = rowToPolygon(row, varargin)
%ROWTOPOLYGON  Create a polygon from a row vector.
%
%   POLY = rowToPolygon(ROW)
%   Convert a 1-by-2*N row vector that concatenates all polygon vertex
%   coordinates into a N-by-2 array of coordinates.
%   Default ordering of coordinates in ROW is:
%   [X1 Y1 X2 Y2 X3 Y3 .... XN YN].
%
%   POLY = rowToPolygon(ROW, METHOD)
%   Specifies the method for concatenating coordinates. METHOS is one of:
%   'interlaced': default method, described above.
%   'packed': the vector ROW has format:
%   [X1 X2 X3 ... XN Y1 Y2 Y3 ... YN].
%
%   POLYS = rowToPolygon(ROWS, ...)
%   When ROWS is a NP-by-NV array containing the vertex coordinates of NP
%   polygons, returns a 1-by-NP cell array containing in each cell the
%   coordinates of the polygon.
%
%
%   Example
%   % Concatenate coordinates of a circle and draw it as a polygon
%     t = linspace (0, 2*pi, 200);
%     row = [cos(t) sin(t)];
%     poly = rowToPolygon(row, 'packed');
%     figure;drawPolygon(poly)
%
%   See also
%   polygons2d, polygonToRow

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

type = 'interlaced';
if ~isempty(varargin)
    type = varargin{1};
end

% number of polygons
nPolys = size(row, 1);

% polygon vertex number
Np = size(row, 2) / 2;

if strcmp(type, 'interlaced')
    % ordering is [X1 Y1 X2 X2... XN YN]
    if nPolys == 1
        poly = reshape(row, [2 Np])';
    else
        poly = cell(1, nPolys);
        for i = 1:nPolys
            poly{i} = reshape(row(i,:), [2 Np])';
        end
    end

elseif strcmp(type, 'packed')
    % ordering is [X1 X2 X3... XN Y1 Y2 Y3... YN]
    if nPolys == 1
        poly = [row(1:Np)' row(Np+1:end)'];
    else
        poly = cell(1, nPolys);
        for i = 1:nPolys
            poly{i} = [row(i, 1:Np)' row(i, Np+1:end)'];
        end
    end
end

author: mattools

repo: matGeom

save_path: github-repos/MATLAB/mattools-matGeom


text:
function [x_preds, y_preds] = load_pred_result(mat_folder_path, img_basename)
   % Load pred mat file for each test image from LSP dataset.
   mat_name = strcat(img_basename, '.mat');
   mat_path = strcat(mat_folder_path, '/', mat_name);
   %disp(mat_path);
   pred_struct = load(mat_path);
   x_preds = pred_struct.joints(:, 2);
   y_preds = pred_struct.joints(:, 1);
   %disp(x_preds);
end

author: Guanghan

repo: GNet-pose

save_path: github-repos/MATLAB/Guanghan-GNet-pose


text:
function varargout = beale(X)
% Beale funcion
%
%   BEALE([x1, x2]) returns the value of the value of the Beale
%   function at the specified points. [x1] and [x2] may be vectors.
%   The search domain is
%
%               -4.5 < x_i < 4.5
%
%   The global minimum is
%
%               f(x1, x2) = f(3, 0.5) = 0.

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

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

    % otherwise, output function value
    else

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

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

        % output function value
        varargout{1} = (1.5 - x1 + x1.*x2).^2 + (2.25 - x1 + x1.*x2.^2).^2 + (2.625 - x1 + x1.*x2.^3).^2;

    end

end


text:
function [x, info] = analyticCenter(f, x, opts)
% x = analyticCenter(f, opts, x)
% compute the analytic center for the domain {Ax=b} intersect the domain of f
%
% Input:
%    f - a ConvexProgram
%    x - a feasible initial point (optional)
%    opts - a structure for options with the following properties (optional)
%       MaxIter - maximum number of iterations
%       Output - maximum number of iterations
%       CentralityTol, FeasibilityTol - stop the following are satisfied
%           ||(A' * lambda - grad f(x)) / sqrt(hess)||_inf < CentralityTol
%           ||A x - b||_inf < FeasibilityTol
%       CollapseDistanceTol -
%           if x_i is CollapseDistanceTol close to some boundary,
%           we assume the block i is tight.
%       SolverIter - the number of iteration in solving linear systems
%       VectorType - the class we use for x (@double, @ddouble, @qdouble)
%
% Output:
%  x - It outputs the analytic center of f
%  info - a structure containing centrality, feasibility and iter.
%  f - the problem f will be modified as we discover collapsed subspace

defaultOpts = struct('MaxIter', 1000, 'Output', @disp, 'CentralityTol', 1e-8, 'FeasibilityTol', 1e-12, ...
            'CollapseDistanceTol', 1e-8, 'SolverIter', 3, 'VectorType', @ddouble);
if nargin >= 3
   opts = setField(defaultOpts, opts);
else
   opts = defaultOpts;
end

if f.feasible == false
   x = []; info = [];
   return;
end

%% prepare the printout
output = TableDisplay('Iter', '5i', 'PredStep', '13.2e', 'CorrStep', '13.2e', 'Centrality', '13.2e', 'Feasibility', '13.2e', 'mu', '13.2e');
output.output = opts.Output;
output.header();

%% initial parameters
f.removeRedundantRows();
if nargin <= 1 || isempty(x) || any(f.distance(x) <= 0)
   % a heuristic initial point
   f.solver.factorize(ones(size(f.A,2),1));
   v = f.A' * f.solver.solve(f.b);
   x = f.barrier.center;
   t = stepSize(x, v-x, 0.95);
   x = x + t * (v-x);
end

A = f.A; At = A'; b = f.b; rx = [];
x = opts.VectorType(x);
y = 0*b;
feasibilityMode = true;
mu = 0;
lastProgress = 0; % record the last iteration making progress
bestCentrality = 1e32; bestFeasibility = 1e32;

%% find the central path
for iter = 1:opts.MaxIter
   % Compute the residual
   if (isempty(rx))
      [rs, rx, h] = updateResidual(x, y, mu);
   end

   % Compute the cholesky decomposition
   cholErr = f.solver.factorize(1./h);

   % Compute the direction
   if (cholErr < f.solver.cholTol)
      v = f.solver.solve([A*(rs./h) rx], [], opts.SolverIter);
      Atv = At * v;
      y = y - v(:,1);

      if (feasibilityMode)
         % corr_dx = (rs - At * (R\(R'\(A*(rs./hess)))))./hess;
         % pred_dx = (At * (R\(R'\rx)))./hess;

         corr_dx = (rs - Atv(:,1))./h;
         pred_dx = -Atv(:,2)./h;

         pred_t = stepSize(x, pred_dx, 1.0);
         dx = corr_dx + pred_t * pred_dx;
         t = stepSize(x, dx, 0.95);
         x = x + t * dx;
      else
         % y = y - (R\(R'\(A*(rs./hess))))
         % dx = (rs + At * (R\(R'\(rx - A*(rs./hess)))))./hess;

         pred_t = 1.0;
         dx = (rs - Atv(:,1))./h - Atv(:,2)./h;
         t = stepSize(x, dx, 0.95);
         x = x + t * dx;
      end

      [rs, rx, h] = updateResidual(x, y, mu);

      centrality = max(abs(rs./sqrt(h)));
      feasibility = max(abs(rx));
      if isempty(feasibility), feasibility = 0; end % Fix the case rx is []

      % Output the error
      o = struct('Iter', iter, 'PredStep', t * pred_t, 'CorrStep', t, 'Centrality', centrality, 'Feasibility', feasibility, 'mu', mu);
      output.row(o);

      if (centrality < 0.9 * bestCentrality)
         bestCentrality = centrality;
         if ~feasibilityMode, lastProgress = iter; end
      end

      if (feasibility < 0.9 * bestFeasibility)
         bestFeasibility = feasibility;
         if feasibilityMode, lastProgress = iter; end
      end

      % Check stop criteria
      if centrality < opts.CentralityTol && feasibility < opts.FeasibilityTol
         if feasibilityMode && (~isempty(f.c) || ~isempty(f.df))
            feasibilityMode = false;
            mu = 1;
            [rs, rx, h] = updateResidual(x, y, mu);
            opts.Output('Switch to optimize mode.');

            bestCentrality = 1e32; bestFeasibility = 1e32;
         else
            break;
         end
      end
   end

   % Perform the collpase if needed
   dist = f.distance(x);
   if any(dist < opts.CollapseDistanceTol / 10)
      blocks = find(dist < opts.CollapseDistanceTol);

      % Update barrier
      block_idx = f.collapse(x, blocks);

      % check feasibility after collapsing
      if f.feasible == false
         x = []; info = [];
         opts.Output('The problem is infeasible.');
         return;
      end

      % Update other variables
      A = f.A; At = A'; b = f.b;
      x(block_idx) = [];
      y = 0*b;
      rx = [];

      lastProgress = iter;
      bestCentrality = 1e32; bestFeasibility = 1e32;
      opts.Output(sprintf('Collapsed %i blocks.', length(blocks)));
      continue;
   end

   if (cholErr >= f.solver.cholTol)
      opts.Output('Failed due to numerical issue.');
      f.feasible = false;
      x = []; info = [];
      return;
   end

   if (iter > lastProgress + 20)
      opts.Output('Stopped due to no progress.');
      break;
   end
end
x = double(x);
info = struct('centrality', centrality, 'feasibility', feasibility, 'iter', iter, 'hess', h);

function t = stepSize(x, dx, factor)
   t = min(double(factor*min(f.barrier.distance(x, dx))),1); % min(Nan,1) = 1 for double
end

function [rs, rx, h] = updateResidual(x, y, mu)
   [gb, hb] = f.barrier.derivatives(x);

   gc = zeros(size(gb), class(gb));
   hc = zeros(size(hb), class(hb));

   if ~isempty(f.c)
      gc = gc + f.c;
   end

   if ~isempty(f.df)
      gc_ = f.df(f.export(x));
      gc = gc + f.scale .* gc_(f.idx);

      hc_ = f.ddf(f.export(x));
      hc = hc + f.scale2 .* hc_(f.idx);
   end

   g = gb + mu * gc; h = hb + mu * hc;

   rs = At * y - g;
   rx = A * x - b;
end

end

author: opencobra

repo: cobratoolbox


text:
%% RATE OF CONVERGENCE OF LINEAR FINITE ELEMENT METHOD
%
% This example is to show the rate of convergence of linear finite element
% approximation of the Poisson equation on the unit square:
%
% $$- \Delta u = f \; \hbox{in } (0,1)^2$$
%
% for the following boundary condition:
%
% # Non-empty Dirichlet boundary condition. $u=g_D \hbox{ on }\Gamma_D, \quad \nabla u\cdot n=g_N \hbox{ on }\Gamma_N. \Gamma _D = \{(x,y): x=0, y\in [0,1]\}, \; \Gamma _N = \partial \Omega \backslash \Gamma _D$.
% # Pure Neumann boundary condition. $\Gamma _N = \partial \Omega$.
% # Robin boundary condition. $g_R u + \nabla u\cdot n=g_N \hbox{ on }\partial \Omega$

%%
clear all; close all;
[node,elem] = cubemesh([0,1,0,1,0,1],0.25);
pde = sincosdata3;
option.L0 = 0;
option.maxIt = 4;
option.printlevel = 1;
option.elemType = 'CR';
option.plotflag = 0;

%% Non-empty Dirichlet boundary condition.
bdFlag = setboundary3(node,elem,'Dirichlet','~(x==0)','Neumann','x==0');
% bdFlag = setboundary3(node,elem,'Dirichlet');
femPoisson3(node,elem,pde,bdFlag,option);

%% Pure Neumann boundary condition.
% pde = sincosNeumanndata;
bdFlag = setboundary3(node,elem,'Neumann');
femPoisson3(node,elem,pde,bdFlag,option);

%% Pure Robin boundary condition.
pdeRobin = sincosRobindata3;
bdFlag = setboundary3(node,elem,'Robin');
femPoisson3(node,elem,pdeRobin,bdFlag,option);

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

repo: Aero-matlab


text:
function relaxPreprocess(dataDir, outDir, mtOffsetFreqs, excludeScans)
%
% relaxPreprocess(dataDir, outDir, mtOffsetFreqs, excludeScans)
%
% This file computes the T1 map (in seconds) of data obtained using
% the TOF sequence. Also saved are:
%   S0: unsaturated reference map
%   PD: a map that includes spin-density (M0), scanner scaling constant
%       (G), and T2*: PD = M0 * G * exp(-TE / T2*).
%
% SEE ALSO:
%
% relaxAlignRefToAnat.m to align the reference image output by this
% function to a structural anatomy template image.
%
% relaxMtFit.m to fit the f and k maps to the output of this function.
%
% HISTORY:
% 2007.02.?? Nikola Stikov wrote it.
% 2007.03.01 RFD: restructured the code and added auto-detection of
% most scan params.

if(~exist('dataDir','var')||isempty(dataDir))
  dataDir = fullfile(pwd,'raw');
end
if(~exist('outDir','var')||isempty(outDir))
  outDir = fileparts(dataDir);
  if(isempty(outDir)) outDir = pwd; end
end
if(~exist('mtOffsetFreqs','var')||isempty(mtOffsetFreqs))
  mtOffsetFreqs = [];
  error('Auto-detection of the mtOffsetFreqs is not implemented yet!');
end
if(~exist('excludeScans','var')||isempty(excludeScans))
  excludeScans = [];
end

showFigs = false;

s = dicomLoadAllSeries(dataDir);

% Exclude scans that don't match most scans imSize
for(ii=1:length(s))
  sz(ii) = prod(size(s(ii).imData));
end
badScans = find(sz~=median(sz));
if(~isempty(badScans))
  disp(['Adding ' num2str(badScans) ' to exclude list due to size mis-match.']);
  excludeScans = [excludeScans badScans];
end

%% Apply exclusion list
%
if(~isempty(excludeScans))
  s(excludeScans) = [];
end

if(showFigs)
    for(ii=1:length(s))
        showMontage(s(ii).imData);
        set(gcf,'name',[s(ii).desc '(' num2str(s(ii).seriesNum) ')']);
    end
end

% Align all the images
disp('Aligning all series to the first...');
bb = [1 1 1; size(s(1).imData)];
for(ii=2:length(s))
    xform = mrAnatRegister(s(ii).imData,s(1).imData);
    s(ii).imData =  mrAnatResliceSpm(s(ii).imData, xform, bb, [1 1 1], [7 7 7 0 0 0], 0);
    s(ii).imData(isnan(s(ii).imData)|s(ii).imData<0) = 0;
end

%% Process the T1 relaxometry scans
%
t1Inds = find([s.mtOffset]==0);
nT1 = length(t1Inds);

meanOfT1s = mean(cat(4,s(t1Inds).imData),4);
%showMontage(mn);
[brainMask,checkSlices] = mrAnatExtractBrain(meanOfT1s,[],0.25);
%figure;image(checkSlices);

% We save as we go
xform = s(t1Inds(1)).imToScanXform;
dtiWriteNiftiWrapper(uint8(brainMask), xform, fullfile(outDir,'brainMask'));
dtiWriteNiftiWrapper(meanOfT1s, xform, fullfile(outDir,'ref'));
clear meanOfT1s;

brainInds = find(brainMask);
nVox = length(brainInds);

%% Process MT scans
%
% For these, we average multiple repeats/
mtInds = find([s.mtOffset]~=0);
nMt = length(mtInds);

% for(ii=1:length(mtInds)), for(jj=1:length(mtInds))
% 	rms(ii,jj) = sqrt(mean((s(mtInds(ii)).imData(brainInds)-s(mtInds(jj)).imData(brainInds)).^2));
% end, end
mtOffsets = unique(mtOffsetFreqs);
for(ii=1:length(mtOffsets))
    mt = mean(cat(4,s(mtInds(mtOffsetFreqs==mtOffsets(ii))).imData),4);
    dtiWriteNiftiWrapper(mt, xform, fullfile(outDir,sprintf('MT_%02dkHz',mtOffsets(ii))));
end
clear mt;

theta = [s(t1Inds).flipAngle]*pi/180;
x = zeros(nVox,nT1);
y = zeros(nVox,nT1);
for(ii=1:nT1)
  y(:,ii) = abs(s(t1Inds(ii)).imData(brainInds)./sin(theta(ii)));
  x(:,ii) = abs(s(t1Inds(ii)).imData(brainInds)./tan(theta(ii)));
end

% Following is a slow loop. Can it be vectorized?
disp('Fitting T1 estimates for each voxel (SLOW!)');
t1 = zeros(nVox,1);
tic;
for(ii=1:nVox)
  if(mod(ii,20000)==0)
    fprintf('  Processed %d of %d voxels (%0.1f %%) in %0.1f secs...\n',ii,nVox,ii/nVox*100,toc);
  end
  if(max(y(ii,:)>0.001))
    d = polyfit (x(ii,:), y(ii,:), 1);
    t1(ii) = d(1);
  end
end

% Why do an abs here? All the negative values seem to be junk.
%t1 = abs(t1);
t1(t1<0) = NaN;
T1 = (-log(t1)/20e-3).^-1;
% Clip to plausible values
T1(isnan(T1)|T1<0.2) = 0.2;
T1(T1>20) = 20;
% Get a PD map for each flip angle
PDmap = zeros(nVox,nT1);
for(ii=1:nT1)
  PDmap(:,ii) = s(t1Inds(ii)).imData(brainInds)./sin(theta(ii)).*((1-cos(theta(ii)).*t1)./(1-t1));
end
% Could do least squares for better PDmap?
PD = mean(PDmap,2);
% WARNING: following clip values are empirically determined with no
% theoretical basis!
PD(PD<0) = 0;
PD(PD>12000) = 12000;

% compute the S0 for the flip angle used in the MT scans
theta = s(mtInds(1)).flipAngle*pi/180;
S0 = PD.*sin(theta).*(1-exp(-32e-3./T1))./(1-cos(theta).*exp(-32e-3./T1));

im=zeros(size(brainMask)); im(brainInds) = T1; T1 = im;
im=zeros(size(brainMask)); im(brainInds) = PD; PD = im;
im=zeros(size(brainMask)); im(brainInds) = S0; S0 = im;

dtiWriteNiftiWrapper(T1, xform, fullfile(outDir,'T1'));
dtiWriteNiftiWrapper(PD, xform, fullfile(outDir,'PD'));
dtiWriteNiftiWrapper(S0, xform, fullfile(outDir,'S0'));

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% EX_MAXWELL_SRC_CUBE: solve Maxwell source problem in the unit cube.

% 1) PHYSICAL DATA OF THE PROBLEM
clear problem_data
% Physical domain, defined as NURBS map given in a text file
problem_data.geo_name = 'geo_cube.txt';

% Type of boundary conditions
problem_data.nmnn_sides   = [3 4 5 6];
problem_data.drchlt_sides = [1 2];

% Physical parameters
problem_data.c_mass  = @(x, y, z) ones(size(x));
problem_data.c_stiff = @(x, y, z) ones(size(x));

% Source and boundary terms
problem_data.f = @(x, y, z) cat(1, ...
                    reshape (sin(y) .* (2*z - exp(x)) - exp(z) .* sin(x), [1, size(x)]), ...
                    zeros ([1, size(x)]), ...
                    reshape (2 * exp(z) .* cos(x), [1, size(x)]));
problem_data.g = @(x, y, z, ind) test_maxwell_cube_g_nmnn (x, y, z, ind);
problem_data.h = @(x, y, z, ind) cat(1, ...
                    reshape (sin(y) .* z, [1, size(x)]), ...
                    reshape (exp(x) .* cos(y), [1, size(x)]), ...
                    reshape (exp(z) .* cos(x), [1, size(x)]));

% Exact solution (optional)
problem_data.uex     = @(x, y, z) cat(1, ...
                    reshape (sin(y) .* z, [1, size(x)]), ...
                    reshape (exp(x) .* cos(y), [1, size(x)]), ...
                    reshape (exp(z) .* cos(x), [1, size(x)]));
problem_data.curluex = @(x, y, z) cat(1, ...
                    zeros ([1, size(x)]), ...
                    reshape (sin(y) + exp(z) .* sin(x), [1, size(x)]), ...
                    reshape (exp(x).*cos(y) - z .* cos(y), [1, size(x)]));

% 2) CHOICE OF THE DISCRETIZATION PARAMETERS
clear method_data
method_data.degree     = [2 2 2];     % Degree of the bsplines
method_data.regularity = [1 1 1];     % Regularity of the splines
method_data.nsub       = [3 3 3];     % Number of subdivisions
method_data.nquad      = [3 3 3];     % Points for the Gaussian quadrature rule

% 3) CALL TO THE SOLVER
[geometry, msh, space, u] = solve_maxwell_src (problem_data, method_data);

% 4) POST-PROCESSING
% 4.1) EXPORT TO PARAVIEW
output_file = 'maxwell_cube_Deg2_Reg1_Sub3';

vtk_pts = {linspace(0, 1, 15), linspace(0, 1, 15), linspace(0, 1, 15)};
fprintf ('The result is saved in the file %s \n \n', output_file);
sp_to_vtk (u, space, geometry, vtk_pts, output_file, 'u')

% 4.2) Comparison with the exact solution
[error_hcurl, error_l2] = ...
    sp_hcurl_error (space, msh, u, problem_data.uex, problem_data.curluex)

%!demo
%! ex_maxwell_src_cube

%!test
%! problem_data.geo_name = 'geo_cube.txt';
%! problem_data.nmnn_sides   = [3 4 5 6];
%! problem_data.drchlt_sides = [1 2];
%! problem_data.c_mass  = @(x, y, z) ones(size(x));
%! problem_data.c_stiff = @(x, y, z) ones(size(x));
%! problem_data.f = @(x, y, z) cat(1, ...
%!                     reshape (sin(y) .* (2*z - exp(x)) - exp(z) .* sin(x), [1, size(x)]), ...
%!                     zeros ([1, size(x)]), ...
%!                     reshape (2 * exp(z) .* cos(x), [1, size(x)]));
%! problem_data.g = @(x, y, z, ind) test_maxwell_cube_g_nmnn (x, y, z, ind);
%! problem_data.h = @(x, y, z, ind) cat(1, ...
%!                     reshape (sin(y) .* z, [1, size(x)]), ...
%!                     reshape (exp(x) .* cos(y), [1, size(x)]), ...
%!                     reshape (exp(z) .* cos(x), [1, size(x)]));
%! problem_data.uex     = @(x, y, z) cat(1, ...
%!                     reshape (sin(y) .* z, [1, size(x)]), ...
%!                     reshape (exp(x) .* cos(y), [1, size(x)]), ...
%!                     reshape (exp(z) .* cos(x), [1, size(x)]));
%! problem_data.curluex = @(x, y, z) cat(1, ...
%!                     zeros ([1, size(x)]), ...
%!                     reshape (sin(y) + exp(z) .* sin(x), [1, size(x)]), ...
%!                     reshape (exp(x).*cos(y) - z .* cos(y), [1, size(x)]));
%! method_data.degree     = [2 2 2];     % Degree of the bsplines
%! method_data.regularity = [1 1 1];     % Regularity of the splines
%! method_data.nsub       = [3 3 3];     % Number of subdivisions
%! method_data.nquad      = [3 3 3];     % Points for the Gaussian quadrature rule
%! [geometry, msh, space, u] = solve_maxwell_src (problem_data, method_data);
%! [error_hcurl, error_l2] = sp_hcurl_error (space, msh, u, problem_data.uex, problem_data.curluex);
%! assert (msh.nel, 27)
%! assert (space.ndof, 300)
%! assert (error_l2, 0.00897255476928207, 1e-14)
%! assert (error_hcurl, 0.0131952570451473, 1e-14)

author: rafavzqz

repo: geopdes

save_path: github-repos/MATLAB/rafavzqz-geopdes


text:
%% contour2logic
% Below is a demonstration of the features of the |contour2logic| function

%%
clear; close all; clc;

%% Syntax
% |[varargout]=contour2logic(M,v,Vcs);|

%% Description
% This function converts contours to a logic or labelled data. The logics
% represent wether voxels are in, on, or outside the contour.
%
% The input consists of:
%%
%
% * A 3D image |M| (or alternatively the size of M).
% * The |vozelSize| a 1x3 vector specifying the size of the voxels in the
% row, column, and slice direction.
% * A cell array |Vcs| containing one or more contours per slice. If the
% image has n slices then Vcs should be an nx1 cell array, i.e. contours
% are defined for each slice.

%% Examples
%
%% Import image data for this demo

defaultFolder = fileparts(fileparts(mfilename('fullpath'))); %Set main folder
pathNameImageData=fullfile(defaultFolder,'data','DICOM','0001_human_calf');
loadNameImageData=fullfile(pathNameImageData,'IMDAT','IMDAT.mat');
IMDAT_struct=load(loadNameImageData); %The image data structure
G = IMDAT_struct.G; %Geometric/spatial information
v=G.v; %The voxel size
M= IMDAT_struct.type_1; %The image data

%%
contourName='imseg_calf_tibia';
pathName=fullfile(defaultFolder,'data','imseg'); %Folder name for contours

%% Compute levelset

loadName=fullfile(pathName,contourName);
load(loadName); %Load segmentation structure
Vcs=saveStruct.ContourSet; %Access the contour data

[logicIn,logicOn,N]=contour2logic(M,v,Vcs);

%%
% Visualize logic image and contours together

%Visualize logic image
sv3(logicIn,v); %Open slice viewer for levelset

%Visualize contours
optionStruct.Color='r';
plotContours({Vcs},optionStruct);  %Plot contours

%Add colorbar labels
[~,hc]=icolorbar;
hc.TickLabels={'Out','In'};

drawnow;

%%
% Visualize labelled image and contours together

%Visualize label image
vizStruct.colormap=viridis(3); %Set colormap for levelset visualization
hf2=sv3(N,v,vizStruct); %Open slice viewer for levelset

%Visualize contours
optionStruct.Color='r';
plotContours({Vcs},optionStruct);  %Plot contours

%Add colorbar labels
[~,hc]=icolorbar;
hc.TickLabels={'Out','On','In'};

drawnow;

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