text:
function nframe=eyematch(frame,thresh,template);

[face_a,skin_region]=face(frame);

im1 = double(rgb2gray(face_a));
img=im1-mean(mean(im1));
template=template-mean(mean(template));

C=normxcorr2(template,img);
Csub = C((size(template,1)-1)/2+1:size(img,1)-1,(size(template,2)-1)/2+1:size(img,2)-1);
figure;
% imshow(Csub, []), pixval;

BW = Csub;
BW(find(BW<thresh))=0;
BW(find(BW>thresh))=1;
% figure, imshow(BW);

[L,n] = bwlabel(BW);
stats = regionprops(L, 'centroid');
% Draw an asterisk

figure(1);imshow(frame);
centroids = cat(1, stats.Centroid);
hold on
plot(centroids(:,1), centroids(:,2), 'r*');
hold off


text:
%% Check rate of convergence for 3D H1 interface problem
%     -div(A grad u)  = f,    x\in \Omega
%      where A is a piecewise constant on Omega^+ and Omega^-.
%
% Domain: Rectangular domain: [xmin,xmax] X [ymin,ymax] X [zmin,zmax]
% Mesh: Cartesian triangular mesh.
% Method: PP-IFE

%% Geometry and Boundary Conditions
% clear
% close all
%clc

domain = [-1,1,-1,1,-1,1];
bc = [1,1,1,1,1,1]; % Dirichelet BC

%% Finite Element Type
femtype = 'P1';
disp(['FEM Type =  Conforming ', femtype]);

%% Initial Partition
nx0 = 10;
ny0 = nx0;
nz0 = nx0;

%% Task
showErr = 0;
showMesh = 0;
computErr = 1;

%% PDE
test = 6;
switch test
    case 0
        pde = poissonNonPoly3D;
    case 1 % circular interface
        r = pi/4; bm = 1; bp = pi/(2*r^2);
        %r = sqrt(pi/2); bm = 1; bp = pi/(2*r^2);domain = [-2,2,-2,2,-2,2];
        x0 = 0; y0 = 0; z0 = 0; rx = r; ry = r; rz = r;
        pde = elli3DcircIntf(bm,bp,r,x0,y0,z0,rx,ry,rz);
    case 2 % orthotorus interface
        domain = [-1.2,1.2,-1.2,1.2,-1.2,1.2];
        bm = 1; bp = 1;
        rx = 1; ry = 0.075; rz = 3;
        pde = elli3DorthocircIntf(bm,bp,rx,ry,rz);
    case 3 % line interface
        bm = 1; bp = 1;
        rx = 1; ry = 10^(-4); rz = -1; cx = 10^(-4); cy = 0; cz = 0;
        cm = 1; cp = 1; a = 1;
        pde = elli3DlinIntf(bm,bp,cx,cy,cz,rx,ry,rz,a,cm,cp);
    case 4 % line interface but with zero boundary conditions
        bm = 1; bp = 10; delta = pi/100; % delta can not be -1,1
        pde = elli3DlinIntf2(bm,bp,delta);
    case 5 % circular interface
        r = pi/4; bm = 1; bp = 10;domain = [-1,1,-1,1,-1,1];
        x0 = 0; y0 = 0; z0 = 0;
        pde = elli3DcircIntf5(bm,bp,r,x0,y0,z0);
    case 6 % circular interface
        bm = 1; bp = 1;domain = [-1.3,1.3,-1.3,1.3,-1.3,1.3];
        x1 = -0.3; y1 = 0; z1 = 0; r11 = pi/5; r12 = 0.2;
        x2 =  0.3; y2 = 0; z2 = 0; r21 = pi/5; r22 = 0.2;
        pde = elli3DtorusTwin(bm,bp,x1,y1,z1,r11,r12,x2,y2,z2,r21,r22);
end
%% PPIFEM Type
PPtype = 'S';
disp(['PPIFEM Type = ', PPtype]);
sig = max(bm,bp)/min(bm,bp);

%% Max Iteration
maxIt = 10;
time = zeros(9,maxIt);

