text:
function tfrqview(tfr,sig,t,method,p1,p2,p3,p4,p5);
%TFRQVIEW Quick visualization of time-frequency representations.
%       TFRQVIEW(TFR,SIG,T,METHOD,P1,P2,P3,P4,P5) allows a quick
%       visualization of a time-frequency representation.
%
%       TFR     : time-frequency representation (MxN).
%       SIG     : signal in time. If unavailable, put sig=[] as input
%                 parameter.                    (default : []).
%       T       : time instants                 (default : 1:N).
%       METHOD  : name of chosen representation (default : 'TYPE1').
%                 See the TFR* files for authorized names.
%                 TYPE1 : the representation TFR goes in normalized
%                       frequency from -0.5 to 0.5 ;
%                 TYPE2 : the representation TFR goes in normalized
%                       frequency from 0 to 0.5.
%       P1...P5 : optional parameters of the representation : run the
%                 file TFRPARAM(METHOD) to know the meaning of P1..P5
%                 for your method.
%
%	When you use the 'save' option in the main menu, you save all your
%	variables as well as two strings, TfrQView and TfrView, in a mat
%	file. If you load this file and do eval(TfrQView), you will restart
%	the display session under tfrqview ; if you do eval(TfrView), you
%	will obtain the exact layout of the screen you had when clicking on
%	the 'save' button.
%
%       Example :
%        sig=fmsin(128); tfr=tfrwv(sig);
%        tfrqview(tfr,sig,1:128,'tfrwv');
%
%       See also TFRVIEW, TFRSAVE, TFRPARAM.

%       F. Auger, September 1994, July 1995
%       O. Lemoine, Oct 1995, May-July 1996.
%       F. Auger, May 1998.
%       Copyright (c) 1996 by CNRS (France).
%
%  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 2 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, write to the Free Software
%  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

comp=computer;  %  so as to know the running computer
MatlabVersion=version; MatlabVersion=str2num(MatlabVersion(1));

% Tests on the input arguments
if nargin<1,
 error('At least one parameter required'); % at least the tfr
end
[tfrrow,tfrcol]=size(tfr);
if nargin==1,
 sig=[]; t=1:tfrcol; method='type1'; % empty signal
elseif nargin==2,
 t=1:tfrcol; method='type1';
elseif nargin==3,
 method='type1';
end
[trow,tcol] = size(t);
if (trow~=1),
 error('T must only have one row'); % t must be a row vector
end;
if (tfrcol~=tcol),
 error('T must have as much elements as tfr has columns');
end;

[Nsig,Ncol]=size(sig);
if Ncol>2,
 error('SIG must have one or two columns');
end

% Computation of Nf2, the number of interresting points in frequency
method=upper(method);
if istfr2(method),
 Nf2=tfrrow;
elseif istfr1(method),
 Nf2=tfrrow/2;
else
 error('Unknown representation. Use type1 or type2');
end;

% Computation of freq (vector of frequency samples)
if istfraff(method),
 freq=eval(['p',num2str(nargin-4)]);   	% last input argument is freqs.
else
 freq=(0.5*(0:Nf2-1)/Nf2);
end

% Initialization of the variables
if exist('options.mat'),
 load options
 colormap(SavedColorMap);
else
 threshold=5.0;  % visualization threshold
 linlogtfr=0;    % tfr visualization scale : 0 for linear 1 for logarithmic
 linlogspec=1;   % spectrum visualization scale
 sigenveloppe=0; % signal enveloppe visualization

 levelnumb=64;   % number of levels in the contour plot
 colmap=1;       % colormap index

 display=2;      % display index

 isgridsig=0;    % grid on signal
 isgridspec=0;   % grid on spectrum
 isgridtfr=0;    % grid on tfr

 issig=0;        % display signal
 isspec=0;       % display spectrum
 iscolorbar=0;   % display colorbar

 fs=1.0;         % sampling frequency (Hz)
 fmin=0.0;       % smallest displayed frequency
 fmax=0.5*fs;    % highest displayed frequency
end;

% Test of analycity
if ~isempty(sig),
 for k=1:Ncol,
  % spec(:,k)=abs(fft(sig(min(t):max(t),k))).^2; Nsp=length(spec);
  % modifications :  F. Auger (fog), 30/11/97
  Lt_fog=max(t)-min(t)+1;
  Nb_tranches_fog = floor(Lt_fog/tfrrow);
  % fprintf('%f \n',Nb_tranches_fog);
  spec(:,k)=zeros(tfrrow,1);
  for Num_tranche_fog=0:Nb_tranches_fog-1,
   % fprintf('%f \n',Num_tranche_fog);
   spec(:,k)=spec(:,k)+abs(fft(sig(min(t)+tfrrow*Num_tranche_fog+(0:tfrrow-1),k))).^2;
  end;

  if (Lt_fog>Nb_tranches_fog*tfrrow),
   spectre_fog=fft(sig(min(t)+tfrrow*Nb_tranches_fog:max(t),k),tfrrow);
   spectre_fog=spectre_fog(:);
   spec(:,k)=spec(:,k)+abs(spectre_fog).^2;
  end;

  % spec1=sum(spec(1:tfrrow/2,k));
  % spec2=sum(spec(tfrrow/2+1:Nsp,k));
  spec1=sum(spec(1:tfrrow/2,k));
  spec2=sum(spec(tfrrow/2+1:tfrrow,k));

  if spec2>spec1/10,
   disp('Be careful : the signal is not analytic!');
  end
 end
end

% Test of reality
if (Ncol==2 & ~isreal(tfr)),
 disp('Cross distribution. As the result is complex, we display the real part.');
 tfr=real(tfr);
end

ChoiceDisplay     =  1; % All the possible values of the choice variable
ChoiceLayout      =  2;
ChoiceSampling    =  3;
ChoiceFreqBounds  =  4;
ChoiceThreshold   =  5;
ChoiceLinlog      =  6;
ChoiceRedraw      =  7;
ChoiceNewFigure   =  8;
ChoiceSaveResults =  9;
ChoiceSaveOptions = 10;
ChoicePrint       = 11;
ChoiceClose       = 12;

CallTfrView = 1; % 1 to call tfrview, 0 not to do it
RefreshFigure=1; % 1 to refresh figure every time, 0 to freeze

