Newer
Older
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# 3D structure tensor - a small example\n",
"A small example demonstrating the use of 3D structure tensor for visualizing and clustering dominant orientation. "
]
},
{
"cell_type": "code",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib notebook\n",
"import numpy as np\n",
"import scipy.io\n",
"import st3d"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reading the data\n",
"\n",
"The data is a small cube from a volumetric image og fibre composite. The cube contains five bundles in layers: UD fibre (0 deg), crossing fibre (45 deg), 90 deg fibre, -45 deg bundle, UD bundle (0 deg)."
]
},
{
"cell_type": "code",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(150, 150, 150)\n"
]
}
],
"source": [
"volume = scipy.io.loadmat('example_data_3D/multi_cube.mat')['vol']\n",
"print(volume.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Computing the structure tensor and the dominant orientation\n",
"Computation of structure tensor requires only two parameters: the noise scale sigma and the integration scale rho. Parameter sigma controls smothing while computing gradientes, and structures smaller than sigma will be removed by smoothing. Parameter rho gives the size over the neighborhood in which the orientation is to be analysed for every volume voxel. Larger rho will result in a smoother orientation field.\n",
"\n",
"Structure tensor is a 3x3 matrix, but as it is symmetrical we only carry values of 6 elements: $s_{xx}$, $s_{yy}$, $s_{zz}$, $s_{xy}$, $s_{xz}$ and $s_{yz}$.\n",
"\n",
"Eigenvalues (val) carry the information about the degree of anisotropy - this is not used or visualized here. Eigenvectors (vec) carry the orientation information, as $x$, $y$, and $z$ component of the orientation vector."
]
},
{
"cell_type": "code",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"The volume has a shape (150, 150, 150), i.e. 3375000 voxels.\n",
"Structure tensor information is carried in a (6, 3375000) array.\n",
"Orientation information is carried in a (3, 3375000) array.\n"
]
}
],
"source": [
"sigma = 0.5; # noise scale\n",
"rho = 2; # integration scale\n",
"S = st3d.structure_tensor(volume, sigma, rho)\n",
"val,vec = st3d.eig_special(S)\n",
"print(f'The volume has a shape {volume.shape}, i.e. {volume.size} voxels.')\n",
"print(f'Structure tensor information is carried in a {S.shape} array.')\n",
"print(f'Orientation information is carried in a {vec.shape} array.')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualizing the dominant orientation\n",
"Here we show only dominant orientation, ignoring shape measures. \n",
"\n",
"Since slicng is in $z$ direction, the arrows show $x$ and $y$ component of orientation vectors, on top of every slice.\n",
"\n",
"To interactively inspect the volume use arrow keys."
]
},
{
"cell_type": "code",
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
"metadata": {},
"outputs": [
{
"data": {
"application/javascript": [
"/* Put everything inside the global mpl namespace */\n",
"window.mpl = {};\n",
"\n",
"\n",
"mpl.get_websocket_type = function() {\n",
" if (typeof(WebSocket) !== 'undefined') {\n",
" return WebSocket;\n",
" } else if (typeof(MozWebSocket) !== 'undefined') {\n",
" return MozWebSocket;\n",
" } else {\n",
" alert('Your browser does not have WebSocket support. ' +\n",
" 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
" 'Firefox 4 and 5 are also supported but you ' +\n",
" 'have to enable WebSockets in about:config.');\n",
" };\n",
"}\n",
"\n",
"mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
" this.id = figure_id;\n",
"\n",
" this.ws = websocket;\n",
"\n",
" this.supports_binary = (this.ws.binaryType != undefined);\n",
"\n",
" if (!this.supports_binary) {\n",
" var warnings = document.getElementById(\"mpl-warnings\");\n",
" if (warnings) {\n",
" warnings.style.display = 'block';\n",
" warnings.textContent = (\n",
" \"This browser does not support binary websocket messages. \" +\n",
" \"Performance may be slow.\");\n",
" }\n",
" }\n",
"\n",
" this.imageObj = new Image();\n",
"\n",
" this.context = undefined;\n",
" this.message = undefined;\n",
" this.canvas = undefined;\n",
" this.rubberband_canvas = undefined;\n",
" this.rubberband_context = undefined;\n",
" this.format_dropdown = undefined;\n",
"\n",
" this.image_mode = 'full';\n",
"\n",
" this.root = $('<div/>');\n",
" this._root_extra_style(this.root)\n",
" this.root.attr('style', 'display: inline-block');\n",
"\n",
" $(parent_element).append(this.root);\n",
"\n",
" this._init_header(this);\n",
" this._init_canvas(this);\n",
" this._init_toolbar(this);\n",
"\n",
" var fig = this;\n",
"\n",
" this.waiting = false;\n",
"\n",
" this.ws.onopen = function () {\n",
" fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n",
" fig.send_message(\"send_image_mode\", {});\n",
" if (mpl.ratio != 1) {\n",
" fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n",
" }\n",
" fig.send_message(\"refresh\", {});\n",
" }\n",
"\n",
" this.imageObj.onload = function() {\n",
" if (fig.image_mode == 'full') {\n",
" // Full images could contain transparency (where diff images\n",
" // almost always do), so we need to clear the canvas so that\n",
" // there is no ghosting.\n",
" fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n",
" }\n",
" fig.context.drawImage(fig.imageObj, 0, 0);\n",
" };\n",
"\n",
" this.imageObj.onunload = function() {\n",
" fig.ws.close();\n",
" }\n",
"\n",
" this.ws.onmessage = this._make_on_message_function(this);\n",
"\n",
" this.ondownload = ondownload;\n",
"}\n",
"\n",
"mpl.figure.prototype._init_header = function() {\n",
" var titlebar = $(\n",
" '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n",
" 'ui-helper-clearfix\"/>');\n",
" var titletext = $(\n",
Loading
Loading full blame...