【图像分割】基于matlab 2D水平集三维医学图像分割【含Matlab源码 584期】

一、获取代码方式

获取代码方式1：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

获取代码方式2：
完整代码已上传我的资源：【图像分割】基于matlab 2D水平集三维医学图像分割【含Matlab源码 584期】

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、图像分割简介

理论知识参考：【基础教程】基于matlab图像处理图像分割【含Matlab源码 191期】

三、部分源代码

unction [phi L_z] = ls_sparse(phi, L_z, h, iter)
  % Whitaker "A Level-Set Approach to 3D Reconstruction..." IJCV 1998
  CFL = 0.15;
  % shift and peek operations
  peekD = @(M,rr,cc) M(safe_sub2ind(size(phi), rr-1, cc));
  peekU = @(M,rr,cc) M(safe_sub2ind(size(phi), rr+1, cc));
  peekR = @(M,rr,cc) M(safe_sub2ind(size(phi), rr,   cc-1));
  peekL = @(M,rr,cc) M(safe_sub2ind(size(phi), rr,   cc+1));
  shiftD = @(M) M([1 1:end-1],:);
  shiftL = @(M) M(:,[2:end end]);
  shiftR = @(M) M(:,[1 1:end-1]);
  shiftU = @(M) M([2:end end],:);
  
  %- gather layers
  Nd = shiftD(phi); Nu = shiftU(phi); Nr = shiftR(phi); Nl = shiftL(phi);
  L_i = find(phi < -.5 & (Nu >= -.5 | Nd >= -.5 | Nl >= -.5 | Nr >= -.5));
  L_o = find(phi >  .5 & (Nu <=  .5 | Nd <=  .5 | Nl <=  .5 | Nr <=  .5));
  %- layer reverse lookup
  L_rev = 2*ones(size(phi), 'int8'); % default: outside
  L_rev(phi < 0) = -2; % inside
  L_rev(L_i) = -1;
  L_rev(L_o) =  1;
  L_rev(L_z) =  0;
  
  for i = 1:iter
    % record sign
    L_z_ = L_z;
    L_z_phi_ = phi(L_z);

    %- 1. update value of zero layer (with uncontrollable force)
    F = h.init_iteration(phi, L_z);
    phi(L_z) = phi(L_z) - CFL*F/max(abs(F)); % control that force

    %- 3. update status of zero layer
    S_i = L_z(phi(L_z) < -.5     );
    S_o = L_z(      .5 < phi(L_z));
    L_z(phi(L_z) < -.5 | .5 < phi(L_z)) = []; % drop
    
    %- 2. update value of nonzero layers
    [rr cc] = ind2sub(size(phi), L_i);
    N = [peekR(phi,rr,cc) peekL(phi,rr,cc) peekU(phi,rr,cc) peekD(phi,rr,cc)]';
    phi(L_i) = max(N) - 1;
    [rr cc] = ind2sub(size(phi), L_o);
    N = [peekR(phi,rr,cc) peekL(phi,rr,cc) peekU(phi,rr,cc) peekD(phi,rr,cc)]';
    phi(L_o) = min(N) + 1;

    %- 4. update status of nonzero layers
    S_zi = L_i(     -.5 <= phi(L_i));
    S_zo = L_o(phi(L_o) <= .5      );
    L_rev(L_i(phi(L_i) < -1.5)) = -2;
    L_rev(L_o(phi(L_o) >  1.5)) =  2;
    L_i(phi(L_i) < -1.5 | -.5 <= phi(L_i)) = [];
    L_o(phi(L_o) <=  .5 | 1.5 <  phi(L_o)) = [];

    %- 5. update reverse lookup
    L_rev(S_i) = -1;
    L_rev(S_o) =  1;
    L_rev([S_zi; S_zo]) = 0;
    
    %- 7. add new neighbors
    Nl = shiftL(L_rev); Nr = shiftR(L_rev);
    Nu = shiftU(L_rev); Nd = shiftD(L_rev);
    is_near = Nl==0 | Nr==0 | Nu==0 | Nd==0;
    S_oo = find(L_rev ==  2 & is_near);
    S_ii = find(L_rev == -2 & is_near);
    L_o = [L_o; S_o; S_oo];
    L_i = [L_i; S_i; S_ii];
    L_z = [L_z; S_zi; S_zo];
    L_rev(S_ii) = -1;
    L_rev(S_oo) =  1;

    %- 8. update values of new neighbors
    [rr cc] = ind2sub(size(phi), S_ii);
    N = [peekR(phi,rr,cc) peekL(phi,rr,cc) peekU(phi,rr,cc) peekD(phi,rr,cc)]';
    phi(S_ii) = max(N) - 1;
    [rr cc] = ind2sub(size(phi), S_oo);
    N = [peekR(phi,rr,cc) peekL(phi,rr,cc) peekU(phi,rr,cc) peekD(phi,rr,cc)]';
    phi(S_oo) = min(N) + 1;
    
    %- 9. update bookkeeping
    L_z_phi = phi(L_z_);
    idx_out = L_z_(L_z_phi_ >= 0  & L_z_phi < 0);
    idx_in  = L_z_(L_z_phi_ < 0 & L_z_phi >= 0);
    h.move_in(idx_in);
    h.move_out(idx_out);
  end
