22、显微镜图像中的粒度分析与木炭量化

显微镜图像中粒度与木炭量化分析

sprite

于 2025-11-23 11:22:43 发布

阅读量9

点赞数

CC 4.0 BY-SA版权

分类专栏： MATLAB地球科学实战指南文章标签：图像分析粒度分析木炭量化

本文链接：https://blog.youkuaiyun.com/sprite/article/details/155208502

MATLAB地球科学实战指南专栏收录该内容

24 篇文章 ¥499.90

订阅专栏¥69.90

会员秒杀 ¥9.9 重磅福利

超级会员免费看

显微镜图像中的粒度分析与木炭量化

1. 图像分析在地质科学中的应用概述

图像分析在地质科学领域有着广泛的应用，经典的应用包括识别、测量和计数图像中的颗粒。具体例子有粒度分析、花粉粒计数以及通过薄片确定岩石的矿物成分等。使用图像分析的主要动机在于能够在短时间内以相对较低的成本对大量样本进行自动化分析。

常见的用于识别和计数图像中物体的方法有三种：
- 基于区域的分割，使用分水岭分割算法。
- 物体检测，使用霍夫变换。
- 阈值处理，利用颜色差异分离物体。

这里将通过两个具体例子来介绍图像处理在识别、测量和计数颗粒方面的应用。

2. 粒度分析中的分水岭分割应用

2.1 数据准备与初步处理

读取图像 ：首先读取一张不同大小的粗岩屑颗粒图像，并将其存储在变量 I1 中。图像大小为 284 × 367 像素，已知宽度为 3 cm，可计算出高度为 3 cm × 284 像素 / 367 像素 = 2.32 cm。

clear
I1 = imread('grainsize.tif');
ix = 3; iy = 284 * 3 / 367;
imshow(I1,'XData',[0 ix],'YData',[0 iy])
title('Original Image')

转换为灰度图像 ：去除图像的颜色信息，使用 rgb2gray 函数将 I1 转换为灰度图像 I2 。

I2 = rgb2gray(I1);
imshow(I2,'XData',[0 ix],'YData',[0 iy])
title('Grayscale Image')

调整图像对比度 ：由于灰度图像 I2 的对比度相对较低，使用 imadjust 函数调整图像强度值，使 1% 的数据在低强度和高强度处饱和，从而增加新图像 I3 的对比度。

I3 = imadjust(I2);
imshow(I3,'XData',[0 ix],'YData',[0 iy])
title('Adjusted Intensity Values')

2.2 背景去除与二值化处理

确定背景并去除 ：使用 imopen 函数确定图像 I3 中低于一定像素大小的物体，生成背景图像，然后得到无背景图像 I4 。

I4 = imopen(I3,strel('disk',1));
imshow(I4,'XData',[0 ix],'YData',[0 iy])
title('No Background')

观察去除的背景 ：通过 imsubtract 函数从原始灰度图像 I3 中减去无背景图像 I4 ，得到被去除的背景图像 I5 。

I5 = imsubtract(I3,I4);
imshow(I5,'XData',[0 ix],'YData',[0 iy])
title('Background')

二值化处理 ：使用 im2bw 函数将无背景图像 I4 转换为二值图像 I6 ，手动调整阈值，发现 0.2 是一个合适的阈值。

I6 = im2bw(I4,0.2);
imshow(I6,'XData',[0 ix],'YData',[0 iy])
title('Binary Image')

2.3 去除边界重叠物体与物体识别

去除边界重叠物体 ：使用 imclearborder 函数消除二值图像 I6 中与图像边界重叠的物体，得到图像 I7 。

I7 = imclearborder(I6);
himage1 = imshow(I6,'XData',[0 ix],'YData',[0 iy]); hold on
set(himage1, 'AlphaData', 0.7);
himage2 = imshow(imsubtract(I6,I7),'XData',[0 ix],'YData',[0 iy]);
set(himage2, 'AlphaData', 0.4);
title('Image Border'), hold off

追踪物体边界 ：使用 bwboundaries 函数追踪二值图像 I7 中物体的边界，选择 noholes 选项抑制孔洞追踪。

[B,L] = bwboundaries(I7,'noholes');
imshow(label2rgb(L,@jet,'w','shuffle'),...
   'XData',[0 ix],'YData',[0 iy])
title('Define Objects')

