diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..af3151e --- /dev/null +++ b/.gitignore @@ -0,0 +1,2 @@ +# Ignore temporary matlab ASV files +*.asv diff --git a/data/BodyParts3D/mat/FMA16586.mat b/data/BodyParts3D/mat/FMA16586.mat index f02bf20..df77966 100644 Binary files a/data/BodyParts3D/mat/FMA16586.mat and b/data/BodyParts3D/mat/FMA16586.mat differ diff --git a/data/BodyParts3D/mat/FMA24474.mat b/data/BodyParts3D/mat/FMA24474.mat index 46ac0ad..cb27c8f 100644 Binary files a/data/BodyParts3D/mat/FMA24474.mat and b/data/BodyParts3D/mat/FMA24474.mat differ diff --git a/data/BodyParts3D/mat/FMA24477.mat b/data/BodyParts3D/mat/FMA24477.mat index bfb5d73..52091e3 100644 Binary files a/data/BodyParts3D/mat/FMA24477.mat and b/data/BodyParts3D/mat/FMA24477.mat differ diff --git a/data/BodyParts3D/mat/FMA24480.mat b/data/BodyParts3D/mat/FMA24480.mat index 35a5362..e83c719 100644 Binary files a/data/BodyParts3D/mat/FMA24480.mat and b/data/BodyParts3D/mat/FMA24480.mat differ diff --git a/data/BodyParts3D/mat/FMA24486.mat b/data/BodyParts3D/mat/FMA24486.mat index 400694b..5450692 100644 Binary files a/data/BodyParts3D/mat/FMA24486.mat and b/data/BodyParts3D/mat/FMA24486.mat differ diff --git a/data/BodyParts3D/mat/FMA7163.mat b/data/BodyParts3D/mat/FMA7163.mat index 73e0bbf..b41bb21 100644 Binary files a/data/BodyParts3D/mat/FMA7163.mat and b/data/BodyParts3D/mat/FMA7163.mat differ diff --git a/data/BodyParts3D/post/skin_leg_right.mat b/data/BodyParts3D/post/skin_leg_right.mat index 32ecdf1..08fa4d1 100644 Binary files a/data/BodyParts3D/post/skin_leg_right.mat and b/data/BodyParts3D/post/skin_leg_right.mat differ diff --git a/docs/docTools/addCodeFooter.m b/docs/docTools/addCodeFooter.m new file mode 100644 index 0000000..7251341 --- /dev/null +++ b/docs/docTools/addCodeFooter.m @@ -0,0 +1,69 @@ +clear; close all; clc; + +%% + +toolboxPath=fileparts(fileparts(fileparts(mfilename('fullpath')))); +docPath=fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'docs'); +libPath=fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'mcode'); +docToolsPath=fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'docs','docTools'); + +pathNames={docToolsPath,libPath,docPath}; + +%% Prepare boilerplate + +licenseBoilerPlate=fullfile(toolboxPath,'licenseBoilerPlate.txt'); +[T]=txtfile2cell(licenseBoilerPlate); + +footerTargetText='% _*SimuLimb footer text*_ '; + +%% + +licenseLink='https://github.com/SimuLimb/simulateAmputation.m/blob/main/LICENSE'; + +%Add comment symbols +for q=1:1:numel(T) + T{q}=['% ',T{q}]; +end + +%Add target header +T=[{'%% '};{footerTargetText};{'% '};{['% License: <',licenseLink,'>']};T(1:end)]; + +%% + +fileExtension='.m'; + +%% + +numPaths=numel(pathNames); +for q_path=1:1:numPaths + pathName=pathNames{q_path}; + files = dir(fullfile(pathName,'*.m')); + files={files(1:end).name}; + files=sort(files(:)); + numFiles=numel(files); + for q_file=1:1:numFiles + fileName=fullfile(pathName,files{q_file}); + [T_now]=txtfile2cell(fileName); + T_now(end+1:end+numel(T))=T; + cell2txtfile(fileName,T_now,0); + end +end + +%% +% _*SimuLimb footer text*_ +% +% License: +% +% Copyright (C) 2006-2021 Kevin Mattheus Moerman and the SimuLimb contributors +% +% Licensed under the Apache License, Version 2.0 (the "License"); +% you may not use this file except in compliance with the License. +% You may obtain a copy of the License at +% +% http://www.apache.org/licenses/LICENSE-2.0 +% +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" BASIS, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. diff --git a/docs/docTools/removeCodeFooter.m b/docs/docTools/removeCodeFooter.m new file mode 100644 index 0000000..5606fa0 --- /dev/null +++ b/docs/docTools/removeCodeFooter.m @@ -0,0 +1,58 @@ +clear; close all; clc; + +%% Path names + +toolboxPath=fileparts(fileparts(fileparts(mfilename('fullpath')))); +docPath=fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'docs'); +libPath=fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'mcode'); +docToolsPath=fullfile(fileparts(fileparts(fileparts(mfilename('fullpath')))),'docs','docTools'); + +pathNames={docToolsPath,libPath,docPath}; + +%% + +fileExtension='.m'; %Extension to remove boilerplate from + +footerTargetStart='% _*SimuLimb footer text*_ '; %Target for removal + +%N.B. Removal is up to end from the start so make sure the target is +%appropriately set!!!!! + +%% Removing footer text + +numPaths=numel(pathNames); +for q_path=1:1:numPaths + pathName=pathNames{q_path}; + files = dir(fullfile(pathName,'*.m')); + files={files(1:end).name}; + files=sort(files(:)); + numFiles=numel(files); + for q_file=1:1:numFiles + fileName=fullfile(pathName,files{q_file}); + [T_now]=txtfile2cell(fileName); + targetStartIndex = find(strcmp(footerTargetStart,T_now)); + if ~isempty(targetStartIndex) + targetStartIndex=targetStartIndex(end); %Keep last occurance + T_now=T_now(1:targetStartIndex-1-1); %NB -1 is used to remove %% above target + cell2txtfile(fileName,T_now,0,0); + end + end +end +%% +% _*SimuLimb footer text*_ +% +% License: +% +% Copyright (C) 2006-2021 Kevin Mattheus Moerman and the SimuLimb contributors +% +% Licensed under the Apache License, Version 2.0 (the "License"); +% you may not use this file except in compliance with the License. +% You may obtain a copy of the License at +% +% http://www.apache.org/licenses/LICENSE-2.0 +% +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" BASIS, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. diff --git a/licenseBoilerPlate.txt b/licenseBoilerPlate.txt new file mode 100644 index 0000000..44d31b8 --- /dev/null +++ b/licenseBoilerPlate.txt @@ -0,0 +1,14 @@ + +Copyright (C) 2006-2021 Kevin Mattheus Moerman and the SimuLimb contributors + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. diff --git a/mcode/A_BodyParts3D_convert_STL_to_MAT.m b/mcode/A_BodyParts3D_convert_STL_to_MAT.m index fc46db4..d248320 100644 --- a/mcode/A_BodyParts3D_convert_STL_to_MAT.m +++ b/mcode/A_BodyParts3D_convert_STL_to_MAT.m @@ -1,69 +1,87 @@ -clear; close all; clc; - -%% Description -% This code parses all STL files in the /stl folder and exports them as MAT -% files to the /mat folder. -% -% This code requires the GIBBON MATLAB toolbox -% - -%% Control parameters - -% Path names -projectFolder = fileparts(fileparts(mfilename('fullpath'))); %Main code path -loadFolder=fullfile(projectFolder,'data','BodyParts3D','stl'); %The STL loading folder -saveFolder=fullfile(projectFolder,'data','BodyParts3D','mat'); %The MAT saving folder - -%% - -%Get list of all STL files -fileList = dir(fullfile(loadFolder,['*','stl'])); -fileList={fileList(1:end).name}; - -Dc=readtable(fullfile(projectFolder,'data','BodyParts3D','FMA_ID_label_obj.csv')); -P=lower(Dc.FMAID); - -%Convert all STL files to MAT -numFiles=numel(fileList); %Number of files to parse -for q=1:1:numFiles %Loop over all files - %Get file name - [~,fileNameClean,~]=fileparts(fileList{q}); %File name without path or extension - - %Prepare STL file name - fileNameLoad=fullfile(loadFolder,fileList{q}); %STL file name - - %Prepare save file name - fileNameSave=fullfile(saveFolder,[fileNameClean,'.mat']); %MAT file name - - %Get FMAID code - FMAID=sscanf(fileNameClean,'FMA%d'); - - %Get preferred name - logicFMAID=P==FMAID; - preferredName=Dc.preferredName{logicFMAID}; - - %Parse STL import - disp(['Parsing file ', num2str(q),' of ',num2str(numFiles),', ',sprintf('%3.0f',round(100*q/numFiles)),'% done, ',fileNameClean,', ',preferredName]) - try %MATLAB's stlread - TR = stlread(fileNameLoad); - v=TR.Points; - f=TR.ConnectivityList; - catch %GIBBON's importer - [stlStruct] = import_STL_bin(fileNameLoad); - v=stlStruct.solidVertices{1}; - f=stlStruct.solidFaces{1}; - end - - %Merge vertices - [f,v]=mergeVertices(f,v); - - %Build model structure - model.sourceName=fileList{q}; - model.preferredName=preferredName; - model.faces=f; - model.vertices=v; - - %Export MAT file - save(fileNameSave,'-struct','model'); -end - +clear; close all; clc; + +%% Description +% This code parses all STL files in the /stl folder and exports them as MAT +% files to the /mat folder. +% +% This code requires the GIBBON MATLAB toolbox +% + +%% Control parameters + +% Path names +projectFolder = fileparts(fileparts(mfilename('fullpath'))); %Main code path +loadFolder=fullfile(projectFolder,'data','BodyParts3D','stl'); %The STL loading folder +saveFolder=fullfile(projectFolder,'data','BodyParts3D','mat'); %The MAT saving folder + +%% + +%Get list of all STL files +fileList = dir(fullfile(loadFolder,['*','stl'])); +fileList={fileList(1:end).name}; + +Dc=readtable(fullfile(projectFolder,'data','BodyParts3D','FMA_ID_label_obj.csv')); +P=lower(Dc.FMAID); + +%Convert all STL files to MAT +numFiles=numel(fileList); %Number of files to parse +for q=1:1:numFiles %Loop over all files + %Get file name + [~,fileNameClean,~]=fileparts(fileList{q}); %File name without path or extension + + %Prepare STL file name + fileNameLoad=fullfile(loadFolder,fileList{q}); %STL file name + + %Prepare save file name + fileNameSave=fullfile(saveFolder,[fileNameClean,'.mat']); %MAT file name + + %Get FMAID code + FMAID=sscanf(fileNameClean,'FMA%d'); + + %Get preferred name + logicFMAID=P==FMAID; + preferredName=Dc.preferredName{logicFMAID}; + + %Parse STL import + disp(['Parsing file ', num2str(q),' of ',num2str(numFiles),', ',sprintf('%3.0f',round(100*q/numFiles)),'% done, ',fileNameClean,', ',preferredName]) + try %MATLAB's stlread + TR = stlread(fileNameLoad); + v=TR.Points; + f=TR.ConnectivityList; + catch %GIBBON's importer + [stlStruct] = import_STL_bin(fileNameLoad); + v=stlStruct.solidVertices{1}; + f=stlStruct.solidFaces{1}; + end + + %Merge vertices + [f,v]=mergeVertices(f,v); + + %Build model structure + model.sourceName=fileList{q}; + model.preferredName=preferredName; + model.faces=f; + model.vertices=v; + + %Export MAT file + save(fileNameSave,'-struct','model'); +end + +%% +% _*SimuLimb footer text*_ +% +% License: +% +% Copyright (C) 2006-2021 Kevin Mattheus Moerman and the SimuLimb contributors +% +% Licensed under the Apache License, Version 2.0 (the "License"); +% you may not use this file except in compliance with the License. +% You may obtain a copy of the License at +% +% http://www.apache.org/licenses/LICENSE-2.0 +% +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" BASIS, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. diff --git a/mcode/B_BodyParts3D_remesh_surfaces.m b/mcode/B_BodyParts3D_remesh_surfaces.m index 10041a2..f8b8f23 100644 --- a/mcode/B_BodyParts3D_remesh_surfaces.m +++ b/mcode/B_BodyParts3D_remesh_surfaces.m @@ -1,76 +1,95 @@ -clear; close all; clc; - -%% Description -% This code parses a selection of MAT files in the /mat folder for surface -% meshes, and remeshes them at a desired point spacing. The remeshed -% surfaces are exported to the /post folder. -% -% This code requires the GIBBON MATLAB toolbox -% - -%% Plotting settings -lineWidth=0.5; - -%% Control parameters - -% Path names -projectFolder = fileparts(fileparts(mfilename('fullpath'))); %Main code path -loadFolder=fullfile(projectFolder,'data','BodyParts3D','mat'); %The MAT loading folder -saveFolder=fullfile(projectFolder,'data','BodyParts3D','post'); %The MAT saving folder for processed data - -fileNames_FMA={'FMA16586','FMA24474','FMA24477','FMA24480','FMA24486'}; - -pointSpacings=[4 4 4 3 3]; - -optionStructRemesh.disp_on=0; % Turn off command window text display - -saveOn=0; - -%% - -for q=1:1:numel(fileNames_FMA) - - %% Import mesh - fileName_FMA=fileNames_FMA{q}; - fileName_mat=fullfile(loadFolder,[fileName_FMA,'.mat']); - - model=load(fileName_mat); - F=model.faces; - V=model.vertices; - - %% Remesh - optionStructRemesh.pointSpacing=pointSpacings(q); %Set desired point spacing - [Fn,Vn]=ggremesh(F,V,optionStructRemesh); - - %% Visualisation - - cFigure; - gtitle([fileName_FMA,' ',model.preferredName]) - subplot(1,2,1); - title('Raw'); - gpatch(F,V,'w','r',1,lineWidth); - axisGeom; - camlight headlight; - - subplot(1,2,2); - title('Remeshed'); - gpatch(Fn,Vn,'w','g',1,lineWidth); - axisGeom; - camlight headlight; - gdrawnow; - - %% Saving - - modelNew.source=model; - modelNew.faces=Fn; - modelNew.vertices=Vn; - modelNew.pointSpacing=pointSpacings(q); - - if saveOn==1 - %Create save name with lowercase letters and underscores instead of spaces - saveNameMesh=regexprep(lower(model.preferredName),' ','_'); - saveName_mat=fullfile(saveFolder,[saveNameMesh,'.mat']); - save(saveName_mat,'-struct','modelNew') - end - -end \ No newline at end of file +clear; close all; clc; + +%% Description +% This code parses a selection of MAT files in the /mat folder for surface +% meshes, and remeshes them at a desired point spacing. The remeshed +% surfaces are exported to the /post folder. +% +% This code requires the GIBBON MATLAB toolbox +% + +%% Plotting settings +lineWidth=0.5; + +%% Control parameters + +% Path names +projectFolder = fileparts(fileparts(mfilename('fullpath'))); %Main code path +loadFolder=fullfile(projectFolder,'data','BodyParts3D','mat'); %The MAT loading folder +saveFolder=fullfile(projectFolder,'data','BodyParts3D','post'); %The MAT saving folder for processed data + +fileNames_FMA={'FMA16586','FMA24474','FMA24477','FMA24480','FMA24486'}; + +pointSpacings=[4 4 4 3 3]; + +optionStructRemesh.disp_on=0; % Turn off command window text display + +saveOn=0; + +%% + +for q=1:1:numel(fileNames_FMA) + + %% Import mesh + fileName_FMA=fileNames_FMA{q}; + fileName_mat=fullfile(loadFolder,[fileName_FMA,'.mat']); + + model=load(fileName_mat); + F=model.faces; + V=model.vertices; + + %% Remesh + optionStructRemesh.pointSpacing=pointSpacings(q); %Set desired point spacing + [Fn,Vn]=ggremesh(F,V,optionStructRemesh); + + %% Visualisation + + cFigure; + gtitle([fileName_FMA,' ',model.preferredName]) + subplot(1,2,1); + title('Raw'); + gpatch(F,V,'w','r',1,lineWidth); + axisGeom; + camlight headlight; + + subplot(1,2,2); + title('Remeshed'); + gpatch(Fn,Vn,'w','g',1,lineWidth); + axisGeom; + camlight headlight; + gdrawnow; + + %% Saving + + modelNew.source=model; + modelNew.faces=Fn; + modelNew.vertices=Vn; + modelNew.pointSpacing=pointSpacings(q); + + if saveOn==1 + %Create save name with lowercase letters and underscores instead of spaces + saveNameMesh=regexprep(lower(model.preferredName),' ','_'); + saveName_mat=fullfile(saveFolder,[saveNameMesh,'.mat']); + save(saveName_mat,'-struct','modelNew') + end + +end + +%% +% _*SimuLimb footer text*_ +% +% License: +% +% Copyright (C) 2006-2021 Kevin Mattheus Moerman and the SimuLimb contributors +% +% Licensed under the Apache License, Version 2.0 (the "License"); +% you may not use this file except in compliance with the License. +% You may obtain a copy of the License at +% +% http://www.apache.org/licenses/LICENSE-2.0 +% +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" BASIS, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. diff --git a/mcode/C_BodyParts3D_isolateLegOuterSkinSurface.asv b/mcode/C_BodyParts3D_isolateLegOuterSkinSurface.asv deleted file mode 100644 index 2a9f754..0000000 --- a/mcode/C_BodyParts3D_isolateLegOuterSkinSurface.asv +++ /dev/null @@ -1,97 +0,0 @@ -clear; close all; clc; - -%% Description -% This code processes the MAT file for the skin surface (found in the /mat -% folder), since it contains both the inner and outer skin surfaces. The -% inner surface is removed and the outer surface is processed to produce a -% single "water tight" mesh. The surface is also remeshed using ggremesh. -% The remeshed surface is exported to the /post folder. -% -% This code requires the GIBBON MATLAB toolbox -% - -%% Plotting settings -lineWidth=0.5; - -%% Control parameters - -% Path names -projectFolder = fileparts(fileparts(mfilename('fullpath'))); %Main code path -loadFolder=fullfile(projectFolder,'data','BodyParts3D','mat'); %The MAT loading folder -saveFolder=fullfile(projectFolder,'data','BodyParts3D','post'); %The MAT saving folder for processed data - -fileName_FMA='FMA7163'; -saveNameMesh='skin_leg_right'; - -pointSpacing=4; -cutHeight=675; -saveOn=1; - -%% - -fileName_mat=fullfile(loadFolder,[fileName_FMA,'.mat']); - -model=load(fileName_mat); -F=model.faces; -V=model.vertices; - -%% - -logicVertices=V(:,1)<0 & V(:,3)<=cutHeight+mean(patchEdgeLengths(F,V)); -logicFaces=any(logicVertices(F),2); -logicFaces=triSurfLogicSharpFix(F,logicFaces); -Fs=F(logicFaces,:); -[Fs,Vs]=patchCleanUnused(Fs,V); - -optionStruct.outputType='label'; -[G,~,groupSize]=tesgroup(Fs,optionStruct); -[~,indMax]=max(groupSize); -logicKeep=G==indMax; -Fs=Fs(logicKeep,:); -[Fs,Vs]=patchCleanUnused(Fs,Vs); - -snapTolerance=mean(patchEdgeLengths(Fs,Vs))/100; -n=vecnormalize([0 0 1]); %Normal direction to plane -P=[0 0 cutHeight]; %Point on plane -[Fc,Vc,~,logicSide]=triSurfSlice(Fs,Vs,[],P,n,snapTolerance); -[Fc,Vc]=patchCleanUnused(Fc(logicSide,:),Vc); - -%% - -shp = alphaShape(Vc,max(patchEdgeLengths(Fc,Vc)),'HoleThreshold',500,'RegionThreshold',1); -[Fa,Va] = boundaryFacets(shp); -[Fa,Va]=patchCleanUnused(Fa,Va); - -optionStruct2.pointSpacing=pointSpacing; %Set desired point spacing -optionStruct2.disp_on=1; % Turn off command window text display -[Fn,Vn]=ggremesh(Fa,Va,optionStruct2); - -%% Visualisation - -cFigure; -subplot(1,2,1); hold on; -gpatch(F,V,'w','none',0.25); -gpatch(Fn,Vn,'bw','none'); -axisGeom; -camlight headlight; - -subplot(1,2,2); hold on; -gpatch(Fn,Vn,'bw','k',1,lineWidth); -axisGeom; -camlight headlight; -gdrawnow; - -%% Store model in structure - -modelNew.source=fileName_FMA; -modelNew.faces=Fn; -modelNew.vertices=Vn; -modelNew.pointSpacing=pointSpacing; -modelNew.cutHeight=cutHeight; - -%% Save model - -if saveOn==1 - saveName_mat=fullfile(saveFolder,[saveNameMesh,'.mat']); - save(saveName_mat,'-struct','modelNew') -end diff --git a/mcode/C_BodyParts3D_isolateLegOuterSkinSurface.m b/mcode/C_BodyParts3D_isolateLegOuterSkinSurface.m index 16fa078..d689767 100644 --- a/mcode/C_BodyParts3D_isolateLegOuterSkinSurface.m +++ b/mcode/C_BodyParts3D_isolateLegOuterSkinSurface.m @@ -1,101 +1,120 @@ -clear; close all; clc; - -%% Description -% This code processes the MAT file for the skin surface (found in the /mat -% folder), since it contains both the inner and outer skin surfaces. The -% inner surface is removed and the outer surface is processed to produce a -% single "water tight" mesh. The surface is also remeshed using ggremesh. -% The remeshed surface is exported to the /post folder. -% -% This code requires the GIBBON MATLAB toolbox -% - -%% Plotting settings -lineWidth=0.5; - -%% Control parameters - -% Path names -projectFolder = fileparts(fileparts(mfilename('fullpath'))); %Main code path -loadFolder=fullfile(projectFolder,'data','BodyParts3D','mat'); %The MAT loading folder -saveFolder=fullfile(projectFolder,'data','BodyParts3D','post'); %The MAT saving folder for processed data - -fileName_FMA='FMA7163'; -saveNameMesh='skin_leg_right'; - -pointSpacing=4; -cutHeight=675; -saveOn=1; - -%% - -fileName_mat=fullfile(loadFolder,[fileName_FMA,'.mat']); - -model=load(fileName_mat); -F=model.faces; -V=model.vertices; - -%% - -logicVertices=V(:,1)<0 & V(:,3)<=cutHeight+mean(patchEdgeLengths(F,V)); -logicFaces=any(logicVertices(F),2); -logicFaces=triSurfLogicSharpFix(F,logicFaces); -Fs=F(logicFaces,:); -[Fs,Vs]=patchCleanUnused(Fs,V); - -optionStruct.outputType='label'; -[G,~,groupSize]=tesgroup(Fs,optionStruct); -[~,indMax]=max(groupSize); -logicKeep=G==indMax; -Fs=Fs(logicKeep,:); -[Fs,Vs]=patchCleanUnused(Fs,Vs); - -snapTolerance=mean(patchEdgeLengths(Fs,Vs))/100; -n=vecnormalize([0 0 1]); %Normal direction to plane -P=[0 0 cutHeight]; %Point on plane -[Fc,Vc,~,logicSide]=triSurfSlice(Fs,Vs,[],P,n,snapTolerance); -[Fc,Vc]=patchCleanUnused(Fc(logicSide,:),Vc); - -%% Construct alpha shape -% This deteriorates surface quality (bridges concave regions) but is a -% termporary "quick-and-dirty" work-around to obtain outer skin surface -% only. - -shp = alphaShape(Vc,max(patchEdgeLengths(Fc,Vc)),'HoleThreshold',500,'RegionThreshold',1); -[Fa,Va] = boundaryFacets(shp); %Get boundary faces of alpha shape -[Fa,Va] = patchCleanUnused(Fa,Va); %Remove unused vertices - -%% Remesh alpha shape using |ggremesh| -optionStructRemesh.pointSpacing=pointSpacing; %Set desired point spacing -optionStructRemesh.disp_on=1; % Turn off command window text display -[Fn,Vn]=ggremesh(Fa,Va,optionStructRemesh); - -%% Visualisation - -cFigure; -subplot(1,2,1); hold on; -gpatch(F,V,'w','none',0.25); -gpatch(Fn,Vn,'bw','none'); -axisGeom; -camlight headlight; - -subplot(1,2,2); hold on; -gpatch(Fn,Vn,'bw','k',1,lineWidth); -axisGeom; -camlight headlight; -gdrawnow; - -%% Store model in structure - -modelNew.source=fileName_FMA; -modelNew.faces=Fn; -modelNew.vertices=Vn; -modelNew.pointSpacing=pointSpacing; -modelNew.cutHeight=cutHeight; - -%% Save model - -if saveOn==1 - saveName_mat=fullfile(saveFolder,[saveNameMesh,'.mat']); - save(saveName_mat,'-struct','modelNew') -end +clear; close all; clc; + +%% Description +% This code processes the MAT file for the skin surface (found in the /mat +% folder), since it contains both the inner and outer skin surfaces. The +% inner surface is removed and the outer surface is processed to produce a +% single "water tight" mesh. The surface is also remeshed using ggremesh. +% The remeshed surface is exported to the /post folder. +% +% This code requires the GIBBON MATLAB toolbox +% + +%% Plotting settings +lineWidth=0.5; + +%% Control parameters + +% Path names +projectFolder = fileparts(fileparts(mfilename('fullpath'))); %Main code path +loadFolder=fullfile(projectFolder,'data','BodyParts3D','mat'); %The MAT loading folder +saveFolder=fullfile(projectFolder,'data','BodyParts3D','post'); %The MAT saving folder for processed data + +fileName_FMA='FMA7163'; +saveNameMesh='skin_leg_right'; + +pointSpacing=4; +cutHeight=675; +saveOn=1; + +%% + +fileName_mat=fullfile(loadFolder,[fileName_FMA,'.mat']); + +model=load(fileName_mat); +F=model.faces; +V=model.vertices; + +%% + +logicVertices=V(:,1)<0 & V(:,3)<=cutHeight+mean(patchEdgeLengths(F,V)); +logicFaces=any(logicVertices(F),2); +logicFaces=triSurfLogicSharpFix(F,logicFaces); +Fs=F(logicFaces,:); +[Fs,Vs]=patchCleanUnused(Fs,V); + +optionStruct.outputType='label'; +[G,~,groupSize]=tesgroup(Fs,optionStruct); +[~,indMax]=max(groupSize); +logicKeep=G==indMax; +Fs=Fs(logicKeep,:); +[Fs,Vs]=patchCleanUnused(Fs,Vs); + +snapTolerance=mean(patchEdgeLengths(Fs,Vs))/100; +n=vecnormalize([0 0 1]); %Normal direction to plane +P=[0 0 cutHeight]; %Point on plane +[Fc,Vc,~,logicSide]=triSurfSlice(Fs,Vs,[],P,n,snapTolerance); +[Fc,Vc]=patchCleanUnused(Fc(logicSide,:),Vc); + +%% Construct alpha shape +% This deteriorates surface quality (bridges concave regions) but is a +% termporary "quick-and-dirty" work-around to obtain outer skin surface +% only. + +shp = alphaShape(Vc,max(patchEdgeLengths(Fc,Vc)),'HoleThreshold',500,'RegionThreshold',1); +[Fa,Va] = boundaryFacets(shp); %Get boundary faces of alpha shape +[Fa,Va] = patchCleanUnused(Fa,Va); %Remove unused vertices + +%% Remesh alpha shape using |ggremesh| +optionStructRemesh.pointSpacing=pointSpacing; %Set desired point spacing +optionStructRemesh.disp_on=1; % Turn off command window text display +[Fn,Vn]=ggremesh(Fa,Va,optionStructRemesh); + +%% Visualisation + +cFigure; +subplot(1,2,1); hold on; +gpatch(F,V,'w','none',0.25); +gpatch(Fn,Vn,'bw','none'); +axisGeom; +camlight headlight; + +subplot(1,2,2); hold on; +gpatch(Fn,Vn,'bw','k',1,lineWidth); +axisGeom; +camlight headlight; +gdrawnow; + +%% Store model in structure + +modelNew.source=fileName_FMA; +modelNew.faces=Fn; +modelNew.vertices=Vn; +modelNew.pointSpacing=pointSpacing; +modelNew.cutHeight=cutHeight; + +%% Save model + +if saveOn==1 + saveName_mat=fullfile(saveFolder,[saveNameMesh,'.mat']); + save(saveName_mat,'-struct','modelNew') +end + +%% +% _*SimuLimb footer text*_ +% +% License: +% +% Copyright (C) 2006-2021 Kevin Mattheus Moerman and the SimuLimb contributors +% +% Licensed under the Apache License, Version 2.0 (the "License"); +% you may not use this file except in compliance with the License. +% You may obtain a copy of the License at +% +% http://www.apache.org/licenses/LICENSE-2.0 +% +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" BASIS, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License. diff --git a/mcode/html/A_BodyParts3D_convert_STL_to_MAT.html b/mcode/html/A_BodyParts3D_convert_STL_to_MAT.html index 3ae1f8b..f7d149e 100644 --- a/mcode/html/A_BodyParts3D_convert_STL_to_MAT.html +++ b/mcode/html/A_BodyParts3D_convert_STL_to_MAT.html @@ -128,77 +128,94 @@ Parsing file 4 of 6, 67% done, FMA24480, Right fibula Parsing file 5 of 6, 83% done, FMA24486, Right patella Parsing file 6 of 6, 100% done, FMA7163, Skin -


