An Overview of MATLAB
focusing on image processing
19/10/2008
rsan Aytekin
Matlab Environment
Every data is kept as matrices or arrays. Operations are based on matrix and element-wise operations.
19/10/2008
rsan Aytekin
Generating Matrices
Four functions generating basic matrices.
zeros ones rand randn [...] All zeros All ones Uniformly distributed random elements Normally distributed random elements Matrix operator
19/10/2008
rsan Aytekin
Generating Matrices
Here are some examples:
Z = zeros(2,4) Z = 0 0 0 0 F = 5*ones(3,3) F = 5 5 5 5 5 5
0 0
0 0
5 5 5
D = [1 3 5; 2 4 6] D = 1 3 5 2 4 6
19/10/2008 rsan Aytekin 4
Accessing Single Elements
Elements of matrices are specified by subscripts (row,col). To reference a particular element in a matrix, specify its row and column number using the following syntax:
A(row, column)
Other ways
Linear indexing Logical indexing
19/10/2008 rsan Aytekin 5
Linear Indexing
MATLAB stores matrices and arrays as a single column of elements. This single column is composed of all of the columns from the matrix, each appended to the last. So, matrix A
A = [2 6 9; 4 2 8; 3 5 1] A = 2 6 9 4 2 8 3 5 1
is actually stored in memory as the sequence
2, 4, 3, 6, 2, 5, 9, 8, 1
The element at row 3, column 2 of matrix A (value = 5) can also be identified as element 6 in the actual storage sequence. To access this element, A(3,2), or A(6) which is referred to as linear indexing.
19/10/2008
rsan Aytekin
Logical Indexing
A logical array index designates the elements of an array A based on their position in the indexing array, B, not their value. Every true element in the indexing array is treated as a positional index into the array being accessed. Logical expressions returns a logical array entries of which are composed of true or false (1 or 0). Examples
A = rand(5); %5-by-5 matrix of random entries B = A > 0.5; %returns a logical matrix, 1 for the %entries satisfying the expression and 0 else B = 0 1 0 1 0 1 0 0 1 A(B) = 5; %the entries corresponding to the locations of 1s are assigned 5. The other ways to do this assignment are as follows: A(A > 0.5) = 5; or B = find(A > 0.5); A(B) = 5; %find is slower than logical indexing so previous ways are recommended
19/10/2008
rsan Aytekin
Some Useful Operations
MATLAB is operations. adapted to vectorial
MATLAB has useful functions and operators to handle vector and matrix operations.
19/10/2008
rsan Aytekin
find
Find indices and values of nonzero elements ind=find(X) locates all nonzero elements of array X, returns the linear indices of those elements in vector ind. logical expressions can be used Example 1 X = [1 0 4 -3 0 0 0 8 6]; indices = find(X) returns linear indices for the nonzero entries of X.
indices =
19/10/2008
9
9
rsan Aytekin
find
Example 2 You can use a logical expression to define X. For example, find(X > 2) returns linear indices corresponding to the entries of X that are greater than 2. ans =
19/10/2008
8
rsan Aytekin
9
10
find
Example 3 The following find command X = [3 2 0; -5 0 7; 0 0 1]; [r,c,v] = find(X) returns a vector of row indices of the nonzero entries of X r = 1 2 1 2 3
a vector of column indices of the nonzero entries of X c = 1 1 2 3 3
and a vector containing the nonzero entries of X. v =
19/10/2008
-5
1
11
rsan Aytekin
ind2sub, sub2ind
The ind2sub command determines the equivalent subscript values corresponding to a single index into an array. [I,J] = ind2sub(siz,IND) The sub2ind command determines the equivalent single index corresponding to a set of subscript values. IND = sub2ind(siz,I,J) IND = sub2ind(siz,I1,I2,...,In)
19/10/2008 rsan Aytekin 12
ind2sub
Example 1 matrix given below.
Two-Dimensional Matrices: The mapping from linear indexes to subscript equivalents for a 3-by-3
This code determines the row and column subscripts in a 3-by-3 matrix, of elements with linear indices 3, 4, 5, 6. IND = [3 4 5 6] s = [3,3]; [I,J] = ind2sub(s,IND)
I = 3 J = 1 2 2 2 1 2 3
19/10/2008
rsan Aytekin
13
ind2sub
Example 2 Three-Dimensional Matrices: The mapping from linear indexes to subscript equivalents for a 2-by-2-by-2 array given below. This code determines the subscript equivalents in a 2-by-2-by-2 array, of elements whose linear indices 3, 4, 5, 6 are specified in the IND matrix. IND = [3 4;5 6]; s = [2,2,2]; [I,J,K] = ind2sub(s,IND)
I = 1 1 J = 2 1 K = 1 2 1 2 2 1 2 2
19/10/2008
rsan Aytekin
14
sub2ind
Example
Create a 3-by-4-by-2 array, A. A = [17 24 1 8; 2 22 7 14; 4 6 13 20]; A(:,:,2) = A - 10 A(:,:,1) = 17 2 4 A(:,:,2) = 7 -8 -6 14 12 -4 -9 -3 3 -2 4 10 24 22 6 1 7 13 8 14 20
The value at row 2, column 1, page 2 of the array is -8. A(2,1,2) ans = -8 To convert A(2,1,2) into its equivalent single subscript, use sub2ind. sub2ind(size(A),2,1,2) ans = 14 You can now access the same location in A using the single subscripting method. A(14) ans = -8
19/10/2008
rsan Aytekin
15
reshape
Reshape array to specified form B = reshape(A,m,n) Returns the m-by-n matrix B whose elements are taken column-wise from A. An error results if A does not have m*n elements.
19/10/2008 rsan Aytekin 16
reshape
Examples Reshape a 3-by-4 matrix into a 2-by-6 matrix. A = 1 2 3 4 5 6 7 8 9 10 11 12
B = reshape(A,2,6) B = 1 3 5 2 4 6 B = reshape(A,2,[]) B = 1 3 5 2 4 6
7 8
9 10
11 12
7 8
9 10
11 12
19/10/2008
rsan Aytekin
17
repmat
Replicate and "tile" the elements of array A in an m-by-n arrangement:
A = [1 2 3; 4 5 6]; B = repmat(A,2,3); B = 1 2 3 1 4 5 6 4 1 2 3 1 4 5 6 4
2 5 2 5
3 6 3 6
1 4 1 4
2 5 2 5
3 6 3 6
19/10/2008
rsan Aytekin
18
circshift
Shift array circularly B = circshift(A,shiftsize) Circularly shifts the values in the array, A, by shiftsize elements. shiftsize is a vector of integer scalars where the nth element specifies the shift amount for the nth dimension of array A.
If positive, the values of A are shifted down (or to the right). If negative, the values of A are shifted up (or to the left). If 0, the values in that dimension are not shifted.
19/10/2008
rsan Aytekin
19
circshift
Examples
Circularly shift first dimension values down by 1. A = [ 1 2 3;4 5 6; 7 8 9] A = 1 2 3 4 5 6 7 8 9 B = circshift(A,1) B =
7 1 4
8 2 5
9 3 6
Circularly shift first dimension values down by 1 and second dimension values to the left by 1. B = circshift(A,[1 -1]); B = 8 9 7 2 3 1 5 6 4
Note that in MATLAB first dimension second dimension third dimension > column > row > depth
19/10/2008
rsan Aytekin
20
Colon (:) Operator
The rules to create regularly spaced vectors:
j:k is the same as [j,j+1,...,k] j:k is empty if j > k j:i:k is the same as [j,j+i,j+2i, ...,k] j:i:k is empty if i == 0, if i > 0 and j > k, or if i < 0 and j < k
where i, j, and k are all scalars.
19/10/2008 rsan Aytekin 21
Colon Operator
A(:,j) is the jth column of A A(i,:) is the ith row of A A(:,:) is the equivalent two-dimensional array. For matrices this is the same as A. A(j:k) is A(j), A(j+1),...,A(k) A(:,j:k)is A(:,j), A(:,j+1),...,A(:,k) A(:,:,k) is the kth page of three-dimensional array A. A(i,j,k,:) is a vector in four-dimensional array A. The vector includes A(i,j,k,1), A(i,j,k,2), A(i,j,k,3), and so on. A(:) is all the elements of A, regarded as a single column. On the left side of an assignment statement, A(:) fills A, preserving its shape from before. In this case, the right side must contain the same number of elements as A.
19/10/2008 rsan Aytekin 22
Colon Operator
Examples Using the colon with integers, D = 1:4 results in D = 1 2 3 4 Using two colons to create a vector with arbitrary real increments between the elements, E = 0:.1:.5 results in E = 0 0.1000 0.2000 0.3000 0.4000 0.5000
19/10/2008
rsan Aytekin
23
Colon Operator
Examples
Selecting a column of a.
a(2,:) a(:,3) a(2:4,:) a(:,3:5) a(2:4,3:5) a([Link],:)
You can put a vector into a row (col.) position within a.
a(:,[1 2 5]) a([2 5],[2 4 5])
19/10/2008 rsan Aytekin 24
Concatenation
The process of joining small matrices to make bigger ones. [] operator B = [A A+32; A+48 A+16] Also cat can be used. Concatenate arrays along specified dimension
A = 1 3 2 4 B = 5 7 6 8
concatenating along different dimensions produces
19/10/2008
rsan Aytekin
25
Deleting Rows and Columns
You can delete rows and columns from a matrix using just a pair of square brackets. Start with X = A; Then, to delete the second column of X, use X(:,2) = [];
19/10/2008 rsan Aytekin 26
Matrix and Arithmetic Operations
Arithmetic operations on arrays are done element by element. MATLAB uses a dot, or decimal point, as part of the notation for multiplicative array operations. The list of operators includes
+ .* ./ .\ .^ .' Addition Subtraction Element-by-element multiplication Element-by-element division Element-by-element left division Element-by-element power Unconjugated array transpose
19/10/2008
rsan Aytekin
27
Image Processing Toolbox
MATLAB stores most images as two-dimensional arrays (i.e., matrices). Each element of the matrix corresponds to a single pixel in the displayed image. The Pixel Coordinate System
19/10/2008
rsan Aytekin
28
Properties of Image Types
Image Type Binary (Also known as a bilevel image) Grayscale (Also known as an intensity, gray scale, or gray level image) Interpretation
Logical array containing only 0s and 1s, interpreted as
black and white, respectively.
Array of class uint8, uint16, int16, single, or double whose pixel values specify intensity values. For single or double arrays, values range from [0, 1]. For uint8, values range from [0,255]. For uint16, values range from [0, 65535]. For int16, values range from [32768, 32767]. m-by-n-by-3 array of class uint8, uint16, single, or double whose pixel values specify intensity values. For single or double arrays, values range from [0, 1]. For uint8, values range from [0, 255]. For uint16, values range from [0, 65535].
Truecolor (Also known as an RGB image )
19/10/2008
rsan Aytekin
29
Binary Images
19/10/2008
rsan Aytekin
30
Gray Scale Images
19/10/2008
rsan Aytekin
31
Truecolor Images
Each pixel is specified by three values;one each for the red, blue, and green components of the pixel's color. Truecolor images are stored as an m-by-n-by-3 data array that defines red, green, and blue color components for each individual pixel. Truecolor to grayscale conversion rgb2gray(I)
19/10/2008 rsan Aytekin 32
Image Arithmetic Saturation Rules
The results of integer arithmetic can easily overflow the data type allotted for storage. For example, the maximum value you can store in uint8 data is 255. Arithmetic operations can also result in fractional values, which cannot be represented using integer arrays. Rules for integer arithmetic:
Values that exceed the range of the integer type are saturated to that range. Fractional values are rounded.
For example, if the data type is uint8, results greater than 255 (including Inf) are set to 255. The following table lists some additional examples.
Result 300 -45 10.5 Class uint8 uint8 uint8 Truncated Value 255 0 11
19/10/2008
rsan Aytekin
33
Reading and Writing Images
Read and Display an Image imread: reads [Link], and stores it in an array named I. Then show the image I = imread('[Link]'); figure, imshow(I) Write the Image to a Working or Specified Disk File You can also specify other formats such as .jpg, .bmp, tiff etc. imwrite (I2, '[Link]'); imwrite (I2, 'D:\EE583\[Link]');
19/10/2008 rsan Aytekin 34
Some Useful Commands for Image Processing
19/10/2008
rsan Aytekin
35
im2col,col2im
im2col Rearrange image blocks into columns B = im2col(A,[m n],block_type) blocktype: distinct or sliding col2im Rearrange matrix columns into blocks A = col2im(B,[m n],[mm nn],'distinct') A = col2im(B,[m n],[mm nn],'sliding') If B = [A11(:) A21(:) A12(:) A22(:)], where each column has length m*n, then A = [A11 A12; A21 A22] where each Aij is m-by-n.
19/10/2008 rsan Aytekin 36
im2col,col2im
Example
B = reshape(uint8(1:25),[5 5])' C = im2col(B,[1 5]) A = col2im(C,[1 5],[5 5],'distinct') B = 1 6 11 16 21 C = 1 2 3 4 5 A = 1 6 11 16 21 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 2 7 12 17 22 3 8 13 18 23 4 9 14 19 24 5 10 15 20 25
19/10/2008
rsan Aytekin
37
imresize
Resize image
B = imresize(A, scale) B = imresize(A, [mrows ncols]) B = imresize(A, scale) returns image B that is scale times the size of A. Interpolation methods: bilinear, bicubic.
The input image A can be a grayscale, RGB, or binary image.
19/10/2008 rsan Aytekin 38
imresize
Example
I = imread('[Link]'); J = imresize(I,1.25); imshow(I) figure, imshow(J)
19/10/2008
rsan Aytekin
39
imcrop
Crop rectangular image
I2 = imcrop(I) I2 = imcrop(I, rect)
rect is a four-element position vector [xmin ymin width height] that specifies the size and position (upper-right corner) of the crop rectangle. Note that xmin corresponds to column, ymin corresponds to row.
I can be a grayscale image, truecolor image, or logical array.
19/10/2008
rsan Aytekin
40
imcrop
Example I = imread('[Link]'); I2 = imcrop(I,[75 68 130 112]); imshow(I), figure, imshow(I2)
19/10/2008
rsan Aytekin
41
Filtering in the Spatial Domain
Filtering: neighborhood operation, in which the value of any given pixel in the output image is determined by applying some algorithm to the values of the pixels in the neighborhood of the corresponding input pixel. Linear filtering: the value of an output pixel is a linear combination of the values of the pixels in the input pixel's neighborhood. MATLAB has two types of convention for linear filtering:
convolution correlation
19/10/2008
rsan Aytekin
42
Convolution
Linear filtering of an image, also called neighborhood operation Each output pixel is the weighted sum of neighboring input pixels. The matrix of weights is called the convolution kernel, also known as the filter. A convolution kernel is a correlation kernel that has been rotated 180 degrees (flip horizontally and then vertically or the other way around).
19/10/2008 rsan Aytekin 43
Convolution
Example Suppose the convolution kernel is h = [8 3 4 1 5 9 6 7 2]
The following figure shows how to compute the (2,4) output pixel 1.2+8.9+15.4+7.7+14.5+16.3+13.6+20.1+22.8 = 575
19/10/2008
rsan Aytekin
44
Correlation
Closely related to convolution. Each output pixel is also computed as a weighted sum of neighboring pixels. The matrix of weights, called correlation kernel, is not rotated during the computation. The Image Processing Toolbox filter design functions return correlation kernels.
19/10/2008 rsan Aytekin 45
Correlation
Example Suppose the correalation kernel is h = [8 3 4 1 5 9 6 7 2]
The (2,4) output pixel from the correlation is 1.8+8.1+15.6+7.3+14.5+16.7+13.4+20.9+22.2 = 585
19/10/2008
rsan Aytekin
46
Performing Linear Filtering of Images Using imfilter
Filtering of images, either by correlation or convolution, by the toolbox function imfilter. Default operation: correlation filtering Example Filter an image with a 5-by-5 filter containing equal weights (averaging filter). I = imread('[Link]'); h = ones(5,5) / 25; I2 = imfilter(I,h); imshow(I), title('Original Image'); figure, imshow(I2), title('Filtered Image')
19/10/2008 rsan Aytekin 47
Performing Linear Filtering of Images Using imfilter
19/10/2008
rsan Aytekin
48
Morphological Operations
Morphology is image processing operations that process images based on shapes. Morphological operations apply a structuring element to an input image, creating an output image of the same size. The most basic morphological operations are dilation and erosion.
19/10/2008 rsan Aytekin 49
Basic Morphological Operations Dilation & Erosion
Dilation
The value of the output pixel is the maximum value of all the pixels in the input pixel's neighborhood. In a binary image, if any of the pixels is set to the value 1, the output pixel is set to 1.
Erosion
The value of the output pixel is the minimum value of all the pixels in the input pixel's neighborhood. In a binary image, if any of the pixels is set to 0, the output pixel is set to 0. dilation > output=1 erosion > output=0
19/10/2008
rsan Aytekin
50
Performance Issues
The idea: MATLAB is
very fast on vector and matrix operations correspondingly slow with loops
Try to avoid loops Try to vectorize your code Try to use block processing
rsan Aytekin 51
19/10/2008
Vectorize Loops
Example Given an=n, and bn =1000n for n=1,...,1000. Calculate Sum{[Link] } for n=1 to 1000. %poor style tic a = 1:1000; b = 1000 - a; ssum=0; for n=1:1000 ssum = ssum +a(n)*b(n); end toc %better style: Note that sum is inner product. tic ssum = a*b'; toc %computation time is almost 2,5 times faster Repmat, reshape, ind2sub, sub2ind, find use this approach. Follow the guides to vectorize on the internet. One link is given below: [Link]
19/10/2008
rsan Aytekin
52
Neighborhood or Block Processing
Certain image processing operations involve processing an image in sections, called blocks or neighborhoods, rather than processing the entire image at once. Instead of using loops employ toolbox block functions Steps of Block processing
19/10/2008
break the input image into blocks or neighborhoods, call the specified function to process each block or neighborhood, reassemble the results into an output image.
rsan Aytekin 53
Neighborhood or Block Processing
Two approaches sliding neighborhood operation, nlfilter
performed a pixel at a time, so overlapping occurs the value of any given pixel in the output image is determined by the application of an algorithm to the values of the corresponding input pixel's neighborhood.
distinct block processing, blkproc
divide a matrix into m-by-n sections. These sections, or distinct blocks, overlay the image matrix starting in the upper left corner, with no overlap.
rsan Aytekin 54
19/10/2008
Neighborhood or Block Processing
Example
I = imread('[Link]'); f = @(x) uint8(round(mean2(x)*ones(size(x)))); I2 = blkproc(I,[8 8],f); imshow(I) figure, imshow(I2);
19/10/2008
rsan Aytekin
55