计算连通物体数量 ：使用 bwlabeln 函数计算二值图像中连通物体的数量，在这个例子中，该方法识别出 192 个颗粒，明显低于手动计数的 236 个（减去 3 个与边界重叠的物体）。

[labeled,numObjects] = bwlabeln(I7,8);
numObjects

2.4 计算颗粒面积与统计信息

计算颗粒面积 ：使用 regionprops 函数计算每个颗粒的面积，然后将其转换为实际面积（单位：cm²）。

graindata = regionprops(labeled,'basic');
grainareas= [graindata(:).Area];
objectareas = 3^2 * grainareas * 367^(-2);

计算统计信息 ：计算图像中所有颗粒的最大、最小和平均面积。

max_area = max(objectareas)
min_area = min(objectareas)
mean_area = mean(objectareas)

绘制直方图 ：绘制所有颗粒面积的直方图。

clf
v = 0.0005 : 0.0005 : 0.15;
hist(objectareas,v)
xlabel('Grain Size in Millimeters^2')
ylabel('Number of Grains')
axis([0 0.1 0 30])

2.5 分水岭分割算法改进

为了部分克服粒度分析中连通颗粒的问题，采用分水岭分割算法。该算法将灰度值解释为数字高程模型中的高程，通过计算最近邻距离来分割物体。
- 计算距离变换 ：由于在这个例子中物体像素值为 1，背景像素值为 0，因此对二值图像 I7 的补集应用 bwdist 函数计算距离变换。

D = bwdist(~I7,'cityblock');

处理距离矩阵 ：对距离变换结果取反，并将不属于物体的像素值设为 -Inf。

D = -D;
D(~I7) = -Inf;

计算分水岭变换 ：对距离矩阵进行分水岭变换，得到标签矩阵 L2 ，并显示结果。

L2 = watershed(D); 
imshow(label2rgb(L2,@jet,'w','shuffle'),...
   'XData',[0 ix],'YData',[0 iy])
title('Watershed Segmentation')

计算物体像素数量与面积 ：计算每个物体的像素数量，然后转换为实际面积（单位：cm²），再次计算最大、最小和平均面积，并绘制直方图。

objects = sortrows(L2(:),1);
max(objects)
clear objectsizes
for i = 2 : max(objects)
    clear individualobject
    individualobject = objects(objects == i);
    objectsizes(i) = length(individualobject);
end
objectsizes = objectsizes';
objectsizes = sortrows(objectsizes,1);
objectsizes = objectsizes(objectsizes~=0);
objectareas = 3^2 * objectsizes * 367^(-2);
max_area = max(objectareas)
min_area = min(objectareas)
mean_area = mean(objectareas)
clf
v = 0.0005 : 0.0005 : 0.15;
hist(objectareas,v)
xlabel('Grain Size in Millimeters^2'), ylabel('Number of Grains')
axis([0 0.1 0 70])

2.6 结果验证

为了验证最终结果，对一个较大颗粒的轮廓进行数字化，存储在变量 data 中，然后使用 polyarea 函数计算该颗粒的面积，计算结果与粒度分析结果大致相符。

figure
imshow(I1,'XData',[0 ix],'YData',[0 iy])
data = ginput;
data(end+1,:) = data(1,:)
imshow(I1,'XData',[0 ix],'YData',[0 iy]), hold on
plot(data(:,1),data(:,2)), hold off
polyarea(data(:,1),data(:,2))

3. 霍夫变换相关资源

对于使用霍夫变换检测物体（如花粉粒）感兴趣的读者，可以参考陶鹏（马里兰大学）的优秀算法。相关的 MATLAB 例程 circularhough_grd.m 和 hough_grd.m 可以从以下链接下载：
http://www.mathworks.com/matlabcentral/fileexchange/

circularhough_grd.m 在 2008 年 5 月 23 日被选为 MATLAB Central 本周精选，该博客展示了如何使用霍夫变换检测图像中的圆圈以计数红细胞：
http://blogs.mathworks.com/pick/2008/05/23/detecting-circles-in-an-image/

4. 粒度分析流程总结

graph TD;
    A[读取图像] --> B[转换为灰度图像];
    B --> C[调整图像对比度];
    C --> D[确定并去除背景];
    D --> E[二值化处理];
    E --> F[去除边界重叠物体];
    F --> G[追踪物体边界与识别物体];
    G --> H[计算颗粒面积与统计信息];
    H --> I[分水岭分割算法改进];
    I --> J[结果验证];

