text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

domain = [-1,1,-1,1,-1,1];
bc = [1,1,1,1,1,1];
nx = 10;
ny = nx;
nz = nx;
r1 = 1; r2 = 0; r3 = 0; x0 = pi/10^4; y0 = 0; z0 = 0;
bm = 1; bp = 1; am = 1; ap = 1;
pde = ConstFun(am,ap,bm,bp,r1,r2,r3,x0,y0,z0);
mesh = genMesh3D(domain, nx, ny, nz);
mesh = enrichMesh3D(mesh,2); % Mesh detail level = 1 (for IFE).
mesh = genIntfMesh3D(mesh,pde.intf);
fem = genNedFEM3D(mesh,bc);

%% 1. use vectorization to generate face and face2elem data structure

% reoder the nodes such that the nodes of interface elements appear at the
% end

IntMeshNode = unique(reshape(mesh.t(mesh.tLoc<0,:),[],1));
MeshNodeID = true(size(mesh.p,1),1);
MeshNodeID(IntMeshNode) = false;
NodeOrderNew = zeros(size(mesh.p,1),1);
NodeOrderNew(MeshNodeID) = 1:sum(MeshNodeID);
NodeOrderNew(~MeshNodeID) = sum(MeshNodeID)+1:size(mesh.p,1);

mesh.t = NodeOrderNew(mesh.t);
mesh.e = NodeOrderNew(mesh.e);
[~,NodeOrderNewSortID] = sort(NodeOrderNew);
mesh.p = mesh.p(NodeOrderNewSortID,:);
mesh.pLoc = mesh.pLoc(NodeOrderNewSortID);

tic

Np = size(mesh.p,1);
node = [mesh.p;mesh.eIntP];
allCutElem = mesh.t(mesh.tLoc<0,:);
isCutElem = (mesh.tLoc<0);
vSign = [mesh.pLoc;zeros(size(mesh.eIntP,1),1)];
isInterfaceNode = false(size(node,1), 1);
isInterfaceNode(allCutElem(:)) = true; % include the vertices of allCutElem
numCutElemV1 = sum(isInterfaceNode);
isInterfaceNode(Np+1:end) = true;    % and the cut points and the aux points
allCutElemReoderTmp = zeros(size(node,1),1);
numCutElemV2 = sum(isInterfaceNode);
allCutElemReoderTmp(isInterfaceNode) = 1:numCutElemV2;
allCutElemReoder = allCutElemReoderTmp(allCutElem);
interfaceNode = node(isInterfaceNode,:);
ntI = -min(mesh.tLoc);
intID = find(mesh.tLoc<0);

tetElem = zeros(12*ntI,4);
tetElemLoc = zeros(12*ntI,1);
idx2cube = zeros(12*ntI,1);
tetcount = 0;

