我写了一个自动识别水中絮体并计算絮体体积、分形维数的算法

貌似直接上代码，虽然显得本人实诚，但是并不被CSDN的算法看中，看来我还要再写一些废话。
好了
这个代码是直接放在matlab编译器里就可以写的
只不过需要你提前存放一张水体图片。
大家就用我这个吧。

clear;clc;
%%xt.jpg、xt.png都是我们的原始絮体图片的名字
img = imread('xt.jpg');%读取jpg格式文件，存为img空间
figure,imshow(img);%在新的一张figure中显示img
img = im2bw(img);%将img转为二进制黑白文件，存为img空间。
original_picture=imread('xt.png');%读取png格式文件，存为original_picture
img2=im2gray(original_picture);%将original_picture转为灰度图，存为img2空间
[m,n] = size(img2);%这里获取灰度图的像素的行列总数
%%
% 进行最大灰度、最小灰度计算。总的思想就是遍历所有img2像素，并两两作比较
max = 0;
min = 255;
sum = 0;
avg = 0;
for i=1:1:m
    for j=1:1:n
        if (img2(i,j) > max)
            max = img2(i,j);%最大灰度
        end
        if (img2(i,j) < min)
            min = img2(i,j);%最小灰度
        end
        sum = sum + double(img2(i, j));
    end
end
avg = sum/(m*n);%平均灰度

figure; imshow(img2);
%创建一个新的窗口，将img2绘制。
figure,imshow(img);
%创建一个新的窗口，将二进制黑白的图像img绘制。
[B,L] = bwboundaries(img);%利用该函数获取絮体的各个轮廓，其中B返回的每个絮体，L是每个絮体的坐标
%%
%对每个絮体进行标记
hold on;
for k = 1:length(B)%对絮体进行遍历
    boundary = B{k};
    plot(boundary(:,2),boundary(:,1),'g','LineWidth',2);
end
[L,N] = bwlabel(img);
img_rgb = label2rgb(L,'hsv',[.5 .5 .5],'shuffle');%给每个絮体生成一个颜色
figure,imshow(img_rgb);%绘制絮体区域
hold on
for k =1:length(B)
    boundary = B{k};
    plot(boundary(:,2),boundary(:,1),'w','LineWidth',2);%绘制絮体的边界：w代表white白色
    text(boundary(1,2),boundary(1,1),num2str(k),'Color','y','Fontsize',14,'FontWeight','bold');%给每个絮体进行编号
end
hold off;
%%
%输出每个标签的面积与周长
for i=1:length(B)%同样遍历每个絮体B
    sumL(i)=0;%初始化面积
    sumC(i)=0;%初始化周长
    %以下初始化x、y方向的长度
    minx_lenth(i)=100000;
    maxx_lenth(i)=0;
    miny_lenth(i)=100000;
    maxy_lenth(i)=0;
    %因为L存储的是絮体的像素坐标，用两层for遍历像素坐标
    for k =1:length(L(:,1))%k是行数
        for j=1:length(L(1,:))%j是列数
            if L(k,j)==i%如果改坐标对应的标签L是i， sumL(i)代表i标签对应的絮体的面积就+1
                sumL(i)=sumL(i)+1;
                %%%%%%%%%%%%%%%%%%%%判断i标签的最大长度，最大长度主要是用来计算分形维数
                if k<minx_lenth(i)
                    minx_lenth(i)=k;
                end
                if k>maxx_lenth(i)
                    maxx_lenth(i)=k;
                end
                if j<miny_lenth(i)
                    miny_lenth(i)=j;
                end
                if j>maxy_lenth(i)
                    maxy_lenth(i)=j;
                end
                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            end
        end
    end
    if maxy_lenth(i)-miny_lenth(i)>maxx_lenth(i)-minx_lenth(i)%判断哪个方向是絮体的长度方向。
        maxL(i)=maxy_lenth(i)-miny_lenth(i);%各个絮体最大长度
    else
        maxL(i)=maxx_lenth(i)-minx_lenth(i);%各个絮体最大长度
    end
    equi_diameter(i)=sqrt(4*sumL(i)/pi);%这里计算的是各个标签i的等效直径
    [row,col]=size(B{i});%获取B的行数、列数
    for t=1:row-1%利用两点的距离完成轮廓的周长计算
        delta_x=abs(B{i}(t,1)-B{i}(t+1,1));
        delta_y=abs(B{i}(t,2)-B{i}(t+1,2));
        sumC(i)=sumC(i)+ sqrt(delta_y.^2+delta_x.^2);%sumC是絮体的周长
    end
end

%%
% 计算分形维数
for i=1:length(B)
    DF1(i)=log(sumL(i))./log(10);
    DF2(i)=log(maxL(i))./log(10);
    DF(i)= DF1(i)./DF2(i);
end
%%
%变量说明：DF(i)—— 分形维数；sumC(i)——絮体周长； maxL(i)——絮体最大长度； sumL(i)————絮体面积；avg
%——图像平均灰度；max ————最大灰度；min ————最小灰度