This weeks assignment is a quick introduction to MATLAB for image analysis. The original file is found under the "Undervisningsplan" on the course web site. The solution to this exercise will be uploaded when it is thirthy minutes left of the group lecture.
If you are on a Red Hat machine, MATLAB can be launched by the command 'matlab' (without the quotes) at the command prompt. In Windows there should be a MATLAB icon in the start-menu or on the desktop.
(Which, when I publish it online, becomes just simple text)
Make a folder named 'INF4300' and set this as the 'current directory' using the UI or using the MATLAB-function 'cd'. When using the CD function and to use MATLAB more efficiently, the autocompletion is often handy. Activate it by hitting 'Tab'. Save a copy of this file in that folder, and open it in MATLAB. To run a part of the code, select the code and hit F9. To run the selected block (the block with the caret), hit 'Ctrl+Enter'. To run an entire file, you can: - type the name of the file, e.g. 'intro', at the command prompt, or - select the file and hit F9, or - open the file and hit F5.
% The matrix | 1 0 0 | % | 0 1 0 | % | 0 0 1 | % may be created like this: someMatrix = [1 0 0 ; 0 1 0 ; 0 0 1] % e.g. [row1 ; row2 ; ... ; rowM], where a row can be defined as % [number1 number2 ... numberN] or [number1, number2, ..., numberN] % MATLAB comes built-in with function for many common operations, % for the case above we could have typed: someMatrix = diag([1 1 1]) % If you wonder how a function/operator works, type 'help <function name>' % or for (a sometimes more thorough) description in a separate window, % type 'doc <function name>'. % It is also possible to simply press F1 while the cursor is at a certain % function to access MATLAB's help for that particular funtion. % MATLAB matrices are dynamic objects, meaning one need not know the exact % size of a matrix in advance. Initially a simple matrix can be made, and % and later rows or colums may be added or removed. If, however, the size % is known in advance and the matrix is large, it is wise to preallocate % memory for the matrix before you start using it. This may greatly improve % performance. % E.g. allocate memory for 100*100 matrix and fill it with 1's: big_1 = ones(100); %<-- Terminating a line with ';' will suppress output. % Or, allocate a 50*40*3 matrix and fill it with 0's: big_0 = zeros([50 40 3]); % A vector can be created using the '<from>:<step>:<to>' syntax: steg1 = 1:2:50; % The <step> is optional. If omitted the step defaults to 1: steg2 = 1:50; % The standard plotting function is PLOT, e.g.: plot(steg1); % An alternative is STEM. Look up how it works using 'help stem'. % You can add a title to the plot by using the title() function: title('Demonstration of plot()'); % Extracting a value, or a series of values, from a matrix is easily % achieved like this: mat = [1 7 4 ; 3 4 5] mat(2,1) % retrives the '3'. NOTE: The first index is 1 in MATLAB! % MATLAB lays out a column at a time in memory, hence the value '7' can % either be retrieved using linear index: mat(3) % or using (row,column)-index: mat(1,2) % A range of elements may retrieved using the <from>:<to> syntax: mat(1:2,1) % First column mat(1:end,1) % The same, since n = 2 mat(:,1) % The same, : is here 'all rows' % We can use the same syntax to set a range of elements: mat(1:2,1:2) = 0 % Note that a matrix can be stored in various formats, e.g. UINT8, INT8 or % DOUBLE. They all have their conversion functions, see e.g. 'help uint8'. % If your result looks fishy, and you have a hard time figuring out why, % a type conversion (of lack of one) may be the reason!
someMatrix = 1 0 0 0 1 0 0 0 1 someMatrix = 1 0 0 0 1 0 0 0 1 mat = 1 7 4 3 4 5 ans = 3 ans = 7 ans = 7 ans = 1 3 ans = 1 3 ans = 1 3 mat = 0 0 4 0 0 5
Here we'll demonstrate how some simple matrix operations is done in MATLAB
% We'ss start mat1 = [2 3; 5 7]; mat2 = [4 5; 1 9]; % Elementwise addition and subtraction. mat1 + mat2 mat1 - mat2 % Transpose. mat1.' % Complex conjugate transpose (same as transpose when the matrix is real). mat1' % But notice the difference for a complex matrix complexMatrix = [0 1i-2; 1i+2 0] % Transpose for complex, see the function transpose complexMatrix.' % Complex conjugate transpose, see the function ctranspose complexMatrix' % Multiplication, division and power. mat1 * mat2 mat1 / mat2 mat1^2 % Elementwise multiplication, division and power % (add '.' in front of the operator). mat1 .* mat2 mat1 ./ mat2 mat1 .^ mat2
ans = 6 8 6 16 ans = -2 -2 4 -2 ans = 2 5 3 7 ans = 2 5 3 7 complexMatrix = 0.0000 + 0.0000i -2.0000 + 1.0000i 2.0000 + 1.0000i 0.0000 + 0.0000i ans = 0.0000 + 0.0000i 2.0000 + 1.0000i -2.0000 + 1.0000i 0.0000 + 0.0000i ans = 0.0000 + 0.0000i 2.0000 - 1.0000i -2.0000 - 1.0000i 0.0000 + 0.0000i ans = 11 37 27 88 ans = 0.4839 0.0645 1.2258 0.0968 ans = 19 27 45 64 ans = 8 15 5 63 ans = 0.5000 0.6000 5.0000 0.7778 ans = 16 243 5 40353607
% While and for loops. % Can often be used instead of matrix operations, but tend to lead to a % slightly inferior perfomance, especially in older versions of MATLAB. % % For some tips on how to write fast MATLAB-code: % http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=5685 % % The reason for introducing performance-tips already is because we will be % storing images as matrices, and these tend to be large. If we also % consider that the algorithms we use are computationally intensive, % writing somewhat optimised code is often worth the extra effort. % E.g. can this simple loop: A = zeros(1,100); for i = 1:100 A(i) = 2; end % be replaced by a single and more efficient command: A = ones(1,100) * 2; % And why write: for i = 1:100 for j = 1:100 stor_1(i,j) = i^2; % MATLAB warning: consider preallocating for speed end end % when you can do this: big_1_temp = repmat((1:100)'.^2, 1,100); % or alternatively: big_1_temp2 = (1:100)'.^2 * ones(1,100); % Note: % If you are not used to 'doing everything with vectors', you'll likely % want to use loops far more often than you need to. When you get the hang % of it, typing vector-code is faster and often less error-prone than the % loop-equivalents, but loops offer better flexibility. % As always, add lots of COMMENTS, especially for complex one-liners! % Logical structures: if-elseif-else. Which number is displayed? if -1 1 elseif -2 2 else 3 end % Logical structures: switch-case. var = 't'; switch var case 't' 1 case 2 2 end
ans = 1 ans = 1
% A very handy method, but may initially seem a bit tricky. % Example: Find elements which are equal to 4. M = [1 2 3 ; 4 4 5 ; 0 0 2] M == 4 % a logical map I = find(M == 4) % the linear indices [II, J] = find(M == 4) % the (row,column)-indices M(M == 4) = 3 % change them to 3 M(I) = 4 % and back to 4
M = 1 2 3 4 4 5 0 0 2 ans = 0 0 0 1 1 0 0 0 0 I = 2 5 II = 2 2 J = 1 2 M = 1 2 3 3 3 5 0 0 2 M = 1 2 3 4 4 5 0 0 2
Now, we'll finally look at some images. This is image alaysis, isn't it?
% Open a couple of images. img1 = imread('football.jpg'); img2 = imread('coins.png'); % Display the first image. imshow(img1) title('The first image'); % A title is nice % Display the second image in a new figure. figure, imshow(img2) title('The second image'); % Another function for displaying the images is: figure, imagesc(img1), colorbar title('Image displayed by imagesc'); % Histogram the image intensity values. H = imhist(img2); plot(H) title('The histogram for the second image'); % Converting a colour image with 3 colour channels to a greyscale image % (the Y in the YIQ colour model). img1 = imread('football.jpg'); img3 = rgb2gray(img1); figure, imshow(img3) title('Grayscale image of football');
Now we have jumped into the curriculum of INF2310. You don't remember, or don't know?! Don't hesitate to ask or read up on last courses materials awailable at http://www.uio.no/studier/emner/matnat/ifi/INF2310/v15/
% 2D correlation: FILTER2 % 2D convolution: CONV2 % With IMFILTER you can choose either 2D correlation (default) and 2D % convolution, and it also has some boundary options. % Filter IMG3 using a 5x5 mean filter. h1 = ones(5,5) / 25; img4 = imfilter(img3,h1); % symmetric filter => correlation = convolution imshow(img3), title('Original image'); figure, imshow(img4), title('Filtered image') % Find the gradient magnitude of IMG3. h2x = [-1 -2 -1 ; 0 0 0 ; 1 2 1] h2y = [-1 0 1 ; -2 0 2 ; -1 0 1] resX = conv2(double(img3), h2x); % NOTE: DOUBLE type conversion resY = conv2(double(img3), h2y); resXY = sqrt(resX.^2 + resY.^2); % Display the gradient magnitude, but not like this: figure, imshow(resXY); title('Not like this'); % because the assumed range of the DOUBLE type is [0,1], % but e.g. one of these ways: figure, imshow(resXY, ); title('But like this,'); figure, imshow(resXY, [min(resXY(:)) max(resXY(:))]); title('this') figure, imshow(resXY/max(resXY(:))); title('this or') figure, imshow(uint8(resXY/max(resXY(:)).*255)); title('this. They are all equal.'); %%%%%%%%%%%%%%%%
h2x = -1 -2 -1 0 0 0 1 2 1 h2y = -1 0 1 -2 0 2 -1 0 1
%%%%%%%%%%%%%%%% % Above, we loaded the image 'football.jpg', converted it to a greyscale % image and applied a 5x5 mean filter, by using the commands: img1 = imread('football.jpg'); img3 = rgb2gray(img1); figure, imshow(img3) h1 = ones(5,5) / 25; img4 = imfilter(img3,h1); imshow(img3), title('Original image'); figure, imshow(img4), title('Filtered image') % The filtered image have a two pixel wide black frame. % a) Use the indexing techniques described above to remove these. % b) Use FILTER2 and CONV2 with the option 'valid' to remove these. % Hint 1: Type convert IMG3 to DOUBLE. % Hint 2: The image should be the second argument to FILTER2, % not the first as it is to CONV2 and IMFILTER. % c) Use the boundary option of IMFILTER to get a same-sized image without % the black frame. %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% % Make a function that returns the same as IMHIST when the parameter is a % 8-bits greyscal image. % % Although it is allowed to use loops, try to avoid using them where % it is possible. One loop should suffice. % % Hint: How to create and use the function % A function is stored in .m files with the same name as the function. % 1. Create an m-file named 'histogram.m', e.g. using 'edit histogram.m' % 2. The first line should be 'function h = histogram(img)' % 3. Write the code to produce a histogram of IMG below the function % declaration. The histogram should be stored in a variable named H. % 4. Save the file. % 5. From another m-file, or from the command line call HISTOGRAM(IMG) % where IMG is a greyscale image. %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% % Above, we loaded the image 'coins.png' using the command: img2 = imread('coins.png'); % % a) % Use the operators >, <, >=, <= to threshold IMG2 using an arbitrary % threshold. % % b) % Image Processing Toolbox (IPT) in MATLAB have a function for computing % the 'optimal' threshold based on Otsu's algorithm. % - Find this function using MATLAB's help system. % - Use it to with the 'optimal' threshold. % - Use the threshold and the function IM2BW to threshold IMG2. % % c) % Compare the binary image resulting from part a with the one from part b % by displaying the images. Do you notice any differences? % Display also the difference between the images. %%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% % Normalize resXY such that max(resXY(:)) = 255 and min(resXY(:)) = 0. % Threshold the result with T = 100. % % What would you do if you wanted to obtain an image containing % only the seam, and the entire seam, of the the ball?