Commit 33ca45b4 authored by Mark van Turnhout's avatar Mark van Turnhout
Browse files

stripped lsm_get_ellipsoid for bas_estIndDir

parent 041d68c3
......@@ -106,7 +106,7 @@ for b = 1:numel(basil)
% OK, let's not optimise at all, and rely on the mmmoi for now
% com and eigenvectors of the mass moment of inertia tensor
[~,com,~,ev,~] = lsm_get_ellipsoid(X(Li),Y(Li),Z(Li)); %
[com,ev] = bas_getIndDir(X(Li),Y(Li),Z(Li)); %
% indenter is mostly in the xy-plane, I hope
idir = ev(:,ev(3,:) == min(ev(3,:)) );
% indenter always points in the positive y-direction, I hope
......
function [com, ev] = bas_getIndDir(Xr,Yr,Zr,m)
%
% [com, ev] = bas_getIndDir(X,Y,Z,m)
%
% Calculates eigenvectors of the mass moment of inertia tensor of the
% (indenter) voxels at the locations [X Y Z]. One eigenvector points along
% the long axis of the indenter.
%
% X, Y and Z are columns of equal length with indenter voxel coordinates.
% The output com is the centre of mass of the indenter voxels, ev is the
% matrix with eigenvectors. Adapted (stripped) from lsm_get_ellipsoid.
%
% See also, lsm_get_ellipsoid, bas_estCylPos, bas_estIndDisp, bas_prepFEM
if nargin == 3, m = ones(size(Zr)); end
if nargin < 3
if nargin == 2
m = Yr;
else
if size(Xr,2) == 3
m = ones(size(Xr,1),1);
elseif size(Xr,2) == 4
m = Zr(:,4);
end
end
Zr = Xr(:,3);
Yr = Xr(:,2);
Xr = Xr(:,1);
end
M = sum(m);
% centre of mass
com = [sum(m.*Xr) sum(m.*Yr) sum(m.*Zr)]/M;
% set com-coordinates to [0, 0, 0]
Xr = Xr - com(1);
Yr = Yr - com(2);
Zr = Zr - com(3);
% I's
Ixx = sum(m.*(Yr.^2+Zr.^2));
% Iyy
Iyy = sum(m.*(Xr.^2+Zr.^2));
% Izz
Izz = sum(m.*(Xr.^2+Yr.^2));
% Ixy
Ixy = -sum(m.*Xr.*Yr);
% Ixz
Ixz = -sum(m.*Xr.*Zr);
% Iyz
Iyz = -sum(m.*Yr.*Zr);
I = zeros(3,3);
I(1,1) = Ixx; I(2,2) = Iyy; I(3,3) = Izz;
I(1,2) = Ixy; I(2,1) = Ixy;
I(1,3) = Ixz; I(3,1) = Ixz;
I(2,3) = Iyz; I(3,2) = Iyz;
[ev, ~] = eig(I);
end
......@@ -138,7 +138,7 @@ BasilLab is very lazy on the masking part.
The function \bas{maskBG}{} loads the voxel set that you wish to segment and shows a maximum intensity projection over the $z$-direction. In this image, you are asked to interactively select the centre and radius of a circle. The circle will serve as the basis of a cylinder-shaped 3D `mask' for \texttt{uiContourSegment}. \bas{maskBG}{} does not make separate mask for skin and bone.
The mask for \basilhome\texttt{/<basilid>/}\texttt{MRI\_slices\_<foo>.mat}\footnote{\texttt{<foo>} can be \texttt{preload} or \texttt{loading} or \texttt{postload} or \texttt{flash} (see table \ref{mritomatout}.}will be saved as \basilhome\texttt{/}\-\texttt{<basilid>}\-\texttt{/MRI\_BG\_<foo>.mat}. To mask the data in the scans before (\texttt{1}), during (\texttt{2}) and after indentation (\texttt{3}) for \texttt{<basilid> = 140611} call:
The mask for \basilhome\texttt{/<basilid>/}\texttt{MRI\_slices\_<foo>.mat}\footnote{\texttt{<foo>} can be \texttt{preload} or \texttt{loading} or \texttt{postload} or \texttt{flash} (see table \ref{mritomatout}).}will be saved as \basilhome\texttt{/}\-\texttt{<basilid>}\-\texttt{/MRI\_BG\_<foo>.mat}. To mask the data in the scans before (\texttt{1}), during (\texttt{2}) and after indentation (\texttt{3}) for \texttt{<basilid> = 140611} call:
\begin{lstlisting}[numbers=none]
>> bas_maskBG(140611, [1 2 3]);
\end{lstlisting}
......@@ -163,9 +163,9 @@ The next step is to map the displacement and rotation of the tibia as a result o
The function \bas{mapBone} uses the segmented contours \texttt{contour\_bone\_<foo>-ncc.mat} to map the movement and displacement of the tibia in the \texttt{loading}-contours and in the \texttt{postload}-contours, both with the \texttt{preload}-contours as reference.
The output of \bas{mapBone}{} consists of a $3\times 3$ matrix for each of the two cases. The first row contains the centre of mass $(x,y,z)$ of the tibia contours (in mm) in the \texttt{preload}-contours, and the second row contains the centre of mass $(x,y,z)$ of the tibia contours in the \texttt{loading}- or \texttt{postload}-contours. The third row contains three `Euler angles' $(\alpha, \beta, \gamma)$ that describe the rotation of the tibia around its centre of mass.
The output of \bas{mapBone}{} consists of a $3\times 3$ matrix for each of the two cases. The first row contains the centre of mass $[x\;y\;z]$ of the tibia contours (in mm) in the \texttt{preload}-contours, and the second row contains the centre of mass of the tibia contours in the \texttt{loading}- or \texttt{postload}-contours. The third row contains three `Euler angles' $[\alpha\;\beta\;\gamma]$ that describe the rotation of the tibia around its centre of mass. The angles are estimated with a minimisation algorithm that runs on the function \bas{estBoneRot}.
These Euler angles follow the \gibbon-convention (see \texttt{euler2DCM.m}): $\alpha$ is the rotation around the $x$-axis, $\beta$ the rotation around the $y$-axis, and $\gamma$ the rotation around the $z$-axis, such that the 3D rotation matrix is given by:
These Euler angles follow the \gibbon-convention (see \texttt{euler2DCM.m}): $\alpha$ is the rotation around the $x$-axis, $\beta$ the rotation around the $y$-axis, and $\gamma$ the rotation around the $z$-axis, such that the 2D $(3\times 3)$ rotation matrix that describes the 3D rotation is given by:
\begin{equation}
\mat{R}_{\alpha\beta\gamma} = \mat{R}_\alpha\cdot\mat{R}_\beta \cdot \mat{R}_\gamma
\end{equation}
......@@ -176,10 +176,12 @@ The $3\times 3$ matrix (centres of mass, Euler angles) is written to \texttt{Bon
>> bas_mapBone([140611, 140613]);
\end{lstlisting}
\subsection{Map indenter movement}\label{mapind}
\subsection{\bas{prepFEM}}
The indenter movement is of course needed as boundary condition in the FEM. The function \bas{estIndDisp}{}
\subsection{Putting it all together: \bas{prepFEM}}
\section{Building the FE model}\label{buildmodel}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment