MATLAB进行图像分割问题

wanglegend

做一个项目，主要是分析污染物在土中的扩散速率，已经拍出很多照片。希望能把图片中的污染物的轮廓分离出来。我已经用matlab做了很长时间，增强对比度，去噪，提取轮廓，但总是得不出很好的结果。matlab的例题和网上的很多资料都是只对他们的照片管用，用在这上面效果不理想。这应该是比较简单的问题，可是临时学图像处理，弄了半天，实在很费劲。请求大侠帮忙指导了！不胜感激！
一张样图见附件。

jzfan

我想这个应该根据污染物在照片中的特性,采用合适的廓提取算法,而不能按照按照书中的算法照搬,比如你这张照片,分析颜色特征就比较好

wanglegend

大侠能否进一步支教一二，感激不尽！

alexqxp

简单做了一下。
1。RGB中G通道分离；
2。平滑；
3。双峰阈值。

我是用PS做的。MATLAB操作类似。

wanglegend

万分感谢！
我我是学水利的，学matlab图像处理才三四天，完全是分析问题需要。
我已经尝试先将其转化为灰度图像，再增强对比度，再滤波，再提取边界，但是增强对比度、滤波后背景仍然有很多无用的图形。呵呵，可能是因为是色弱的缘故，没有想起来从开头就考虑颜色的影响，只是按照例题先转为灰度图像了。
我需要处理几十张相似图片，因为要求其扩散速率。而且原始图像还有周围的大背景，如附件所示。那么我这么计算行不行，先检测污染物目标的R、G、B各自的大致范围，然后对图像进行处理，将象素点R、G、B值在这个范围内的保留，去掉所有在区间外的象素点。这样实现可以吗？能否得到污染物目标的图像？实现起来会不会很难呢？
多谢！

alexqxp

你的思路可行，具体实现参考一下

http://www.simwe.com/forum/viewthread.php?tid=612042

大头

如果你在拍摄图像的时候，背景是一致的，可以尝试先拍摄初始图像，然后将随后拍摄的图像与之相比较，这样的方法提取图像。

wanglegend

说到这我只想哭。这是我师兄做的实验，自己花了100多块买了个摄像头放在离心机中间随着试验箱一起转，试验做了几天，拍了上百张照片，早上、中午、夜里的照片都有，背景光照很不一样，摄像头又很糙。所以照片的背景还是有一些区别的。简单的相减不是很容易。
导师给的任务是希望能作出即时分析的软件来，就这试验条件，真是能折腾人。

wanglegend

感谢“大头”网友的提示。本来觉得由于照片数量非常多，且试验时间很长，照片背景光照很不一致，以为无法用普通的图片相减的方法取得大致的目标区域的图像。今天仔细看了全部图片，发现虽然光照不一致，但图像范围是一致的，即拍摄过程中摄像头没有移动，照片中的物体的位置没有变化，突然想起PS中的mask，于是用PS做了一个mask，目标大致区域RGB都取0，即黑色，背景区域RGB都取255，即白色，将该mask与其他试验图片相加，即可得到范围小了很多的目标区域。再接着往下分析应该会简单些。如附图。
        现在的问题是，这样取mask能否在matlab中实现，我做的例子用椭圆截取，条件放宽一些，如果我给定原始图片中一个矩形的四个角点坐标，能否只取出这个矩形内的图像，而删除周围所有的象素，并且只保存该矩形内的图像。因为如附图所示，虽然周围图像没了，却仍然留下了大片空白，怎样只留中间的那个矩形呢？
        望高手指点！

wanglegend

mask

wanglegend

最后结果

wanglegend

那么能否直接得到如下图片呢。
如果我给定四个角的坐标。

alexqxp

imcrop

wanglegend