Published with MATLAB® R2021a

\ No newline at end of file diff --git a/mcode/html/B_BodyParts3D_remesh_surfaces.html b/mcode/html/B_BodyParts3D_remesh_surfaces.html index 422e55f..85479b4 100644 --- a/mcode/html/B_BodyParts3D_remesh_surfaces.html +++ b/mcode/html/B_BodyParts3D_remesh_surfaces.html @@ -116,83 +116,102 @@ save(saveName_mat,'-struct','modelNew') end 
end
-


Published with MATLAB® R2021a

\ No newline at end of file diff --git a/mcode/html/B_BodyParts3D_remesh_surfaces.png b/mcode/html/B_BodyParts3D_remesh_surfaces.png index c7a1bbc..e50295d 100644 Binary files a/mcode/html/B_BodyParts3D_remesh_surfaces.png and b/mcode/html/B_BodyParts3D_remesh_surfaces.png differ diff --git a/mcode/html/B_BodyParts3D_remesh_surfaces_01.png b/mcode/html/B_BodyParts3D_remesh_surfaces_01.png index 6eb12eb..a37c78a 100644 Binary files a/mcode/html/B_BodyParts3D_remesh_surfaces_01.png and b/mcode/html/B_BodyParts3D_remesh_surfaces_01.png differ diff --git a/mcode/html/B_BodyParts3D_remesh_surfaces_02.png b/mcode/html/B_BodyParts3D_remesh_surfaces_02.png index fbc88a0..6226cef 100644 Binary files a/mcode/html/B_BodyParts3D_remesh_surfaces_02.png and b/mcode/html/B_BodyParts3D_remesh_surfaces_02.png differ diff --git a/mcode/html/B_BodyParts3D_remesh_surfaces_03.png b/mcode/html/B_BodyParts3D_remesh_surfaces_03.png index 5160a0b..693012a 100644 Binary files a/mcode/html/B_BodyParts3D_remesh_surfaces_03.png and b/mcode/html/B_BodyParts3D_remesh_surfaces_03.png differ diff --git a/mcode/html/B_BodyParts3D_remesh_surfaces_04.png b/mcode/html/B_BodyParts3D_remesh_surfaces_04.png index a631182..8b9d3f4 100644 Binary files a/mcode/html/B_BodyParts3D_remesh_surfaces_04.png and b/mcode/html/B_BodyParts3D_remesh_surfaces_04.png differ diff --git a/mcode/html/B_BodyParts3D_remesh_surfaces_05.png b/mcode/html/B_BodyParts3D_remesh_surfaces_05.png index e9418eb..75a8dbe 100644 Binary files a/mcode/html/B_BodyParts3D_remesh_surfaces_05.png and b/mcode/html/B_BodyParts3D_remesh_surfaces_05.png differ diff --git a/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface.html b/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface.html index cbfb3d7..5cdef84 100644 --- a/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface.html +++ b/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface.html @@ -111,9 +111,9 @@ [Fn,Vn]=ggremesh(Fa,Va,optionStructRemesh); 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
------->  Geogram/vorpalite for resmeshing  <------ 01-Sep-2021 17:26:26
-# Export mesh input file.                          01-Sep-2021 17:26:26
-# Run Geomgram/vorpalite.                          01-Sep-2021 17:26:26
+------>  Geogram/vorpalite for resmeshing  <------ 01-Sep-2021 17:42:02
+# Export mesh input file.                          01-Sep-2021 17:42:02
+# Run Geomgram/vorpalite.                          01-Sep-2021 17:42:02
  ______________________________________________________________________________ 
 |                                                                              |
 | o-[config      ] Configuration file name:geogram.ini                         |