choice=ChoiceSampling;
while choice~=ChoiceClose,                       % while not close
 if RefreshFigure & CallTfrView,                 % Call to tfrview
  linlog=linlogtfr+2*linlogspec+4*sigenveloppe;
  isgrid=isgridsig+2*isgridspec+4*isgridtfr;
  layout=issig+isspec*2+iscolorbar*4+1;
  param = [display, linlog, threshold, levelnumb, Nf2, layout,...
           fs, isgrid, fmin, fmax];
  if (nargin<=4),
   tfrview(tfr,sig,t,method,param);
  elseif (nargin==5),
   tfrview(tfr,sig,t,method,param,p1);
  elseif (nargin==6),
   tfrview(tfr,sig,t,method,param,p1,p2);
  elseif (nargin==7),
   tfrview(tfr,sig,t,method,param,p1,p2,p3);
  elseif (nargin==8),
   tfrview(tfr,sig,t,method,param,p1,p2,p3,p4);
  elseif (nargin==9),
   tfrview(tfr,sig,t,method,param,p1,p2,p3,p4,p5);
  end;
 end;

 if (linlogtfr==0),                              % Lin/log scale of the tfr
  linlogstr='Change to a logarithmic scale';
 else
  linlogstr='Change to a linear scale';
 end;

 if (RefreshFigure==1),
  redrawstr='Don''t redraw yet';
 else
  redrawstr='Redraw now';
 end;

 % Main menu
 choice=menu ('TFRQVIEW MENU :',...
              'Change the display mode',...       % ChoiceDisplay
              'Change the display layout',...     % ChoiceLayout
              'Change the sampling frequency',... % ChoiceSampling
              'Change the frequency bounds',...   % ChoiceFreqBounds
              'Change the threshold',...          % ChoiceThreshold
              linlogstr,...                       % ChoiceLinlog
              redrawstr,...                       % ChoiceRedraw
              'New figure',...                    % ChoiceNewFigure
              'Save results',...                  % ChoiceSaveResults
              'Save options',...                  % ChoiceSaveOptions
              'Print',...                         % ChoicePrint
              'Close');                           % ChoiceClose

 if (choice==ChoiceDisplay),                      % Change the display mode

  OldDisplay=display;
  display=menu('DISPLAY MODE :',...
               'contour',...                      % 1
               'imagesc',...                      % 2
               'pcolor',...                       % 3
               'surf',...                         % 4
               'mesh',...                         % 5
               'change the color map',...         % 6
               'change the number of colors or levels',...  % 7
               'cancel');                         % 8

  if (display>=1)&(display<=5),
   CallTfrView=1;
  elseif (display==6),
   if MatlabVersion>=5,
    colmap=menu('COLOR MAP :',...
                'hsv','jet','cool','bone','gray','hot','prism',...
                'pink','colorcube','autumn','winter','spring','summer',...
                'brighten','darken','permute','spin','cancel');
    if     colmap== 1,  colormap(hsv(levelnumb));
    elseif colmap== 2,  colormap(jet(levelnumb));
    elseif colmap== 3,  colormap(cool(levelnumb));
    elseif colmap== 4,  colormap(bone(levelnumb));
    elseif colmap== 5,  colormap(gray(levelnumb));
    elseif colmap== 6,  colormap(hot(levelnumb));
    elseif colmap== 7,  colormap(prism(levelnumb));
    elseif colmap== 8,  colormap(pink(levelnumb));
    elseif colmap== 9,  colormap(colorcube(levelnumb));
    elseif colmap==10,  colormap(autumn(levelnumb));
    elseif colmap==11,  colormap(winter(levelnumb));
    elseif colmap==12,  colormap(spring(levelnumb));
    elseif colmap==13,  colormap(summer(levelnumb));
    elseif colmap==14,  brighten(+0.20);
    elseif colmap==15,  brighten(-0.10);
    elseif colmap==16,  MyMap = colormap; colormap(flipud(MyMap));
    elseif colmap==17,  spinmap;
    end
   else
    colmap=menu('COLOR MAP :',...
                'hsv','jet','cool','bone','gray','hot','prism',...
                'brighten','darken','permute','spin','cancel');
    if     colmap== 1,  colormap(hsv(levelnumb));
    elseif colmap== 2,  colormap(jet(levelnumb));
    elseif colmap== 3,  colormap(cool(levelnumb));
    elseif colmap== 4,  colormap(bone(levelnumb));
    elseif colmap== 5,  colormap(gray(levelnumb));
    elseif colmap== 6,  colormap(hot(levelnumb));
    elseif colmap== 7,  colormap(prism(levelnumb));
    elseif colmap== 8,  brighten(+0.25);
    elseif colmap== 9,  brighten(-0.25);
    elseif colmap==10,  MyMap = colormap; colormap(flipud(MyMap));
    elseif colmap==11,  spinmap;
    end
   end
   display=OldDisplay; CallTfrView=0;

  elseif (display==7),
   fprintf(' Old number of levels: %f\n',levelnumb); levelold=levelnumb;
   levelnumb=input(' New number of levels: ');
   if isempty(levelnumb),
    levelnumb=levelold; CallTfrView=0;
   else
    if levelnumb<levelold,
     CallTfrView=1;
     MyMap = colormap; MyMap=MyMap(1:levelnumb,:); colormap(MyMap);
    elseif levelnumb>levelold,
     CallTfrView=1;
     MyMap = ones(levelnumb, 3); MyMap(1:levelold,:)=colormap;
     fprintf('warning : the colormap size has been increased by identical vectors\n');
     fprintf('You should redefine the colormap\n');
    else
     CallTfrView=0;
    end
   end
   display=OldDisplay;

  elseif (display==8),
   display=OldDisplay; CallTfrView=0;
  end;

 elseif (choice==ChoiceLayout),                 % Change the display layout

  layout=1;
  if issig==0,
   SignalStr= 'display signal';
  else
   SignalStr='remove signal';
  end;

  if isspec==0,
   SpectrumStr= 'display spectrum';
  else
   SpectrumStr='remove spectrum';
  end;

  if ~issig & ~isspec,
   if isgridtfr,
    GridStr='Remove the grid';
   else
    GridStr='Add a grid';
   end;
  else
   GridStr='Grids';
  end;

  if iscolorbar==0,
   ColorbarStr='display colorbar';
  else
   ColorbarStr='remove colorbar';
  end;

  layout=menu('DISPLAY LAYOUT',...
               SignalStr,...
               SpectrumStr,...
               GridStr,...
               ColorbarStr,...
              'cancel');

  if layout==1,
   issig=~issig;
   if issig==1,
    if isempty(sig),
     disp('Impossible action : the signal is unavailable'); issig=0; CallTfrView=0;
    else
     sigenveloppe=menu('SIGNAL REPRESENTATION','signal only','signal with enveloppe')-1;
     CallTfrView=1;
    end;
   else
    isgridsig=0;
   end;
  elseif layout==2,
   isspec=~isspec;
   if isspec==1,
    if isempty(sig),
     disp('Impossible action : the signal is unavailable'); isspec=0; CallTfrView=0;
    else
     linlogspec=menu('FREQUENCY REPRESENTATION','linear scale','log scale')-1;
     CallTfrView=1;
    end;
   else
    isgridspec=0;
   end;

  elseif layout==3,

   if ~issig & ~isspec,	                 % No signal and no spectrum
    isgridtfr=1-isgridtfr;
    CallTfrView=1;

   elseif issig & ~isspec,               % A signal, no spectrum
    Grid=1;
    if ~isgridsig,
     gridsigstr='add a grid on the signal';
    else
     gridsigstr='remove the grid on the signal';
    end

    if ~isgridtfr,
     gridtfrstr='add a grid on the TFR';
    else
     gridtfrstr='remove the grid on the TFR';
    end

    Grid=menu('GRID MENU :',gridsigstr,gridtfrstr,'cancel');
    if Grid==1,
     isgridsig=1-isgridsig; CallTfrView=1;
    elseif Grid==2,
     isgridtfr=1-isgridtfr; CallTfrView=1;
    else
     CallTfrView=0;
    end

   elseif ~issig & isspec,               % No signal, a spectrum
    Grid=1;
    if ~isgridspec,
     gridspestr='add a grid on the spectrum';
    else
     gridspestr='remove the grid on the spectrum';
    end
    if ~isgridtfr,
     gridtfrstr='add a grid on the TFR';
    else
     gridtfrstr='remove the grid on the TFR';
    end
    Grid=menu('GRID MENU :',gridspestr,gridtfrstr,'Close');
    if Grid==1,
     isgridspec=1-isgridspec; CallTfrView=1;
    elseif Grid==2,
     isgridtfr=1-isgridtfr; CallTfrView=1;
    else CallTfrView=0;
    end

  else                                  % A signal and a spectrum
   Grid=1;
   if ~isgridsig,
    gridsigstr='add a grid on the signal';
   else
    gridsigstr='remove the grid on the signal';
   end
   if ~isgridspec,
    gridspestr='add a grid on the spectrum';
   else
    gridspestr='remove the grid on the spectrum';
   end
   if ~isgridtfr,
    gridtfrstr='add a grid on the TFR';
   else
    gridtfrstr='remove the grid on the TFR';
   end
   Grid=menu('GRID MENU :',gridsigstr,gridspestr,gridtfrstr,'cancel');
   if Grid==1,
    isgridsig=1-isgridsig; CallTfrView=1;
   elseif Grid==2,
    isgridspec=1-isgridspec; CallTfrView=1;
   elseif Grid==3,
    isgridtfr=1-isgridtfr; CallTfrView=1;
   else CallTfrView=0;
   end
  end

  elseif layout==4,
   iscolorbar=~iscolorbar; CallTfrView=1;
  elseif layout==5,
   CallTfrView=0;
  end;

 elseif (choice==ChoiceSampling),                   % Change the sampling frequency

  fprintf(' Old sampling frequency: %f\n',fs);
  fsold=fs; fs=input(' New sampling frequency: ');
  if isempty(fs),
   fs=fsold; CallTfrView=0;
  else
   CallTfrView=1;
  end;
 elseif (choice==ChoiceFreqBounds),                 % Change the frequency bounds
  CallTfrView=0;

     fprintf(' Old smallest normalized frequency : %f\n',fmin); fminold=fmin;
  fmin=input(' New smallest normalized frequency : ');
  if isempty(fmin),
   fmin=fminold;
  elseif fmin>0.5,
   fprintf('normalized frequency desired ! value unmodified\n');
   fmin=fminold;
  else
   CallTfrView=1;
  end;

     fprintf(' Old highest normalized frequency  : %f\n',fmax); fmaxold=fmax;
  fmax=input(' New highest normalized frequency  : ');
  if isempty(fmax),
   fmax=fmaxold;
  elseif fmax>0.5,
   fprintf('normalized frequency desired ! value unmodified\n');
   fmax=fmaxold;
  else
   CallTfrView=1;
  end;

 elseif (choice==ChoiceThreshold),                  % Change the threshold

  fprintf(' Old threshold: %f\n', threshold); throld=threshold;
  threshold=input(' New threshold: ');
  if isempty(threshold),
   threshold=throld; CallTfrView=0;
  else
   CallTfrView=1;
  end

 elseif (choice==ChoiceLinlog),                     % Change the lin/log scale of tfr

  linlogtfr=1-linlogtfr;

 elseif (choice==ChoiceRedraw),                           % redraw ?
  RefreshFigure=1-RefreshFigure;
  if RefreshFigure==1, CallTfrView=1; end;

 elseif (choice==ChoiceNewFigure),                        % new figure
  figure; CallTfrView=1;

 elseif (choice==ChoiceSaveResults),                      % Save the results

  f=freq*fs;
  Nmethod=length(method);
  if comp(1:2)=='PC',
   DefaultName=[method(4:Nmethod),num2str(Nsig),'.mat'];
   [name,PathWorkDir] = uiputfile(DefaultName, 'Save As');
  else
   DefaultName=[method(4:Nmethod),num2str(Nsig)];
   nameStr=[' Name of the MAT file [',DefaultName,'] : '];
   name=input(nameStr,'s');
   while (length(name)>8),
    disp(' The name must have less than 8 characters');
    name=input(nameStr,'s');
   end
   if isempty(name),
    name=DefaultName;
   end
   PathWorkDir='';
  end
  linlog=linlogtfr+2*linlogspec+4*sigenveloppe;
  isgrid=isgridsig+2*isgridspec+4*isgridtfr;
  param = [display,linlog,threshold,levelnumb,Nf2,layout,fs,isgrid,fmin,fmax];
  SavedColorMap=colormap;
  if (nargin<=4),
   TfrQView=['colormap(SavedColorMap); tfrqview(tfr,sig,t,method)'];
   TfrView =['clf;colormap(SavedColorMap); tfrview(tfr,sig,t,method,param)'];
   eval(['save ',PathWorkDir,name,...
         ' tfr sig t f fs method param SavedColorMap TfrView TfrQView']);
  elseif (nargin==5),
   TfrQView=['colormap(SavedColorMap); tfrqview(tfr,sig,t,method,p1)'];
   TfrView =['clf; colormap(SavedColorMap); tfrview(tfr,sig,t,method,param,p1)'];
   eval(['save ',PathWorkDir,name, ...
         ' tfr sig t f fs method param p1 SavedColorMap TfrView TfrQView']);
  elseif (nargin==6),
   TfrQView=['colormap(SavedColorMap); tfrqview(tfr,sig,t,method,p1,p2)'];
   TfrView =['clf; colormap(SavedColorMap); tfrview(tfr,sig,t,method,param,p1,p2)'];
   eval(['save ',PathWorkDir,name,...
         ' tfr sig t f fs method param p1 p2 SavedColorMap TfrView TfrQView']);
  elseif (nargin==7),
   TfrQView=['colormap(SavedColorMap); tfrqview(tfr,sig,t,method,p1,p2,p3)'];
   TfrView =['clf; colormap(SavedColorMap); tfrview(tfr,sig,t,method,param,p1,p2,p3)'];
   eval(['save ',PathWorkDir,name,...
         ' tfr sig t f fs method param p1 p2 p3 SavedColorMap TfrView TfrQView']);
  elseif (nargin==8),
   TfrQView=['colormap(SavedColorMap); tfrqview(tfr,sig,t,method,p1,p2,p3,p4)'];
   TfrView =['clf; colormap(SavedColorMap); tfrview(tfr,sig,t,method,param,p1,p2,p3,p4)'];
   eval(['save ',PathWorkDir,name,...
         ' tfr sig t f fs method param p1 p2 p3 p4 SavedColorMap TfrView TfrQView']);
  elseif (nargin==9),
   TfrQView=['colormap(SavedColorMap); tfrqview(tfr,sig,t,method,p1,p2,p3,p4,p5)'];
   TfrView =['clf; colormap(SavedColorMap); tfrview(tfr,sig,t,method,param,p1,p2,p3,p4,p5)'];
   eval(['save ',PathWorkDir,name,...
         ' tfr sig t f fs method param p1 p2 p3 p4 p5 SavedColorMap TfrView TfrQView']);
  end;
  disp(' ');
  fprintf('The file is saved in the directory %s\n',PathWorkDir);
  fprintf('under the name %s\n',name);

  fprintf('If you want to find again the exact layout of this screen, do\n');
  fprintf('load %s%s; eval(TfrView);\n\n', PathWorkDir,name);
  fprintf('If you want to restart the display session under tfrqview, do\n');
  fprintf('load %s%s; eval(TfrQView);\n',PathWorkDir,name);
  CallTfrView=0;

 elseif (choice==ChoiceSaveOptions),                 % Save options
  SavedColorMap=colormap;
  save options fs fmin fmax threshold linlogtfr linlogspec levelnumb  ...
       display layout colmap SavedColorMap iscolorbar  ...
       isgridsig isgridspec isgridtfr issig sigenveloppe isspec;
  fprintf('\n Options saved\n');
  CallTfrView=0;

 elseif (choice==ChoicePrint),	                     % Print the current figure

  Nmethod=length(method);
  TFTBDevice = MENU('Choose a device',...
                    '-deps','-depsc','-deps2','-depsc2','-djpeg','-dtiff','other','cancel');
  if TFTBDevice==1,
   TFTBDeviceName='-deps '  ; TFTBExtension='.eps';
  elseif TFTBDevice==2,
   TFTBDeviceName='-depsc ' ; TFTBExtension='.eps';
  elseif TFTBDevice==3,
   TFTBDeviceName='-deps2 ' ; TFTBExtension='.eps';
  elseif TFTBDevice==4,
   TFTBDeviceName='-depsc2 '; TFTBExtension='.eps';
  elseif TFTBDevice==5,
   TFTBDeviceName='-djpeg ' ; TFTBExtension='.jpg';
  elseif TFTBDevice==6,
   TFTBDeviceName='-dtiff ' ; TFTBExtension='.tif';
  elseif TFTBDevice==7,
   TFTBDeviceName=input('device option : ','s'); ; TFTBDeviceName=[TFTBDeviceName,' '];
   TFTBExtension =input('file extension : ','s'); ;
  end;
  if TFTBDevice~=8,
   if comp(1:2)=='PC',
    DefaultName=[method(4:Nmethod),num2str(Nsig),TFTBExtension];
    [name,PathWorkDir] = uiputfile(DefaultName, 'Save As');
   else
    DefaultName=[method(4:Nmethod),num2str(Nsig),TFTBExtension];
    nameStr=[' file name [',DefaultName,'] : '];
    name=input(nameStr,'s');
    while (length(name)>8),
     disp('The name must have less than 8 characters');
     name=input(nameStr,'s');
    end
    if isempty(name),
     name=DefaultName;
    end
   end
   % ['print ', TFTBDeviceName, PathWorkDir, name]
   eval(['print ', TFTBDeviceName, PathWorkDir, name]);
   fprintf(' The file is saved in the directory %s\n',PathWorkDir);
   fprintf('under the name %s\n',name);
  end;
  CallTfrView=0;

 end;

end;

% good bye. I hope that everything happened fine.

author: HeLiangHIT


text:
function imgOut = VanHaterenTMO(img, pupil_area)
%
%
%        imgOut = VanHaterenTMO(img, pupil_area)
%
%
%        Input:
%           -img: input HDR image
%           -pupil_area:
%
%        Output:
%           -imgOut: tone mapped image
%
%     This is the stable version of the Van Hateren 2006 algorithm, this is
%     not suitable for HDR videos.
%
%     Copyright (C) 2010-17  Francesco Banterle
%
%     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/>.
%
%     The paper describing this technique is:
%     "Encoding of High Dynamic Range Video with a Model of Human Cones"
% 	  by J. Hans Van Hateren
%     in ACM Transaction on Graphics 2006
%