感谢各位网友的建议。
      现在我已经把目标范围缩小到一个小矩形内。
       本来按照上面网友的建议将图像的G通道单独抽出作分析，抽出来的是一个灰度图，那么与RGB2GRAY函数将原图转为灰度图并没有本质上的区别。现在需要将目标从背景中提取出来。我看了几本书，上面都提到“基于灰度分割”，但是只提到两种算法。一种是简单地通过imhist命令目测阈值，再进行分割，一般只能对一张图适用，换图像又得改；一种是下面所示的算法，但如程序所示，th值仍然靠自己取定，而且一般情况下也只对一张图适用。我用该程序结算，后面几十张照片效果还好，但前面几十张th参数必须更改才能得出结果。
      使用imadjust函数增强图像对比度时可以使用stretchlim传递给imadjust灰度范围，为什么“基于灰度分割”既没有相应的函数，也没有像
stretchlim这样的参数自动获取明显有双峰的灰度图的阈值呢？
PS：matlab能否循环调用图片进行处理（图片可以按自然数命名，如1、2、3……），可是imread('')似乎只能读一张图片，不能循环调用啊？
多谢各位！

i2=imread('1.bmp');
rect =[250 65 109 150];
subplot(221);imshow(i2);title('origina');
a2=imcrop(i2,rect);
subplot(222);imshow(a2);title('cut');
greena2=a2(:,:,2);%抽取绿通道灰度图，事实上也可以直接使用原图像的灰度图
subplot(223);imshow(greena2);title('green gray');
low_high=stretchlim(greena2);%增强图像，似乎也不是一定需要
greena2=imadjust(greena2,low_high,[]);
subplot(224);imshow(greena2);title('after adjust');

count=imhist(greena2);
[r,t]=size(greena2);
n=r*t;
l=256;
count=count/n;%各级灰度出现的概率
for i=2:l
    if count(i)~=0
        st=i-1;
        break
    end
end
%以上循环语句实现寻找出现概率不为0的最小灰度值
for i=l:-1:1
    if count(i)~=0;
        nd=i-1;
        break
    end
end
%实现找出出现概率不为0的最大灰度值
f=count(st+1:nd+1);
p=st;q=nd-st;%p和分别是灰度的起始和结束值
u=0;
for i=1:q;
    u=u+f(i)*(p+i-1);
    ua(i)=u;
end
%计算图像的平均灰度值
for i=1:q;
    w(i)=sum(f(1:i));
end
%计算出选择不同k的时候，A区域的概率
d=(u*w-ua).^2./(w.*(1-w));%求出不同k值时类间方差
[y,tp]=max(d);%求出最大方差对应的灰度级
th=tp+p;
if th<=105
    th=tp+p;
else
    th=105;       %根据具体情况适当修正门限，那意味着岂不是每张图都得单独指定？
end
y1=zeros(r,t);
for i=1:r
    for j=1:t
        x1(i,j)=double(greena2(i,j));
    end
end
for i=1:r
    for j=1:t
        if (x1(i,j)>=th)
            y1(i,j)=x1(i,j);
        else
            y1(i,j)=0;
        end
    end
end
%上面一段代码实现分割
figure,imshow(y1);
title('灰度门限分割的图像');

alexqxp

1.分割的方法很多，比如用聚类分割；
2.双峰法也可以通过程序实现自动分割；
3.可以进行图像批处理，论坛里面有帖子，查查看。
下面是聚类分割，全自动的。

wanglegend

谢谢，我再查查看。希望能多找集中方法。
以下是我通过选定目标图像的RGB范围进行分割的程序。算是与大家交流吧。欢迎大家提出意见。不胜感激！

clear,clc;
i=imread('15_37_55.bmp');%读入照片14_3_55
rect =[250 65 109 150];%取定一个范围
subplot(121);imshow(i);title('origina');%显示原来的照片
a2=imcrop(i,rect);%按照范围截取图片
a3=imcrop(a2);

a31=a3(:,:,1);%取R
r1=size(a31,1);
c1=size(a31,2);
a31=reshape(a31,r1*c1,1);
a31mean=mean(a31),%求平均值
a31std=std(double(a31)),%求方差

