ECE565 Project 2 repot PDF

Preview

Summary

Download ECE565 Project 2 repot PDF

Description

ECE565 Project #2 Name: Phanimithra Gangalakurti ID: A20340643 I declare that all the work performed for this project and presented in the document below are what I have worked on and developed myself and have not taken from any other/older work.

Question 1: Write a global thresholding program in which the threshold is estimated automatically using the procedure discussed in Section 10.3.2. The output of your program should be a segmented (binary) image. Use your program to segment “noisy_fingerprint.tiff” and produce a segmented image. The global thresholding algorithm generated threshold is 125. Here is the thresholded image:

Matlab Source Code: close all; clear all; f = imread('noisy_fingerprint.tif');

g = global_thresholding(f,1);

Functions: function [g Threshold] = global_thresholding(f,enable_plots); L = 256; if enable_plots figure; imshow(f); end % Compute histogram of input image h = imhist(f); if enable_plots figure; stem(h); end % Find min and max intensity of the image f_min = double(min(min(f(find(f~=0))))) f_max = double(max(max(f(find(f~=0))))) % Initial threshold T = (f_min + f_max)/2 threshold_limit = 2; delta_T = L; % Iteratively determine the threshold while ( delta_T > threshold_limit ) mean_class1 = ([0:(T(end)-1)]*h(1:T(end)))/sum(h(1:T(end))) mean_class2 = ([T(end):L-1]*h(T(end)+1:end))/sum(h(T(end)+1:end)); new_threshold = 0.5*(mean_class1+mean_class2) T = [T new_threshold]; delta_T = abs(T(end) - T(end-1)) end Threshold = round(T(end)); g = f; g(find(fThreshold)) = 255; if enable_plots figure; imshow(g); figure; imhist(g); end end

Question 2: a. Implement Otsu’s optimum thresholding algorithm from scratch by yourself. Use the algorithm to segment provided image: Here is the source code (function) for Otsu’s optimum thresholding: Matlab source code: close all; clear all; f = imread('polymersomes.tif'); [g T_otsus] = otsus_thresholding(f,1); [g_global T_global] = global_thresholding(f,1);

