Transcript Introduction to Matlab

Computers in Medicine
Bone visualization and analyses
- An introduction to Matlab -
Presented by: Kathryn Stok
ETH Zurich
[email protected]
Prepared by: Steven Boyd
University of Calgary
Objectives (I)

Provide an overview of Matlab
– What is Matlab?
– Basic commands
– Accessing online help
– Writing M-files
– Programming style
– Create plots
Objectives (II)

Work through examples
– Reading an image file
– Probe the data
– Filtering
– Create a histogram
– Thresholding
– Counting voxels
– Calculating 2nd moment of area
Why bone visualization?

Carries the body
→ it should be strong

But in many cases, bone
strength is impaired.
– congenital deformities
– osteoporosis
– fracture
– fatigue (stress fracture)
Bone architecture
Bone structure – femoral head

Young Normal

Osteoporotic
Bone structure – lumbar spine

Young Normal

Osteoporotic
Why bone visualization?

Bone (structure) analyses and visualization are of
great importance in
– Diagnosis
– Treatment
– Prevention
Why Matlab?

High-performance language for technical computing.
– Matrix Laboratory

Advantages
– Easy to use:
– No dimensioning required.
– No compiling
– Good for prototyping, testing
– Functionality is greatly expanded by toolboxes.
– Many fields of application.

Disadvantages:
– Uses much memory (slow for large applications)
– Not good for final software design
What is matlab?

Matlab language
– matrix/array language
– control flow (if, for, while)
– functions

Working environment
– workspace for computing, importing, and
exporting data
– text editor for creating M-files
This Matlab course

It covers:
• Using the development environment.
• Syntax of the language.
• Basic matrix operations.
• Basic graphics features.
• Writing Matlab scripts and functions (m-files)

It does not cover:
• Any of the toolboxes including Simulink.
• c-mex and f-mex files, Matlab compilers.
Possibilities…

Develop tools for:
– Image processing and basic analysis.
– Visualization: 2D and 3D.
– Simulation: Simulate bone loss.
Matlab Environment
Matrices

Magic square
Albrecht Dürer
Matrices

Entering matrices
– explicitly
semicolon
Matrices

Sum, transpose
Matrices

Diagonal (diag)
Matrices

Subscripts: A(i,j)
row
column
Matrices

Colon operator ' : '
Matrices

Colon operator ' : '
Matlab is 1-based
Matrix multiplication ( * )
Matrix multiplication ( .* )
Variables

Scalar, vector, matrix
Variables

Check workspace
– command line
– window
Numbers

Integer
-3

Real
0.0001

Exponential
1.60210e-20

Imaginary
5 + 2i
Operators: Arithmetic

+
*
/
^

'

()




Addition
Subtraction
Multiplication (or .* element-wise)
Division
Power (or .^ element-wise)
Transpose
Brackets for evaluation order
A = [ 4 1
A.^2 = [16
A.*2 = [ 8
3]
1 9]
2 6]
 16 4 12 


A'* A   4 1 3 
 12 3 9 
Operators: Relational

<
Less than

<= Less than or equal to

>

>= Greater than or equal to

== Equal to

~= Not equal to
Greater than
Operators: Logical

&
AND

|
OR

~
NOT
Type help ops to see:
the arithmetic,relational and logical functions.
Matrices

Entering matrices
– explicitly
– load from external file
– your own M-files
– built-in function
Albrecht Dürer
Generating matrices

Load
data file
Generating matrices

M-files
Generating matrices

Built-in functions
Help

Help Window
– pop-up window

Type help
Online help
Online help
Online help
Useful commands

(up arrow)
– recall previous command
– p recalls previous command starting with ‘p’

clear
– clear workspace
– clear A Z Q clears only A, Z, and Q.

size
– size(A) returns size of matrix A
Multidimensional matrices
column j, Y
page k,
Z
row
i, X
Multidimensional matrices

Generate a 3D array manually
Advanced indexing

Matrices are always stored as
columns
– subscripts (i,j,k)
– array dimension [d1 d2 d3]
i
j
3
2
7
k
–
if A is of dimension 3 x 3 x 3
then A(3,3,3) = A(27) = 8
5 4
1 6
1 8
6 1 3
2 4 6
6 9 7
5 9 1
6 5 2
7 6 8
3
2
7
5
1
1
4
6
8
6
2
6
1
4
9
3
6
7
5
6
7
9
5
6
1
2
8
Characters and Text

s = 'myimage'
– 7 character array (string)

num2str() - convert a number to string
User input

Keyboard input
a_number = input('Enter number');
a_string = input('Enter filename', 's');

Using GUI
[filename,path] = uigetfile('*.*');
Graphics: plot
Graphics: plot(s)
‘hold on’
allows more plots
to be added
Graphics: contours
Graphics: subplots
Graphics: labels and titles
Graphics: labels and titles
… or use the Figure window !
Graphics: mesh
Flow control

if, else, and elseif

switch

while

for
if, else, and elseif
if logical_expression
statements
end
if (rem(a,2) == 0)
disp(‘a is even’)
end
if (n<0)
disp(‘Input must be positive’)
elseif (rem(n,2) == 0)
disp(‘a is even’)
else
disp(‘a is odd’)
end
switch
switch expression
case value1
statements
case value2
statements
otherwise
statements
end
switch (input_num)
case (-1)
disp(‘negative one’)
case (0)
disp(‘zero’)
case (1)
disp(‘positive one’)
otherwise
disp(‘othervalue’)
end
while, for
while expression
for index = start:inc:end
end
end
n = 1;
while (prod(1:n) < 1e10)
n = n + 1;
end
for i = 2:6
x(i) = 2*x(i-1);
end
statements
statements
for i = 1:m
for j = 1:n
A(i,j) = 1/(i + j - 1);
end
end
Vectorization