check13Color(img);

checkNegative(img);

if(~exist('pupil_area', 'var'))
    pupil_area = -1.0;
end

if(pupil_area <= 0.0)
    pupil_area = 10; %fixed pupil area 10 mm^2
end

k_beta = 1.6e-4; % td/ms
a_C = 9e-2;
C_beta = 2.8e-3; % 1/ms

%Calculate Ios,max
polIosMax = [a_C, 1, 0, 0, 0, -1 / C_beta];
maxIos = max(real(roots(polIosMax)));

%Luminance channel
Lori = lum(img);

%conversion from cd/m^2 to trolands (tr)
L = Lori * pupil_area;

%Range reduction
tmpI = - 1 ./ (C_beta + k_beta * L);
[r, c] = size(L);
n = r * c;
Ios = zeros(size(tmpI));
base = [a_C, 1, 0, 0, 0];

for i = 1:n
    tmp = [base, tmpI(i)];
    Ios(i) = max(real(roots(tmp)));
end

Ld = ClampImg(1 - Ios / maxIos, 0, 1);

%Changing luminance
imgOut = ChangeLuminance(img, Lori, Ld);

warning('The image does not require gamma correction.');

end

author: banterle


text:
function [J, grad] = linearRegCostFunction(X, y, theta, lambda)
%LINEARREGCOSTFUNCTION Compute cost and gradient for regularized linear
%regression with multiple variables
%   [J, grad] = LINEARREGCOSTFUNCTION(X, y, theta, lambda) computes the
%   cost of using theta as the parameter for linear regression to fit the
%   data points in X and y. Returns the cost in J and the gradient in grad

% Initialize some useful values
m = length(y); % number of training examples

% You need to return the following variables correctly
J = 0;
grad = zeros(size(theta));
J = 1/(2*m)*sum((X*theta-y).^2)+lambda/(2*m)*(sum(theta.^2)-theta(1).^2);

grad = 1/m*X'*(X*theta-y)+lambda/m*theta;
grad(1) = grad(1) -lambda/m*theta(1);
% ====================== YOUR CODE HERE ======================
% Instructions: Compute the cost and gradient of regularized linear
%               regression for a particular choice of theta.
%
%               You should set J to the cost and grad to the gradient.
%

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

grad = grad(:);

end

author: Ayatans

repo: Machine-Learning-homework


text:
function tapas_validate_data(data)
%% Check that data complies with the inferface.
%
% aponteeduardo@gmail.com
% copyright (C) 2016
%

if ~isstruct(data)
    error('tapas:validate:data', 'data is not a structure')
end

if ~isfield(data, 'y')
    error('tapas:validate:data', 'data lacks field y');
end

if ~isfield(data, 'u')
    error('tapas:validate:data', 'data lacks field u');
end

end

author: translationalneuromodeling

repo: tapas


text:
function varargout = process_pac_sur_zscore( varargin )
% PROCESS_PAC_ANALYSIS: Further analysis of tpac maps.
%
% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Soheila Samiee, 2018
%   - 1.0: SS. Oct. 2018
%
%
eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'PAC: Z-score with surrogate';
    sProcess.FileTag     = '';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Standardize';%{'Frequency','Time-resolved Phase-Amplitude Coupling'};
    sProcess.Index       = 1020;
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'timefreq'};
    sProcess.OutputTypes = {'timefreq'};
    sProcess.nInputs     = 2;
    sProcess.nMinFiles   = 1;
    sProcess.isPaired    = 1;
    sProcess.isSeparator = 0;

%     % === ANALYSIS TO BE DONE
%     sProcess.options.label.Comment = 'Analysis:';
%     sProcess.options.label.Type    = 'label';
%     sProcess.options.analyze_type.Comment = {'Mean (Over sources)', ...
%         'Median (Over sources)','Z-score on time (If no negative time, on total recording)', ...
%         'Mean (Over time)'};
%     sProcess.options.analyze_type.Type    = 'radio';
%     sProcess.options.analyze_type.Value   = 1;

    % === Label
    sProcess.options.label.Comment = 'Mean on sources: Files A Z-scored with respect to Files B

';
    sProcess.options.label.Type    = 'label';
    % === Using phase
    sProcess.options.usePhase.Comment = 'Use phase in averaging (mean)';
    sProcess.options.usePhase.Type    = 'checkbox';
    sProcess.options.usePhase.Value   = 0;

end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function sOutput = Run(sProcess, sInputsA, sInputsB)  %#ok<DEFNU>
%     sOutput = sInputsA;
%     % Colormap
%     sOutput.ColormapType = 'stat2';
%     % Time-frequency: Change the measure type
%     if strcmpi(sInputsA(1).FileType, 'timefreq')
%         sOutput.Measure = 'other';
%     end

    usePhase = sProcess.options.usePhase.Value;
    tag = '|sur-zscore';

%     % Load first TF file
%     tpacMat = in_bst_timefreq(sInputs(1).FileName, 0);
%     % Error
%     if isempty(tpacMat)
%         bst_report('Error', sProcess, sInput, Messages);
%         return;
%     end

    N = length(sInputsA);
    for iFile = 1:N
        tpac = in_bst_timefreq(sInputsA(iFile).FileName, 0);
        tpac_sur = in_bst_timefreq(sInputsB(iFile).FileName, 0);
        if usePhase
            tmp1 = abs(mean(tpac.sPAC.DynamicPAC.*exp(1i*tpac.sPAC.DynamicPhase),1));
            tmp2 = abs(mean(tpac_sur.sPAC.DynamicPAC.*exp(1i*tpac_sur.sPAC.DynamicPhase),1));
        else
            tmp1 = abs(mean(tpac.sPAC.DynamicPAC,1));
            tmp2 = abs(mean(tpac_sur.sPAC.DynamicPAC,1));
        end

        tpac.sPAC.DynamicPAC  = (tmp1 - mean(tmp2,4))./std(tmp2, [],4);
        if usePhase
            tpac.sPAC.DynamicPhase  = angle(mean(tpac.sPAC.DynamicPAC.*exp(1i*tpac.sPAC.DynamicPhase),1));
        else
            tpac.sPAC.DynamicPhase  = [];
        end

        tpac.TF = tpac.sPAC.DynamicPAC;
        %Saving the files
        tpac.Comment = [tpac.Comment, ' ', tag];
        % Output filename: add file tag
        FileTag = strtrim(strrep(tag, '|', ''));
        pathName = file_fullpath(sInputsA(iFile).FileName);
        sOutput{1} = strrep(pathName, '.mat', ['_' FileTag '.mat']);
        sOutput{1} = file_unique(sOutput{1});
        % Save file
        bst_save(sOutput{1}, tpac, 'v6');
        % Add file to database structure
        db_add_data(sInputsA(iFile).iStudy, sOutput{1}, tpac);
    end

end

author: brainstorm-tools


text:
function planC = ImportDVH(command, planC)
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%global planC;

indexS = planC{end};

if isempty(planC), return; end;
try cd(command);
catch
end;
% read DVH
[FileName,PathName] = uigetfile({'*.*','All files (*.*)';}, sprintf('Loading DVH files'), 'MultiSelect', 'on');
if isequal(FileName,0), return ; end;

scSize = get(0,'ScreenSize');
x = 220; y = 120;
setupDlg = figure('Name','DVH Data Type','NumberTitle','off', ...
    'Position',[(scSize(3)-x)/2, (scSize(4)-y)/2, x, y], 'windowstyle', 'modal');
h = uibuttongroup();
u0 = uicontrol('Style','Radio','String','differential', 'pos',[10 80 100 30],'parent',h,'HandleVisibility','off');
u1 = uicontrol('Style','Radio','String','cumulative', 'pos',[10 40 100 30],'parent',h,'HandleVisibility','off');
u2 = uicontrol('Style','pushbutton','String','Continue', 'pos',[110 10 100 20],'parent',setupDlg,'HandleVisibility','off', ...
    'callback', @continue_callback);
uiwait(setupDlg);
    function continue_callback(source, eventdata)
        so = get(h,'selectedObject');
        if isequal(so, u0)
            setappdata(findobj('tag', 'CERRStartupFig'), 'dvhVolumeType', get(u0, 'string'));
        else
            setappdata(findobj('tag', 'CERRStartupFig'), 'dvhVolumeType', get(u1, 'string'));

        end
        close;
    end

volumeType = getappdata(findobj('tag', 'CERRStartupFig'), 'dvhVolumeType');

if iscell(FileName)
    fileCount = length(FileName);
else
    fileCount = 1;
end;

DVHIndex = length(planC{indexS.DVH}) + 1;

for i = 1:fileCount
    if iscell(FileName)
        dvhFileName = fullfile(PathName, FileName{i});
    else
        dvhFileName = fullfile(PathName, FileName);
    end

    [pathstr, name, ext] = fileparts(dvhFileName);
    name = [name,'.',ext];
    strName = ext(2:end);
    indName = strfind(name, '.');
    doseName = name((indName(end-1)+1):end);
    doseIndex = 1;
    for j=1:length(planC{indexS.dose})
        if isequal(doseName, planC{indexS.dose}(j).fractionGroupID)
            doseIndex = j;
            break;
        end
    end

    try
        [doseBinsV, doseHistV] = readDVHfile(dvhFileName);
    catch
        continue;
    end

    %     doseIndex = stateS.doseSet;
    planC{indexS.DVH}(DVHIndex).volumeType      = volumeType;
    planC{indexS.DVH}(DVHIndex).doseType        = 'plunc';
    planC{indexS.DVH}(DVHIndex).doseUnits        = 'cGy';
    planC{indexS.DVH}(DVHIndex).structureName   = strName;
    planC{indexS.DVH}(DVHIndex).doseIndex       = doseIndex;

    strNum = '';
    for j = 1:length(planC{indexS.structures})
        if  strcmp(upper(strName), upper(planC{indexS.structures}(j).structureName))
            strNum = j;
            break;
        end
    end
    planC{indexS.DVH}(DVHIndex).assocStrUID=planC{indexS.structures}(strNum).strUID;
    planC{indexS.DVH}(DVHIndex).dvhUID = createUID('dvh');
    planC{indexS.DVH}(DVHIndex).assocDoseUID = planC{indexS.dose}(doseIndex).doseUID;

    planC{indexS.DVH}(DVHIndex).fractionGroupID = planC{indexS.dose}(doseIndex).fractionGroupID;

    %Store computational results
    % planC = saveDVHMatrix(DVHIndex, doseBinsV, doseHistV, planC);

    planC{indexS.DVH}(DVHIndex).DVHMatrix(:,1) = doseBinsV;
    planC{indexS.DVH}(DVHIndex).DVHMatrix(:,2) = doseHistV;

    DVHIndex = DVHIndex +1;
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [source] = estimate_fwhm2(source, maxdist)

% ESTIMATE_FWHM2(SOURCE, MAXDIST)
%
% This function computes the Gaussian fwhm of the spatial filters, according to
% least-squares Gaussian fit including data points up until MAXDIST from the
% locations of interest.
%
% This function can only deal with scalar filters.

if nargin<2, maxdist = 2.5; end % maxdist should be in units of the pos in source
if isempty(maxdist), maxdist = inf; end

if islogical(source.inside)
  inside = find(source.inside);
else
  inside = source.inside;
end
ninside = numel(inside);

if ~isfield(source.avg, 'filter')
  ft_error('the input should contain spatial filters in');
end

nchan   = size(source.avg.filter{inside(1)},2);
ndir    = size(source.avg.filter{inside(1)},1);
if ndir~=1,
  ft_error('only scalar filters are allowed as input');
end

%get filters and positions
filter = cat(1,source.avg.filter{inside});
pos    = source.pos(inside,:);