Functions: function [g Otsus_k] = otsus_thresholding(f,enable_plots) % Otsu's method of maximum between class variance thresholding % clear all; % close all; if enable_plots figure; imshow(f); end f_double = double(f); L = 255; for idx = 0:L p(idx+1) = length(find(f_double==idx)); end % Probability density function of the image intensities p = p/prod(size(f)); % Threshold iteration for k = 0:L idx = k+1; P1(idx) = sum(p(1:idx)); P2(idx) = 1-P1(idx); m1(idx) = (p(1:idx) * [0:k]')/P1(idx); m2(idx) = (p(idx+1:end) * [k+1:L]')/P2(idx); mG(idx) = P1(idx) * m1(idx) + P2(idx) * m2(idx); % Between class variance Sigma_B_sq(idx) = P1(idx)*(m1(idx)-mG(idx))^2 + P2(idx)*(m2(idx)mG(idx))^2; end

[value Otsus_k] = max(Sigma_B_sq); Otsus_k g = f; g(find(fOtsus_k)) = 255; if enable_plots figure;imshow(g) end end

Here is the image thresholded using Otsu’s binary threshold:

b. Here is the global thresholded image:

Question 3: a. Generate smoothed image g using 9x9 averaging filter: Here is the smoothed image with a 9x9 averaging window using imfilter built in function

b. Generate a binary image gB by thresholding image g obtained in part (a): The image is thresholded with 60% of max intensity of the image. Here is the thresholded image:

c. Extract the outer boundary of gB and display the results as a binary image Boundary extraction routine using boundary tracing algorithm has been implemented (code included after results)

d. Subsample the boundary obtained in Part.C into a grid whose lines are separated by 50 pixels. Subsampling has been done using the bsubsamp function provided in the appendix of the project question paper.

Image of subsampled image points connected by line segments:

e. Write a program that computes the Freeman chain code ‘c’ of the boundary b with the code connectivity specified by CONN (i.e. c = fchcode(b,CONN)).

Results from the function fchcode: c = struct with fields: fcc: [0 0 0 0 6 0 diff: [0 0 0 6 2 6 mm: [0 0 0 0 6 0 diffmm: [0 0 0 0 0 0 x0y0: [154 460]

6 0 6 0

6 0 6 6

6 0 6 0

6 0 6 0

6 0 6 0

6 0 6 0

6 0 6 0

6 6 6 6

4 0 4 2

4 0 4 6

4 0 4 0

4 0 4 0

4 0 4 0

4 6 4 0

2 2 2 6

4 6 4 2

2 0 2 0

2 0 2 6

2 0 2 2

2 0 2 6

2 6 2 0

0 2 0 0

2 0 2 0

2 6 2 6

0 2 0 2

2] 6] 2] 6]

Matlab source code: clear all; close all; f = imread('circular_stroke.tif'); figure; imshow(f); title('Input image'); % Smoothing the image w = ones(9); g = imfilter(double(f),w,'replicate'); figure; imshow(g,[]); title('Smoothed image to remove specular noise'); % Threhsolding the image T = max(g(:)); gB = zeros(size(g)); gB = g >= 0.6*T; figure; imshow(gB); title('Thresholded binary image, specular noise removed'); % Boundary tracing (trace boundary and generates correpsonding coordinates) % Did not use bwboundaries, extracted boundary by the boundary tracing algorithm [gBound gBound_coord] = outer_boundary_trace(gB); figure; imshow(gBound); title('Outer boundary traced image'); % Verify if generated coordinates are correct with reference function boundIm = bound2im(gBound_coord,size(gB,1),size(gB,2)); figure; imshow(boundIm); title('Validated boundary co-ordinates computed by the function'); err = sum(sum(gBound-boundIm)); if err error('Problem with generated coordinates'); else disp('Phew...finally coordinates and boundary image match...!!! :)'); end gS_coord = bsubsamp(gBound_coord,50); gS = bound2im(gS_coord,size(gB,1),size(gB,2)); figure; imshow(gS); title('Boundaries sub-sampled with a 50 pixel grid separation'); gSc_coord = connectpoly(gS_coord(:,1),gS_coord(:,2)); gSc = bound2im(gSc_coord,size(gB,1),size(gB,2)); figure; imshow(gSc); title('Sub-sampled image co-ordinates connected by line segments'); CONN = 8; c = fchcode(gS_coord,CONN)

Functions: Function1: function [gBound,gBound_coord] = outer_boundary_trace(f); % Expects a binary input image % Finds the boundary of the binary image. Handles only images that do not % connect to image edges. Computes binary image gBound of the same % dimension as f and traces the outer boundary of image in f % gBound_coord contains NumBoundaryPixels x 2 dimensioned coordinates of % the boundary pixels starting from topmost leftmost boundary pixel % Find uppermost-leftmost point [y x] = find(f'==1); b0 = [x(1) y(1)]; b = b0; c = b-[1 0]; x_coord = [b(1)]; y_coord = [b(2)]; gBound = zeros(size(f)); gBound(b(1),b(2)) = 1; % Map of pixels in the 3x3 bounding box in clockwise direction starting % from top left pixel_extract_idx = [1 4 7 8 9 6 3 2]; offset = 1; boundary_closed = 0; while ~boundary_closed next_fg_pixel_found = 0; wind_raw = f(b(1)+[-1:1],b(2)+[-1:1]); wind_vect_raw = reshape(wind_raw,1,[]); [y x] = meshgrid(-1:1,-1:1); % x represents rows and y columns in images x_vect = reshape(x,1,[]); y_vect = reshape(y,1,[]); x_vect_cs y_vect_cs

= x_vect(circshift(pixel_extract_idx,-offset)); = y_vect(circshift(pixel_extract_idx,-offset));

x_vect = x_vect(pixel_extract_idx); y_vect = y_vect(pixel_extract_idx); pixel_idx

= 0;

while ~next_fg_pixel_found pixel_idx = pixel_idx + 1; if f(b(1)+x_vect_cs(pixel_idx),b(2)+y_vect_cs(pixel_idx)) next_fg_pixel_found = 1; coord_step = [x_vect_cs(pixel_idx) y_vect_cs(pixel_idx)];

if pixel_idx == 1 pixel_idx_m1 = 8; else pixel_idx_m1 = pixel_idx - 1; end coord_step_m1 = [x_vect_cs(pixel_idx_m1) y_vect_cs(pixel_idx_m1)]; c = b + coord_step_m1; new_b = b + coord_step; b = new_b; x_coord = [x_coord; b(1)]; y_coord = [y_coord; b(2)]; gBound(b(1),b(2)) = 1; a = (c-b); for offset = 1:8 if (a(1)==x_vect(offset) && a(2)==y_vect(offset)) break else continue end end end end if (b0(1)==b(1) && b0(2)==b(2)) boundary_closed = 1; end end gBound_coord = [x_coord y_coord]; end

Function2: function c = fchcode(b,conn); % inp_co_ord = gS_coord; % conn = 8; inp_co_ord = b; b = [inp_co_ord(end,:); inp_co_ord]; % Find out phase changes interms of signs (magnitude does not matter) coordinate_diff = diff(b); % The 4/8 connectivity vectors have zero pointing in horizontal right % direction but that is the direction of y. Hence swapped real/imag

dir_vector = coordinate_diff(:,2)+j*coordinate_diff(:,1); % Prefixed a 0 to account for the freeman_chain_code = [mod(-angle(dir_vector)/(2*pi/conn),conn)]; c.fcc = freeman_chain_code'; % Computing first diff first_diff = mod(circshift(freeman_chain_code,-1) - freeman_chain_code,conn); c.diff = first_diff';

% Finding the integer of minimum magnitude for Freeman Chain Code sequence seq_to_check = freeman_chain_code; non_zero_coords = [1; find(seq_to_check ~= 0)]; dist_between_non_zeros = diff(non_zero_coords); [max_val max_zeros_prec_idx] = max(dist_between_non_zeros); c.mm = circshift(seq_to_check,-(non_zero_coords(max_zeros_prec_idx)-1))'; % Finding integer of minimum mangnitude for First diff sequence seq_to_check = first_diff; non_zero_coords = [1; find(seq_to_check ~= 0)]; dist_between_non_zeros = diff(non_zero_coords); [max_val max_zeros_prec_idx] = max(dist_between_non_zeros); max_zeros_start_loc = non_zero_coords(max_zeros_prec_idx)+1; c.diffmm = circshift(seq_to_check,-(max_zeros_start_loc-1))'; c.x0y0 = inp_co_ord(max_zeros_start_loc,:); return;

Question 4: a. Extract the boundary of the chromosome and display the result as a binary image: I’ve used the function outer_boundary_trace copied above for extracting the boundary of the chromosome. Here is the binary image:

b. Write a program to compute the Fourier descriptors of a boundary S (i.e. z = fourierdescp(s)). The input S is an np x 2 sequence of ordered coordinates describing the boundary and the output z is a sequence of Fourier descriptors obtained. Compute the Fourier descriptors for the boundary obtained in Part.a Matlab source code: close all; clear all; f = imread('chromosome.tif'); figure; imshow(f); title('Input image'); [gB,gB_coord] = outer_boundary_trace(f); figure; imshow(gB); title('Boundary of the Chromosome extracted from input image'); gBreconstruct = bound2im(gB_coord,size(gB,1),size(gB,2)); figure; imshow(gBreconstruct); title('Validation of boundary description coordinates'); % Compute Fourier descriptors z = fourierdescp(gB_coord); % Reconstruct image from Fourier descriptors nd = size(gB_coord,1); gB_coord_nd = ifourierdescp(z,nd); gBreconstruct_nd = bound2im(gB_coord_nd,size(gB,1),size(gB,2)); figure; imshow(gBreconstruct_nd); title('Binary image reconstructed from Fourier descriptors'); % Reconstruct image considering significant 50% Fourier descriptors nd = size(gB_coord,1)*0.5; gB_coord_50pct = ifourierdescp(z,nd); gBreconstruct_50pct = bound2im(gB_coord_50pct,size(gB,1),size(gB,2)); figure; imshow(gBreconstruct_50pct); title('Binary image reconstructed considering 50% of Fourier descriptors'); % Reconstruct image considering significant 1% Fourier descriptors nd = floor(size(gB_coord,1)*0.01); nd = nd + mod(nd,2); gB_coord_1pct = ifourierdescp(z,nd); gBreconstruct_1pct = bound2im(gB_coord_1pct,size(gB,1),size(gB,2)); figure; imshow(gBreconstruct_1pct); title('Binary image reconstructed considering 1% of Fourier descriptors');

Functions: Function1: function z = fourierdescp(f); s = f(:,1) + j * f(:,2); % Performing fftshift on DFT output so that the spectrum is from [-pi:pi] z = fftshift(fft(s,size(f,1))); return;

c. Write a program to compute the inverse Fourier descriptors (i.e. s = ifourierdescp(z,nd)). The input z is a sequence of Fourier descriptors and nd is the number of descriptors used to compute the inverse. nd must be an even integer no greater than length(z). Reconstruct the boundary using 50% of the total possible descriptors and display the result as a binary image. Then, reconstruct the boundary using 1% of the total possible descriptors and display the result as a binary image.

Matlab source code (function): function g = ifourierdescp(z,nd); if nd>size(z,1) error('nd should be less than or equal to length of Fourier descriptor array'); elseif mod(nd,2) error('nd should be an even number'); end z_processed = zeros(size(z)); z_processed(size(z,1)/2+[-nd/2+1:nd/2]) = z(size(z,1)/2+[-nd/2+1:nd/2]); s = ifft(ifftshift(z_processed),size(z,1)); g = [real(s) imag(s)]; return;

Image of the chromosome reconstructed from 50% Fourier descriptors:

Image reconstructed using 1% Fourier descriptors:...