matlab将孔隙空间充满毛细血管算法实现

二维情况

• 假设：毛细血管横截面是圆形
• 给定的图像是灰度图

代码实现

• 预处理

  %% binarize the origin image
  a = imread('1.jpg');
  b = im2bw(a, mean(a(:))/255);
  d = ones(4);
  c = imopen(b, d);
  imwrite(c, '2.bmp')
  

  原图和二值化后经过开运算处理的图片如下：
  

  

• 产生所有的在孔隙内部可能的圆

  %% generate circles that can be in the area
  tic
  img = imread('2.bmp');
  edgeI = edge(img);
  [rows, cols] = find(img == 0 & edgeI ~= 1 );
  [edgeR, edgeC] = find(edgeI == 1);
  edges = [ edgeR, edgeC ];
  resu = [];
  for ii = 1:length(rows)
      nowP = [rows(ii), cols(ii)];
      tmp = repmat( nowP, size(edges, 1), 1 );
      dis = sqrt( sum( (tmp-edges).^2, 2) );
      r = floor(min(dis(:)));
      resu = [resu; [nowP, r ]];
  end
  toc
  

• 删除重复或者不可取的圆

  对于有些互相相交的圆或者包围在其他圆内部的圆，是需要对其进行删除的。

  %% delete circles that are unuseful
  tic
  circles = resu;
  circles = sortrows(circles, -3);
  tplot = [];
  lpC = size(circles,1);
  for ii = 1:lpC-1
      if circles(ii,3) == 0        
          continue;
      end
      ct = repmat( circles(ii,1:2), size(circles,1)-ii, 1 );
      otherP = circles(ii+1:end, 1:2);
      dis = sqrt( sum( (ct-otherP).^2, 2) );
      iRows = find(dis < (circles(ii,3)+circles(ii+1:end,3)));
      circles(iRows+ii, :) = repmat([inf inf 0], length(iRows), 1);
  end
  [rows, cols] = find(circles(:,3) == 0 || circles(:,3) == 1); %删除不可取的圆
  circles(rows, :) = [];
  
  toc
  

• 可视化，显示孔隙和最后的圆，即毛细血管横截面

  将圆和二值化图像放在一起叠加显示

  %% draw circles to visualize the result
  cToDraw  =circles;
  imshow(img)
  hold on
  for ii = 1:size(cToDraw, 1)
      rectangle('Position',[cToDraw(ii,2)-cToDraw(ii,3),cToDraw(ii,1)-cToDraw(ii,3),2*cToDraw(ii,3),2*cToDraw(ii,3)],'Curvature',[1,1],  'EdgeColor','g')
  end
  hold off
  axis equal
  

  结果图如下：
  

• 思路来源的论文链接
  http://cdmd.cnki.com.cn/Article/CDMD-10358-2007097246.htm