a32=a3(:,:,2);%取G
r2=size(a32,1);
c2=size(a32,2);
a32=reshape(a32,r1*c1,1);
a32mean=mean(a32),%求平均值
a32std=std(double(a32)),%求方差

a33=a3(:,:,3);%取G
r3=size(a33,1);
c3=size(a33,2);
a33=reshape(a33,r1*c1,1);
a33mean=mean(a33),%求平均值
a33std=std(double(a33)),%求方差

%subplot(122);imshow(a2);title('cut');%截取图像的一部分

%a2(50,50,2)
for i=1:(rect(4)+1)
    for j=1:(rect(3)+1)
        if a2(i,j,1)>(a31mean-3*a31std)&a2(i,j,1)<(a31mean+3*a31std)&a2(i,j,2)>(a32mean-3*a32std)&a2(i,j,2)<(a32mean+3*a32std)&a2(i,j,3)>(a33mean-3*a33std)&(a2(i,j,3)<a33mean+3*a33std)
            a2(i,j,1)=0;%目标图像部分令其为黑色
            a2(i,j,2)=0;
            a2(i,j,3)=0;
        else
             a2(i,j,1)=255;%其他部分为白色
            a2(i,j,2)=255;
            a2(i,j,3)=255;
        end
    end
end
figure,imshow(a2);
a2=RGB2GRAY(a2);
title('色素分割的图像');

i_filter3=medfilt2(a2,[7 7]);%中值滤波，滤去噪音点
figure,imshow(i_filter3);title('窗口大小为7×7');%窗口大小为

bw1=edge(i_filter3,'canny',0.2);%提取边缘
figure,imshow(bw1);%显示边缘

wanglegend

对不起。我已经看了4本MATLAB图像处理的书，没有发现双峰法的自动实现方法（有需要人工确定阈值的）。请问能多给点指导吗？谢谢！

alexqxp

你可以参考下面的：
双峰法

双峰法的原理及其简单：它认为图像由前景和背景组成，在灰度直方图上，前后二景都形成高峰，在双峰之间的最低谷处就是图像的阈值所在。根据这一原理，我们给出了它的实现，部分代码如下（Pascal语言描述，以下同）：

//intPeak、intPeak2、intValley：峰值和直方图值

//intIndx:：相应的灰度值

intPeak,intIndx,intPeak2,intIndx2,intValley,intValleyIndx:integer;

//初始双峰值

    intPeak:=0;

    intPeak2:=0;

//取得第一峰值

    for intLoop:=0 to 255 do

      if intPeak<=intGrayLevel[intLoop] then

      begin

        intPeak:=intGrayLevel[intLoop];

        intIndx:=intLoop;

      end;

//取得第二峰值

    for intLoop:=0 to 255 do

    Begin

      if (intPeak2<=intGrayLevel[intLoop]) and (intLoop<>intIndx) then

      begin

        intPeak2:=intGrayLevel[intLoop];

        intIndx2:=intLoop;

      end

    end;

//取得双峰之间的谷值

    intValley:=intSize;

    if intIndx2<intIndx then

      for intLoop:=intIndx2 to intIndx do

        if intValley>intGrayLevel[intLoop] then

        begin

          intValley:=intGrayLevel[intLoop];

          intValleyIndx:=intLoop;

        end;

从分割的效果来看，当前后景的对比较为强烈时，分割效果较好；否则基本无效。

*********************************************************
或者用这个试试。结果就是我放上去的那张图。

clear,clc;
he=imread('15_37_55.bmp');%读入照片14_3_55
imshow(he), title('H&E image');
cform = makecform('srgb2lab');
lab_he = applycform(he,cform);
ab = double(lab_he(:,:,2:3));
nrows = size(ab,1);
ncols = size(ab,2);
ab = reshape(ab,nrows*ncols,2);

