diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 00000000..f7dbcad5 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "STAPLE/external/PelvicLandmarkIdentification"] + path = STAPLE/external/PelvicLandmarkIdentification + url = https://github.com/RWTHmediTEC/PelvicLandmarkIdentification diff --git a/Example_use_Fischer2019.m b/Example_use_Fischer2019.m new file mode 100644 index 00000000..62ffe5db --- /dev/null +++ b/Example_use_Fischer2019.m @@ -0,0 +1,109 @@ +%-------------------------------------------------------------------------% +% Copyright (c) 2024 MCM Fischer. % +% Author: Maximilian C. M. Fischer, 2024 % +% email: N/A % +% ----------------------------------------------------------------------- % +% This example demonstrates how to setup a STAPLE workflow using an % +% algorithm for the morphological analyses of the pelvis published by % +% Fischer et al. https://doi.org/10.1038/s41598-019-49573-4 % +% ----------------------------------------------------------------------- % +clear; clc; close all +addpath(genpath('./STAPLE')); + +%----------% +% SETTINGS % +%----------% +output_models_folder = 'opensim_models_examples'; + +% folder where the various datasets (and their geometries) are located. +datasets_folder = 'bone_datasets'; + +% datasets that you would like to process +dataset_set = {'JIA_MRI'}; + +% anthropometrics +% H = ?; %cm +subj_mass = 76.5; %kg + +% Cell array with the bone geometries that you would like to process +% NOTE: The complete pelvis mesh including both hip bones and the sacrum is +% required for the algorithm of [Fischer 2019]. Consider the requirements +% for the mesh of the pelvis described in the document header of the +% fuction pelvicLandmarkID. +bones_list = {'pelvis','femur_r'}; + +% visualization geometry format (options: 'stl' or 'obj') +vis_geom_format = 'obj'; + +% choose the definition of the joint coordinate systems (see documentation) +% options: 'Modenese2018' or 'auto2020' +workflow = 'auto2020'; +%-------------------------------------- + +tic + +% create model folder if required +if ~isfolder(output_models_folder); mkdir(output_models_folder); end + +for n_d = 1:numel(dataset_set) + + % current dataset being processed + cur_dataset = dataset_set{n_d}; + + % folder from which triangulations will be read + tri_folder = fullfile(datasets_folder, cur_dataset,'tri'); + + % create geometry set structure for all 3D bone geometries in the dataset + triGeom_set = createTriGeomSet(bones_list, tri_folder); + + % get the body side (can also be specified by user as input to funcs) + side = inferBodySideFromAnatomicStruct(triGeom_set); + + % model and model file naming + model_name = ['auto_',dataset_set{n_d},'_',upper(side)]; + model_file_name = [model_name, '.osim']; + + % log printout + log_file = fullfile(output_models_folder, [model_name, '.log']); + logConsolePrintout('on', log_file); + + % create bone geometry folder for visualization + geometry_folder_name = [model_name, '_Geometry']; + geometry_folder_path = fullfile(output_models_folder,geometry_folder_name); + + % convert geometries in chosen format (30% of faces for faster visualization) + writeModelGeometriesFolder(triGeom_set, geometry_folder_path, vis_geom_format,0.3); + + % initialize OpenSim model + osimModel = initializeOpenSimModel(model_name); + + % create bodies + osimModel = addBodiesFromTriGeomBoneSet(osimModel, triGeom_set, geometry_folder_name, vis_geom_format); + + % process bone geometries (compute joint parameters and identify markers) + % pelvis femur tibia + [JCS, BL, CS] = processTriGeomBoneSet(triGeom_set, side, 'Fischer2019','Kai2014','Kai2014'); + + % create joints + createOpenSimModelJoints(osimModel, JCS, workflow); + + % update mass properties to those estimated using a scale version of + % gait2392 with COM based on Winters's book. + osimModel = assignMassPropsToSegments(osimModel, JCS, subj_mass); + + % add markers to the bones + addBoneLandmarksAsMarkers(osimModel, BL); + + % finalize connections + osimModel.finalizeConnections(); + + % print + osimModel.print(fullfile(output_models_folder, model_file_name)); + + % inform the user about time employed to create the model + disp('-------------------------') + disp(['Model generated in ', num2str(toc)]); + disp(['Saved as ', fullfile(output_models_folder, model_file_name),'.']); + disp('-------------------------') + logConsolePrintout('off'); +end \ No newline at end of file diff --git a/README.md b/README.md index 1649c0d3..377d7cac 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,8 @@ # STAPLE: Shared Tools for Automatic Personalised Lower Extremity modelling -[![License: CC BY-NC 4.0](https://img.shields.io/badge/License-CC%20BY--NC%204.0-lightgrey.svg)](https://creativecommons.org/licenses/by-nc/4.0/) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4275285.svg)](https://doi.org/10.5281/zenodo.4275285) ![visitors](https://visitor-badge.glitch.me/badge?page_id=modenaxe.msk-STAPLE) +[![License: CC BY-NC 4.0](https://img.shields.io/badge/License-CC%20BY--NC%204.0-lightgrey.svg)](https://creativecommons.org/licenses/by-nc/4.0/) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4428103.svg)](https://doi.org/10.5281/zenodo.4428103) ![visitors](https://visitor-badge.laobi.icu/badge?page_id=modenaxe.msk-STAPLE) + +### ⚠️ STAPLE is released under a [non-commercial license](#license) and it is free to use for academic purposes only. For any other use please contact the authors. # Table of contents - [What is STAPLE?](#what-is-staple) @@ -14,7 +16,7 @@ - [Detailed steps to setup a STAPLE workflow](#detailed-steps-to-setup-a-staple-workflow) - [Algorithms for bone morphological analysis](#algorithms-for-bone-morphological-analysis) - [Datasets available for testing](#datasets-available-for-testing) - - [Further notes on STAPLE](#further-notes-on-staple) + - [STAPLE variables and conventions](#STAPLE-variables-and-conventions) - [Troubleshooting your workflow](#troubleshooting-your-workflow) - [Does STAPLE work only with OpenSim?](#does-staple-work-only-with-opensim) - [What are the differences between the STAPLE toolbox, NMSBuilder and the MAP Client for generating OpenSim models?](#what-are-the-differences-between-the-staple-toolbox-nmsbuilder-and-the-map-client-for-generating-opensim-models) @@ -62,6 +64,10 @@ For example, models of hip, knee and ankle joints can be created as individual m ![articular_surfaces](./images/artic_surfaces.png) +* **Transformations to local reference systems**: STAPLE performs morphological analyses using the medical images reference system where the bone geometries are segmented, but then it can transform both the bone geometries and model objects to standard segment-based biomechanical reference systems. This is an important feature if you are using STAPLE in OpenSim, as joint reactions are expressed in body reference systems. See the provided [example of this functionality](./Example_BCS_leg_right.m). + +![reference_systems](./images/reference_systems.png) + * **Basic identification of bony landmarks**: certain bony landmarks can be easily identified following the morphological analysis of the bone surfaces. These landmarks are intended as first guess for registration with gait analysis data. * **Merging subject-specific and generic models**: STAPLE includes some basic utilities to merge partial subject-specific skeletal models (obtainable from localized MRI scans) with generic musculoskeletal models. See the provided [advanced examples](./advanced_examples). @@ -286,11 +292,11 @@ The BCS is a MATLAB structure with the following fields: When running a morphological analysis, the joint coordinate system (JCS) associated with the analysed geometry is also returned, together with the BCS. The JCS is a MATLAB structure with the following fields: -1. ****: +1. **joint name**: the name of the joint if more than one are stored in the variable. 2. **Origin** [3x1] vector: the origin of the joint in the bone being analysed 3. **V** [3x3] matrix: the pose matrix of the joint 4. **parent_location** [1x3] vector: the origin of the parent coordinate system in OpenSim -5 **parent_orientation** [1x3] vector: three angles describing the pose of the parent reference system in OpenSim, obtained by matrix decomposition with XYZ rotation order. +5. **parent_orientation** [1x3] vector: three angles describing the pose of the parent reference system in OpenSim, obtained by matrix decomposition with XYZ rotation order. 6. **child_location** [1x3] vector: the origin of the child coordinate system in OpenSim 7. **child_orientation** [1x3] vector: three angles describing the pose of the parent reference system in OpenSim, obtained by matrix decomposition with XYZ rotation order. @@ -315,6 +321,14 @@ Note that the algorithms used on a single bone normally cannot define a joint co Bone landmarks are stored in MATLAB structures as [3x1] vectors. Each field name of the structure corresponds to the name of a landmark. +| Bone geometry | Landmarks | +| --- | --- | +| pelvis | | +| femur | | +| tibia | | +| patella | RLOW/LLOW: most distal point of patella | +| talus | None | +| foot | | ## Troubleshooting your workflow @@ -323,9 +337,11 @@ Before informing us as suggested in [the contributing guidelines](/CONTRIBUTING. - [ ] ensure that the entire `STAPLE` folder is on your MATLAB path - [ ] if you are using one of the workflows from the examples, ensure that you are using the correct names for the bone geometries, e.g. `pelvis_no_sacrum`, `femur_r`, `tibia_r`, etc. - [ ] ensure that the quality of your bone surface geometries is sufficient for running your selected algorithm. The GIBOC algorithms, in particular, require relatively good quality surface meshes. You can have an idea of what we mean by "good mesh" consulting the table that describes the datasets provided with STAPLE. If necessary, use the filters and tools available on software like [MeshLab](https://www.meshlab.net/) to improve your dataset. +- [ ] processing the meshes so that they are triangular and contain a reasonable number of vertices (decimation can be applied) also helps with the OpenSim visualization. See [this related tweet](https://twitter.com/JohnRHutchinson/status/1455490276298997766). - [ ] if possible, try running alternative algorithms in order to establish if your bone mesh is processable in the first place or if you have encounter a proper bug in the algorithm. You can specify the processing algorithms for each bone as input to the `processTriGeomBoneSet.m` function. For bad quality meshes, we recommend using the `STAPLE` algorithm at the pelvis and `Kai2014` algorithms for femur and tibia. Keep an eye on issue #78 as we will develop an example on how to do this. - [ ] verify that processing of your dataset is not failing because of the [current limitations](#current-limitations) of the STAPLE toolbox. + ## Does STAPLE work only with OpenSim? The algorithms collected in the STAPLE toolbox were proposed in publications that did not have modelling focus, and can be applied in broader contexts, such as reference system definition for orthopaedics applications or modelling in other platforms. Similarly, the outputs of a STAPLE workflow, e.g. from `processTriGeomBoneSet.m`, include all the necessary information to create a kinematic or kinetic model in any biomechanical modelling workflow. All our modelling examples, however, rely on OpenSim. diff --git a/STAPLE/GIBOC-core/SubFunctions/FittingFun/LSSLFitPatellaRidge.m b/STAPLE/GIBOC-core/SubFunctions/FittingFun/LSSLFitPatellaRidge.m index 5eae29b9..d5a0aef5 100644 --- a/STAPLE/GIBOC-core/SubFunctions/FittingFun/LSSLFitPatellaRidge.m +++ b/STAPLE/GIBOC-core/SubFunctions/FittingFun/LSSLFitPatellaRidge.m @@ -1,72 +1,86 @@ -function [ U, Uridge ,LowestPoints_End ] = LSSLFitPatellaRidge( TR,U,nbSlice,StartDist,EndDist, debug_plot) %Least Square Straight Line Fit on Patellar Ridge %To obtain the direction of the ridge a Least Square Straight Line is %fitted on the lowest points on slices of normal U, and the normal U is %updated until convergence (of U or number of iterations > 100) +% U is the vector ~aligned with the ridge +%-------------------------------------------------------------------------% +% Author: Luca Modenese & Jean-Baptiste Renault. +% Copyright 2020 Luca Modenese & Jean-Baptiste Renault +%-------------------------------------------------------------------------% +function [U, Uridge, LowestPoints] = LSSLFitPatellaRidge( patellaTri,... + U,... + nbSlice,... + startDist,... + endDist,... + debug_plots) +% initialization +if nargin<3; nbSlice = 50; end +if nargin<4; startDist = 0.025*range(patellaTri.Points*U); end %Offset distance from first points +if nargin<5; endDist = startDist; end%Offset distance from last points +if nargin<6; debug_plots=0; end -%% inputs Tests -if nargin<3 - StartDist = 0.025*range(TR.Points*U); %Offset distance from first points - EndDist = StartDist; %Offset distance from last points - nbSlice = 50; -end - -if nargin<4 - StartDist = 0.025*range(TR.Points*U); %Offset distance from first points - EndDist = StartDist; %Offset distance from last points -end - -if nargin<6 - debug_plot=0; -end -%% Code +% initialize fitting U0 = U; -reslts = [0;0;0]; -k=0; +U_old = [0;0;0]; +k=0; test = true; while test && k < 100 % Cut at 100 iterations (convergence generally occuring before 10) k=k+1; - Alt = linspace( min(TR.Points*U)+StartDist , max(TR.Points*U)-EndDist, nbSlice); + % points along U where the slices will be performed + Alt = linspace( min(patellaTri.Points*U)+startDist , max(patellaTri.Points*U)-endDist, nbSlice); LowestPoints = zeros(length(Alt),3); + % slice patella along U i=0; for d = -Alt i=i+1; - [ Curves , ~ , ~ ] = TriPlanIntersect( TR, U , d ); + [ Curves , ~ , ~ ] = TriPlanIntersect( patellaTri, U , d ); EdgePts = vertcat(Curves(:).Pts); [~,lowestPointID] = min(EdgePts(:,3)); LowestPoints(i,:) = EdgePts(lowestPointID(1),:); + if debug_plots == 1 + title('Slicing and LowerPoints on AS') + quickPlotTriang(patellaTri); + plot3(Curves(:).Pts(:,1), Curves(:).Pts(:,2), Curves(:).Pts(:,3),'k'); hold on + plot3(LowestPoints(:,1), LowestPoints(:,2), LowestPoints(:,3),'r'); hold on + axis equal + end end + % updating direction of U usign the variance of the patellar ridge path [Vridge_all,~] = eig(cov(LowestPoints)); - U = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3); - U(3) = 0; U = normalizeV( U ); - reslts(1:3,end+1) = U; - test = reslts(1:3,end-1)'*reslts(1:3,end) > (1 - 100*eps); -end - - -if nargout == 2 -Uridge = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3); - -elseif nargout == 3 - Uridge = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3); - LowestPoints_End = LowestPoints; + % vector describing the direction of the ridge + Uridge = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3); + + % the vector will be in the XY plane where the patella is a "circle" + % third component = 0 +% U = Uridge; +% U(3) = 0; + U = normalizeV([Uridge(1:2); 0]); + + % how much the vector is changing in direction + U_old(1:3,end+1) = U; + test = U_old(1:3,end-1)'*U_old(1:3,end) > (1 - 100*eps); end -%% Plots -if debug_plot == 1 +% debug plots +if debug_plots == 1 figure() - trisurf(TR,'facecolor','red','faceAlpha',0.5) - hold on - grid off - axis equal + title('Final vector of ridge') + trisurf(patellaTri,'facecolor','red','faceAlpha',0.5) + hold on; grid off; axis equal; + LocalCS = struct(); + LocalCS.Origin = [0 0 0]'; + LocalCS.V = eye(3); + quickPlotRefSystem(LocalCS); + + % plot trace of points on the proximal surface pl3t(LowestPoints,'gs') plotArrow( U, 1.5, mean(LowestPoints), 30, 1, 'k') - Uridge = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3); plotArrow( Uridge, 1.5, mean(LowestPoints), 50, 1, 'c') hold off end + end diff --git a/STAPLE/GIBOC-core/SubFunctions/FittingFun/fitQuadriTalus.m b/STAPLE/GIBOC-core/SubFunctions/FittingFun/fitQuadriTalus.m index a6676d3c..c8ff8590 100644 --- a/STAPLE/GIBOC-core/SubFunctions/FittingFun/fitQuadriTalus.m +++ b/STAPLE/GIBOC-core/SubFunctions/FittingFun/fitQuadriTalus.m @@ -17,78 +17,105 @@ plotOn = 0; end +% First project the mesh on the plane perpendicular to X0 (1st inertia axis) +Pts_inertia = (V_all'*Tr.Points')'; -%% First project the mesh on the plan perpendicular to X0 (1st inertia axis) -YZ0 = [V_all(:,2),V_all(:,3)]; -PX0 = YZ0*inv(YZ0'*YZ0)*YZ0'; %Projection matrix -Pts_proj = V_all'*PX0*Tr.Points'; +if plotOn + figure('color','w','numbertitle','off', ... + 'name', ['Debug Figure: ' mfilename '.m: Projection on X0 (1st inertia axis)']); + axis equal tight; hold on + plot3(Pts_inertia(:,1),Pts_inertia(:,2),Pts_inertia(:,3),'g.') + xlabel('X'); ylabel('Y'); zlabel('Z') + view(90,0) +end -% Keep only the 2 and 3rd coordinates -Pts_proj_2D = Pts_proj(2:3,:)'; +% Keep only the 2nd (Y0) and 3rd (Z0) dimension +Pts_proj_2D = Pts_inertia(:,2:3); +assert(sum(sum(Pts_inertia(:,2:3)-Pts_proj_2D)) < 1e-8) %Get the convex hull of the point cloud K = convhull(Pts_proj_2D,'simplify', false); -% Find the max edge -maxDist = -1; -for i=1:length(K) - for j=1:length(K) - diff1 = Pts_proj_2D(K(j),1)-Pts_proj_2D(K(i),1); - diff2 = Pts_proj_2D(K(j),2)-Pts_proj_2D(K(i),2); - dist = sqrt(diff1^2+diff2^2); - if dist>maxDist - maxDist=dist; - pointPair1=[K(i),K(j)]; - U1=[diff1;diff2]; - end - end +if plotOn + projFig = figure('color','w','numbertitle','off', ... + 'name', ['Debug Figure: ' mfilename '.m: Convex hull of the projection on X0']); + projAxes = axes(); + axis equal tight; hold on + plot(Pts_proj_2D(:,1),Pts_proj_2D(:,2),'b.') + plot(Pts_proj_2D(K,1),Pts_proj_2D(K,2),'k-') + xlabel('Y'); ylabel('Z') +end + +% Find the two points of the convex hull with max. distance to each other +% assuming that it is the 1st diagonal of the quadrilateral +pwDist = squareform(pdist(Pts_proj_2D(K,:))); +[~, maxDistIdx] = max(pwDist(:)); +[rowPP1,colPP1] = ind2sub(size(pwDist),maxDistIdx); +pointPair1 = [K(rowPP1),K(colPP1)]; +U1 = (Pts_proj_2D(pointPair1(2),:)-Pts_proj_2D(pointPair1(1),:))'; + +if plotOn + maxDist1D_H = plot(Pts_proj_2D(pointPair1,1),Pts_proj_2D(pointPair1,2),'r-o'); + legend(maxDist1D_H, '1st diagonal') end % Hacky way to find the second diagonal of the quadrilateral by shrinking -% along the point cloud in the direction (U1) of the first found diagonal +% along the point cloud in the direction (U1) of the 1st diagonal found U1 = normalizeV(U1); Ru = [U1(1),-U1(2);U1(2),U1(1)]; -Pts_proj_2D_shr = transpose(Ru*[0.5,0;0,1]*Ru'*Pts_proj_2D'); -% plot(Pts_proj_2D_shr(:,1),Pts_proj_2D_shr(:,2),'b.') - - -% Find the max edge on the shrink points cloud -maxDist = -1; -for i=1:length(K) - for j=1:length(K) - diff1 = Pts_proj_2D_shr(K(j),1)-Pts_proj_2D_shr(K(i),1); - diff2 = Pts_proj_2D_shr(K(j),2)-Pts_proj_2D_shr(K(i),2); - dist = sqrt(diff1^2+diff2^2); - if dist>maxDist - maxDist=dist; - pointPair2=[K(i),K(j)]; - U2=[diff1;diff2]; - end - end +Pts_proj_2D_shr = (Ru*[0.5,0;0,1]*Ru'*Pts_proj_2D')'; + +if plotOn + figure('color','w','numbertitle','off', ... + 'name', ['Debug Figure: ' mfilename '.m: Projection on X0 shrinked along 1st diagonal']); + axis equal tight; hold on + plot(Pts_proj_2D_shr(:,1),Pts_proj_2D_shr(:,2),'g.') + line(Pts_proj_2D_shr(K,1),Pts_proj_2D_shr(K,2)) + xlabel('Y'); ylabel('Z') end -% Get the quadrilateral vertices -quadriV = [pointPair1(1); pointPair2(1);... - pointPair1(2); pointPair2(2);... - pointPair1(1)]; +% Find the two points of the convex hull with max. distance to each other +% assuming that it is the 2nd diagonal of the quadrilateral +pwDist = squareform(pdist(Pts_proj_2D_shr(K,:))); +[~, maxDistIdx] = max(pwDist(:)); +[rowPP2,colPP2] = ind2sub(size(pwDist),maxDistIdx); +pointPair2 = [K(rowPP2),K(colPP2)]; + +if plotOn + maxDist2D_H = plot(Pts_proj_2D_shr(pointPair2,1),Pts_proj_2D_shr(pointPair2,2),'r--o'); + legend(maxDist2D_H, '2nd diagonal') +end + +% Sort the quadrilateral vertices in counter-clockwise direction along the +% convex hull +quadriVIdx = sort([rowPP1,colPP1,rowPP2,colPP2]); +quadriV = K([quadriVIdx,quadriVIdx(1)]); -% Get the length of the quadri edges +if plotOn + maxDist2D_H = plot(projAxes, Pts_proj_2D(pointPair2,1),Pts_proj_2D(pointPair2,2),'r--o'); + qLitVerts_H = text(projAxes, Pts_proj_2D(quadriV(1:4),1),Pts_proj_2D(quadriV(1:4),2),{'1','2','3','4'},... + 'HorizontalAlignment','right', 'FontSize',14, 'FontWeight','bold'); + legend([maxDist1D_H,maxDist2D_H], {'1st diagonal','2nd diagonal'}) +end + +% Get the length of the edges of the quadrilateral edgesLength = zeros(4,1); -for i=1:4 - diff1 = Pts_proj_2D(quadriV(i+1),1)-Pts_proj_2D(quadriV(i),1); - diff2 = Pts_proj_2D(quadriV(i+1),2)-Pts_proj_2D(quadriV(i),2); - edgesLength(i) = sqrt(diff1^2+diff2^2); +for i=1:4 + dy = Pts_proj_2D(quadriV(i+1),1)-Pts_proj_2D(quadriV(i),1); + dz = Pts_proj_2D(quadriV(i+1),2)-Pts_proj_2D(quadriV(i),2); + edgesLength(i) = sqrt(dy^2+dz^2); end - [~,Imax] = max(edgesLength); -% Indices of start vertices of quadrilateral -I_V_sup = mod(Imax+2,4); -% Edge corresponding to the superior part of the bone is assumed to be the -% one opposing the largest one +% The edge corresponding to the superior part of the bone is assumed to be +% the one opposing the largest one +I_V_sup = mod(Imax+2,4); % Index of the start vertex of the superior edge +if I_V_sup == 0 + I_V_sup = 4; +end Edge_sup = quadriV(I_V_sup:I_V_sup+1); -% Get the direction of the edge : +% Get the direction of the edge U_edge_sup = [Pts_proj_2D(Edge_sup(2),1)-Pts_proj_2D(Edge_sup(1),1);... Pts_proj_2D(Edge_sup(2),2)-Pts_proj_2D(Edge_sup(1),2)]; %Z0_proj is normal to edge direction @@ -103,18 +130,14 @@ Y0 = cross(Z0,V_all(:,1)); if plotOn - figure() - plot(Pts_proj_2D(:,1),Pts_proj_2D(:,2),'b.') - hold on - plot(Pts_proj_2D(K,1),Pts_proj_2D(K,2),'k-') - axis equal - plot(Pts_proj_2D(pointPair1,1),Pts_proj_2D(pointPair1,2),'r-*') - plot(Pts_proj_2D(pointPair2,1),Pts_proj_2D(pointPair2,2),'r-*') - plot(Pts_proj_2D(quadriV,1),Pts_proj_2D(quadriV,2),'g-s',... - 'linewidth',2) - plot(Pts_proj_2D(Edge_sup,1),Pts_proj_2D(Edge_sup,2),'m-o',... - 'linewidth',3) -end - + qLitEdge_H = plot(projAxes, Pts_proj_2D(quadriV,1),Pts_proj_2D(quadriV,2),'g-s',... + 'linewidth',2); + supEdge_H = plot(projAxes, Pts_proj_2D(Edge_sup,1),Pts_proj_2D(Edge_sup,2),'m-o',... + 'linewidth',3); + legend([maxDist1D_H,maxDist2D_H,qLitEdge_H,supEdge_H],... + {'1st diagonal','2nd diagonal','Quadrilateral edges','Superior edge'}) + uistack(qLitVerts_H, 'top'); + figure(projFig) end +end \ No newline at end of file diff --git a/STAPLE/GIBOC-core/SubFunctions/GeometricFun/PtsOnCondylesFemur.m b/STAPLE/GIBOC-core/SubFunctions/GeometricFun/PtsOnCondylesFemur.m index a28554cc..a7399935 100644 --- a/STAPLE/GIBOC-core/SubFunctions/GeometricFun/PtsOnCondylesFemur.m +++ b/STAPLE/GIBOC-core/SubFunctions/GeometricFun/PtsOnCondylesFemur.m @@ -75,21 +75,24 @@ new_ver_line = R*ver_line; new_horz_line = R*horz_line; -% figure() -% plot(Pts_Epiphysis(:,3),Pts_Epiphysis(:,1),'g.') -% hold on -% axis equal -% % plot(Pts_Epiphysis(OUT_Elps,3),Pts_Epiphysis(OUT_Elps,1),'c*') -% plot( rotated_ellipse(:,1),rotated_ellipse(:,2),'r' ); -% plot(Pts_Epiphysis(I_kept,3),Pts_Epiphysis(I_kept,1),'rs') -% plot( new_ver_line(1,:),new_ver_line(2,:),'r' ); -% plot( new_horz_line(1,:),new_horz_line(2,:),'r' ); -% quiver(Elps.X0_in,Elps.Y0_in,50*cos(-Elps.phi),50*sin(-Elps.phi)); -% quiver(Elps.X0_in,Elps.Y0_in,50*sin(Elps.phi),50*cos(Elps.phi)); -% plot(mean(Pts_Epiphysis(:,3)),mean(Pts_Epiphysis(:,1)),'ks') -% plot(Pts_Epiphysis(I,3),Pts_Epiphysis(I,1),'ks') -% plot(PtsCondyle_0(K,3),PtsCondyle_0(K,1),'k-') -% plot(PtsCondyle_0(:,3),PtsCondyle_0(:,1),'kd') -% plot(PtsCondyle_0(:,3),PtsCondyle_0(:,1),'k*') +debug_plots = 0; +if debug_plots + figure() + plot(Pts_Epiphysis(:,3),Pts_Epiphysis(:,1),'g.') + hold on + axis equal + % plot(Pts_Epiphysis(OUT_Elps,3),Pts_Epiphysis(OUT_Elps,1),'c*') + plot( rotated_ellipse(:,1),rotated_ellipse(:,2),'r' ); + plot(Pts_Epiphysis(I_kept,3),Pts_Epiphysis(I_kept,1),'rs') + plot( new_ver_line(1,:),new_ver_line(2,:),'r' ); + plot( new_horz_line(1,:),new_horz_line(2,:),'r' ); + quiver(Elps.X0_in,Elps.Y0_in,50*cos(-Elps.phi),50*sin(-Elps.phi)); + quiver(Elps.X0_in,Elps.Y0_in,50*sin(Elps.phi),50*cos(Elps.phi)); + plot(mean(Pts_Epiphysis(:,3)),mean(Pts_Epiphysis(:,1)),'ks') + plot(Pts_Epiphysis(I,3),Pts_Epiphysis(I,1),'ks') + plot(PtsCondyle_0(K,3),PtsCondyle_0(K,1),'k-') + plot(PtsCondyle_0(:,3),PtsCondyle_0(:,1),'kd') + plot(PtsCondyle_0(:,3),PtsCondyle_0(:,1),'k*') +end end diff --git a/STAPLE/GIBOC-core/SubFunctions/TriangulationFun/TriFillPlanarHoles.m b/STAPLE/GIBOC-core/SubFunctions/TriangulationFun/TriFillPlanarHoles.m index ab64de28..5f32f19c 100644 --- a/STAPLE/GIBOC-core/SubFunctions/TriangulationFun/TriFillPlanarHoles.m +++ b/STAPLE/GIBOC-core/SubFunctions/TriangulationFun/TriFillPlanarHoles.m @@ -27,7 +27,7 @@ Curves(i).NodesID(j+1)=Segments(1,2); Curves(i).FB(end+1,:) = Segments(1,:); - % Remove the semgents added to Curves(i) from the segments matrix + % Remove the segments added to Curves(i) from the segments matrix Segments(1,:)=[]; j=j+1; @@ -64,9 +64,16 @@ i=i+1; end -for i = 1 : length(Curves) - Curves(i).NodesID(Curves(i).NodesID==0) = []; - Curves(i).Pts = Pts(Curves(i).NodesID,:); +if isempty(fields(Curves)) + % warning + disp('No holes on triangulation.'); + TRout = TRin; + return +else + for i = 1 : length(Curves) + Curves(i).NodesID(Curves(i).NodesID==0) = []; + Curves(i).Pts = Pts(Curves(i).NodesID,:); + end end diff --git a/STAPLE/algorithms/Fischer2019_pelvis.m b/STAPLE/algorithms/Fischer2019_pelvis.m new file mode 100644 index 00000000..05ae0ccc --- /dev/null +++ b/STAPLE/algorithms/Fischer2019_pelvis.m @@ -0,0 +1,116 @@ +function [BCS, JCS, PelvisBL] = Fischer2019_pelvis(Pelvis, side_raw, ... + result_plots, debug_plots, in_mm) +%FISCHER2019_PELVIS Calculation of the pelvic CS based on [Fischer 2019]. +% +% This is a wrapper function to integrate the MATLAB function: +% https://github.com/RWTHmediTEC/PelvicLandmarkIdentification +% +% The method is described in detail in: +% Fischer, M. C. M. et al. A robust method for automatic +% identification of landmarks on surface models of the pelvis. +% Sci. Rep. 9, 13322 (2019) +% https://doi.org/10.1038/s41598-019-49573-4 +% +% ATTENTION: Consider the requirements for the mesh of the pelvis +% described in the document header of the fuction pelvicLandmarkID. +% +% AUTHOR: Maximilian C. M. Fischer +% VERSION: 1.0.0 +% DATE: 2024-08-19 +% COPYRIGHT (C) 2024 Maximilian C. M. Fischer +% LICENSE: CC BY-NC 4.0 +% + +if nargin<2; side_raw = 'r'; end +if nargin<3; result_plots = 1; end +if nargin<4; debug_plots = 0; end +if nargin<5; in_mm = 1; end +if in_mm == 1; dim_fact = 0.001; else; dim_fact = 1; end + +% Get side identifier correspondent to the body side (used for hip joint +% parent). No need for the sign because left and right reference systems +% are identical +[~, side_low] = bodySide2Sign(side_raw); + +disp('---------------------') +disp(' FISCHER - PELVIS '); +disp('---------------------') +disp(['* Hip Joint : ', upper(side_low)]); +disp(['* Method : ', 'https://doi.org/10.1038/s41598-019-49573-4']); +disp(['* Result Plots: ', convertBoolean2OnOff(result_plots)]); +disp(['* Debug Plots: ', convertBoolean2OnOff(debug_plots)]); +disp(['* Triang Units: ', 'mm']); +disp('---------------------') +disp('Initializing method...') + +% Convert from MATLAB triangulation to mesh struct. +pelvisMesh.vertices = Pelvis.Points; +pelvisMesh.faces = Pelvis.ConnectivityList; + +disp('Landmarking...') + +% Run [Fischer 2019] +[~, Landmarks] = pelvicLandmarkID(pelvisMesh, ... + 'visu',debug_plots, 'CS','APP' ,'debug',debug_plots); + +% Extract the landmarks from the landmarks struct +LASIS = Landmarks.ASIS(1,:); +RASIS = Landmarks.ASIS(2,:); +LPSIS = Landmarks.PSIS(1,:); +RPSIS = Landmarks.PSIS(2,:); +SYMP = Landmarks.PS; + +% Sanity check if superior iliac spine (SIS) landmarks were correctly identified +if norm(RASIS-LASIS)(min(Alt) + 1.3*dd)); -ElmtsEpi = find(DistFem.incenter*Z0<(min(Alt) + 1.3*dd)); - -% create diaphysis triangulation -DiaFem = TriReduceMesh( DistFem, ElmtsDia ); -DiaFem = TriFillPlanarHoles(DiaFem); - -% get inertial properties from diaphysis -[V_DiaFem, DiaFem_Center] = TriInertiaPpties( DiaFem ); -Zdia = V_DiaFem(:,1); - -% create distal epiphysis triangulation -EpiFem = TriReduceMesh( DistFem, ElmtsEpi ); - -% debug plot -quickPlotTriang(EpiFem, 'b') - -%% Find Pt1 described in their method -% Pt1 : Closest epiphysis point to the diaphysis 1st inertia axis - -% Get the vectors of the epiphysis points to a point on the -% 1st inertial axis (the centroid of the diaphysis) -LinePtNodes = bsxfun(@minus, EpiFem.Points, DiaFem_Center'); - -% Get the vectors connecting the epiphysis point to their othogonal -% projection point on the 1st inertia axis -CP = (cross(repmat(Zdia',length(LinePtNodes),1),LinePtNodes)); - -% Find the closest point to the diaphysis 1st inertia axis -Dist = sqrt(sum(CP.^2,2)); -[~,IclosestPt] = min(Dist); -Pt1 = EpiFem.Points(IclosestPt,:); - -% debug plot -plot3(Pt1(1),Pt1(2),Pt1(3),'o','LineWidth',4) - -%% Find Pt2 -% curve at second location (max + 1/2 area) -[ Curves , ~, ~ ] = TriPlanIntersect(DistFem, coeff*Z0 , min(Alt) + dd ); - -% debug plot -plot3(Curves.Pts(:,1), Curves.Pts(:,2), Curves.Pts(:,3),'k'); hold on; axis equal - -% moving curve in inertial axes ref system -NewPts = Curves.Pts*V_all; -% finding the centre of the bounding box -Center_BBox = mean([min(NewPts); max(NewPts)]); -% TODO: find points posterior in axial ref syst -% post_curve = NewPts(:,3)>Center_BBox(3); -% [~, ind] = min(abs(NewPts(post_curve,2)-Center_BBox(2))); -% figure; plot3(NewPts(:,1), NewPts(:,2), NewPts(:,3),'-k'); grid on; axis equal;hold on -% plot3(Center_BBox(:,1), Center_BBox(:,2), Center_BBox(:,3),'ok'); grid on; axis equal -% plot3(NewPts(ind,1), NewPts(ind,2), NewPts(ind,3),'*k') - -CenterCS = V_all*Center_BBox'; - -%======================== -% THIS IS A MESS AND NEEDS TO BE REWRITTEN - -% Identify the posterior part of the cross section : -% assume that the closest point of the cross section -% border is posterior -% getting the nearest neighbour in Curves for the point in CenterCS -IDX = knnsearch(Curves.Pts,CenterCS'); - -% Get a raw Anterior to Posterior vector : Uap -Uap = normalizeV(Curves.Pts(IDX,:)-CenterCS'); - -% From the cross section border keep only the points that are -% posterior to the center point of the bounding box -PosteriorPts = Curves.Pts(Curves.Pts*Uap>CenterCS'*Uap,:); - -% Find the The Point P2 -% Get the orthogonal distance of the points on the posterior -% part of the cross section border to the 3rd inertia axis of -% the whole distal femur -LinePtNodes = bsxfun(@minus, PosteriorPts, CenterCS'); -CP = (cross(repmat(X0',length(LinePtNodes),1),LinePtNodes)); -Dist = sqrt(sum(CP.^2,2)); - -% Identify the closest point which is Pt2 -[~,IclosestPt] = min(Dist); -Pt2 = PosteriorPts(IclosestPt,:); -%============================================== - -% debug plot -plot3(Pt2(1),Pt2(2),Pt2(3),'ro','LineWidth',4) - -%% Define first plan iteration -npcs = normalizeV(cross( Pt1-Pt2, Y0)); - -if (Center0'-Pt1)*npcs > 0 - npcs = -npcs; -end - -ElmtsDPCs = find(EpiFem.incenter*npcs > Pt1*npcs); - -% first iteration fitted geometry -PCsFem = TriReduceMesh( EpiFem, ElmtsDPCs ); - -% debug plot -quickPlotTriang(PCsFem, 'm') - -%% First Cylinder Fit -% Initialize axis, radius and center of the cylinder -Axe0 = Y0'; -Radius0 = 0.5*(max(PCsFem.Points*npcs)-min(PCsFem.Points*npcs)); -Center0 = mean(PCsFem.Points)' - 2*npcs; -% Get a subset of the points for speed -CylPtsSS = PCsFem.Points(1:3:end,:) - -% Fit the 1st cylinder -[x0n, an, rn, d] = lscylinder(CylPtsSS, Center0, Axe0, Radius0, 0.001, 0.001); - -plotCylinder( an, rn, x0n, 15, 1, 'b') - -%% Define second plan iteration -% TODO: double check: second iteration is very similar to first -% The normal of the cut plan separting the posterior condyles is -% updated with the just identified cylinder axis -npcs = cross( Pt1-Pt2, an); -npcs = npcs'/norm(npcs); - -% Make sure the plan normal npcs points posteriorly and distally -if (Center0'-Pt1)*npcs > 0 - npcs = -npcs; -end - -% Select the points of the epiphysis lying posteriorly to the plan -ElmtsDPCs = find(EpiFem.incenter*npcs > Pt1*npcs); -PCsFem = TriReduceMesh( EpiFem, ElmtsDPCs ); - -% Second and last Cylinder Fit -Axe1 = normalizeV(an); -Radius1 = rn; -Center1 = X0n; -% Get a subset of the points for speed -CylPtsSS = PCsFem.Points(1:3:end,:); - -% Fit the 2nd cylinder -[x0n, an, rn, d] = lscylinder(CylPtsSS, Center1, Axe1, Radius1, 0.001, 0.001); - -%% Get the origin of the CS : the centroid of the fitted cylinder -% Here we take it as the middle point on the cylinder axis between the most -% lateral point projected on the axis and the most medial point projected on the -% cylinder axis. -EpiPtsOcyl_tmp = bsxfun(@minus,PCsFem.Points,x0n'); - -CylStart = min(EpiPtsOcyl_tmp*an)*an' + x0n'; -CylStop = max(EpiPtsOcyl_tmp*an)*an' + x0n'; - -CylCenter = 1/2*(CylStart + CylStop); - -% debug plot -quickPlotTriang(PCsFem, 'g') -plotCylinder( an, rn, x0n, norm(CylStart - CylStop), 1, 'r') - -CSs.Yend_Miranda = an; -CSs.Xend_Miranda = cross(an,Zdia); -CSs.Xend_Miranda = CSs.Xend_Miranda / norm(CSs.Xend_Miranda); -CSs.Zend_Miranda = cross(CSs.Xend_Miranda,CSs.Yend_Miranda); -CSs.CenterKnee_Miranda = CylCenter; - -end diff --git a/advanced_examples/upper_limb/GIBOC_humerus.m b/advanced_examples/upper_limb/GIBOC_humerus.m new file mode 100644 index 00000000..2c8ce410 --- /dev/null +++ b/advanced_examples/upper_limb/GIBOC_humerus.m @@ -0,0 +1,97 @@ +% GIBOC_HUMERUS Automatically define a reference system based on bone +% morphology by fitting analytical shapes to the articular surfaces. +% The algorithm leverages GIBOC_femur.m, of which it shared the inputs. +% The GIBOC algorithm can also extract articular +% surfaces. +%-------------------------------------------------------------------------% +% Author: Luca Modenese +% Copyright 2021 Luca Modenese & Jean-Baptiste Renault +%-------------------------------------------------------------------------% +function [BCS, JCS, HumerusBL, ArtSurf, AuxCSInfo] = GIBOC_humerus(humerusTri,... + side_raw,... + fit_method,... + result_plots,... + debug_plots,... + in_mm) + +% NOTE THAT fit_method = 'cylinder' is HARDCODED +fit_method = 'cylinder'; + +% defaults +if nargin<2; side_raw = 'r'; end +if nargin<4; result_plots = 1; end +if nargin<5; debug_plots = 0; end +if nargin<6; in_mm = 1; end +if in_mm == 1; dim_fact = 0.001; else; dim_fact = 1; end %placeholder + +% get sign correspondent to body side +[~, side_low] = bodySide2Sign(side_raw); + +% run a GIBOC_femur-cylinder analysis without plotting anything +[BCS, JCS, ~, ArtSurf, AuxCSInfo] = GIBOC_femur(humerusTri,... + side_raw,... + fit_method,... + 0,... + debug_plots,... + in_mm); + +% landmark bone according to CS but for the humerus (overwrite femoral BL) +HumerusBL = landmarkBoneGeom(humerusTri, BCS, ['humerus_', side_low]); + +% field names used in script +hip_name = ['hip_', side_low]; +knee_name = ['knee_', side_low]; +glenohumeral_name = ['glenohumeral_', side_low]; +elbow_name = ['elbow_', side_low]; + +% change structs to elbow +JCS.(glenohumeral_name) = JCS.(hip_name); +JCS.(elbow_name) = JCS.(knee_name); + +% upd ArtSurf struct +ArtSurf.(glenohumeral_name) = ArtSurf.(hip_name); +ArtSurf.(['dist_humerus_', side_low]) = ArtSurf.(['dist_femur_', side_low]); + +% remove femur structs +JCS = rmfield(JCS, {hip_name, knee_name}); +ArtSurf = rmfield(ArtSurf, {hip_name, ['dist_femur_', side_low]}); + +%% plotting + +% cut humerus for plotting +[ProxHumTri, DistHumTri] = cutLongBoneMesh(humerusTri, BCS.V(:,3)); + +label_switch=1; +if result_plots == 1 + figure('Name', ['GIBOC | bone: humerus | fit: ',fit_method,' | side: ', side_low]) + alpha = 0.5; + + % plot full femur and final JCSs + subplot(2,2,[1,3]); + plotTriangLight(humerusTri, BCS, 0) + quickPlotRefSystem(JCS.(glenohumeral_name)); + quickPlotRefSystem(JCS.(elbow_name)); + % plot humerus condyle surface + quickPlotTriang(ArtSurf.condyles_r, 'r', alpha); + + % add markers and labels + plotBoneLandmarks(HumerusBL, label_switch); + + % glenohumeral head + subplot(2,2,2); + plotTriangLight(ProxHumTri, BCS, 0); hold on + quickPlotRefSystem(JCS.(glenohumeral_name)); + plotSphere(AuxCSInfo.CenterFH_Renault, AuxCSInfo.RadiusFH_Renault, 'g', alpha); + + % distal humerus + subplot(2,2,4); + plotTriangLight(DistHumTri, BCS, 0); hold on + quickPlotRefSystem(JCS.(elbow_name)); + % plot fitting method + plotCylinder( AuxCSInfo.Cyl_Y, AuxCSInfo.Cyl_Radius, AuxCSInfo.Cyl_Pt, AuxCSInfo.Cyl_Range*1.1, alpha, 'g') + +end +grid off + +end + diff --git a/advanced_examples/upper_limb/README.md b/advanced_examples/upper_limb/README.md new file mode 100644 index 00000000..01822c8e --- /dev/null +++ b/advanced_examples/upper_limb/README.md @@ -0,0 +1,16 @@ +## Using STAPLE algorithms for upper limb applications + +### Applying `GIBOC_femur.m` to the humerus +This folder includes a simple example of how the STAPLE algorithms focussed on the lower limb can be adapted for processing also upper limb bones. +A `GIBOC_humerus.m` algorithm is drafted based on `GIBOC_femur.m` and tested by processing a humerus bone geometry (from the LHDL dataset). +The outputs of `GIBOC_femur.m` have been modified in order to: +1. obtain the correct fields on the JCS and ArtSurf structures +2. perform an appropriate landmarking of the bone +3. plot the bone and morphological items + +### Limitations +Please note that: +1. the `GIBOC_humerus.m` file has still the same console printout of `GIBOC_femur.m`. +2. this example is a proof-of-concept and it has not been tested on an extensive dataset, so we cannot ensure the robustness of `GIBOC_humerus.m`. + +![processed-humerus](../../images/adv_example_humerus.png) \ No newline at end of file diff --git a/advanced_examples/upper_limb/bone_geom/tri/humerus_r.mat b/advanced_examples/upper_limb/bone_geom/tri/humerus_r.mat new file mode 100644 index 00000000..59ab16f5 Binary files /dev/null and b/advanced_examples/upper_limb/bone_geom/tri/humerus_r.mat differ diff --git a/advanced_examples/upper_limb/process_humerus_r.m b/advanced_examples/upper_limb/process_humerus_r.m new file mode 100644 index 00000000..467326f8 --- /dev/null +++ b/advanced_examples/upper_limb/process_humerus_r.m @@ -0,0 +1,40 @@ +%-------------------------------------------------------------------------% +% Copyright (c) 2020 Modenese L. % +% Author: Luca Modenese, 2020 % +% email: l.modenese@imperial.ac.uk % +% ----------------------------------------------------------------------- % +% Script representing a proof of concept of the extensibility of STAPLE to +% the upper limb based on similarity of the humerus with the femur (long +% bone with spherical articular surface on one extremity and condilar +% structure on the other extremity. +% ----------------------------------------------------------------------- % + +clear; clc; close all +addpath(genpath('../../STAPLE')); + +%----------% +% SETTINGS % +%----------% +% folder where the various datasets (and their geometries) are located. +bones_folder = 'bone_geom'; + +% names of the bones to process with STAPLE +bones_list = {'humerus_r'}; +%-------------------------------------- + +% folder of the bone geometries in MATLAB format ('tri'/'stl') +tri_folder = fullfile(bones_folder, 'tri'); + +% create geometry set structure +geom_set = createTriGeomSet(bones_list, tri_folder); + +% body side from TriGeomSet +side = 'r'; + +% process humerus using an algorithm derived by GIBOC_femur +[HumerusCS, HumerusJCS, ~, ArtSurfFem] = GIBOC_humerus(geom_set.humerus_r, side); + +% a plot will be shown + +% remove paths +rmpath(genpath('STAPLE')); \ No newline at end of file diff --git a/bone_datasets/LHDL_CT/tri/toes_r.mat b/bone_datasets/LHDL_CT/tri/toes_r.mat new file mode 100644 index 00000000..9b93ec32 Binary files /dev/null and b/bone_datasets/LHDL_CT/tri/toes_r.mat differ diff --git a/convert_stl2triang.m b/convert_stl2triang.m index 162360bf..ce983e74 100644 --- a/convert_stl2triang.m +++ b/convert_stl2triang.m @@ -17,7 +17,7 @@ % dataset_folder = './bone_datasets/TLEM2'; % dataset_folder = './bone_datasets/ICL_MRI'; % dataset_folder = './bone_datasets/JIA_MRI'; -dataset_folder = './bone_datasets/JIA_ANKLE_MRI2'; +dataset_folder = './bone_datasets/__upper'; % folder where to store the resulting triangulations triang_folder = [dataset_folder, filesep, 'tri']; diff --git a/convert_triang2STL.m b/convert_triang2STL.m new file mode 100644 index 00000000..b39b5e73 --- /dev/null +++ b/convert_triang2STL.m @@ -0,0 +1,76 @@ +%-------------------------------------------------------------------------% +% Copyright (c) 2020 Modenese L. % +% % +% Author: Luca Modenese % +% email: l.modenese@imperial.ac.uk % +% ----------------------------------------------------------------------- % +% This script should be run when the bone geometries stored as MATLAB file +% are desired as STL files. It converts the triangulation objects stored as +% the .mat files (minimal size) as binary STL files (larger but usable in +% software like MeshLab or NMSBuilder). +% ----------------------------------------------------------------------- % +clear; clc +addpath(genpath('./STAPLE')) + +%---------------- USER'S SETTINGS ------------------ +% folder where to look for triangulation files +dataset_name = 'TLEM2_CT'; + +% main dataset folder for bone geometries +dataset_folder = ['./bone_datasets',filesep, dataset_name]; + +% folder where to store the resulting triangulations +stl_out_folder = [dataset_folder, filesep, 'stl']; +%--------------------------------------------------- + +% if triang_folder unspecified, then it is set to stl_folder plus '_tri' +if isempty(stl_out_folder) || strcmp(stl_out_folder,'') + stl_out_folder = [dataset_folder, filesep, 'stl']; +end + +% check if triang folder exists, if it doesn't then create it +if ~isfolder(stl_out_folder); mkdir(stl_out_folder); end + +% getting list of stl files from the specified folder +tri_input_folder = fullfile(dataset_folder, 'tri'); +list_of_tri = dir([tri_input_folder, filesep, '*.mat']); + +% total number of stl files +N_tri = size(list_of_tri,1); + +% looping through all the stl files found in folder +for n = 1:N_tri + + % check if directory is empty + if isempty(list_of_tri)==1 + error(['There are no tri geometry files in specified folder: ',tri_input_folder]); + end + + % read in stl file + curr_tri_name = list_of_tri(n).name; + [~, name, ext] = fileparts(curr_tri_name); + disp('-------------------') + disp(['GEOMETRY ', curr_tri_name, ' (', num2str(n),'/',num2str(N_tri),')']) + disp('-------------------') + disp(['Processing file : ', fullfile(tri_input_folder, curr_tri_name)]); + + % read traingulation + temp_tri = load(fullfile(tri_input_folder, curr_tri_name)); + + % standardize name for triangulation + tri_field_name = fields(temp_tri); + triang_geom = temp_tri.(tri_field_name{1}); + + % save it as triangulation matlab file in specified folder + curr_stl_name = [name, '.stl']; + curr_triang_file = fullfile(stl_out_folder, curr_stl_name); + disp(['Saving triangulation as STL: ', fullfile(stl_out_folder, curr_stl_name)]); + stlwrite(triang_geom, curr_triang_file) +end + +% final print +disp('-------------------') +disp('DONE') + +% remove paths +rmpath(genpath('./STAPLE')) diff --git a/images/adv_example_humerus.png b/images/adv_example_humerus.png new file mode 100644 index 00000000..90b8d838 Binary files /dev/null and b/images/adv_example_humerus.png differ diff --git a/images/reference_systems.png b/images/reference_systems.png new file mode 100644 index 00000000..1c6f0877 Binary files /dev/null and b/images/reference_systems.png differ diff --git a/tests/test_all_methods.m b/tests/test_all_methods.m index 6e049d8f..b068984a 100644 --- a/tests/test_all_methods.m +++ b/tests/test_all_methods.m @@ -54,19 +54,19 @@ %---- FEMUR ----- if isfield(geom_set, femur_name) % [FemurCS0, JCS0] = Miranda2010_buildfACS(geom_set.(femur_name)); - [FemurCS1, FemurJCS1] = Kai2014_femur(geom_set.(femur_name), side); - [FemurCS2, FemurJCS2] = GIBOC_femur(geom_set.(femur_name), side, 'spheres'); - [FemurCS3, FemurJCS3] = GIBOC_femur(geom_set.(femur_name), side, 'ellipsoids'); - [FemurCS4, FemurJCS4] = GIBOC_femur(geom_set.(femur_name), side, 'cylinder'); + [FemurCS1, FemurJCS1, FemurBL1] = Kai2014_femur(geom_set.(femur_name), side); + [FemurCS2, FemurJCS2, FemurBL2] = GIBOC_femur(geom_set.(femur_name), side, 'spheres'); + [FemurCS3, FemurJCS3, FemurBL3] = GIBOC_femur(geom_set.(femur_name), side, 'ellipsoids'); + [FemurCS4, FemurJCS4, FemurBL4] = GIBOC_femur(geom_set.(femur_name), side, 'cylinder'); end %---- TIBIA ----- if isfield(geom_set, tibia_name) % [TibiaCS0, JCS0] = Miranda2010_buildtACS(geom_set.(tibia_name)); - [TibiaCS1, TibiaJCS1] = Kai2014_tibia(geom_set.(tibia_name), side); - [TibiaCS2, TibiaJCS2] = GIBOC_tibia(geom_set.(tibia_name), side, 'plateau'); - [TibiaCS3, TibiaJCS3] = GIBOC_tibia(geom_set.(tibia_name), side, 'ellipse'); - [TibiaCS4, TibiaJCS4] = GIBOC_tibia(geom_set.(tibia_name), side, 'centroids'); + [TibiaCS1, TibiaJCS1, TibiaBL1] = Kai2014_tibia(geom_set.(tibia_name), side); + [TibiaCS2, TibiaJCS2, TibiaBL2] = GIBOC_tibia(geom_set.(tibia_name), side, 'plateau'); + [TibiaCS3, TibiaJCS3, TibiaBL3] = GIBOC_tibia(geom_set.(tibia_name), side, 'ellipse'); + [TibiaCS4, TibiaJCS4, TibiaBL4] = GIBOC_tibia(geom_set.(tibia_name), side, 'centroids'); end % %---- PATELLA ----- @@ -84,7 +84,7 @@ % % %---- CALCANEUS/SUBTALAR ----- if isfield(geom_set, calcn_name) - [FootCS, FootJCS] = STAPLE_foot(geom_set.(calcn_name), side); + [FootCS, FootJCS, FootBL] = STAPLE_foot(geom_set.(calcn_name), side); end % close all