Possible bug in GIBOC_femur_processEpiPhysis.m line 67

[RiadAkhundov]

Issue: Possible bug in GIBOC_femur_processEpiPhysis.m line 67

64 % Are CSs.Y0 and U_Axes pointing in the same direction? 65 Orientation = round(mean(sign(U_Axes*CSs.Y0))); 66 % Make all the axes point in the Laterat -> Medial direction 67 U_Axes = Orientation*U_Axes;

Logic: If CSs.Y0 and U_Axes are not pointing in the same direction Orientation*U_Axes = 0 and prevents the script from progressing. If they do point in the same direction Orientation*U_Axes = 1*U_Axes, meaning line 67 is either superfluous, or breaks the script.

Fix: 67 %U_Axes = Orientation*U_Axes; Commenting out line 67 fixed the issue for me, but I may be wrong (and this might be an underlying issue with my .stl, although the script runs smoothly after commenting out that line), so please double-check it :)

[modenaxe]

Hi Riad, thank you for reporting the issue! if you can share the stl that caused the issue that would really help debugging. @renaultJB this seem a GIBOC inherited issue.

[RiadAkhundov]

Hi Luca,

Here is the link to my .stl files and the hip model example with the settings I used:

https://mega.nz/file/gdAHSIYJ#03SDkiGnCApGg2sIRLW68fHIJrV-IzGTABVcO5WW72Y

Cheers, Riad

---

From: Luca Modenese notifications@github.com Sent: Wednesday, 17 February 2021 3:29 AM To: modenaxe/msk-STAPLE msk-STAPLE@noreply.github.com Cc: Riad Akhundov Riad.Akhundov@uon.edu.au; Author author@noreply.github.com Subject: Re: [modenaxe/msk-STAPLE] Possible bug in GIBOC_femur_processEpiPhysis.m line 67 (#87)

Hi Riad, thank you for reporting the issue! if you can share the stl that caused the issue that would really help debugging. @renaultJBhttps://github.com/renaultJB this seem a GIBOC inherited issue.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/modenaxe/msk-STAPLE/issues/87#issuecomment-779996071, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AMANBH4HMZ2VY57KP4QLESDS7KTPPANCNFSM4XVSCNVA.

[modenaxe]

@RiadAkhundov I think you are right.

@renaultJB in this function we are identifying the lines between condyles that are grown into the condyles:

64 % Are CSs.Y0 and U_Axes pointing in the same direction? 65 Orientation = round(mean(sign(U_Axes*CSs.Y0))); 66 % Make all the axes point in the Laterat -> Medial direction 67 U_Axes = Orientation*U_Axes;

I think the comments (probably mine) are misleading and L67 indeed doesn't do much: as sign(U_Axes*CSs.Y0) yields axes in both directions, L67 simply flips them.

@renaultJB shall we comment it out?

[renaultJB]

@RiadAkhundov Thank you for letting us know.

Orientation is supposed to be -1 or +1, but you are right nothing prevents it to be equal to 0. But Orientation being 0 is a bad signal in my opinion. It means that dot products of U_Axes vectors with CSs.Y0 are sometimes positive and sometimes negative, whileU_Axes vectors should all point to the same rough direction. Direction, that should be roughly paralllel to CSs.Y0.

U_Axes is actually a collection of vectors that are roughly directed from lateral to medial or from medial to lateral. So, I think the error come from the unnecessary reorientation operation that could have been dealt with earlier in the code. We know that CSs.Y0 is directed from lateral to medial. Line 57 we remove the duplicate of each axes in the opposite direction. We should just do instead:

% Remove duplicate Axes that are not directed from Lateral to Medial (CSs.Y0)
I_Axes_duplicate = find(Axes*CSs.Y0 < 0);

IdCdlPts(I_Axes_duplicate,:)=[];
Axes(I_Axes_duplicate,:)=[];

%Normalize Axes to get unitary vectors
U_Axes = Axes./repmat(sqrt(sum(Axes.^2,2)),1,3);

This part is not useful anymore

% Are CSs.Y0 and U_Axes pointing in the same direction?
Orientation = round(mean(sign(U_Axes*CSs.Y0)));
% Make all the axes point in the Laterat -> Medial direction
U_Axes = Orientation*U_Axes;
% Axes = Orientation*Axes;

However Orientation is used at the end of the function to switch values within an outputed array, so we should either modify the end part of the function or leave Orientation=round(mean(sign(U_Axes*CSs.Y0))); and raise en error or a warning if Orientation is different than 1.

@modenaxe What do you think ?

@RiadAkhundov I'll try to test this modification over the WE, but feel free to test it and let us know.

[renaultJB]

There is a lot small U_Axes (black lines) roughly orthogonals to the expected Lateral -> Medial axis this is what caused the error I think.

Code used to make the figure (within the function)

figure;
trisurf(EpiFem,'facealpha',0.4,'facecolor','y',...
        'edgecolor','k'); hold on
plot3(EpiFem.Points(IdCdlPts(:,1),1), EpiFem.Points(IdCdlPts(:,1),2), EpiFem.Points(IdCdlPts(:,1),3),'r*');
plot3(EpiFem.Points(IdCdlPts(:,2),1), EpiFem.Points(IdCdlPts(:,2),2), EpiFem.Points(IdCdlPts(:,2),3),'b*');
for i = 1:length(IdCdlPts(:,1))
    plot3(EpiFem.Points(IdCdlPts(i,:),1), EpiFem.Points(IdCdlPts(i,:),2), EpiFem.Points(IdCdlPts(i,:),3),'k-','LineWidth',3);
end
axis equal;

[renaultJB]

After applying the proposed changes, and displaying the same figure at the end of the function we obtain what we are expecting.

After testing on Riad dataset. I've made a pull request, If you can have a look @modenaxe before merging.

And @RiadAkhundov you can either copy the code in the pull request or clone again the repo once it has been merged with master. Thank you again for pointing out this issue, the code is now more robust thanks to you.

[modenaxe]

@renaultJB thank you! I have double checked and run the current tests - all seems good. I was initially concerned about right/left leg pointing axes, but it is actually ok. I have also added to the GIBOC_femur_processEpiPhysis.m the plotting code, so next time it will be easier to debug similar issues using the debug_plots. @RiadAkhundov please git pull and let us know if you have further issues!