workflow.tex 33.9 KB
Newer Older
1
2
% !TeX root = BasilLab.tex
\begin{savequote}
Mark van Turnhout's avatar
Mark van Turnhout committed
3
Modelling is a scientific technique that requires a good deal of intuition and insight to be really effective. Nevertheless, it is the most powerful tool available to the experimenter from the instrumentarium. \qauthor{Rom Harr{\'e} \cite{Harre2009}}
4
\end{savequote}
Mark van Turnhout's avatar
Mark van Turnhout committed
5
\chapter{Work flow}
Mark van Turnhout's avatar
Mark van Turnhout committed
6

Mark van Turnhout's avatar
Mark van Turnhout committed
7
8
9
10
11
%Before you can use BasilLab, you wil need to take two steps. First, add the directory that contains the (Matlab) scripts, possibly <basillab> to your Matlab search path, e.g. go (in Matlab) to that directory and issue \texttt{addpath(pwd);} on the command line. 
%
%Second, tell your local copy of BasilLab where the raw MRI-data is to be found, and where to store and for look processed MRI-data. You should use the functions \texttt{bas\_setMRIroot} and \texttt{bas\_setDataRoot} to do that (once). 
%
%The directory with the raw MRI-data (\texttt{bas\_setMRIroot}) remains untouched. The MRI-data is read and copied to the directory \texttt{<basilid>} which is created in the directory \texttt{<basilhome>} (\texttt{bas\_setDataRoot}). The link between the raw MRI-data folder inside  \texttt{<MRI-root>} and  \texttt{<basilid>} is \textsl{hard-coded} in \texttt{bas\_getMeta.m}. 
Mark van Turnhout's avatar
Mark van Turnhout committed
12

Mark van Turnhout's avatar
Mark van Turnhout committed
13
14
15
After initialisation of BasilLab, the work flow can be roughly divided into four steps:
\begin{enumerate}[start=0]
\item Initialise BasilLab (once)
Mark van Turnhout's avatar
Mark van Turnhout committed
16
17
18
19
20
21
\item Pre-processing of (raw) MRI-data. Including reading and copying (raw) MRI-data, segmentation of MR images,and processing (segmented) MRI-data. Requires Matlab (and human intervention).
\item Building the FE model. Requires Matlab and Abaqus.
\item Running the FE model, Requires Abaqus.
\item Post-processing of FE results. Requires Abaqus.
\end{enumerate}

Mark van Turnhout's avatar
Mark van Turnhout committed
22
23
24
25
26
27
28
29
30
31
32
33

\setcounter{section}{-1}

\section{Initialise BasilLab}\label{basilini}

Matlab needs to be aware of the location of the BasilLab-scripts,  so you first have to add the directory that contains the (Matlab) scripts, possibly \basilhome{} to your Matlab search path, e.g.\ go (in Matlab) to that directory and issue \texttt{addpath(pwd);} on the command line. 

As a matter of principle and for philosophical reasons, BasilLab leaves the raw MRI-data completely untouched. You have to initialise BasilLab \textsl{once} by defining the path to the raw MRI-data (\MRIroot{}), and the path to the BasilLab work directory (\basilhome{}, see also section \ref{dirstruct}). The necessary MRI (meta) data will be copied from \MRIroot{} to a unique folder \texttt{<basilid>} inside \basilhome{} for each animal. All (model) information that is extracted from the MRI-data by BasilLab will be written to, and sought in, that folder: \basilhome{}\texttt{/<basilid>}.\\

\noindent \MRIroot{} is defined (hard-coded) in \bas{getMRIroot}. This script is called by BasilLab when it needs to know the path to the \MRIroot{}-folder. You can edit this script manually, or you can use \bas{setMRIroot} to update this script to your local situation. When you call  \bas{setMRIroot}, the \bas{getMRIroot}-script is overwritten with a new version that provides the \MRIroot{} that you provided as an argument to \bas{setMRIroot}. 

Similarly, \basilhome{} is defined in \bas{getDataRoot} and updated with \bas{setDataRoot}. The argument for \bas{setDataRoot} or \bas{setMRIroot} is the \textsl{full} path to the folder that you wish to define, e.g.:
Mark van Turnhout's avatar
Mark van Turnhout committed
34
\begin{lstlisting}[numbers=none, language=basillab]
Mark van Turnhout's avatar
Mark van Turnhout committed
35
36
37
38
>> bas_setMRIroot('/home/mark/BrukerData');
>> bas_setDataRoot(pwd);
\end{lstlisting}

39
40
\section{Pre-processing of MRI-data}

41
42
All steps that are required to prepare the model data can be performed with \bas{prepFEM}. Briefly, these steps are
\begin{enumerate}
Mark van Turnhout's avatar
Mark van Turnhout committed
43
\item Copy MRI-data from \MRIroot{}\texttt{/<experiment-folder>} to \basilhome\texttt{/<basilid>}
44
45
46
47
48
49
50
51
\item Segment MRI-data to obtain skin- and bone contours
\item Determine bone movement and  rotation due to indentation from segmented bone contours
\item Determine indenter orientation and displacements from MRI-data
\end{enumerate}
Only step 2. requires user interaction, the rest is performed 'under the hood', automagically.

\subsection{Copying MRI data}

Mark van Turnhout's avatar
Mark van Turnhout committed
52
53
54
55
Extracting and copying MRI (meta) data is done by \bas{MRItoMat}. This function calls a number of other functions: \bas{getMRIroot} for the path to the raw MRI-data, \bas{getDataRoot} for the path to the BasilLab work directory, and \bas{getMeta} for the experiment specific sub-folders within these two directories.  \bas{getMeta} also provides information on the scans that are to be processed in the raw MRI-data.

The actual parsing of raw MRI-data is done by a call to \bas{readMRI}, which in turn calls \bas{readMRImeta} to parse the MRI meta data. Note that MRI data parsing in BasilLab is copied (and slightly adapted) from older scripts that existed long before BasilLab: \texttt{jcampread.m} by B.C.\ Hamans (2003) and Edwin Heijman (update 2005), \texttt{load2seq.m} by Edwin Heijman (July 2005), and \texttt{PV3read.m} by B.J.\ van Nierop (March 2008).\\

56
\noindent The first argument of \bas{MRItoMat} is a list of \texttt{<basilid>}'s for which to parse and copy the MRI-data. The second argument can be used to skip copying and parsing of data that already exists in \basilhome\texttt{/<basilid>}. When the list of \texttt{<basilid>}'s is omitted, all data defined in \bas{getMeta} is parsed. When the last argument is \texttt{0}, data that is already present in  \basilhome\texttt{/<basilid>} is not processed; when it is  \texttt{1}, all data present in  \basilhome\texttt{/<basilid>} will be overwritten.  If \bas{MRItoMat} is called without arguments, all data defined in \bas{getMeta} is parsed and all existing data is overwritten.
Mark van Turnhout's avatar
Mark van Turnhout committed
57
58

To process the data for \texttt{<basilid> = 140611} and \texttt{<basilid> = 140613} and skip processing of data that is already present  in  \basilhome\texttt{/<basilid>}, call:
Mark van Turnhout's avatar
Mark van Turnhout committed
59
\begin{lstlisting}[numbers=none, language=basillab]
Mark van Turnhout's avatar
Mark van Turnhout committed
60
61
62
63
>> bas_MRItoMat([140611, 140613], 0);
\end{lstlisting}

The output of \bas{MRItoMat} consists of the files listed in table \ref{mritomatout}. These files are written to \basilhome\texttt{/<basilid>} as explained before. The bold entries are actually needed by BasilLab to prepare the model data. The other (meta) data is written for your convenience and backward compatibility.
Mark van Turnhout's avatar
Mark van Turnhout committed
64
\begin{table}[p!]
Mark van Turnhout's avatar
Mark van Turnhout committed
65
\center
66
\caption{Output of \bas{MRItoMat}: \texttt{`\_preload'} is from the scans before indentation, \texttt{`\_loading'} is from the scans during indentation, and \texttt{`\_postload'} is from the scans after removal of the indentation;  \texttt{`\_flash'} is from the flash scans. \textbf{Bold} entries are \textbf{required} to prepare model data, \textsl{slanted} entries are used for parameters that should facilitate \textsl{post-processing} of the FE data.}\label{mritomatout}
Mark van Turnhout's avatar
Mark van Turnhout committed
67
68
\begin{tabulary}{\linewidth}{l R}
\hline\noalign{\smallskip}
Mark van Turnhout's avatar
Mark van Turnhout committed
69
70
file in \basilhome\texttt{/<basilid>/}& description\\  
\noalign{\smallskip}\hline\noalign{\smallskip}
71
72
73
74
75
76
77
78
79
\textbf{MRI\_slices\_preload.mat\index{MRI\_slices\_preload.mat@\texttt{MRI\_slices\_preload.mat}}} &   \multirow{4}*{3D matrix $(256\times 256\times 20)$ with MRI-slice voxel values}\\
\textbf{MRI\_slices\_loading.mat\index{MRI\_slices\_loading.mat@\texttt{MRI\_slices\_loading.mat}}} & \\
\textsl{MRI\_slices\_postload.mat\index{MRI\_slices\_postload.mat@\texttt{MRI\_slices\_postload.mat}}} & \\ 
\texttt{MRI\_slices\_flash.mat\index{MRI\_slices\_flash.mat@\texttt{MRI\_slices\_flash.mat}}} & \\ \hline\noalign{\smallskip}
\texttt{IMDAT\_preload.mat\index{IMDAT\_<foo>.mat@\texttt{IMDAT\_<foo>.mat}}} & as \texttt{MRI\_slices\_preload.mat\index{MRI\_slices\_preload.mat@\texttt{MRI\_slices\_preload.mat}}} \\
\texttt{IMDAT\_loading.mat\index{IMDAT\_<foo>.mat@\texttt{IMDAT\_<foo>.mat}}} & as \texttt{MRI\_slices\_loading.mat\index{MRI\_slices\_loading.mat@\texttt{MRI\_slices\_loading.mat}}} \\
\texttt{IMDAT\_postload.mat\index{IMDAT\_<foo>.mat@\texttt{IMDAT\_<foo>.mat}}} & as \texttt{MRI\_slices\_postload.mat\index{MRI\_slices\_postload.mat@\texttt{MRI\_slices\_postload.mat}}} \\
\texttt{IMDAT\_flash.mat\index{IMDAT\_<foo>.mat@\texttt{IMDAT\_<foo>.mat}}} & as \texttt{MRI\_slices\_flash.mat\index{MRI\_slices\_flash.mat@\texttt{MRI\_slices\_flash.mat}}} \\ \hline\noalign{\smallskip}
\texttt{MRI\_subject\_preload.mat} &  \multirow{4}*{`subject' from MRI meta data\index{MRI\_subject\_<foo>.mat@\texttt{MRI\_subject\_<foo>.mat}}}\\
Mark van Turnhout's avatar
Mark van Turnhout committed
80
81
82
\texttt{MRI\_subject\_loading.mat} & \\
\texttt{MRI\_subject\_postload.mat} &\\ 
\texttt{MRI\_subject\_flash.mat} &\\\hline\noalign{\smallskip}
83
\texttt{MRI\_acqp\_preload.mat} &  \multirow{4}*{`acqp' (aquisition parameters) from MRI meta\index{MRI\_acqp\_<foo>.mat@\texttt{MRI\_acqp\_<foo>.mat}}}\\
Mark van Turnhout's avatar
Mark van Turnhout committed
84
85
86
\texttt{MRI\_acqp\_loading.mat} & \\
\texttt{MRI\_acqp\_postload.mat} & \\ 
\texttt{MRI\_acqp\_flash.mat} & \\ \hline
87
\texttt{MRI\_method\_preload.mat}&\multirow{4}*{`method' from MRI meta data\index{MRI\_method\_<foo>.mat@\texttt{MRI\_method\_<foo>.mat}}} \\
Mark van Turnhout's avatar
Mark van Turnhout committed
88
89
90
\texttt{MRI\_method\_loading.mat} & \\
\texttt{MRI\_method\_postload.mat} &  \\ 
\texttt{MRI\_method\_flash.mat} &  \\ \hline\noalign{\smallskip} 
91
\texttt{MRI\_imnd\_preload.mat} & \multirow{4}*{`imnd' from MRI meta data\index{MRI\_imnd\_<foo>.mat@\texttt{MRI\_imnd\_<foo>.mat}}} \\
Mark van Turnhout's avatar
Mark van Turnhout committed
92
93
94
\texttt{MRI\_imnd\_loading.mat} &  \\
\texttt{MRI\_imnd\_postload.mat} & \\
\texttt{MRI\_imnd\_flash.mat} & \\  \hline\noalign{\smallskip}
95
\texttt{MRI\_reco\_preload.mat} & \multirow{4}*{`reco' from MRI meta data\index{MRI\_reco\_<foo>.mat@\texttt{MRI\_reco\_<foo>.mat}}} \\
Mark van Turnhout's avatar
Mark van Turnhout committed
96
97
98
\texttt{MRI\_reco\_loading.mat} & \\
\texttt{MRI\_reco\_postload.mat} & \\ 
\texttt{MRI\_reco\_flash.mat} & \\ \hline\noalign{\smallskip}
99
\texttt{MRI\_d3proc\_preload.mat} &  \multirow{4}*{`d3proc' from MRI meta data\index{MRI\_d3proc\_<foo>.mat@\texttt{MRI\_d3proc\_<foo>.mat}}} \\
Mark van Turnhout's avatar
Mark van Turnhout committed
100
101
102
\texttt{MRI\_d3proc\_loading.mat} &  \\
\texttt{MRI\_d3proc\_postload.mat} &  \\ 
\texttt{MRI\_d3proc\_flash.mat} &  \\ \hline\noalign{\smallskip}
103
\texttt{MRI\_visupars\_preload.mat} & \multirow{4}*{`visupars' (visualisation parameters) from MRI meta data\index{MRI\_visupars\_<foo>.mat@\texttt{MRI\_visupars\_<foo>.mat}}} \\
Mark van Turnhout's avatar
Mark van Turnhout committed
104
105
106
\texttt{MRI\_visupars\_loading.mat} & \\
\texttt{MRI\_visupars\_postload.mat} &  \\ 
\texttt{MRI\_visupars\_flash.mat} &  \\\hline\noalign{\smallskip}
107
\textbf{MRI\_voxel\_preload.txt} & \multirow{4}*{row with the voxel sizes in $[x y z]$  in mm from MRI meta data\index{MRI\_voxel\_<foo>.mat@\texttt{MRI\_voxel\_<foo>.mat}}}\\
Mark van Turnhout's avatar
Mark van Turnhout committed
108
109
110
\textbf{MRI\_voxel\_loading.txt} & \\
\textsl{MRI\_voxel\_postload.txt} & \\
\texttt{MRI\_voxel\_flash.txt} & \\
Mark van Turnhout's avatar
Mark van Turnhout committed
111
112
113
\noalign{\smallskip}\hline 
\end{tabulary}
\end{table}
114
115
116


\subsection{Segment MRI-data}
Mark van Turnhout's avatar
Mark van Turnhout committed
117

118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
We need to segment five contours in the MR images (table \ref{contours}): skin contours in de scans before and during indentation, and bone contours in the scans before, during and after indentation.
\begin{table}[t!]
\center
\caption{Contours that need to be segmented in the MRI-data}\label{contours}
\begin{tabulary}{\linewidth}{l c R l}
\hline\noalign{\smallskip}
contour & scans & purpose & section \\  
\noalign{\smallskip}\hline\noalign{\smallskip}
\multirow{2}*{skin} & before indentation & geometry for the FEM & \ref{buildmodel}\\
 & during indentation & determine end-point of the indenter, boundary condition for the FEM & \ref{mapind}\\ \hline\noalign{\smallskip}
\multirow{3}*{bone} & before indentation & geometry for the FEM & \ref{buildmodel}\\
& during indentation & determine tibia movement as a result of the application of indentation, boundary condition for the FEM & \ref{mapbone}\\
& after indentation & determine tibia movement as a result of the release of indentation, used in post-processing of the results of the FEM & \ref{mapbone}\\
\noalign{\smallskip}\hline 
\end{tabulary}
\end{table}

Mark van Turnhout's avatar
Mark van Turnhout committed
135
For the segmentation of the MR image data we rely on the \gibbon{}-function \gib{uiContourSegment}. This function requires a (3D) `mask' to go with the 3D MR voxel data that is to be segmented. The mask should have the same size as the voxel data, and can change between (2D) MRI-slices. 
136
137
138

BasilLab is very lazy on the masking part.

Mark van Turnhout's avatar
Mark van Turnhout committed
139
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 \gib{uiContourSegment}.  \bas{maskBG}{} does not make separate mask for skin and bone.
140

141
The mask for \basilhome\texttt{/<basilid>/}\texttt{MRI\_slices\_<foo>.mat\index{MRI\_slices\_<foo>.mat@\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\index{MRI\_BG\_<foo>.mat@\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:
Mark van Turnhout's avatar
Mark van Turnhout committed
142
\begin{lstlisting}[numbers=none, language=basillab]
143
144
145
>> bas_maskBG(140611, [1 2 3]);
\end{lstlisting}

146
\noindent With the masks, we can call \bas{segmentMRI} which in turn calls \gibbon{} for the actual segmentation. The output of \gib{uiContourSegment} is a Matlab structure called \texttt{Vcs} that contains a cell with contour points $(x, y,z)$ for each 2D MRI-slice. \bas{segmentMRI} saves this structure `as is'.  The \texttt{Vcs}-structure for the skin in \basilhome\texttt{/<basilid>/}\texttt{MRI\_slices\_<foo>.mat\index{MRI\_slices\_<foo>.mat@\texttt{MRI\_slices\_<foo>.mat}}} will be saved as  \basilhome\texttt{/<basilid>/}\texttt{contour\_skin\_<foo>.mat\index{contour\_skin\_<foo>.mat@\texttt{contour\_skin\_<foo>.mat}}}, that of the tibia as \basilhome\texttt{/<basilid>}\-\texttt{/contour\_bone\_<foo>.mat\index{contour\_bone\_<foo>.mat@\texttt{contour\_bone\_<foo>.mat}}}.
147

148
Note that \bas{segmentMRI} loaded \basilhome\texttt{/<basilid>/}\texttt{MRI\_voxel\_<foo>.txt\index{MRI\_voxel\_<foo>.mat@\texttt{MRI\_voxel\_<foo>.mat}}}  in the process. The dimensions of the saved contours are in mm.  To segment the skin  and bone data in the scans before (\texttt{1,2}), during (\texttt{3,4}), and the bone contour after indentation (\texttt{6}) for \texttt{<basilid> = 140611} call:
Mark van Turnhout's avatar
Mark van Turnhout committed
149
\begin{lstlisting}[numbers=none,language=basillab]
150
151
152
>> bas_segmentMRI(140611, [1:4, 6]);
\end{lstlisting}

153
\noindent Finally, we need to check for missing contours in de segmented data. The  \texttt{Vcs}-structure of  \gib{uiContourSegment} may contain empty cells, and missing contour data may mess op BasilLab and post-processing. The function \bas{checkVsize} loads five contour-sets \texttt{contour\_skin\_preload.mat\index{contour\_skin\_preload.mat@\texttt{contour\_skin\_preload.mat}}}, \texttt{contour\_skin\_loading.mat\index{contour\_skin\_loading.mat@\texttt{contour\_skin\_loading.mat}}}, \texttt{contour\_bone\_preload.mat}, \texttt{contour\_bone\_loading.mat\index{contour\_bone\_loading.mat@\texttt{contour\_bone\_loading.mat}}}, and \texttt{contour\_bone\_postload.mat\index{contour\_bone\_postload.mat@\texttt{contour\_bone\_postload.mat}}} from \basilhome\texttt{/<basilid>/}, and deletes a contour that is missing from \textsl{any} set from \textsl{all} sets. 
154

155
The (new) contour sets are saved with \texttt{-ncc} appended to their file name, i.e.\  the corrected  set for \texttt{contour\_skin\_preload.mat\index{contour\_skin\_preload.mat@\texttt{contour\_skin\_preload.mat}}} will be saved as  \texttt{contour\_skin\_preload-ncc.mat\index{contour\_skin\_preload-ncc.mat@\texttt{contour\_skin\_preload-ncc.mat}}}. The original data will remain untouched. To post-process the contours  for \texttt{<basilid> = 140611} and \texttt{<basilid> = 140613} call:
Mark van Turnhout's avatar
Mark van Turnhout committed
156
\begin{lstlisting}[numbers=none,language=basillab]
157
158
159
160
161
162
163
164
165
>> bas_checkVsize([140611, 140613]);
\end{lstlisting}

\subsection{Map bone movement}\label{mapbone}

The next step is to map the displacement and rotation of the tibia as a result of the application and release of the indentation. The movement of the tibia as a result of the indentation is modelled in the FEM in the loading phase and is therefore needed for its boundary conditions. The movement of the tibia as a result of the release of the indentation is used in post-processing of the FEM output.

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. 

166
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}.
167

Mark van Turnhout's avatar
Mark van Turnhout committed
168
These Euler angles follow the \gibbon-convention (see \gib{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:
169
170
171
\begin{equation}
\mat{R}_{\alpha\beta\gamma} = \mat{R}_\alpha\cdot\mat{R}_\beta \cdot \mat{R}_\gamma
\end{equation}
Turnhout, M.C. van's avatar
Turnhout, M.C. van committed
172
This kind of angle multiplication is not commutative. Abaqus uses a different definition for Euler angles. These three Euler angles from  \warning\gibbon{} can not be directly used in Abaqus. It will wreak havoc.
173

174
The $3\times 3$ matrix (centres of mass, Euler angles) is written to \texttt{Bone\_preload\_loading.txt\index{Bone\_preload\_loading.txt@\texttt{Bone\_preload\_loading.txt}}} or \texttt{Bone\_preload\_postload.txt\index{Bone\_preload\_postload.txt@\texttt{Bone\_preload\_postload.txt}}}, depending on the case, in  \basilhome\texttt{/<basilid>/}. To map the tibia movements for \texttt{<basilid> = 140611} and \texttt{<basilid> = 140613} call:
Mark van Turnhout's avatar
Mark van Turnhout committed
175
\begin{lstlisting}[numbers=none,language=basillab]
176
177
>> bas_mapBone([140611, 140613]);
\end{lstlisting}
Mark van Turnhout's avatar
Mark van Turnhout committed
178

179
\subsection{Map indenter movement}\label{mapind}
Mark van Turnhout's avatar
Mark van Turnhout committed
180

181
The indenter movement is of course needed as boundary condition in the FEM. Mapping of the indenter comprises of three steps: (1) The function \bas{estCylPos} finds the 3D orientation and the centre of mass of the indenter with the help from \texttt{MRI\_slices\_loading.mat\index{MRI\_slices\_loading.mat@\texttt{MRI\_slices\_loading.mat}}} and \texttt{contour\_skin\_loading-ncc.mat\index{contour\_skin\_loading-ncc.mat@\texttt{contour\_skin\_loading-ncc.mat}}}; next (2) \bas{estIndDisp}{} uses \texttt{contour\_skin\_preload-ncc.mat\index{contour\_skin\_preload-ncc.mat@\texttt{contour\_skin\_preload-ncc.mat}}} to find the start position of the indenter (for the FEM); and finally (3)  \bas{estIndDisp}{} uses \texttt{contour\_skin\_loading-ncc.mat\index{contour\_skin\_loading-ncc.mat@\texttt{contour\_skin\_loading-ncc.mat}}} to find the end position of the indenter.\\
182

Mark van Turnhout's avatar
Mark van Turnhout committed
183
\noindent First, we find the orientation of the indenter in the MRI-data during indentation. The skin contours during indentation are loaded from \texttt{contour\_skin\_loading-ncc.mat\index{contour\_skin\_loading-ncc.mat@\texttt{contour\_skin\_loading-ncc.mat}}} and used to dismiss voxel data that belongs to the soft tissue in the \texttt{MRI\_slices\_loading.mat\index{MRI\_slices\_loading.mat@\texttt{MRI\_slices\_loading.mat}}}-MR-data. We are left with the very bright indenter and some noise in these images and from these we can threshold the voxels that belong to the indenter. A call to \bas{getIndDir} gives us the eigenvectors of the mass moment of inertia of the indenter voxels, and the centre of mass of these voxels.
Mark van Turnhout's avatar
Mark van Turnhout committed
184
185
186
187
188

We use the eigenvector that is closest to the image $\vec{x}\vec{y}$-plane as the (main) orientation of the indenter. This orientation vector $\vec{v}$ is transformed to a rotation matrix $\mat{R}$ by the \gibbon-script \gib{vecAngle2Rot} such that 
\begin{equation}
\vec{v}^\text{T} = \begin{bmatrix} 0& 0& 1\end{bmatrix}\cdot \mat{R}^\text{T}
\end{equation} 
Mark van Turnhout's avatar
Mark van Turnhout committed
189
A first estimate of the indenter tip position during indentation is calculated from the centre of mass and the indenter orientation. The results are written to disk as a $4\times 3$ matrix in \basilhome\texttt{/<basilid>}\texttt{/Indenter\_tipNR.txt\index{Indenter\_tipNR.txt@\texttt{Indenter\_tipNR.txt}}}. The first row of this matrix contains the estimate of the indenter tip position $[x y z]$ (in mm), and the second to last rows are the $3 \times 3$ rotation matrix $\mat{R}$ from  \gib{vecAngle2Rot}.\\
Mark van Turnhout's avatar
Mark van Turnhout committed
190

Mark van Turnhout's avatar
Mark van Turnhout committed
191
\noindent Next, in \bas{estIndDisp} we use the skin contours before indentation \texttt{contour\_skin\_preload-ncc.mat\index{contour\_skin\_preload-ncc.mat@\texttt{contour\_skin\_preload-ncc.mat}}} and the indenter orientation and tip position from \bas{estCylPos} to find the start position of the indenter for the FEM:  the indenter is moved `backwards' along its (main) orientation until the indenter tip is 1 mm outside the undeformed skin contours. \\
Mark van Turnhout's avatar
Mark van Turnhout committed
192

Mark van Turnhout's avatar
Mark van Turnhout committed
193
\noindent Finally,  we use the skin contours during indentation \texttt{contour\_skin\_loading-ncc.mat\index{contour\_skin\_loading-ncc.mat@\texttt{contour\_skin\_loading-ncc.mat}}} and the indenter orientation and tip position from \bas{estCylPos} to find the end position of the indenter. The indenter is moved along its (main) orientation until the indenter tip starts to penetrate the (deformed) skin contours. This is a minimalisation problem that uses \bas{estIndEnd} as anonymous function to minimize.\\ 
Mark van Turnhout's avatar
Mark van Turnhout committed
194

Mark van Turnhout's avatar
Mark van Turnhout committed
195
\noindent The results of  \bas{estIndDisp} are written to  \basilhome\texttt{/<basilid>}\texttt{/Indenter\_StartEnd.txt\index{Indenter\_StartEnd.txt@\texttt{Indenter\_StartEnd.txt}}} in a $2 \times 3$ matrix. The first row contains the coordinates of the indenter tip  $[x y z]$ (in mm) before indentation, the second row contains the coordinates of the indenter tip during indentation.
196
197

\subsection{Putting it all together: \bas{prepFEM}}
Mark van Turnhout's avatar
Mark van Turnhout committed
198

Mark van Turnhout's avatar
Mark van Turnhout committed
199
200
201
202
203
204
205
206
207
If all went well, you should now have the files listed in table \eqref{prepfemfiles} in your  \basilhome\texttt{/<basilid>/}-directory (in addition to the files listed in table \ref{mritomatout}). The function \bas{prepFEM} checks for the presence of these files and runs the appropriate functions (column two in table \ref{prepfemfiles}) when files are missing.  If the \texttt{contour}-files are missing,  \bas{prepFEM} calls \bas{segmentMRI} to prompt the user to (manually) segment
\begin{enumerate}
\item the skin contours before indentation;
\item the bone contours before indentation;
\item the skin contours during indentation;
\item the bone contours during indentation; and
\item the bone contours after indentation; 
\end{enumerate}
in that order. Other files are produced automagically and do not require user interaction. To prepare your MRI-data for building of the FEM for \texttt{<basilid> = 140611} and \texttt{<basilid> = 140613} call:
Mark van Turnhout's avatar
Mark van Turnhout committed
208
\begin{lstlisting}[numbers=none,language=basillab]
Mark van Turnhout's avatar
Mark van Turnhout committed
209
210
211
212
213
214
215
216
217
218
>> bas_prepFEM([140611, 140613]);
\end{lstlisting}

\begin{table}[b!]
\center
\caption{Files in \basilhome\texttt{/<basilid>/} after processing of MRI-data (in addition to those listed in table \ref{mritomatout}). The function \bas{prepFEM} checks for the existence of these files, and runs the functions in column two when necessary (i.e.\ when files are missing). }\label{prepfemfiles}
\begin{tabulary}{\linewidth}{l l R}
\hline\noalign{\smallskip}
file in \basilhome\texttt{/<basilid>/}& written by & needed for \\  
\noalign{\smallskip}\hline\noalign{\smallskip}
219
\texttt{contour\_skin\_preload.mat\index{contour\_skin\_preload.mat@\texttt{contour\_skin\_preload.mat}}} & \bas{segmentMRI} & \bas{checkVsize}\\
Mark van Turnhout's avatar
Mark van Turnhout committed
220
\texttt{contour\_bone\_preload.mat} & \bas{segmentMRI} & \bas{checkVsize}\\
221
222
223
224
225
226
227
228
229
230
231
\texttt{contour\_skin\_loading.mat\index{contour\_skin\_loading.mat@\texttt{contour\_skin\_loading.mat}}} & \bas{segmentMRI} &  \bas{checkVsize}\\
\texttt{contour\_bone\_loading.mat\index{contour\_bone\_loading.mat@\texttt{contour\_bone\_loading.mat}}} & \bas{segmentMRI} &  \bas{checkVsize}\\
\texttt{contour\_bone\_postload.mat\index{contour\_bone\_postload.mat@\texttt{contour\_bone\_postload.mat}}} & \bas{segmentMRI} &  \bas{checkVsize}\\
\texttt{contour\_skin\_preload-ncc.mat\index{contour\_skin\_preload-ncc.mat@\texttt{contour\_skin\_preload-ncc.mat}}} & \bas{checkVsize} & \bas{estIndDisp}, \bas{buildFEM}\\
\texttt{contour\_bone\_preload-ncc.mat\index{contour\_bone\_preload-ncc.mat@\texttt{contour\_bone\_preload-ncc.mat}}} & \bas{checkVsize} & \bas{mapBone}, \bas{buildFEM}\\
\texttt{contour\_skin\_loading-ncc.mat\index{contour\_skin\_loading-ncc.mat@\texttt{contour\_skin\_loading-ncc.mat}}} & \bas{checkVsize} & \bas{estIndDisp}\\
\texttt{contour\_bone\_loading-ncc.mat\index{contour\_bone\_loading-ncc.mat@\texttt{contour\_bone\_loading-ncc.mat}}} & \bas{checkVsize} & \bas{mapBone}\\
\texttt{contour\_bone\_postload-ncc.mat\index{contour\_bone\_postload-ncc.mat@\texttt{contour\_bone\_postload-ncc.mat}}} & \bas{checkVsize} & \bas{mapBone}\\
\texttt{Bone\_preload\_loading.txt\index{Bone\_preload\_loading.txt@\texttt{Bone\_preload\_loading.txt}}} & \bas{mapBone} & \bas{buildFEM}\\
\texttt{Bone\_preload\_postload.txt\index{Bone\_preload\_postload.txt@\texttt{Bone\_preload\_postload.txt}}} & \bas{mapBone} & post-processing\\
\texttt{Indenter\_tipNR.txt\index{Indenter\_tipNR.txt@\texttt{Indenter\_tipNR.txt}}} & \bas{estCylPos} & \bas{estIndDisp}, \bas{buildFEM}\\
Mark van Turnhout's avatar
Mark van Turnhout committed
232
\texttt{Indenter\_Lc.mat\index{Indenter\_Lc.mat@\texttt{Indenter\_Lc.mat}}} & \bas{estCylPos} & backward compatibility\\
233
\texttt{Indenter\_StartEnd.txt\index{Indenter\_StartEnd.txt@\texttt{Indenter\_StartEnd.txt}}} & \bas{indDisp} & \bas{buildFEM}\\
Mark van Turnhout's avatar
Mark van Turnhout committed
234
235
236
237
\noalign{\smallskip}\hline 
\end{tabulary}
\end{table}

238
\section{Building the FE model}\label{buildmodel}
Mark van Turnhout's avatar
Mark van Turnhout committed
239

Mark van Turnhout's avatar
Mark van Turnhout committed
240
The FEM building process is described in detail in the next chapter. Briefly, we use the script \bas{buildFEM} to read the necessary model data from \basilhome\texttt{/<basilid>/} and to write an Abaqus python script \textsl{in the current working directory}\footnote{Note:  \textsl{in the current working directory}, these files are \textsl{not} written to \basilhome\texttt{/<basilid>/}.}. With this python script, Abaqus builds (and optionally: runs) the FE model: geometry, material parameters, boundary conditions, contact, etc. To write the Abaqus model scripts  for \texttt{<basilid> = 140611} and \texttt{<basilid> = 140613} call:
Mark van Turnhout's avatar
Mark van Turnhout committed
241
\begin{lstlisting}[numbers=none,language=basillab]
Mark van Turnhout's avatar
Mark van Turnhout committed
242
243
>> bas_buildFEM([140611, 140613], param);
\end{lstlisting}
244
where \texttt{param} is a structure with modelling options (see chapter \ref{buildfem}). 
Mark van Turnhout's avatar
Mark van Turnhout committed
245

Mark van Turnhout's avatar
Mark van Turnhout committed
246
The (main) output of \bas{buildFEM} is the python script \texttt{basil<basilid>\_job.py\index{basil<basilid>\_job.py@\texttt{basil<basilid>\_job.py}}}, and the FEM can be build by calling Abaqus with
Mark van Turnhout's avatar
Mark van Turnhout committed
247
248
249
\begin{lstlisting}[numbers=none, language=bash]
mark@telab11:~$ abaqus cae noGUI=basil<basilid>_job.py
\end{lstlisting}
Mark van Turnhout's avatar
Mark van Turnhout committed
250
The (Abaqus) output will be the input-file that Abaqus needs to run the model (again: in the \textsl{working directory}):  \texttt{basil<basilid>.inp\index{basil<basilid>.inp@\texttt{basil<basilid>.inp}}}
Mark van Turnhout's avatar
Mark van Turnhout committed
251

Mark van Turnhout's avatar
Mark van Turnhout committed
252

Mark van Turnhout's avatar
Mark van Turnhout committed
253
\section{Running the FE model}
254

Mark van Turnhout's avatar
Mark van Turnhout committed
255
256
257
258
The \texttt{basil<basilid>.inp\index{basil<basilid>.inp@\texttt{basil<basilid>.inp}}} file can be submitted to Abaqus to run the job\footnote{With \texttt{<jobid>}, you tell Abaqus to write the results of the simulation to \texttt{<jobid>.odb}. In this manual, it is assumed that  \texttt{basil<basilid>.inp\index{basil<basilid>.inp@\texttt{basil<basilid>.inp}}} gets \texttt{<jobid> = basil<basilid>} and that the results are written to  \texttt{basil<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}}.}:
\begin{lstlisting}[numbers=none, language=bash]
mark@telab11:~$ abaqus job=<jobid> input=basil<basilid>.inp
\end{lstlisting}
259

Mark van Turnhout's avatar
Mark van Turnhout committed
260
The Abaqus input file  \texttt{basil<basilid>.inp\index{basil<basilid>.inp@\texttt{basil<basilid>.inp}}} is human readable but visually not very pleasing. With \texttt{basil<basilid>\_job.py\index{basil<basilid>\_job.py@\texttt{basil<basilid>\_job.py}}}, you can ask Abaqus to (also) write  \texttt{basil<basilid>.cae\index{basil<basilid>.cae@\texttt{basil<basilid>.cae}}} which can be opened with the Abaqus pre-processing GUI. Then you can have a look at the model, and run the model from the Abaqus GUI. 
261

262
Abaqus writes the simulation results to a (binary) file \texttt{basil<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}} (\texttt{odb} = output database). 
263

Mark van Turnhout's avatar
Mark van Turnhout committed
264
265
\section{Post-processing of FE results}

Mark van Turnhout's avatar
Mark van Turnhout committed
266
267
BasilLab does not provide any post-processing functionality beyond conversion of a selection of parameters from the (binary)  Abaqus output database  \texttt{basil<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}} to human and Matlab readable text-files; and some basic plotting of the FEM results.

268
269
As with building and running the FE model, post-processing works in the current working directory, not on files in \basilhome\texttt{/<basilid>/}.

Mark van Turnhout's avatar
Mark van Turnhout committed
270
271
\subsection{Extracting data from the \texttt{odb}-file}

272
273
Since \texttt{basil<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}} is only readable for Abaqus, BasilLab provides some Abaqus python scripts to extract simulation data from this file. The data that is read from the \texttt{odb}-file is written to a text-file that can easily be fed to your favorite post-processing software (Matlab). 

274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
Note that only results for the soft tissue in the model are extracted, results for the indenter or the cast are not processed.\\

\noindent The python script \texttt{extract\_topology.py\index{extract\_topology.py@\texttt{extract\_topology.py}}} writes the mesh toplogy for \texttt{basil\-<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}}  to the file  \texttt{basil<basilid>.top\index{basil<basilid>.top@\texttt{basil<basilid>.top}}}. This file contains a row for each element in the soft tissue set with the format: 
\begin{framed}
\texttt{element number}, \texttt{node number 1}, \texttt{node number 2}, \texttt{node number 3}, \dots,  \texttt{node number n}
\end{framed}
The length of each row depends on the number of nodes \texttt{n} for the specific element. \\

\noindent  The python script \texttt{extract\_COORD.py\index{extract\_COORD.py@\texttt{extract\_COORD.py}}} writes the node coordinates for \texttt{basil\-<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}} to the file  \texttt{basil<basilid>.xyz\index{basil<basilid>.xyz@\texttt{basil<basilid>.xyz}}}. This file contains a row for each node at the start of the simulation and at the end of the simulation with the format
\begin{framed}
\texttt{time}, \texttt{node number}, \texttt{x-coordinate}, \texttt{y-coordinate}, \texttt{z-coordinate}
\end{framed}
The \texttt{time} at the start of the simulation is \texttt{0}, the \texttt{time} at the end of the simulation is the \textsl{step} time of the last (indentation) step in the FEM\footnote{With the default settings (see chapter \ref{buildfem}) this step lasts 100\,s. If you find any other \texttt{time} than that step duration, the simulation did not converge.}. Coordinates are in mm.\\

\noindent  The python script \texttt{extract\_EVOL.py\index{extract\_EVOL.py@\texttt{extract\_EVOL.py}}} writes the element volumes for \texttt{basil\-<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}} to the file  \texttt{basil<basilid>.evol\index{basil<basilid>.evol@\texttt{basil<basilid>.evol}}}. This file contains a row for each element at the start of the simulation and at the end of the simulation with the format
\begin{framed}
\texttt{time}, \texttt{element number}, \texttt{element volume}
\end{framed}
The \texttt{time} at the start of the simulation is \texttt{0}, the \texttt{time} at the end of the simulation is the \textsl{step} time of the last (indentation) step in the FEM. Volumes are in mm$^3$.\\

\noindent  The python script \texttt{extract\_ELSE.py\index{extract\_ELSE.py@\texttt{extract\_ELSE.py}}} writes the strain energy per element for \texttt{basil\-<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}} to the file  \texttt{basil<basilid>.sen\index{basil<basilid>.sen@\texttt{basil<basilid>.sen}}}. This file contains a row for each element at the end of the simulation with the format
\begin{framed}
\texttt{time}, \texttt{element number}, \texttt{strain energy}
\end{framed}
The \texttt{time} at the end of the simulation is the \textsl{step} time of the last (indentation) step in the FEM. Strain energies are in J.\\

\noindent  The python script \texttt{extract\_ESEDEN.py\index{extract\_ESEDEN.py@\texttt{extract\_ESEDEN.py}}} writes the strain energy density per element for \texttt{basil\-<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}} to the file  \texttt{basil<basilid>.sed\index{basil<basilid>.sed@\texttt{basil<basilid>.sed}}}. This file contains a row for each element at the end of the simulation with the format
\begin{framed}
\texttt{time}, \texttt{element number}, \texttt{strain energy density}
\end{framed}
The \texttt{time} at the end of the simulation is the \textsl{step} time of the last (indentation) step in the FEM. Strain energies are in J/mm$^3$.\\

\noindent  The python script \texttt{extract\_LE.py\index{extract\_LE.py@\texttt{extract\_LE.py}}} writes the (unique) logarithmic strain tensor components per node for \texttt{basil\-<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}} to the file  \texttt{basil<basilid>.le\index{basil<basilid>.le@\texttt{basil<basilid>.le}}}. This file contains a row for each node at the end of the simulation with the format
\begin{framed}
\texttt{time}, \texttt{node number}, \texttt{LE11}, \texttt{LE22}, \texttt{LE33}, \texttt{LE12}, \texttt{LE13}, \texttt{LE23}
\end{framed}
The \texttt{time} at the end of the simulation is the \textsl{step} time of the last (indentation) step in the FEM. Strains are dimensionless.\\


\noindent To extract data from  \texttt{basil\-<basilid>.odb\index{basil<basilid>.odb@\texttt{basil<basilid>.odb}}} with these Abaqus python scripts, use 
314
\begin{lstlisting}[numbers=none, language=bash]
315
mark@telab11:~$ abaqus python extract_<var>.py -odb basil<basilid>
316
317
\end{lstlisting}

Mark van Turnhout's avatar
Mark van Turnhout committed
318
319

\subsection{Basic plotting of FEM results}