%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name: gmac = global minimization of the active contour model
% Description: see paper "Fast Global Minimization of the Active Contour/Snake Model" in JMIV07
% Author: Xavier Bresson (xbresson@math.ucla.edu)
% Lastest version: 07-09-21
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function gmac(Im0,std_Gb,beta,region_descriptor,theta,lambda,dt,std_pwm,NbIterUpdate)
close all;
% RD = Region Descriptor
RD_MEAN = 0; % RD = Mean Intensity
RD_PDF = 1; % RD = Probability Density Function (PDF)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Processing input image
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Ny,Nx,Nc] = size(Im0);
if Nc>1; Im0=rgb2gray(Im0); end; % Convert color image to gray-scale image
Im0 = double(Im0);
Im0 = Im0/ max(Im0(:)); % Intensity range = [0 1]
maxIm = 255; % Maximum intensity value
Im0 = 1 + maxIm* Im0; % Intensity range = [1 256]
figure(1); imagesc(Im0); colormap(gray); title('Input image'); colorbar;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computing the edge detector function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Smoothing input image
Ng=ceil(6*std_Gb)+1; Gaussian = fspecial('gaussian',[Ng Ng],std_Gb);
Im0s = conv2(Im0,Gaussian,'same');
% Norm of gradient of the smoothed image
[Gx,Gy] = gradient(Im0s);
NormGrad = sqrt(Gx.^2 + Gy.^2);
% Edge detector function
Gb = 1./ (1 + beta* NormGrad.^2);
figure(2); imagesc(Gb); colormap(gray); title('The stopping function'); colorbar;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Computing 1-D Gaussian Function for PDF
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if region_descriptor == RD_PDF
Ng=ceil(6*std_pwm); if (mod(Ng,2)==0); Ng=Ng+1; end; Nc=round(Ng/2);
z=1:Ng; G1D = exp(-(z-Nc).^2/(2*std_pwm^2))/ sqrt(2*pi*std_pwm^2);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pxu = zeros(size(Im0));
pyu = zeros(size(Im0));
regionTermIn = zeros(size(Im0));
regionTermOut = zeros(size(Im0));
u = Im0/max(Im0(:)); % u between 0 and 1.
v = u;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Iterative segmentation process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for iteration=1:10000
iteration
% Update u
Divp = ( BackwardX(pxu) + BackwardY(pyu) );
Term = Divp - v/ theta;
TermX = ForwardX(Term);
TermY = ForwardY(Term);
Norm = sqrt(TermX.^2 + TermY.^2);
Denom = Gb + dt*Norm;
pxu = Gb.* (pxu + dt*TermX)./Denom;
pyu = Gb.* (pyu + dt*TermY)./Denom;
u = v - theta* Divp;
% Computing region-based term every "NbIterUpdate" iterations
if iteration==1 | rem(iteration,NbIterUpdate)==0
% MEAN
if region_descriptor == RD_MEAN
MeanIn = mean(mean(Im0(u>=0.5)));
MeanOut = mean(mean(Im0(u<0.5)));
RegionTerm = ((Im0-MeanIn).^2 - (Im0-MeanOut).^2);
MinRegionTerm = lambda* theta *min(RegionTerm(:))
MaxRegionTerm = lambda* theta *max(RegionTerm(:))
end
% PDF
if region_descriptor == RD_PDF
% Computing pdf inside the active contour
ImInside=Im0; ImInside(u>=0.5)=maxIm+1; NbPtIn=sum(sum(u<0.5));
H=reshape(ImInside,1,Nx*Ny); X=0:maxIm; H=histc(H,X); H=H/(Nx*Ny*NbPtIn);
PdfIn=conv(H,G1D); PdfIn=PdfIn(Nc:Nc+maxIm);
% Compute pdf outside the active contour
ImOutside=Im0; ImOutside(u<0.5)=maxIm+1; NbPtOut=sum(sum(u>=0.5));
H=reshape(ImOutside,1,Nx*Ny); X=0:maxIm; H=histc(H,X); H=H/(Nx*Ny*NbPtOut);
PdfOut=conv(H,G1D); PdfOut=PdfOut(Nc:Nc+maxIm);
for ct=0:maxIm;
regionTermIn(Im0==ct)=PdfIn(ct+1);
regionTermOut(Im0==ct)=PdfOut(ct+1);
end;
% Compute region term = difference between inside and outside pdfs
RegionTerm = regionTermIn - regionTermOut;
MinRegionTerm = lambda* theta *min(RegionTerm(:))
MaxRegionTerm = lambda* theta *max(RegionTerm(:))
end
end
% Update v
v = u - lambda* theta * RegionTerm;
v = max(v,0);
v = min(v,1);
% Displaying result every "NbIterUpdate" iterations
if rem(iteration,NbIterUpdate)==0
if region_descriptor == RD_PDF
figure(10); clf; plot(X,PdfIn,'r'); hold on; plot(X,PdfOut,'b'); title('Red plot: Inside PDF, Blue plot: Outside PDF.');
end
figure(11); clf;
subplot(1,3,1); imagesc(Im0); colormap(gray); title('Given Image'); hold on;
[c,h] = contour(u,[0.5 0.5],'w-'); set(h,'LineWidth',2);
subplot(1,3,2); imagesc(u); colormap(gray); title('u'); hold on;
[c,h] = contour(u,[0.5 0.5],'m-'); set(h,'LineWidth',2);
subplot(1,3,3); imagesc(v); colormap(gray); title('v');
s=1.5; truesize(11,[round(s*Ny) round(s*Nx)]);
pause(0.1)
%pause
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of main function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sub-functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dx]=BackwardX(u);
[Ny,Nx] = size(u);
dx = u;
dx(2:Ny-1,2:Nx-1)=( u(2:Ny-1,2:Nx-1) - u(2:Ny-1,1:Nx-2) );
dx(:,Nx) = -u(:,Nx-1);
function [dy]=BackwardY(u);
[Ny,Nx] = size(u);
dy = u;
dy(2:Ny-1,2:Nx-1)=( u(2:Ny-1,2:Nx-1) - u(1:Ny-2,2:Nx-1) );
dy(Ny,:) = -u(Ny-1,:);
function [dx]=ForwardX(u);
[Ny,Nx] = size(u);
dx = zeros(Ny,Nx);
dx(1:Ny-1,1:Nx-1)=( u(1:Ny-1,2:Nx) - u(1:Ny-1,1:Nx-1) );
function [dy]=ForwardY(u);
[Ny,Nx] = size(u);
dy = zeros(Ny,Nx);
dy(1:Ny-1,1:Nx-1)=( u(2:Ny,1:Nx-1) - u(1:Ny-1,1:Nx-1) );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of sub-function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
快速全局最小化活动轮廓图像分割源代码by--Xavier Bresson
4星 · 超过85%的资源 需积分: 10 74 浏览量
2010-08-20
17:09:56
上传
评论
收藏 53KB ZIP 举报
gloriachan824
- 粉丝: 0
- 资源: 3
最新资源
- 基于qt+c++实现ddos小工具可用于网站压测等性能测试+源码(期末大作业&课设&项目开发)
- docker-compose-linux-x86-64-v2.26.0.zip
- (5积分下载)学生信息管理系统-C++ 容器版本
- 3341L-VB一款P-Channel沟道SOT23的MOSFET晶体管参数介绍与应用说明
- 三菱PLC例程源码灯数任意控制FX2n
- Android算法部署-在Android平台基于NCNN部署YOLOv5目标检测算法-优质项目实战.zip
- python 链接ms sqlserver的通用文件
- 人工智能-项目实践基于TensorFlow的CNN-RNN中文文本分类源码+说明文档.zip
- 双系统安装02--在已有win10基础上安装统信UOS-CS (1).mhtml
- 体感技术kinectV2.0开发体感游戏资料,人机互动例子,手势识别,骨骼绑定,手势翻书,语音识别,包含Unity例子
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