for i = 1:ntI
    tID = intID(i);
    t_e = mesh.t_e(tID,:);
    pLocK = mesh.pLoc(mesh.t(tID,:));
    vert0 = mesh.p(mesh.t(tID,:),:);
    nodeid = [allCutElemReoder(i,:),numCutElemV1-mesh.eLoc(t_e(mesh.eLoc(t_e)<0))'];
    vert2 = vert0(pLocK>0,:); % plus domain
    vert1 = vert0(pLocK<0,:); % minus domain
    intpt0 = mesh.eIntP(-mesh.eLoc(t_e(mesh.eLoc(t_e)<0)),:);
    id1 = [find(pLocK<0)', 5:4+size(intpt0,1)];
    id2 = [find(pLocK>0)', 5:4+size(intpt0,1)];
    p1 = [vert1;intpt0]; DT = delaunayTriangulation(p1); t1 = DT.ConnectivityList;
    p2 = [vert2;intpt0]; DT = delaunayTriangulation(p2); t2 = DT.ConnectivityList;
    tetElem(tetcount+1:tetcount+size(t1,1),:) = nodeid(id1(t1));
    idx2cube(tetcount+1:tetcount+size(t1,1)) = tID;
    tetElemLoc(tetcount+1:tetcount+size(t1,1)) = 1; % inside subdomain
    tetcount = tetcount+size(t1,1);
    tetElem(tetcount+1:tetcount+size(t2,1),:) = nodeid(id2(t2));
    idx2cube(tetcount+1:tetcount+size(t2,1)) = tID;
    tetElemLoc(tetcount+1:tetcount+size(t2,1)) = 2; % outside subdomain
    tetcount = tetcount+size(t2,1);
end
tetElem(tetcount+1:end,:) = [];
tetElemLoc(tetcount+1:end,:) = [];
idx2cube(tetcount+1:end,:) = [];
[tetElem, ortidx, volume] = fixorder3(interfaceNode, tetElem);
% there are some tets which are coplane, need to get rid of them
idPlane = (volume<=10^(-12));
tetElem(idPlane,:) = [];
volume(idPlane,:) = [];
tetElemLoc(idPlane,:) = [];
idx2cube(idPlane,:) = [];
localidx2globalidx = find(isInterfaceNode); % tetElem points to interfaceNode
tetElem = localidx2globalidx(tetElem);
% reorder tetElem such that the index of its nodes are local, i.e.,
% starting from 1,2,...

% there are some tets which are on the interface, but contained in the surrounding tets
% need to get rid of them (but the following algorithm may not be robust)
tetSign = vSign(tetElem);
idSliver = (sum(abs(tetSign),2)==0);
tetElem(idSliver,:) = [];
volume(idSliver,:) = [];
tetElemLoc(idSliver,:) = [];
idx2cube(idSliver,:) = [];
PolyVolume = accumarray([idx2cube,tetElemLoc],volume);
PolyVolume = PolyVolume(mesh.tLoc<0,:);
h = (PolyVolume(:,1) + PolyVolume(:,2)).^(1/3);
% the first component contains the volume of the polyhedron on the
% subdomain 1, the first component contains the volume of the polyhedron
% on the subdomain 2

toc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Get triangular faces for interrior elements and interface

tic
localFace = [2 3 4; 1 4 3; 1 2 4; 1 3 2];
NT = size(tetElem,1);
tface = zeros(4*NT, 3);
tface2elem = zeros(4*NT, 1);
iface = zeros(4*NT, 3);
iface2elem = zeros(4*NT, 1);
% find the interior tet elements
isIntTet = min(vSign(tetElem),[], 2) == -1; % can not be == -1 as there is sliver
intTet = tetElem(isIntTet,:);
% find the corresponding cube indices
intIdx2cube = idx2cube(isIntTet);
% find triangular interface
T = auxstructure3(intTet);
neighbor = T.neighbor; % if a face is on the boundary, then the neighbor element index is itself
clear T;
tmp = (1:size(intTet, 1))';
ct = 0;
ci = 0;
for i = 1:4
    face = intTet(:, localFace(i,:));
    % find the triangle faces of polyhedron
    % 1. face and its neighbor associated to different cubes, and
    % 2. face is not on the boundary of all cut elements (which are squares)
    isPolyTriFace = ((neighbor(:, i) == tmp) | (intIdx2cube ~= intIdx2cube(neighbor(:,i)))) &...
        (sum(abs(vSign(face)), 2) > 0);% & (sum(abs(vSign(face)), 2) > 0);
    c2 = ct + sum(isPolyTriFace);
    tface((ct+1):c2,:) = face(isPolyTriFace,:);
    tface2elem((ct+1):c2,:) = intIdx2cube(isPolyTriFace); % the indices of the polyhedron
    ct = c2;
    % note that only interior elements are treated

    % find the triangle faces on interface with normal points to exterior
    % 1. face is on the boundary and
    % 2. all vertices are on the interface
    isInterface = (sum(abs(vSign(face)), 2) == 0);

    % add to interface
    c4 = ci + sum(isInterface);
    iface((ci+1):c4,:) = face(isInterface,:); % the interface tri faces.
    iface2elem((ci+1):c4,:) = intIdx2cube(isInterface); % the indices of the polyhedron
    ci = c4;
end
iface((ci+1):end,:) = [];
iface2elem((ci+1):end,:) = [];
c2old = c2;
% plot to check
%trisurf(tface(1:c2,:),node(:,1),node(:,2),node(:,3))
%trisurf(iface(1:ci,:),node(:,1),node(:,2),node(:,3))

% Find the triangular faces for exterior elements
extTet = tetElem(~isIntTet, :);
extIdx2cube = idx2cube(~isIntTet);

T = auxstructure3(extTet);
neighbor = T.neighbor;
clear T;
tmp = (1:size(extTet,1))';
for i = 1:4
    face = extTet(:, localFace(i,:));
    % find the triangle faces of polyhedron
    % 1. face and its neighbor associated to different cubes, and
    % 2. face is not on the boundary of all cut elements (which are squares)
    isPolyTriFace = ((neighbor(:, i) == tmp) | (extIdx2cube ~= extIdx2cube(neighbor(:,i)))) &...
        (sum(abs(vSign(face)), 2) > 0);
    c2 = ct + sum(isPolyTriFace);
    tface((ct+1):c2,:) = face(isPolyTriFace,:);
    tface2elem((ct+1):c2,:) = extIdx2cube(isPolyTriFace);
    % index of exterior polyhedron is append to the end of elem
    ct = c2;
end
toc

tface((ct+1):end,:) = [];    % remove empty meomory
tface2elem((ct+1):end) = [];
face2elemLoc = [ones(c2old,1);2*ones(c2-c2old,1)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% generate edge structure of elements
% reorder the DoFs: put the DoFs of the edges of non-interface elements (sans those of interface elements)
% at the begining, then cancate the new edges including interface edges and
% non-interface edges

IntNodeID = [unique(reshape(mesh.t(mesh.tLoc<0,:),[],1));...
    (size(mesh.p,1)+1:size(mesh.p,1) + size(mesh.eIntP,1))'];
AllNodeID = zeros(size(node,1),1);
AllNodeID(IntNodeID) = 1:length(IntNodeID);
tfaceNew = AllNodeID(tface);
[face2EDoF,Newedge,face2EDoFSign] = dof3edge(tfaceNew);
% IMPORTANT: the Newedge contains those inside elements (4 intersection points)
% the edge orientation of dof3edge is from 1->2, 2->3, 3->1 of nodes of tfacenew
% face2EDoFSign assigns a sign such that orientation is from small to large
Newedge = IntNodeID(Newedge); % transfer the node index of edges to the global ones
% the new index of non-interface edges of interface elements (starting from 1)
NewedgeNintID = find(abs(sum(vSign(Newedge),2))==2);

OldedgeID = unique(reshape(mesh.t_e(mesh.tLoc<0,:),[],1));
Oldedge = mesh.e(OldedgeID,:);
OldEtmp = abs(sum(vSign(Oldedge),2))==2;
OldedgeNintID = OldedgeID(OldEtmp);
% the old index of non-interface edges of interface elements
OldedgeNint = Oldedge(OldEtmp,:);
[OldedgeNintSort,OldedgeNintSortID] = sortrows(OldedgeNint);
% ~ should be exactly the same as Newedge(NewedgeNintID)
% therefore OldedgeNintSortID should be the new index of all the
% non-interface edges of interface elements
BasicEdge = ones(size(mesh.e,1),1); BasicEdge(OldedgeID) = false;
BasicEdge = logical(BasicEdge);
% basically the non-interface edges sans edges of interface elements
BasicEdgeNum = sum(BasicEdge);
TotalOldEdge = zeros(size(mesh.e,1),1);
TotalOldEdge(BasicEdge) = 1:BasicEdgeNum;
TotalOldEdge(OldedgeNintID(OldedgeNintSortID)) = BasicEdgeNum + NewedgeNintID;
% map OldedgeNintSortID to global ones and put them at the location of
% the original non-interface edges
% NOTE: at the location of interface edges, TotalOldEdge still contains zero

% generate new global edge DoF structure
NumEdge = size(mesh.e,1) - size(Oldedge,1) + size(Newedge,1);
% NumEdge should equal size(mesh.e,1) + sum(mesh.eLoc<0) + 2*sum(mesh.fLoc<0)
gdof = zeros(NumEdge,2);
gdof(1:BasicEdgeNum,:) = mesh.e(BasicEdge,:);
gdof(BasicEdgeNum+1:end,:) = Newedge;
NintElem = mesh.tLoc>0;
g2ldofNint = TotalOldEdge(mesh.t_e(NintElem,:));
face2EDoF = face2EDoF + BasicEdgeNum;
% each row contains the DoFs of each element

% generate orientation structure for each edge
t_e_orit = zeros(size(g2ldofNint));
e_ind = [1,2; 1,3; 1,4; 2,3; 2,4; 3,4];
for i = 1:6
    d = mesh.p(mesh.t(NintElem,e_ind(i,2)),:) - mesh.p(mesh.t(NintElem,e_ind(i,1)),:);
    %d_crect = mesh.p(mesh.e(mesh.t_e(NintElem,i),2),:) - mesh.p(mesh.e(mesh.t_e(NintElem,i),1),:);
    d_crect = node(gdof(g2ldofNint(:,i),2),:) - node(gdof(g2ldofNint(:,i),1),:);
    t_e_orit(:,i) = sign(sum(d.*d_crect,2));
end

meshI.tface = tface;
meshI.tface2elem = tface2elem;
meshI.iface = iface;
meshI.iface2elem = iface2elem;
meshI.face2elemLoc = face2elemLoc;
meshI.PolyVolume = PolyVolume;
% meshI.vSign = vSign;
meshI.node = node;
meshI.tetElem = tetElem;
meshI.tetElemLoc = tetElemLoc;
meshI.idx2cube = idx2cube;
meshI.tetVolume = volume;
meshI.face2EDoF = face2EDoF;
meshI.face2EDoFSign = face2EDoFSign;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2.use face, face2elem and face2EDoF data structure to compute ife functions

face2elem = meshI.tface2elem;
tface = meshI.tface;
N = size(node, 1); NEdof = NumEdge;
NF = size(tface,1); NC = size(mesh.t, 1);

NP = max(face2elem);
isCutPoly = false(NP, 1);
isCutPoly(face2elem) = true;
% this is the same as:
% isCutPoly = (mesh.tLoc<0);
% isCutPoly = isCutPoly(1:max(face2elem));
cutPolyIdx = zeros(NP, 1);
isExtPoly = false(NP,1);
isExtPoly(NC+1:NP) = true;

NP = sum(isCutPoly); % this new NP is max(-mesh.tLoc(mesh.tLoc<0));
cutPolyIdx(isCutPoly) = 1:NP;
% this is the same as:
% cutPolyIdx = zeros(size(mesh.tLoc,1), 1);
% cutPolyIdx(mesh.tLoc<0) = -mesh.tLoc(mesh.tLoc<0);
% cutPolyIdx = cutPolyIdxNew(1:max(face2elem));
face2elem = cutPolyIdx(face2elem);
% the old face2elem is the global index of interface elements
% the new one is only the index of interface elements
isExtPoly = isExtPoly(isCutPoly);
%%%clear isCutPoly cutPolyIdx

% prepare new data structure to compute projections
poly2node = sparse(face2elem(:)*ones(1,3), tface(:), 1, NP,N);
poly2node = (poly2node > 0);
NV = poly2node*ones(N, 1);

poly2edge = sparse(face2elem(:)*ones(1,3), face2EDoF(:), 1, NP,NEdof);
poly2edge = (poly2edge > 0);
NE = poly2edge*ones(NEdof, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% generate local IFE basis functions which are not associated with any DoFs
% IFEbasis(:,i,:,:) for minus subdomain (i=1) or plus domain (i=2)
iface = meshI.iface;
iface2elem = meshI.iface2elem;
[iface2elemReduce, uniID] = unique(iface2elem);
% length(uniID) should be NP
ifaceReduce = iface(uniID,:);
[IntFnormal,IntFarea,unitNormal,unittgt1,unittgt2] = facenormaltgt(node,ifaceReduce);
[Fnormal,Farea,FunitNormal,Funittgt1,Funittgt2] = facenormaltgt(node,tface);
% IFEbasis contains the associated IFE function on each face
% IFEbasisA is for curl term and IFEbasisB is for u term
IFEbasisA = zeros(size(face2elem,1),3,3);
% for two pieces of an ife function, only the normal vector is different
% on the minus subdomain (1) ife = n; on the plus subdomain (2) ife = n*bm/bp
IFEbasisA(:,1,:) = unitNormal(face2elem,:);
IFEbasisA(:,2,:) = unittgt1(face2elem,:);
IFEbasisA(:,3,:) = unittgt2(face2elem,:);
IFEbasisB = IFEbasisA;
piece1tmp = (face2elemLoc==1);
IFEbasisA(~piece1tmp,2,:) = IFEbasisA(~piece1tmp,2,:)*am/ap;
IFEbasisA(~piece1tmp,3,:) = IFEbasisA(~piece1tmp,3,:)*am/ap;
IFEbasisB(~piece1tmp,1,:) = IFEbasisB(~piece1tmp,1,:)*bm/bp;
% Xmf the centroid of a triangular face on the interface of each interface
% element but assigned to each face
% Xf the centroid of each triangular face
Xmf = (node(ifaceReduce(face2elem,1),:) + node(ifaceReduce(face2elem,2),:) +...
    node(ifaceReduce(face2elem,3),:))/3; % the Xm point on the approximate interface plane
Xf = (node(tface(:,1),:) + node(tface(:,2),:) +...
    node(tface(:,3),:))/3;

%% generate boundary integral for computing projections
% for the curl term, need to compute surface curl (rot) on each face
BIFEA =  zeros(length(face2elem),3);
BIFEB =  zeros(length(face2elem),3,3);
% BIFEA contains  (alpha*(grad vi).(xf -xm))for i=1,2,3 (vi is the test function)
BIFEA(:,1) = sum(squeeze(IFEbasisA(:,1,:)).*(Xf - Xmf),2);
BIFEA(:,2) = sum(squeeze(IFEbasisA(:,2,:)).*(Xf - Xmf),2);
BIFEA(:,3) = sum(squeeze(IFEbasisA(:,3,:)).*(Xf - Xmf),2);
% no need multiply area for this case as computing curl generates |F|^(-1)
BIFEA(piece1tmp,:) = am*BIFEA(piece1tmp,:);
BIFEA(~piece1tmp,:) = ap*BIFEA(~piece1tmp,:);
% BIFEB contains  wh=(beta*(curl vi)times(xf -xm))times n for i=1,2,3
% position to match the structure of mytimes
BIFEB(:,:,1) = -cross(cross(squeeze(IFEbasisB(:,1,:)),(Xf - Xmf)),FunitNormal)/2;
BIFEB(:,:,2) = -cross(cross(squeeze(IFEbasisB(:,2,:)),(Xf - Xmf)),FunitNormal)/2;
BIFEB(:,:,3) = -cross(cross(squeeze(IFEbasisB(:,3,:)),(Xf - Xmf)),FunitNormal)/2;
% devided by 2 is because curl(c times x) = 2c
BIFEB(piece1tmp,:) = bm*BIFEB(piece1tmp,:);
BIFEB(~piece1tmp,:) = bp*BIFEB(~piece1tmp,:);
% project the 3D vector onto each face
FNT = zeros(length(face2elem),3,3);
FNT(:,1,:) = Funittgt1;
FNT(:,2,:) = Funittgt2;
FNT(:,3,:) = FunitNormal;
BIFEB = mytimes(FNT,BIFEB);
BIFEB = BIFEB(:,1:2,:);
%%%%%%%%%%%%%%%%
% mutiply the inverse the matrix on the left hand side of each interface
% element for computing projection
% this is a blocal diagonal matrix
% for BIFEA (curl):
VdiagTmp1 = (am*PolyVolume(:,1) + ap*PolyVolume(:,2)).^(-1);
VdiagTmp2 = (am*PolyVolume(:,1) + am^2/ap*PolyVolume(:,2)).^(-1);
VdiagTmp3 = (am*PolyVolume(:,1) + am^2/ap*PolyVolume(:,2)).^(-1);
Vdiag1 = VdiagTmp1(face2elem);
Vdiag2 = VdiagTmp2(face2elem);
Vdiag3 = VdiagTmp3(face2elem);
Mdiag = zeros(size(Vdiag1,1),3,3);
Mdiag(:,1,1) = Vdiag1; Mdiag(:,2,2) = Vdiag2; Mdiag(:,3,3) = Vdiag3;
BIFEA = mytimes(Mdiag,BIFEA);
%%%%%%%%%%%%%%%%
% for BIFEB (u):
VdiagTmp1 = (bm*PolyVolume(:,1) + bm^2/bp*PolyVolume(:,2)).^(-1);
VdiagTmp2 = (bm*PolyVolume(:,1) + bp*PolyVolume(:,2)).^(-1);
VdiagTmp3 = (bm*PolyVolume(:,1) + bp*PolyVolume(:,2)).^(-1);
Vdiag1 = VdiagTmp1(face2elem);
Vdiag2 = VdiagTmp2(face2elem);
Vdiag3 = VdiagTmp3(face2elem);
Mdiag = zeros(size(Vdiag1,1),3,3);
Mdiag(:,1,1) = Vdiag1; Mdiag(:,2,2) = Vdiag2; Mdiag(:,3,3) = Vdiag3;
%%%%%%%%%%%%%%%%
% compute [-(ym-y1),(xm-x1)]
VecXm2Pt = (Xf - node(tface(:,1),:))/2;
VecXm2Pt = mytimes(FNT,VecXm2Pt);
VecXm2Pt = VecXm2Pt(:,[2,1]); VecXm2Pt(:,1) = -VecXm2Pt(:,1);
% compute [-(ym-y1),(xm-x1)].wh
BIFEB1 = zeros(size(VecXm2Pt));
for i = 1:3
    BIFEB1(:,i) = sum(VecXm2Pt.*squeeze(BIFEB(:,:,i)),2);
end
BIFEB1 = mytimes(Mdiag,BIFEB1);
% compute inv([t1;t2]) where t1 and t2 are for the 1st and 3rd edges
MatT = zeros(length(face2elem),3,2);
MatT(:,:,1) = node(tface(:,2),:) - node(tface(:,1),:);
MatT(:,:,2) = node(tface(:,1),:) - node(tface(:,3),:);
MatT = mytimes(FNT,MatT);
MatT = MatT(:,1:2,:); MatTInv = zeros(size(MatT));
MatTInv(:,1,1) = MatT(:,2,2)./(-2*Farea); % generate inverse matrix
MatTInv(:,2,2) = MatT(:,1,1)./(-2*Farea);
MatTInv(:,1,2) = -MatT(:,1,2)./(-2*Farea);
MatTInv(:,2,1) = -MatT(:,2,1)./(-2*Farea);
clear MatT
BIFEB2 = mytimes(MatTInv,BIFEB,[2,2]);  % T^(-1)w_h
BIFEB2 = mytimes(Mdiag,permute(BIFEB2,[1,3,2]));
BIFEB2 = permute(BIFEB2,[1,3,2]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute projections to IFE spaces and form the matrices
face2EDoF = meshI.face2EDoF;
face2EDoFSign = meshI.face2EDoFSign;
IFNT = zeros(length(ifaceReduce),3,3);
IFNT(:,:,1) = unitNormal;
IFNT(:,:,2) = unittgt1;
IFNT(:,:,3) = unittgt2;

nnz = sum(NE.^2);
iiP = zeros(nnz, 1);
jjP = zeros(nnz, 1);
ssKC = zeros(nnz, 1);
ssMU = zeros(nnz, 1);
indexP = 0;
iiF = zeros(nnz, 1);
jjF = zeros(nnz, 1);
ssSC = zeros(nnz, 1);
ssSU = zeros(nnz, 1);
indexF = 0;
une = unique(NE);
hF = h(face2elem);
b = zeros(NEdof, 1);

for kk = 1:length(une) % group polys according to their # of nodes

    tic
    %% generate IFE functions including the IFE basis functions on each interface
    %% elements which are not associated with any DoFs, and projections of gradients
    %% and the IFE functions with matching face average
    ne = une(kk);
    % some data related to element
    isCurrentFace = (NE(face2elem) == ne);
    CurrentFace = tface(isCurrentFace,:);
    CurrentFaceE = face2EDoF(isCurrentFace,:);
    CurrentFaceS = face2EDoFSign(isCurrentFace,:);
    CNF = size(CurrentFaceE,1);
    Currentfaceloc = face2elemLoc(isCurrentFace);
    Currentface2elem = face2elem(isCurrentFace);
    CurrentFnormal = FunitNormal(isCurrentFace,:);
    FareaCurrent = Farea(isCurrentFace);
    % some data related to face
    isCurrentElem = false(NP, 1);
    isCurrentElem(face2elem(isCurrentFace)) = true;
    CurrentPolyVolume = PolyVolume(isCurrentElem,:);
    IFNTCurrent = IFNT(isCurrentElem,:,:);
    % Current elem to node matrix
    currentPoly2edge = poly2edge(isCurrentElem, :);
    CNP = size(currentPoly2edge, 1);
    currentPolyLocalIdx = zeros(NP, 1);
    currentPolyLocalIdx(isCurrentElem) = 1:CNP;

    % currentElem(i,:) contains the edge index of the (current) i-th element
    [I, ~] = find(currentPoly2edge');
    currentElem = reshape(I, ne, [])';
    % localIdx(i,j) contains the (current) i-th element having the
    % node(dof=1,2,3,4) on the j-th location (j is the global node index)
    localIdx = sparse(repmat((1:CNP)', 1, ne), currentElem, ones(CNP, 1)*(1:ne), CNP, NEdof);
    clear currentPoly2edge

    %% Deal with current triangle face cases
    %% (1) for projection of curl u
    subs1 = currentPolyLocalIdx(face2elem(isCurrentFace));
    subs2 = [ones(CNF, 1); 2*ones(CNF, 1); 3*ones(CNF, 1)];
    val = BIFEA(isCurrentFace,:);

    % compute the projection of gradients
    GPA1 = zeros(CNP,3,ne);
    % coefficient of n t1 and t2 for the subelement on the subdomain 1
    for m = 1:3
        subs3 = full(localIdx((CurrentFaceE(:, m) - 1)*CNP + subs1));
        SignCurrentE = repmat(CurrentFaceS(:,m),3,1);
        % subs1: new order (range in 1:CNP) of the element index
        % subs2: three components in the vector
        % subs3: the local index (DoF) (w.r.t. element) of the m-th edge of this face
        GPA1 = GPA1 + accumarray([repmat(subs1,3,1), subs2, repmat(subs3, 3, 1)],...
            SignCurrentE.*val(:), [CNP, 3, ne]);
    end
    GPA2 = GPA1; GPA2(:,2,:) = am/ap*GPA2(:,2,:); GPA2(:,3,:) = am/ap*GPA2(:,3,:);
    GPA1 = mytimes(IFNTCurrent,GPA1);
    GPA2 = mytimes(IFNTCurrent,GPA2);
    % the updated GP contains three component values of the projection vector
    % generate IFE vectors on each face
    % clear currentPolyLocalIdx localIdx
    % clear subs1 subs2 subs3

    % generate the data for computing stabilization matrices
    CurrentFace2DoF = zeros(NP,ne);
    CurrentFace2DoF(isCurrentElem,:) = currentElem;
    CurrentFace2DoF = CurrentFace2DoF(Currentface2elem,:);
    % CurrentFace2DoF contains the DoFs of the element associated with each face

    % compute curl uh.n on each face (uh is an IFE function)
    GPA1tmp = zeros(NP,3,ne); GPA2tmp = zeros(NP,3,ne);
    GPA1tmp(isCurrentElem,:,:) = GPA1; GPA2tmp(isCurrentElem,:,:) = GPA2;
    GAface = zeros(CNF,3,ne);
    GAface(Currentfaceloc==1,:,:) = GPA1tmp(Currentface2elem(Currentfaceloc==1),:,:);
    GAface(Currentfaceloc==2,:,:) = GPA2tmp(Currentface2elem(Currentfaceloc==2),:,:);
    clear GPA1tmp GPA2tmp
    GAfaceP = zeros(size(GAface,1),ne);
    for n = 1:ne
        GAfaceP(:,n) =  sum(squeeze(GAface(:,:,n)).*CurrentFnormal,2);
    end

    % compute rot_F uh on each face according to DoFs
    TrueCurlTmp = CurrentFaceS./FareaCurrent;
    TrueCurl = zeros(CNF,ne);
    for i = 1:ne
        for j = 1:3
            IDij = (CurrentFaceE(:,j) == CurrentFace2DoF(:,i));
            TrueCurl(IDij,i) = TrueCurlTmp(IDij,j);
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% (2) for projection of u
    val = BIFEB1(isCurrentFace,:);
    GPB1 = zeros(CNP,3,ne);
    % coefficient of n t1 and t2 for the subelement on the subdomain 1
    for m = 1:3
        subs3 = full(localIdx((CurrentFaceE(:, m) - 1)*CNP + subs1));
        SignCurrentE = repmat(CurrentFaceS(:,m),3,1);
        % subs1: new order (range in 1:CNP) of the element index
        % subs2: three components in the vector
        % subs3: the local index (DoF) (w.r.t. element) of the m-th edge of this face
        GPB1 = GPB1 + accumarray([repmat(subs1,3,1), subs2, repmat(subs3, 3, 1)],...
            SignCurrentE.*val(:), [CNP, 3, ne]);
    end

    edgetmp = [1,3];
    for mm = 1:2
        val = squeeze(BIFEB2(isCurrentFace,mm,:));
        m = edgetmp(mm);
        subs3 = full(localIdx((CurrentFaceE(:, m) - 1)*CNP + subs1));
        SignCurrentE = repmat(CurrentFaceS(:,m),3,1);
        Areatmp = repmat(FareaCurrent,3,1);
        % subs1: new order (range in 1:CNP) of the element index
        % subs2: three components in the vector
        % subs3: the local index (DoF) (w.r.t. element) of the m-th edge of this face
        GPB1 = GPB1 + accumarray([repmat(subs1,3,1), subs2, repmat(subs3, 3, 1)],...
            SignCurrentE.*val(:).*Areatmp, [CNP, 3, ne]);
    end
    GPB2 = GPB1; GPB2(:,2,:) = am/ap*GPB2(:,2,:); GPB2(:,3,:) = am/ap*GPB2(:,3,:);
    GPB1 = mytimes(IFNTCurrent,GPB1);
    GPB2 = mytimes(IFNTCurrent,GPB2);
    % the updated GP contains three component values of the projection vector
    % generate IFE vectors on each face
    clear currentPolyLocalIdx localIdx
    clear subs1 subs2 subs3

    % compute uh^(\tau) (projection onto each face)
    GPB1tmp = zeros(NP,3,ne); GPB2tmp = zeros(NP,3,ne);
    GPB1tmp(isCurrentElem,:,:) = GPB1; GPB2tmp(isCurrentElem,:,:) = GPB2;
    GBface = zeros(CNF,3,ne);
    GBface(Currentfaceloc==1,:,:) = GPB1tmp(Currentface2elem(Currentfaceloc==1),:,:);
    GBface(Currentfaceloc==2,:,:) = GPB2tmp(Currentface2elem(Currentfaceloc==2),:,:);
    clear GPA1tmp GPA2tmp
    GBfaceP = zeros(size(GBface));
    for n = 1:ne
        GBfaceP(:,:,n) = GBface(:,:,n) - sum(GBface(:,:,n).*CurrentFnormal,2).*CurrentFnormal;
    end
    clear GBface
    GBfaceP = mytimes(FNT(isCurrentFace,:,:),GBfaceP);
    GBfaceP = GBfaceP(:,1:2,:);

    % compute u^(tau)(xm) uh on each face according to DoFs
    TrueUtTmp = zeros(CNF,2,ne);
    TrueUt = zeros(CNF,2,ne);
    MatTInvCurrent = MatTInv(isCurrentFace,:,:);
    for i = 1:ne
        TrueUt(:,:,i) = TrueUt(:,:,i) + TrueCurl(:,i).*VecXm2Pt(isCurrentFace,:);
        for j = 1:2
            jj = edgetmp(j);
            IDij = (CurrentFaceE(:,jj) == CurrentFace2DoF(:,i));
            TrueUtTmp(IDij,:,i) = MatTInvCurrent(IDij,j,:).*CurrentFaceS(IDij,jj);
        end
    end
    TrueUt = TrueUt + TrueUtTmp;
    clear TrueUtTmp

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% generate stiffness and mass matrices
    for n = 1:ne
        for m = 1:ne
            iiP(indexP+1:indexP + CNP) = currentElem(:, n);
            jjP(indexP+1:indexP + CNP) = currentElem(:, m);
            ssKC(indexP+1:indexP + CNP) = am*dot(GPA1(:,:,n), GPA1(:,:,m),2).*CurrentPolyVolume(:,1)+...
                ap*dot(GPA2(:,:,n), GPA2(:,:,m),2).*CurrentPolyVolume(:,2);
            ssMU(indexP+1:indexP + CNP) = bm*dot(GPB1(:,:,n), GPB1(:,:,m),2).*CurrentPolyVolume(:,1)+...
                bp*dot(GPB2(:,:,n), GPB2(:,:,m),2).*CurrentPolyVolume(:,2);
            indexP = indexP + CNP;

            iiF(indexF+1:indexF + CNF) = CurrentFace2DoF(:, n);
            jjF(indexF+1:indexF + CNF) = CurrentFace2DoF(:, m);
            ssSC(indexF+1:indexF + CNF) = dot((TrueCurl(:,m) - GAfaceP(:,m)),...
                (TrueCurl(:,n) - GAfaceP(:,n)),2).*...
                Farea(isCurrentFace).*hF(isCurrentFace);
            ssSU(indexF+1:indexF + CNF) = ( dot((TrueUt(:,1,m) - GBfaceP(:,1,m)),...
                (TrueUt(:,1,n) - GBfaceP(:,1,n)),2) + dot((TrueUt(:,2,m) - GBfaceP(:,2,m)),...
                (TrueUt(:,2,n) - GBfaceP(:,2,n)),2) ).*Farea(isCurrentFace).*hF(isCurrentFace);
            indexF = indexF + CNF;
        end
    end

    %% form the right-hand side vector
    idx2cubeNew = -mesh.tLoc(idx2cube);
    TetID = isCurrentElem(idx2cubeNew);
    currentvolume = volume(TetID);
    %%% DoFs of elements for tet
    currentTetDoF = zeros(size(isCurrentElem,1),ne);
    currentTetDoF(isCurrentElem,:) = currentElem;
    currentTetDoF = currentTetDoF(idx2cubeNew(TetID),:);
    %%% quadrature points
    X1 = node(tetElem(TetID,1),:);
    X2 = node(tetElem(TetID,2),:);
    X3 = node(tetElem(TetID,3),:);
    X4 = node(tetElem(TetID,4),:);
    ng = 1;
    [gx, gy, gz] = gaussPtetra(X1, X2, X3, X4, ng);
    gw = gaussWtetra(ng);
    Xm = (node(ifaceReduce(isCurrentElem,1),:) + node(ifaceReduce(isCurrentElem,2),:) +...
        node(ifaceReduce(isCurrentElem,3),:))/3;
    Xmtet = zeros(size(isCurrentElem,1),3);
    Xmtet(isCurrentElem,:) = Xm;
    Xmtet = Xmtet(idx2cubeNew(TetID),:);
    %%% projection of u
    Bas1 = zeros(size(isCurrentElem,1),3,ne); Bas2 = Bas1;
    Bas1(isCurrentElem,:,:) = GPB1;
    Bas2(isCurrentElem,:,:) = GPB2;
    Bas1 = Bas1(idx2cubeNew(TetID),:,:); Bas2 = Bas2(idx2cubeNew(TetID),:,:);
    Bas = zeros(sum(TetID),3,ne);
    piecetmp = tetElemLoc(TetID);
    Bas(piecetmp==1,:,:) = Bas1(piecetmp==1,:,:);
    Bas(piecetmp==2,:,:) = Bas2(piecetmp==2,:,:);

    ft1 = pde.f1(gx,gy,gz);
    ft2 = pde.f2(gx,gy,gz);
    ft3 = pde.f3(gx,gy,gz);
    Currentb = zeros(size(ft1));
    cnt = size(ft1,1);
    bii = ne*cnt;
    count = 0 ;
    for i = 1:ne
        Basi = squeeze(Bas(:,:,i));
        Currentb(:,i) = sum(gw*(ft1.*Basi(:,1)+ft2.*Basi(:,2)+ft3.*Basi(:,3)),2).*...
            currentvolume;
        bii(count+1:count+cnt) = currentTetDoF(:,i);
        count = count+cnt;
    end
    b = b + sparse(bii,1,reshape(Currentb,[],1),NEdof,1);

    toc

end

Kcurl = sparse(iiP, jjP, ssKC, NEdof, NEdof);
Scurl = sparse(iiF, jjF, ssSC, NEdof, NEdof);
MU = sparse(iiP, jjP, ssMU, NEdof, NEdof);
SU = sparse(iiF, jjF, ssSU, NEdof, NEdof);

MatI = Kcurl + Scurl + MU +SU;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% generate matrix on non-interface elements (for test)
feEvalBas1 = @EvalNed1Bas3D;
feEvalBas2 = @EvalNed1Bas3D;

dof1 = 6; dof2 = 6; nloc = dof1*dof2;
ntID = find(mesh.tLoc > 0); ntN = length(ntID);

%% 1. Matrix on noninterface elements
MatInd = 1;
AN = fem.area(ntID);
gxN = fem.gx(ntID,:); gyN = fem.gy(ntID,:); gzN = fem.gz(ntID,:); gw = fem.gw;
XN = zeros(nloc*ntN, 1);

coefN = feval(pde.A,gxN,gyN,gzN);
Ibasx = cell(dof1,1); Ibasy = cell(dof1,1); Ibasz = cell(dof1,1);
Jbasx = cell(dof2,1); Jbasy = cell(dof2,1); Jbasz = cell(dof2,1);

for i = 1:dof1
    Ibasx{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 1).*t_e_orit(:,i);
    Ibasy{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 2).*t_e_orit(:,i);
    Ibasz{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 3).*t_e_orit(:,i);
end
for j = 1:dof2
    Jbasx{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 1).*t_e_orit(:,j);
    Jbasy{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 2).*t_e_orit(:,j);
    Jbasz{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 3).*t_e_orit(:,j);
end

IN = reshape(repmat(g2ldofNint(:,1:6),6,1),nloc*ntN,1);
JN = repmat(reshape(g2ldofNint(:,1:6),dof2*ntN,1),6,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XN(ind+1:ind+ntN) = AN.*(sum(((Ibasx{i}.*(coefN.*Jbasx{j})).*gw'),2) + ...
            sum(((Ibasy{i}.*(coefN.*Jbasy{j})).*gw'),2) + ...
            sum(((Ibasz{i}.*(coefN.*Jbasz{j})).*gw'),2));
        ind = ind + ntN;
    end
end
ID = find(XN~=0);
SN = sparse(IN(ID),JN(ID),XN(ID),NEdof,NEdof);

StiffCurl = SN + Kcurl + Scurl;

%%%%%%

MatInd = 0;
AN = fem.area(ntID);
gxN = fem.gx(ntID,:); gyN = fem.gy(ntID,:); gzN = fem.gz(ntID,:); gw = fem.gw;
XN = zeros(nloc*ntN, 1);

coefN = feval(pde.A,gxN,gyN,gzN);
Ibasx = cell(dof1,1); Ibasy = cell(dof1,1); Ibasz = cell(dof1,1);
Jbasx = cell(dof2,1); Jbasy = cell(dof2,1); Jbasz = cell(dof2,1);

for i = 1:dof1
    Ibasx{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 1).*t_e_orit(:,i);
    Ibasy{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 2).*t_e_orit(:,i);
    Ibasz{i} = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, MatInd, 3).*t_e_orit(:,i);
end
for j = 1:dof2
    Jbasx{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 1).*t_e_orit(:,j);
    Jbasy{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 2).*t_e_orit(:,j);
    Jbasz{j} = feEvalBas2(fem.bas, ntID, gxN, gyN, gzN, j, MatInd, 3).*t_e_orit(:,j);
end

IN = reshape(repmat(g2ldofNint(:,1:6),6,1),nloc*ntN,1);
JN = repmat(reshape(g2ldofNint(:,1:6),dof2*ntN,1),6,1);
ind = 0;
for i = 1:dof1
    for j = 1:dof2
        XN(ind+1:ind+ntN) = AN.*(sum(((Ibasx{i}.*(coefN.*Jbasx{j})).*gw'),2) + ...
            sum(((Ibasy{i}.*(coefN.*Jbasy{j})).*gw'),2) + ...
            sum(((Ibasz{i}.*(coefN.*Jbasz{j})).*gw'),2));
        ind = ind + ntN;
    end
end
ID = find(XN~=0);
MN = sparse(IN(ID),JN(ID),XN(ID),NEdof,NEdof);

MassU = MN + MU + SU;

%%%%%%%%%%

dof = 6;  nloc = dof;
X = zeros(nloc*ntN, 1);

fN1 = feval(pde.exactu1,gxN,gyN,gzN);
fN2 = feval(pde.exactu2,gxN,gyN,gzN);
fN3 = feval(pde.exactu3,gxN,gyN,gzN);
ind = 0;
I = reshape(g2ldofNint(:,1:6),nloc*ntN,1);
for i = 1:dof
    ibas1 = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, 0, 1);
    ibas2 = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, 0, 2);
    ibas3 = feEvalBas1(fem.bas, ntID, gxN, gyN, gzN, i, 0, 3);
    X(ind+1:ind+ntN) = AN.*sum((ibas1.*fN1+ibas2.*fN2+ibas3.*fN3).*gw',2).*fem.t_e_orit(ntID,i);
    ind = ind + ntN;
end
rhsN = sparse(I,1,X,NEdof,1);

rhs = rhsN + b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[gew,gex,gey,gez] = gaussPedge(node(gdof(:,1),:),node(gdof(:,2),:),1);

tu1 = sum(feval(pde.exactu1,gex,gey,gez).*gew,2);
tu2 = sum(feval(pde.exactu2,gex,gey,gez).*gew,2);
tu3 = sum(feval(pde.exactu3,gex,gey,gez).*gew,2);

tgt = node(gdof(:,2),:) - node(gdof(:,1),:);
tgt = tgt./sum(tgt.^2,2).^(1/2);
tu = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);

Atotal = StiffCurl + MassU;
%b = ones(size(MassU,1),1);
x0 = zeros(size(b));

bcind = fem.bcind;
[bc,mapper] = boundaryEdge3D(node,gdof,bcind);
bdidx = zeros(NEdof,1);
isBdEdge = true(NEdof,1);
isBdEdge(mapper) = false;
bdidx(isBdEdge) = 1;
Tbd = spdiags(bdidx,0,NEdof,NEdof);
T = spdiags(1-bdidx,0,NEdof,NEdof);
A = T*Atotal*T + Tbd;

ub = tu;
ub(mapper) = 0;
rhsB = Atotal*ub;
f = rhs - rhsB;
f(isBdEdge) = tu(isBdEdge);

option.outsolver = 'cg';
alpha = am*ones(size(gdof,1),1);
% alpha(femI.eLoc>0) = ap;
% alpha(femI.eLoc==0) = (am+ap)/2;
beta = bm*ones(size(gdof,1),1);
% beta(femI.eLoc>0) = bp;
% beta(femI.eLoc==0) = (bm+bp)/2;
option.alpha = alpha;
option.beta = beta;
option.solver = 'amg';
edge = gdof;
option.isBdEdge = isBdEdge;
option.smoother = 'BD';
D = diag(A); D = D(1:BasicEdgeNum);
M = sparse(1:BasicEdgeNum,1:BasicEdgeNum,D,NEdof,NEdof);
M(BasicEdgeNum+1:end,BasicEdgeNum+1:end) = M(BasicEdgeNum+1:end,BasicEdgeNum+1:end) +...
    A(BasicEdgeNum+1:end,BasicEdgeNum+1:end);
option.M = M;
%option.M = A(BasicEdgeNum+1:end,BasicEdgeNum+1:end);
option.blkId = BasicEdgeNum;
[x,info] = amgMaxwell(A,f,node,edge,option);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% generate exact solution

% [gew,gex,gey,gez] = gaussPedge(node(gdof(:,1),:),node(gdof(:,2),:),1);
%
%
% tu1 = sum(feval(pde.exactu1,gex,gey,gez).*gew,2);
% tu2 = sum(feval(pde.exactu2,gex,gey,gez).*gew,2);
% tu3 = sum(feval(pde.exactu3,gex,gey,gez).*gew,2);
%
% tgt = node(gdof(:,2),:) - node(gdof(:,1),:);
% tgt = tgt./sum(tgt.^2,2).^(1/2);
% tu = tu1.*tgt(:,1) + tu2.*tgt(:,2) + tu3.*tgt(:,3);
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %% test
% DC = full(StiffCurl*tu);
% DU = full(SU*tu);

% for i = 1:CNP
%     vv = sum(tu(currentElem(1,:)).*(squeeze(GPA1(1,:,:))'),1);
%     if sum(abs(vv)) > 10^(-8)
%         stp
%     end
% end

repo: ifem


text:
function a = r8col_permute ( m, n, a, p )

%*****************************************************************************80
%
%% R8COL_PERMUTE permutes an R8COL in place.
%
%  Discussion:
%
%    An R8COL is an M by N array of R8's, regarded
%    as an array of N columns of length M.
%
%    This routine permutes an array of real "objects", but the same
%    logic can be used to permute an array of objects of any arithmetic
%    type, or an array of objects of any complexity.  The only temporary
%    storage required is enough to store a single object.  The number
%    of data movements made is N + the number of cycles of order 2 or more,
%    which is never more than N + N/2.
%
%  Example:
%
%    Input:
%
%      M = 2
%      N = 5
%      P = (   2,    4,    5,    1,    3 )
%      A = ( 1.0,  2.0,  3.0,  4.0,  5.0 )
%          (11.0, 22.0, 33.0, 44.0, 55.0 )
%
%    Output:
%
%      A    = (  2.0,  4.0,  5.0,  1.0,  3.0 )
%             ( 22.0, 44.0, 55.0, 11.0, 33.0 ).
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 January 2007
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the size of the objects.
%
%    Input, integer N, the number of objects.
%
%    Input, real A(M,N), the array to be permuted.
%
%    Input, integer P(N), the permutation.  P(I) = J means
%    that the I-th element of the output array should be the J-th
%    element of the input array.  P must be a legal permutation
%    of the integers from 1 to N, otherwise the algorithm will
%    fail catastrophically.
%
%    Output, real A(M,N), the permuted array.
%

%
%  Search for the next element of the permutation that has not been used.
%
  for istart = 1 : n

    if ( p(istart) < 0 )

      continue

    elseif ( p(istart) == istart )

      p(istart) = - p(istart);
      continue

    else

      a_temp(1:m) = a(1:m,istart);
      iget = istart;
%
%  Copy the new value into the vacated entry.
%
      while ( 1 )

        iput = iget;
        iget = p(iget);

        p(iput) = - p(iput);

        if ( iget < 1 || n < iget )
          fprintf ( 1, '\n' );
          fprintf ( 1, 'R8COL_PERMUTE - Fatal error!\n' );
          fprintf ( 1, '  A permutation index is out of range.\n' );
          fprintf ( 1, '  P(%d) = %d\n', iput, iget );
          error ( 'R8COL_PERMUTE - Fatal error!' );
        end

        if ( iget == istart )
          a(1:m,iput) = a_temp(1:m)';
          break
        end

        a(1:m,iput) = a(1:m,iget);

      end

    end

  end
%
%  Restore the signs of the entries.
%
% p(1:n) = -p(1:n);

  return
end

author: johannesgerer

repo: jburkardt-m


text:
classdef TestBrightEdges
    %TestBrightEdges

    methods (Static)
        function test_1
            img = imread(fullfile(mexopencv.root(),'test','balloon.jpg'));
            out = cv.BrightEdges(img, 'Contrast',1);
            validateattributes(out, {'uint8'}, ...
                {'size',[size(img,1) size(img,2)]});
        end
    end

end

author: kyamagu

repo: mexopencv

save_path: github-repos/MATLAB/kyamagu-mexopencv


text:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%	Implemetation of the tracker described in paper
%	"MEEM: Robust Tracking via Multiple Experts using Entropy Minimization",
%   Jianming Zhang, Shugao Ma, Stan Sclaroff, ECCV, 2014
%
%	Copyright (C) 2014 Jianming Zhang
%
%	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/>.
%
%	If you have problems about this software, please contact: jmzhang@bu.edu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function expertsDo(I_vf,lambda,sigma)
global sampler;
global svm_tracker;
global experts;
global config

roi_reg = sampler.roi; roi_reg(3:4) = sampler.roi(3:4)-sampler.roi(1:2);

feature_map = imresize(I_vf,config.ratio,'nearest'); %
ratio_x = size(I_vf,2)/size(feature_map,2);
ratio_y = size(I_vf,1)/size(feature_map,1);
% tmp_mask = zeros(sampler.template_size(1:2));
% tmp_mask(1:2:end,1:2:end) = 1;
% tmp_mask = repmat(tmp_mask,[1,1,size(I_vf,3)]);
patterns = im2colstep(feature_map,[sampler.template_size(1:2), size(I_vf,3)],[1, 1, size(I_vf,3)]);
% patterns = patterns(tmp_mask(:)>0,:); % columnwise

x_sz = size(feature_map,2)-sampler.template_size(2)+1;
y_sz = size(feature_map,1)-sampler.template_size(1)+1;
[X Y] = meshgrid(1:x_sz,1:y_sz);
temp = repmat(svm_tracker.output,[numel(X),1]);
temp(:,1) = (X(:)-1)*ratio_x + sampler.roi(1);
temp(:,2) = (Y(:)-1)*ratio_y + sampler.roi(2);
state = temp;

%% select expert

label_prior = fspecial('gaussian',[y_sz,x_sz],sigma);
label_prior_neg = ones(size(label_prior))/numel(label_prior);

% compute log likelihood and entropy
n = numel(experts);
score_temp = zeros(n,1);
rect_temp = zeros(n,4);

if config.debug
    loglik_vec=[];
    ent_vec=[];
    figure(3)
end

kernel_size = sampler.template_size(1:2)*0.5;%half template size;
rad = 0.5*min(sampler.template_size(1:2));

mask_temp = zeros(y_sz,x_sz);
idx_temp = [];
svm_scores = [];
svm_score = {};
svm_density = {};
peaks_collection = {};
peaks = zeros(n,2);
peaks_pool = [];

[X Y] = meshgrid(1:round(rad):x_sz,1:round(rad):y_sz);

for i = 1:n
    % find the highest peak
    svm_score{i} = -(experts{i}.w*patterns+experts{i}.Bias);
    svm_density{i} = normcdf(svm_score{i},0,1).*label_prior(:)';
    [val idx] = max(svm_density{i});
    best_rect = state(idx,:);
    rect_temp(i,:) = best_rect;
    svm_scores(i) = svm_score{i}(idx);
    idx_temp(i) = idx;
    [r c] = ind2sub(size(mask_temp),idx);
    peaks(i,:) = [r c];

    % find the possible peaks

    density_map = reshape(svm_density{i},y_sz,[]);
    density_map = (density_map - min(density_map(:)))/(max(density_map(:)) - min(density_map(:)));
    mm = (imdilate(density_map,strel('square',round(rad))) == density_map) & density_map > 0.9;
    [rn cn] = ind2sub(size(mask_temp),find(mm));
    peaks_pool = cat(1,peaks_pool,[rn cn]);
    peaks_collection{i} = [rn cn];
%     mask_temp(r,c) = 1;
end
peaks_orig = peaks;
% merg peaks
peaks = mergePeaks(peaks,rad);
peaks_pool = mergePeaks(peaks_pool,rad);
mask_temp(sub2ind(size(mask_temp),round(peaks(:,1)),round(peaks(:,2)))) = 1;

%%
for i = 1:n

    dis = pdist2(peaks_pool,peaks_collection{i});
    [rr cc] = ind2sub([size(peaks_pool,1),size(peaks_collection{i},1)],find(dis < rad));
    [C,ia,ic] = unique(cc);
    peaks_temp = peaks_pool;
    peaks_temp(rr(ia),:) = peaks_collection{i}(cc(ia),:);
    mask = zeros(size(mask_temp));
    mask(sub2ind(size(mask_temp),round(peaks_temp(:,1)),round(peaks_temp(:,2)))) = 1;
    mask = mask>0;

    [loglik ent] = getLogLikelihoodEntropy(svm_score{i}(mask(:)),label_prior(mask(:)),label_prior_neg(mask(:)));
    if config.debug
        loglik_vec(end+1) = loglik;
        ent_vec(end+1) = ent;
        subplot(2,4,i)
        imagesc(reshape(svm_score{i},y_sz,[]));
        colorbar
        subplot(2,4,i+4)
        imagesc(reshape(mask,y_sz,[]))
    end

    experts{i}.score(end+1) =  loglik - lambda*ent;
    score_temp(i) = sum(experts{i}.score(max(end+1-config.entropy_score_winsize,1):end));
end

%%
svm_tracker.best_expert_idx = numel(score_temp);
if numel(score_temp) >= 2 && config.use_experts
    [val idx] = max(score_temp(1:end-1));
    if score_temp(idx) > score_temp(end) && size(peaks,1) > 1%svm_scores(idx) > config.svm_thresh
        % recover previous version
%         output = svm_tracker.output;
%         experts{end}.snapshot = svm_tracker;
        experts{end}.score = experts{idx}.score;
        svm_tracker = experts{idx}.snapshot;
%         svm_tracker.output = rect_temp(idx,:);
        svm_tracker.best_expert_idx = idx;
%         experts([idx end]) = experts([end idx]);
    end
end
svm_tracker.output = rect_temp(svm_tracker.best_expert_idx,:);
svm_tracker.confidence = svm_scores(svm_tracker.best_expert_idx);
svm_tracker.output_exp = rect_temp(end,:);
svm_tracker.confidence_exp = svm_scores(end);

% svm_tracker.w = experts{svm_tracker.best_expert_idx}.w;
% svm_tracker.Bias = experts{svm_tracker.best_expert_idx}.Bias;

if config.debug

    for i = 1:n
        subplot(2,4,i)
        if i == svm_tracker.best_expert_idx
            color = [1 0 0];
        else
            color = [1 1 1];
        end
        text(0,1,num2str(experts{i}.score(end)),'BackgroundColor',color);
        text(15,1,num2str(score_temp(i)),'BackgroundColor',color);
        text(0,3,num2str(loglik_vec(i)),'BackgroundColor',color);
        text(15,3,num2str(ent_vec(i)),'BackgroundColor',color);
    end
    figure(2)
    imagesc(mask_temp)
    figure(1)
end

%% update training sample
% approximately 200 training samples
step = round(sqrt((y_sz*x_sz)/120));
mask_temp = zeros(y_sz,x_sz);
mask_temp(1:step:end,1:step:end) = 1;
mask_temp = mask_temp > 0;
sampler.patterns_dt = patterns(:,mask_temp(:))';
sampler.state_dt = state(mask_temp(:),:);
sampler.costs = 1 - getIOU(sampler.state_dt,svm_tracker.output);
if min(sampler.costs)~=0
    sampler.state_dt = [sampler.state_dt; rect_temp(svm_tracker.best_expert_idx,:)];
    sampler.patterns_dt = [sampler.patterns_dt; patterns(:,idx_temp(svm_tracker.best_expert_idx))'];
    sampler.costs = [sampler.costs;0];
end

% better localization and add the predicted state and pattern
% output = svm_tracker.output;
% output(1:2) = output(1:2) - sampler.roi(1:2) + 1;
% [shift pattern_loc] = localize(I_vf,...
%     -reshape(svm_tracker.w,config.template_sz(1),config.template_sz(2),[]),...
%     output,5);
% svm_tracker.output(1:2) = svm_tracker.output(1:2) + shift;

end

function merged_peaks = mergePeaks(peaks, rad)

dis_mat = pdist2(peaks,peaks) + diag(inf*ones(size(peaks,1),1));
while min(dis_mat(:)) < rad && size(peaks,1) > 1
    [val idx] = min(dis_mat(:));
    [id1 id2] = ind2sub(size(dis_mat),idx);
    merged_peak = 0.5*(peaks(id1,:) + peaks(id2,:));
    peaks([id1 id2],:) = [];
    peaks = [peaks;merged_peak];
    dis_mat = pdist2(peaks,peaks) + diag(inf*ones(size(peaks,1),1));
end

merged_peaks = peaks;

end

author: flyers

repo: drone-tracking


text:
%% Homework 10
% We Wish to solve the following partial differential equation:
%
% $\frac{\partial T}{\partial t}+u \frac{\partial T}{\partial x}=\alpha \frac{\partial^2 T}{\partial x^2}$
%
%% Part 2: Convection and Diffusion (Assume $\alpha=0.001$)
% Now using the same initial boundary conditions in Part 1, Solve the
% convection-diffusion equation. Repeat part(a): i and ii with the
% addition of the second-order central difference for the diffusion term.
%
% Following the condition: $cfl < 1$
%  where $cfl = u \frac{\Delta t}{\Delta x}$
%
% And now define $\alpha*\frac{\Delta t}{(\Delta x)^2}$ in the range [0,1]

clc
clear all
close all

% Initial Parameters
cfl    = 0.7;
u      = 0.08;
dx     = 1/50; % at least 51 grid points
dt     = cfl*dx/u;
alpha  = 0.001;
beta   = alpha*dt/dx^2;
t_end  = 8;

% Grid
x = 0:dx:1;
t = 0:dt:t_end;
X = length(x);
N = length(t);
T = zeros(N,X);

% Initial Condition
T(1,1:12)  = 1.-(10.*x(1:12)-1).^2;
T(1,12:51) = 0;

% Plot options
xmin=1; xmax=51; ymin=-0.2; ymax=1.7;

%% Barckward Euler time advancement
% a) Using Explicit Euler time advancement and second-order central
% difference for the spatial derivate.

% Numerical Scheme:
T(:,1) = 0; T(:,X) = 0;
for n = 1:N-1
    for j = 2:X-1
    T(n+1,j) = T(n,j) - 1/2*cfl*(T(n,j+1) - T(n,j-1))...
        + beta*(T(n,j+1) - 2*T(n,j) + T(n,j-1));
    end
end

% Plot Figures:
figure
Axis([xmin xmax ymin ymax])
Title(['Backward Euler Scheme',', CFL = ',num2str(cfl),', Beta = ',num2str(beta)])
hold on
plot(T(1,:),'ro')
plot(T(floor(N/2),:),'b+')
plot(T(N,:),'gs')
legend('@ t = 0','@ t = 4','@ t = 8')
hold off

%% Leapfrog time advancement
% b) Using Leapfrog time advancement and second-order central
% difference for the spatial derivate and a central difference for the
% diffusion term.
%
% This is a multistep method that is we need to use the first two lines on
% the forward euler to start our new scheme.

T2 = zeros(N,X);
T2(1,:) = T(1,:);
T2(2,:) = T(2,:);
T2(1,1) = 0; T2(1,X) = 0;

% Numerical Scheme:
T2(:,1) = 0; T2(:,X) = 0;
for n = 2:N-1
    for j = 2:X-1
    T2(n+1,j) = T2(n-1,j) - cfl*(T2(n,j+1) - T2(n,j-1))...
        + 2*beta*(T2(n,j+1) - 2*T2(n,j) + T2(n,j-1));
    end
end

% Plot Figures:
figure
%Axis([xmin xmax ymin ymax])
Title(['Leapfrog Scheme',', CFL = ',num2str(cfl),', Beta = ',num2str(beta)])
hold on
plot(T2(1,:),'-.ro')
plot(T2(floor(N/2),:),'-.b+')
plot(T2(N,:),'-.gs')
legend('@ t = 0','@ t = 4','@ t = 8')
hold off

%% Leapfrog time advancement with the diffusion term lagged in time.
% c) Using Leapfrog time advancement and second-order central
% difference for the spatial derivate and a central difference for the
% diffusion term. But this time the diffusion term would be evaluated at
% step "n-1" rather than "n".
%
% This is a multistep method that is we need to use the first two lines on
% the forward euler to start our new scheme.

T3 = zeros(N,X);
T3(1,:) = T(1,:);
T3(2,:) = T(2,:);
T3(1,1) = 0; T3(1,X) = 0;

% Numerical Scheme:
T3(:,1) = 0; T3(:,X) = 0;
for n = 2:N-1
    for j = 2:X-1
    T3(n+1,j) = T3(n-1,j) - cfl*(T3(n,j+1) - T3(n,j-1))...
        + 2*beta*(T3(n-1,j+1) - 2*T3(n-1,j) + T3(n-1,j-1));
    end
end

% Plot Figures:
figure
Axis([xmin xmax ymin ymax])
Title(['Leapfrog with diff-term lagged Scheme',', CFL = ',num2str(cfl),', Beta = ',num2str(beta)])
hold on
plot(T3(1,:),'ro')
plot(T3(floor(N/2),:),'b+')
plot(T3(N,:),'gs')
legend('@ t = 0','@ t = 4','@ t = 8')
hold off

%% Discusion
% How the presence of the diffusion term affected the physical behavior of
% the solution and stability properties of the numerical solution?

repo: Aero-matlab


text:
function [sDataCT] = correctPatient51_HN(pathDATA)
% -------------------------------------------------------------------------
% function [sDataCT] = correctPatient51_HN(pathDATA)
% -------------------------------------------------------------------------
% DESCRIPTION:
% Reason of existence of the function: the RTstruct of the CT scan of
% patient 51 is incomplete (most probably it was previously incorrectly
% processed by MIM). The problem is solved by copying and processing the
% 'contour' structure of the PET scan.
% -------------------------------------------------------------------------
% INPUTS:
% - pathDATA: Full path to the 'DATA' folder of the HN workspace.
% -------------------------------------------------------------------------
% OUTPUTS:
% - sData: Corrected sData.
% -------------------------------------------------------------------------
% 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(pathDATA)
load('Patient51_CT'), sDataCT = sData;
load('Patient51_PET'), sDataPET = sData;

% ADJUSTING CT CONTOURS ACCORDING TO PET CONTOURS.
sDataCT{2}.scan.contour = sDataPET{2}.scan.contour;
nContour = length(sDataPET{2}.scan.contour);
downF = sDataCT{2}.scan.pixelW/sDataPET{2}.scan.pixelW;
offsetStart = round(abs(sDataPET{3}(1).ImagePositionPatient(2)-sDataCT{3}(1).ImagePositionPatient(2))*sDataCT{2}.scan.pixelW);
offsetEnd = round((size(sDataCT{2}.scan.volume,1)*sDataCT{2}.scan.pixelW - size(sDataPET{2}.scan.volume,1)*sDataPET{2}.scan.pixelW)/sDataCT{2}.scan.pixelW)-offsetStart;
for i = 1:nContour
    sDataCT{2}.scan.contour(i).boxBound(1:2,1) = round((sDataPET{2}.scan.contour(i).boxBound(1:2,1)/downF + offsetStart)); sDataCT{2}.scan.contour(i).boxBound(1:2,2) = round((sDataPET{2}.scan.contour(i).boxBound(1:2,2)/downF + offsetEnd));
    boxBound = sDataCT{2}.scan.contour(i).boxBound;
    nSlices = size(sDataCT{2}.scan.contour(i).boxMask,3);
    sDataCT{2}.scan.contour(i).boxMask = zeros(boxBound(1,2)-boxBound(1,1)+1,boxBound(2,2)-boxBound(2,1)+1,nSlices);
    for j = 1:nSlices
        sDataCT{2}.scan.contour(i).boxMask(:,:,j) = imresize(sDataPET{2}.scan.contour(i).boxMask(:,:,j),[boxBound(1,2)-boxBound(1,1)+1,boxBound(2,2)-boxBound(2,1)+1],'nearest');
    end
end

cd(startpath)
end

author: mvallieres

repo: radiomics


text:
function varargout = ghostscript(cmd)
%GHOSTSCRIPT  Calls a local GhostScript executable with the input command
%
% Example:
%   [status result] = ghostscript(cmd)
%
% Attempts to locate a ghostscript executable, finally asking the user to
% specify the directory ghostcript was installed into. The resulting path
% is stored for future reference.
%
% Once found, the executable is called with the input command string.
%
% This function requires a Ghostscript installation on your system.
% You can download Ghostscript from http://ghostscript.com (Windows/Linux)
% or http://pages.uoregon.edu/koch (MacOS).
%
% IN:
%   cmd - Command string to be passed into ghostscript.
%
% OUT:
%   status - 0 iff command ran without problem.
%   result - Output from ghostscript.

% Copyright: Oliver Woodford, 2009-2015, Yair Altman 2015-
%{
% Thanks to Jonas Dorn for the fix for the title of the uigetdir window on Mac OS.
% Thanks to Nathan Childress for the fix to default location on 64-bit Windows systems.
% 27/04/11 - Find 64-bit Ghostscript on Windows. Thanks to Paul Durack and
%            Shaun Kline for pointing out the issue
% 04/05/11 - Thanks to David Chorlian for pointing out an alternative
%            location for gs on linux.
% 12/12/12 - Add extra executable name on Windows. Thanks to Ratish
%            Punnoose for highlighting the issue.
% 28/06/13 - Fix error using GS 9.07 in Linux. Many thanks to Jannick
%            Steinbring for proposing the fix.
% 24/10/13 - Fix error using GS 9.07 in Linux. Many thanks to Johannes
%            for the fix.
% 23/01/14 - Add full path to ghostscript.txt in warning. Thanks to Koen
%            Vermeer for raising the issue.
% 27/02/15 - If Ghostscript croaks, display suggested workarounds
% 30/03/15 - Improved performance by caching status of GS path check, if ok
% 14/05/15 - Clarified warning message in case GS path could not be saved
% 29/05/15 - Avoid cryptic error in case the ghostscipt path cannot be saved (issue #74)
% 10/11/15 - Custom GS installation webpage for MacOS. Thanks to Andy Hueni via FEX
% 15/01/20 - Various message cleanups/fixes in case of errors
%}

    try
        % Call ghostscript
        [varargout{1:nargout}] = system([gs_command(gs_path()) cmd]);
    catch err
        % Display possible workarounds for Ghostscript croaks
        url1 = 'https://github.com/altmany/export_fig/issues/12#issuecomment-61467998';  % issue #12
        url2 = 'https://github.com/altmany/export_fig/issues/20#issuecomment-63826270';  % issue #20
        hg2_str = ''; if using_hg2, hg2_str = ' or Matlab R2014a'; end
        fprintf(2, 'Ghostscript error. Rolling back to GS 9.10%s may possibly solve this:\n * %s ', hg2_str, hyperlink(url1));
        if using_hg2
            fprintf(2, '(GS 9.10)\n * %s (R2014a)', hyperlink(url2));
        end
        fprintf('\n\n');
        if ismac || isunix
            url3 = 'https://github.com/altmany/export_fig/issues/27';  % issue #27
            fprintf(2, 'Alternatively, this may possibly be due to a font path issue:\n * %s\n\n', hyperlink(url3));
            % issue #20
            % TODO: in Unix/Mac, find a way to automatically determine whether to use "export" (bash) or "setenv" (csh/tcsh)
            if isdeployed
                url = [mfilename '.m'];
            else
                fpath = which(mfilename);
                if isempty(fpath), fpath = [mfilename('fullpath') '.m']; end
                url = ['' fpath ''];
            end
            fprintf(2, 'Alternatively, if you are using csh, modify shell_cmd from "export ..." to "setenv ..."\nat the bottom of %s\n\n', url);
        end
        rethrow(err);
    end
end

function path_ = gs_path
    % Return a valid path
    % Start with the currently set path
    path_ = user_string('ghostscript');
    % Check the path works
    if check_gs_path(path_)
        return
    end
    % Check whether the binary is on the path
    if ispc
        bin = {'gswin32c.exe', 'gswin64c.exe', 'gs'};
    else
        bin = {'gs'};
    end
    for a = 1:numel(bin)
        path_ = bin{a};
        if check_store_gs_path(path_)
            return
        end
    end
    % Search the obvious places
    if ispc
        default_location = 'C:\Program Files\gs\';
        dir_list = dir(default_location);
        if isempty(dir_list)
            default_location = 'C:\Program Files (x86)\gs\'; % Possible location on 64-bit systems
            dir_list = dir(default_location);
        end
        executable = {'\bin\gswin32c.exe', '\bin\gswin64c.exe'};
        ver_num = 0;
        % If there are multiple versions, use the newest
        for a = 1:numel(dir_list)
            ver_num2 = sscanf(dir_list(a).name, 'gs%g');
            if ~isempty(ver_num2) && ver_num2 > ver_num
                for b = 1:numel(executable)
                    path2 = [default_location dir_list(a).name executable{b}];
                    if exist(path2, 'file') == 2
                        path_ = path2;
                        ver_num = ver_num2;
                    end
                end
            end
        end
        if check_store_gs_path(path_)
            return
        end
    else
        executable = {'/usr/bin/gs', '/usr/local/bin/gs'};
        for a = 1:numel(executable)
            path_ = executable{a};
            if check_store_gs_path(path_)
                return
            end
        end
    end
    % Ask the user to enter the path
    while true
        if strncmp(computer, 'MAC', 3) % Is a Mac
            % Give separate warning as the uigetdir dialogue box doesn't have a
            % title on MacOS
            uiwait(warndlg('Ghostscript installation not found - please locate the program.', 'Ghostscript'))
            base = uigetdir('/', 'Ghostcript program location');
        else
            base = uigetdir('/', 'Ghostcript program not found - please locate it');
        end
        if isequal(base, 0)
            % User hit cancel or closed window
            break;
        end
        base = [base filesep]; %#ok<AGROW>
        bin_dir = {'', ['bin' filesep], ['lib' filesep]};
        for a = 1:numel(bin_dir)
            for b = 1:numel(bin)
                path_ = [base bin_dir{a} bin{b}];
                if exist(path_, 'file') == 2
                    if check_store_gs_path(path_)
                        return
                    end
                end
            end
        end
    end
    if ismac
        url = 'http://pages.uoregon.edu/koch';
    else
        url = 'http://ghostscript.com';
    end
    error('Ghostscript:NotFound', 'Ghostscript not found. Have you installed it from %s ?', hyperlink(url));
end

function good = check_store_gs_path(path_)
    % Check the path is valid
    good = check_gs_path(path_);
    if ~good
        return
    end
    % Update the current default path to the path found
    if ~user_string('ghostscript', path_)
        %filename = fullfile(fileparts(which('user_string.m')), '.ignore', 'ghostscript.txt');
        [unused, filename] = user_string('ghostscript'); %#ok<ASGLU>
        warning('Ghostscript:path', 'Path to ghostscript installation could not be saved in %s (perhaps a permissions issue). You can manually create this file and set its contents to %s, to improve performance in future invocations (this warning is safe to ignore).', filename, path_);
        return
    end
end

function good = check_gs_path(path_)
    persistent isOk
    if isempty(path_)
        isOk = false;
    elseif ~isequal(isOk,true)
        % Check whether the path is valid
        [status, message] = system([gs_command(path_) '-h']); %#ok<ASGLU>
        isOk = status == 0;
    end
    good = isOk;
end

function cmd = gs_command(path_)
    % Initialize any required system calls before calling ghostscript
    % TODO: in Unix/Mac, find a way to automatically determine whether to use "export" (bash) or "setenv" (csh/tcsh)
    shell_cmd = '';
    if isunix
        shell_cmd = 'export LD_LIBRARY_PATH=""; '; % Avoids an error on Linux with GS 9.07
    end
    if ismac
        shell_cmd = 'export DYLD_LIBRARY_PATH=""; ';  % Avoids an error on Mac with GS 9.07
    end
    % Construct the command string
    cmd = sprintf('%s"%s" ', shell_cmd, path_);
end

author: altmany


text:
function x = vec(X)
% x = vec(X)
%
% Y = VEC(x)  Given an m x n matrix x, this produces the vector Y of length
%   m*n that contains the columns of the matrix x, stacked below each other.
%
% See also mat.

% This file is part of SeDuMi 1.1 by Imre Polik and Oleksandr Romanko
% Copyright (C) 2005 McMaster University, Hamilton, CANADA  (since 1.1)
%
% Copyright (C) 2001 Jos F. Sturm (up to 1.05R5)
%   Dept. Econometrics & O.R., Tilburg University, the Netherlands.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% Affiliation SeDuMi 1.03 and 1.04Beta (2000):
%   Dept. Quantitative Economics, Maastricht University, the Netherlands.
%
% Affiliations up to SeDuMi 1.02 (AUG1998):
%   CRL, McMaster University, Canada.
%   Supported by the Netherlands Organization for Scientific Research (NWO).
%
% 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 Street, Fifth Floor, Boston, MA
% 02110-1301, USA

x = reshape(X,numel(X),1);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function X = cs_qleft (V, Beta, p, Y)
%CS_QLEFT apply Householder vectors on the left.
%   X = cs_qleft(V,Beta,p,Y) computes X = Hn*...*H2*H1*P*Y = Q'*Y where Q is
%   represented by the Householder vectors V, coefficients Beta, and
%   permutation p.  p can be [], which denotes the identity permutation.
%
%   Example:
%       Prob = UFget ('HB/well1033') ; A = Prob.A ; [m n] = size (A) ;
%       b = rand (m,1) ;
%       [V,beta,p,R] = cs_qr (A) ; % QR factorization of A(p,:)
%       b1 = cs_qleft (V, beta, p, b) ;
%       x1 = R (1:n,1:n) \ b1 (1:n) ;
%       x2 = A\b ;
%       norm (x1-x2)
%
%   See also CS_QR, CS_QRIGHT.

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

[m2 n] = size (V) ;
[m ny] = size (Y) ;
X = Y ;
if (m2 > m)
    if (issparse (Y))
        X = [X ; sparse(m2-m,ny)] ;
    else
        X = [X ; zeros(m2-m,ny)] ;
    end
end
if (~isempty (p))
    X = X (p,:) ;
end
for k = 1:n
    X = X - V (:,k) * (Beta (k) * (V (:,k)' * X)) ;
end

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
% Anchored Neighborhood Regression for Fast Example-Based Super-Resolution
% Example code
%
% March 22, 2013. Radu Timofte, VISICS @ KU Leuven
%
% Revised version: (includes all [1] methods)
% October 3, 2013. Radu Timofte, CVL @ ETH Zurich
%
% Updated version: (adds A+ methods [2])
% September 5, 2014. Radu Timofte, CVL @ ETH Zurich
% %
% Please reference to both:
% [1] Radu Timofte, Vincent De Smet, Luc Van Gool.
% Anchored Neighborhood Regression for Fast Example-Based Super-Resolution.
% International Conference on Computer Vision (ICCV), 2013.
%
% [2] Radu Timofte, Vincent De Smet, Luc Van Gool.
% A+: Adjusted Anchored Neighborhood Regression for Fast Super-Resolution.
% Asian Conference on Computer Vision (ACCV), 2014.
%
% For any questions, email me by timofter@vision.ee.ethz.ch
%

clear;

p = pwd;
addpath(fullfile(p, '/methods'));  % the upscaling methods

addpath(fullfile(p, '/ksvdbox')) % K-SVD dictionary training algorithm

addpath(fullfile(p, '/ompbox')) % Orthogonal Matching Pursuit algorithm

imgscale = 1; % the scale reference we work with
flag = 1;       % flag = 0 - only GR, ANR, A+, and bicubic methods, the other get the bicubic result by default
                % flag = 1 - all the methods are applied

upscaling = 4; % the magnification factor x2, x3, x4...

input_dir = 'Set5'; % Directory with input images from Set5 image dataset
%input_dir = 'Set14'; % Directory with input images from Set14 image dataset

pattern = '*.bmp'; % Pattern to process

dict_sizes = [2 4 8 16 32 64 128 256 512 1024 2048 4096 8192 16384 32768 65536];
neighbors = [1:1:12, 16:4:32, 40:8:64, 80:16:128, 256, 512, 1024];
%d = 7
%for nn=1:28
%nn= 28

clusterszA = 2048; % neighborhood size for A+

disp('The experiment corresponds to the results from Table 2 in the referenced [1] and [2] papers.');

disp(['The experiment uses ' input_dir ' dataset and aims at a magnification of factor x' num2str(upscaling) '.']);
if flag==1
    disp('All methods are employed : Bicubic, Yang et al., Zeyde et al., GR, ANR, NE+LS, NE+NNLS, NE+LLE, A+ (0.5 mil), A+, A+ (16 atoms).');
else
    disp('We run only for Bicubic, GR, ANR and A+ methods, the other get the Bicubic result by default.');
end

fprintf('\n\n');

for d=10    %1024
    %d = 9; % 512
    %d = 8; %256
    %d = 7; %128
    %d = 6; % 64
    %d = 5; % 32
    %d=4;  %16
    %d=3;  %8
    %d=2; %4
    %d=1; %2

    tag = [input_dir '_x' num2str(upscaling) '_' num2str(dict_sizes(d)) 'atoms'];

    disp(['Upscaling x' num2str(upscaling) ' ' input_dir ' with Zeyde dictionary of size = ' num2str(dict_sizes(d))]);

    mat_file = ['conf_Zeyde_' num2str(dict_sizes(d)) '_finalx' num2str(upscaling)];

    if exist([mat_file '.mat'],'file')
        disp(['Load trained dictionary...' mat_file]);
        load(mat_file, 'conf');
    else
        disp(['Training dictionary of size ' num2str(dict_sizes(d)) ' using Zeyde approach...']);
        % Simulation settings
        conf.scale = upscaling; % scale-up factor
        conf.level = 1; % # of scale-ups to perform
        conf.window = [3 3]; % low-res. window size
        conf.border = [1 1]; % border of the image (to ignore)

        % High-pass filters for feature extraction (defined for upsampled low-res.)
        conf.upsample_factor = upscaling; % upsample low-res. into mid-res.
        O = zeros(1, conf.upsample_factor-1);
        G = [1 O -1]; % Gradient
        L = [1 O -2 O 1]/2; % Laplacian
        conf.filters = {G, G.', L, L.'}; % 2D versions
        conf.interpolate_kernel = 'bicubic';

        conf.overlap = [1 1]; % partial overlap (for faster training)
        if upscaling <= 2
            conf.overlap = [1 1]; % partial overlap (for faster training)
        end

        startt = tic;
        conf = learn_dict(conf, load_images(...
            glob('CVPR08-SR/Data/Training', '*.bmp') ...
            ), dict_sizes(d));
        conf.overlap = conf.window - [1 1]; % full overlap scheme (for better reconstruction)
        conf.trainingtime = toc(startt);
        toc(startt)

        save(mat_file, 'conf');

        % train call
    end

    if dict_sizes(d) < 1024
        lambda = 0.01;
    elseif dict_sizes(d) < 2048
        lambda = 0.1;
    elseif dict_sizes(d) < 8192
        lambda = 1;
    else
        lambda = 5;
    end

    %% GR
    if dict_sizes(d) < 10000
        conf.ProjM = inv(conf.dict_lores'*conf.dict_lores+lambda*eye(size(conf.dict_lores,2)))*conf.dict_lores';
        conf.PP = (1+lambda)*conf.dict_hires*conf.ProjM;
    else
        % here should be an approximation
        conf.PP = zeros(size(conf.dict_hires,1), size(conf.V_pca,2));
        conf.ProjM = [];
    end

    conf.filenames = glob(input_dir, pattern); % Cell array

    conf.desc = {'Original', 'Bicubic', 'Yang et al.', ...
        'Zeyde et al.', 'Our GR', 'Our ANR', ...
        'NE+LS','NE+NNLS','NE+LLE','Our A+ (0.5mil)','Our A+', 'Our A+ (16atoms)'};
    conf.results = {};

    %conf.points = [1:10:size(conf.dict_lores,2)];
    conf.points = [1:1:size(conf.dict_lores,2)];

    conf.pointslo = conf.dict_lores(:,conf.points);
    conf.pointsloPCA = conf.pointslo'*conf.V_pca';

    % precompute for ANR the anchored neighborhoods and the projection matrices for
    % the dictionary

    conf.PPs = [];
    if  size(conf.dict_lores,2) < 40
        clustersz = size(conf.dict_lores,2);
    else
        clustersz = 40;
    end
    D = abs(conf.pointslo'*conf.dict_lores);

    for i = 1:length(conf.points)
        [vals idx] = sort(D(i,:), 'descend');
        if (clustersz >= size(conf.dict_lores,2)/2)
            conf.PPs{i} = conf.PP;
        else
            Lo = conf.dict_lores(:, idx(1:clustersz));
            conf.PPs{i} = 1.01*conf.dict_hires(:,idx(1:clustersz))*inv(Lo'*Lo+0.01*eye(size(Lo,2)))*Lo';
        end
    end

    ANR_PPs = conf.PPs; % store the ANR regressors

    save([tag '_' mat_file '_ANR_projections_imgscale_' num2str(imgscale)],'conf');

    %% A+ computing the regressors
    Aplus_PPs = [];

    fname = ['Aplus_x' num2str(upscaling) '_' num2str(dict_sizes(d)) 'atoms' num2str(clusterszA) 'nn_5mil.mat'];

    if exist(fname,'file')
       load(fname);
    else
        %%
       disp('Compute A+ regressors');
       ttime = tic;
       tic
       [plores phires] = collectSamplesScales(conf, load_images(...
        glob('CVPR08-SR/Data/Training', '*.bmp')), 12, 0.98);

        if size(plores,2) > 5000000
            plores = plores(:,1:5000000);
            phires = phires(:,1:5000000);
        end
        number_samples = size(plores,2);

        % l2 normalize LR patches, and scale the corresponding HR patches
        l2 = sum(plores.^2).^0.5+eps;
        l2n = repmat(l2,size(plores,1),1);
        l2(l2<0.1) = 1;
        plores = plores./l2n;
        phires = phires./repmat(l2,size(phires,1),1);
        clear l2
        clear l2n

        llambda = 0.1;

        for i = 1:size(conf.dict_lores,2)
            D = pdist2(single(plores'),single(conf.dict_lores(:,i)'));
            [~, idx] = sort(D);
            Lo = plores(:, idx(1:clusterszA));
            Hi = phires(:, idx(1:clusterszA));
            Aplus_PPs{i} = Hi*inv(Lo'*Lo+llambda*eye(size(Lo,2)))*Lo';
%Aplus_PPs{i} = Hi*(inv(Lo*Lo'+llambda*eye(size(Lo,1)))*Lo)';
        end
        clear plores
        clear phires

        ttime = toc(ttime);
        save(fname,'Aplus_PPs','ttime', 'number_samples');
        toc
    end

    %% A+ (0.5mil) computing the regressors with 0.5 milion training samples
    Aplus05_PPs = [];

    fname = ['Aplus_x' num2str(upscaling) '_' num2str(dict_sizes(d)) 'atoms' num2str(clusterszA) 'nn_05mil.mat'];

    if exist(fname,'file')
       load(fname);
    else
        %%
       disp('Compute A+ (0.5 mil) regressors');
       ttime = tic;
       tic
       [plores phires] = collectSamplesScales(conf, load_images(...
        glob('CVPR08-SR/Data/Training', '*.bmp')), 1,1);

        if size(plores,2) > 500000
            plores = plores(:,1:500000);
            phires = phires(:,1:500000);
        end
        number_samples = size(plores,2);

        % l2 normalize LR patches, and scale the corresponding HR patches
        l2 = sum(plores.^2).^0.5+eps;
        l2n = repmat(l2,size(plores,1),1);
        l2(l2<0.1) = 1;
        plores = plores./l2n;
        phires = phires./repmat(l2,size(phires,1),1);
        clear l2
        clear l2n

        llambda = 0.1;

        for i = 1:size(conf.dict_lores,2)
            D = pdist2(single(plores'),single(conf.dict_lores(:,i)'));
            [~, idx] = sort(D);
            Lo = plores(:, idx(1:clusterszA));
            Hi = phires(:, idx(1:clusterszA));
            Aplus05_PPs{i} = Hi*inv(Lo'*Lo+llambda*eye(size(Lo,2)))*Lo';
        end
        clear plores
        clear phires

        ttime = toc(ttime);
        save(fname,'Aplus05_PPs','ttime', 'number_samples');
        toc
    end

    %% load the A+ (16 atoms) for comparison results
    conf16 = [];
    fname = ['Aplus_x' num2str(upscaling) '_16atoms' num2str(clusterszA) 'nn_05mil.mat'];
    fnamec = ['Set14_x' num2str(upscaling) '_16atoms_conf_Zeyde_16_finalx' num2str(upscaling) '_ANR_projections_imgscale_' num2str(imgscale) '.mat'];
    if exist(fname,'file') && exist(fnamec,'file')
       kk = load(fnamec);
       conf16 = kk.conf;
       kk = load(fname);
       conf16.PPs = kk.Aplus05_PPs;
       clear kk
    end
    %%
    conf.result_dirImages = qmkdir([input_dir '/results_' tag]);
    conf.result_dirImagesRGB = qmkdir([input_dir '/results_' tag 'RGB']);
    conf.result_dir = qmkdir(['Results-' datestr(now, 'YYYY-mm-dd_HH-MM-SS')]);
    conf.result_dirRGB = qmkdir(['ResultsRGB-' datestr(now, 'YYYY-mm-dd_HH-MM-SS')]);

    %%
    t = cputime;

    conf.countedtime = zeros(numel(conf.desc),numel(conf.filenames));

    res =[];
    for i = 1:numel(conf.filenames)
        f = conf.filenames{i};
        [p, n, x] = fileparts(f);
        [img, imgCB, imgCR] = load_images({f});
        if imgscale<1
            img = resize(img, imgscale, conf.interpolate_kernel);
            imgCB = resize(imgCB, imgscale, conf.interpolate_kernel);
            imgCR = resize(imgCR, imgscale, conf.interpolate_kernel);
        end
        sz = size(img{1});

        fprintf('%d/%d\t"%s" [%d x %d]\n', i, numel(conf.filenames), f, sz(1), sz(2));

        img = modcrop(img, conf.scale^conf.level);
        imgCB = modcrop(imgCB, conf.scale^conf.level);
        imgCR = modcrop(imgCR, conf.scale^conf.level);

            low = resize(img, 1/conf.scale^conf.level, conf.interpolate_kernel);
            if ~isempty(imgCB{1})
                lowCB = resize(imgCB, 1/conf.scale^conf.level, conf.interpolate_kernel);
                lowCR = resize(imgCR, 1/conf.scale^conf.level, conf.interpolate_kernel);
            end

        interpolated = resize(low, conf.scale^conf.level, conf.interpolate_kernel);
        if ~isempty(imgCB{1})
            interpolatedCB = resize(lowCB, conf.scale^conf.level, conf.interpolate_kernel);
            interpolatedCR = resize(lowCR, conf.scale^conf.level, conf.interpolate_kernel);
        end

        res{1} = interpolated;

        if (flag == 1) && (dict_sizes(d) == 1024) && (upscaling==3)
            startt = tic;
            res{2} = {yima(low{1}, upscaling)};
            toc(startt)
            conf.countedtime(2,i) = toc(startt);
        else
            res{2} = interpolated;
        end

        if (flag == 1)
            startt = tic;
            res{3} = scaleup_Zeyde(conf, low);
            toc(startt)
            conf.countedtime(3,i) = toc(startt);
        else
            res{3} = interpolated;
        end

        %if flag == 1
            startt = tic;
            res{4} = scaleup_GR(conf, low);
            toc(startt)
            conf.countedtime(4,i) = toc(startt);
        %else
            %res{4} = interpolated;
        %end

        startt = tic;
        conf.PPs = ANR_PPs;
        res{5} = scaleup_ANR(conf, low);
        toc(startt)
        conf.countedtime(5,i) = toc(startt);

        if flag == 1
            startt = tic;
            if 12 < dict_sizes(d)
                res{6} = scaleup_NE_LS(conf, low, 12);
            else
                res{6} = scaleup_NE_LS(conf, low, dict_sizes(d));
            end
            toc(startt)
            conf.countedtime(6,i) = toc(startt);
        else
            res{6} = interpolated;
        end

        if flag == 1
            startt = tic;
            if 24 < dict_sizes(d)
                res{7} = scaleup_NE_NNLS(conf, low, 24);
            else
                res{7} = scaleup_NE_NNLS(conf, low, dict_sizes(d));
            end
            toc(startt)
            conf.countedtime(7,i) = toc(startt);
        else
            res{7} = interpolated;
        end

        if flag == 1
            startt = tic;
            if 24 < dict_sizes(d)
                res{8} = scaleup_NE_LLE(conf, low, 24);
            else
                res{8} = scaleup_NE_LLE(conf, low, dict_sizes(d));
            end
            toc(startt)
            conf.countedtime(8,i) = toc(startt);
        else
            res{8} = interpolated;
        end

        % A+ (0.5 mil)
        if flag == 1 && ~isempty(Aplus05_PPs)
            fprintf('A+ (0.5mil)\n');
            conf.PPs = Aplus05_PPs;
            startt = tic;
            res{9} = scaleup_ANR(conf, low);
            toc(startt)
            conf.countedtime(9,i) = toc(startt);
        else
            res{9} = interpolated;
        end

        % A+
        if ~isempty(Aplus_PPs)
            fprintf('A+\n');
            conf.PPs = Aplus_PPs;
            startt = tic;
            res{10} = scaleup_ANR(conf, low);
            toc(startt)
            conf.countedtime(10,i) = toc(startt);
        else
            res{10} = interpolated;
        end
        % A+ 16atoms
        if flag == 1 && ~isempty(conf16)
            fprintf('A+ 16atoms\n');
            startt = tic;
            res{11} = scaleup_ANR(conf16, low);
            toc(startt)
            conf.countedtime(11,i) = toc(startt);
        else
            res{11} = interpolated;
        end

        result = cat(3, img{1}, interpolated{1}, res{2}{1}, res{3}{1}, ...
            res{4}{1}, res{5}{1}, res{6}{1}, res{7}{1}, res{8}{1}, ...
            res{9}{1}, res{10}{1}, res{11}{1});

        result = shave(uint8(result * 255), conf.border * conf.scale);

        if ~isempty(imgCB{1})
            resultCB = interpolatedCB{1};
            resultCR = interpolatedCR{1};
            resultCB = shave(uint8(resultCB * 255), conf.border * conf.scale);
            resultCR = shave(uint8(resultCR * 255), conf.border * conf.scale);
        end

        conf.results{i} = {};
        for j = 1:numel(conf.desc)
            conf.results{i}{j} = fullfile(conf.result_dirImages, [n sprintf('[%d-%s]', j, conf.desc{j}) x]);
            imwrite(result(:, :, j), conf.results{i}{j});

            conf.resultsRGB{i}{j} = fullfile(conf.result_dirImagesRGB, [n sprintf('[%d-%s]', j, conf.desc{j}) x]);
            if ~isempty(imgCB{1})
                rgbImg = cat(3,result(:,:,j),resultCB,resultCR);
                rgbImg = ycbcr2rgb(rgbImg);
            else
                rgbImg = cat(3,result(:,:,j),result(:,:,j),result(:,:,j));
            end

            imwrite(rgbImg, conf.resultsRGB{i}{j});
        end
        conf.filenames{i} = f;
    end
    conf.duration = cputime - t;

    % Test performance
    scores = run_comparison(conf);
    process_scores_Tex(conf, scores,length(conf.filenames));

    run_comparisonRGB(conf); % provides color images and HTML summary
    %%
    save([tag '_' mat_file '_results_imgscale_' num2str(imgscale)],'conf','scores');
end
%

author: thomas-koehler

repo: SupER


text:
function c = linkage(c,varargin);

% c = LINKAGE(c) Creates a hierarchical cluster tree for all waveforms in
% the correlation object. It reads the CORR field (must be filled) and
% fills in the LINK field. This iuse is identical to c =
% LINKAGE(c,'average').
%
% c = LINKAGE(c,...) is just a wrapper program which calls the linkage
% function in the matlab stats toolbox.  Usage is the same as the stats
% toolbox version except that the first argument and the returned value are
% correlation objects. In most cases these will be the same object. See
% HELP LINKAGE for details and usage. All options and user controls are the
% same.
%
% ** IMPORTANT NOTE **: There is one significant difference from the
% built-in linkage routine. The native linkage command operates on
% *dissimilarity* information. That is, the opposite of correlation values.
% When linkage is called on a correlation object, it expects max
% correlation values as input. Temporary dissimiliarity information is
% created internally and passed to the built in linkage routine.
%
% Common uses include:
%   c = LINKAGE(c)           % returns clusters of similar waveforms
% 	c = LINKAGE(c,'average') % same as first use
% 	c = LINKAGE(c,'single')  % useful for evolving waveforms

% Author: Michael West, Geophysical Institute, Univ. of Alaska Fairbanks
% $Date$
% $Revision$

if nargin <= 0
    error('Not enough inputs');
end

if ~isa(c,'correlation')
    error('First input must be a correlation object');
end;

if get(c,'TRACES')<2
    error('correlationLinkageTooFewTraces','Correlation object must contain at least two traces to use the LINKAGE function');
end;

% if isempty(get(c,'LAG'))
%     error('LAG field must be filled in input object');
%     error('See correlation/linkage function');
% end;

K = 1.001 - c.C;			% create dissimilarity matrix
K = K - diag(diag(K));			% remove diagonal (required format)
Y = squareform(K);          % transform to "pdist" vector format

if nargin == 1
    c.link = linkage(Y,'average');
else
    c.link = linkage(Y,varargin{:});
end

author: geoscience-community-codes

repo: GISMO