text:
%% DEMO_aorta_build_passive_01
% Below is a demonstration for:
%
% * Building geometry for a subject-specific aorta
% * Creating the FEA mesh
% * Coding the abaqus structure

%% Keywords
%
% * ABAQUS
% * aorta
% * hexahedral elements, hex8

%%

clear; close all; clc;

%% Directory

saveOn=0;

%Define working directory
defaultFolder = fileparts(fileparts(mfilename('fullpath')));
savePath=fullfile(defaultFolder,'data','MAT');
fileName='Concannon_aorta_segmentation.mat';

%Define save names
abaqusInpFileNamePart='tempModel';
abaqusInpFileName=fullfile(savePath,[abaqusInpFileNamePart,'.inp']); %INP file name
matfileSaveName=fullfile(savePath,[abaqusInpFileNamePart,'.mat']); %INP file name

%Load data as structure
dataStruct=load(fullfile(savePath,fileName));

%% Control parameters
%Define smoothing Parameters
pointSpacing=1.7; %
numThickenSteps=2;
smoothFactorCentreLine=0.01; %Cubic smooth spline parameter [0-1] use empty to turn off
smoothFactorSegments=0.01; %Cubic smooth spline parameter [0-1], 0=straight line, 1=cubic
trunkSegmentReduceFactor=1; %

numSmoothTrunk=100; %
numSmoothPants_LAP=100; %Number of Laplacian smoothing iterations for Iliacs
numSmoothPants_HC=100; %Number of HC smoothing iterations for Iliacs
numSmoothBranchAttachments=100; %Number of smoothing iterations for branch origins (50)
numSmoothBranches=100; %Number of smoothing iterations for rest of branches
numOffsetSteps=1; %Number of steps taken to offset mesh by wall thickness
numSmoothOffset=2; %
methodSmoothOffset='LAP'; %Smoothing method Laplacian=agressive
qualityMetricSmoothOffset=120; %
distSmoothGrowth=10; %
numHexSplit=1; %Split elements by factor
t=1; %Time value corresponding to phase in cardiac cycle
%Compressibility factor (Nolan, McGarry)
kfactor=8.56;
%Define thickness information
dataStruct.WallThickness=[1.425	0.9	1	1.025	0.833333333	0.891666667	0.95	0.975	0.9	0.825];
thicknessIndexTransitionRegion=10; %Index of thickness at Iliac bifurcation
thicknessIndexBifurcation=10; %Index of thickness distal to Iliac bifurcation
thicknessIndicesBranches=[8 8 7 7 7 7 7]; %R_Renal_Ori,L_Renal_Ori,SMA_Ori,COE_Ori,LSA_Ori,LCCA_Ori,V_BCA_Ori
wallThickness=dataStruct.WallThickness; %Raw data for wall thickness as a function of location
numThicknessLevels=numel(wallThickness);
%Define material parameter information (from segment fitting [1-10])
numMaterials=100;
materialIndicesSelect=linspace(1,10,10); %Plane indices to interpolate between
dataStruct.mu=0.02; %Constant
dataStruct.k1=0.0001; %Constant
dataStruct.k2=1; %Constant
dataStruct.kappa=0; %Constant
dataStruct.theta1=90; %Circumferential
dataStruct.theta2=-90; %Circumferential
dataStruct.sigact=0; %Constant
dataStruct.ke=0.1; %Constant
dataStruct.Ee=[0.82 0.74 0.63 0.5 0.45 0.4 0.35 0.28 0.27 0.27];
dataStruct.thetaE1=90; %Circumferential
dataStruct.thetaE2=-90; %Circumferential
dataStruct.iswitch=1; %[1=local; 2=global] %Constant
dataStruct.xeps=[0.34 0.31 0.28 0.28 0.25 0.20 0.18 0.20 0.20 0.20];
dataStruct.Vcol=[0.20 0.31 0.28 0.30 0.44 0.44 0.44 0.44 0.45 0.52];
dataStruct.Vela=[0.30 0.28 0.25 0.22 0.22 0.22 0.24	0.20 0.18 0.16];
dataStruct.Vsmc=[0.20 0.20 0.20	0.23 0.24 0.25 0.26	0.26 0.27 0.27];
%Acces data
mu_data=dataStruct.mu;
k1_data=dataStruct.k1;
k2_data=dataStruct.k2;
kappa_data=dataStruct.kappa;
theta1_data=dataStruct.theta1;
theta2_data=dataStruct.theta2;
sigact_data=dataStruct.sigact;
ke_data=dataStruct.ke;
Ee_data=dataStruct.Ee;
thetaE1_data=dataStruct.thetaE1;
thetaE2_data=dataStruct.thetaE2;
xeps_data=dataStruct.xeps;
Vcol_data=dataStruct.Vcol;
Vela_data=dataStruct.Vela;
Vsmc_data=dataStruct.Vsmc;
iswitch=dataStruct.iswitch;
%Define plot color (black or white background)
colorMode=1;
switch colorMode
    case 1
        figStruct.ColorDef='white';
        figStruct.Color='w';
    case 2
        figStruct.ColorDef='black';
        figStruct.Color='k';
end

%% Access data structure components
V_cent=dataStruct.Cent; %Centroid list
segmentCell=dataStruct.Points; %Lumen boundary coordinates

%% Resampling aorta section contours
%Resample boundary points so each plane has same number of points for lofting
%Find number of points to use based on biggest circumference
d=zeros(size(segmentCell,2),1);
for indNow=1:1:size(segmentCell,2)
    Vs_1=segmentCell{t,indNow}';
    d(indNow)=max(pathLength(Vs_1));
end
nSegment=round(max(d)/pointSpacing);
%Resample
segmentCellSmooth=segmentCell;
segmentCellMean=segmentCell;
w=ones(size(V_cent,1),1); %Cubic smoothing spline weights
indexPlanePoints_V_cent=zeros(1,size(segmentCell,2)); %Indices of centre line points at sections
for indNow=1:1:size(segmentCell,2)
    %Resample section contour
    Vs_1=segmentCell{t,indNow}'; %Current contour
    Vs_1_smooth=evenlySampleCurve(Vs_1,nSegment,smoothFactorSegments,1); %Resample evenly
    Vs_1_mean=mean(Vs_1_smooth,1);
    segmentCellSmooth{t,indNow}=Vs_1_smooth';
    segmentCellMean{t,indNow}=Vs_1_mean;
    %Prepare for center line smoothing by setting weight vector
    [~,indVertex_1]=min(sqrt(sum((V_cent-Vs_1_mean(ones(size(V_cent,1),1),:)).^2,2))); %Index closest to section
    w(indVertex_1)=1e9; %Heigh weight at contour sections
    indexPlanePoints_V_cent(indNow)=indVertex_1; %Store index of closets
end