@@ -137,11 +137,11 @@
    _______________________
  _/ =====[remeshing]===== \____________________________________________________
 |                                                                              |
-||| o-[Newton      ] Elapsed time: 1.32s                                         |
+||| o-[Newton      ] Elapsed time: 1.29s                                         |
 | o-[Remesh      ] Computing RVD...                                            |
 | o-[Validate    ] (FP64) nb_v:22222 nb_e:0 nb_f:44440 nb_b:0 tri:1 dim:3      |
 |                  Attributes on vertices: point[3]                            |
-| o-[Remesh      ] Elapsed time: 1.65 s                                        |
+| o-[Remesh      ] Elapsed time: 1.64 s                                        |
    ____________________________
  _/ =====[postprocessing]===== \_______________________________________________
 |                                                                              |
@@ -159,11 +159,11 @@
 |                  Attributes on vertices: point[3]                            |
 | o-[I/O         ] Saving file /mnt/data/MATLAB/GIBBON/data/temp/temp_out.obj. |
 |                  ..                                                          |
-| o-[Total time  ] Elapsed time: 3.74 s                                        |
+| o-[Total time  ] Elapsed time: 3.73 s                                        |
 \______________________________________________________________________________/
-# Importing remeshed geometry.                     01-Sep-2021 17:26:30
-# Removing temporary files.                        01-Sep-2021 17:26:31
-# Done!                                            01-Sep-2021 17:26:31
+# Importing remeshed geometry.                     01-Sep-2021 17:42:06
+# Removing temporary files.                        01-Sep-2021 17:42:07
+# Done!                                            01-Sep-2021 17:42:07