%get the filter correlation matrix
Cmat   = filter*filter';
Cmat   = abs(Cmat)./sqrt(diag(Cmat)*diag(Cmat)');

fwhm    = zeros(size(source.pos,1),1);
onesvec = ones(ninside,1);
for k = 1:ninside
 d   = sqrt(sum( (pos-pos(k*onesvec,:)).^2, 2));
 sel = d<=maxdist;
 s   = gaussfit(Cmat(sel,k)',d(sel)');
 fwhm(inside(k)) = s;
end
source.fwhm = fwhm;

% fwhm    = zeros(size(source.pos,1),3,3);
% onesvec = ones(ninside,1);
% for k = 1:ninside
%   dpos = pos-pos(k*onesvec,:);
%   d   = sqrt(sum(dpos.^2, 2));
%   sel = d<=maxdist;
%   [s,dum] = gaussfit3D(dpos(sel,:), Cmat(sel,k));
%   fwhm(inside(k),:,:) = s;
% end
% source.fwhm = fwhm;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fwhm] = gaussfit(dat, design)

%this function performs a least-squares gaussian fit

%create independent variable
ivar = design.^2;
dat  = log(dat+eps);
beta = dat*ivar'*pinv(ivar*ivar');

sigma = sqrt(-0.5./beta(:,1));
fwhm  = 2.*sqrt(2.*log(2)).*sigma;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [invsigma, R] = gaussfit3D(x, y)

% function for a linear least squares fit of a Gaussian function in 3D
% (with a peak of 1 at (0,0,0)).
%
% use as
%   sigma = gaussfit3D(x, y)
%
% x = Nx3 and y = Nx1, sigma is the covariance describing the gaussian

% create design matrix based in the independent variable x
design = [x(:,1).^2 x(:,2).^2 x(:,3).^2 2.*(x(:,1).*x(:,2)) 2.*(x(:,1).*x(:,3)) 2.*(x(:,2).*x(:,3))];
%design = design-repmat(mean(design,1),[size(design,1) 1]);
design = cat(2, design, ones(size(design,1),1));

% log-transform the dependent variable y
dat = -2.*log(y+eps);

% regression
beta = design\dat;

% residuals
res  = dat - design*beta;
R    = 1-sum(res.^2)./sum(dat.^2);

% create output
invsigma = [beta(1) beta(4) beta(5);beta(4) beta(2) beta(6);beta(5) beta(6) beta(3)];
%sigma = inv(invsigma);

% by construction the exponent to the gaussian looks like this
%
% exp( -0.5.*(-x'*invsigma*x) )
%
% -x'*invsigma*x = [x1 x2 x3]*[a d e [x1
%                              d b f  x2
%                              e f c] x3] =
%
% a*x1^2+b*x2^2+c*x3^2+2d*(x1x2)+2e*(x1x3)+2f*(x2x3)
%
% the variables a through f correspond with the ordered beta weights

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% POWER_BOUNDED computes the power cells about the points (x,y) inside
% the bounding box (must be a rectangle or a square) crs.  If crs is not supplied, an
% axis-aligned box containing (x,y) is used.
% It is optimised to work fast on large number of sites (e.g. 10000 sites or more)
% Input:
%   * x, y: coordinate of the Voronoi point (numPoints x 1)
%   * wts: weights of each point (numPoints x 1)
%   * crs: vortices of the bounding box in cw order (numVert x 2)
% Output:
%   * V: x,y-coordinate of vertices of the power cells
%   * C: indices of the Voronoi cells from V
% See Matlab's voronoin for more information about the output
% Made by: Aaron Becker, atbecker@uh.edu, and Muhammad Kasim, muhammad.kasim@wolfson.ox.ac.uk

function [V,C] = power_bounded(x,y, wts, crs)
    bnd=[min(x) max(x) min(y) max(y)]; %data bounds
    if nargin < 3
        crs=double([bnd(1) bnd(4);bnd(2) bnd(4);bnd(2) bnd(3);bnd(1) bnd(3);bnd(1) bnd(4)]);
    end

    rgx = max(crs(:,1))-min(crs(:,1));
    rgy = max(crs(:,2))-min(crs(:,2));
    rg = max(rgx,rgy);
    midx = (max(crs(:,1))+min(crs(:,1)))/2;
    midy = (max(crs(:,2))+min(crs(:,2)))/2;

    % add 4 additional edges
    xA = [x; midx + [0;0;-5*rg;+5*rg]];
    yA = [y; midy + [-5*rg;+5*rg;0;0]];

    if (all(wts == 0))
        [vi,ci] = voronoin([xA,yA]);
    else
        [vi,ci] = powerDiagram2([xA,yA], [wts;zeros(4,1)]);
    end

    % remove the last 4 cells
    C = ci(1:end-4);
    V = vi;
    % use Polybool to crop the cells
    %Polybool for restriction of polygons to domain.

    maxX = max(crs(:,1)); minX = min(crs(:,1));
    maxY = max(crs(:,2)); minY = min(crs(:,2));
    for ij=1:length(C)
        % thanks to http://www.mathworks.com/matlabcentral/fileexchange/34428-voronoilimit
        Cij = C{ij};
        if (length(Cij) == 0) continue; end;

        % first convert the contour coordinate to clockwise order:
        pts = V(Cij,:);
        K = convhull(pts);
        K = K(end-1:-1:1);
        C{ij} = Cij(K);
        X2 = pts(K,1);
        Y2 = pts(K,2);

        % if all points are inside the bounding box, then skip it
        if (all((X2 <= maxX) & (X2 >= minX) & (Y2 <= maxY) & (Y2 >= minY))) continue; end;

        [xb, yb] = clip_polygons(crs(:,1),crs(:,2),X2,Y2);
        % xb = xb'; yb = yb';
        ix=nan(1,length(xb));
        for il=1:length(xb)
            if any(V(:,1)==xb(il)) && any(V(:,2)==yb(il))
                ix1=find(V(:,1)==xb(il));
                ix2=find(V(:,2)==yb(il));
                for ib=1:length(ix1)
                    if any(ix1(ib)==ix2)
                        ix(il)=ix1(ib);
                    end
                end
                if isnan(ix(il))==1
                    lv=length(V);
                    V(lv+1,1)=xb(il);
                    V(lv+1,2)=yb(il);
                    ix(il)=lv+1;
                end
            else
                lv=length(V);
                V(lv+1,1)=xb(il);
                V(lv+1,2)=yb(il);
                ix(il)=lv+1;
            end
        end
        C{ij} = ix;
    end
end

author: optimaltransport


text:
function value = r8_cinh ( x )

%*****************************************************************************80
%
%% R8_CINH: alternate hyperbolic cosine integral Cinh of an R8 argument.
%
%  Discussion:
%
%    Cinh ( x ) = Integral ( 0 <= t <= x ) ( cosh ( t ) - 1 ) dt / t
%
%    The original text of this program had a mistake:
%      y = x * x / 9.0 - 1.0
%    has been corrected to
%      y = x * x / 4.5 - 1.0
%    JVB, 27 March 2010
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    28 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the hyperbolic cosine integral Cinh
%    evaluated at X.
%
  persistent cinhcs
  persistent eul
  persistent ncinh
  persistent xmin
  persistent xsml

  eul = 0.57721566490153286060651209008240;

  if ( isempty ( ncinh ) )

    cinhcs = [ ...
      0.1093291636520734431407425199795917, ...
      0.0573928847550379676445323429825108, ...
      0.0028095756978830353416404208940774, ...
      0.0000828780840721356655731765069792, ...
      0.0000016278596173914185577726018815, ...
      0.0000000227809519255856619859083591, ...
      0.0000000002384484842463059257284002, ...
      0.0000000000019360829780781957471028, ...
      0.0000000000000125453698328172559683, ...
      0.0000000000000000663637449497262300, ...
      0.0000000000000000002919639263594744, ...
      0.0000000000000000000010849123956107, ...
      0.0000000000000000000000034499080805, ...
      0.0000000000000000000000000094936664, ...
      0.0000000000000000000000000000228291, ...
      0.0000000000000000000000000000000484 ]';

    ncinh = r8_inits ( cinhcs, 16, 0.1 * r8_mach ( 3 ) );
    xsml = sqrt ( r8_mach ( 3 ) );
    xmin = 2.0 * sqrt ( r8_mach ( 1 ) );

  end

  absx = abs ( x );

  if ( x == 0.0 )
    value = 0.0;
  elseif ( absx <= xmin )
    value = 0.0;
  elseif ( x <= xsml )
    y = - 1.0;
    value = x * x * ( 0.25 + r8_csevl ( y, cinhcs, ncinh ) );
  elseif ( x <= 3.0 )
    y = x * x / 4.5 - 1.0;
    value = x * x * ( 0.25 + r8_csevl ( y, cinhcs, ncinh ) );
  else
    value = r8_chi ( absx ) - eul - log ( absx );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function sump=sum(p,dim)
% sympoly/sum: sum a sympoly array along a given dimension
% usage: sump=sum(p);
% usage: sump=sum(p,dim);
%
% arguments:
%    p - sympoly array object
%  dim - (OPTIONAL) dimension to sum over
%        DEFAULT ==1, unless sp is a row vector
%
%  sump - sympoly object containing the sum reduced array

s = size(p);
np = length(s);

% default for dim is 1, UNLESS p is a row vector.
if (nargin<2) || isempty(dim)
  if (s(1) == 1) && (np==2)
    % a row vector
    dim = 2;
  else
    % any other shape array
    dim = 1;
  end
end

% for dim == 1 or dim == 2, do the sum using a dot product
if (dim == 1) && (np == 2)
  % sum down rows
  sump = ones(1,s(1))*p;

elseif (dim == 2) && (np == 2)
  % sum across columns
  sump = p*ones(s(2),1);

else
  % its an n-d array
  ss = s;
  ss(dim) = 1;

  si = cell(1,np);
  for i = 1:np
    si{i} = 1:s(i);
  end

  if any(ss~=1)
    sump = repmat(sympoly(0),ss);
  else
    sump = sympoly(0);
  end
  for i = 1:s(dim)
    si{dim} = i;
    sump = sump + p(si{:});
  end

end

author: Sable

repo: mcbench-benchmarks


text:
%PNMTEST Test program for the PNM Toolbox.
%   PNMTEST is a test program for the PNM Toolbox. It is not really
%   intended for end users. It checks the most important read/write
%   functionality to see if the PNM Toolbox behaves as expected.
%
%   See also PNMDEMO.

%   Author:      Peter J. Acklam
%   Time-stamp:  1998-05-11 17:15:29
%   E-mail:      jacklam@math.uio.no (Internet)
%   URL:         http://www.math.uio.no/~jacklam

more off;
echo on;
home; clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PPM test
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% Read PPM image as an RGB image.
%

[ r, g, b ] = ppmread( 'logo5' );

%pause

%
% Simplify color map and display image to see if it appears to be OK.
%

[ x, map ] = im2ind( r, g, b );

figure( gcf ); clf;
colormap( map ); image( x ); axis image;

pause

%
% Write an ascii and a binary encoded version from indexed data.
%

ppmwrite( x, map, 'logo5_asc', 'ascii' );
ppmwrite( x, map, 'logo5_bin', 'binary' );

%pause

%
% Read them as indexed images and compare them.
%

[ x1, map1 ] = ppmread( 'logo5_asc' );
[ x2, map2 ] = ppmread( 'logo5_bin' );

if any( any( map1(x1(:),:) ) ~= any( map2(x2(:),:) ) )
   error( 'Test failed.' );
end

%pause

%
% Read them as RGB images and compare them.
%

[ r1, g1, b1 ] = ppmread( 'logo5_asc' );
[ r2, g2, b2 ] = ppmread( 'logo5_bin' );

if any( r1(:) ~= r2(:) ) | any( g1(:) ~= g2(:) ) | any( b1(:) ~= b2(:) )
   error( 'Test failed.' );
end

%pause

%
% The indexed and RGB images should be identical.
%

if any( any( map1(x1(:),:) ~= [ r(:) g(:) b(:) ] ) )
   error( 'Test failed.' );
end

%pause

%
% Write an ascii and a binary encoded version from RGB data.
%

ppmwrite( r, g, b, 'logo5_asc', 'ascii' );
ppmwrite( r, g, b, 'logo5_bin', 'binary' );

%pause

%
% Read them as indexed images and compare them.
%

[ x1, map1 ] = ppmread( 'logo5_asc' );
[ x2, map2 ] = ppmread( 'logo5_bin' );

if any( any( map1(x1(:),:) ~= map2(x2(:),:) ) )
   error( 'Test failed.' );
end

%pause

%
% Read them as RGB images and compare them.
%