end
function [phi L_z] = ls_sparse3(phi, L_z, h, iter,img)
  % Whitaker "A Level-Set Approach to 3D Reconstruction..." IJCV 1998
  CFL = .4;
  % shift and peek operations
  peekD = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr-1, cc,   dd));
  peekU = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr+1, cc,   dd));
  peekR = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr,   cc-1, dd));
  peekL = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr,   cc+1, dd));
  peekF = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr,   cc,   dd+1));
  peekB = @(M,rr,cc,dd) M(safe_sub2ind(size(phi), rr,   cc,   dd-1));
  shiftD = @(M) M([1 1:end-1],:,:);
  shiftL = @(M) M(:,[2:end end],:);
  shiftR = @(M) M(:,[1 1:end-1],:);
  shiftU = @(M) M([2:end end],:,:);
  shiftF = @(M) M(:,:,[2:end end]);
  shiftB = @(M) M(:,:,[1 1:end-1]);
  
  %- gather layers
  Nd = shiftD(phi); Nu = shiftU(phi); Nr = shiftR(phi); Nl = shiftL(phi); Nf = shiftF(phi); Nb = shiftB(phi);
  L_i = find(phi < -.5 & (Nu >= -.5 | Nd >= -.5 | Nl >= -.5 | Nr >= -.5 | Nf >= -.5 | Nb >= -.5));
  L_o = find(phi >  .5 & (Nu <=  .5 | Nd <=  .5 | Nl <=  .5 | Nr <=  .5 | Nf <=  .5 | Nb <=  .5));
  %- layer reverse lookup
  L_rev = 2*ones(size(phi), 'int8'); % default: outside
  L_rev(phi < 0) = -2; % inside
  L_rev(L_i) = -1;
  L_rev(L_o) =  1;
  L_rev(L_z) =  0;
  
  for i = 1:iter
    % record sign
    L_z_ = L_z;
    L_z_phi_ = phi(L_z);

    %- 1. update value of zero layer (with uncontrollable force)
    F = h.init_iteration(phi, L_z,img);%speed1%speed2
%     F = h.init_iteration(phi, L_z);%mean throshhold..
    phi(L_z) = phi(L_z) - CFL*F/max(abs(F)); % control that force

    %- 3. update status of zero layer
    S_i = L_z(phi(L_z) < -.5     );
    S_o = L_z(      .5 < phi(L_z));
    L_z(phi(L_z) < -.5 | .5 < phi(L_z)) = []; % drop
    
    %- 2. update value of nonzero layers
    [rr cc dd] = ind2sub(size(phi), L_i);
    N = [peekR(phi,rr,cc,dd) peekL(phi,rr,cc,dd) ...
         peekU(phi,rr,cc,dd) peekD(phi,rr,cc,dd) ...
         peekF(phi,rr,cc,dd) peekB(phi,rr,cc,dd)];
    phi(L_i) = max(N,[],2) - 1;
    [rr cc dd] = ind2sub(size(phi), L_o);
    N = [peekR(phi,rr,cc,dd) peekL(phi,rr,cc,dd) ...
         peekU(phi,rr,cc,dd) peekD(phi,rr,cc,dd) ...
         peekF(phi,rr,cc,dd) peekB(phi,rr,cc,dd)];
    phi(L_o) = min(N,[],2) + 1;

    %- 4. update status of nonzero layers
    S_zi = L_i(     -.5 <= phi(L_i));
    S_zo = L_o(phi(L_o) <= .5      );
    L_rev(L_i(phi(L_i) < -1.5)) = -2;
    L_rev(L_o(phi(L_o) >  1.5)) =  2;
    L_i(phi(L_i) < -1.5 | -.5 <= phi(L_i)) = [];
    L_o(phi(L_o) <=  .5 | 1.5 <  phi(L_o)) = [];

    %- 5. update reverse lookup
    L_rev(S_i) = -1;
    L_rev(S_o) =  1;
    L_rev([S_zi; S_zo]) = 0;
    
    %- 7. add new neighbors
    Nl = shiftL(L_rev) == 0; Nr = shiftR(L_rev) == 0; Nu = shiftU(L_rev) == 0;
    Nd = shiftD(L_rev) == 0; Nf = shiftF(L_rev) == 0; Nb = shiftB(L_rev) == 0;
    is_near = Nl | Nr | Nu | Nd | Nf | Nb;
    S_oo = find(L_rev ==  2 & is_near);
    S_ii = find(L_rev == -2 & is_near);
    L_o = [L_o; S_o; S_oo];
    L_i = [L_i; S_i; S_ii];
    L_z = [L_z; S_zi; S_zo];
    L_rev(S_ii) = -1;
    L_rev(S_oo) =  1;

    %- 8. update values of new neighbors
    [rr cc dd] = ind2sub(size(phi), S_ii);
    N = [peekR(phi,rr,cc,dd) peekL(phi,rr,cc,dd) ...
         peekU(phi,rr,cc,dd) peekD(phi,rr,cc,dd) ...
         peekF(phi,rr,cc,dd) peekB(phi,rr,cc,dd)];
    phi(S_ii) = max(N,[],2) - 1;
    [rr cc dd] = ind2sub(size(phi), S_oo);
    N = [peekR(phi,rr,cc,dd) peekL(phi,rr,cc,dd) ...
         peekU(phi,rr,cc,dd) peekD(phi,rr,cc,dd) ...
         peekF(phi,rr,cc,dd) peekB(phi,rr,cc,dd)];
    phi(S_oo) = min(N,[],2) + 1;
    
    %- 9. update bookkeeping
    L_z_phi = phi(L_z_);
    idx_out = L_z_(L_z_phi_ > 0  & L_z_phi <= 0);
    idx_in  = L_z_(L_z_phi_ <= 0 & L_z_phi > 0);
    h.move_in(idx_in);
    h.move_out(idx_out);
  end
end

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.
[5]赵勇,方宗德,庞辉,王侃伟.基于量子粒子群优化算法的最小交叉熵多阈值图像分割[J].计算机应用研究. 2008,(04)

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/115023106