Matlab图像分割(隐Markov随机场模型+期望最大算法)

% Note % ---- % The code is an implementation which refers to the following paper and his % codes. % % Quan Wang. HMRF-EM-image: Implementation of the Hidden Markov Random Field % Model and its Expectation-Maximization Algorithm. % arXiv:1207.3510 [cs.CV], 2012. % % His code on Github: https://github.com/wq2012/HMRF-EM-image classdef HMRF_EM properties (Access=public) % Max iter step of MAP procedure. Default is 10. MAP_iter = 10 % Max iter step of EM(expectation-maximization) procedure. Default is 10. EM_iter = 10 % Number of clusters generated by k-means algorithm. Default is 2. kmeans_num = 2 % Arguments cell for `edge` function to extract edge of image. The cell % fills in the second and subsequent parameters of the `edge` function % in sequence. edgeArgs = {'canny', 0.75} % Arguments of image filter. The cell fills in the second and subsequent % parameters of the `imfilter` function in sequence. imfilterArgs end methods (Access=public) function newcls = HMRF_EM(kmeans_num, MAP_iter, EM_iter) % Initialization of `HMRF_EM` class. % >>> Parameters % kmeans_num: Number of clusters generated by k-means algorithm. % MAP_iter: Max iter step of MAP procedure. % EM_iter: Max iter step of EM(expectation-maximization) procedure. % >>> Return % newcls: The `HMRF_EM` object. % Check environment runInOctave = isOctaveEnv(); if runInOctave==1 % Load `image` and `statistics` package when the code is running in % GNU Octave. pkg load statistics; pkg load image; end % Use `validateattributes` to check input switch nargin case 0 MAP_iter = 10; EM_iter = 10; kmeans_num = 2; case 1 MAP_iter = 10; EM_iter = 10; validateattributes(kmeans_num,{'numeric'},{'scalar','positive'}); case 2 EM_iter = 10; validateattributes(kmeans_num,{'numeric'},{'scalar','positive'}); validateattributes(MAP_iter,{'numeric'},{'scalar','positive'}); otherwise validateattributes(kmeans_num,{'numeric'},{'scalar','positive'}); validateattributes(MAP_iter,{'numeric'},{'scalar','positive'}); validateattributes(EM_iter,{'numeric'},{'scalar','positive'}); end kmeans_num = int64(kmeans_num); MAP_iter = int64(MAP_iter); EM_iter = int64(EM_iter); % Apply initialization parameters. newcls.MAP_iter = MAP_iter; newcls.EM_iter = EM_iter; newcls.kmeans_num = kmeans_num; newcls.imfilterArgs = {fspecial('gaussian', 3*ceil(3)+1, 3),... 'replicate'}; end function [labelIMG, mu, sigma] = segmentation(obj, img) %Main function of HMRF-EM image segmentation algorithm. % [labelIMG, mu, sigma] = obj.segmentation(img) % >>> Parameters % obj: The instance of `HMRF_EM` class itself % img: File name of the image for segmentation, or a uint8(:,:,3) % image matrix. % >>> Returns % labelIMG: Final matrix with 2D labels. % mu: Final vector of means. % sigma: Final vector of standard deviations. % Read image if isa(img,'char') IMG = imread(img); elseif isa(img, 'uint8') if size(img,3)==3 IMG = img; else error('Input matrix is not an RGB image matrix!'); end else error('Invalid input!'); end % Generate gray image grayIMG = rgb2gray(IMG); % Generate edge image edgeIMG = edge(grayIMG, obj.edgeArgs{:}); % imwrite(edgeIMG*255, 'edge_image.png'); % Use filter to blur `grayIMG`. % NOTE: The new `grayIMG` is a double matrix. Use % `imshow(uint8(grayIMG))` to show it. grayIMG = double(grayIMG); grayIMG = imfilter(grayIMG, obj.imfilterArgs{:}); % imwrite(uint8(grayIMG), 'blurred_gray_image.png'); % K-means algorithm. [labelIMG, mu, sigma] = obj.imgKmeans(grayIMG); % imwrite(uint8(labelIMG*ceil(255/obj.kmeans_num)), ... % 'init_label_image.png'); % sum_U : record value of U in each iteration sum_U = zeros(obj.EM_iter, 1); % Main loop for HMRF-EM algorithm. for iternum = 1:obj.EM_iter % Update `labelIMG` [labelIMG, sum_U(iternum)] = obj.MRF_MAP(labelIMG, grayIMG, edgeIMG, ... mu, sigma); % Calculate P_lyi P_lyi = obj.get_P_lyi(labelIMG, grayIMG, edgeIMG, mu, sigma); % Update mu and sigma for L=1:obj.kmeans_num sum_P_lyi = sum(P_lyi(:,L)); mu(L) = sum(P_lyi(:,L) .* grayIMG(:)) ./ sum_P_lyi; sigma(L) = sqrt(sum(P_lyi(:,L).*(grayIMG(:)-mu(L,1)).^2)/sum_P_lyi); end % Judge whether the loop can be ended. if iternum>=3 && std(sum_U(iternum-2:iternum))/sum_U(iternum)<1e-4 break; end end % imwrite(uint8(labelIMG .* ceil(255/obj.kmeans_num)), 'final_label_image.png'); end end methods (Access=protected) function [labelIMG, sum_U] = MRF_MAP(obj, labelIMG, grayIMG, edgeIMG, ... mu, sigma) % Algorithm for the MAP estimation. % >>> Parameters % obj: The `HMRF_EM` class itself % labelIMG: Matrix with 2D labels. % grayIMG : A doubly-type matrix. Use `uint8(grayIMG)` to convert % into gray image. % edgeIMG : A uint8-type matrix, generated by `edge` function. % mu: Vector of means. `mu(i,1)` is the mean value of the i-th label. % sigma: Vector of standard deviations. `sigma(i,1)` is the standard % deviations of the i-th label. % >>> Returns % labelIMG: Matrix with 2D labels. % sum_U : Sum of prior energy function. % X, Y, Z: labelIMG, grayIMG, edgeIMG % Get size of image as [m, n] [m, n] = size(grayIMG); % Iterations k = obj.kmeans_num; sum_U_MAP = zeros(obj.MAP_iter, 1); % `U` is a matrix which records prior energy function of each cluster. % Value of each column means energy function of each cluster. U = zeros(m*n, k); % Calculate U1: a part of U % $U1 = U(y|x, \Theta) = \sum_{i} U(y_i|x_i, \Theta) = \sum_{i} % [\frac{(y_i - \mu_{xi})^2}{2 \sigma_{xi}^2} - ln(\simga_{xi})]$ U1 = zeros(m*n, k); Y = grayIMG(:); % Flatten `grayIMG` for L = 1:k U1(:,L) = (Y - mu(L,1)).^2./(2 * sigma(L,1)^2) + log(sigma(L,1)^2); end for iternum = 1:obj.MAP_iter U2 = zeros(m*n, k); % Inner iteration for all labels. for L = 1:k % Loop for all pixels and calculate U2 % $U2 = U(x) = \sum_{c \in C} V_c(x)$ % $V_c(x)$ is the clique potential and C is the set of all possible % cliques. Assume that one pixel has at most 4 neighbors. Then: % $V_{c}(x_i,x_j) = 0.5 * \delta_{ij}$ % where: % $\delta_{ij} = 0 (x_i = x_j)$ and $\delta_{ij} = 1 (x_i \neq x_j)$ for i = 1:m for j = 1:n ind = (j-1)*m+i; % Flatten index.(Column priority) % A temporary variable to calculate sum of U2. u2 = 0; % Value of labelIMG(i,j) label_ij = labelIMG(i,j); % calculate U2 at (i,j) if (i-1>=1) && edgeIMG(i-1,j)==0 u2 = u2 + (label_ij ~= labelIMG(i-1,j))/2; end if (i+1<=m) && edgeIMG(i+1,j)==0 u2 = u2 + (label_ij ~= labelIMG(i+1,j))/2; end if (j-1>=1) && edgeIMG(i,j-1)==0 u2 = u2 + (label_ij ~= labelIMG(i,j-1))/2; end if (j+1<=n) && edgeIMG(i,j+1)==0 u2 = u2 + (label_ij ~= labelIMG(i,j+1))/2;