Efficiency!
– Example: Table of logarithms
– Find log10 of numbers between 0 and 10.
x = 0;
for k = 1:1001
y(k) = log10(x);
x = x + 0.01;
end
x = 0:0.01:10;
y = log10(x);
Looped code
Vectorized code
AVOID LOOPS!
Scripts and functions

Scripts
–
No input arguments nor
output arguments.
–
Operate on data in the
workspace
–
Useful for automating a
series of steps that will
be performed many
times

Functions
–
Can accept input
arguments and return
output arguments
–
Variables are local to
the function
–
Useful for extending the
MATLAB language
beyond the built-in
functions
Script example

Create m-files:
– type edit
– use the interface
Comment
lines
Computations
Graphics
output
Script example
Hit return to advance plot
Script Example
Function example
Function
definition
Help
comments
Error check
Calculations
function [y] = average(x)
input argument
function name
output argument
keyword
Function example
semi-colon
M-files

Tips:
– use descriptive variable names
– e.g. number_students = 20
– functions end with “_fcn”
– e.g. an_example_fcn(number_students)
– indent if, for, while statements 3-4 spaces
– add comments using ‘%’
– develop with smaller portions of data
– debug by pasting line-by-line from M-file
Online tutorial
http://www.imrtweb.ethz.ch/matlab/
Or google “mathworks tutorial”
Examples
Exercises
Fundamentals
Part II: Practical Programming

Image analysis on 2D section
– Read image data (micro-CT)
– Plot the image
– Examine image
– Filter the image
– Create a histogram
– Threshold (distinguish between bone and marrow)
– Calculate porosity
– Centroid, principal axes
Read 2D Image

Want to import the image from micro-CT

Write a program to:
– clear workspace
– set filename
– load file
–
Image data to visualize:
/Phalanx/hp2d_34-95_380.bmp
Set up main program (a script file)
comments
clear
workspace
close all
figure windows
file to
load
load it
Note: You need to write the 2D image reading function.
Read in 2D image
Read the file
Convert data
Display input data
Display input data
Display input data
Examine image

min, max, mean
Examine image

min, max, mean
Image Processing

Filter

Thresholding
– histogram

Porosity = 1 - (bone area / total area)

Moment of area
– First moment of area
– Centroid
– Second moment of area
Filtering – example 1
Filtering – example 2
original image
filtered original image
noisy image
filtered noisy image
Filtering the data

Apply Gaussian filtration
–
Shape of filter depends on:
window size [n1 n2] and s
Filtering
Filter
Window
s
=3
= 1.2
Filtering
Filter
Window
s
=5
= 1.2
Filtering
Filter
Window
s
= 15
= 1.2
Filtering
Filter
Window
s
= 15
=7
Filtering
Filter
Window
s
= 15
= 100
Filtering: Convolution
fx  f y  fz
Filtering: Edge Effects
edge data nulled
• Important to begin
with image dimensions
large enough to
account for edge
effects
Filtering

Use built-in functions:
– fspecial
– imfilter
– (smooth3)
Thresholding

Segmentation of bone
– divides bone from other tissue
• Good to try different thresholds to see the effect.
Thresholding
threshold2D = zeros matrix
Find all indices (i,j) where element(i,j) > threshold
threshold2D(i,j) = 1
Thresholding
Thresholding
Threshold = 90
Threshold = 120
Histogram
–
lookfor and help
Histogram

hist works on vectors, not 2D matrices
– need reshape function to vectorize
Porosity

Porosity = 1 - (Bone volume divided by total volume)

Need total bone volume (area)
– Create mask:
– filter image to remove holes
– threshold filtered image
– Count total “mask” pixels
– Count total “bone” pixels
within mask
Porosity

Create mask:
original
image
filtered
image
threshold
filtered image
1
2
3
Image 3 super-imposed
on image 1
Porosity
Porosity = 1 -
Bone pixels within mask
Mask pixels
Thresholded
image (bone)
Image of bone
within the mask
Mask image
Image analysis – determine porosity
Segmenting the bone
Filtered
Thresholded
Find bone within mask
Original
Creating the mask
Filtered
Thresholded
Porosity = 1 -
Areabone
Areamask
Erosion, dilation
Erode
Erode = 2
Erode = 4
Dilate
Dilate = 2
Dilate = 4
Part of MAIN program
Stress in Pure Bending

Mechanical stiffness
– First moment of area
– Centroid of area
– Second moment of area
First moment of an area

First moment of area with respect to the X axis:
(mm3)

First moment of area with respect to the Y axis:
(mm3)
Can be positive, negative, or zero
depending on the axis location
Centroid of an area

Centroid of area defined as:

Knowing the first moment:
First moment of area
Second moment of area

Stress resulting from bending: σ = My/I

Resistance against bending.
– Take Ix, Iy, Ixy about centroid C
– Or use the parallel axis theorem: Ix = Ixc + Ad2
Second moment of area - rotation

Transformation of 2nd moments:
Second moment of area - principal axes

Two values of q for which Ixy = 0.

Can calculate the max/min 2nd moments:
maximum associated
with Ix
minimum associated
with Iy
Principal axes & Transformation
2D: Pompe Disease
Radius
Tibia
Mouse knees
CTRL
NTRT