nColors = 3;
% repeat the clustering 3 times to avoid local minima
[cluster_idx cluster_center] = kmeans(ab,nColors,'distance','sqEuclidean', ...
                                      'Replicates',3)

pixel_labels = reshape(cluster_idx,nrows,ncols);
imshow(pixel_labels,[]), title('image labeled by cluster index');
segmented_images = cell(1,3);
rgb_label = repmat(pixel_labels,[1 1 3]);

for k = 1:nColors
    color = he;
    color(rgb_label ~= k) = 0;
    segmented_images{k} = color;
end

imshow(segmented_images{2}), title('objects in cluster 2');

wanglegend

bow!
太谢谢了。我得好好研究研究。
谢谢指点。我已经找到依次调用多个文件处理的方法。也放在这里。
（文件名分别是my1.bmp，my2.bmp，my3.bmp）

for j=1:3
A=['my',num2str(j),'.bmp'];
i=imread(A);
figure,imshow(i),
end

如你所说。当前景色与背景色区别较明显时，我们的方法效果都很好。但是当它们的对比不是那么强烈或是边缘较模糊时，结果真的还是不太好。
再次感谢。

wanglegend

这两天又好好研究了阈值分割的方法
刚刚弄出来，程序如下，仍然是对目标和背景区别较大的图像有用。请问
alexqxp 能提供你给的聚类方法分割程序中的
makecform和applycform等函数的代码吗？实在惭愧，研究了一段时间，如果没有你编的函数，实在还是看不懂你的程序。惭愧惭愧。
另：我查资料说可以利用imregionalmin函数来搜索直方图中的谷值，然后进行阈值分割。可是我的matlab好像有问题，居然没有imregionalmin的help，请问有人知道这个函数的用法吗？bow!

function th=thresh_md(a);
%该函数实现最大方差法计算分割门限。
%输入参数为灰度图像，输出为灰度门限。
count=imhist(a);
%返回图像矩阵a各个灰度等级像素个数
[m,n]=size(a);
N=m*n-sum(sum(find(a==0),1));
L=256;
%指定图像灰度等级为256级。
%d=std2((double(a)+1)/N)
%err=0.003
count=count/N;
%计算出各灰度出现的概率。
for i=2:L
    if count(i)~=0
        st=i-1;
        break
    end
end
%找出出现概率不为0的最小灰度
for i=L:-1:1
    if count(i)~=0;
        nd=i-1;
        break;
    end
end
%找出出现概率不为0的最大灰度
f=count(st+1:nd+1);
p=st;  q=nd-st;
%p和q分别为灰度起始值和结束值
u=0;
for i=1:q;
    u=u+f(i)*(p+i-1);
    ua(i)=u;
end;
%计算图像的平均灰度
for i=1:q;
    w(i)=sum(f(1:i));
end;
%计算出选择不同k值时，A区域的概率。
d=(u*w-ua).^2./(w.*(1-w));
%求出不同k值时类间的方差。
[y,tp]=max(d);
%求出最大方差对应的灰度级。
th=tp+p;

主程序:

close all,clear,clc;
i=imread('13_59_55.bmp');
rect =[250 65 109 150];
subplot(231);imshow(i);title('origina');
image=imcrop(i,rect);
subplot(232);imshow(image);title('after cropping');
imagegreen=image(:,:,2);%抽取绿通道灰度图，事实上也可以直接使用原图像的灰度图
subplot(233);imshow(imagegreen);title('green gray');

th=thresh_md(imagegreen);
imagebw=im2bw(imagegreen,th/255);
subplot(234);imshow(imagebw);title('after effect');

imagebw=(~imagebw);

i_filter3=medfilt2(imagebw,[7 7]);%中值滤波，滤去噪音点
subplot(235);imshow(i_filter3);title('中值滤波后');%窗口大小为7×7

bw1=edge(i_filter3,'canny',0.2);%提取边缘
subplot(236);imshow(bw1);title('提取的图像边缘');%显示边缘