[ r1, g1, b1 ] = ppmread( 'logo5_asc' );
[ r2, g2, b2 ] = ppmread( 'logo5_bin' );

if any( r1(:) ~= r2(:) ) | any( g1(:) ~= g2(:) ) | any( b1(:) ~= b2(:) )
   error( 'Test failed.' );
end

%pause

%
% The indexed and RGB images should be identical.
%

if any( any( map1(x1(:),:) ~= [ r(:) g(:) b(:) ] ) )
   error( 'Test failed.' );
end

%pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PGM test
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pgmwrite( x, map, 'logo5', 'binary' );
i = pgmread( 'logo5' );

%pause

%
% Simplify color map and display image to see if it appears to be OK.
%

[ x, map ] = im2ind( i );

figure( gcf ); clf;
colormap( map ); image( x ); axis image;

pause

%
% Write an ascii and a binary encoded version from indexed data.
%

pgmwrite( x, map, 'logo5_asc', 'ascii' );
pgmwrite( x, map, 'logo5_bin', 'binary' );

%pause

%
% Read them as indexed images and compare them.
%

[ x1, map1 ] = pgmread( 'logo5_asc' );
[ x2, map2 ] = pgmread( 'logo5_bin' );

if any( any( map1(x1(:),:) ~= map2(x2(:),:) ) )
   error( 'Test failed.' );
end

%pause

%
% Read them as intensity images and compare them.
%

i1 = pgmread( 'logo5_asc' );
i2 = pgmread( 'logo5_bin' );

if any( i1(:) ~= i2(:) )
   error( 'Test failed.' );
end

%pause

%
% The indexed and intensity images should be identical.
%

if any( any( map1(x1(:),:) ~= [ i(:) i(:) i(:) ] ) )
   error( 'Test failed.' );
end

%pause

%
% Write an ascii and a binary encoded version from intensity data.
%

pgmwrite( i, 'logo5_asc', 'ascii' );
pgmwrite( i, 'logo5_bin', 'binary' );

%pause

%
% Read them as indexed images and compare them.
%

[ x1, map1 ] = pgmread( 'logo5_asc' );
[ x2, map2 ] = pgmread( 'logo5_bin' );

if any( any( map1(x1(:),:) ~= map2(x2(:),:) ) )
   error( 'Test failed.' );
end

%pause

%
% Read them as intensity images and compare them.
%

i1 = pgmread( 'logo5_asc' );
i2 = pgmread( 'logo5_bin' );

if any( i1(:) ~= i2(:) )
   error( 'Test failed.' );
end

%pause

%
% The indexed and intensity images should be identical.
%

if any( any( map1(x1(:),:) ~= [ i(:) i(:) i(:) ] ) )
   error( 'Test failed.' );
end

%pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PBM test
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

pbmwrite( x, map, 'logo5', 'binary' );
i = pbmread( 'logo5' );

%pause

%
% Simplify color map and display image to see if it appears to be OK.
%

[ x, map ] = im2ind( i );

figure( gcf ); clf;
colormap( map ); image( x ); axis image;

pause

%
% Write an ascii and a binary encoded version from indexed data.
%

pbmwrite( x, map, 'logo5_asc', 'ascii' );
pbmwrite( x, map, 'logo5_bin', 'binary' );

%pause

%
% Read them as indexed images and compare them.
%

[ x1, map1 ] = pbmread( 'logo5_asc' );
[ x2, map2 ] = pbmread( 'logo5_bin' );

if any( any( map1(x1(:),:) ~= map2(x2(:),:) ) )
   error( 'Test failed.' );
end

%pause

%
% Read them as intensity images and compare them.
%

i1 = pbmread( 'logo5_asc' );
i2 = pbmread( 'logo5_bin' );

if any( i1(:) ~= i2(:) )
   error( 'Test failed.' );
end

%pause

%
% The indexed and intensity images should be identical.
%

if any( any( map1(x1(:),:) ~= [ i(:) i(:) i(:) ] ) )
   error( 'Test failed.' );
end

%pause

%
% Write an ascii and a binary encoded version from intensity data.
%

pbmwrite( i, 'logo5_asc', 'ascii' );
pbmwrite( i, 'logo5_bin', 'binary' );

%pause

%
% Read them as indexed images and compare them.
%

[ x1, map1 ] = pbmread( 'logo5_asc' );
[ x2, map2 ] = pbmread( 'logo5_bin' );

if any( any( map1(x1(:),:) ~= map2(x2(:),:) ) )
   error( 'Test failed.' );
end

%pause

%
% Read them as intensity images and compare them.
%

i1 = pbmread( 'logo5_asc' );
i2 = pbmread( 'logo5_bin' );

if any( i1(:) ~= i2(:) )
   error( 'Test failed.' );
end

%pause

%
% The indexed and intensity images should be identical.
%

if any( any( map1(x1(:),:) ~= [ i(:) i(:) i(:) ] ) )
   error( 'Test failed.' );
end

%pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Miscellaneous testing.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
% PPMREAD should be able to read a PGM image as an indexed image.
%

[ x1, map1 ] = ppmread( 'logo5_bin.pgm' );
[ x2, map2 ] = pgmread( 'logo5_bin.pgm' );

if any( x1(:) ~= x2(:) ) | any( map1(:) ~= map2(:) )
   error( 'Test failed.' );
end

%
% PPMREAD should be able to read a PGM image as RGB image.
%

[ r, g, b ] = ppmread( 'logo5_bin.pgm' );
i           = pgmread( 'logo5_bin.pgm' );

if any( r(:) ~= i(:) ) | any( g(:) ~= i(:) ) | any( b(:) ~= i(:) )
   error( 'Test failed.' );
end

%
% PPMREAD should be able to read a PBM image as indexed image.
%

[ x1, map1 ] = ppmread( 'logo5_bin.pbm' );
[ x2, map2 ] = pbmread( 'logo5_bin.pbm' );

if any( x1(:) ~= x2(:) ) | any( map1(:) ~= map2(:) )
   error( 'Test failed.' );
end

%
% PPMREAD should be able to read a PBM image as an RGB image.
%

[ r, g, b ] = ppmread( 'logo5_bin.pbm' );
i           = pbmread( 'logo5_bin.pbm' );

if any( r(:) ~= i(:) ) | any( g(:) ~= i(:) ) | any( b(:) ~= i(:) )
   error( 'Test failed.' );
end

%
% PGMREAD should be able to read a PBM image as an indexed image.
%

[ x1, map1 ] = pgmread( 'logo5_bin.pbm' );
[ x2, map2 ] = pbmread( 'logo5_bin.pbm' );

if any( x1(:) ~= x2(:) ) | any( map1(:) ~= map2(:) )
   error( 'Test failed.' );
end

%
% PGMREAD should be able to read PBM a image as an intensity image.
%

i1 = pgmread( 'logo5_bin.pbm' );
i2 = pbmread( 'logo5_bin.pbm' );

if any( i1(:) ~= i2(:) )
   error( 'Test failed.' );
end

%
% That should be all for now, so clean up.
%

delete logo5.pgm logo5.pbm
delete logo5_asc.ppm logo5_asc.pgm logo5_asc.pbm
delete logo5_bin.ppm logo5_bin.pgm logo5_bin.pbm

echo off

fprintf( '\nPassed all tests.\n\n' );


text:
function len = s_len_trim ( s )

%*****************************************************************************80
%
%% S_LEN_TRIM returns the length of a string to the last nonblank.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 June 2003
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, string S, the string to be measured.
%
%    Output, integer LEN, the length of the string up to the last nonblank.
%
  len = length ( s );

  while ( 0 < len )
    if ( s(len) ~= ' ' )
      return
    end
    len = len - 1;
  end

  return
end


text:
% This is material illustrating the methods from the book
% Financial Modelling  - Theory, Implementation and Practice with Matlab
% source
% Wiley Finance Series
% ISBN 978-0-470-74489-5
%
% Date: 02.05.2012
%
% Authors:  Joerg Kienitz
%           Daniel Wetterau
%
% Please send comments, suggestions, bugs, code etc. to
% kienitzwetterau_FinModelling@gmx.de
%
% (C) Joerg Kienitz, Daniel Wetterau
%
% Since this piece of code is distributed via the mathworks file-exchange
% it is covered by the BSD license
%
% This code is being provided solely for information and general
% illustrative purposes. The authors will not be responsible for the
% consequences of reliance upon using the code or for numbers produced
% from using the code.

function result = cvalue(x1, x2, a, b, N, V, ...
                        model, t, r, q, varargin)

Nstrike = length(x1);