%% Smooth center line
%Fit smoothing spline through centreline points for loft
if ~isempty(smoothFactorCentreLine)
    V_cent_original=V_cent;
    d=pathLength(V_cent);
    V_cent = csaps(d,V_cent_original',smoothFactorCentreLine,d,w)'; %Smoothed
    cFigure(figStruct); hold on;
    h1=plotV(V_cent_original,'g.-','LineWidth',3,'MarkerSize',35);
    h2=plotV(V_cent,'r-','LineWidth',2,'MarkerSize',15);
    h3=plotV(V_cent(w==max(w),:),'b.','MarkerSize',40);
    legend([h1 h2 h3],{'Original','New','At contours'})
    axisGeom;
    drawnow;
end

%% Offsetting section curves outward if thickening is inward
segmentCellSmooth_pre=segmentCellSmooth;

for q=1:size(segmentCellSmooth,2)
    Vc=segmentCellSmooth{t,q}'; %Current curve vertices
    segmentCellSmooth{t,q}=curveOffset(Vc,wallThickness(q))';
end

%% Visualize offset curves
cFigure(figStruct); hold on;
plotV(V_cent,'b-','LineWidth',2,'MarkerSize',15);
for q=1:size(segmentCellSmooth,2)
    plotV(segmentCellSmooth_pre{t,q}','r.-','LineWidth',3,'MarkerSize',15);
    plotV(segmentCellSmooth{t,q}','g.-','LineWidth',3,'MarkerSize',15);
end
axisGeom;
drawnow;

%% Perform main trunk loft
%

% Initialize figure with center line
hf1=cFigure(figStruct); hold on;
plotV(V_cent,'b.-','LineWidth',3,'markerSize',25);
axisGeom;
camlight headlight;
drawnow;
hp=gobjects(11,1);
F_main_cell=cell(1,size(segmentCell,2)-1); %Faces cell
V_main_cell=cell(1,size(segmentCell,2)-1); %Vertex cell
C_main_gradient_cell=cell(1,size(segmentCell,2)-1); %Color/label dataa
segmentCurve_cell=cell(1,size(segmentCell,2)); %
% Eb_main_cell=cell(1,size(segmentCell,2)-1);
logicBoundary=false(0,0);
maxC=0;
indAdd=0;
%Loft from one plane to next in loop along aorta
for indNow=1:1:(size(segmentCell,2)-1)
    Vs_1=segmentCell{t,indNow}';
    Vs_1=Vs_1(1:end-1,:);
    Vs_2=segmentCell{t,indNow+1}';
    Vs_2(1:end-1,:);
    Vs_1_smooth=segmentCellSmooth{t,indNow}';
    Vs_2_smooth=segmentCellSmooth{t,indNow+1}';
    Vs_1_mean=segmentCellMean{t,indNow};
    Vs_2_mean=segmentCellMean{t,indNow+1};
    d=pathLength(V_cent);
    %Get curve part for lofting
    [~,indVertex_1]=min(sqrt(sum((V_cent-Vs_1_mean(ones(size(V_cent,1),1),:)).^2,2))); %Start
    [~,indVertex_2]=min(sqrt(sum((V_cent-Vs_2_mean(ones(size(V_cent,1),1),:)).^2,2))); %End
    V_cent_part=V_cent(indVertex_1:indVertex_2,:);
    nPart=round(max(pathLength(V_cent_part))/pointSpacing);
    [V_cent_part_smooth] = evenlySampleCurve(V_cent_part,nPart,'spline',0);
    Vs_1_smooth=Vs_1_smooth-Vs_1_mean(ones(size(Vs_1_smooth,1),1),:);
    Vs_1_smooth=Vs_1_smooth+V_cent_part_smooth(ones(size(Vs_1_smooth,1),1),:);
    Vs_2_smooth=Vs_2_smooth-Vs_2_mean(ones(size(Vs_2_smooth,1),1),:);
    Vs_2_smooth=Vs_2_smooth+V_cent_part_smooth(size(V_cent_part_smooth,1)*ones(size(Vs_2_smooth,1),1),:);
    v1=vecnormalize(V_cent_part_smooth(2,:)-V_cent_part_smooth(1,:));
    [Q]=pointSetPrincipalDir(Vs_1_smooth-Vs_1_mean(ones(size(Vs_1_smooth,1),1),:));
    n1=Q(:,3)';
    if dot(v1,n1)<0
        n1=-n1;
    end
    v2=vecnormalize(V_cent_part_smooth(end,:)-V_cent_part_smooth(end-1,:));
    [Q]=pointSetPrincipalDir(Vs_2_smooth-Vs_2_mean(ones(size(Vs_2_smooth,1),1),:));
    n2=Q(:,3)';
    if dot(v2,n2)<0
        n2=-n2;
    end
    plotOn=0;
    nLoft=round(nPart/trunkSegmentReduceFactor);
    [Fs,Vs,Cs]=sweepLoft(Vs_1_smooth,Vs_2_smooth,n1,n2,V_cent_part_smooth,nLoft,0,plotOn);
    indLoftBottom=1:nLoft:size(Vs,1);
    indLoftTop=(nLoft:nLoft:size(Vs,1));
    segmentCurve_cell{indNow}=indLoftBottom+indAdd;
    if indNow==(size(segmentCell,2)-1)
        segmentCurve_cell{indNow+1}=indLoftTop+indAdd;
    end
    indAdd=indAdd+size(Vs,1);
    Eb=patchBoundary(Fs,Vs);
    logicBoundaryNow=false(size(Vs,1),1);
    logicBoundaryNow(unique(Eb(:)))=1;
    F_main_cell{indNow}=Fs;
    V_main_cell{indNow}=Vs;
    %Store color gradient information
    C_main_gradient_cell{indNow}=Cs+maxC;
    maxC=maxC+max(Cs);
    logicBoundary=[logicBoundary; logicBoundaryNow];

    %%
    % Plot Loft of main trunk
    figure(hf1);
    gpatch(Fs,Vs,'kw','kw',0.85);
    plotV(V_cent_part_smooth,'m.-','LineWidth',2,'markerSize',15);
    plotV(Vs_1_smooth,'r-','LineWidth',2);
    plotV(V_cent(indVertex_1,:),'r+','markerSize',25);
    plotV(Vs_2_smooth,'r-','LineWidth',2);
    axis off;
    xlim([140 200]); ylim([120 200]); zlim([20 370]);
    drawnow;
end

%% Material model
%linearly interpolate parameters between planes
% Get center line distance metric for each plane
V_cent_crop=V_cent(indexPlanePoints_V_cent(1):indexPlanePoints_V_cent(end),:);
indexPlanePoints_V_cent_crop=(indexPlanePoints_V_cent-indexPlanePoints_V_cent(1))+1;
d=pathLength(V_cent_crop);
d=d./max(d(:)); %Normalised curve length
curveLengthPlanePoints=d(indexPlanePoints_V_cent_crop);

%mu_data=interp1(curveLengthPlanePoints(materialIndicesSelect),mu_data(materialIndicesSelect),curveLengthPlanePoints,'linear');
Ee_data=interp1(curveLengthPlanePoints(materialIndicesSelect),Ee_data(materialIndicesSelect),curveLengthPlanePoints,'linear');
xeps_data=interp1(curveLengthPlanePoints(materialIndicesSelect),xeps_data(materialIndicesSelect),curveLengthPlanePoints,'linear');
Vcol_data=interp1(curveLengthPlanePoints(materialIndicesSelect),Vcol_data(materialIndicesSelect),curveLengthPlanePoints,'linear');
Vela_data=interp1(curveLengthPlanePoints(materialIndicesSelect),Vela_data(materialIndicesSelect),curveLengthPlanePoints,'linear');

%% Main trunk
% Join element sets to form main trunk and form colour information
C_main=[];
for q=1:1:numel(F_main_cell)
    C_main=[C_main; q*ones(size(F_main_cell{q},1),1)];
end
[F_main,V_main,C_path]=joinElementSets(F_main_cell,V_main_cell,C_main_gradient_cell);
%Scale color gradient for wall thickness interpolation
C_path = rescale(C_path,1,numThicknessLevels);
F_main=fliplr(F_main); %Invert orientation
[F_main,V_main,ind1,indFix]=mergeVertices(F_main,V_main);
for q=1:1:numel(segmentCurve_cell)
    segmentCurve_cell{q}=indFix(segmentCurve_cell{q});
end
logicBoundary=logicBoundary(ind1);

% Perform Smoothing on main trunk
controlParameter.n=numSmoothTrunk;
controlParameter.Method='HC';
controlParameter.RigidConstraints=find(logicBoundary); %Ensure smoothing cannot change coordinates of planes of interest
[V_main]=patchSmooth(F_main,V_main,[],controlParameter);
% Plot
cFigure(figStruct);
subplot(1,2,1); hold on;
title('Segments')
plotV(V_cent,'k.-','LineWidth',3,'markerSize',25);
gpatch(F_main,V_main,C_main,'k');
patchNormPlot(F_main,V_main,1,'v'); %Check normals of elements all facing outward
for q=1:1:numel(segmentCurve_cell)
    plotV(V_main(segmentCurve_cell{q},:),'b.-','LineWidth',5,'markerSize',35);
end
plotV(V_main(logicBoundary,:),'r.','markerSize',25);
axisGeom;
colormap(gca,gjet(250)); icolorbar;
camlight headlight;
lighting gouraud;
subplot(1,2,2); hold on;
title('Continuous')
plotV(V_cent,'k.-','LineWidth',3,'markerSize',25);
gpatch(F_main,V_main,C_path,'k');
plotV(V_main(logicBoundary,:),'r.','markerSize',25);
% patchNormPlot(F_main,V_main,1,'v');
axisGeom;
colormap(gca,gjet(250)); colorbar;
camlight headlight;
lighting gouraud;
drawnow;

%%

E_rings=[F_main(:,[1 2]); fliplr(F_main(:,[3 4]))];
E_rings=uniqueIntegerRow(E_rings);

%% Load Iliac Branch Data
%Right Origin
V_R_ori=dataStruct.R_Ori(1:end-1,:);
V_R_ori=resampleCurve(V_R_ori,pointSpacing,1);
V_R_ori(:,1)=V_R_ori(:,1);
V_R_ori=curveOffset(V_R_ori,wallThickness(end));

%Left Origin
V_L_ori=dataStruct.L_Ori(1:end-1,:);
V_L_ori=resampleCurve(V_L_ori,pointSpacing,1);
V_L_ori(:,1)=V_L_ori(:,1)-1;
V_L_ori=curveOffset(V_L_ori,wallThickness(end));

%Left End
V_L_Iliac=dataStruct.L_Iliac(1:end-1,:);
V_L_Iliac=resampleCurve(V_L_Iliac,pointSpacing,1);
V_L_Iliac=curveOffset(V_L_Iliac,wallThickness(end));

%Right End
V_R_Iliac=dataStruct.R_Iliac(1:end-1,:);
V_R_Iliac=resampleCurve(V_R_Iliac,pointSpacing,1);
V_R_Iliac=curveOffset(V_R_Iliac,wallThickness(end));

%Bifurcation Data
V_bifurc=dataStruct.bifurc;
V_bifurc=curveOffset(V_bifurc,wallThickness(end));

%% Resample curve to have same number of points as main trunk for loft
Eb_all=patchBoundary(F_main,V_main);
Eb_lowerSegment=patchBoundary(F_main(C_main==max(C_main),:),V_main);
Eb_lower=Eb_lowerSegment(all(ismember(Eb_lowerSegment,Eb_all),2),:);
indLowerCurve=edgeListToCurve(Eb_lower);
indLowerCurve=indLowerCurve(1:end-1);
indLowerCurve=indLowerCurve(:);
V_main_lowerCurve=V_main(indLowerCurve,:);
if isPolyClockwise(V_bifurc)~=isPolyClockwise(V_main_lowerCurve)
    V_bifurc=flipud(V_bifurc);
end
V_bifurc=evenlySampleCurve(V_bifurc,size(V_main_lowerCurve,1),'spline',1); %resample
[~,indMin]=minDist(V_main_lowerCurve(1,:),V_bifurc);
if indMin>1
    V_bifurc=[V_bifurc(indMin:end,:); V_bifurc(1:indMin-1,:)];
end
% Mesh transition region
cPar.closeLoopOpt=1;
cPar.patchType='quad';
[F_add,V_add,indStart,indEndBifurc]=polyLoftLinear(V_main_lowerCurve,V_bifurc,cPar);
F_add=fliplr(F_add);
C_add=ones(size(F_add,1),1);
F_main=[F_main; F_add+size(V_main,1)];
indEndBifurc=indEndBifurc+size(V_main,1);
V_main=[V_main; V_add];
C_main=[C_main; C_add+max(C_main)];
C_path=[C_path; thicknessIndexTransitionRegion.*ones(size(C_add))];
[F_main,V_main,~,indFix]=mergeVertices(F_main,V_main);
indLowerCurve=indFix(indLowerCurve);
indEndBifurc=indFix(indEndBifurc);
E_rings=indFix(E_rings);
for q=1:1:numel(segmentCurve_cell)
    segmentCurve_cell{q}=indFix(segmentCurve_cell{q});
end

%% Perform Smoothing on Transition Region
indTouch=indLowerCurve(:);
numSmoothGrowSteps=3;
for q=1:1:numSmoothGrowSteps
    logicFacesTouch=any(ismember(F_main,indTouch),2);
    indTouch=F_main(logicFacesTouch,:);
end
smoothControlParameters.n=numSmoothPants_HC;
smoothControlParameters.Method='HC';
smoothControlParameters.RigidConstraints=[unique(F_main(~logicFacesTouch,:)); indLowerCurve];
[V_main]=patchSmooth(F_main,V_main,[],smoothControlParameters);

%%
% Plot Main and Transition Region
cFigure(figStruct); hold on;
gpatch(F_main,V_main,logicFacesTouch,'none');
patchNormPlot(F_main,V_main);
plotV(V_main(indLowerCurve,:),'r.-','markerSize',25);
plotV(V_main(indEndBifurc,:),'g.-','markerSize',25);
for q=1:1:numel(segmentCurve_cell)
    plotV(V_main(segmentCurve_cell{q},:),'b.-','LineWidth',5,'markerSize',35);
end
axisGeom;
colormap gjet;
camlight headlight;
lighting flat;
drawnow;

%% Pinched ellipse to two circle split
ns=ceil(abs(mean(V_main(indEndBifurc,3))-mean(V_R_ori(:,3)))/pointSpacing)+1;
V_cell={V_main(indEndBifurc,:),V_R_ori,V_L_ori};
patchType='quad';
splitMethod='nearMid';
controlParSmooth.Method='LAP';
controlParSmooth.n=numSmoothPants_LAP;
[F_split,V_split,curveIndices,C_split]=splitCurveSetMesh(V_cell,ns,patchType,controlParSmooth,splitMethod,1);

%%
% Plot
cFigure(figStruct); hold on;
% gpatch(F_main,V_main,'kw','k');
gpatch(F_split,V_split,C_split,'k');
patchNormPlot(F_split,V_split);
plotV(V_main(indEndBifurc,:),'r.-','markerSize',25);
% plotV(V_split(curveIndices{1},:),'g.-','markerSize',25);
axisGeom;
colormap gjet;
camlight headlight;
drawnow;

%% Ensure all points for curves are clockwise for loft
indEndBifurc_split=curveIndices{1}; indEndBifurc_split=indEndBifurc_split(:);
indBranch11=curveIndices{2}; indBranch11=indBranch11(:);
indBranch21=curveIndices{3}; indBranch21=indBranch21(:);
V_branch12=evenlySampleCurve(V_L_Iliac,numel(indBranch11),'spline',1);
V_branch22=evenlySampleCurve(V_R_Iliac,numel(indBranch21),'spline',1);
V_branch11=V_split(indBranch11,:);
if ~isPolyClockwise(V_branch11)
    indBranch11=flipud(indBranch11);
    V_branch11=V_split(indBranch11,:);
end
V_branch21=V_split(indBranch21,:);
if ~isPolyClockwise(V_branch21)
    indBranch21=flipud(indBranch21);
    V_branch21=V_split(indBranch21,:);
end
if ~isPolyClockwise(V_branch22)
    V_branch22=flipud(V_branch22);
end
if ~isPolyClockwise(V_branch12)
    V_branch12=flipud(V_branch12);
end
[~,indMin]=min(V_branch12(:,1));
if indMin>1
    V_branch12=[V_branch12(indMin:end,:); V_branch12(1:indMin-1,:)];
end
[~,indMin]=min(V_branch22(:,1));
if indMin>1
    V_branch22=[V_branch22(indMin:end,:); V_branch22(1:indMin-1,:)];
end

%% Define points for iliac extensions and Loft
V_loft_cell1{1}=V_branch11;
V_loft_cell2{1}=V_branch12;
V_loft_cell1{2}=V_branch21;
V_loft_cell2{2}=V_branch22;
V_loft_cent_cell{1}=dataStruct.Cent_I_R;
V_loft_cent_cell{2}=dataStruct.Cent_I_L;
F_iliac_cell=cell(1,2);
V_iliac_cell=cell(1,2);
for q=1:1:2
    V1=V_loft_cell1{q};
    V2=V_loft_cell2{q};
    Vc=V_loft_cent_cell{q};
    [~,indMin]=minDist(mean(V1,1),Vc);
    if indMin>1
        Vc=Vc(indMin:end,:);
    end
    Vc(1,:)=mean(V1,1); %Overide first point with mean of segment
    np=ceil(max(pathLength(Vc))/pointSpacing)+1;
    Vc = evenlySampleCurve(Vc,np,'spline',0);

    %%
    % plot loft paths and profile
    cFigure(figStruct); hold on;
    plotV(V1,'b.-','markerSize',25);
    plotV(V2,'r.-','markerSize',25);
    plotV(Vc,'g.-','markerSize',25);
    axisGeom;
    colormap gjet;
    camlight headlight;
    lighting gouraud;
    drawnow;

    %%
    V2 = evenlySampleCurve(V2,size(V1,1),'spline',1);
    v1=vecnormalize(Vc(2,:)-Vc(1,:));
    [Q]=pointSetPrincipalDir(V1-Vc(ones(size(V1,1),1),:));
    n1=Q(:,3)';
    if dot(v1,n1)<0
        n1=-n1;
    end
    v2=vecnormalize(Vc(end,:)-Vc(end-1,:));
    [Q]=pointSetPrincipalDir(V2-Vc(size(Vc,1)*ones(size(V2,1),1),:));
    n2=Q(:,3)';
    if dot(v2,n2)<0
        n2=-n2;
    end
    pointSpacingNow=mean(diff(pathLength(V1)));
    np=ceil(max(pathLength(Vc))/pointSpacingNow);
    Vc = evenlySampleCurve(Vc,np,'spline',0);
    [F_loft,V_loft,C_loft]=sweepLoft(V1,V2,n1,n2,Vc,size(Vc,1),0,0);
    F_iliac_cell{q}=F_loft;
    V_iliac_cell{q}=V_loft;
end
F_branch1=F_iliac_cell{1};
V_branch1=V_iliac_cell{1};
F_branch2=F_iliac_cell{2};
V_branch2=V_iliac_cell{2};

%%
% Plot Iliac extension Loft
cFigure(figStruct); hold on;
% gpatch(F_main,V_main,C_main,'k');
% gpatch(F_main(C_main==max(C_main),:),V_main,'kw','none');
gpatch(F_split,V_split,C_split,'k');
gpatch(F_branch1,V_branch1,'gw');
patchNormPlot(F_branch1,V_branch1);
gpatch(F_branch2,V_branch2,'bw');
patchNormPlot(F_branch2,V_branch2);
plotV(V_main(indEndBifurc,:),'r.-','markerSize',25);
plotV(V_branch11,'g.-','markerSize',25);
plotV(V_branch12,'g.-','markerSize',25);
plotV(V_branch21,'b.-','markerSize',25);
plotV(V_branch22,'b.-','markerSize',25);
axisGeom;
colormap gjet;
camlight headlight;
lighting gouraud;
drawnow;

%% Join Iliac extensions to Bifurcation
[Fp,Vp,Cp]=joinElementSets({F_split,F_branch1,F_branch2},{V_split,V_branch1,V_branch2});
[Fp,Vp,~,indFix]=mergeVertices(Fp,Vp);
Eb_p=patchBoundary(Fp,Vp);
indEndBifurc_split=indFix(indEndBifurc_split);
indBranch11=indFix(indBranch11);
indBranch21=indFix(indBranch21);
%%Smooth Iliacs
indTouch=[indBranch11(:); indBranch21(:)];
numSmoothGrowSteps=3;
for q=1:1:numSmoothGrowSteps
    logicFacesTouch=any(ismember(Fp,indTouch),2);
    indTouch=Fp(logicFacesTouch,:);
end
smoothControlParameters.n=numSmoothPants_HC;
smoothControlParameters.Method='HC';
smoothControlParameters.RigidConstraints=[unique(Fp(~logicFacesTouch,:)); indEndBifurc_split];
[Vp]=patchSmooth(Fp,Vp,[],smoothControlParameters);

%%
% Plot Iliacs
cFigure(figStruct); hold on;
gpatch(Fp,Vp,logicFacesTouch,'k');
patchNormPlot(Fp,Vp);
plotV(Vp(indBranch11,:),'r-','LineWidth',3)
plotV(Vp(indBranch21,:),'r-','LineWidth',3)
axisGeom;
colormap gjet;
camlight headlight;
lighting gouraud;
drawnow;

%%
% Add to main trunk
F_main=[F_main; Fp+size(V_main,1)];
indBranch11=indBranch11+size(V_main,1);
indBranch21=indBranch21+size(V_main,1);
V_main=[V_main; Vp];
C_main=[C_main; Cp+max(C_main)];
C_path=[C_path; thicknessIndexBifurcation.*ones(size(Cp))];
C_path_mat=C_path; %Copy over path color data for material assignments
[F_main,V_main,~,indFix]=mergeVertices(F_main,V_main);
indLowerCurve=indFix(indLowerCurve);
indBranch11=indFix(indBranch11);
indBranch21=indFix(indBranch21);
E_rings=indFix(E_rings);

for q=1:1:numel(segmentCurve_cell)
    segmentCurve_cell{q}=indFix(segmentCurve_cell{q});
end
%
indTouch=[indLowerCurve(:); indBranch11(:); indBranch21(:)];
nGrowthSteps=round(distSmoothGrowth/pointSpacing);
for q=1:1:numSmoothGrowSteps
    logicFacesTouch=any(ismember(F_main,indTouch),2);
    indTouch=F_main(logicFacesTouch,:);
end
ind1=F_main(~logicFacesTouch,:);
indRigid=unique([ind1(:); indBranch11(:); indBranch21(:)]);
smoothControlParameters.Tolerance=0.01;
smoothControlParameters.n=numSmoothPants_HC;
smoothControlParameters.Method='HC';
smoothControlParameters.RigidConstraints=indRigid;
[V_main]=patchSmooth(F_main,V_main,[],smoothControlParameters);
%
cFigure(figStruct); hold on;
gpatch(F_main,V_main,logicFacesTouch,'k');
% patchNormPlot(F_main,V_main);
plotV(V_main(indLowerCurve,:),'r.-','markerSize',25);
plotV(V_main(indBranch11,:),'r.-','markerSize',25);
plotV(V_main(indBranch21,:),'r.-','markerSize',25);
for q=1:1:numel(segmentCurve_cell)
    plotV(V_main(segmentCurve_cell{q},:),'b.-','LineWidth',5,'markerSize',35);
end
axisGeom;
colormap gjet;
camlight headlight;
lighting gouraud;
drawnow;

%% Load Branch data

%Other branches
V_L_Renal_Ori=dataStruct.L_Renal_Ori(1:end-1,:);
V_L_Renal_Ori=resampleCurve(V_L_Renal_Ori,pointSpacing,1);
V_L_Renal_Ori=curveOffset(V_L_Renal_Ori,wallThickness(7));
V_L_Renal=dataStruct.L_Renal(1:end-1,:);
V_L_Renal=resampleCurve(V_L_Renal,pointSpacing,1);
V_L_Renal=curveOffset(V_L_Renal,wallThickness(7));
V_Cent_L_Renal=dataStruct.Cent_Renal_L;
V_Cent_L_Renal=resampleCurve(V_Cent_L_Renal,pointSpacing,0);

V_R_Renal_Ori=dataStruct.R_Renal_Ori(1:end-1,:);
V_R_Renal_Ori=resampleCurve(V_R_Renal_Ori,pointSpacing,1);
V_R_Renal_Ori=curveOffset(V_R_Renal_Ori,wallThickness(7));

V_R_Renal=dataStruct.R_Renal(1:end-1,:);
V_R_Renal=resampleCurve(V_R_Renal,pointSpacing,1);
V_R_Renal=curveOffset(V_R_Renal,wallThickness(7));
V_Cent_R_Renal=dataStruct.Cent_Renal_R;
V_Cent_R_Renal=resampleCurve(V_Cent_R_Renal,pointSpacing,0);

V_SMA_Ori=dataStruct.SMA_O(1:end-1,:);
V_SMA_Ori=resampleCurve(V_SMA_Ori,pointSpacing,1);
V_SMA_Ori=curveOffset(V_SMA_Ori,wallThickness(7));

V_SMA=dataStruct.SMA(1:end-1,:);
V_SMA=resampleCurve(V_SMA,pointSpacing,1);
V_SMA=curveOffset(V_SMA,wallThickness(7));
V_Cent_SMA=dataStruct.Cent_SMA;
V_Cent_SMA=resampleCurve(V_Cent_SMA,pointSpacing,0);
V_COE_Ori=dataStruct.COE_O(1:end-1,:);
V_COE_Ori=resampleCurve(V_COE_Ori,pointSpacing,1);
V_COE_Ori=curveOffset(V_COE_Ori,wallThickness(7));
V_COE=dataStruct.COE(1:end-1,:);
V_COE=resampleCurve(V_COE,pointSpacing,1);
V_COE=curveOffset(V_COE,wallThickness(7));
V_Cent_COE=dataStruct.Cent_Coeliac;
V_Cent_COE=resampleCurve(V_Cent_COE,pointSpacing,0);
V_BCA_Ori=dataStruct.BCA_O(1:end-1,:);
V_BCA_Ori=resampleCurve(V_BCA_Ori,pointSpacing,1);
V_BCA_Ori=curveOffset(V_BCA_Ori,wallThickness(2));
V_BCA=dataStruct.BCA(1:end-1,:);
V_BCA=resampleCurve(V_BCA,pointSpacing,1);
V_BCA=curveOffset(V_BCA,wallThickness(2));
V_Cent_BCA=dataStruct.Cent_BCA;
V_Cent_BCA=resampleCurve(V_Cent_BCA,pointSpacing,0);
V_LCCA_Ori=dataStruct.LCCA_O(1:end-1,:);
V_LCCA_Ori=resampleCurve(V_LCCA_Ori,pointSpacing,1);
V_LCCA_Ori=curveOffset(V_LCCA_Ori,wallThickness(2));
V_LCCA=dataStruct.LCCA(1:end-1,:);
V_LCCA=resampleCurve(V_LCCA,pointSpacing,1);
V_LCCA=curveOffset(V_LCCA,wallThickness(2));
V_Cent_LCCA=dataStruct.Cent_LCCA;
V_Cent_LCCA=resampleCurve(V_Cent_LCCA,pointSpacing,0);
V_LSA_Ori=dataStruct.LSA_O(1:end-1,:);
V_LSA_Ori=resampleCurve(V_LSA_Ori,pointSpacing,1);
V_LSA_Ori=curveOffset(V_LSA_Ori,wallThickness(2));
V_LSA=dataStruct.LSA(1:end-1,:);
V_LSA=resampleCurve(V_LSA,pointSpacing,1);
V_LSA=curveOffset(V_LSA,wallThickness(2));
V_Cent_LSA=dataStruct.Cent_LSA;
V_Cent_LSA=resampleCurve(V_Cent_LSA,pointSpacing,0);

% Plot
cFigure(figStruct); hold on;
gpatch(F_main,V_main,C_main,0.65);
plotV(V_R_Renal_Ori,'r.-','markerSize',5);
plotV(V_R_Renal,'w.-','markerSize',5);
plotV(V_L_Renal_Ori,'g.-','markerSize',5);
plotV(V_L_Renal,'w.-','markerSize',5);
plotV(V_COE_Ori,'w.-','markerSize',5);
plotV(V_COE,'w.-','markerSize',5);
plotV(V_SMA_Ori,'w.-','markerSize',5);
plotV(V_SMA,'w.-','markerSize',5);
plotV(V_BCA_Ori,'w.-','markerSize',5);
plotV(V_BCA,'w.-','markerSize',5);
plotV(V_LCCA_Ori,'w.-','markerSize',5);
plotV(V_LCCA,'w.-','markerSize',5);
plotV(V_LSA_Ori,'w.-','markerSize',5);
plotV(V_LSA,'w.-','markerSize',5);
xlim([140 200]); ylim([120 200]); zlim([20 370]);
axisGeom;
colormap jet;
camlight headlight;
lighting gouraud;
drawnow;

%% Perform Extrude cut
% Find Centroid of Ori and shoot vector in direction of nearest Centreline
% point. Apply same vector to each point on Ori and every element on the
% main trunk wall it touches is deleted. Loft then from deleted elements of
% main to each branch ori
V_cut=V_R_Renal_Ori;
V_cut_cell={V_R_Renal_Ori,V_L_Renal_Ori,V_SMA_Ori,V_COE_Ori,...
    V_LSA_Ori,V_LCCA_Ori,V_BCA_Ori};
%Loop over connections
smoothControlParameters.n=numSmoothBranchAttachments;
smoothControlParameters.Method='HC';
hw = waitbar(0,'Please wait...');
numSteps=numel(V_cut_cell);
V_endCurve_cell=cell(1,numSteps);
C_main_max=max(C_main)+1;
for q=1:1:numSteps
    waitbar(q/numSteps,hw,['Processing case ',num2str(q),' of ',num2str(numSteps)]);
    VF=patchCentre(F_main,V_main); %Face centre coordinates
    [~,indMin]=minDist(mean(V_cut_cell{q},1),VF);
    c=C_path(indMin);
    [F_main,V_main,C_main,indEnd,logicRemoveFaces,segmentCurve_cell,E_rings]=circleCutExtrude(F_main,V_main,C_main,V_cent,V_cut_cell{q},pointSpacing,0,smoothControlParameters,segmentCurve_cell,E_rings);
    numFacesNewFeature=nnz(C_main==max(C_main));
    C_path=[C_path(~logicRemoveFaces); thicknessIndicesBranches(q)*ones(numFacesNewFeature,1)];
    C_path_mat=[C_path_mat(~logicRemoveFaces); c*ones(numFacesNewFeature,1)];
    V_endCurve_cell{q}=V_main(indEnd,:);
end
close(hw);
% plot
cFigure(figStruct); hold on;
gpatch(F_main,V_main,C_path,'k',1);
for q=1:1:numel(V_endCurve_cell)
    plotV(V_endCurve_cell{q},'m.-','markerSize',25);
end
axisGeom;
colormap gjet; colorbar;
camlight headlight;
drawnow;
% plot
cFigure(figStruct); hold on;
gpatch(F_main,V_main,C_path_mat,'k',1);
axisGeom;
colormap gjet; colorbar;
camlight headlight;
drawnow;

% Resample branch ends and sweep
V_loft_cell1=V_endCurve_cell;
V_loft_cell2={V_R_Renal,V_L_Renal,V_SMA,V_COE,...
    V_LSA,V_LCCA,V_BCA};
V_loft_cent_cell={V_Cent_R_Renal,V_Cent_L_Renal,V_Cent_SMA,V_Cent_COE,V_Cent_LSA,V_Cent_LCCA,V_Cent_BCA};
F_branch_cell=cell(1,numel(V_endCurve_cell));
V_branch_cell=cell(1,numel(V_endCurve_cell));
indBranchTop_cell=cell(1,numel(V_endCurve_cell));
indBranchBottom_cell=cell(1,numel(V_endCurve_cell));
for q=1:1:numel(V_endCurve_cell)
    V1=V_loft_cell1{q};
    pointSpacingNow=mean(diff(pathLength(V1)));
    V2=V_loft_cell2{q};
    Vc=V_loft_cent_cell{q};
    np=ceil(max(pathLength(Vc))/pointSpacingNow);
    Vc = evenlySampleCurve(Vc,np,'spline',0);
    V2 = evenlySampleCurve(V2,size(V1,1),'spline',1);
    v1=vecnormalize(Vc(2,:)-Vc(1,:));
    [Q]=pointSetPrincipalDir(V1-Vc(ones(size(V1,1),1),:));
    n1=Q(:,3)';
    if dot(v1,n1)<0
        n1=-n1;
    end
    v2=vecnormalize(Vc(end,:)-Vc(end-1,:));
    [Q]=pointSetPrincipalDir(V2-Vc(size(Vc,1)*ones(size(V2,1),1),:));
    n2=Q(:,3)';
    if dot(v2,n2)<0
        n2=-n2;
    end
    b2=vecnormalize(V2(2,:)-V2(1,:));
    a2=vecnormalize(V2(1,:)-mean(V2,1));
    c2=cross(b2,a2);
    if dot(n2,c2)<0
        V2=flipud(V2);
    end
    [F_loft,V_loft,C_loft]=sweepLoft(V1,V2,n1,n2,Vc,size(Vc,1),0,0);
    F_loft=fliplr(F_loft);
    indLoftBottom=1:size(Vc,1):size(V_loft,1);
    indLoftTop=(size(Vc,1):size(Vc,1):size(V_loft,1));
    F_branch_cell{q}=F_loft;
    V_branch_cell{q}=V_loft;
    if q>1
        sizAll=cellfun(@(x) size(x,1),V_branch_cell);
        indBranchTop_cell{q}=indLoftTop+sum(sizAll(1:q-1));
        indBranchBottom_cell{q}=indLoftBottom+sum(sizAll(1:q-1));
    else
        indBranchTop_cell{q}=indLoftTop;
        indBranchBottom_cell{q}=indLoftBottom;
    end
    % The below plot highlights each loft path in a loop
    %     cFigure(figStruct); hold on;
    %     title(num2str(q))
    %     gpatch(F_main,V_main,'kw','none',0.5);
    %     gpatch(F_loft,V_loft,'rw','r',1);
    %     plotV(V_loft(indLoftBottom,:),'b.-','LineWidth',3);
    %     plotV(V_loft(indLoftTop,:),'g.-','LineWidth',3);
    % %     plotV(V1,'r.-','LineWidth',3);
    % %     plotV(V2,'b.-','LineWidth',3);
    % %     plotV(Vc,'k.-','LineWidth',3);
    % %     quiverVec(Vc(1,:),n1,10,'k');
    % %     quiverVec(Vc(end,:),n2,10,'k');
    %     axisGeom;
    %     camlight headlight;
    %     drawnow;
end
%Join branches together
VF=patchCentre(F_main,V_main);
%Joining branches
[F_branch,V_branch,C_branch]=joinElementSets(F_branch_cell,V_branch_cell);
%Adding branches to main thing
numVerticesInitial=size(V_main,1);
F_main=[F_main; F_branch+numVerticesInitial];
V_main=[V_main; V_branch];
C_main=[C_main; C_branch+max(C_main)];
C_path_branch=thicknessIndicesBranches(C_branch);
C_path=[C_path; C_path_branch(:)];

%%
% plot
cFigure(figStruct); hold on;
gpatch(F_main,V_main,'kw','none',0.5);
axisGeom;
colormap gjet; colorbar;
camlight headlight;
drawnow;
for q=1:1:numel(V_cut_cell)
    V_mean_now=mean(V_cut_cell{q},1);
    plotV(V_mean_now,'r.','markerSize',25);
    [~,indMin]=minDist(V_mean_now,VF);
    plotV(VF(indMin,:),'b.','markerSize',25);
    c=C_path_mat(indMin);
    C_path_mat=[C_path_mat; c*ones(size(F_branch_cell{q},1),1)];
end

%%

%Fix curve indices for joining sets
for q=1:1:numel(indBranchBottom_cell)
    indBranchTop_cell{q}=indBranchTop_cell{q}+numVerticesInitial;
    indBranchBottom_cell{q}=indBranchBottom_cell{q}+numVerticesInitial;
end
%merging nodes together
[F_main,V_main,~,indFix]=mergeVertices(F_main,V_main);
%Fix curve indices for merging
for q=1:1:numel(indBranchBottom_cell)
    indBranchTop_cell{q}=indFix(indBranchTop_cell{q})';
    indBranchBottom_cell{q}=indFix(indBranchBottom_cell{q})';
end
for q=1:1:numel(segmentCurve_cell)
    segmentCurve_cell{q}=indFix(segmentCurve_cell{q});
end

E_rings=indFix(E_rings);

% Snapping material color data to number of materials
C_path_index=C_path_mat;
C_path_mat=round(rescale(C_path_mat,1,numMaterials));

%%
% plot
cFigure(figStruct); hold on;
gpatch(F_main,V_main,C_path,'k',1);
axisGeom;
colormap gjet; colorbar;
camlight headlight;
drawnow;

%%
%Smoothing of Branches
indTouch=[indBranchBottom_cell{:}];
numSmoothGrowSteps=3;
for q=1:1:numSmoothGrowSteps
    logicFacesTouch=any(ismember(F_main,indTouch),2);
    indTouch=F_main(logicFacesTouch,:);
end

smoothControlParameters.n=numSmoothBranches;
smoothControlParameters.Method='HC';
smoothControlParameters.RigidConstraints=unique(F_main(~logicFacesTouch,:));
[V_main]=patchSmooth(F_main,V_main,[],smoothControlParameters);

% Interpolating thicknesses
thicknessData=dataStruct.WallThickness; %Thickness data
indexData=1:1:numel(thicknessData); %Index data for x-axis for interpolation
C_thickness=interp1(indexData,thicknessData,C_path,'spline');

%%
% plot
hf=cFigure(figStruct); hold on;
gtitle('Wall Thickness')
gpatch(F_main,V_main,C_thickness,'k',1);
axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;
drawnow;

%% Inverting offset direction
%Converting face thickness data to vertex data
nodalThickness=faceToVertexMeasure(F_main,V_main,C_thickness);
[~,~,N]=patchNormal(F_main,V_main);

%%
% plot
cFigure(figStruct); hold on;
gpatch(F_main,V_main,'kw','kw',0.5);
title('Mesh Offset')
quiverVec(V_main,-N.*nodalThickness,[],nodalThickness);
axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;
drawnow;

%% Thicken to create hexahedral elements
[ET,VT,Fp1,Fp2]=patchThick(F_main,V_main,-1,nodalThickness,numThickenSteps);
CT=repmat((1:1:size(F_main,1))',numThickenSteps,1);

ET=ET(:,[5:8 1:4]);
FT_inner=Fp2;
FT_outer=Fp1;

[~,logicPositive]=hexVol(ET,VT);
% logicPositive
if any(logicPositive==0)
    error('Negative hex volume found');
end

C_ET_path_mat_index=C_path_mat;
C_ET_path_index=C_path_index;
indicesNodesInner=unique(FT_inner(:));
% Find elements touching the branch ends
indBranchTopAll=[indBranchTop_cell{:}]; %+size(V_main,1)*numThickenSteps;
logicBranchEndElement=any(ismember(ET,indBranchTopAll),2);

%% Retrieve segment curve indices
segmentCurve_cell_outer=segmentCurve_cell;
for q=1:1:numel(segmentCurve_cell)
    segmentCurve_cell{q}=segmentCurve_cell{q}+size(V_main,1)*numThickenSteps;
end

[FT,CFT]=element2patch(ET,CT);

%%
% plot
cFigure(figStruct); hold on;
gpatch(FT,VT,CFT,'k',0.5);
gpatch(FT_inner,VT,'g','k',1);

for q=1:1:numel(segmentCurve_cell)
    plotV(VT(segmentCurve_cell{q},:),'g.-','LineWidth',5,'markerSize',15);
end

for q=1:1:numel(segmentCurve_cell)
    plotV(VT(segmentCurve_cell_outer{q},:),'r.-','LineWidth',5,'markerSize',15);
end

axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;
drawnow;

%% Grouping ring edges and offset so they are on the inside

%Offset indices inward
E_rings=E_rings+size(V_main,1)*numThickenSteps;

%Group indices to form rings
optionStruct.outputType='label';
G_rings=tesgroup(E_rings,optionStruct);

%Compose ring cell
ringCurve_cell=cell(1,max(G_rings));
for q=1:1:max(G_rings)
    indList=edgeListToCurve(E_rings(G_rings==q,:)); %1=Ring number
    indList=indList(1:end-1);
    ringCurve_cell{q}=indList;
end

%%
% plot
cFigure(figStruct); hold on;
gpatch(FT,VT,'kw','none',0.5);

plotColors=gjet(max(G_rings));
for q=1:1:max(G_rings)
    hp=plotV(VT(ringCurve_cell{q},:),'b.-','LineWidth',5,'markerSize',15);
%     hp.Color=plotColors(q,:);
end

axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;
drawnow;

%% Create color data for hex elements

C_ET_path_mat_index=C_ET_path_mat_index(CT);
C_ET_path_index=C_ET_path_index(CT);
logicBranchEndElement=logicBranchEndElement(CT); %Colors for original element indices (and sweeping steps)
%%Define Inner Surface Set
logicElementsInner=any(ismember(ET,indicesNodesInner),2);
indicesElementsInner=find(logicElementsInner);

%%
% plot
cFigure(figStruct); hold on;
gpatch(FT_inner,VT,'gw','none',0.5);
gpatch(FT_outer,VT,'rw','none',0.5);
patchNormPlot(FT_inner,VT)
axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;
drawnow;

%% Material Properties
% Assign element material parameters
% Derive interpolatable parameters [Vcol,Vela,Ee,xeps,mu]
indexData=1:1:numel(xeps_data);
indexDataInterp=linspace(1,numel(xeps_data),numMaterials);
indexData_ET=1:1:numMaterials;

xeps_vector=interp1(indexData,xeps_data,indexDataInterp,'spline');
xeps_ET=interp1(indexData_ET,xeps_vector,C_ET_path_mat_index,'spline');
Ee_vector=interp1(indexData,Ee_data,indexDataInterp,'spline');
Ee_ET=interp1(indexData_ET,Ee_vector,C_ET_path_mat_index,'spline');
Vela_vector=interp1(indexData,Vela_data,indexDataInterp,'spline');
Vela_ET=interp1(indexData_ET,Vela_vector,C_ET_path_mat_index,'spline');
Vcol_vector=interp1(indexData,Vcol_data,indexDataInterp,'spline');
Vcol_ET=interp1(indexData_ET,Vcol_vector,C_ET_path_mat_index,'spline');
Vsmc_vector=interp1(indexData,Vsmc_data,indexDataInterp,'spline');
Vsmc_ET=interp1(indexData_ET,Vsmc_vector,C_ET_path_mat_index,'spline');

%% Define Branch Ends Set
[FT,CFT]=element2patch(ET,logicElementsInner);
[~,CFT_logicBranchEndElement]=element2patch(ET,logicBranchEndElement);
[~,CFT_path_mat]=element2patch(ET,C_ET_path_mat_index);
[~,xeps_FT]=element2patch(ET,xeps_ET);
[~,Ee_FT]=element2patch(ET,Ee_ET);
[~,Vela_FT]=element2patch(ET,Vela_ET);
[~,Vcol_FT]=element2patch(ET,Vcol_ET);
[~,Vsmc_FT]=element2patch(ET,Vsmc_ET);
indBoundary=tesBoundary(FT,VT);
Fb=FT(indBoundary,:);
F1=sort(ET(:,[1 2 3 4]),2);
F2=sort(ET(:,[5 6 7 8]),2);
FT_boundary=FT(indBoundary,:);
sizVirt=size(VT,1)*ones(1,size(FT,2));
indVirt_F1=sub2indn(sizVirt,F1);
indVirt_F2=sub2indn(sizVirt,F2);
indVirt_FT_boundary=sub2indn(sizVirt,sort(FT_boundary,2));
indVirt_FT=sub2indn(sizVirt,sort(FT,2));
CFT_logicBranchEndElement=CFT_logicBranchEndElement & ismember(indVirt_FT,indVirt_FT_boundary); %Only keep boundary members
CFT_logicBranchEndElement=CFT_logicBranchEndElement & ~ismember(indVirt_FT,indVirt_F1); %Cant be member of top
CFT_logicBranchEndElement=CFT_logicBranchEndElement & ~ismember(indVirt_FT,indVirt_F2); %Cant be member of bottom

%Nodes for boundary conditions
indNodesFix=FT(CFT_logicBranchEndElement,:);
indNodesFix=unique(indNodesFix(:));

%%
% plot full mesh
cFigure(figStruct); hold on;
gpatch(FT,VT,CFT_logicBranchEndElement,'r',1);
xlim([140 200]); ylim([120 200]); zlim([20 370]);
axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;
drawnow;

%%
% plot interpolated parameters

cFigure(figStruct);
subplot(2,3,1);hold on;
title('Indexing color');
gpatch(FT(indBoundary,:),VT,CFT_path_mat(indBoundary,:),'none',1);
axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;

subplot(2,3,2);hold on;
title('xeps');
gpatch(FT(indBoundary,:),VT,xeps_FT(indBoundary,:),'none',1);
axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;

subplot(2,3,3);hold on;
title('Ee');
gpatch(FT(indBoundary,:),VT,Ee_FT(indBoundary,:),'none',1);
axisGeom;
colormap(gjet(250)); colorbar;
caxis([min(Ee_data) max(Ee_data)]);
camlight headlight;

subplot(2,3,4);hold on;
title('Vcol');
gpatch(FT(indBoundary,:),VT,Vcol_FT(indBoundary,:),'none',1);
axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;
%axis off

subplot(2,3,5);hold on;
title('Vela');
gpatch(FT(indBoundary,:),VT,Vela_FT(indBoundary,:),'none',1);
axisGeom;
colormap(gjet(250)); colorbar;
camlight headlight;
%axis off
drawnow;

%End of Model Generation steps

%% Create ABAQUS structure
% Setup structure to define an Abaqus inp file

%%--> Heading
abaqus_spec.Heading.COMMENT{1}='Job name: AORTA';

%%--> Preprint
abaqus_spec.Preprint.ATTR.echo='NO';
abaqus_spec.Preprint.ATTR.model='NO';
abaqus_spec.Preprint.ATTR.history='NO';
abaqus_spec.Preprint.ATTR.contact='NO';

%--> Part

% Node
nodeIds=(1:1:size(VT,1))';
abaqus_spec.Part.COMMENT='This section defines the part geometry in terms of nodes and elements';
abaqus_spec.Part.ATTR.name='Aorta';
abaqus_spec.Part.Node={nodeIds,VT};

% Element
elementIds=(1:1:size(ET,1))';
abaqus_spec.Part.Element{1}.ATTR.type='C3D8';
abaqus_spec.Part.Element{1}.VAL={elementIds,ET};

% Element sets
for q=1:1:numMaterials
    elementIdsSetNow=find(C_ET_path_mat_index==q);
    abaqus_spec.Part.Elset{q}.ATTR.elset=['MatSet-',num2str(q)];
    abaqus_spec.Part.Elset{q}.VAL=elementIdsSetNow';
end

surfaceElementSetName='elementSetInnerSurface';
abaqus_spec.Part.Elset{numMaterials+1}.ATTR.elset=surfaceElementSetName;
abaqus_spec.Part.Elset{numMaterials+1}.ATTR.internal=''; %Remains hidden uppon import
abaqus_spec.Part.Elset{numMaterials+1}.VAL=indicesElementsInner(:)';

% Surfaces
sidePick=1;
abaqus_spec.Part.Surface{1}.ATTR.type='ELEMENT';
abaqus_spec.Part.Surface{1}.ATTR.name=[surfaceElementSetName,'_side',num2str(sidePick)];
abaqus_spec.Part.Surface{1}.VAL={surfaceElementSetName,['S',num2str(sidePick)]};

% Sections
for q=1:1:numMaterials
    elementIdsSetNow=find(C_ET_path_mat_index==q);
    abaqus_spec.Part.Solid_section{q}.ATTR.elset=['MatSet-',num2str(q)];
    abaqus_spec.Part.Solid_section{q}.ATTR.material=['Mat_',num2str(q)];
end

%--> Assembly
abaqus_spec.Assembly.ATTR.name='Assembly-1';
abaqus_spec.Assembly.Instance.ATTR.name='Aorta-assembly';
abaqus_spec.Assembly.Instance.ATTR.part='Aorta';

abaqus_spec.Assembly.Nset{1}.ATTR.nset='NSet_Inner';
abaqus_spec.Assembly.Nset{1}.ATTR.instance=abaqus_spec.Assembly.Instance.ATTR.name;
abaqus_spec.Assembly.Nset{1}.VAL=indicesNodesInner;

%Add segment curve node sets
for q=1:1:numel(segmentCurve_cell)
    indNow=numel(abaqus_spec.Assembly.Nset)+1;
    abaqus_spec.Assembly.Nset{indNow}.ATTR.nset=['Segment',num2str(q)];
    abaqus_spec.Assembly.Nset{indNow}.ATTR.instance=abaqus_spec.Assembly.Instance.ATTR.name;
    abaqus_spec.Assembly.Nset{indNow}.VAL=segmentCurve_cell{q};
end

%Add ring curve node sets
for q=1:1:numel(ringCurve_cell)
    indNow=numel(abaqus_spec.Assembly.Nset)+1;
    abaqus_spec.Assembly.Nset{indNow}.ATTR.nset=['Ring',num2str(q)];
    abaqus_spec.Assembly.Nset{indNow}.ATTR.instance=abaqus_spec.Assembly.Instance.ATTR.name;
    abaqus_spec.Assembly.Nset{indNow}.VAL=ringCurve_cell{q};
end

%Add fix node set
% indNow=numel(abaqus_spec.Assembly.Nset)+1;
% setNameFix=['Set-',num2str(indNow)];
% abaqus_spec.Assembly.Nset{indNow}.ATTR.nset=setNameFix;
% abaqus_spec.Assembly.Nset{indNow}.ATTR.instance=abaqus_spec.Assembly.Instance.ATTR.name;
% abaqus_spec.Assembly.Nset{indNow}.VAL=indNodesFix';

%%--> Material
for q=1:1:numMaterials
    abaqus_spec.Material{q}.ATTR.name=['mat_',num2str(q)];
    abaqus_spec.Material{q}.Depvar.VAL=13;
    abaqus_spec.Material{q}.User_Material.ATTR.constants=17;
    %Define material parameters
    mu_now=mu_data;
    k_now=mu_now*kfactor;
    D_now=2/k_now;
    k1_now=k1_data;
    k2_now=k2_data;
    kappa_now=kappa_data;
    theta1_now=theta1_data;
    theta2_now=theta2_data;
    sigact_now=sigact_data;
    ke_now=ke_data;
    Ee_now=Ee_vector(q);
    thetaE1_now=thetaE1_data;
    thetaE2_now=thetaE2_data;
    iswitch_now=iswitch;
    xeps_now=xeps_vector(q);
    Vcol_now=Vcol_vector(q);
    Vela_now=Vela_vector(q);
    Vsmc_now=Vsmc_vector(q);
    t=vec2strIntDouble([mu_now/2 D_now k1_now k2_now kappa_now theta1_now theta2_now sigact_now ke_now Ee_now thetaE1_now thetaE2_now iswitch_now xeps_now Vcol_now Vela_now Vsmc_now],'%6.7e');
    t=strwrap(t,8,', '); %Wrap to max width of 8 entries
    %abaqus_spec.User_Material{q}.VAL=t;
    abaqus_spec.Material{q}.User_Material.VAL=t;
end
%%--> Step
abaqus_spec.Step.ATTR.name='Step-1';
abaqus_spec.Step.ATTR.nlgeom='YES';
abaqus_spec.Step.Static=[0.01 1 1e-6 0.01];

% Boundary
% setNameFix=abaqus_spec.Assembly.Nset{indNow}.ATTR.nset;
% abaqus_spec.Step.Boundary{1}.VAL={setNameFix,[1,1]};
% abaqus_spec.Step.Boundary{2}.VAL={setNameFix,[2,2]};
% abaqus_spec.Step.Boundary{3}.VAL={setNameFix,[3,3]};

%Output
% abaqus_spec.Step.Restart.ATTR.write='';
% abaqus_spec.Step.Restart.ATTR.frequency=0;
%
% abaqus_spec.Step.Output{1}.ATTR.field='';
% abaqus_spec.Step.Output{1}.ATTR.variable='PRESELECT';
% abaqus_spec.Step.Output{2}.ATTR.history='';
% abaqus_spec.Step.Output{2}.ATTR.variable='PRESELECT';
% abaqus_spec.Step.Node_print.ATTR.nset='all';
% abaqus_spec.Step.Node_print.ATTR.frequency = 1;
% abaqus_spec.Step.Node_print.VAL='COORD';
% abaqus_spec.Step.El_print{1}.VAL='S';
% abaqus_spec.Step.El_print{2}.VAL='E';

% Creating the INP file
% You can use |abaqusStruct2inp| to write the structure data to a file.

%%
if saveOn==1
    % Export INP file
    abaqusStruct2inp(abaqus_spec,abaqusInpFileName);

    saveStruct.ET=ET;
    saveStruct.FT=FT;
    saveStruct.VT=VT;
    saveStruct.Fb=Fb;
    saveStruct.segmentCurve_cell=segmentCurve_cell;
    saveStruct.abaqus_spec=abaqus_spec;
    save(matfileSaveName,'-struct','saveStruct');

    disp('inp file write complete')
end

%% FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Fs,Vs,Cs,indEnd,logicRemoveFaces,segmentCurve_cell,E_rings]=circleCutExtrude(F,V,C,V_cent,V_cut,pointSpacing,plotOn,smoothControlParameters,segmentCurve_cell,E_rings)

%% Conver to triangles
Ft=[F(:,[1 2 3]);F(:,[3 4 1])];

%% Resample center line
np=round(max(pathLength(V_cent))/pointSpacing);
V_cent = evenlySampleCurve(V_cent,np,'spline',0);

%% Get mean of cut curve, project to surface

V_cut_mean=mean(V_cut,1);

rMax=max(sqrt(sum((V_cut-V_cut_mean(ones(size(V_cut,1),1),:)).^2,2)));

[~,indMin]=minDist(V_cut_mean,V_cent);

v=V_cent(indMin,:)-V_cut_mean;

optStruct.eps      = 1e-6;
optStruct.triangle = 'one sided';
optStruct.ray      = 'ray';
optStruct.border   = 'normal';

[V_intersect,L_intersect,~] = triangleRayIntersection(V_cut_mean(ones(size(Ft,1),1),:),v(ones(size(Ft,1),1),:),V,Ft,optStruct);

V_intersect=V_intersect(L_intersect,:);

D=minDist(V,V_intersect);

logicRemoveVertices=D<=rMax;

logicRemoveFaces=any(logicRemoveVertices(F),2);

Eb_cut=patchBoundary(F(logicRemoveFaces,:),V);
indCurveCut=edgeListToCurve(Eb_cut);
indCurveCut=indCurveCut(1:end-1);
indCurveCut=indCurveCut(:);

logicFlip=~any((Eb_cut(:,1)==indCurveCut(1))&(Eb_cut(:,2)==indCurveCut(2)));
if logicFlip
    indCurveCut=flipud(indCurveCut);
end

V_curve_cut=V(indCurveCut,:);

%Resample so they have the same amount of points
np=size(V_curve_cut,1);
V_cut = evenlySampleCurve(V_cut,np,'spline',1);

%Get appropriate start point
[~,indMin]=minDist(V_cut(1,:),V_curve_cut);
if indMin>1
    V_curve_cut=[V_curve_cut(indMin:end,:); V_curve_cut(1:indMin-1,:)];
end

V_curve_cut_mean=mean(V_curve_cut,1);
V_cut_mean=mean(V_cut,1);

% Check order again
nq=round(0.25*np);
a1=dot(V_cut(nq,:)-V_cut_mean,V_curve_cut(nq,:)-V_curve_cut_mean);

V_cut_test=flipud(V_cut);
a2=dot(V_cut_test(nq,:)-V_cut_mean,V_curve_cut(nq,:)-V_curve_cut_mean);

if a2>a1
    V_cut=V_cut_test;
end

%Get appropriate start point
[~,indMin]=minDist(V_cut(1,:),V_curve_cut);
if indMin>1
    V_curve_cut=[V_curve_cut(indMin:end,:); V_curve_cut(1:indMin-1,:)];
end

%Minimize twist
SSQD=zeros(np,1);
for q=1:1:np
    if q>1
        indSort=[q:np 1:q-1];
    else
        indSort=1:np;
    end
    V_curve_cut_test=V_curve_cut(indSort,:);
    SSQD(q)=sum(sqrt(sum((V_curve_cut_test-V_cut).^2,2)).^2);
end

%Get sort order
[~,q]=min(SSQD);
if q>1
    indSort=[q:np 1:q-1];
    V_curve_cut=V_curve_cut(indSort,:);
end

%%

F=F(~logicRemoveFaces,:);
C=C(~logicRemoveFaces);

%% Loft

cPar.closeLoopOpt=1;
cPar.patchType='quad';

[F_merge,V_merge,~,indEnd]=polyLoftLinear(V_curve_cut,V_cut,cPar);
C_merge=max(C(:))+ones(size(F_merge,1),1);

%%

% if plotOn==1
%     cFigure(figStruct); hold on;
%     gpatch(F,V,C,'k',0.5);
%     gpatch(F_merge,V_merge,'g','k',1);
%     plotV(V_cent,'g.-','LineWidth',3,'markerSize',25);
%     plotV(V_cut,'r.-','LineWidth',3,'markerSize',25);
%     plotV(V_cut_mean,'r.','markerSize',50);
%     plotV(V_curve_cut,'b.-','LineWidth',3);
%     plotV(V_intersect,'b.','markerSize',50);
%     quiverVec(V_cut_mean,v,[],'r');
%     axisGeom;
%     colormap gjet;
%     camlight headlight;
%     drawnow;
% end

%%

[Fs,Vs,Cs]=joinElementSets({F,F_merge},{V,V_merge},{C,C_merge});
indEnd=indEnd+size(V,1);

%% Constrained smoothing

nGrowthSteps=2;
logicSmooth=any(ismember(Fs,indCurveCut),2);
for q=1:1:nGrowthSteps-1
    indTouch=unique(Fs(logicSmooth,:));
    logicSmooth=any(ismember(Fs,indTouch),2);
end

logicSmooth=logicSmooth | Cs==max(Cs(:));

[Fs,Vs,~,indFix]=mergeVertices(Fs,Vs); %Merging points
indEnd=indFix(indEnd);

for q=1:1:numel(segmentCurve_cell)
    indSegment=segmentCurve_cell{q};
    indSegment=indSegment(indSegment>0);
    segmentCurve_cell{q}=indFix(indSegment);
end
E_rings=indFix(E_rings);

[Fs,Vs,indFix]=patchCleanUnused(Fs,Vs); %removing unused at hole
indEnd=indFix(indEnd);
E_rings=indFix(E_rings);
E_rings=E_rings(all(E_rings>0,2),:);

for q=1:1:numel(segmentCurve_cell)
    indSegment=segmentCurve_cell{q};
    indSegment=indSegment(indSegment>0);
    segmentCurve_cell{q}=indFix(indSegment);
end

indRigid=unique(Fs(~logicSmooth,:));
indRigid=unique([indRigid(:);indEnd(:)]);

Eb=patchBoundary(Fs,Vs);
smoothControlParameters.RigidConstraints=unique([indRigid(:);Eb(:)]);
[Vs]=patchSmooth(Fs,Vs,[],smoothControlParameters);

%%

if plotOn==1
    cFigure(figStruct); hold on;
    gpatch(Fs,Vs,logicSmooth,'k',1);
    patchNormPlot(Fs,Vs);
    plotV(Vs(indRigid,:),'b.','markerSize',25);
    axisGeom;
    colormap gjet;
    camlight headlight;
    drawnow;
end

end

%%
function V=resampleCurve(V,pointSpacing,closeLoopOpt)

np=ceil(max(pathLength(V))/pointSpacing);
V = evenlySampleCurve(V,np,'spline',closeLoopOpt);

end

%%

function [Vn]=curveOffset(V,wallThickness)
p1=mean(V,1); %Curve center

vf=-vecnormalize(V-[V(2:end,:);V(1,:)]); %Allong curve path vectors
vb=vecnormalize(V-[V(end,:);V(1:end-1,:)]); %Allong curve path vectors
v=(vf+vb)/2;

r=vecnormalize(V-p1); %Position vector wrt mean
v1=vecnormalize(cross(v,r)); %perimeter quasi Z-vectors
n=vecnormalize(cross(v1(ones(size(v,1),1),:),v)); %Outward normal vectors
Vn=(V+n*wallThickness); %Offset to create new curve

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

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function [out] = evap_1(S,Ep,dt)
%evap_1

% Copyright (C) 2019, 2021 Wouter J.M. Knoben, Luca Trotter
% This file is part of the Modular Assessment of Rainfall-Runoff Models
% Toolbox (MARRMoT).
% MARRMoT is a free software (GNU GPL v3) and distributed WITHOUT ANY
% WARRANTY. See <https://www.gnu.org/licenses/> for details.

% Flux function
% ------------------
% Description:  Evaporation at the potential rate
% Constraints:  f <= S/dt
% @(Inputs):    S    - current storage [mm]
%               Ep   - potential evaporation rate [mm/d]
%               dt   - time step size

out = min(S/dt,Ep);

end

author: wknoben

repo: MARRMoT

save_path: github-repos/MATLAB/wknoben-MARRMoT


text:
% T in K
% es in Pa
function [T, T_approx] = calculate_saturation_vapor_pressure_liquid_inv (...
es, method, tol_T)
    if (nargin < 2) || isempty(method),  method = '';  end % [] ~= ''
    if (nargin < 3) || isempty(tol_T),  tol_T = eps;  end

    switch method
    case {'Bolton1980', 'approx'}
        % this also serves as a good starting point for more elaborate methods.
        es2 = es ./ 1e3;
        temp = log(es2 ./ 0.6112) ./ 17.67;
        t = 243.5 .* temp ./ (1 - temp);
        T = t + 273.15;
        if (nargout > 1),  T_approx = T;  end
    otherwise
        T_approx = calculate_saturation_vapor_pressure_liquid_inv (es,'approx');
        f = @(T_) calculate_saturation_vapor_pressure_liquid (T_, method);
        temp = warning('off', ...
            'calculate_saturation_vapor_pressure_liquid:outRange');
        T = inv_func2 (f, es, T_approx, tol_T);
        warning(temp);
        % we disable the warning because fzero may rightesously probe
        % outside the valid range, in the search for the root.
    end
end

%!test
%! method = 'Bolton1980';
%! T = (123:10:332)';
%! es = calculate_saturation_vapor_pressure_liquid (T, method);
%! [T2, T3] = calculate_saturation_vapor_pressure_liquid_inv (es, method);
%! %[T, T2, T3, T2-T, T3-T]  % DEBUG
%! %max(abs(T2-T))  % DEBUG
%! tol = eps;
%! myassert(T2, T, -tol)

%!test
%! method = 'Murphy&Koop2005';
%! T = (123:10:332)';
%! es = calculate_saturation_vapor_pressure_liquid (T, method);
%! [T2, T3] = calculate_saturation_vapor_pressure_liquid_inv (es, method);
%! %[T, T2, T3, T2-T, T3-T]  % DEBUG
%! %max(abs(T2-T))  % DEBUG
%! tol = nthroot(eps, 3);
%! myassert(T2, T, -tol)

%!test
%! % empty method should work.
%! method = [];
%! T = (123:10:332)';
%! es = calculate_saturation_vapor_pressure_liquid (T, method);
%! [T2, T3] = calculate_saturation_vapor_pressure_liquid_inv (es, method);

author: Sable

repo: mcbench-benchmarks


text:
function x = p04_root ( i )

%*****************************************************************************80
%
%% P04_ROOT returns a root for problem 4.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 May 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer I, the index of the requested root.
%
%    Output, real X, the value of the root.
%
  if ( i == 1 )
    x = 0.09534461720025875;
  else
    fprintf ( 1, '\n' );
    fprintf ( 1, 'P04_ROOT - Fatal error!\n' );
    fprintf ( 1, '  Illegal root index = %d\n', i );
    error ( 'P04_ROOT - Fatal error!' );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
clear all
        % Declararea variabilei simbolice k
syms k
        % Calcularea primei sume
S1=symsum(k,1,k)
        % Afisarea rezultatului nesimplificat
pretty(S1)
        % Simplificarea rezultatului
S1=simple(S1)
        % Afisarea rezultatului simplificat
pretty(S1)
        % Calcularea sumei a doua
S2=symsum(1/(k*(k+1)),1,inf)


text:
function [Y,FS,NBITS,OPTS] = mp3read(FILE,N,MONO,DOWNSAMP,DELAY)
% MP3READ   Read MP3 audio file via use of external binaries.
%   Y = MP3READ(FILE) reads an mp3-encoded audio file into the
%     vector Y just like wavread reads a wav-encoded file (one channel
%     per column).  Extension ".mp3" is added if FILE has none.
%     Also accepts other formats of wavread, such as
%   Y = MP3READ(FILE,N) to read just the first N sample frames (N
%     scalar), or the frames from N(1) to N(2) if N is a two-element vector.
%   Y = MP3READ(FILE,FMT) or Y = mp3read(FILE,N,FMT)
%     with FMT as 'native' returns int16 samples instead of doubles;
%     FMT can be 'double' for default behavior (to exactly mirror the
%     syntax of wavread).
%
%   [Y,FS,NBITS,OPTS] = MP3READ(FILE...) returns extra information:
%     FS is the sampling rate,  NBITS is the bit depth (always 16),
%     OPTS.fmt is a format info string; OPTS has multiple other
%     fields, see WAVREAD.
%
%   SIZ = MP3READ(FILE,'size') returns the size of the audio data contained
%     in the file in place of the actual audio data, returning the
%     2-element vector SIZ=[samples channels].
%
%   [Y...] = MP3READ(FILE,N,MONO,DOWNSAMP,DELAY) extends the
%     WAVREAD syntax to allow access to special features of the
%     mpg123 engine:  MONO = 1 forces output to be mono (by
%     averaging stereo channels); DOWNSAMP = 2 or 4 downsamples by
%     a factor of 2 or 4 (thus FS returns as 22050 or 11025
%     respectively for a 44 kHz mp3 file); DELAY controls how many
%     "warm up" samples to drop at the start of the file; the
%     default value of 2257 makes an mp3write/mp3read loop for a 44
%     kHz mp3 file be as close as possible to being temporally
%     aligned; specify as 0 to prevent discard of initial samples.
%
%   Example:
%   To read an mp3 file as doubles at its original width and sampling rate:
%     [Y,FS] = mp3read('piano.mp3');
%   To read the first 1 second of the same file, downsampled by a
%   factor of 4, cast to mono, using the default filename
%   extension:
%     [Y,FS4] = mp3read('piano', FS/4, 1, 4);
%
%   Note: Because the mp3 format encodes samples in blocks of 26 ms (at
%   44 kHz), and because of the "warm up" period of the encoder,
%   the file length may not be exactly what you expect.
%
%   Note: requires external binaries mpg123 and mp3info; you
%   can find binaries for several platforms at:
%     http://labrosa.ee.columbia.edu/matlab/mp3read.html
%
%   See also mp3write, wavread.

% 2003-07-20 dpwe@ee.columbia.edu  This version calls mpg123.
% 2004-08-31 Fixed to read whole files correctly
% 2004-09-08 Uses mp3info to get info about mp3 files too
% 2004-09-18 Reports all mp3info fields in OPTS.fmt; handles MPG2LSF sizes
%            + added MONO, DOWNSAMP flags, changed default behavior.
% 2005-09-28 Fixed bug reading full-rate stereo as 1ch (thx bjoerns@vjk.dk)
% 2006-09-17 Chop off initial 2257 sample delay (for 44.1 kHz mp3)
%            so read-write loop doesn't get progressively delayed.
%            You can suppress this with a 5th argument of 0.
% 2007-02-04 Added support for FMT argument to match wavread
%            Added automatic selection of binary etc. to allow it
%            to work cross-platform without editing prior to
%            submitting to Matlab File Exchange
% 2007-07-23 Tweaks to 'size' mode so it exactly agrees with read data.

% find our baseline directory
path = fileparts(which('mp3read'));

% %%%%% Directory for temporary file (if needed)
% % Try to read from environment, or use /tmp if it exists, or use CWD
tmpdir = getenv('TMPDIR');
if isempty(tmpdir) || exist(tmpdir,'file')==0
  tmpdir = '/tmp';
end
if exist(tmpdir,'file')==0
  tmpdir = '';
end
% ensure it exists
%if length(tmpdir) > 0 && exist(tmpdir,'file')==0
%  mkdir(tmpdir);
%end

%%%%%% Command to delete temporary file (if needed)
rmcmd = 'rm';

%%%%%% Location of the binaries - attempt to choose automatically
%%%%%% (or edit to be hard-coded for your installation)
ext = lower(computer);
if ispc
  ext = 'exe';
  rmcmd = 'del';
end
mpg123 = fullfile(path,['mpg123.',ext]);
mp3info = fullfile(path,['mp3info.',ext]);

%%%%% Process input arguments
if nargin < 2
  N = 0;
end

% Check for FMT spec (per wavread)
FMT = 'double';
if ischar(N)
  FMT = lower(N);
  N = 0;
end

if length(N) == 1
  % Specified N was upper limit
  N = [1 N];
end
if nargin < 3
  forcemono = 0;
else
  % Check for 3rd arg as FMT
  if ischar(MONO)
    FMT = lower(MONO);
    MONO = 0;
  end
  forcemono = (MONO ~= 0);
end
if nargin < 4
  downsamp = 1;
else
  downsamp = DOWNSAMP;
end
if downsamp ~= 1 && downsamp ~= 2 && downsamp ~= 4
  error('DOWNSAMP can only be 1, 2, or 4');
end
if nargin < 5
  mpg123delay44kHz = 2257;  % empirical delay of lame/mpg123 loop
  delay = round(mpg123delay44kHz/downsamp);
else
  delay = DELAY;
end

if strcmp(FMT,'native') == 0 && strcmp(FMT,'double') == 0 && ...
      strcmp(FMT,'size') == 0
  error(['FMT must be ''native'' or ''double'' (or ''size''), not ''',FMT,'''']);
end

%%%%%% Constants
NBITS=16;

%%%%% add extension if none (like wavread)
[path,file,ext] = fileparts(FILE);
if isempty(ext)
  FILE = [FILE, '.mp3'];
end

%%%%%% Probe file to find format, size, etc. using "mp3info" utility
cmd = ['"',mp3info, '" -r m -p "%Q %u %b %r %v * %C %e %E %L %O %o %p" "', FILE,'"'];
% Q = samprate, u = #frames, b = #badframes (needed to get right answer from %u)
% r = bitrate, v = mpeg version (1/2/2.5)
% C = Copyright, e = emph, E = CRC, L = layer, O = orig, o = mono, p = pad
w = mysystem(cmd);
% Break into numerical and ascii parts by finding the delimiter we put in
starpos = findstr(w,'*');
nums = str2num(w(1:(starpos - 2)));
strs = tokenize(w((starpos+2):end));

SR = nums(1);
nframes = nums(2);
nchans = 2 - strcmp(strs{6}, 'mono');
layer = length(strs{4});
bitrate = nums(4)*1000;
mpgv = nums(5);
% Figure samples per frame, after
% http://board.mp3-tech.org/view.php3?bn=agora_mp3techorg&key=1019510889
if layer == 1
  smpspfrm = 384;
elseif SR < 32000 && layer ==3
  smpspfrm = 576;
  if mpgv == 1
    error('SR < 32000 but mpeg version = 1');
  end
else
  smpspfrm = 1152;
end

OPTS.fmt.mpgBitrate = bitrate;
OPTS.fmt.mpgVersion = mpgv;
% fields from wavread's OPTS
OPTS.fmt.nAvgBytesPerSec = bitrate/8;
OPTS.fmt.nSamplesPerSec = SR;
OPTS.fmt.nChannels = nchans;
OPTS.fmt.nBlockAlign = smpspfrm/SR*bitrate/8;
OPTS.fmt.nBitsPerSample = NBITS;
OPTS.fmt.mpgNFrames = nframes;
OPTS.fmt.mpgCopyright = strs{1};
OPTS.fmt.mpgEmphasis = strs{2};
OPTS.fmt.mpgCRC = strs{3};
OPTS.fmt.mpgLayer = strs{4};
OPTS.fmt.mpgOriginal = strs{5};
OPTS.fmt.mpgChanmode = strs{6};
OPTS.fmt.mpgPad = strs{7};
OPTS.fmt.mpgSampsPerFrame = smpspfrm;

if SR == 16000 && downsamp == 4
  error('mpg123 will not downsample 16 kHz files by 4 (only 2)');
end

if downsamp == 1
  downsampstr = '';
else
  downsampstr = [' -',num2str(downsamp)];
end
FS = SR/downsamp;

if forcemono == 1
  nchans = 1;
  chansstr = ' -m';
else
  chansstr = '';
end

% Size-reading version
if strcmp(FMT,'size') == 1
   Y = [floor(smpspfrm*nframes/downsamp)-delay, nchans];
else

  % Temporary file to use
  tmpfile = fullfile(tmpdir, ['tmp',num2str(round(1000*rand(1))),'.wav']);

  skipx = 0;
  skipblks = 0;
  skipstr = '';
  sttfrm = N(1)-1;

  % chop off transcoding delay?
  %sttfrm = sttfrm + delay;  % empirically measured
  % no, we want to *decode* those samples, then drop them
  % so delay gets added to skipx instead

  if sttfrm > 0
    skipblks = floor(sttfrm*downsamp/smpspfrm);
    skipx = sttfrm - (skipblks*smpspfrm/downsamp);
    skipstr = [' -k ', num2str(skipblks)];
  end
  skipx = skipx + delay;

  lenstr = '';
  endfrm = -1;
  decblk = 0;
  if length(N) > 1
    endfrm = N(2);
    if endfrm > sttfrm
      decblk = ceil((endfrm+delay)*downsamp/smpspfrm) - skipblks + 10;
      % we read 10 extra blks (+10) to cover the case where up to 10 bad
      % blocks are included in the part we are trying to read (it happened)
      lenstr = [' -n ', num2str(decblk)];
      % This generates a spurious "Warn: requested..." if reading right
      % to the last sample by index (or bad blks), but no matter.
    end
 end

  % Run the decode
  cmd=['"',mpg123,'"', downsampstr, chansstr, skipstr, lenstr, ...
       ' -q -w "', tmpfile,'"  "',FILE,'"'];
  %w =
  mysystem(cmd);

  % Load the data
  Y = wavread(tmpfile);

%  % pad delay on to end, just in case
%  Y = [Y; zeros(delay,size(Y,2))];
%  % no, the saved file is just longer

  if decblk > 0 && length(Y) < decblk*smpspfrm/downsamp
    % This will happen if the selected block range includes >1 bad block
    disp(['Warn: requested ', num2str(decblk*smpspfrm/downsamp),' frames, returned ',num2str(length(Y))]);
  end

  % Delete tmp file
  mysystem([rmcmd,' "', tmpfile,'"']);

  % debug
%  disp(['sttfrm=',num2str(sttfrm),' endfrm=',num2str(endfrm),' skipx=',num2str(skipx),' delay=',num2str(delay),' len=',num2str(length(Y))]);

  % Select the desired part
  if skipx+endfrm-sttfrm > length(Y)
      endfrm = length(Y)+sttfrm-skipx;
  end

  if endfrm > sttfrm
    Y = Y(skipx+(1:(endfrm-sttfrm)),:);
  elseif skipx > 0
    Y = Y((skipx+1):end,:);
  end

  % Convert to int if format = 'native'
  if strcmp(FMT,'native')
    Y = int16((2^15)*Y);
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function w = mysystem(cmd)
% Run system command; report error; strip all but last line
[s,w] = system(cmd);
if s ~= 0
  error(['unable to execute ',cmd,' (',w,')']);
end
% Keep just final line
w = w((1+max([0,findstr(w,10)])):end);
% Debug
%disp([cmd,' -> ','*',w,'*']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = tokenize(s)
% Break space-separated string into cell array of strings
% 2004-09-18 dpwe@ee.columbia.edu
a = [];
p = 1;
n = 1;
l = length(s);
nss = findstr([s(p:end),' '],' ');
for ns = nss
  % Skip initial spaces
  if ns == p
    p = p+1;
  else
    if p <= l
      a{n} = s(p:(ns-1));
      n = n+1;
      p = ns+1;
    end
  end
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
%KGAUSS Gaussian kernel
%
% K = KGAUSS(SIGMA) is a 2-dimensional Gaussian kernel of standard deviation
% SIGMA, and  centred within the matrix K whose half-width is H=2xSIGMA and
% W=2xH+1.
%
% K = KGAUSS(SIGMA, H) as above but the half-width H is specified.
%
% Notes::
% - The volume under the Gaussian kernel is one.
%
% See also KDGAUSS, KDOG, KLOG, ICONV.

% Copyright (C) 1993-2011, by Peter I. Corke
%
% This file is part of The Machine Vision Toolbox for Matlab (MVTB).
%
% MVTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MVTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with MVTB.  If not, see <http://www.gnu.org/licenses/>.

function m = kgauss(sigma, w)

    if nargin == 1,
        w = ceil(3*sigma);
    end
    ww = 2*w + 1;

    [x,y] = meshgrid(-w:w, -w:w);

    m = 1/(2*pi*sigma^2) * exp( -(x.^2 + y.^2)/2/sigma^2);

    % area under the curve should be 1, but the discrete case is only
    % an approximation, correct it
    %m = m / sum(m(:));

author: petercorke

repo: machinevision-toolbox-matlab


text:
function t_load2disp(quiet)
%T_LOAD2DISP  Tests for LOAD2DISP.

%   MATPOWER
%   Copyright (c) 2016, Power Systems Engineering Research Center (PSERC)
%   by Ray Zimmerman, PSERC Cornell
%
%   This file is part of MATPOWER.
%   Covered by the 3-clause BSD License (see LICENSE file for details).
%   See https://matpower.org for more info.

if nargin < 1
    quiet = 0;
end

num_tests = 72;

t_begin(num_tests, quiet);

[PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
    VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
    MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
    TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
    ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

casefile = 't_case9_opf';
if quiet
    verbose = 0;
else
    verbose = 0;
end
if have_feature('octave')
    if have_feature('octave', 'vnum') >= 4
        file_in_path_warn_id = 'Octave:data-file-in-path';
    else
        file_in_path_warn_id = 'Octave:load-file-in-path';
    end
    s1 = warning('query', file_in_path_warn_id);
    warning('off', file_in_path_warn_id);
end

mpopt = mpoption('opf.violation', 1e-6, 'mips.gradtol', 1e-8, ...
        'mips.comptol', 1e-8, 'mips.costtol', 1e-9);
mpopt = mpoption(mpopt, 'out.all', 0, 'verbose', verbose, 'opf.ac.solver', 'MIPS');

%% set up indices
ib_data     = [1:BUS_AREA BASE_KV:VMIN];
ib_voltage  = [VM VA];
ib_lam      = [LAM_P LAM_Q];
ib_mu       = [MU_VMAX MU_VMIN];
ig_data     = [GEN_BUS QMAX QMIN MBASE:APF];
ig_disp     = [PG QG VG];
ig_mu       = (MU_PMAX:MU_QMIN);
ibr_data    = (1:ANGMAX);
ibr_flow    = (PF:QT);
ibr_mu      = [MU_SF MU_ST];
ibr_angmu   = [MU_ANGMIN MU_ANGMAX];

%% get solved AC OPF case from MAT-file
load soln9_opf;     %% defines bus_soln, gen_soln, branch_soln, f_soln

%% load2disp
mpc0 = loadcase(casefile);
mpc1 = load2disp(mpc0);
ng = size(mpc0.gen, 1);
ib = find(mpc0.bus(:, PD));
ig = ng+(1:length(ib));

t = 'mpc = load2disp(mpc) : ';
t_is(size(mpc1.gen, 1), size(mpc0.gen, 1)+length(ib), 12, [t 'number of gens']);
t_is(mpc1.bus(:, PD), 0, 12, [t 'PD is zero']);
t_is(mpc1.bus(:, QD), 0, 12, [t 'QD is zero']);
t_is(mpc1.gen(ig, PG), -mpc0.bus(ib, PD), 12, [t 'PG = -PD']);
t_is(mpc1.gen(ig, PMIN), -mpc0.bus(ib, PD), 12, [t 'PMIN = -PD']);
t_is(mpc1.gen(ig, PMAX), 0, 12, [t 'PMAX = 0']);
t_is(mpc1.gen(ig, QMIN), -mpc0.bus(ib, QD), 12, [t 'QMIN = -QD']);
t_is(mpc1.gen(ig, QMAX), 0, 12, [t 'QMAX = 0']);
t_is(mpc1.gen(ig, VG), mpc0.bus(ib, VM), 12, [t 'VG']);
t_is(mpc1.gencost(ig, MODEL), 2, 12, [t 'gencost(:, MODEL)']);
t_is(mpc1.gencost(ig, NCOST), 2, 12, [t 'gencost(:, NCOST)']);
t_is(mpc1.gencost(ig, COST), 5000, 12, [t 'gencost(:, COST)']);

t = 'mpc = load2disp(mpc, '''', idx, voll) : ';
idx = [5;7];
ig = ng+(1:length(idx));
mpc2 = load2disp(mpc0, '', idx, 1000);
t_is(size(mpc2.gen, 1), size(mpc0.gen, 1)+length(idx), 12, [t 'number of gens']);
t_is(mpc2.bus(:, PD), [zeros(8,1); 125], 12, [t 'PD']);
t_is(mpc2.bus(:, QD), [zeros(8,1);  50], 12, [t 'QD']);
t_is(mpc2.gen(ig, PG), -mpc0.bus(idx, PD), 12, [t 'PG = -PD']);
t_is(mpc2.gen(ig, PMIN), -mpc0.bus(idx, PD), 12, [t 'PMIN = -PD']);
t_is(mpc2.gen(ig, PMAX), 0, 12, [t 'PMAX = 0']);
t_is(mpc2.gen(ig, QMIN), -mpc0.bus(idx, QD), 12, [t 'QMIN = -QD']);
t_is(mpc2.gen(ig, QMAX), 0, 12, [t 'QMAX = 0']);
t_is(mpc2.gen(ig, VG), mpc0.bus(idx, VM), 12, [t 'VG']);
t_is(mpc2.gencost(ig, MODEL), 2, 12, [t 'gencost(:, MODEL)']);
t_is(mpc2.gencost(ig, NCOST), 2, 12, [t 'gencost(:, NCOST)']);
t_is(mpc2.gencost(ig, COST), 1000, 12, [t 'gencost(:, COST)']);

%% run OPF
t = 'fixed load OPF soln : ';
r0 = runopf(mpc0, mpopt);
t_ok(r0.success, [t 'success']);
t_is(r0.f, f_soln, 3, [t 'f']);
t_is(r0.bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
t_is(r0.bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
t_is(r0.bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
t_is(r0.bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
t_is(r0.gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
t_is(r0.gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
t_is(r0.gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
t_is(r0.branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
t_is(r0.branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
t_is(r0.branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

t = 'disp load OPF soln : ';
r1 = runopf(mpc1, mpopt);
t_ok(r0.success, [t 'success']);
t_is(r0.f, f_soln, 3, [t 'f']);
t_is(r0.bus(:,ib_data   ),    bus_soln(:,ib_data   ), 10, [t 'bus data']);
t_is(r0.bus(:,ib_voltage),    bus_soln(:,ib_voltage),  3, [t 'bus voltage']);
t_is(r0.bus(:,ib_lam    ),    bus_soln(:,ib_lam    ),  3, [t 'bus lambda']);
t_is(r0.bus(:,ib_mu     ),    bus_soln(:,ib_mu     ),  2, [t 'bus mu']);
t_is(r0.gen(:,ig_data   ),    gen_soln(:,ig_data   ), 10, [t 'gen data']);
t_is(r0.gen(:,ig_disp   ),    gen_soln(:,ig_disp   ),  3, [t 'gen dispatch']);
t_is(r0.gen(:,ig_mu     ),    gen_soln(:,ig_mu     ),  3, [t 'gen mu']);
t_is(r0.branch(:,ibr_data  ), branch_soln(:,ibr_data  ), 10, [t 'branch data']);
t_is(r0.branch(:,ibr_flow  ), branch_soln(:,ibr_flow  ),  3, [t 'branch flow']);
t_is(r0.branch(:,ibr_mu    ), branch_soln(:,ibr_mu    ),  2, [t 'branch mu']);

t = 'mpc = load2disp(solved_mpc) + gentype/genfuel : ';
r0.genfuel = {'ng'; 'coal'; 'hydro'};
r0.gentype = {'GT'; 'ST'; 'HY'};
mpc3 = load2disp(r0);
ib = find(mpc0.bus(:, PD));
ig = ng+(1:length(ib));
t_is(size(mpc3.gen, 1), size(mpc0.gen, 1)+length(ib), 12, [t 'number of gens']);
t_is(mpc3.bus(:, PD), 0, 12, [t 'PD is zero']);
t_is(mpc3.bus(:, QD), 0, 12, [t 'QD is zero']);
t_is(mpc3.gen(ig, PG), -mpc0.bus(ib, PD), 12, [t 'PG = -PD']);
t_is(mpc3.gen(ig, PMIN), -mpc0.bus(ib, PD), 12, [t 'PMIN = -PD']);
t_is(mpc3.gen(ig, PMAX), 0, 12, [t 'PMAX = 0']);
t_is(mpc3.gen(ig, QMIN), -mpc0.bus(ib, QD), 12, [t 'QMIN = -QD']);
t_is(mpc3.gen(ig, QMAX), 0, 12, [t 'QMAX = 0']);
t_is(mpc3.gen(ig, VG), r0.bus(ib, VM), 12, [t 'VG']);
t_is(mpc3.gencost(ig, MODEL), 2, 12, [t 'gencost(:, MODEL)']);
t_is(mpc3.gencost(ig, NCOST), 2, 12, [t 'gencost(:, NCOST)']);
t_is(mpc3.gencost(ig, COST), 5000, 12, [t 'gencost(:, COST)']);
t_ok(isfield(mpc3, 'gentype'), [t 'results.gentype']);
t_ok(isfield(mpc3, 'genfuel'), [t 'results.genfuel']);
if isfield(mpc3, 'gentype')
    t_ok(isequal(mpc3.gentype, {r0.gentype{:}, 'DL', 'DL', 'DL'}'), [t 'results.results.gentype']);
else
    t_ok(0, [t 'results.results.gentype']);
    got = mpc3.gentype
    expected = {r0.gentype{:}, 'DL', 'DL', 'DL'}'
end
if isfield(mpc3, 'genfuel')
    t_ok(isequal(mpc3.genfuel, {r0.genfuel{:}, 'dl', 'dl', 'dl'}'), [t 'results.results.genfuel']);
else
    t_ok(0, [t 'results.results.genfuel']);
    got = mpc3.genfuel
    expected = {r0.genfuel{:}, 'dl', 'dl', 'dl'}'
end

t = 'loadshed(mpc) : ';
shed = loadshed(mpc3.gen);
t_is(size(shed), [3 1], 12, [t 'size']);
t_is(shed, [0;0;0], 12, [t 'all zeros']);

t = 'loadshed(mpc, ild) : ';
ild = [6; 5];
shed = loadshed(mpc3.gen, ild);
t_is(size(shed), [2 1], 12, [t 'size']);
t_is(shed, [0;0], 12, [t 'all zeros']);

ild = find(isload(mpc3.gen));
mpc3.gen(ild, PG) = [-80; -100; -120];

t = 'loadshed(mpc) : ';
shed = loadshed(mpc3.gen);
t_is(size(shed), [3 1], 12, [t 'size']);
t_is(shed, [10;0;5], 12, [t 'all zeros']);

t = 'loadshed(mpc, ild) : ';
ild = [6; 5];
shed = loadshed(mpc3.gen, ild);
t_is(size(shed), [2 1], 12, [t 'size']);
t_is(shed, [5;0], 12, [t 'all zeros']);

if have_feature('octave')
    warning(s1.state, file_in_path_warn_id);
end

t_end;

author: MATPOWER

repo: matpower

save_path: github-repos/MATLAB/MATPOWER-matpower


text:
function sig_out = wsola_time_scaling(sig_in, fs, scale_factor, simil_method)
% sig_in - input signal
% fs - sampling frequency
% scale_factor - by what ratio to strech/compress the audio.
% simil_method - optional parameter either 'xcorr' or 'amdf' (default)
%
% from the paper by Werner VERHELST and Marc ROELANDS
% found here:
% http://www.etro.vub.ac.be/research/dssp/PUB_FILES/int_conf/ICASSP-1993.pdf

if nargin<4
	simil_method = 'amdf';
end

sf_split_margin = 3;
orig_len = length(sig_in);
% recursively treat scale factors too big or too small
if (scale_factor > sf_split_margin) || (scale_factor < 1/sf_split_margin)
    num_recur = ceil(abs((log(scale_factor)./log(sf_split_margin))));

    scale_factor_recur = scale_factor.^((num_recur-1)/num_recur);

    sig_in = wsola_time_scaling(sig_in, fs, scale_factor_recur, simil_method);

    scale_factor = scale_factor./(length(sig_in)./orig_len);
end

% consts
win_time = 0.020; %sec
overlap_ratio = 0.5;
max_err = min(0.005, win_time*overlap_ratio./scale_factor/2);

% lengths
win_len = ceil(win_time.*fs);
max_err_len = ceil(max_err.*fs);
step_len = floor(overlap_ratio.*win_len);

% vectors
win = hann(win_len);
% orig_scale = 1:length(sig_in);
new_scale = (1:round((length(sig_in).*scale_factor*2)))';
sig_out = zeros(size(new_scale));
win_out = zeros(size(new_scale));

cursor_in = 1;
cursor_out = 1;

while cursor_in<(length(sig_in)-win_len-max(step_len, step_len./scale_factor) - 2*max_err_len)...
                                    && cursor_out<(length(sig_out)-win_len)
    % input segments
    new_seg = sig_in(cursor_in:(cursor_in+win_len-1)).*win;
    new_seg_neihbour = sig_in((cursor_in+step_len):(cursor_in+step_len+win_len-1)).*win;
    % overlap add
    sig_out(cursor_out:(cursor_out+win_len-1)) = ...
            sig_out(cursor_out:(cursor_out+win_len-1))+...
            new_seg;
    % overlap add window normalization vec
    win_out(cursor_out:(cursor_out+win_len-1)) = ...
            win_out(cursor_out:(cursor_out+win_len-1))+...
            win;
    % move cursors
    cursor_out = cursor_out + step_len;
    cursor_in = cursor_in + round(step_len./scale_factor);
    % new candidate
    new_seg_cand = sig_in((cursor_in):(cursor_in+win_len-1)).*win;
    % similarity calc. Note: matlab apears to be using FFT for xcorr, and
    % so computes the whole thing (instead of just the center).
    % Writing own version of xcorr would probably make it faster.
    if strcmp(simil_method, 'xcorr')
        shift = max_xcorr_similarity(new_seg_neihbour, new_seg_cand, max_err_len);
    elseif strcmp(simil_method, 'amdf')
        shift = min_amdf_similarity(new_seg_neihbour, new_seg_cand, max_err_len);
    else
        return
    end
    % adjust cursor place
    cursor_in = cursor_in - shift;
end

% remove slack
sig_out(cursor_out:end) = [];
win_out(cursor_out:end) = [];

sig_out = sig_out./(win_out+eps); %normalize to remove possible modulations

end

    % cross correlation based similarity calculation
    function shift = max_xcorr_similarity(seg1, seg2, max_lag)
        [~, max_i] = max(xcorr(seg1, seg2, max_lag, 'unbiased'));
        shift = max_i - max_lag;
    end

    % amdf based similarity calculation
    function shift = min_amdf_similarity(seg1, seg2, max_lag)
        n = length(seg1);
        amdf = ones(1,2*max_lag-1);
        for lag=-max_lag:max_lag
            amdf(lag+max_lag+1) = sum(abs(seg2(max(1,(-lag+1)):min(n,(n-lag)))-...
                           seg1(max(1,(lag+1)):min(n,(n+lag))) ))/n;
        end
        [~, min_i] = min(amdf);
        shift = min_i - max_lag;
    end


text:
function view = roiLoadHdf5(view, select, clr)
% got rid of absPathFlag,defaultPathFlag
%
% view = roiLoadHdf5(view,  [select], [color])

%% Analagous to old getROIfilename.m

global MRFILES;

pathStr=['/', view.subdir, '/ROIs'];
MRFILES = mrFilesSet(MRFILES, 'path', pathStr);

groups = mrFilesGet(MRFILES, 'groups');
% there should be no datasets here
datasets = mrFilesGet(MRFILES, 'datasets');

if ~isempty(datasets)
    error('mrFiles:roiLoadHdf5:Malformed', ...
              ['found dataset under' pathStr]);
end

if isempty(groups)
  myErrorDlg('No ROIs found')
end

roi_offset = length(pathStr)+2;
group_names = {};
for i = 1:length(groups)
    group_names = {group_names{:}, groups{i}(roi_offset:end)};
end

% right now we just return the max / most recent ROI from each group.  we
% should offer a full range!
[sel,ok] = listdlg('PromptString','Select an ROI',...
       'SelectionMode','multiple',...
       'ListString',group_names);

if ok
    groups = groups(sel);
    group_names = group_names(sel);
else
    groups = {};
    group_names = {};
end

%% Analagous to loadROI.m
if ~exist('select', 'var'),            select=1;           end

for i=1:length(groups)
    % This is a little inefficient, but not so bad
    MRFILES = mrFilesSet(MRFILES, 'path', groups{i});

    MRFILES = mrfPos(MRFILES, 'max');

    disp(['loading ', groups{i}]);
    thisROI = mrfLoadHdf5(MRFILES, 'coords');

    thisROI.name = group_names{i};

    % SHOULD THIS BE DONE?
    % this is guaranteed by the fact that we're using view.subdir
    % using view.subdir directly might even be better, but I'm guessing
    % these things are separate for a reason...
    thisROI.viewType = view.viewType;

    % Check color is set well
    if exist('clr','var') % from the function call
        thisROI.ROI.color=clr;
    end
    if ~isfield(thisROI,'color')
        thisROI.color='b';
    end

    view = addROI(view,thisROI,select);
end

disp('Done loading');

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function test_ft_volumebiascorrect

% MEM 8gb
% WALLTIME 00:60:00
% DEPENDENCY ft_volumebiascorrect

mri = [];
mri.anatomy = randn(181,217,181);
mri.dim = [181 217 181];
mri.transform = eye(4);
mri.unit = 'cm';
mri.coordsys = 'ctf';

cfg = [];
mriout = ft_volumebiascorrect(cfg, mri);

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [favar]=favar_fevd(gamma_record,It,Bu,n,IRFperiods,FEVDband,favar,IRFt,strctident)

% function [fevd_estimates]=bear.olsfevd(irf_estimates,IRFperiods,gamma,n,endo,datapath)
% computes and displays fevd values for the OLS VAR model
% inputs:  - cell 'irf_estimates': lower bound, point estimates, and upper bound for the IRFs
%          - integer 'IRFperiods': number of periods for IRFs
%          - matrix 'gamma': structural disturbance variance-covariance matrix (defined p 48 of technical guide)
%          - integer 'n': number of endogenous variables in the BVAR model (defined p 7 of technical guide)
%          - cell 'endo': list of endogenous variables of the model
%          - string 'datapath': user-supplied path to excel data spreadsheet
% outputs: - cell 'fevd_estimates': lower bound, point estimates, and upper bound for the FEVD

% preliminary tasks
npltX=favar.npltX;

% load IRFs
favar_irf_record=favar.IRF.favar_irf_record;

if IRFt==1||IRFt==2||IRFt==3
    identified=n; % fully identified
elseif IRFt==4 || IRFt==6 %if the model is identified by sign restrictions or sign restrictions (+ IV)
    identified=size(strctident.signreslabels,1); % count the labels provided in the sign res sheet (+ IV)
elseif IRFt==5
    identified=1; % one IV shock
end

% create the first cell
temp=cell(npltX,identified+1);

% start by filling the first column of every Tij matrix in the cell
% loop over rows of temp
for jj=1:npltX
    % loop over columns of temp
    for ii=1:identified
        % square each element
        temp{jj,ii}(:,1)=favar_irf_record{jj,ii}(:,1).^2;
    end
end
% fill all the other entries of the Tij matrices
% loop over rows of temp
for jj=1:npltX
    % loop over columns of temp
    for ii=1:identified
        % loop over remaining columns
        for kk=2:IRFperiods
            % define the column as the square of the corresponding column in orthogonalised_irf_record
            % additioned to the value of the preceeding columns, which creates the cumulation
            temp{jj,ii}(:,kk)=favar_irf_record{jj,ii}(:,kk).^2+temp{jj,ii}(:,kk-1);
        end
    end
end

% reshape gamma for loop
gamma_gibbs=reshape(gamma_record,n,n,It-Bu);

% recall L from the sampling process in this case, analogue to beta and sigma
Lgibbs=reshape(favar.L_gibbs,size(favar.L,1),size(favar.L,2),It-Bu);
%relevant loadings of restricted information variables
Lgibbs=Lgibbs(favar.plotX_index,:,:);
% R2 for scaling
R2_gibbs=reshape(favar.R2_gibbs,size(favar.plotX_index,1),1,It-Bu);

% multiply each matrix in the cell by the variance of the structural shocks
% to do so, loop over simulations (rows of the Tij matrices)
for kk=1:It-Bu
    % recover the covariance matrix of structural shocks gamma for this iteration
    gamma=squeeze(gamma_gibbs(:,:,kk));
    L=squeeze(Lgibbs(:,:,kk));
    R2=squeeze(R2_gibbs(:,:,kk));

    % scale gamma, irf estimates are already scaled,do we need this step????
    for ii=1:npltX
        for ll=1:identified
            favar_gamma{ii}(:,ll)=L(ii,ll)*gamma(:,ll);
        end
    end

    % loop over rows of temp
    for ii=1:npltX
        % loop over columns of temp
        for jj=1:identified
            % multiply column jj of the matrix by the variance of the structural shock
            temp{ii,jj}(1,:)=temp{ii,jj}(1,:)*favar_gamma{ii}(jj,jj);
        end
    end
end

% obtain now the values for Ti, the (n+1)th matrix of each row
% loop over rows of temp
for ii=1:npltX
    % start the summation over Tij matrices
    temp{ii,identified+1}=temp{ii,1};
    % sum over remaining columns
    for jj=2:identified
        temp{ii,identified+1}=temp{ii,identified+1}+temp{ii,jj};
    end
end

% create the output cell fevd_record, scale the shocks with R2 in spirit of BBE (2005)
favar_fevd_record=cell(npltX,identified);

% fill the cell
% loop over rows of fevd_estimates
for ii=1:npltX
    % load the R2 to determine the "true" share of variance explained by
    % the common component
    scale=R2(ii);
    shocks=[];
    % loop over columns of fevd_estimates
    for jj=1:identified
        % define the matrix Vfij as the division (pairwise entry) of Tfij by Tfj
        shock=(temp{ii,jj}./temp{ii,identified+1})*scale;
        favar_fevd_record{ii,jj}=shock; % changed to abs(shock) from shock here
        % save shocks to compute residual
        shocks(:,:,jj)=shock;
    end
    % finally add the residual, reflecting the share explained by the idiosyncratic component (residual)
    favar_fevd_record{ii,jj+1}=1-sum(shocks,3);
end

%% create the FEVD estimates output
% create first the cell that will contain the estimates
favar_fevd_estimates=cell(npltX,identified+1);

% for each variable and each variable contribution along with each period, compute the median, lower and upper bound from the Gibbs sampler records
% consider variables in turn
for ii=1:npltX
    % consider contributions in turn
    for jj=1:identified+1
        % consider periods in turn
        for kk=1:IRFperiods
            % compute first the lower bound
            favar_fevd_estimates{ii,jj}(1,kk)=quantile(favar_fevd_record{ii,jj}(:,kk),(1-FEVDband)/2);
            % then compute the median
            favar_fevd_estimates{ii,jj}(2,kk)=quantile(favar_fevd_record{ii,jj}(:,kk),0.5);
            % finally compute the upper bound
            favar_fevd_estimates{ii,jj}(3,kk)=quantile(favar_fevd_record{ii,jj}(:,kk),1-(1-FEVDband)/2);
        end
    end
end

%save output
favar.FEVD.favar_fevd_estimates=favar_fevd_estimates;

author: european-central-bank

repo: BEAR-toolbox


text:
function printH(H)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=size(H,1);
m=size(H,2);
fprintf('H     =');
fprintf('\n');

for i=1:n
    for j=1:m
        if H(i,j)==0,continue;end
        fprintf('%10.5f',H(i,j));
    end
    fprintf('\n');
end

return;

repo: GINav


text:
function [elem2edge,edge,elem2edgeSign] = dof3edge(elem)
%% DOF3EDGE dof structure for the lowest order edge elements.
%
% [elem2edge,edge] = DOF3EDGE(elem) constructs data structure
% for the lowest order edge element. elem is the connectivity matrix for a
% 3-D triangulation. elem2edge is the elementwise pointer from elem to dof
% indices. edge is the edge matrix. The orientation of edge is
% from the node with smaller index to bigger one. The orientation of local
% edges of a tetrahedron formed by [1 2 3 4] is lexicographic order:
%           [1 2], [1 3], [1 4], [2 3], [2 4], [3 4]
%
% [elem2edge,edge,elem2edgeSign] = DOF3EDGE(elem) also output elem2edgeSign
% which records the consistency of the local and global edge orientation.
% If elem is ascend ordered, then elem2edgeSign is 1 and do not needed.
%
% See also dof3edge.
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

NT = size(elem,1);
totalEdge = int32([elem(:,[1 2]); elem(:,[1 3]); elem(:,[1 4]); ...
                    elem(:,[2 3]); elem(:,[2 4]); elem(:,[3 4])]);
sortedTotalEdge = sort(totalEdge,2);
matlabversion = version;
if str2double(matlabversion(end-5:end-2)) > 2012
    [edge, tempvar, j] = unique(sortedTotalEdge,'rows','legacy'); %#ok<ASGLU>
else
    [edge, tempvar, j] = unique(sortedTotalEdge,'rows'); %#ok<ASGLU>
end
elem2edge = uint32(reshape(j,NT,6));
direction = ones(6*NT,1);
idx = (totalEdge(:,1)>totalEdge(:,2));
direction(idx) = -1;
elem2edgeSign = reshape(direction,NT,6);

repo: Aero-matlab


text:
% pdfbdenoiseimagethmt.m
% written by: Duncan Po
% Date: January 14, 2004
% Clean the noisy image using the model and state probabilities of the
% noisy image. The image is assumed to be square.
% Usage: cleanimage = pdfbdenoiseimage( model, stateprob, nvar, nc, imname, imformat)
% Inputs:   model       - model of the noisy image
%           stateprob   - state probabilities of the noisy image
%           nvar        - noise variance, normalized to the image range
%                           (ie. should be between 0 and 1)
%           nc          - noise variance in contourlet domain. Input '' if
%                           unknown
%           imname      - name of the image file
%           imformat    - format of the image file (e.g. 'gif')
% Output:   cleanimage  - the denoised image

function cleanimage = pdfbdenoiseimage( model, stateprob, nvar, nc, imname, imformat)

pyrfilter = '9-7';
dirfilter = 'pkva';
nlevel = model{1}.nlevels;
for lev = 1:nlevel
    levndir(lev) = log2(length(model{1}.stdv{lev})*length(model));
end;

if nargin == 6
    coef = contourlet(pyrfilter, dirfilter, levndir, imname, imformat);
elseif nargin == 5
    coef = contourlet(pyrfilter, dirfilter, levndir, imname);
end;

% Verify image is square and obtain image dimensions
[nrow, ncol] = size(coef{2}{1});
imagestats = imfinfo(imname, imformat);
imdim = imagestats.Width;
if imdim ~= imagestats.Height
    error('Image must be square.');
end

for state = 1:size(stateprob{1}{1},2)
    for lev = 1:nlevel+1
        newstateprob{state}{lev} = [];
        covariance{state}{lev} = [];
    end;
end;

% process stateprob, and covariance so that it has the same structure as coef
for state = 1:size(stateprob{1}{1},2)
    for dir = 1:2.^levndir(1)
        for scale = 1:length(stateprob{dir})
            tempsp{scale} = stateprob{dir}{scale}(:,state);
        end;
        tempstateprob = tree2contourlet(tempsp, dir, levndir, nrow, ncol);

        for lev = 1:nlevel
            newstateprob{state}{lev+1} = [newstateprob{state}{lev+1} tempstateprob{lev}];
            for subdir = 1:length(model{dir}.stdv{lev})
               covariance{state}{lev+1} = [covariance{state}{lev+1} ...
                     {ones(size(coef{lev+1}{(dir-1)*(2.^(levndir(lev)-levndir(1)))+subdir},1),...
                        size(coef{lev+1}{(dir-1)*(2.^(levndir(lev)-levndir(1)))+subdir},2))*...
                        model{dir}.stdv{lev}{subdir}(state)*model{dir}.stdv{lev}{subdir}(state)}];
            end;
        end;
    end;
end;

% calculates the noise variance in contourlet domain
if isempty(nc)
    nc = contournc(nvar, pyrfilter, dirfilter, levndir, imdim);
end;

% computes the attenuator for each coefficient and state
for state = 1:length(covariance)
   for lev = 2:nlevel+1
      for dir = 1:length(covariance{state}{lev})
         [covrow, covcol] = size(covariance{state}{lev}{dir});
         covariance{state}{lev}{dir}=max(covariance{state}{lev}{dir}...
            -nc{lev}{dir}, zeros(covrow, covcol))./covariance{state}{lev}{dir};
      end;
   end;
end;

% multiplies each contourlet coefficient with its corresponding state probabilities
for state = 1:length(newstateprob)
   for lev = 2:nlevel+1
      for dir = 1:length(newstateprob{state}{lev})
         newcoef{state}{lev}{dir} = newstateprob{state}{lev}{dir}...
            .*coef{lev}{dir};
      end;
   end;
end;

% we don't denoise the scaling function since its SNR is high
cleancoef{1} = coef{1};

% initialize the clean coefficients
for lev = 2:nlevel+1
    for dir = 1:length(newstateprob{state}{lev})
        cleancoef{lev}{dir} = zeros(size(coef{lev}{dir},1), size(coef{lev}{dir},2));
    end;
end;

% denoise
for state = 1:length(newstateprob)
   for lev = 2:nlevel+1
      for dir = 1:length(newstateprob{state}{lev})
         cleancoef{lev}{dir}=cleancoef{lev}{dir}+covariance{state}{lev}{dir}...
            .*newcoef{state}{lev}{dir};
      end;
   end;
end;

cleanimage = pdfbrec(cleancoef, pyrfilter, dirfilter);
figure;
imshow(uint8(cleanimage));


text:
function [kept,added]= ...
  setDifferencesScoreFeatures(scoreFeaturesOld,scoreFeaturesNew)
% Given an old and new scoreFeatures, computes a boolean array that tells
% which elements of the old one are retained in the new one (kept), and a
% boolean array indicating which elements of the new one are not present in
% the old one (added).

% Get the list of (absolute) files names for each scoreFeatures array
fileNameListOld={scoreFeaturesOld(:).classifierfile};
fileNameListNew={scoreFeaturesNew(:).classifierfile};

% define a local function that checks whether a file name is in the new file
% name list
inFileNameListNew=@(fileName)(ismember(fileName,fileNameListNew));

% mark the ones in the old list that are not in the new list (these are
% ones that will be deleted)
kept=cellfun(inFileNameListNew,fileNameListOld);

% define a local function that checks whether a file name is in the old file
% name list
notInFileNameListOld=@(fileName)(~ismember(fileName,fileNameListOld));

% mark the ones in the new list that were not in the old list (added
% ones)
added=cellfun(notInFileNameListOld,fileNameListNew);

end

author: kristinbranson

repo: JAABA


text:
function  FRIEDEL= Tpm_FR(P,X,G,DH,ERH)
% TPM_FR  Friedel two-phase multiplier
% TPM_FR(P,X,G,DH,ERH) Returns the Friedel two-phase
% multipier for a steam-water system
%  Called function: h2o_rhof(P), h2o_rhog(P), h2o_muf(P),
%                   h2o_mug(P), h2o_sigma(P), h2o_rhotp(P,X),
%                   ffcw(RE,DH,ERH),
%  Required Inputs are: P  - pressure (kPa)
%                       X  - quality (fraction)
%                       G  - mass flux (kg/m^2s)
%                       DH - hydraulic diameter (m)
%                       ERH - equivalent roughness height (m)
% ---------------------------------------------------------------
% The MATLAB function was created by Tibor Balint, December 1998
% TBoreal Research Corporation, Toronto, Ont. Canada
% (tibor@netcom.ca) and also, University of Warwick, UK
% ---------------------------------------------------------------

format long g;  % set the format of the calculations

% set the density and viscosity properties
RF=h2o_rhof(P); %sat. fluid density
RG=h2o_rhog(P); %sat. steam density
MUF=h2o_muf(P); %dynamic viscosity of sat. fluid
MUG=h2o_mug(P); %dynamic viscosity of sat. steam
SIGMA=h2o_sigma(P); %surface tension
VFG=(1/RG)-(1/RF); %diff. between steam and fluid spec.volumes
VF=1/RF; %saturated fluid specific volume

% the two-phase specific density
RHOTP=h2o_rhotp(P,X);

% the two-phase specific volume
VTP=1/RHOTP;

% the Froude number
FR=(G^2)*(VTP^2)/(9.81*DH);

% the Weber number
WE=(G^2)*DH*VTP/SIGMA;

% the Reynolds number FOR LIQUID
REF=G*DH/MUF;

% the Reynolds number FOR STEAM
REG=G*DH/MUG;

% the Colebrook-White friction factor FOR LIQUID
FCWF=ffcw(REF,DH,ERH);

% the Colebrook-White friction factor FOR STEAM
FCWG=ffcw(REG,DH,ERH);

% Calculate the Friedel correlation
% the first part of the Friedel correlation
A=((1-X)^2)+(X^2)*(RF*FCWG)/(RG*FCWF);

%  THE FRIEDEL TWO-PHASE FRICTION MULTIPLIER
FRIEDEL=A+3.21*(X^0.78)*((1-X)^0.224)*((RF/RG)^0.91)*((MUG/MUF)^0.19)*...
       ((1-(MUG/MUF))^0.7)/((FR^0.0454)*(WE^0.035));

return


text:
function processing = FeaturePipe2GPU(processing)

for i=1:length(processing)
    if isfield(processing{i}, 'transform')
        processing{i}.transform = gpuArray(processing{i}.transform);
    end
    if isfield(processing{i}, 'bias')
        processing{i}.bias = gpuArray(processing{i}.bias);
    end
end

author: singaxiong

repo: SignalGraph


text:
function [newnode,newelem,newface]=meshrefine(node,elem,varargin)
%
% [newnode,newelem,newface]=meshrefine(node,elem,face,opt)
%
% refine a tetrahedral mesh by adding new nodes or constraints
%
% author: Qianqian Fang (fangq<at> nmr.mgh.harvard.edu)
%
% input parameters:
%      node: existing tetrahedral mesh node list
%      elem: existing tetrahedral element list
%      face: (optional) existing tetrahedral mesh surface triangle list
%      opt:  options for mesh refinement:
%        if opt is a Nx3 array, opt is treated as a list of new nodes to
%          be inserted into the mesh (must be located on the surface or inside)
%        if opt is a struct, it can have the following fields:
%          opt.newnode: same as setting opt to an Nx3 array
%          opt.reratio: radius-edge ratio, by default, iso2mesh uses 1.414
%          opt.maxvol: maximum element volume
%
% outputs:
%      newnode: node coordinates of the tetrahedral mesh
%      newelem: element list of the tetrahedral mesh
%      newface: mesh surface element list of the tetrahedral mesh
%             the last column denotes the boundary ID
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

exesuff=getexeext;
exesuff=fallbackexeext(exesuff,'tetgen');

newpt=[];
opt=struct;
if(length(varargin)==1)
	face=[];
	if(isstruct(varargin{1}))
		opt=varargin{1};
    else
		newpt=varargin{1};
	end
elseif(length(varargin)>=2)
    face=varargin{1};
    if(isstruct(varargin{2}))
        opt=varargin{2};
    else
        newpt=varargin{2};
    end
else
	error('meshrefine requires at least 3 inputs');
end
if(isstruct(opt) && isfield(opt,'newnode'))
    newpt=opt.newnode;
end

% call tetgen to create volumetric mesh
deletemeshfile(mwpath('pre_refine.*'));
deletemeshfile(mwpath('post_refine.*'));

moreopt='';
setquality=0;
if(isstruct(opt) && isfield(opt,'reratio'))
	moreopt=[moreopt sprintf(' -q %.10f ',opt.reratio)];
	setquality=1;
end
if(isstruct(opt) && isfield(opt,'maxvol'))
    moreopt=[moreopt sprintf(' -a %.10f ',opt.maxvol)];
end

if(~isempty(newpt))
	savetetgennode(newpt,mwpath('pre_refine.1.a.node'));
	moreopt=' -i ';
end
if(size(elem,2)==3 && setquality==0)
    if(~isempty(newpt))
        error('inserting new point can not be used for surfaces');
    end
    nedge=savegts(node, elem,mwpath('pre_refine.gts'));
    exesuff=fallbackexeext(getexeext,'gtsrefine');
elseif(size(elem,2)==3)
    savesurfpoly(node,elem,[],[],[],[],mwpath('pre_refine.poly'));
else
    savetetgennode(node, mwpath('pre_refine.1.node'));
    savetetgenele (elem, mwpath('pre_refine.1.ele'));
end

fprintf(1,'refining the input mesh ...\n');

if(size(elem,2)==3 && setquality==0)
    if(isstruct(opt) && isfield(opt,'scale'))
        moreopt=sprintf('%s -n %d ',moreopt,round(nedge*opt.scale));
    else
        error('you must give opt.scale value for refining a surface');
    end
end
if(isstruct(opt) && isfield(opt,'moreopt'))
	moreopt=[moreopt opt.moreopt];
end

if(size(elem,2)==3 && setquality==0)
    system([' "' mcpath('gtsrefine') exesuff '" ' moreopt ' < "' ...
          mwpath('pre_refine.gts') '" > "' mwpath('post_refine.gts') '"']);
    [newnode,newelem]=readgts(mwpath('post_refine.gts'));
    newface=newelem;
elseif(size(elem,2)==3)
    system([' "' mcpath('tetgen') exesuff '" ' moreopt ' -p -A "' mwpath('pre_refine.poly') '"']);
    [newnode,newelem,newface]=readtetgen(mwpath('pre_refine.1'));
else
    system([' "' mcpath('tetgen') exesuff '" ' moreopt ' -r "' mwpath('pre_refine.1') '"']);
    [newnode,newelem,newface]=readtetgen(mwpath('pre_refine.2'));
end

% read in the generated mesh

fprintf(1,'mesh refinement is complete\n');

repo: REST


text:
function zz=lpcss2zz(ss)
%LPCSS2ZZ Convert s-place poles to z-plane poles ZZ=(SS)
%the s-plane is in units of Normalized Hz and so the imaginary part
% of each ss() value is in the range +-0.5
%
% If you multiply ss by the sample frequency, a formant with
% frequency f and bandwidth b will give an s-plane pole-pair
% of approximately -b/2 +- j f

%      Copyright (C) Mike Brookes 1997
%      Version: $Id: lpcss2zz.m 8210 2016-07-20 20:58:39Z dmb $
%
%   VOICEBOX is a MATLAB toolbox for speech processing.
%   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   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 can obtain a copy of the GNU General Public License from
%   http://www.gnu.org/copyleft/gpl.html or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

zz=exp(2*pi*ss);

author: jtkim-kaist

repo: VAD

save_path: github-repos/MATLAB/jtkim-kaist-VAD