for i = 1:maxIt

    %% 1. Generate Mesh
    tic
    nx = nx0 + 10*(i-1); h = (domain(2) - domain(1))/nx;
    ny = ny0 + 10*(i-1);
    nz = nz0 + 10*(i-1);
    time(1,i) = nx;
    disp(' ')
    disp('**************************************************************************************')
    disp(['Partition =  ',int2str(nx),' X ',int2str(ny),' X ',int2str(nz)]);
    disp(' ')

    mesh = genMesh3D(domain, nx, ny, nz);
    mesh = enrichMesh3D(mesh,2); % Mesh detail level = 2 (for PPIFE).
    disp(['number of element =  ', int2str(length(mesh.t))]);
    mesh = genIntfMesh3D(mesh,pde.intf);
    disp(['number of interface element =  ', int2str(-min(mesh.tLoc)),...
        ', is ', num2str(100*-min(mesh.tLoc)/length(mesh.t)), '% of all elements']);
    if showMesh == 1
        tetramesh(mesh.t,mesh.p,ones(size(mesh.t)));
    end
    time(2,i) = toc;

    %% 2. Generate FEM DoF
    tic
    fem = genFEM3D(mesh,femtype,bc);
    disp(['number of DoF =  ', int2str(length(fem.p))]);
    time(3,i) = toc;
    tic
    femI = genP1IFEM3D(mesh,fem,bm,bp);
    time(4,i) = toc;
    tic
    femIF = genP1IFEM3DFace(mesh,fem,femI);
    time(5,i) = toc;

    %% 3. Assemble Matrix
    tic
    disp(' '); disp('Start Assembling Matrix');
    matrix = genMatEllPPIFE3D(pde,mesh,fem,femI,femIF,PPtype,sig);
    time(6,i) = toc;

    %% 4. Solve the linear system Au = f
    tic
    disp(' ')
    if strcmp(PPtype,'S')
        disp('Start Solving Linear System: using PCG with ichol precond');
        L = ichol(matrix.A);
        u = pcg(matrix.A,matrix.f,1e-8,300,L,L');
    elseif strcmp(PPtype,'N') || strcmp(PPtype,'I')
        disp('Start Solving Linear System: using GMRES with ilu precond');
        [L,U] = ilu(matrix.A,struct('milu','row'));
        u = gmres(matrix.A,matrix.f,[],1e-8,300,L,U);
    end
    time(7,i) = toc;
    tu = pde.exactu(fem.p(:,1),fem.p(:,2),fem.p(:,3));
    uh = tu; uh(fem.mapper) = u;

    %% 5. Postprocess: Calculating Errors
    tic
    errND = max(abs(uh - tu)); % Error on nodes
    err.nd = errND; err.inf = 0; err.l2 = 0; err.h1 = 0;
    if computErr == 1
        disp(' ');  disp('Start computing error in Inf norm');
        [errInfN,errInfI] = getErrInfIFE3D(uh,pde.exactu,mesh,femI,fem,[0,0,0]);
        eNorm = 'L2'; disp(['Start computing error in ',eNorm,'  norm']);
        [errL2,errL2K] = getErrIFE3D(uh, pde, mesh, fem, femI,eNorm);
        eNorm = 'H1x'; disp(['Start computing error in ',eNorm,' norm']);
        [errH1x,errH1xK] = getErrIFE3D(uh, pde, mesh, fem, femI, eNorm);
        eNorm = 'H1y'; disp(['Start computing error in ',eNorm,' norm']);
        [errH1y,errH1yK] = getErrIFE3D(uh, pde, mesh, fem, femI, eNorm);
        eNorm = 'H1z'; disp(['Start computing error in ',eNorm,' norm']);
        [errH1z,errH1zK] = getErrIFE3D(uh, pde, mesh, fem, femI, eNorm);

        disp(' ')
        disp(['Max error on Regular Elements   = ', num2str(errInfN)])
        disp(['Max error on Interface Elements = ', num2str(errInfI)])
        err.inf = max([errND,errInfN,errInfI]);
        err.l2 = errL2;
        err.h1 = sqrt(errH1x^2+errH1y^2+errH1z^2);
    end
    time(8,i) = toc;
    time(9,i) = 1e6*sum(time(2:8,i))/length(mesh.t);

    %% 6: Display Output
    disp(' ')
    disp('Errors')
    disp('Node        Inf norm    L2 norm     H1 norm')
    formatSpec = '%6.4e  %6.4e  %6.4e  %6.4e\n';
    fprintf(formatSpec, err.nd, err.inf, err.l2, err.h1)

    if i > 1
        format short
        rNd = log(err0.nd/err.nd)./log(h0/h);
        rInf = log(err0.inf/err.inf)./log(h0/h);
        rL2 = log(err0.l2/err.l2)./log(h0/h);
        rH1 = log(err0.h1/err.h1)./log(h0/h);
        disp(' ')
        disp('Convergence Rate')
        disp('Node        Inf norm    L2 norm     H1 norm')
        formatSpec = '%6.4f      %6.4f      %6.4f      %6.4f\n';
        fprintf(formatSpec, rNd, rInf, rL2, rH1)
    end
    err0 = err; h0 = h;

    disp(' '); disp('CPU Time')
    disp('   N     Mesh     FEM      FEMI    FEMIF    Matrix   Solve    Error    Time/1M cell')
    formatSpec = '%4i  %7.2f  %7.2f  %7.2f  %7.2f  %7.2f  %7.2f  %7.2f   %7.2f\n';
    fprintf(formatSpec, time(:,1:i))

    %% 6. Plot Solution and Error
    if showErr == 1
    end

end

repo: ifem


text:
function quadmom_test ( )

%*****************************************************************************80
%
%% QUADMOM_TEST tests the QUADMOM library.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    05 October 2013
%
%  Author:
%
%    John Burkardt
%
  timestamp ( );
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADMOM_TEST\n' );
  fprintf ( 1, '  MATLAB version\n' );
  fprintf ( 1, '  Test the QUADMOM library.\n' );

  quadmom_test01 ( );
  quadmom_test02 ( );
  quadmom_test03 ( );
  quadmom_test04 ( );
  quadmom_test05 ( );
  quadmom_test06 ( );
  quadmom_test07 ( );
  quadmom_test08 ( );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'QUADMOM_TEST\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( );

 return
end

author: johannesgerer

repo: jburkardt-m


text:
% Gnufft_test0.m
% Basic tests of small Gnufft object

if ~isvar('A') % create Gnufft class object
	omega = linspace(0, 10*2*pi, 101)'; % crude spiral:
	omega = pi*[cos(omega) sin(omega)].*omega(:,[1 1])/max(omega);
	if 1 % 2d
		N = [16 14];
		J = [8 6];
		testadj = @(sys) test_adjoint(sys, 'complex', 1);
	else % 3d
		N = [16 14 10];
		J = [8 6 4];
		omega(:,3) = linspace(-pi, pi, size(omega,1));
		testadj = @(sys) test_adjoint(sys, 'big', 1, 'complex', 1);
	end
	K = 2*N;
	mask = true(N);
	mask(:,end) = false;
	cl = 'Fatrix';
	cl = 'fatrix2';
	A = Gnufft(cl, mask, {omega, N, J, K});

	if 1
		fatrix2_tests(A, 'complex', 1, 'tol_gram', 7e-5)
		testadj(A);
	end
end

if ~isvar('W')
	wi = [1:size(omega,1)]';
	W = Gdiag(wi);
end

if 1, printm 'Gnufft gram'
	T = build_gram(A, wi);
	if isa(A, 'Fatrix')
		Fatrix_test_basic(T, mask, 'complex', 1)
	else
		fatrix2_tests(T, 'complex', 1)
	end
	if length(N) == 3
		testadj(T);
	else
		[t0 t1] = test_adjoint(T, 'complex', 1);
		im plc 1 3, im(1, t0), im(2, t1), im(3, t0 - t1')
	end
prompt
end

if length(N) == 2, printm 'T vs A''WA' % slow in 3D
	T2 = T(:,:);
	Af = A(:,:);
	T1 = Af' * diag(wi) * Af;
	max_percent_diff T1 T2
	im plc 1 3, im(1, T1), im(2, T2), im(3, T1 - T2)
%	equivs(y1, y2)
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
% Compute h-step Weisfeiler-Lehman shortest path delta kernel for a set of graphs
% Author: Nino Shervashidze - nino.shervashidze@tuebingen.mpg.de
% Copyright 2010 Nino Shervashidze
% Input: Graphs - a 1xN array of graphs
% 		  Graphs(i).am is the adjacency matrix of the i'th graph,
% 		  Graphs(i).al is the adjacency list of the i'th graph,
%                 Graphs(i).nl.values is a column vector of node
%                 labels for the i'th graph.
%                 Graphs(i).sp (may be missing) is the shortest
%                 path length matrix of the i'th graph
%                 Graphs(i) may have other fields, but they will not be
%                 used here.
%	 h - a natural number: number of iterations of WL
%	 nl - a boolean: 1 if we want to use original node labels, 0 otherwise
%        spprec - a boolean: 1 if shortest paths are precomputed, 0
%        otherwise. Default: 0
% Output: K - a h+1-element cell array of NxN kernel matrices K for
%             each iter = 0,...,h
%         runtime - scalar (total runtime in seconds)

function [K,runtime] = WLspdelta(Graphs,h,nl,spprec)
% THIS SOURCE CODE IS SUPPLIED "AS IS" WITHOUT WARRANTY OF ANY KIND, AND
% ITS AUTHOR AND THE JOURNAL OF MACHINE LEARNING RESEARCH (JMLR) AND
% JMLR'S PUBLISHERS AND DISTRIBUTORS, DISCLAIM ANY AND ALL WARRANTIES,
% INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES OF MERCHANTABILITY
% AND FITNESS FOR A PARTICULAR PURPOSE, AND ANY WARRANTIES OR NON
% INFRINGEMENT. THE USER ASSUMES ALL LIABILITY AND RESPONSIBILITY FOR
% USE OF THIS SOURCE CODE, AND NEITHER THE AUTHOR NOR JMLR, NOR JMLR'S
% PUBLISHERS AND DISTRIBUTORS, WILL BE LIABLE FOR DAMAGES OF ANY KIND
% RESULTING FROM ITS USE. Without limiting the generality of the foregoing,
% neither the author, nor JMLR, nor JMLR's publishers and distributors,
% warrant that the Source Code will be error-free, will operate without
% interruption, or will meet the needs of the user.

N=size(Graphs,2);
Lists = cell(1,N);
K = cell(1,h+1);
if nargin<4 spprec=0; end
n_nodes=0;
% compute adjacency lists and n_nodes, the total number of nodes in the dataset
for i=1:N
  Lists{i}=Graphs(i).al;
  n_nodes=n_nodes+size(Graphs(i).am,1);
end
% compute/copy shortest path matrices
maxpath = 0;
if ~spprec
  disp('computing shortest paths...');
  for i=1:N
    Ds{i}=floydwarshall(Graphs(i).am);
    aux=max(Ds{i}(~isinf(Ds{i})));
    if aux > maxpath
      maxpath=aux;
    end
  end
else
  for i=1:N
    Ds{i}=Graphs(i).sp;
    aux=max(Ds{i}(~isinf(Ds{i})));
    if aux > maxpath
      maxpath=aux;
    end
  end
end
t=cputime; % for measuring runtime

%%% INITIALISATION
% initialize the node labels for each graph with their labels or
% with degrees (for unlabeled graphs).
label_lookup=containers.Map();
label_counter=uint32(1);
% label_lookup is an associative array, which will contain the
% mapping from multiset labels (strings) to short labels (integers)
if nl==1
  for i=1:N
    % the type of labels{i} is uint32, meaning that it can only handle
    % 2^32 labels and compressed labels over all iterations. If
    % more is needed, switching (all occurences of uint32) to
    % uint64 is a possibility
    labels{i}=zeros(size(Graphs(i).nl.values,1),1,'uint32');
    for j=1:length(Graphs(i).nl.values)
      str_label=num2str(Graphs(i).nl.values(j));
      % str_label is the node label of the current node of the
      % current graph converted into a string
      if ~isKey(label_lookup, str_label)
        %str_label
        %label_counter
        label_lookup(str_label)=label_counter;
        labels{i}(j)=label_counter;
        label_counter=label_counter+1;
      else
        labels{i}(j)=label_lookup(str_label);
      end
    end
  end
else
  for i=1:N
    labels{i}=uint32(full(sum(Graphs(i).am,2)));
    for j=1:length(labels{i})
      str_label=num2str(labels{i}(j));
      % str_label is the node label of the current node of the
      % current graph converted into a string
      if ~isKey(label_lookup, str_label)
        label_lookup(str_label)=label_counter;
        labels{i}(j)=label_counter;
        label_counter=label_counter+1;
      else
        labels{i}(j)=label_lookup(str_label);
      end
    end
  end
end
L=double(label_counter)-1;
clear Graphs;
disp(['Number of original labels: ',num2str(L)]);
disp(['Number of potential shortest path features: ',num2str((maxpath+1)*L*(L+1)/2)]);
sp=sparse((maxpath+1)*L*(L+1)/2,N);
for i=1:N
  labels_aux=repmat(double(labels{i}),1,length(labels{i}));
  a=min(labels_aux, labels_aux');
  b=max(labels_aux, labels_aux');
  I=triu(~(isinf(Ds{i})));
  Ind=Ds{i}(I)*L*(L+1)/2+(a(I)-1).*(2*L+2-a(I))/2+b(I)-a(I)+1;
  minind=min(Ind);
  diff=max(Ind)-minind;
  aux=accumarray(Ind,ones(nnz(I),1),[],[],[],(minind > 5000 || diff > 3000));
  % sparse of full accumarray depending on the range of values in Ind
  % (and based on empirical observations on the speed of accumarray)
  sp(Ind,i)=aux(Ind);
end
sp=sp(sum(sp,2)~=0,:);
K{1}=full(sp'*sp);

%%% MAIN LOOP
iter=1;
new_labels=labels;
while iter<=h
  disp(['iter=',num2str(iter)]);
  % create an empty lookup table
  label_lookup=containers.Map();
  label_counter=uint32(1);
  for i=1:N
    for v=1:length(Lists{i})
      % form a multiset label of the node v of the i'th graph
      % and convert it to a string
      long_label=[labels{i}(v), sort(labels{i}(Lists{i}{v}))'];
      long_label_2bytes=typecast(long_label,'uint16');
      long_label_string=char(long_label_2bytes);
      % if the multiset label has not yet occurred, add it to the
      % lookup table and assign a number to it
      if ~isKey(label_lookup, long_label_string)
        label_lookup(long_label_string)=label_counter;
        new_labels{i}(v)=label_counter;
        label_counter=label_counter+1;
      else
        new_labels{i}(v)=label_lookup(long_label_string);
      end
    end
  end
  L=double(label_counter)-1;
  disp(['Number of compressed labels: ',num2str(L)]);
  disp(['Number of potential shortest path features: ',num2str((maxpath+1)*L*(L+1)/2)]);
  labels=new_labels;
  sp=sparse((maxpath+1)*L*(L+1)/2,N);
  for i=1:N
    labels_aux=repmat(double(labels{i}),1,length(labels{i}));
    a=min(labels_aux, labels_aux');
    b=max(labels_aux, labels_aux');
    I=triu(~(isinf(Ds{i})));
    Ind=Ds{i}(I)*L*(L+1)/2+(a(I)-1).*(2*L+2-a(I))/2+b(I)-a(I)+1;
    minind=min(Ind);
    diff=max(Ind)-minind;
    aux=accumarray(Ind,ones(nnz(I),1),[],[],[],(minind > 5000 || diff > 3000));
    % sparse of full accumarray depending on the range of values in Ind
    % (and based on empirical observations on the speed of accumarray)
    sp(Ind,i)=aux(Ind);
  end
  sp=sp(sum(sp,2)~=0,:);
  K{iter+1}=K{iter}+full(sp'*sp);
  iter=iter+1;
end
runtime=cputime-t; % computation time of K
end

function [D] = floydwarshall(A, sym, w)
% Input: A - nxn adjacency matrix,
%           sym - boolean, 1 if A and w symmetric
%	    w - nxn weight matrix
% Output: D - nxn distance matrix

n = size(A,1); % number of nodes
D=zeros(n,n);

if nargin<2 % if the graph is not weighted and we have no information about sym, then
  sym=1;
  w=A;
end

if nargin<3 % if the graph is not weighted, then
  w=A;
end

D=w.*A;
D(A+diag(repmat(Inf,n,1))==0)=Inf;
D=full(D.*(ones(n)-eye(n))); % set the diagonal to zero

if sym % then it is a bit faster
  for k=1:n
    Daux=repmat(full(D(:,k)),1,n);
    Sumdist=Daux+Daux';
    D(Sumdist<D)=Sumdist(Sumdist<D);
  end
else
  for k=1:n
    Daux1=repmat(full(D(:,k)),1,n);
    Daux2=repmat(full(D(k,:)),n,1);
    Sumdist=Daux1+Daux2;
    D(Sumdist<D)=Sumdist(Sumdist<D);
  end
end
end

author: muhanzhang

repo: DGCNN

save_path: github-repos/MATLAB/muhanzhang-DGCNN


text:
function calcAllFusedTextures_HNissam(pathWORK,nPatient,roiNumb,CTinv_cell,CTweight_mat,R_mat,scale_cell,algo_cell,Ng_mat,fusWeights)
% -------------------------------------------------------------------------
% function calcAllFusedTextures_HNissam(pathWORK,nPatient,roiNumb,CTinv_cell,CTweight_mat,R_mat,scale_cell,algo_cell,Ng_mat)
% -------------------------------------------------------------------------
% DESCRIPTION:
% This function computes FUSED texture features for all patients, for all
% different combinations of the following texture extraction parameters:
% - CT Inv.: Inversion of CT intensities in the PET/CT fusion process.
% - CT weight: Weight given to MRI wavelet low-pass sub-bands in the
%               PET/CT fusion process.
% - WPBF ratio 'R': Ratio specifying the ratio of weight given to band-pass
%                   coefficients to the weigth given to the rest of
%                   coefficients (HHH and LLL) in the wavelet band-pass
%                   filtering process.
% - Scale: Resolution at which the ROI is isotropically resampled.
% - Quant. Algo.: Quantization algorithm used on analyzed ROI.
% - Ng: Number of gray-levels in the quantization process.
%
% Different extraction parameters are passed as arrays or cells in the
% function in order to test all possible combinations. This function is
% used for FUSED scans specifically. See Ref. [1,2] and 'prepareVolume.m'
% for more details.
%
% Texture features are computed for all head and neck (HN) DICOM  imaging
% data downloaded from The Cancer Imaging Archive (TCIA) website at:
% <http://dx.doi.org/10.7937/K9/xxxxxxxxxxxxxxxxxxx, and first organized
% in a 'DATA' directory using the function readAllDICOM_HN.m.  Results are
% then saved in a folder 'TEXTURES' in the HN WORKSPACE.
% -------------------------------------------------------------------------
% REFERENCE:
% [1] Vallieres, M. et al. (2015). FDG-PET/CT radiomics models for the
%     early prediction of different tumour outcomes in head and neck cancer.
%     The Journal of Nuclear Medicine, aa(bb), xxx-yyy.
%     doi:
% [2] Vallieres, M. et al. (2015). A radiomics model from joint FDG-PET and
%     MRI texture features for the prediction of lung metastases in soft-tissue
%     sarcomas of the extremities. Physics in Medicine and Biology, 60(14),
%     5471-5496. doi:10.1088/0031-9155/60/14/5471
% -------------------------------------------------------------------------
% INPUTS:
% - pathWORK: Full path to the HN WORKSPACE directory.
% - nPatient: Number of patients to read.
% - roiNumb: Array of size [N X 2] specifying the contour number
%            corresponding to the GTV in sData{2}.scan.contour for all
%            patients. First column is for the PET data, second column for
%            CT data.
% - CTinv_cell: Cell vector specifying the different 'CT Inv.' values to test.
% - CTweight_mat: Array vector specifying the different 'CT weight' values to test.
% - R_mat: Array vector specifying the different 'R' ratio values to test.
% - scale_cell: Cell vector specifying the different 'Scale' values to test.
% - algo_cell: Cell vector specifying the different 'Quant. Algo.' values to test.
% - Ng_mat: Array vector specifying the different 'Ng' values to test.
% -------------------------------------------------------------------------
% EXAMPLE:
% CTinv_cell = {'NoInv','Inv'};
% CTweight_mat = [1/4,1/3,1/2,2/3,3/4];
% R_mat = [1/2,2/3,1,3/2,2];
% scale_cell = {'pixelW',1,2,3,4,5};
% algo_cell = {'Uniform','Equal','Lloyd'};
% Ng_mat = [8,16,32,64];
% -------------------------------------------------------------------------
% AUTHOR(S): Martin Vallieres <mart.vallieres@gmail.com>
% -------------------------------------------------------------------------
% HISTORY:
% - Creation: July 2015
%--------------------------------------------------------------------------
% STATEMENT:
% This file is part of <https://github.com/mvallieres/radiomics/>,
% a package providing MATLAB programming tools for radiomics analysis.
% --> Copyright (C) 2015  Martin Vallieres
%
%    This package 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 package 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 package.  If not, see <http://www.gnu.org/licenses/>.
% -------------------------------------------------------------------------

startpath=pwd;
cd(pathWORK), cd('TEXTURES'), pathTEXT = pwd;
cd .., cd('DATA'), pathDATA = pwd;

% INITIALIZATION
scans = {'CT'};
nScans = numel(scans);
suffix = ['_TEXT'];

% COMPUTATION
fprintf('\n')
for i = 1:nPatient
    load(['Patient',num2str(i),'_PET']), sDataPET = sData;
    for j = 1:nScans
        nameLoad = ['Patient',num2str(i),'_CT'];
        load(nameLoad), sDataCT = sData;
        tStart = tic;
        fprintf(['\n*********************** COMPUTING TEXTURES: PATIENT %u, FUSED PET/',scans{j},' SCAN ***********************'],i)
        [textures] = calcPatientFusText_HNissam(sDataPET,sDataCT,roiNumb(i,:),CTinv_cell,CTweight_mat,R_mat,scale_cell,algo_cell,Ng_mat,fusWeights);
        cd(pathTEXT), save(['Patient',num2str(i),'_PET_',scans{j},suffix],'textures')
        cd(pathDATA)
        time = toc(tStart);
        fprintf('TOTAL TIME: %.2f seconds\n',time)
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
% Author: Jai Juneja (adapted from video lecture series by Joan Sola)
% Date: 12/02/2013
%
% Given a range-bearing measurement y, obtain position p of observed
% landmark in the robot's frame.
%
% Inputs:
%   y = [d a]       :   Laser range and bearing measurement
%
% Outputs:
%   p_r = [p_x p_y]':   Cartesian position of observation in robot's frame
%   Optional:
%   P_y             :   Jacobian of p wrt. y

function [p_r, P_y] = getInvMeasurement(y)
    d = y(1, :);           % Range measurement
    a = y(2, :);           % Bearing measurement

    p_x = d .* cos(a);   % x-position of observation in robot's frame
    p_y = d .* sin(a);   % y-position of observation in robot's frame

    p_r = [p_x; p_y];

    if nargout > 1      % Compute Jacobian (only works for single measurement)
        P_y = [...
            cos(a)  -d*sin(a)
            sin(a)  d*cos(a)];
    end
end

author: jaijuneja

repo: ekf-slam-matlab


text:
function dist = line_par_point_dist_2d ( f, g, x0, y0, p )

%*****************************************************************************80
%
%% LINE_PAR_POINT_DIST_2D: distance ( parametric line, point ) in 2D.
%
%  Discussion:
%
%    The parametric form of a line in 2D is:
%
%      X = X0 + F * T
%      Y = Y0 + G * T
%
%    We normalize by choosing F*F+G*G=1 and 0 <= F.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 February 2005
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Adrian Bowyer and John Woodwark,
%    A Programmer's Geometry,
%    Butterworths, 1983.
%
%  Parameters:
%
%    Input, real F, G, X0, Y0, the parametric line parameters.
%
%    Input, real P(2,1), the point whose distance from the line is
%    to be measured.
%
%    Output, real DIST, the distance from the point to the line.
%
  dx =   g * g * ( p(1,1) - x0 ) - f * g * ( p(2,1) - y0 );
  dy = - f * g * ( p(1,1) - x0 ) + f * f * ( p(2,1) - y0 );

  dist = sqrt ( dx * dx + dy * dy ) / ( f * f + g * g );

  return
end


text:
function filesize = GetFileSize(fileobject)

if ischar(fileobject)
    if exist(fileobject, 'file') == 2
        fd = fopen(fileobject);
    else
        filesize = -1;
        return;
    end
else
    fd = fileobject;
end

% Since we're changing file position, save the initial one, then restore
% it before exiting
p0 = ftell(fd);

fseek(fd, 0, 'eof');     % move the cursor to the end of the file
filesize = ftell(fd);

% Restore original file position
fseek(fd, p0, 'bof');

if ischar(fileobject)
    fclose(fd);
end

author: BUNPC


text:
function y = barycentricInterpolate(Y,Idx,W)
% y = barycentricInterpolate(Y,Idx,W)
%
% This function interpolates a point y, based on the data in Y and the
% weights (W) on the verticies Y(Idx)
%
% INPUTS:
%   Y = [N, 1] data set.
%   W = [d+1, n] weight to apply to each index
%   Idx = [d+1, n] linear index corresponding to each weight
%
% OUTPUTS:
%   y = [n,d] = function value at query points
%
% NOTES:
%   --> y = f(x), where x is d-dimensional, and y is scalar
%   --> Y = Y(x1,x2,...,xd)
%   --> n = size(Y);
%   --> y(X(:,i)) = W(i,:)*Y(Idx(:,i));
%
% See Also: barycentricWeights

n = size(Idx,2);
y = zeros(n,1);

for query=1:n
    y(query) = dot(W(:,query),Y(Idx(:,query)));
end

end

author: MatthewPeterKelly

repo: dscTutorials


text:
function writechanloc(sqdfile,locfile,format,channum)
% WRITECHANLOC(SQD_FILE, LOCATION_FILE, FORMAT, CHANNEL_LIST); Writes the location of
% the sensors of the channels from the CHANNEL_LIST in the SQD_FILE to LOC_FILE
% following the FORMAT structure
% where FORMAT can be:
%     'xyz' : (Default) Matlab/EEGLAB Cartesian coordinates (NOT EGI Cartesian).
%             x is toward nose; y is toward left ear; z is toward vertex
% Default values:
%   LOCATION_FILE = './testsqd.xyz'
%   FORMAT = 'xyz'
%   CHANNEL_LIST = 0:156
%
% This is useful for analysing the data in EEGLAB toolbox
% Please see  http://sccn.ucsd.edu/eeglab for more details of the toolbox.
%
% Examples:
% If the path for the sqd-toolbox has been set properly and all of the files present
% (see README), there should be a sample raw sqd-file "testsqd.sqd" for testing. The
% sqd-file "testsqd.sqd" has data for 192 channels, 10 seconds of data and sampled
% at 1000Hz.
% 1. To read in channel location from "testsqd.sqd":
%   locationfile = 'testlocation.xyz';  % Filename for the location file to be created
%   sqdfilename  = 'testsqd.sqd';       % sqd-file to read the locations from
%   writechanloc(sqdfilename,locationfile); % Write location file
%   type(locationfile);                 % To verify
% 2. To read in channel location from "testsqd.sqd" for channels N1,N2,N3:
%   locationfile = 'testlocation.xyz';  % Filename for the location file to be created
%   sqdfilename  = 'testsqd.sqd';       % sqd-file to read the locations from
%   N1 = 1;N2 = 100;N3 = 101;           % Init channels to read
%   writechanloc(sqdfilename,locationfile,'xyz',[N1,N2,N3]); % Write location file
%   type(locationfile);                 % To verify
%
% See also,
% @chanhandle/writechanloc

% Note:
% This function is just a wrapper around @chanhandle/writechanloc
% Error correction is taken care by respective functions

% Init input arguments
switch nargin
case 0
    error('Atleast one input - SQDFILENAME required');
case 1
    locfile = 'testsqd.xyz';
    format = 'xyz';
    channum = 0:156;
case 2
    format = 'xyz';
    channum = 0:156;
case 3
    channum = 0:156;
end;

chant = chanhandle(sqdfile,channum);    % get chanhandle
writechanloc(chant,locfile,format);     % Write file

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
%% Example
% real time control of the KUKA iiwa 7 R 800
% Moving first joint of the robot, using a sinisoidal function
% control is executed at joint velocity level

% An example script, it is used to show how to use the different
% functions of the KUKA Sunrise matlab toolbox

% First run the following script in Matlab
% Then start the client on the KUKA iiwa controller

% Mohammad SAFEEA, 8th-August-2018

close all,clear all;clc;

warning('off')

ip='172.31.1.147'; % The IP of the controller
% start a connection with the server
global t_Kuka;
t_Kuka=net_establishConnection( ip );

if ~exist('t_Kuka','var') || isempty(t_Kuka) || strcmp(t_Kuka.Status,'closed')
  disp('Connection could not be establised, script aborted');
  return;
end

%% Go to initial position

jPos={0,0,0,0,0,0,0};

setBlueOn(t_Kuka); % turn on blue light

relVel=0.25;
movePTPJointSpace( t_Kuka , jPos, relVel); % move to initial configuration

%% Start direct servo in joint space
realTime_startVelControlJoints( t_Kuka );

w=2; % motion constants, frequency rad/sec
A=pi/6; % motion constants, amplitude of motion
counter=0;
jvel={0,0,0,0,0,0,0};
%% Initiate timing variables
dt=0;
tic;
t0=toc; % calculate initial time
%% Control loop
try
    while(dt<(8*pi/w))
        %% Perform velocity command calculation here
        time=toc;
        dt=time-t0;
        jvel{1}=A*w*sin(w*dt);
        counter=counter+1;
        %% Send joint velocties to robot
        sendJointsVelocities( t_Kuka ,jvel);
    end
    tstart=t0;
    tend=time;
    rate=counter/(tend-tstart);
    %% Stop the direct servo motion
    realTime_stopVelControlJoints( t_Kuka );
    fprintf('\nTotal execution time is %f: \n',tend-t0 );
    fprintf('\nThe rate of command update per second is: \n');
    disp(rate);
    fprintf('\n')
    pause(2);
    %% turn off light
    setBlueOff(t_Kuka);
    net_turnOffServer( t_Kuka )
    disp('Direct servo motion completed successfully')
    warning('on')
    return;
catch
    % in case of error turn off the server
    net_turnOffServer( t_Kuka );
    disp('Error during execution the direct servo motion')
    fclose(t_Kuka);
    warning('on')
end

repo: KST-Kuka-Sunrise-Toolbox


text:
% VL_ERFILL  Fill extremal region
%  MEMBERS=VL_ERFILL(I,ER) returns the list MEMBERS of the pixels which
%  belongs to the extremal region represented by the pixel ER.
%
%  The selected region is the one that contains pixel ER and of
%  intensity I(ER).
%
%  I must be of class UINT8 and ER must be a (scalar) index of the
%  region representative point.
%
%  See also: VL_MSER(), VL_HELP().

% Copyright (C) 2007-12 Andrea Vedaldi and Brian Fulkerson.
% All rights reserved.
%
% This file is part of the VLFeat library and is made available under
% the terms of the BSD license (see the COPYING file).

author: yihui-he

repo: panorama

save_path: github-repos/MATLAB/yihui-he-panorama


text:
function [imgLabel, labels, labels_shift] = computeConnectedComponents(img, mode)
%
%
%       [imgLabel, labels, labels_shift] = computeConnectedComponents(img, mode)
%
%
%        Input:
%           -img: an integer grayscale image
%           -mode: 4 or 8 neighbors; this is the connetion type
%
%        Output:
%           -imgLabel: labelled image
%           -labels: all labels in imgLabel
%           -labels_shift:
%
%     Copyright (C) 2011-2020  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/>.
%

r = size(img, 1);
c = size(img, 2);
imgLabel = zeros(r, c);

labels = unique(img);
n = length(labels);
totLabels = 0;

if(~exist('type', 'var'))

end

if(mode ~= 4 || mode ~= 8)
    mode = 4;
end

labels_shift = 0;

for i=1:n

    indx = find(img == labels(i));
    if(~isempty(indx))
        imbBin_i = zeros(r,c); %create a binary image
        imbBin_i(indx) = 1;
        segmented_i = bwlabel(logical(imbBin_i), mode);

        imgLabel = imgLabel + (segmented_i + totLabels);
        totLabels = totLabels + max(segmented_i(:)) + 1;

        labels_shift = [labels_shift, totLabels];
    end
end

end

author: banterle