5. 相关软件资源

RADIUS 软件：由美因茨大学的 Klemens Seelos 开发的基于 MATLAB 的软件，实现了上述粒度分析和木炭量化的应用。该软件是一种基于薄片数字图像评估的粒度测量技术，具有亚毫米级的样本分辨率，可用于研究中粉砂至粗砂粒径范围内的沉积过程。其 MATLAB 代码可从以下链接下载：
http://www.particle-analysis.info/

通过以上内容，我们详细介绍了图像处理在粒度分析中的应用，包括基本的图像处理步骤、分水岭分割算法的应用以及相关软件和资源。下一部分将介绍如何使用阈值处理方法量化显微镜图像中的木炭含量。

显微镜图像中的粒度分析与木炭量化

6. 显微镜图像中木炭的量化

在地质科学中，量化物质的组成是一项重要任务，例如确定岩石薄片中的矿物成分，或者测量筛分沉积物样本中的木炭含量。图像处理方法为这些任务提供了便利，其中阈值处理是一种简单有效的分割图像中具有不同颜色或灰度值物体的方法。

6.1 图像读取与预处理

以肯尼亚纳库鲁湖的筛分湖底沉积物样本图像为例，该图像尺寸为 1500×1500 像素，假设图像为正方形，宽和高均为 1 厘米。

clear
I1 = imread('lakesediment.jpg');
ix = 1; iy = 1;
imshow(I1,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeter'), ylabel('Centimeter')
title('Original Image')

将 RGB 彩色图像转换为灰度图像，以简化后续处理。

I2 = rgb2gray(I1);
imshow(I2,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('Grayscale')

由于图像对比度较低，使用 imadjust 函数调整图像强度值，使 1% 的数据在低强度和高强度处饱和，从而增强图像对比度。

I3 = imadjust(I2);
imshow(I3,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('Better Contrast')

6.2 背景去除

使用 imopen 函数确定图像中低于一定像素大小的物体，以去除背景纹理，得到无背景图像 I4 。

I4 = imopen(I3,strel('disk',5));
imshow(I4,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('W/O Background')

通过 imsubtract 函数从原始灰度图像 I3 中减去无背景图像 I4 ，可以观察到被去除的背景图像 I5 。

I5 = imsubtract(I3,I4);
imshow(I5,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('Background')

6.3 阈值处理分离木炭颗粒

使用 im2bw 函数将无背景图像 I4 转换为二值图像，通过手动调整阈值，发现 0.03 是分离黑色木炭颗粒的合适阈值。

I6 = im2bw(I4,0.03);
imshow(I6,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('Only Charcoal')

此时，可以通过计数这些颗粒的像素数量来估计图像中木炭的总量。

6.4 确定所有颗粒并计算木炭百分比

使用另一个阈值 0.6，将图像中的所有物体（包括木炭和其他物质）分离出来，得到二值图像 I7 。

I7 = im2bw(I4,0.6);
imshow(I7,'XData',[0 ix],'YData',[0 iy]), axis on
xlabel('Centimeters'), ylabel('Centimeters')
title('All Objects')

我们关注的是样本中木炭的百分比，通过以下代码计算：

100*sum(sum(I6==0))/sum(sum(I7==0))

计算结果表明，筛分样本中约 13% 为木炭。通过选择不同的阈值，还可以进一步量化样本中的其他成分，如介形虫或矿物颗粒。

7. 木炭量化流程总结

graph TD;
    A[读取图像] --> B[转换为灰度图像];
    B --> C[调整图像对比度];
    C --> D[确定并去除背景];
    D --> E[第一次阈值处理分离木炭];
    E --> F[第二次阈值处理分离所有物体];
    F --> G[计算木炭百分比];

8. 阈值处理方法总结

处理步骤	操作函数	说明
读取图像	`imread`	读取指定路径的图像文件
转换为灰度图像	`rgb2gray`	将 RGB 彩色图像转换为灰度图像
调整图像对比度	`imadjust`	调整图像强度值，增强对比度
确定并去除背景	`imopen` 、 `imsubtract`	使用结构元素确定背景并去除
二值化处理	`im2bw`	通过设定阈值将图像转换为二值图像
计算木炭百分比	`sum`	通过计数像素数量计算木炭在样本中的百分比