exp2 = exp( 1i .* repmat((1:N)',1,Nstrike) * diag((x2 - a) ./ (b - a)) .* pi );    % init
exp1 = exp( 1i .* repmat((1:N)',1,Nstrike) * diag((x1 - a) ./ (b - a)) .* pi );    % init

m = zeros(3*N-1, Nstrike);                                        % init base

m(N,:) = 1i * pi * (x2 - x1) ./ (b - a);
m(N+1:2*N,:) = 1 ./ repmat((1:N)',1,Nstrike) .* ( exp2 - exp1 );
m(1:N-1,:) = - conj(flipud(m(N+1:2*N-1, :)));
m(2*N+1:3*N-1,:) = ( exp2(1:N-1,:) * diag(exp2(N, :)) - exp1(1:N-1,:) ...
    * diag(exp1(N,:)) ) ./ ( repmat((N+1:2*N-1)',1,Nstrike) );

Grid_j = (0:N-1)';                                          % fix grid

% compute u values
u = exp(feval(@CF, model,pi*repmat(Grid_j,1,Nstrike)*diag(1./(b-a)), t,r,q,varargin{:})) .* V;
u(1,:) = 0.5*u(1,:);

m_s = [m(N:-1:1, :); zeros(1,Nstrike); m(2*N-1:-1:N+1, :)];
u_s = [u; zeros(N, Nstrike)];
m_c = m(3*N-1:-1:N, :);

shortCut = 1;

% apply fft five times
if shortCut == 1
    zeta = -ones(2*N, Nstrike);
    zeta(2 .* (1:N)' - 1,:) = 1;

    fft_u_s = fft(u_s);
    xi_s = ifft((fft(m_s)) .* fft_u_s);
    xi_c = ifft((fft(m_c)) .* (zeta .* fft_u_s));

    result = exp(-r * t) / pi .* imag( xi_s(1:N,:) + flipud(xi_c(1:N,:)) );
else
    M_c = zeros(N, N);
    M_s = zeros(N, N);

    result = zeros(N,Nstrike);

    for k = 1:Nstrike
        for n = 0:N-1
            M_c(:, n+1) = m(N+n:2*N-1+n,k);
            M_s(:, n+1) = m(N+n:-1:1+n,k);
        end

        result(:,k) = exp(-r*t) / pi .* imag((M_c + M_s) * u(:,k));
    end

end


text:
function [ n_data, n, x, fx ] = lobatto_polynomial_values ( n_data )

%*****************************************************************************80
%
%% LOBATTO_POLYNOMIAL_VALUES returns values of the completed Lobatto polynomials.
%
%  Discussion:
%
%    In Mathematica, the function can be evaluated by:
%
%      n * LegendreP [ n - 1, x ] - n * x * LegendreP [ n, x ]
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 May 2013
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Milton Abramowitz and Irene Stegun,
%    Handbook of Mathematical Functions,
%    US Department of Commerce, 1964.
%
%    Stephen Wolfram,
%    The Mathematica Book,
%    Fourth Edition,
%    Wolfram Media / Cambridge University Press, 1999.
%
%  Parameters:
%
%    Input/output, integer N_DATA.  The user sets N_DATA to 0 before the
%    first call.  On each call, the routine increments N_DATA by 1, and
%    returns the corresponding data; when there is no more data, the
%    output value of N_DATA will be 0 again.
%
%    Output, integer N, the order of the function.
%
%    Output, real X, the point where the function is evaluated.
%
%    Output, real FX, the value of the function.
%
  n_max = 31;

  fx_vec = [ ...
    0.9375000000000000, ...
    0.7031250000000000, ...
   -0.9667968750000000, ...
   -1.501464843750000, ...
    0.3639221191406250, ...
    2.001914978027344, ...
    0.6597948074340820, ...
   -1.934441328048706, ...
   -1.769941113889217, ...
    1.215243665501475, ...
    0.000000000000000, ...
    0.8692500000000000, ...
    1.188000000000000, ...
    1.109250000000000, ...
    0.7680000000000000, ...
    0.2812500000000000, ...
   -0.2520000000000000, ...
   -0.7507500000000000, ...
   -1.152000000000000, ...
   -1.410750000000000, ...
   -1.500000000000000, ...
   -1.410750000000000, ...
   -1.152000000000000, ...
   -0.7507500000000000, ...
   -0.2520000000000000, ...
    0.2812500000000000, ...
    0.7680000000000000, ...
    1.109250000000000, ...
    1.188000000000000, ...
    0.8692500000000000, ...
    0.000000000000000 ];

  n_vec = [ ...
     1,  2, ...
     3,  4,  5, ...
     6,  7,  8, ...
     9, 10,  3, ...
     3,  3,  3, ...
     3,  3,  3, ...
     3,  3,  3, ...
     3,  3,  3, ...
     3,  3,  3, ...
     3,  3,  3, ...
     3,  3 ];

  x_vec = [ ...
    0.25, ...
    0.25, ...
    0.25, ...
    0.25, ...
    0.25, ...
    0.25, ...
    0.25, ...
    0.25, ...
    0.25, ...
    0.25, ...
   -1.00, ...
   -0.90, ...
   -0.80, ...
   -0.70, ...
   -0.60, ...
   -0.50, ...
   -0.40, ...
   -0.30, ...
   -0.20, ...
   -0.10, ...
    0.00, ...
    0.10, ...
    0.20, ...
    0.30, ...
    0.40, ...
    0.50, ...
    0.60, ...
    0.70, ...
    0.80, ...
    0.90, ...
    1.00 ];

  if ( n_data < 0 )
    n_data = 0;
  end

  n_data = n_data + 1;

  if ( n_max < n_data )
    n_data = 0;
    n = 0;
    x = 0.0;
    fx = 0.0;
  else
    n = n_vec(n_data);
    x = x_vec(n_data);
    fx = fx_vec(n_data);
  end

  return
end


text:
function y = in2cm(x)
%IN2CM converts inches to centimeters
%
%   Usage: y = in2cm(x)
%
%   Input parameters:
%       x     - input / inches
%
%   Output parameters:
%       y     - output / cm
%
%   See also: in2px, cm2px, cm2in, px2in, px2cm

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
% IMPLIED, INCLUDING BUT  NOT LIMITED TO THE  WARRANTIES OF MERCHANTABILITY, *
% FITNESS  FOR A PARTICULAR  PURPOSE AND  NONINFRINGEMENT. IN NO EVENT SHALL *
% THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
% LIABILITY, WHETHER  IN AN  ACTION OF CONTRACT, TORT  OR OTHERWISE, ARISING *
% FROM,  OUT OF  OR IN  CONNECTION  WITH THE  SOFTWARE OR  THE USE  OR OTHER *
% DEALINGS IN THE SOFTWARE.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input parameters ====================================
nargmin = 1;
nargmax = 1;
narginchk(nargmin,nargmax);

%% ===== Computation =====================================================
y = x * 2.54;

author: sfstoolbox

repo: sfs-matlab


text:
% --------------------------
% XML TOOLBOX FOR MATLAB 3.1
% --------------------------
%
% FUNCTIONS:
%  xml_format       converts a Matlab variable into XML string
%  xml_formatany    converts a Matlab variable into XML string with attributes
%  xml_parse        parses & converts XML string into Matlab structure
%  xml_parseany     parses & converts XML string into Matlab structure with attributes
%  xml_save         saves a Matlab variable/structure in XML format in a file
%  xml_load         loads the .xml file written with xml_save back into a variable
%  xml_help         this file, displays info about available xml_* commands
%  tests/xml_tests  tests the xml toolbox by writing/reading a number of xml test files
%  strsplit         utility function which splits a string at specified characters
%  unblank          utility function which removes leading and trailing spaces
%
% FILES:
%  doc/xml_matlab   documentation containing info on installation, usage, implementation, etc.
%  matlab.xsd       contains the Schema to validate xml files for the toolbox
%                   (if not present, look at http://www.geodise.org/matlab.xsd)
%
% RELATED:
%  xmlread, xmlwrite (shipped with Matlab from version 6.5)
%
% Further information can be obtained by using the help command on
% a specific function, e.g. help xml_format or by visiting
% http://www.geodise.org

% Copyright (C) 2002-2005, University of Southampton
% Author: Dr Marc Molinari <m.molinari@soton.ac.uk>
% $Revision: 1.1 $ $Date: 2005/04/15 17:12:14 $ $Tag$

author: Cloud-CV

repo: object-proposals


text:
function y = wshift(type,x,p)
%WSHIFT Shift Vector or Matrix.
%   Y = WSHIFT(TYPE,X,P) with TYPE = {1,'1','1d' or '1D'}
%   performs a P-circular shift of vector X.
%   The shift P must be an integer, positive for right to left
%   shift and negative for left to right shift.
%
%   Y = WSHIFT(TYPE,X,P) with TYPE = {2,'2','2d' or '2D'}
%   performs a P-circular shift of matrix X.
%   The shifts P must be integers. P(1) is the shift for rows
%   and P(2) is the shift for columns.
%
%   WSHIFT('1D',X) is equivalent to WSHIFT('1D',X,1)
%   WSHIFT('2D',X) is equivalent to WSHIFT('2D',X,[1 1])
%
%   Example 1D:
%     x = [1 2 3 4 5]
%     wshift('1D',x,1)  = [2 3 4 5 1]
%     wshift('1D',x,-1) = [5 1 2 3 4]
%
%   Example 2D:
%     x = [1 2 3;5 6 7]
%     wshift('2D',x,[1 1])  = [6 7 5;2 3 1]
%     wshift('2D',x,[-1,0]) = [5 6 7;1 2 3]

%   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 01-Dec-97.
%   Last Revision: 01-May-1998.
%   Copyright (c) 1995-98 by The MathWorks, Inc.
%   $Revision: 1.2 $  $Date: 1998/07/16 14:33:59 $

if isempty(x) | all(p==0) , y = x; return; end
switch type
  case {1,'1','1d','1D'}
    if nargin<3 , p = 1; end
    L = length(x);
    p = rem(p,L);
    if p<0 , p = L+p; end
    y = x([p+1:L,1:p]);

  case {2,'2','2d','2D'}
    if nargin<3 , p = [1 1]; end
    S = size(x);
    p = rem(p,S);
    k = (p<0);
    p(k) = S(k)+p(k);
    y = x([p(1)+1:S(1),1:p(1)],[p(2)+1:S(2),1:p(2)]);
end


text:
function [vw,newMeshNum] = meshBuild(vw,hemisphere)
% Build a 3D mesh for visualization and analysis
%
%   [vw,newMeshNum] = meshBuild(vw,[hemisphere]);
%
%  Using mrVista data, build a mesh, save it in a file in the anatomy
%  directory, and add the mesh to the 3D Control Window pull down
%  options.
%
% vw:          A VISTASOFT view structure
% hemisphere:  left, right or both.  Default: left
%
% A mrMesh window is opened as well, showing the computed mesh.
%
% Example:
%    [VOLUME{1},newMeshNum] = meshBuild(VOLUME{1},'left');
%    VOLUME{1} = viewSet(VOLUME{1},'currentmeshn',newMeshNum);
%
% See also:  meshBuildFromClass, meshBuildFromNiftiClass, meshSmooth,
%            meshColor
%
% 11/05 ras: also saves mesh path.
%
% (c) Stanford VISTA Team 2008

% Programming TODO.  Check this!
%   We have (or had?) trouble for building 'both' meshes.  We need a new
%   procedure.

% Be sure anatomy is loaded (we need it for the mmPerVox field)
if isempty(vw.anat), vw = loadAnat(vw); end
if ieNotDefined('hemisphere'), hemisphere = 'left'; end

newMeshNum = viewGet(vw,'nmesh') + 1;

% Parameters we establish in this routine
[meshName,numGrayLayers,hemiNum] = readParams(newMeshNum,hemisphere);
if isempty(meshName), newMeshNum = newMeshNum - 1; return; end   % User pressed cancel.

% mmPerVox = viewGet(vw,'mmPerVoxel');

wbar = mrvWaitbar(0.1, ...
    sprintf('meshBuild: Combining white and gray matter...'));

% Load left, right, or both hemispheres.
if (hemiNum==1)
    [voxels,vw] = meshGetWhite(vw, 'left', numGrayLayers);
elseif (hemiNum==2)
    [voxels,vw] = meshGetWhite(vw, 'right', numGrayLayers);
elseif (hemiNum == 0)
    [voxels,vw] = meshGetWhite(vw, 'left', numGrayLayers);
    [voxels,vw] = meshGetWhite(vw, 'right', numGrayLayers,voxels);
end

% host = 'localhost';
% windowID = -1;

% We build a smoothed (mesh) and an unsmoothed mesh (tenseMesh) with these calls
mrvWaitbar(0.35,wbar,sprintf('Building mesh'));
[newMesh, tenseMesh] = mrmBuild(voxels,viewGet(vw,'mmPerVox'),1);

% Must have a name
newMesh   = meshSet(newMesh,'name',meshName);
tenseMesh = meshSet(tenseMesh,'name',sprintf('%s-tense',meshName));

% mrvWaitbar(0.65,wbar,sprintf('meshBuild: Unsmoothed mesh vertex to gray mapping'));
initVertices = meshGet(tenseMesh,'vertices');
newMesh = meshSet(newMesh,'initialvertices',initVertices);
vertexGrayMap = mrmMapVerticesToGray(...
    initVertices, ...
    viewGet(vw,'nodes'), ...
    viewGet(vw,'mmPerVox'),...
    viewGet(vw,'edges'));

newMesh = meshSet(newMesh,'vertexGrayMap',vertexGrayMap);
newMesh = meshSet(newMesh,'name',meshName);
newMesh = meshSet(newMesh,'nGrayLayers',numGrayLayers);

mrvWaitbar(0.9,wbar,sprintf('meshBuild: Saving mesh file %s',meshGet(newMesh,'name')));

% Save mesh file
[newMesh newMesh.path] = mrmWriteMeshFile(newMesh);

mrvWaitbar(1,wbar,sprintf('meshBuild: Done'));
pause(0.5);
close(wbar);

% Now refresh the UI
vw = viewSet(vw,'add and select mesh',newMesh);

return;

%---------------------------------------
function classFile = verifyClassFile(vw,hemisphere)

classFile =  viewGet(vw,'classFileName',hemisphere);
str = sprintf('Class %s',classFile);

r=questdlg(str);
if ~strcmp(r,'Yes')
    switch hemisphere
        case 'left'
            vw = viewSet(vw,'leftClassFileName',[]);
        case 'right'
            vw = viewSet(vw,'rightClassFileName',[]);
    end
    classFile =  viewGet(vw,'classFileName',hemisphere);
end

return;

%---------------------------------------
function voxels = classExtractWhite(voxels,data,voi,whiteValue)
%

% ras 05/07: the indexing of data seems off to me -- is this correct?
voxels(voi(1):voi(2), voi(3):voi(4), voi(5):voi(6)) = ...
    voxels(voi(1):voi(2), voi(3):voi(4), voi(5):voi(6)) ...
    | (data(voi(1):voi(2), voi(3):voi(4), voi(5):voi(6)) == whiteValue);

return;

%----------------------------------------
function [meshName,numGrayLayers,hemiNum,alpha,restrictVOI,relaxIterations] = ...
    readParams(newMeshNum,hemisphere)
%
%  readParams
%
% Internal routine to read the parameters for meshBuild
%
meshName = sprintf('%sSmooth',hemisphere);
numGrayLayers = 0;
switch hemisphere
    case 'left'
        hemiNum = 1;
    case 'right'
        hemiNum = 2;
    case 'both'
        hemiNum = 0;
end

% transparency level (transparency is off by default, but if it gets turned
% on, this alpha parameter will have an effect).
alpha = 200;
restrictVOI = 1;
relaxIterations = 0.2;

prompt = {'Mesh Name:',...
        'Number of Gray Layers (0-4):',...
        'Hemisphere (0=both, 1=left, 2=right):',...
        % 'Default alpha (0-255):',...
        % 'Inflation (0=none, 1=lots):',...
        % 'Restrict to class VOI (0|1):'};
        };
defAns = {meshName,...
        num2str(numGrayLayers),...
        num2str(hemiNum),...
        % num2str(alpha),...
        % num2str(relaxIterations),...
        % num2str(restrictVOI)};
        };

resp = inputdlg(prompt, 'meshBuild Parameters', 1, defAns);

if(~isempty(resp))
    meshName = resp{1};
    numGrayLayers = str2num(resp{2});
    hemiNum = str2num(resp{3});
    % alpha = str2num(resp{4});
    % relaxIterations = round(str2num(resp{5})*160);  % Arbitrary choice, scales iters [0,160]
    % restrictVOI = str2num(resp{6});
else
    meshName = [];
    numGrayLayers = [];
    hemiNum = [];
    % alpha = [];
    % relaxIterations = [];  % Arbitrary choice, scales iters [0,160]
    % restrictVOI = [];
end

return;

%---------------------------------
function [voxels,vw] = meshGetWhite(vw, hemiName, numGrayLayers, voxels)
%
%
%

if ieNotDefined('vw'), error('You must send in a volume vw'); end
if ieNotDefined('hemiName'), error('You must define right,left or both'); end
if ieNotDefined('numGrayLayers'), numGrayLayers = 0; end

classFile = verifyClassFile(vw,hemiName);
if isempty(classFile),
    close(wbar); newMeshNum = -1;
    voxels = [];
    return;
end
classFileParam = [hemiName,'ClassFile'];
vw       = viewSet(vw,classFileParam,classFile);

classData = viewGet(vw,'classdata',hemiName);
if ieNotDefined('voxels'),
    voxelsOld = uint8(zeros(classData.header.xsize, ...
        classData.header.ysize, ...
        classData.header.zsize));
else
    voxelsOld = voxels;
end
voxels = zeros(classData.header.xsize, ...
        classData.header.ysize, ...
        classData.header.zsize);

% Restrict the white matter volume to a size equal to the ROI in which it
% was selected
voxels = classExtractWhite(voxels,...
    classData.data,classData.header.voi,classData.type.white);
%   msh = meshColor(meshSmooth(meshBuildFromClass(voxels,[1 1 1])));
%   meshVisualize(msh);

% Add the gray matter
if(numGrayLayers>0)

    [nodes,edges,classData] = mrgGrowGray(classData,numGrayLayers);
    voxels = ...
        uint8( (classData.data == classData.type.white) | ...
        (classData.data == classData.type.gray));
end

voxels = uint8(voxels | voxelsOld);

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function out = SC_fastdfa(y)
% SC_fastdfa   Matlab wrapper for Max Little's ML_fastdfa code
%
% Measures the scaling exponent of the time series using a fast implementation
% of detrended fluctuation analysis (DFA).
%
%---INPUT:
% y, the input time series, is fed straight into the fastdfa script.

% The original fastdfa code is by Max A. Little and publicly-available at
% http://www.maxlittle.net/software/index.php
% ------------------------------------------------------------------------------
% Copyright (C) 2020, Ben D. Fulcher <ben.d.fulcher@gmail.com>,
% <http://www.benfulcher.com>
%
% If you use this code for your research, please cite the following two papers:
%
% (1) B.D. Fulcher and N.S. Jones, "hctsa: A Computational Framework for Automated
% Time-Series Phenotyping Using Massive Feature Extraction, Cell Systems 5: 527 (2017).
% DOI: 10.1016/j.cels.2017.10.001
%
% (2) B.D. Fulcher, M.A. Little, N.S. Jones, "Highly comparative time-series
% analysis: the empirical structure of time series and their methods",
% J. Roy. Soc. Interface 10(83) 20130048 (2013).
% DOI: 10.1098/rsif.2013.0048
%
% This function 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/>.
% ------------------------------------------------------------------------------

if size(y,2) > size(y,1);
    y = y'; % Ensure input time series is a column vector
end

out = ML_fastdfa(y);

end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
%% FUNCTION Least_Weight2FGLasso
%   Weighted Fused Group Lasso and Weighted Lasso with Least Squares Loss .
%
%% OBJECTIVE
%   argmin_W { \sum_i^t (0.5 * norm (Y{i} - X{i}' * W(:, i))^2)
%              + \sum_{i=1}^d rho1{i} * \|W(i, :)\|_1
%              + \sum_{i=1}^d rho2{i} \sum_{j=1}^{t-1}* \|W(i, j) - W(i, j+1)\|_1
%          }
%   rho1: weighted Lasso sparse.
%   rho2: weighted Fused Lasso.
%   R encodes fused structure relationship [1 -1 0 ...; 0 1 -1 ...; ...]
%      R=zeros(t,t-1);R(1:(t+1):end)=1;R(2:(t+1):end)=-1;
%
%% INPUT
% X: {d * n} * t - input matrix
% Y: {n * 1} * t - output matrix
% rho1: d * 1    - weighted Lasso.
% rho2: d * 1    - weighted fused Lasso.
%
%% OUTPUT
% W: model: d * t
% funcVal: function value vector.
%
%% LICENSE
%   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/>.
%
%   Copyright (C) 2011 - 2012 Jiayu Zhou, Jun Liu and Jieping Ye
%
%   You are suggested to first read the Manual.
%   For any problem, please contact with Jiayu Zhou via jiayu.zhou@asu.edu
%
%   Last modified on June 3, 2012.
%
%% RELATED PAPERS
%
%   [1] Zhou, J., Jun, L., Narayan, A. V. and Ye, J.  Modeling Disease
%   Progression via Fused Sparse Group Lasso. KDD 2012
%
%% RELATED FUNCTIONS
%   Least_NCFGLasso, init_opts

function [W, funcVal] = Least_Weight2FGLasso(X, Y, rho1, rho2, opts)

if nargin <5
    opts = [];
end

% initialize options.
opts=init_opts(opts);

task_num  = length (X);
dimension = size(X{1}, 1);

if length(rho1) ~= dimension
    error('Size of rho1 is not correct! Should be equivalent to dimension')
end

if length(rho2) ~= dimension
    error('Size of rho2 is not correct! Should be equivalent to dimension')
end

funcVal = [];

% Relation
R=zeros(task_num,task_num-1);
R(1:(task_num+1):end)=1;
R(2:(task_num+1):end)=-1;
R = R';

% initialize a starting point
if opts.init==2
    W0 = zeros(dimension, task_num);
elseif opts.init == 0
    W0 = W0_prep;
else
    if isfield(opts,'W0')
        W0=opts.W0;
        if (nnz(size(W0)-[dimension, task_num]))
            error('\n Check the input .W0');
        end
    else
        W0=W0_prep;
    end
end

bFlag=0; % this flag tests whether the gradient step only changes a little

Wz= W0;
Wz_old = W0;

t = 1;
t_old = 0;

iter = 0;
gamma = 1;
gamma_inc = 2;

while iter < opts.maxIter
    alpha = (t_old - 1) /t;

    Ws = (1 + alpha) * Wz - alpha * Wz_old;

    % compute function value and gradients of the search point
    gWs  = gradVal_eval(Ws);
    Fs   = funVal_eval (Ws);

    while true
        Wzp = ReweightFGLasso_projection(Ws - gWs/gamma, rho1 ./ gamma, rho2 ./ gamma);
        Fzp = funVal_eval  (Wzp);

        delta_Wzp = Wzp - Ws;
        r_sum = norm(delta_Wzp, 'fro')^2;
        %         Fzp_gamma = Fs + trace(delta_Wzp' * gWs)...
        %             + gamma/2 * norm(delta_Wzp, 'fro')^2;
        Fzp_gamma = Fs + sum(sum(delta_Wzp.* gWs))...
            + gamma/2 * r_sum;

        if (Fzp <= Fzp_gamma)
            break;
        else
            gamma = gamma * gamma_inc;
        end
    end

    Wz_old = Wz;
    Wz = Wzp;

    funcVal = cat(1, funcVal, Fzp + nonsmooth_eval(Wz, rho1, rho2));

    if (bFlag)
        % fprintf('\n The program terminates as the gradient step changes the solution very small.');
        break;
    end

    % test stop condition.
    switch(opts.tFlag)
        case 0
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <= opts.tol)
                    break;
                end
            end
        case 1
            if iter>=2
                if (abs( funcVal(end) - funcVal(end-1) ) <=...
                        opts.tol* funcVal(end-1))
                    break;
                end
            end
        case 2
            if ( funcVal(end)<= opts.tol)
                break;
            end
        case 3
            if iter>=opts.maxIter
                break;
            end
    end

    iter = iter + 1;
    t_old = t;
    t = 0.5 * (1 + (1+ 4 * t^2)^0.5);

end

W = Wzp;

% private functions

    function [Wp] = ReweightFGLasso_projection (W, lambda_1, lambda_2)
        % solve it in row wise, since that
        % \sum_i^d rho1{i} * \|W(i, :)\|_1 is row coupled.
        % for each row we need to solve the proximal opterator
        % W(i, :) = argmin_w { 0.5 \|w - v\|_2^2
        %            + lambda_1{i} * \|w\|_1
        %            + lambda_2{i} * \|R * w\|_1 }
        % NOTE: Here the R is t-1 * t, the outside R is t * t-1, and
        Wp = zeros(size(W));

        for i = 1 : size(W, 1)
            v = W(i, :);

            w0 = zeros(length(v)-1, 1);
            w = flsa(v, w0,  lambda_1(i), lambda_2(i), length(v), 1000, 1e-9, 1, 6);

            Wp(i, :) = w';
        end
    end

% smooth part gradient.
    function [grad_W] = gradVal_eval(W)
        if opts.pFlag
            grad_W = zeros(size(W));
            parfor i = 1:task_num
                grad_W(:, i) = X{i}*(X{i}' * W(:,i)-Y{i});
            end
        else
            grad_W = [];
            for i = 1:task_num
                grad_W = cat(2, grad_W, X{i}*(X{i}' * W(:,i)-Y{i}) );
            end
        end
    end

% smooth part gradient.
    function [funcVal] = funVal_eval (W)
        funcVal = 0;
        if opts.pFlag
            parfor i = 1: task_num
                funcVal = funcVal + 0.5 * norm (Y{i} - X{i}' * W(:, i))^2;
            end
        else
            for i = 1: task_num
                funcVal = funcVal + 0.5 * norm (Y{i} - X{i}' * W(:, i))^2;
            end
        end
    end

    function [non_smooth_value] = nonsmooth_eval(W, rho_1, rho_2)
        non_smooth_value = 0;
        for i = 1 : size(W, 1)
            w = W(i, :);
            non_smooth_value = non_smooth_value ...
                + rho_1(i) * norm(w, 1) + rho_2(i) * norm(R * w', 1);
        end
    end

end

author: jiayuzhou

repo: MALSAR

save_path: github-repos/MATLAB/jiayuzhou-MALSAR


text:
classdef Absolutepath_Test < TestCase
    methods
        function this = Absolutepath_Test(name)
            this = this@TestCase(name);
        end

        function testAbsolutepath(~)
            tmp = pwd;
            if ispc; tmp(1) = upper(tmp(1)); end;
            assertEqual(tmp(1:find(tmp==filesep,1,'last')-1),  absolutepath('..'));
            assertEqual(tmp,  absolutepath('.'));
            assertEqual([tmp filesep 'runSmallTests.m'],  absolutepath('./runSmallTests.m'));
            assertEqual([tmp filesep 'runSmallTests.m'],  absolutepath('../simulation/./runSmallTests.m'));
            assertEqual([tmp filesep 'runSmallTests.m'],  absolutepath('../simulation/../simulation/./runSmallTests.m'));
            assertEqual([tmp filesep 'output' filesep 'runSmallTests'],  absolutepath('../simulation/../simulation/./output/runSmallTests'));
            assertEqual([tmp filesep 'output' filesep 'runSmallTests'],  absolutepath('..//simulation/../simulation/./output/runSmallTests'));
            tic; absolutepath('../simulation/../simulation/./output/runSmallTests'); toc
        end
    end
end

author: CovertLab

repo: WholeCell

save_path: github-repos/MATLAB/CovertLab-WholeCell


text:
function [ boxlayout,surface_labels ] = getspatiallayout(imdir,imagename,workspcdir)
% GETSPATIALLAYOUT Given an image, estimate its spatial layout, consisting
%   of a boxlayout and pixel labels of different surfaces. Please refer to
%   the readme file provided with this software for detailed meaning of
%   these outputs.
%
% USAGE: [ boxlayout,surface_labels ] =
%               getspatiallayout(imdir,imagename,workspcdir)
%
% INPUT:
%   imdir - directory containing the original image to be processed. Use '/'
%       (not '\') to separate directories.
%   imagename - original image name
%   workspcdir - directory to hold intermediate results. Internally, two
%       independent directories are created inside this directory to hold
%       visualizations and data files.
%
% OUTPUT:
%   binlayout - structure containing the estimated approximation of the
%       indoor scene as a 3D box, with the following fields -
%   .polyg - cell array of size (n x 5), each containing coordinates of
%       planar surfaces of the room (left, right, middle walls, and, floor
%       and ceiling). The (i,j)-th cell entry corresponds to i-th box
%       layout hypothesis, and j-th plane of that hypothesis.
%   .init - An array of size (n x 2). The i-th row (score,index) represents
%       score and index of a hypothesis in the polyg cell array above, as
%       estimated by our algorithm.
%   .reestimated - An array of size (n x 2). The i-th row (score,index) represents
%       score and index of a hypothesis in the polyg cell array above, as
%       estimated by the second iteration of our algorithm.
%   surface_labels - structure containing the per pixel likelihood of each of
%       the surfaces, and objects, with following fields -
%   .init - An array of size (m x 7), where m is th total number of
%       superpixels detected in the image. The (i,j)-th entry represents
%       the likelihood of the i-th superpixel being labeled as j-th
%       surface.
%
% Copyright (c) Varsha Hedau, University of Illinois at Urbana Champaign,
% 2011.
% To use this code or any part of the package, for any purpose, you agree
% to abide by the restrictions set forth in the file Readme.txt

outimgdir = [workspcdir 'Images/'];
if ~exist(outimgdir,'dir')
    mkdir(outimgdir);
end
workspcdir = [workspcdir 'data/'];
if ~exist(workspcdir,'dir')
    mkdir(workspcdir);
end

% tempdir='./sptiallayouttempworkspace/';
%Compute vps
boxlayout=[];
surface_labels=[];

[vp p All_lines]=getVP(imdir,imagename,0,workspcdir);

img=imread([imdir imagename]);
[h w kk]=size(img);
VP=vp;
if numel(VP)<6
    return;
end

vp=[VP(1) VP(2);VP(3) VP(4);VP(5) VP(6)];
[vp P]=ordervp(vp,h,w,p);
[vv linemem]=max(P,[],2);
vpdata.vp=vp;
vpdata.lines=All_lines;
vpdata.linemem=linemem;
vpdata.dim=[h w];
%visvp(vpdata,img)
%Get segmentation and GC surface confidence maps
% sigma=num2str(0.8);
% k1=num2str(100);
% min1=num2str(100);
% inputim=[workspcdir imagename(1:end-4) '.ppm'];
% outputim=[workspcdir imagename(1:end-4) '.pnm'];
% [s w]=system(['./segment' ' ' sigma ' '  k1  ' ' min1 ' ' inputim ' ' outputim ]);

segext='pnm';
nsegments=[5 15 25 35 40 60 80 100];
fn=['../Imsegs/' imagename(1:end-4) '.' segext];
imseg = processSuperpixelImage(fn);

tic
imdata = mcmcComputeImageData(im2double(img), imseg);% made changes here
toc

load(fullfile('../LabelClassifier/', 'Classifiers_gc.mat'));

spfeatures = mcmcGetAllSuperpixelData(imdir, imseg);
[efeatures, adjlist] = mcmcGetAllEdgeData(spfeatures, imseg(1));

nsegments=[5 15 25 35 40 60 80 100];

pE{1} = test_boosted_dt_mc(eclassifier, efeatures{1});
pE{1} = 1 ./ (1+exp(ecal(1)*pE{1}+ecal(2)));
smaps{1} = msCreateMultipleSegmentations(pE{1}, adjlist{1}, ...
    imseg(1).nseg, nsegments);

for k = 1:numel(nsegments)
    if max(smaps{1}(:, k))>0
        segfeatures{1, k} = mcmcGetSegmentFeatures(imseg, ...
            spfeatures{1}, imdata, smaps{1}(:, k), (1:max(smaps{1}(:, k))));
    end
end

%Get surface label confidences initial from GC
normalize = 1;
pg=zeros(imseg.nseg,7);%7 labels
%get P(L/I)
pg = msTest(imseg, segfeatures, smaps, ...
    labelclassifier, segclassifier,normalize);

filename=fullfile(workspcdir, [imagename(1:end-4) '_lc_gc.mat' ]);
save(filename,'pg');

%visualize

%Compute intergral images for features
tic
[integData]=getIntegralimages([imagename],vpdata,imseg,500,workspcdir,imdir);
toc

%Get candidate layouts and their features
[polyg, Features] = getcandboxlayout( vpdata.vp,vpdata.dim(1),vpdata.dim(2),integData);

%Get initial estimate
Features1=Features;

load ../LearntClassifiers/pf_i.mat

lay_score=[];
numL=size(Features,1);
% change features
if(numel(weights)==14)
    tmpf1=sum(Features(:,1:5).*Features(:,11:15),2);
    tmpf2=sum(Features(:,1:5).*Features(:,16:20),2);
    tmpf3=sum(Features(:,6:10).*Features(:,11:15),2);
    tmpf4=sum(Features(:,6:10).*Features(:,16:20),2);
    Features=[Features(:,1:10) tmpf1 tmpf2 tmpf3 tmpf4];

end
%evaluate
score=repmat(weights,[numL,1]).*Features;
score=sum(score,2);
[vv ii]=sort(score,'descend');

boxlayout.polyg=polyg;
boxlayout.init=[vv ii];

load ../LearntClassifiers/pf_il.mat
Features=Features1;
if(numel(weights)==59)
    tmpf1=sum(Features(:,1:5).*Features(:,11:15),2);
    tmpf2=sum(Features(:,1:5).*Features(:,16:20),2);
    tmpf3=sum(Features(:,6:10).*Features(:,11:15),2);
    tmpf4=sum(Features(:,6:10).*Features(:,16:20),2);
    Features=[Features(:,1:10) tmpf1 tmpf2 tmpf3 tmpf4 Features(:,31:75) ];%Features(:,114:122)];

end
score=repmat(weights,[numL,1]).*Features;
score=sum(score,2);
[vv ii]=sort(score,'descend');

boxlayout.reestimated=[vv ii];

lay_scores=[vv ii];
save([workspcdir imagename(1:end-4) '_layres.mat'],'polyg','lay_scores');%,'avg_pg');

% figure(101);
% drawnow;
% for lay=1:25
%     layoutid=ii(lay);
%     Polyg=[];
%     for fie=1:5
%         Polyg{fie}=[];
%         if size(polyg{layoutid,fie})>0
%             Polyg{fie}=polyg{layoutid,fie};
%         end
%     end
%
%     tempimg=displayout(Polyg,w,h,img);
%     subplot(5,5,lay);imshow(uint8(tempimg),[]);title(num2str(vv(lay)));
% end
% saveas(101,[outimgdir imagename(1:end-4) '_boxlayouts.png']);
%
% Polyg=[];
% for fie=1:5
%     Polyg{fie}=[];
%     if size(polyg{ii(1),fie})>0
%         Polyg{fie}=polyg{ii(1),fie};
%     end
% end
% ShowGTPolyg(img,Polyg,103);
% saveas(103,[outimgdir imagename(1:end-4) '_boxlayout.png']);

%Re-Compute Surface labels (GC+box layout features)
load(['../LabelClassifier/' 'Classifiers_stage2.mat']);
tic
xspfields=[];%save per image
numSup=imseg.nseg;
for lay=1:numL %each layout
    xspf=[];
    for  supno=1:numSup
        tempbndy=imdata.tracedbndy{supno}{1};
        if size(tempbndy,1) > 30

            YY1=tempbndy(1:10:end,1);XX1=tempbndy(1:10:end,2);
        else

            YY1=tempbndy(:,1);XX1=tempbndy(:,2);
        end

        for fi=1:5  %each field
            xarea=0;
            if size(polyg{lay,fi},1)>0

                XX2=polyg{lay,fi}(:,1);
                YY2=polyg{lay,fi}(:,2);

                [in on]=inpolygon(XX1,YY1,[XX2;XX2(1)],[YY2;YY2(1)]);

                if numel(find(in==1))==length(in)
                    X0=XX1;Y0=YY1;
                    xarea=polyarea([X0;X0(1)],[Y0;Y0(1)]);
                elseif  numel(find(in==1))==0

                    X0=[];Y0=[];
                    xarea=0;
                else
                    xarea=polyintarea(XX1,YY1,XX2,YY2,0);

                end

            end
            xspf(supno,fi)=xarea;
        end

    end
    xspfields{lay}=xspf;
end
toc

nsp=imseg.nseg; %num of suppixels
smap = [1:nsp];
smaps{1} = smap(:);

clear segfeatures
for k = 1: 1%numel(nsegments)
    features = mcmcGetSegmentFeatures(imseg, ...
        spfeatures{1}, imdata, smaps{1}(:, k), (1:max(smaps{1}(:, k))));

    for lay=1:numL
        tempfeatures=features;
        [fieldfeatures]=segfieldsfeat_sup(imseg,imseg.nseg,smaps{1}(:,k),xspfields{lay});
        tempfeatures(:,95:100)=fieldfeatures(1:size(features,1),:);
        tempfeatures(:,101:106)=pg{1}(:,1:6);
        segfeatures{lay,k}  = tempfeatures;
    end
end

if numel(vv)>100
    A = [vv(1) 1;vv(100) 1];
else
    A = [vv(1) 1;vv(end) 1];
end
b = [0;log(49)];
params = inv(A)*b;

%     conf = 1./(1+exp(params(1)*vv(1:5)+params(2)));
conf = 1./(1+exp(params(1)*vv(1)+params(2)));
conf=conf./sum(conf);

avg_pg=zeros(imseg.nseg,7);%7 labels

tic
for j=1:1 %take best 5 scored layouts
    lay=ii(j);

    %get P(L/F,I)
    pg = msTest(imseg, segfeatures(lay, :), smaps, ...
        labelclassifier);%, segclassifier,normalize);
    avg_pg=avg_pg+pg{1}.*conf(j);

end
toc

filename=fullfile(workspcdir,[imagename(1:end-4) '_lc_st2.mat' ]);
save(filename,'avg_pg');

surface_labels.restimated={avg_pg};
surface_labels.init=pg;

rl=[255  255 0  0 0 255 0];
gl=[255 0 255  255 0 0 0];
bl=[0  0 255 0 255 255 0];
rl1=[255 202 132 255 0 255 0 191 0 0];
gl1=[236 255 112 255 191 0 255 21 0 0];
bl1=[139 112 255  255 255 0 0 133 255 0];
cimages = msPg2confidenceImages(imseg,{avg_pg});
[aa indd]=max(cimages{1}(:,:,1:6),[],3);

figure(102);
clear mask_color;
mask_r = rl(indd);
mask_g = gl(indd);
mask_b = bl(indd);
mask_color(:,:,1) = mask_r;
mask_color(:,:,2) = mask_g;
mask_color(:,:,3) = mask_b;

hsvmask=rgb2hsv(mask_color);
hsvmask(:,:,3)=aa*255;
%     hsvmask(:,:,2)=aa;
mask_color=hsv2rgb(hsvmask);

tempimg = double(img)*0.5 + mask_color*0.5;
%         tempimg =  mask_color;
imshow(uint8(tempimg));
saveas(102,[outimgdir imagename(1:end-4) '_surfacelabels.png']);

end

author: zouchuhang

repo: LayoutNet


text:
function res = exp_same(a,b)
% Check whether two expressions are structurally identical.
% Result = exp_same(Expression-A, Expression-B)
%
% Differs from isequal() in the following ways:
%  a) ignores the values associated with impure expressions
%  b) treats @x as equal to exp_symbol('x'); therefore, @x is a shortcut notation for exp_symbol('x').
%
% See also:
%   exp_match
%
%                                Christian Kothe, Swartz Center for Computational Neuroscience, UCSD
%                                2010-04-19

if ~exp_beginfun('symbolic') return; end

res = utl_same(a,b);

exp_endfun;

repo: REST