Visualisation

cFigure;
 subplot(1,2,1); hold on;
 gpatch(F,V,'w','none',0.25);
@@ -185,109 +185,126 @@
     saveName_mat=fullfile(saveFolder,[saveNameMesh,'.mat']);
     save(saveName_mat,'-struct','modelNew')
 end
-


Published with MATLAB® R2021a

\ No newline at end of file diff --git a/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface.png b/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface.png index 466077c..fd1bf6d 100644 Binary files a/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface.png and b/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface.png differ diff --git a/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface_01.png b/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface_01.png index 2e5ec7c..f252856 100644 Binary files a/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface_01.png and b/mcode/html/C_BodyParts3D_isolateLegOuterSkinSurface_01.png differ diff --git a/mcode/processAmputation.m b/mcode/processAmputation.m index 4aea731..8453312 100644 --- a/mcode/processAmputation.m +++ b/mcode/processAmputation.m @@ -1,556 +1,575 @@ -clear; close all; clc; - -%% Description -% This code processes skin and bone surface meshes for the lower limb to -% simulate amputation. Amputation is "simulated" in the sense that a 3D -% surface mesh is obtained with a shape that may be expected from a typical -% amputation. The current code is for transtibial amputation. -% -% This code requires the GIBBON MATLAB toolbox -% - -%% Control parameters - -projectFolder = fileparts(fileparts(mfilename('fullpath'))); -loadFolder=fullfile(projectFolder,'data','BodyParts3D','post'); -saveFolder=loadFolder; -saveNameGeom='BodyParts3D_right_leg_transtibial_amp'; - -fileNames={'right_femur','right_tibia','right_fibula','right_patella','right_leg_skin'}; - -saveOn=0; - -amputationDistanceTibia=180; -amputationDistanceFibula=amputationDistanceTibia-30; -amputationDistanceSkin=amputationDistanceTibia-25; -distalExcess=25; -boneRoundFactor=1; -radiusEnd=distalExcess+(amputationDistanceTibia-amputationDistanceSkin); -taperHeigth=100; -taperThreshold=15; %Regions more distance than threshold are tappered -taperFraction=0.3; -amputationDistances=[amputationDistanceTibia amputationDistanceFibula amputationDistanceSkin]; -tibiaDistalCutRadialFactor=0.5; -numBezierPoints=50; -nCut=vecnormalize([0 0 1]); %Normal direction to plane -topCropOffset=10; - -%% Loading surfaces into cell array - -%Allocate cell arrays -FT=cell(1,numel(fileNames)); -VT=cell(1,numel(fileNames)); -CT=cell(1,numel(fileNames)); - -%Loop over all surfaces -for q=1:1:numel(fileNames) - - %Import mesh - fileName=fileNames{q}; - fileName_mat=fullfile(loadFolder,[fileName,'.mat']); - model=load(fileName_mat); - F=model.faces; %Faces - V=model.vertices; %Vertices - C=q*ones(size(F,1),1); %Color label - - %Store mesh data in cell array - FT{q}=F; - VT{q}=V; - CT{q}=C; -end - -%% Compute landmarks - -P_patella_centroid=triSurfCentroid(FT{4},VT{4}); %Centroid of the patella - -%% -% Visualize surface and landmarks - -cFigure; hold on; -gpatch(FT,VT,CT,'none',0.5); -plotV(P_patella_centroid,'k.','MarkerSize',35); -axisGeom; camlight headlight; -colormap gjet; icolorbar; -gdrawnow; - -%% Cut surfaces - -%Allocate cell arrays -FT_amp=FT; -VT_amp=VT; -CT_amp=CT; -c=[2 3 5]; %Indices (=labels) for surfaces to cut -for q=1:numel(c) - cutLevelNow=P_patella_centroid(:,3)-amputationDistances(q); - - %Get surface - F=FT_amp{c(q)}; - V=VT_amp{c(q)}; - - %Use triSurfSlice to process cut - snapTolerance=mean(patchEdgeLengths(F,V))/100; - P=[0 0 cutLevelNow]; %Point on plane - [Fc,Vc,~,logicSide]=triSurfSlice(F,V,[],P,nCut,snapTolerance); - [Fc,Vc]=patchCleanUnused(Fc(~logicSide,:),Vc); %Clean-up mesh - Cc=c(q)*ones(size(Fc,1),1); %Adjusted color data (shorter list after cut) - - %Store processed mesh data in cell arrays - FT_amp{c(q)}=Fc; - VT_amp{c(q)}=Vc; - CT_amp{c(q)}=Cc; -end - -%% -% Visualization - -cFigure; hold on; -title('Cut features'); -% gpatch(FT,VT,'w','none',0.1); -gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); -axisGeom; camlight headlight; -colormap gjet; icolorbar; -gdrawnow; - -%% Set up skin taper parameterisation - -%Get tibia end curve -Ft=FT_amp{2}; -Vt=VT_amp{2}; -Ebt=patchBoundary(Ft,Vt); -indBt=edgeListToCurve(Ebt); -indBt=indBt(1:end-1); -P_tibia_end_centroid=mean(Vt(indBt,:),1); - -Fs=FT_amp{5}; -Vs=VT_amp{5}; -D=Vs(:,3)-(min(Vs(:,3))+taperHeigth); -D(D>0)=0; -D=abs(D); -D=D./max(D); -D=D.^2; - -[Dp,indMin]=minDist(Vs(:,[1 2]),Vt(indBt,:)); -Dp(Dp<=taperThreshold)=taperThreshold; -Dp=Dp-min(Dp); -Dp=Dp./max(Dp); - -R=Vt(indBt(indMin),[1 2])-Vs(:,[1 2]); -R=R.*D.*taperFraction.*Dp; - -%% Process skin taper - -Vs1=Vs; %Initialize as original -Vs(:,[1 2])=Vs(:,[1 2])+R; %Push to produce taper - -%Overide skin surface with tapered surface -FT_amp{5}=Fs; -VT_amp{5}=Vs; - -%% -% Visualisation - -% MR=sqrt(sum(R.^2,2)); - -cFigure; hold on; -title('taper morphing'); -gpatch(FT_amp,VT_amp,'w','none',0.5); -% gpatch(Fs,Vs,'w','none',0.5); -quiverVec(Vs1,R); -colormap gjet; colorbar; -axisGeom; camlight headlight; -gdrawnow; - -%% Process skin distal end closure - -P_end=P_tibia_end_centroid-[0 0 distalExcess]; -Ebs=patchBoundary(Fs,Vs); - -[Fsb,Vsb,indBs,Nd,P1,P2,P3,XB,YB,ZB]=smoothBezierClose(Fs,Vs,Ebs,P_end,numBezierPoints); - -%% -% Visualization - -cFigure; hold on; -title('Distal end closure'); -gpatch(FT_amp,VT_amp,'w','none',1); -gpatch(Fs,Vs,'w','none',0.25); -hp=gpatch(Fsb,Vsb,Vsb(:,3),'none',1); hp.FaceColor='interp'; -% gpatch(Ebs,Vs,'none','b',1,3); - -quiverVec(Vs(indBs,:),Nd(indBs,:),radiusEnd/4,'k'); - -plotV(P1,'r.-','MarkerSize',15,'LineWidth',2); -plotV(P2,'g.-','MarkerSize',15,'LineWidth',2); -plotV(P3,'b.-','MarkerSize',15,'LineWidth',2); -for q=1:size(XB,2) - plotV([XB(:,q) YB(:,q) ZB(:,q)],'k.-','LineWidth',0.5,'MarkerSize',5); -end - -plotV(P_end,'k.','MarkerSize',25); - -axisGeom; camlight headlight; -gdrawnow; - -%% -% % Visualization -% -% cFigure; hold on; -% title('Distal end closure'); -% -% hp=gpatch(Fsb,Vsb,Vsb(:,3),'k',1); hp.FaceColor='interp'; -% % gpatch(Ebs,Vs,'none','b',1,3); -% -% axisGeom; camlight headlight; -% gdrawnow; -% -% dfasfa - -%% Merge skin components and remesh - -pointSpacing=mean(patchEdgeLengths(Fs,Vs)); -[Fs,Vs]=joinElementSets({Fs,Fsb},{Vs,Vsb}); -[Fs,Vs]=mergeVertices(Fs,Vs); - -optionStructRemesh.pointSpacing=pointSpacing; %Set desired point spacing -optionStructRemesh.disp_on=0; % Turn off command window text display -[Fs,Vs]=ggremesh(Fs,Vs,optionStructRemesh); - -%% - -FT_amp{5}=Fs; -VT_amp{5}=Vs; -CT_amp{5}=5*ones(size(Fs,1),1); - -%% - -% [FT,VT,CT]=joinElementSets(FT,VT); -% [FT_amp,VT_amp,CT_amp]=joinElementSets(FT_amp,VT_amp); - -%% - -cFigure; -subplot(1,2,1); hold on; -% gpatch(FT,VT,'w','none',0.25); -gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); -axisGeom; camlight headlight; -colormap gjet; icolorbar; - -subplot(1,2,2); hold on; -gpatch(FT_amp,VT_amp,'w','k',1); -axisGeom; camlight headlight; -gdrawnow; - -%% - -Ft=FT_amp{2}; -Vt=VT_amp{2}; - -Ebs=patchBoundary(Ft,Vt); -indBs=edgeListToCurve(Ebs); -indBs=indBs(1:end-1); - -P_mid=mean(Vt(indBs,:),1); - -radiusEnd=boneRoundFactor*mean(minDist(Vt(indBs,:),P_mid))/2; -P_end=P_mid-[0 0 radiusEnd]; - -[Fsb,Vsb,indBs,Nd,P1,P2,P3,XB,YB,ZB]=smoothBezierClose(Ft,Vt,Ebs,P_end,numBezierPoints); - -%% -% Visualization -cFigure; hold on; -gpatch(Ft,Vt,'w','none',0.25); -gpatch(Ebs,Vt,'none','b',1,3); -hp=gpatch(Fsb,Vsb,'w','k',1,1); -% hp.FaceColor='interp'; -plotV(P_tibia_end_centroid,'k.','MarkerSize',35); - -quiverVec(Vt(indBs,:),Nd(indBs,:),radiusEnd/4,'k'); - -plotV(P1,'r.-','MarkerSize',15,'LineWidth',2); -plotV(P2,'g.-','MarkerSize',15,'LineWidth',2); -plotV(P3,'b.-','MarkerSize',15,'LineWidth',2); -for q=1:size(XB,2) - plotV([XB(:,q) YB(:,q) ZB(:,q)],'k.-','LineWidth',0.5,'MarkerSize',5); -end - -axisGeom; camlight headlight; -gdrawnow; - -%% Merge bone components and remesh - -pointSpacing=mean(patchEdgeLengths(Ft,Vt)); -[Ft,Vt]=joinElementSets({Ft,Fsb},{Vt,Vsb}); -[Ft,Vt]=mergeVertices(Ft,Vt); - -optionStructRemesh.pointSpacing=pointSpacing; %Set desired point spacing -optionStructRemesh.disp_on=0; % Turn off command window text display -[Ft,Vt]=ggremesh(Ft,Vt,optionStructRemesh); - -%% - -FT_amp{2}=Ft; -VT_amp{2}=Vt; -CT_amp{2}=2*ones(size(Ft,1),1); - -%% - -cFigure; -subplot(1,2,1); hold on; -% gpatch(FT,VT,'w','none',0.25); -gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); -axisGeom; camlight headlight; -colormap gjet; icolorbar; - -subplot(1,2,2); hold on; -gpatch(FT_amp,VT_amp,'w','none',0.5); -gpatch(Ft,Vt,'w','k',1); -axisGeom; camlight headlight; -gdrawnow; - -%% FIBULA - -Ft=FT_amp{3}; -Vt=VT_amp{3}; - -Ebs=patchBoundary(Ft,Vt); -indBs=edgeListToCurve(Ebs); -indBs=indBs(1:end-1); - -P_mid=mean(Vt(indBs,:),1); - -radiusEnd=boneRoundFactor*mean(minDist(Vt(indBs,:),P_mid))/2; -P_end=P_mid-[0 0 radiusEnd]; - -[Fsb,Vsb,indBs,Nd,P1,P2,P3,XB,YB,ZB]=smoothBezierClose(Ft,Vt,Ebs,P_end,numBezierPoints); - -%% -% Visualization -cFigure; hold on; -gpatch(Ft,Vt,'w','none',0.25); -gpatch(Ebs,Vt,'none','b',1,3); -hp=gpatch(Fsb,Vsb,'w','k',1,1); -% hp.FaceColor='interp'; -plotV(P_tibia_end_centroid,'k.','MarkerSize',35); - -quiverVec(Vt(indBs,:),Nd(indBs,:),radiusEnd/4,'k'); - -plotV(P1,'r.-','MarkerSize',15,'LineWidth',2); -plotV(P2,'g.-','MarkerSize',15,'LineWidth',2); -plotV(P3,'b.-','MarkerSize',15,'LineWidth',2); -for q=1:size(XB,2) - plotV([XB(:,q) YB(:,q) ZB(:,q)],'k.-','LineWidth',0.5,'MarkerSize',5); -end - -axisGeom; camlight headlight; -gdrawnow; - -%% Merge bone components and remesh - -pointSpacing=mean(patchEdgeLengths(Ft,Vt)); -[Ft,Vt]=joinElementSets({Ft,Fsb},{Vt,Vsb}); -[Ft,Vt]=mergeVertices(Ft,Vt); - -optionStructRemesh.pointSpacing=pointSpacing; %Set desired point spacing -optionStructRemesh.disp_on=0; % Turn off command window text display -[Ft,Vt]=ggremesh(Ft,Vt,optionStructRemesh); - -%% - -FT_amp{3}=Ft; -VT_amp{3}=Vt; -CT_amp{3}=3*ones(size(Ft,1),1); - -%% - -cFigure; -subplot(1,2,1); hold on; -% gpatch(FT,VT,'w','none',0.25); -gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); -axisGeom; camlight headlight; -colormap gjet; icolorbar; - -subplot(1,2,2); hold on; -gpatch(FT_amp,VT_amp,'w','none',0.5); -gpatch(Ft,Vt,'w','k',1); -axisGeom; camlight headlight; -gdrawnow; - -%% - -%Get surface -F=FT_amp{5}; -V=VT_amp{5}; -cutLevelNow=max(V(:,3))-topCropOffset; - -%Use triSurfSlice to process cut -snapTolerance=mean(patchEdgeLengths(F,V))/100; - -P=[0 0 cutLevelNow]; %Point on plane -[Fc,Vc,~,logicSide]=triSurfSlice(F,V,[],P,-nCut,snapTolerance); -[Fc,Vc]=patchCleanUnused(Fc(~logicSide,:),Vc); %Clean-up mesh - -% Remesh using ggremesh -optionStruct3.pointSpacing=mean(patchEdgeLengths(F,V)); -optionStruct3.disp_on=0; % Turn off command window text display -optionStruct3.pre.max_hole_area=100; %Max hole area for pre-processing step -optionStruct3.pre.max_hole_edges=0; %Max number of hole edges for pre-processing step - -[Fc,Vc]=ggremesh(Fc,Vc,optionStruct3); - -Ebs=patchBoundary(Fc,Vc); -indBs=edgeListToCurve(Ebs); -indBs=indBs(1:end-1); - - - -Fcs=Fc; -Vcs=Vc; -indBss=indBs; - -%% - -%Get surface -F=FT_amp{1}; -V=VT_amp{1}; - -%Use triSurfSlice to process cut -snapTolerance=mean(patchEdgeLengths(F,V))/100; -P=[0 0 cutLevelNow]; %Point on plane -[Fc,Vc,~,logicSide]=triSurfSlice(F,V,[],P,-nCut,snapTolerance); -[Fc,Vc]=patchCleanUnused(Fc(~logicSide,:),Vc); %Clean-up mesh - -% Remesh using ggremesh -optionStruct3.pointSpacing=mean(patchEdgeLengths(F,V)); -optionStruct3.disp_on=0; % Turn off command window text display -optionStruct3.pre.max_hole_area=100; %Max hole area for pre-processing step -optionStruct3.pre.max_hole_edges=0; %Max number of hole edges for pre-processing step - -[Fc,Vc]=ggremesh(Fc,Vc,optionStruct3); - -Ebs=patchBoundary(Fc,Vc); -indBs=edgeListToCurve(Ebs); -indBs=indBs(1:end-1); - -%% - -z=0.5*(mean(Vc(indBs,3))+mean(Vc(indBs,3))); - -Vcs(indBss,3)=z; -Vc(indBs,3)=z; - -%% - -pointSpacing=mean(patchEdgeLengths(Fcs,Vcs)); -[Ftt,Vtt]=regionTriMesh3D({Vcs(indBss,:),Vc(indBs,:)},pointSpacing,0,'linear'); - -%% - -cFigure; -gpatch(Fcs,Vcs,'gw','k',1); -gpatch(Fc,Vc,'rw','k',1); -gpatch(Ftt,Vtt,'bw','k',0.5); -colormap gjet; -axisGeom; camlight headlight; -gdrawnow; - -%% -FT_amp{1}=fliplr(Fc); %invert femur normals -VT_amp{1}=Vc; -CT_amp{1}=1*ones(size(Fc,1),1); - -FT_amp{5}=Fcs; -VT_amp{5}=Vcs; -CT_amp{5}=5*ones(size(Fcs,1),1); - -FT_amp{6}=Ftt; -VT_amp{6}=Vtt; -CT_amp{6}=6*ones(size(Ftt,1),1); - -%% - -cFigure; hold on; -gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); -% patchNormPlot(FT_amp,VT_amp); -axisGeom; camlight headlight; -colormap gjet; icolorbar; -gdrawnow; - -%% Save model - -if saveOn==1 - saveName_mat=fullfile(saveFolder,[saveNameGeom,'.mat']); - save(saveName_mat,'FT_amp','VT_amp','CT_amp'); -end - -%% - -function [Fsb,Vsb,indBs,Nd,P1,P2,P3,XB,YB,ZB]=smoothBezierClose(Fs,Vs,Ebs,P_end,numPoints) - -indBs=edgeListToCurve(Ebs); -indBs=indBs(1:end-1); - -[~,~,Ns]=patchNormal(Fs,Vs); - -[~,~,NEb]=edgeVec(Ebs,Vs); -NEb=vecnormalize(NEb); -Nd=vecnormalize(cross(NEb,Ns)); - -cPar.n=25; -Nd=patchSmooth(Ebs,Nd,[],cPar); - -%Bezier point set 1 -P1=Vs(indBs,:); -P1_mid=mean(P1,1); - -distanceEnd=sqrt(sum((P_end-P1_mid).^2,2)); - -%Bezier point set 2 -f=1./abs(Nd(indBs,3)); %Extend factor for direction vectors -P2=P1+distanceEnd/2*f(:,ones(1,3)).*Nd(indBs,:); %Position half-way - -%Bezier point set 3 -P3=P2; -P3(:,3)=P_end(:,3); %Shift to bottom -P3=P3-mean(P3,1); %Shift on own centre -[t,r]=cart2pol(P3(:,1),P3(:,2)); %Polar coordinates -[P3(:,1),P3(:,2)]=pol2cart(t,distanceEnd/2*ones(size(r))); %Force constant radius -P3=P3+P_end; %Place centered on end - -%Bezier point set 4 -P4=P_end; - -XB=zeros(numPoints,size(P1,1)); -YB=zeros(numPoints,size(P1,1)); -ZB=zeros(numPoints,size(P1,1)); -for q=1:1:size(P1,1) - p=[P1(q,:); P2(q,:); P3(q,:); P4]; %Control points - V_bezier=bezierCurve(p,numPoints); %Compute bezier curve - XB(:,q)=V_bezier(:,1); - YB(:,q)=V_bezier(:,2); - ZB(:,q)=V_bezier(:,3); -end - -[Fsb,Vsb]=meshToPatch(XB,YB,ZB,1); -[Fsb,Vsb]=quad2tri(Fsb,Vsb); -[Fsb,Vsb]=mergeVertices(Fsb,Vsb); -[Fsb,Vsb]=patchCleanUnused(Fsb,Vsb); - -end - -%% - -function [F,V]=meshToPatch(X,Y,Z,closeSection) - -%Create quad patch data -[F,V] = surf2patch(X,Y,Z); - -%Close section if required -if closeSection==1 - I=[(1:size(Z,1)-1)' (1:size(Z,1)-1)' (2:size(Z,1))' (2:size(Z,1))' ]; - J=[size(Z,2).*ones(size(Z,1)-1,1) ones(size(Z,1)-1,1) ones(size(Z,1)-1,1) size(Z,2).*ones(size(Z,1)-1,1)]; - F_sub=sub2ind(size(Z),I,J); - F=[F;F_sub]; -end - -end +clear; close all; clc; + +%% Description +% This code processes skin and bone surface meshes for the lower limb to +% simulate amputation. Amputation is "simulated" in the sense that a 3D +% surface mesh is obtained with a shape that may be expected from a typical +% amputation. The current code is for transtibial amputation. +% +% This code requires the GIBBON MATLAB toolbox +% + +%% Control parameters + +projectFolder = fileparts(fileparts(mfilename('fullpath'))); +loadFolder=fullfile(projectFolder,'data','BodyParts3D','post'); +saveFolder=loadFolder; +saveNameGeom='BodyParts3D_right_leg_transtibial_amp'; + +fileNames={'right_femur','right_tibia','right_fibula','right_patella','right_leg_skin'}; + +saveOn=0; + +amputationDistanceTibia=180; +amputationDistanceFibula=amputationDistanceTibia-30; +amputationDistanceSkin=amputationDistanceTibia-25; +distalExcess=25; +boneRoundFactor=1; +radiusEnd=distalExcess+(amputationDistanceTibia-amputationDistanceSkin); +taperHeigth=100; +taperThreshold=15; %Regions more distance than threshold are tappered +taperFraction=0.3; +amputationDistances=[amputationDistanceTibia amputationDistanceFibula amputationDistanceSkin]; +tibiaDistalCutRadialFactor=0.5; +numBezierPoints=50; +nCut=vecnormalize([0 0 1]); %Normal direction to plane +topCropOffset=10; + +%% Loading surfaces into cell array + +%Allocate cell arrays +FT=cell(1,numel(fileNames)); +VT=cell(1,numel(fileNames)); +CT=cell(1,numel(fileNames)); + +%Loop over all surfaces +for q=1:1:numel(fileNames) + + %Import mesh + fileName=fileNames{q}; + fileName_mat=fullfile(loadFolder,[fileName,'.mat']); + model=load(fileName_mat); + F=model.faces; %Faces + V=model.vertices; %Vertices + C=q*ones(size(F,1),1); %Color label + + %Store mesh data in cell array + FT{q}=F; + VT{q}=V; + CT{q}=C; +end + +%% Compute landmarks + +P_patella_centroid=triSurfCentroid(FT{4},VT{4}); %Centroid of the patella + +%% +% Visualize surface and landmarks + +cFigure; hold on; +gpatch(FT,VT,CT,'none',0.5); +plotV(P_patella_centroid,'k.','MarkerSize',35); +axisGeom; camlight headlight; +colormap gjet; icolorbar; +gdrawnow; + +%% Cut surfaces + +%Allocate cell arrays +FT_amp=FT; +VT_amp=VT; +CT_amp=CT; +c=[2 3 5]; %Indices (=labels) for surfaces to cut +for q=1:numel(c) + cutLevelNow=P_patella_centroid(:,3)-amputationDistances(q); + + %Get surface + F=FT_amp{c(q)}; + V=VT_amp{c(q)}; + + %Use triSurfSlice to process cut + snapTolerance=mean(patchEdgeLengths(F,V))/100; + P=[0 0 cutLevelNow]; %Point on plane + [Fc,Vc,~,logicSide]=triSurfSlice(F,V,[],P,nCut,snapTolerance); + [Fc,Vc]=patchCleanUnused(Fc(~logicSide,:),Vc); %Clean-up mesh + Cc=c(q)*ones(size(Fc,1),1); %Adjusted color data (shorter list after cut) + + %Store processed mesh data in cell arrays + FT_amp{c(q)}=Fc; + VT_amp{c(q)}=Vc; + CT_amp{c(q)}=Cc; +end + +%% +% Visualization + +cFigure; hold on; +title('Cut features'); +% gpatch(FT,VT,'w','none',0.1); +gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); +axisGeom; camlight headlight; +colormap gjet; icolorbar; +gdrawnow; + +%% Set up skin taper parameterisation + +%Get tibia end curve +Ft=FT_amp{2}; +Vt=VT_amp{2}; +Ebt=patchBoundary(Ft,Vt); +indBt=edgeListToCurve(Ebt); +indBt=indBt(1:end-1); +P_tibia_end_centroid=mean(Vt(indBt,:),1); + +Fs=FT_amp{5}; +Vs=VT_amp{5}; +D=Vs(:,3)-(min(Vs(:,3))+taperHeigth); +D(D>0)=0; +D=abs(D); +D=D./max(D); +D=D.^2; + +[Dp,indMin]=minDist(Vs(:,[1 2]),Vt(indBt,:)); +Dp(Dp<=taperThreshold)=taperThreshold; +Dp=Dp-min(Dp); +Dp=Dp./max(Dp); + +R=Vt(indBt(indMin),[1 2])-Vs(:,[1 2]); +R=R.*D.*taperFraction.*Dp; + +%% Process skin taper + +Vs1=Vs; %Initialize as original +Vs(:,[1 2])=Vs(:,[1 2])+R; %Push to produce taper + +%Overide skin surface with tapered surface +FT_amp{5}=Fs; +VT_amp{5}=Vs; + +%% +% Visualisation + +% MR=sqrt(sum(R.^2,2)); + +cFigure; hold on; +title('taper morphing'); +gpatch(FT_amp,VT_amp,'w','none',0.5); +% gpatch(Fs,Vs,'w','none',0.5); +quiverVec(Vs1,R); +colormap gjet; colorbar; +axisGeom; camlight headlight; +gdrawnow; + +%% Process skin distal end closure + +P_end=P_tibia_end_centroid-[0 0 distalExcess]; +Ebs=patchBoundary(Fs,Vs); + +[Fsb,Vsb,indBs,Nd,P1,P2,P3,XB,YB,ZB]=smoothBezierClose(Fs,Vs,Ebs,P_end,numBezierPoints); + +%% +% Visualization + +cFigure; hold on; +title('Distal end closure'); +gpatch(FT_amp,VT_amp,'w','none',1); +gpatch(Fs,Vs,'w','none',0.25); +hp=gpatch(Fsb,Vsb,Vsb(:,3),'none',1); hp.FaceColor='interp'; +% gpatch(Ebs,Vs,'none','b',1,3); + +quiverVec(Vs(indBs,:),Nd(indBs,:),radiusEnd/4,'k'); + +plotV(P1,'r.-','MarkerSize',15,'LineWidth',2); +plotV(P2,'g.-','MarkerSize',15,'LineWidth',2); +plotV(P3,'b.-','MarkerSize',15,'LineWidth',2); +for q=1:size(XB,2) + plotV([XB(:,q) YB(:,q) ZB(:,q)],'k.-','LineWidth',0.5,'MarkerSize',5); +end + +plotV(P_end,'k.','MarkerSize',25); + +axisGeom; camlight headlight; +gdrawnow; + +%% +% % Visualization +% +% cFigure; hold on; +% title('Distal end closure'); +% +% hp=gpatch(Fsb,Vsb,Vsb(:,3),'k',1); hp.FaceColor='interp'; +% % gpatch(Ebs,Vs,'none','b',1,3); +% +% axisGeom; camlight headlight; +% gdrawnow; +% +% dfasfa + +%% Merge skin components and remesh + +pointSpacing=mean(patchEdgeLengths(Fs,Vs)); +[Fs,Vs]=joinElementSets({Fs,Fsb},{Vs,Vsb}); +[Fs,Vs]=mergeVertices(Fs,Vs); + +optionStructRemesh.pointSpacing=pointSpacing; %Set desired point spacing +optionStructRemesh.disp_on=0; % Turn off command window text display +[Fs,Vs]=ggremesh(Fs,Vs,optionStructRemesh); + +%% + +FT_amp{5}=Fs; +VT_amp{5}=Vs; +CT_amp{5}=5*ones(size(Fs,1),1); + +%% + +% [FT,VT,CT]=joinElementSets(FT,VT); +% [FT_amp,VT_amp,CT_amp]=joinElementSets(FT_amp,VT_amp); + +%% + +cFigure; +subplot(1,2,1); hold on; +% gpatch(FT,VT,'w','none',0.25); +gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); +axisGeom; camlight headlight; +colormap gjet; icolorbar; + +subplot(1,2,2); hold on; +gpatch(FT_amp,VT_amp,'w','k',1); +axisGeom; camlight headlight; +gdrawnow; + +%% + +Ft=FT_amp{2}; +Vt=VT_amp{2}; + +Ebs=patchBoundary(Ft,Vt); +indBs=edgeListToCurve(Ebs); +indBs=indBs(1:end-1); + +P_mid=mean(Vt(indBs,:),1); + +radiusEnd=boneRoundFactor*mean(minDist(Vt(indBs,:),P_mid))/2; +P_end=P_mid-[0 0 radiusEnd]; + +[Fsb,Vsb,indBs,Nd,P1,P2,P3,XB,YB,ZB]=smoothBezierClose(Ft,Vt,Ebs,P_end,numBezierPoints); + +%% +% Visualization +cFigure; hold on; +gpatch(Ft,Vt,'w','none',0.25); +gpatch(Ebs,Vt,'none','b',1,3); +hp=gpatch(Fsb,Vsb,'w','k',1,1); +% hp.FaceColor='interp'; +plotV(P_tibia_end_centroid,'k.','MarkerSize',35); + +quiverVec(Vt(indBs,:),Nd(indBs,:),radiusEnd/4,'k'); + +plotV(P1,'r.-','MarkerSize',15,'LineWidth',2); +plotV(P2,'g.-','MarkerSize',15,'LineWidth',2); +plotV(P3,'b.-','MarkerSize',15,'LineWidth',2); +for q=1:size(XB,2) + plotV([XB(:,q) YB(:,q) ZB(:,q)],'k.-','LineWidth',0.5,'MarkerSize',5); +end + +axisGeom; camlight headlight; +gdrawnow; + +%% Merge bone components and remesh + +pointSpacing=mean(patchEdgeLengths(Ft,Vt)); +[Ft,Vt]=joinElementSets({Ft,Fsb},{Vt,Vsb}); +[Ft,Vt]=mergeVertices(Ft,Vt); + +optionStructRemesh.pointSpacing=pointSpacing; %Set desired point spacing +optionStructRemesh.disp_on=0; % Turn off command window text display +[Ft,Vt]=ggremesh(Ft,Vt,optionStructRemesh); + +%% + +FT_amp{2}=Ft; +VT_amp{2}=Vt; +CT_amp{2}=2*ones(size(Ft,1),1); + +%% + +cFigure; +subplot(1,2,1); hold on; +% gpatch(FT,VT,'w','none',0.25); +gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); +axisGeom; camlight headlight; +colormap gjet; icolorbar; + +subplot(1,2,2); hold on; +gpatch(FT_amp,VT_amp,'w','none',0.5); +gpatch(Ft,Vt,'w','k',1); +axisGeom; camlight headlight; +gdrawnow; + +%% FIBULA + +Ft=FT_amp{3}; +Vt=VT_amp{3}; + +Ebs=patchBoundary(Ft,Vt); +indBs=edgeListToCurve(Ebs); +indBs=indBs(1:end-1); + +P_mid=mean(Vt(indBs,:),1); + +radiusEnd=boneRoundFactor*mean(minDist(Vt(indBs,:),P_mid))/2; +P_end=P_mid-[0 0 radiusEnd]; + +[Fsb,Vsb,indBs,Nd,P1,P2,P3,XB,YB,ZB]=smoothBezierClose(Ft,Vt,Ebs,P_end,numBezierPoints); + +%% +% Visualization +cFigure; hold on; +gpatch(Ft,Vt,'w','none',0.25); +gpatch(Ebs,Vt,'none','b',1,3); +hp=gpatch(Fsb,Vsb,'w','k',1,1); +% hp.FaceColor='interp'; +plotV(P_tibia_end_centroid,'k.','MarkerSize',35); + +quiverVec(Vt(indBs,:),Nd(indBs,:),radiusEnd/4,'k'); + +plotV(P1,'r.-','MarkerSize',15,'LineWidth',2); +plotV(P2,'g.-','MarkerSize',15,'LineWidth',2); +plotV(P3,'b.-','MarkerSize',15,'LineWidth',2); +for q=1:size(XB,2) + plotV([XB(:,q) YB(:,q) ZB(:,q)],'k.-','LineWidth',0.5,'MarkerSize',5); +end + +axisGeom; camlight headlight; +gdrawnow; + +%% Merge bone components and remesh + +pointSpacing=mean(patchEdgeLengths(Ft,Vt)); +[Ft,Vt]=joinElementSets({Ft,Fsb},{Vt,Vsb}); +[Ft,Vt]=mergeVertices(Ft,Vt); + +optionStructRemesh.pointSpacing=pointSpacing; %Set desired point spacing +optionStructRemesh.disp_on=0; % Turn off command window text display +[Ft,Vt]=ggremesh(Ft,Vt,optionStructRemesh); + +%% + +FT_amp{3}=Ft; +VT_amp{3}=Vt; +CT_amp{3}=3*ones(size(Ft,1),1); + +%% + +cFigure; +subplot(1,2,1); hold on; +% gpatch(FT,VT,'w','none',0.25); +gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); +axisGeom; camlight headlight; +colormap gjet; icolorbar; + +subplot(1,2,2); hold on; +gpatch(FT_amp,VT_amp,'w','none',0.5); +gpatch(Ft,Vt,'w','k',1); +axisGeom; camlight headlight; +gdrawnow; + +%% + +%Get surface +F=FT_amp{5}; +V=VT_amp{5}; +cutLevelNow=max(V(:,3))-topCropOffset; + +%Use triSurfSlice to process cut +snapTolerance=mean(patchEdgeLengths(F,V))/100; + +P=[0 0 cutLevelNow]; %Point on plane +[Fc,Vc,~,logicSide]=triSurfSlice(F,V,[],P,-nCut,snapTolerance); +[Fc,Vc]=patchCleanUnused(Fc(~logicSide,:),Vc); %Clean-up mesh + +% Remesh using ggremesh +optionStruct3.pointSpacing=mean(patchEdgeLengths(F,V)); +optionStruct3.disp_on=0; % Turn off command window text display +optionStruct3.pre.max_hole_area=100; %Max hole area for pre-processing step +optionStruct3.pre.max_hole_edges=0; %Max number of hole edges for pre-processing step + +[Fc,Vc]=ggremesh(Fc,Vc,optionStruct3); + +Ebs=patchBoundary(Fc,Vc); +indBs=edgeListToCurve(Ebs); +indBs=indBs(1:end-1); + + + +Fcs=Fc; +Vcs=Vc; +indBss=indBs; + +%% + +%Get surface +F=FT_amp{1}; +V=VT_amp{1}; + +%Use triSurfSlice to process cut +snapTolerance=mean(patchEdgeLengths(F,V))/100; +P=[0 0 cutLevelNow]; %Point on plane +[Fc,Vc,~,logicSide]=triSurfSlice(F,V,[],P,-nCut,snapTolerance); +[Fc,Vc]=patchCleanUnused(Fc(~logicSide,:),Vc); %Clean-up mesh + +% Remesh using ggremesh +optionStruct3.pointSpacing=mean(patchEdgeLengths(F,V)); +optionStruct3.disp_on=0; % Turn off command window text display +optionStruct3.pre.max_hole_area=100; %Max hole area for pre-processing step +optionStruct3.pre.max_hole_edges=0; %Max number of hole edges for pre-processing step + +[Fc,Vc]=ggremesh(Fc,Vc,optionStruct3); + +Ebs=patchBoundary(Fc,Vc); +indBs=edgeListToCurve(Ebs); +indBs=indBs(1:end-1); + +%% + +z=0.5*(mean(Vc(indBs,3))+mean(Vc(indBs,3))); + +Vcs(indBss,3)=z; +Vc(indBs,3)=z; + +%% + +pointSpacing=mean(patchEdgeLengths(Fcs,Vcs)); +[Ftt,Vtt]=regionTriMesh3D({Vcs(indBss,:),Vc(indBs,:)},pointSpacing,0,'linear'); + +%% + +cFigure; +gpatch(Fcs,Vcs,'gw','k',1); +gpatch(Fc,Vc,'rw','k',1); +gpatch(Ftt,Vtt,'bw','k',0.5); +colormap gjet; +axisGeom; camlight headlight; +gdrawnow; + +%% +FT_amp{1}=fliplr(Fc); %invert femur normals +VT_amp{1}=Vc; +CT_amp{1}=1*ones(size(Fc,1),1); + +FT_amp{5}=Fcs; +VT_amp{5}=Vcs; +CT_amp{5}=5*ones(size(Fcs,1),1); + +FT_amp{6}=Ftt; +VT_amp{6}=Vtt; +CT_amp{6}=6*ones(size(Ftt,1),1); + +%% + +cFigure; hold on; +gpatch(FT_amp,VT_amp,CT_amp,'none',0.5); +% patchNormPlot(FT_amp,VT_amp); +axisGeom; camlight headlight; +colormap gjet; icolorbar; +gdrawnow; + +%% Save model + +if saveOn==1 + saveName_mat=fullfile(saveFolder,[saveNameGeom,'.mat']); + save(saveName_mat,'FT_amp','VT_amp','CT_amp'); +end + +%% + +function [Fsb,Vsb,indBs,Nd,P1,P2,P3,XB,YB,ZB]=smoothBezierClose(Fs,Vs,Ebs,P_end,numPoints) + +indBs=edgeListToCurve(Ebs); +indBs=indBs(1:end-1); + +[~,~,Ns]=patchNormal(Fs,Vs); + +[~,~,NEb]=edgeVec(Ebs,Vs); +NEb=vecnormalize(NEb); +Nd=vecnormalize(cross(NEb,Ns)); + +cPar.n=25; +Nd=patchSmooth(Ebs,Nd,[],cPar); + +%Bezier point set 1 +P1=Vs(indBs,:); +P1_mid=mean(P1,1); + +distanceEnd=sqrt(sum((P_end-P1_mid).^2,2)); + +%Bezier point set 2 +f=1./abs(Nd(indBs,3)); %Extend factor for direction vectors +P2=P1+distanceEnd/2*f(:,ones(1,3)).*Nd(indBs,:); %Position half-way + +%Bezier point set 3 +P3=P2; +P3(:,3)=P_end(:,3); %Shift to bottom +P3=P3-mean(P3,1); %Shift on own centre +[t,r]=cart2pol(P3(:,1),P3(:,2)); %Polar coordinates +[P3(:,1),P3(:,2)]=pol2cart(t,distanceEnd/2*ones(size(r))); %Force constant radius +P3=P3+P_end; %Place centered on end + +%Bezier point set 4 +P4=P_end; + +XB=zeros(numPoints,size(P1,1)); +YB=zeros(numPoints,size(P1,1)); +ZB=zeros(numPoints,size(P1,1)); +for q=1:1:size(P1,1) + p=[P1(q,:); P2(q,:); P3(q,:); P4]; %Control points + V_bezier=bezierCurve(p,numPoints); %Compute bezier curve + XB(:,q)=V_bezier(:,1); + YB(:,q)=V_bezier(:,2); + ZB(:,q)=V_bezier(:,3); +end + +[Fsb,Vsb]=meshToPatch(XB,YB,ZB,1); +[Fsb,Vsb]=quad2tri(Fsb,Vsb); +[Fsb,Vsb]=mergeVertices(Fsb,Vsb); +[Fsb,Vsb]=patchCleanUnused(Fsb,Vsb); + +end + +%% + +function [F,V]=meshToPatch(X,Y,Z,closeSection) + +%Create quad patch data +[F,V] = surf2patch(X,Y,Z); + +%Close section if required +if closeSection==1 + I=[(1:size(Z,1)-1)' (1:size(Z,1)-1)' (2:size(Z,1))' (2:size(Z,1))' ]; + J=[size(Z,2).*ones(size(Z,1)-1,1) ones(size(Z,1)-1,1) ones(size(Z,1)-1,1) size(Z,2).*ones(size(Z,1)-1,1)]; + F_sub=sub2ind(size(Z),I,J); + F=[F;F_sub]; +end + +end + +%% +% _*SimuLimb footer text*_ +% +% License: +% +% Copyright (C) 2006-2021 Kevin Mattheus Moerman and the SimuLimb contributors +% +% Licensed under the Apache License, Version 2.0 (the "License"); +% you may not use this file except in compliance with the License. +% You may obtain a copy of the License at +% +% http://www.apache.org/licenses/LICENSE-2.0 +% +% Unless required by applicable law or agreed to in writing, software +% distributed under the License is distributed on an "AS IS" BASIS, +% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +% See the License for the specific language governing permissions and +% limitations under the License.