workflow.tex 33.9 KB
 Mark van Turnhout committed Dec 11, 2015 1 2 % !TeX root = BasilLab.tex \begin{savequote}  Mark van Turnhout committed Dec 30, 2015 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}}  Mark van Turnhout committed Dec 11, 2015 4 \end{savequote}  Mark van Turnhout committed Dec 21, 2015 5 \chapter{Work flow}  Mark van Turnhout committed Dec 11, 2015 6   Mark van Turnhout committed Dec 21, 2015 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 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{} which is created in the directory \texttt{} (\texttt{bas\_setDataRoot}). The link between the raw MRI-data folder inside \texttt{} and \texttt{} is \textsl{hard-coded} in \texttt{bas\_getMeta.m}.  Mark van Turnhout committed Dec 11, 2015 12   Mark van Turnhout committed Dec 21, 2015 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 committed Dec 11, 2015 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 committed Dec 21, 2015 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{} 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{/}.\\ \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 committed Jan 03, 2016 34 \begin{lstlisting}[numbers=none, language=basillab]  Mark van Turnhout committed Dec 21, 2015 35 36 37 38 >> bas_setMRIroot('/home/mark/BrukerData'); >> bas_setDataRoot(pwd); \end{lstlisting}  Mark van Turnhout committed Dec 13, 2015 39 40 \section{Pre-processing of MRI-data}  Mark van Turnhout committed Dec 21, 2015 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 committed Dec 21, 2015 43 \item Copy MRI-data from \MRIroot{}\texttt{/} to \basilhome\texttt{/}  Mark van Turnhout committed Dec 21, 2015 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 committed Dec 21, 2015 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).\\  Mark van Turnhout committed Dec 22, 2015 56 \noindent The first argument of \bas{MRItoMat} is a list of \texttt{}'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{/}. When the list of \texttt{}'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{/} is not processed; when it is \texttt{1}, all data present in \basilhome\texttt{/} 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 committed Dec 21, 2015 57 58  To process the data for \texttt{ = 140611} and \texttt{ = 140613} and skip processing of data that is already present in \basilhome\texttt{/}, call:  Mark van Turnhout committed Jan 03, 2016 59 \begin{lstlisting}[numbers=none, language=basillab]  Mark van Turnhout committed Dec 21, 2015 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{/} 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 committed Dec 22, 2015 64 \begin{table}[p!]  Mark van Turnhout committed Dec 21, 2015 65 \center  Mark van Turnhout committed Dec 22, 2015 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 committed Dec 21, 2015 67 68 \begin{tabulary}{\linewidth}{l R} \hline\noalign{\smallskip}  Mark van Turnhout committed Dec 22, 2015 69 70 file in \basilhome\texttt{//}& description\\ \noalign{\smallskip}\hline\noalign{\smallskip}  Mark van Turnhout committed Dec 27, 2015 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\_.mat@\texttt{IMDAT\_.mat}}} & as \texttt{MRI\_slices\_preload.mat\index{MRI\_slices\_preload.mat@\texttt{MRI\_slices\_preload.mat}}} \\ \texttt{IMDAT\_loading.mat\index{IMDAT\_.mat@\texttt{IMDAT\_.mat}}} & as \texttt{MRI\_slices\_loading.mat\index{MRI\_slices\_loading.mat@\texttt{MRI\_slices\_loading.mat}}} \\ \texttt{IMDAT\_postload.mat\index{IMDAT\_.mat@\texttt{IMDAT\_.mat}}} & as \texttt{MRI\_slices\_postload.mat\index{MRI\_slices\_postload.mat@\texttt{MRI\_slices\_postload.mat}}} \\ \texttt{IMDAT\_flash.mat\index{IMDAT\_.mat@\texttt{IMDAT\_.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\_.mat@\texttt{MRI\_subject\_.mat}}}\\  Mark van Turnhout committed Dec 22, 2015 80 81 82 \texttt{MRI\_subject\_loading.mat} & \\ \texttt{MRI\_subject\_postload.mat} &\\ \texttt{MRI\_subject\_flash.mat} &\\\hline\noalign{\smallskip}  Mark van Turnhout committed Dec 27, 2015 83 \texttt{MRI\_acqp\_preload.mat} & \multirow{4}*{acqp' (aquisition parameters) from MRI meta\index{MRI\_acqp\_.mat@\texttt{MRI\_acqp\_.mat}}}\\  Mark van Turnhout committed Dec 22, 2015 84 85 86 \texttt{MRI\_acqp\_loading.mat} & \\ \texttt{MRI\_acqp\_postload.mat} & \\ \texttt{MRI\_acqp\_flash.mat} & \\ \hline  Mark van Turnhout committed Dec 27, 2015 87 \texttt{MRI\_method\_preload.mat}&\multirow{4}*{method' from MRI meta data\index{MRI\_method\_.mat@\texttt{MRI\_method\_.mat}}} \\  Mark van Turnhout committed Dec 22, 2015 88 89 90 \texttt{MRI\_method\_loading.mat} & \\ \texttt{MRI\_method\_postload.mat} & \\ \texttt{MRI\_method\_flash.mat} & \\ \hline\noalign{\smallskip}  Mark van Turnhout committed Dec 27, 2015 91 \texttt{MRI\_imnd\_preload.mat} & \multirow{4}*{imnd' from MRI meta data\index{MRI\_imnd\_.mat@\texttt{MRI\_imnd\_.mat}}} \\  Mark van Turnhout committed Dec 22, 2015 92 93 94 \texttt{MRI\_imnd\_loading.mat} & \\ \texttt{MRI\_imnd\_postload.mat} & \\ \texttt{MRI\_imnd\_flash.mat} & \\ \hline\noalign{\smallskip}  Mark van Turnhout committed Dec 27, 2015 95 \texttt{MRI\_reco\_preload.mat} & \multirow{4}*{reco' from MRI meta data\index{MRI\_reco\_.mat@\texttt{MRI\_reco\_.mat}}} \\  Mark van Turnhout committed Dec 22, 2015 96 97 98 \texttt{MRI\_reco\_loading.mat} & \\ \texttt{MRI\_reco\_postload.mat} & \\ \texttt{MRI\_reco\_flash.mat} & \\ \hline\noalign{\smallskip}  Mark van Turnhout committed Dec 27, 2015 99 \texttt{MRI\_d3proc\_preload.mat} & \multirow{4}*{d3proc' from MRI meta data\index{MRI\_d3proc\_.mat@\texttt{MRI\_d3proc\_.mat}}} \\  Mark van Turnhout committed Dec 22, 2015 100 101 102 \texttt{MRI\_d3proc\_loading.mat} & \\ \texttt{MRI\_d3proc\_postload.mat} & \\ \texttt{MRI\_d3proc\_flash.mat} & \\ \hline\noalign{\smallskip}  Mark van Turnhout committed Dec 27, 2015 103 \texttt{MRI\_visupars\_preload.mat} & \multirow{4}*{visupars' (visualisation parameters) from MRI meta data\index{MRI\_visupars\_.mat@\texttt{MRI\_visupars\_.mat}}} \\  Mark van Turnhout committed Dec 22, 2015 104 105 106 \texttt{MRI\_visupars\_loading.mat} & \\ \texttt{MRI\_visupars\_postload.mat} & \\ \texttt{MRI\_visupars\_flash.mat} & \\\hline\noalign{\smallskip}  Mark van Turnhout committed Dec 27, 2015 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\_.mat@\texttt{MRI\_voxel\_.mat}}}\\  Mark van Turnhout committed Dec 22, 2015 108 109 110 \textbf{MRI\_voxel\_loading.txt} & \\ \textsl{MRI\_voxel\_postload.txt} & \\ \texttt{MRI\_voxel\_flash.txt} & \\  Mark van Turnhout committed Dec 21, 2015 111 112 113 \noalign{\smallskip}\hline \end{tabulary} \end{table}  Mark van Turnhout committed Dec 21, 2015 114 115 116  \subsection{Segment MRI-data}  Mark van Turnhout committed Dec 11, 2015 117   Mark van Turnhout committed Dec 22, 2015 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 committed Dec 26, 2015 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.  Mark van Turnhout committed Dec 22, 2015 136 137 138  BasilLab is very lazy on the masking part.  Mark van Turnhout committed Dec 26, 2015 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.  Mark van Turnhout committed Dec 22, 2015 140   Mark van Turnhout committed Dec 27, 2015 141 The mask for \basilhome\texttt{//}\texttt{MRI\_slices\_.mat\index{MRI\_slices\_.mat@\texttt{MRI\_slices\_.mat}}}\footnote{\texttt{} can be \texttt{preload} or \texttt{loading} or \texttt{postload} or \texttt{flash} (see table \ref{mritomatout}).}will be saved as \basilhome\texttt{/}\-\texttt{}\-\texttt{/MRI\_BG\_.mat\index{MRI\_BG\_.mat@\texttt{MRI\_BG\_.mat}}}. To mask the data in the scans before (\texttt{1}), during (\texttt{2}) and after indentation (\texttt{3}) for \texttt{ = 140611} call:  Mark van Turnhout committed Jan 03, 2016 142 \begin{lstlisting}[numbers=none, language=basillab]  Mark van Turnhout committed Dec 22, 2015 143 144 145 >> bas_maskBG(140611, [1 2 3]); \end{lstlisting}  Mark van Turnhout committed Dec 27, 2015 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{//}\texttt{MRI\_slices\_.mat\index{MRI\_slices\_.mat@\texttt{MRI\_slices\_.mat}}} will be saved as \basilhome\texttt{//}\texttt{contour\_skin\_.mat\index{contour\_skin\_.mat@\texttt{contour\_skin\_.mat}}}, that of the tibia as \basilhome\texttt{/}\-\texttt{/contour\_bone\_.mat\index{contour\_bone\_.mat@\texttt{contour\_bone\_.mat}}}.  Mark van Turnhout committed Dec 22, 2015 147   Mark van Turnhout committed Dec 27, 2015 148 Note that \bas{segmentMRI} loaded \basilhome\texttt{//}\texttt{MRI\_voxel\_.txt\index{MRI\_voxel\_.mat@\texttt{MRI\_voxel\_.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{ = 140611} call:  Mark van Turnhout committed Jan 03, 2016 149 \begin{lstlisting}[numbers=none,language=basillab]  Mark van Turnhout committed Dec 22, 2015 150 151 152 >> bas_segmentMRI(140611, [1:4, 6]); \end{lstlisting}  Mark van Turnhout committed Dec 27, 2015 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{//}, and deletes a contour that is missing from \textsl{any} set from \textsl{all} sets.  Mark van Turnhout committed Dec 22, 2015 154   Mark van Turnhout committed Dec 27, 2015 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{ = 140611} and \texttt{ = 140613} call:  Mark van Turnhout committed Jan 03, 2016 156 \begin{lstlisting}[numbers=none,language=basillab]  Mark van Turnhout committed Dec 22, 2015 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\_-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.  Mark van Turnhout committed Dec 22, 2015 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}.  Mark van Turnhout committed Dec 22, 2015 167   Mark van Turnhout committed Dec 26, 2015 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:  Mark van Turnhout committed Dec 22, 2015 169 170 171  \mat{R}_{\alpha\beta\gamma} = \mat{R}_\alpha\cdot\mat{R}_\beta \cdot \mat{R}_\gamma  Turnhout, M.C. van committed Jul 08, 2017 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.  Mark van Turnhout committed Dec 22, 2015 173   Mark van Turnhout committed Dec 27, 2015 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{//}. To map the tibia movements for \texttt{ = 140611} and \texttt{ = 140613} call:  Mark van Turnhout committed Jan 03, 2016 175 \begin{lstlisting}[numbers=none,language=basillab]  Mark van Turnhout committed Dec 22, 2015 176 177 >> bas_mapBone([140611, 140613]); \end{lstlisting}  Mark van Turnhout committed Dec 22, 2015 178   Mark van Turnhout committed Dec 22, 2015 179 \subsection{Map indenter movement}\label{mapind}  Mark van Turnhout committed Dec 11, 2015 180   Mark van Turnhout committed Jan 05, 2016 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.\\  Mark van Turnhout committed Dec 22, 2015 182   Mark van Turnhout committed Dec 27, 2015 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 committed Dec 26, 2015 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 \vec{v}^\text{T} = \begin{bmatrix} 0& 0& 1\end{bmatrix}\cdot \mat{R}^\text{T}  Mark van Turnhout committed Dec 27, 2015 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{/}\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 committed Dec 26, 2015 190   Mark van Turnhout committed Dec 27, 2015 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 committed Dec 26, 2015 192   Mark van Turnhout committed Dec 27, 2015 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 committed Dec 26, 2015 194   Mark van Turnhout committed Dec 27, 2015 195 \noindent The results of \bas{estIndDisp} are written to \basilhome\texttt{/}\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.  Mark van Turnhout committed Dec 22, 2015 196 197  \subsection{Putting it all together: \bas{prepFEM}}  Mark van Turnhout committed Dec 11, 2015 198   Mark van Turnhout committed Dec 26, 2015 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{//}-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{ = 140611} and \texttt{ = 140613} call:  Mark van Turnhout committed Jan 03, 2016 208 \begin{lstlisting}[numbers=none,language=basillab]  Mark van Turnhout committed Dec 26, 2015 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{//} 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{//}& written by & needed for \\ \noalign{\smallskip}\hline\noalign{\smallskip}  Mark van Turnhout committed Dec 27, 2015 219 \texttt{contour\_skin\_preload.mat\index{contour\_skin\_preload.mat@\texttt{contour\_skin\_preload.mat}}} & \bas{segmentMRI} & \bas{checkVsize}\\  Mark van Turnhout committed Dec 26, 2015 220 \texttt{contour\_bone\_preload.mat} & \bas{segmentMRI} & \bas{checkVsize}\\  Mark van Turnhout committed Dec 27, 2015 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 committed Dec 27, 2015 232 \texttt{Indenter\_Lc.mat\index{Indenter\_Lc.mat@\texttt{Indenter\_Lc.mat}}} & \bas{estCylPos} & backward compatibility\\  Mark van Turnhout committed Dec 27, 2015 233 \texttt{Indenter\_StartEnd.txt\index{Indenter\_StartEnd.txt@\texttt{Indenter\_StartEnd.txt}}} & \bas{indDisp} & \bas{buildFEM}\\  Mark van Turnhout committed Dec 26, 2015 234 235 236 237 \noalign{\smallskip}\hline \end{tabulary} \end{table}  Mark van Turnhout committed Dec 22, 2015 238 \section{Building the FE model}\label{buildmodel}  Mark van Turnhout committed Dec 11, 2015 239   Mark van Turnhout committed Dec 26, 2015 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{//} 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{//}.}. 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{ = 140611} and \texttt{ = 140613} call:  Mark van Turnhout committed Jan 03, 2016 241 \begin{lstlisting}[numbers=none,language=basillab]  Mark van Turnhout committed Dec 26, 2015 242 243 >> bas_buildFEM([140611, 140613], param); \end{lstlisting}  Mark van Turnhout committed Dec 29, 2015 244 where \texttt{param} is a structure with modelling options (see chapter \ref{buildfem}).  Mark van Turnhout committed Dec 26, 2015 245   Mark van Turnhout committed Dec 27, 2015 246 The (main) output of \bas{buildFEM} is the python script \texttt{basil\_job.py\index{basil\_job.py@\texttt{basil\_job.py}}}, and the FEM can be build by calling Abaqus with  Mark van Turnhout committed Dec 26, 2015 247 248 249 \begin{lstlisting}[numbers=none, language=bash] mark@telab11:~$abaqus cae noGUI=basil_job.py \end{lstlisting}  Mark van Turnhout committed Dec 27, 2015 250 The (Abaqus) output will be the input-file that Abaqus needs to run the model (again: in the \textsl{working directory}): \texttt{basil.inp\index{basil.inp@\texttt{basil.inp}}}  Mark van Turnhout committed Dec 26, 2015 251   Mark van Turnhout committed Dec 11, 2015 252   Mark van Turnhout committed Dec 27, 2015 253 \section{Running the FE model}  Mark van Turnhout committed Dec 27, 2015 254   Mark van Turnhout committed Dec 27, 2015 255 256 257 258 The \texttt{basil.inp\index{basil.inp@\texttt{basil.inp}}} file can be submitted to Abaqus to run the job\footnote{With \texttt{}, you tell Abaqus to write the results of the simulation to \texttt{.odb}. In this manual, it is assumed that \texttt{basil.inp\index{basil.inp@\texttt{basil.inp}}} gets \texttt{ = basil} and that the results are written to \texttt{basil.odb\index{basil.odb@\texttt{basil.odb}}}.}: \begin{lstlisting}[numbers=none, language=bash] mark@telab11:~$ abaqus job= input=basil.inp \end{lstlisting}  Mark van Turnhout committed Dec 27, 2015 259   Mark van Turnhout committed Dec 27, 2015 260 The Abaqus input file \texttt{basil.inp\index{basil.inp@\texttt{basil.inp}}} is human readable but visually not very pleasing. With \texttt{basil\_job.py\index{basil\_job.py@\texttt{basil\_job.py}}}, you can ask Abaqus to (also) write \texttt{basil.cae\index{basil.cae@\texttt{basil.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.  Mark van Turnhout committed Dec 27, 2015 261   Mark van Turnhout committed Dec 28, 2015 262 Abaqus writes the simulation results to a (binary) file \texttt{basil.odb\index{basil.odb@\texttt{basil.odb}}} (\texttt{odb} = output database).  Mark van Turnhout committed Dec 27, 2015 263   Mark van Turnhout committed Dec 11, 2015 264 265 \section{Post-processing of FE results}  Mark van Turnhout committed Dec 27, 2015 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.odb\index{basil.odb@\texttt{basil.odb}}} to human and Matlab readable text-files; and some basic plotting of the FEM results.  Mark van Turnhout committed Dec 28, 2015 268 269 As with building and running the FE model, post-processing works in the current working directory, not on files in \basilhome\texttt{//}.  Mark van Turnhout committed Dec 27, 2015 270 271 \subsection{Extracting data from the \texttt{odb}-file}  Mark van Turnhout committed Dec 28, 2015 272 273 Since \texttt{basil.odb\index{basil.odb@\texttt{basil.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).  Mark van Turnhout committed Dec 28, 2015 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\-.odb\index{basil.odb@\texttt{basil.odb}}} to the file \texttt{basil.top\index{basil.top@\texttt{basil.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\-.odb\index{basil.odb@\texttt{basil.odb}}} to the file \texttt{basil.xyz\index{basil.xyz@\texttt{basil.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\-.odb\index{basil.odb@\texttt{basil.odb}}} to the file \texttt{basil.evol\index{basil.evol@\texttt{basil.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\-.odb\index{basil.odb@\texttt{basil.odb}}} to the file \texttt{basil.sen\index{basil.sen@\texttt{basil.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\-.odb\index{basil.odb@\texttt{basil.odb}}} to the file \texttt{basil.sed\index{basil.sed@\texttt{basil.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\-.odb\index{basil.odb@\texttt{basil.odb}}} to the file \texttt{basil.le\index{basil.le@\texttt{basil.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\-.odb\index{basil.odb@\texttt{basil.odb}}} with these Abaqus python scripts, use  Mark van Turnhout committed Dec 28, 2015 314 \begin{lstlisting}[numbers=none, language=bash]  Mark van Turnhout committed Dec 28, 2015 315 mark@telab11:~\$ abaqus python extract_.py -odb basil  Mark van Turnhout committed Dec 28, 2015 316 317 \end{lstlisting}  Mark van Turnhout committed Dec 27, 2015 318 319  \subsection{Basic plotting of FEM results}