医学影像大数据matlab
学生实验报告
| 学生姓名 | 景涛 | 学号 | B20040301002 | 组员： |
|---|---|---|---|---|
| 实验项目 | 医学图像基础 | |||
| □必修 □选修 | □演示性实验 □验证性实验 □操作性实验 □综合性实验 | |||
| 实验地点 | E405 | 实验仪器台号 | ||
| 指导教师 | 刘巧红 | 实验日期及节次 | 2023.2-2023.6 |
一、实验目的及要求
1、目的
通过Matlab编程实现医学影像的基本读取与查看功能教学，并帮助学生全面掌握灰度直方图绘制与分析方法、各种几何运算及其应用以及傅里叶变换的基本原理与实现过程。
2、内容及要求
1、Matlab 读取DICOM
2、DICOM文件灰度直方图的显示及调整
3﹑医学影像的灰度运算与代数运算
4**、医学影像的几何运算**
5、医学影像的傅里叶变换
透彻理解上述核心知识点，并熟练运用Matlab进行基本运算。(1)学习并掌握相关理论知识及其实现代码；(2)通过实践完成灰度变换等相关技术的学习；(3)撰写基础实验报告并完成相关内容的实践作业。
二、实验工具
| 仪器名称 | 规格/型号 | 数量 | 备注 |
|---|---|---|---|
| 计算机 | Windows10 | 1 | 有网络环境 |
| Matlab | 2016b及以上 | 1 | |
| Python | 3.6及以上 | 1 |
三、实验步骤及结果
(一) Matlab 读取DICOM
（利用matlab读取dicom格式的文件，并进行多种灰度级别的显示）
clc; %clc;清空命令行窗口 clear; %clear;清空工作区 memory? close all graphics;"%关闭所有的图形界面" I = dicomread('D:\data\1\1\1\1\1\1.dcm') ; %% 解析该 DICOM 格式文件中的图像 figure( 'Name', '读取后的图像显示' );
imshow(I); %显示图片****
%设置不同的灰度级来显示图像figure; %新建画布**imshow(I,[100 2048]); %指定其灰度范围以便显示图像
figure; %新建画布****imshow(I,[512 2048]); %设置灰度级为512到2048
figure; %新建画布imshow(I,[1024 2048]); %设置灰度级为1024到2048
使用 subplot 函数创建 x×y 个子图窗口，并生成一个包含 x 行 y 列共 x*y 个子图的窗口；设置第一个子图的位置为 2 行 2 列的第一格；调用 imshow 函数显示图像 I；设置第二个子图的位置为 2 行 2 列的第二格；调用 imshow 函数并指定图像范围 [100, 2048]；设置第三个子图的位置为 2 行 2 列的第三格；调用 imshow 函数并指定图像范围 [512, 2048]；设置第四个子图的位置为 2 行 2 列的第四格；调用 imshow 函数并指定图像范围 [1024, 2048]；
imwrite(uint8(I),'1.jpg');
X = imread('1.jpg')
figure;
imshow(X);
% 查看该DICOM影像的具体信息fprintf('DICOM图像具体信息');
info = dicominfo('000000.dcm');
name = info.PatientName;
modality=info.Modality;
fprintf(modality);
（二） DICOM文件灰度分布图像的展示与调整（借助MATLAB软件对获取的医学影像数据生成灰度分布图像并实现其自动调整）
**clc;clear;I_raw = dicomread('000000.dcm');info = dicominfo('000000.dcm');%% 将DICOM格式转换成常见的图像格式% I = imadjust(I_raw,[0.5;0.55],[0;1]); %调整灰度范围I = double(I_raw); %将灰度级映射到0~255low = min(min(I)); % 先取每一列的最小值 再取一次生成的向量的最小值 也就是矩阵的最小值****high = max(max(I));maxgray = high-low; %计算窗宽rate = 256/maxgray;I = I * rate;I = I + abs(min(min(I))); %加窗I = uint8(I); %转化为8位的位图数据格式imwrite(I,'000000.jpg');****figure;****imshow(I);**title('原始DICOM图像');
该代码行用于生成原始DICOM格式医学图像的概率分布直方图，并对其进行显示。其中imhist(I)函数用于计算并绘制图像I的概率分布直方图；title('原始DICOM图像的直方图')设置了图形窗口标题以反映当前处理的内容；注释内容详细说明了对直方图均衡化的具体实现方式及其默认参数设置：即实现直方图均衡化，默认将第二个参数设置为输出图像所定义的灰度级别，并规定其取值被设定为64
Y = histeq(I);
**figure;**imshow(Y);
**figure;****imhist(Y);****title('直方图均衡化后的DICOM图像直方图');****Y1 = histeq(I,256);****figure;****imhist(Y1);****title('量化级为256的直方图均衡化后的DICOM图像直方图');****Y2 = histeq(I,2);****figure;****imhist(Y2);****title('量化级为2的直方图均衡化后的DICOM图像直方图');****J = imadjust(I,[0.45 0.75],[0 1]);****figure;****imshow(J);****title('灰度变换后的DICOM图像');****figure;****imhist(J);**title('灰度变换后的DICOM图像直方图');
(三) 医学影像的灰度运算与代数运算
**clc;clear;close all;****I = imread('spine.tif');****imshow(I);title('原始脊柱影像','fontsize',20);J1 = imnoise(I,'gaussian',0,0.01); % 随机噪声-高斯噪声，均值0，标准差0.01figure;imshow(J1);****J2 = imnoise(I,'gaussian',0,0.01);****J3 = imnoise(I,'gaussian',0,0.01);****J4 = imnoise(I,'gaussian',0,0.01);**avarage1 = 0.5J1 + 0.5J2;avarage2 = 0.25J1 + 0.25J2 + 0.25J3 + 0.25J4;figure;imshow(avarage1);figure;imshow(avarage2);% 100次求平均K = zeros(512,512);for i = 1:100J = imnoise(I,'gaussian',0,0.02);****J = im2double(J);K= K + J;endK = K/100;**figure;imshow(K);
(四)医学影像的几何运算
clc;clear;close all;I = imread('spine.tif');imshow(I);title('原始脊柱影像','fontsize',20);%% 图像平移% se = translate(strel(1),[150 150]); % 向下移动150像素，向右移动150像素**% Y = imdilate(I,se); % 将影像按照se进行变换****% figure;****% imshow(Y);****Y0 = imtranslate(I,[150 150]);figure;imshow(Y0);se1 = translate(strel(1),[-150 -150]); % 向上移动150像素，向左移动150像素Y1 = imdilate(I,se1); % 将影像按照se进行变换figure;imshow(Y1);%% 图像旋转B1 = imrotate(I,-40,'bilinear');****B2 = imrotate(I,-40,'bilinear','crop');****B3 = imrotate(I,90,'bilinear');****B4 = imrotate(I,-40);****figure;****subplot(2,2,1);****imshow(B1);title('顺时针旋转40度，未裁切');****subplot(2,2,2);****imshow(B2);title('顺时针旋转40度，裁切');****subplot(2,2,3);****imshow(B3);title('逆时针旋转90度，未裁切');subplot(2,2,4);imshow(B4);title('顺时针旋转40度，最近邻');%% 图像缩放C = imresize(I,0.5);****figure;imshow(C);title('缩小影像');****D = imresize(I,2);****figure;imshow(D);title('放大影像');E = imresize(I,[200,200]);figure;imshow(E);title('规定尺寸缩放影像');%% 图像镜像F = flipud(I);****G = fliplr(I);****figure;****subplot(1,3,1);****imshow(I);title('原始影像');****subplot(1,3,2);****imshow(F);title('水平镜像');****subplot(1,3,3);**imshow(G);title('垂直镜像');
(五)医学影像的傅里叶变换
清除命令窗口、清空工作空间并关闭所有图形界面；读取图片文件至变量f；创建一个新的图形窗口；显示图像到当前图形窗口；执行二维快速傅里叶变换（FFT）以获取频谱数据；计算频谱幅度：通过求取每个复数元素的模值来获得频谱图；显示计算出的频率幅值分布图；将频谱中心移动至图像中心位置；对频谱进行对数缩放处理以减小动态范围；执行逆向快速傅里叶变换并取其实部得到还原后的图像；显示还原后的图像
四、讨论与结论
对于图像中的代数运算与几何运算的基本原理感到困惑，在深入理解其相关知识时会遇到较大的困难。为了更好地理解和记忆这些内容，请尝试通过形象化的手段进行学习。
加法运算主要应用:
- 通过计算所有相同位置的像素均值来减少叠加性噪声，在同一景物连续摄取多幅带有互不相关且均值为零的高斯噪声的情况下特别有效；
- 将配准后的同一幅图作为图层覆盖至另一张图片上以改善视觉效果；
- 在高光谱成像技术中使用逐通道的线性组合扩展可见光范围，在绿色与红色通道成像结果结合下可获得接近全色的空间分布信息；
- 常用于合成测试场景并拼接边缘区域以提升整体效果
减法运算的主要应用领域包括：
1. 比较两幅图像间的差异并识别同一场景下的变化情况。例如，在运动目标检测中使用背景减法技术来识别动态物体；此外，在视频监控系统中也可通过此方法检测镜头边缘的变化；
2. 用于消除不必要的叠加元素（如缓慢变化的背景阴影或周期性的噪声）以及在每个像素处已知附加污染的情况（如电视制作中的蓝屏技术）；
3. 在图像分割中被广泛应用；
**4. 用于生成合成图像并增强视觉效果。
学生实验报告
| 学生姓名 | 景涛 | 学号 | B20040301002 | 组员： |
|---|---|---|---|---|
| 实验项目 | 医学图像增强 | |||
| □必修 □选修 | □演示性实验 □验证性实验 □操作性实验 □综合性实验 | |||
| 实验地点 | 实验日期 |
一、实验目的及要求
1**、目的**
通过Matlab编程实现方式，在实验环节中帮助学生深入理解并熟练掌握医学影像增强技术的基本原理与方法。具体包括学习线性空间变换理论及其应用、探索高斯噪声与椒盐噪声去除算法，并实践空间锐化处理技术等核心内容。
2**、**内容及要求：
1、线性空间变换
2、高斯噪声和椒盐噪声的去除
3、空间锐化
二、实验工具
| 仪器名称 | 规格/型号 | 数量 | 备注 |
|---|---|---|---|
| 计算机 | Windows10 | 1 | 有网络环境 |
| Matlab | 2016b及以上 | 1 | |
| Python | 3.6及以上 | 1 |
三、实验步骤及结果
1、线性空间变换
使用clc指令清空命令行窗口；调用clear函数删除工作区中的变量；使用close函数关闭所有未保存的图形窗口；将图像文件'breast_xray.tif'读入变量f中
f = im2double(f); %转化为浮点图像
[M,N] = size(f);figure;imshow(f);figure;imhist(f);%原图的直方图figure;
通过调用图像归一化处理后的直方图计算函数imhist对输入图像f进行处理，并将其分为64个等分区间来生成相应的频率分布曲线；然后使用茎状图形展示各区间对应的频率值；其中x存储了各个灰度区间对应的索引值；而H则表示每个区间内的像素数量；通过这种方式可以直观地观察到输入图像f在不同灰度范围内的分布情况
%增加对比度a1= 2;b1 =-55;f1= a1.*f+ b1/255;figure;imshow(f1);title('a=2 b=-55增加对比度');
%减小对比度a2=0.5;b2 =-55;f2= a2.*f+ b2/255;figure;imshow(f2);title('a=0.5 b=-55 减小对比度');
%线性平移增加亮度a3= 1;b3= 55;f3= a3.*f+ b3/255;figure;imshow(f3);title('a=1 b=55线性平移增加亮度');
%反相显示a4 =-1;b4 = 255;f4 = a4.*f + b4/255;figure;imshow(f4);title('a=-1 b=255反相显示');
figure;subplot(1,5,1);
adjust灰度变换：
clc;clear;f = imread('breast_xray.tif');figure;imshow(f);
g1 = imadjust(f,[0 1],[1 0]); % 反转图像，负片图像figure;imshow(g1);
g2 = imadjust(f,[0.5 0.75],[0 1]); % 将0.5-0.75之间的灰度级扩展到0-1，突出亮度带figure;imshow(g2);
g3 = imadjust(f,[],[],2); % 压缩灰度级的低端并扩展灰度级的高端figure;imshow(g3);
g4 = im2bw(f); % 二值化，默认阈值取0.5figure;imshow(g4);
2、高斯噪声和椒盐噪声的去除
清除命令窗口；清空工作区内容；关闭所有图形界面；读取 knee1.dcm 格式的 DICOM 图像到变量 I0 中；创建并显示 I0 的图像；读取 knee.jpg 格式的 JPEG 图像到变量 I 中；创建并显示 I 的图像；
% 用3*3模板相关滤波 w = [1 1 1;1 1 1;1 1 1]/9;g = imfilter(I,w);figure;imshow(g);
% 添加高斯噪声I1 = imnoise(I,'Gaussian',0,0.1);figure;imshow(I1);
% 添加椒盐噪声I2 = imnoise(I,'salt & pepper',0.1);figure;imshow(I2);
对原始图像应用不同尺寸的平均模板进行平滑滤波处理：w1 = fspecial('average', 3)；k1 = imfilter(I, w1)；随后使用5×5尺寸的平均模板生成k2；接着应用7×7尺寸的平均模板得到k3；最后在同一个窗口中展示原始图像及其三种不同尺寸的平滑后结果图
对椒盐噪声图像进行处理并应用不同尺寸的中值滤波器。通过使用medfilt2函数实现对输入图像I的不同尺寸中值滤波。figure; subplot(…)；显示原始图像及其经过不同尺寸中值滤波后的效果。具体而言：使用[3×3]大小模板生成f1；[5×5]大小模板生成f₂；[7×7]大小模板生成f₃；以及[9×9]大小模板生成f₄等结果图。
3、空间锐化
clc;clear;close all;I = imread('knee.jpg');figure;imshow(I);
以下是按照要求对原文的内容进行同义改写的文本：
% sobel算子进行锐化处理h1 = fspecial('sobel');s1 = imfilter(I,h1);figure;imshow(s1);
% 使用Prewitt算子执行图像锐化
h2 = fspecial('prewitt');
对输入图像I应用预定义的滤镜h1以获得锐化后的图像s2
figure;
imshow(s2);
% prewitt算子进行边缘提取Y = edge(I,'prewitt',0.05);figure;imshow(Y);
四、讨论与结论
近年来，在科技进步的推动下，图像处理技术不断进步，在帮助医生完成诊断与治疗的过程中发挥了重要作用；特别是在临床医学领域的发展中产生了深远的影响。然而，在摄取与传输过程中容易受到多种因素的影响，并最终导致后续处理与分析的效果受到较大程度的影响；尤其是在一些定量与定性分析方面表现尤为突出。
除了在日常生活与科学研究中得到广泛应用之外，在医学领域也发挥着不可或缺的作用。作为图像处理的核心技术之一，同时也是解决许多实际问题的关键技术。该技术涵盖了一维与二维信号能量分布不均程度的量化评估以及增强等技术手段等。
医学图像具有高分辨率、较低的噪声干扰以及在病灶部位有更加清晰显示的特点；这些特点使得医学图像处理面临更为复杂的挑战，并且存在一系列亟待解决的问题，例如评估医学图像的质量等。
学生实验报告
| 学生姓名 | 景涛 | 学号 | B20040301002 | 组员： |
|---|---|---|---|---|
| 实验项目 | 医学图像分割 | |||
| □必修 □选修 | □演示性实验 □验证性实验 □操作性实验 □综合性实验 | |||
| 实验地点 | E405、E503 | 实验日期 | 2023.2-2023.6 |
一、实验目的及要求
目的：
通过Matlab编程技术实现医学影像分割的学习与实践教学模式中，引导学生深入理解并熟练掌握各种基础分割方法。
内容及要求：
1、边缘检测算法
2、阈值分割算法
3、区域增长算法
4、聚类分割算法
二、实验工具
| 仪器名称 | 规格/型号 | 数量 | 备注 |
|---|---|---|---|
| 计算机 | Windows10 | 1 | 有网络环境 |
| Matlab | 2016b及以上 | 1 | |
| Python | 3.6及以上 | 1 |
三、实验步骤及结果
1、边缘检测算法
clc; % 使用"delete"函数代替重复的clc指令
clear; % 将工作区中的所有变量清空
close all; % 通过一次操作关闭所有未保存的工作区窗口
nii = load_nii('brats20_training_001_flair.nii'); // 该代码读取并加载了特定类型的MRI数据文件
%img = nii.img; %[n1,n2,n3] = size(img); % 但原本意图似乎是想获取图像尺寸信息
%i = imrotate(img(:,:,80),90); % 对第80个切片旋转90度
i = imread('headct.tif'); % 读取头部ct扫描数据
imshow(i,[]); % 显示图像，并自动缩放亮度范围
注：
由于第一条代码清除命令行窗口，则之前在命令行中定义的所有变量都将被移除。第二条代码清除工作区，则当前的工作区中没有任何变量留存。接下来的部分是一个处理医学图像数据的程序。但可能由于某种原因未能实现。
而最后两行代码捕获了名为'headct.tif'的CT图像文件，并将其呈现出来。imshow函数通过[]参数自动调节亮度范围，在呈现图像时确保整个颜色映射范围被覆盖（即黑白）。
% 读取图像
I = imread('image.jpg');
% 使用 Roberts 算子检测边缘
b1 = edge(I,'roberts');
% 使用 Prewitt 算子检测边缘
b2 = edge(I,'prewitt');
% 使用 Sobel 算子检测边缘
b3 = edge(I,'sobel');
% 使用 Laplacian of Gaussian（LoG）算法检测边缘
b4 = edge(I,'log');
% 使用 Canny 算子检测边缘
b5 = edge(I,'canny');
生成一个新的图像窗口，并将其中包含原始图像以及所有边缘检测结果的内容以2×3网格的形式进行展示
figure;
subplot(2,3,1); imshow(I,[]); title('原始图像');
subplot(2,3,2); imshow(b1,[]); title('Roberts');
subplot(2,3,3); imshow(b2,[]); title('Prewitt');
subplot(2,3,4); imshow(b3,[]); title('Sobel');
subplot(2,3,5); imshow(b4,[]); title('LoG');
subplot(2,3,6); imshow(b5,[]); title('Canny');
注：
该代码执行了对输入图像进行边缘检测的任务，并将不同算法产生的结果以2×3网格的形式展示在同一幅图中。具体而言，在该代码中分别采用了roberts、prewitt、sobel、log和canny等五种经典的边缘检测算法进行运算。其中edge( )函数是matlab中专门用于完成边缘检测操作的重要函数，在其调用过程中第一个参数指定待处理的对象图像信息，第二个参数则决定了采用哪种具体的边缘检测算法实现目标。随后利用imshow( )函数生成相应的图形窗口以直观展示计算结果。最终呈现给用户的是一幅包含原始图像及其五种不同算法处理后效果的整体图集。
% 设置参数提高边缘检测性能
[gs,ts] = edge(i,'sobel'); % 默认参数的 sobel 边缘检测
[gl,tl] = edge(i,'log'); % 默认参数的 log 边缘检测
[gc,tc] = edge(i,'canny'); % 默认参数的 canny 边缘检测
gsb = edge(i,'sobel',0.1); % 更改 sobel 边缘检测门限
glb = edge(i,'log',0.0065,2.25); % 更改 log 边缘检测参数
gcb = edge(i,'canny',[0.06 0.1],1.5); % 更改 canny 边缘检测参数
figure; % 创建新窗口展示结果
subplot(2,2,1);imshow(i,[]);title('original');
subplot(2,2,2);imshow(gsb,[]);title('sobel');
subplot(2,2,3);imshow(glb,[]);title('log');
subplot(2,2,4);imshow(gcb,[]);title('canny');
注：
该代码被用来完成图像边缘检测任务。它采用了不同设置来执行Sobel、Laplacian of Gaussian（LOG）以及Canny等三种边缘检测方法，并在子图中展示了这些方法的结果。通过调整算法中的相关参数可以调控边缘检测的敏感程度。
2、阈值分割算法
将图像I从'bacteria.bmp'读取为矩阵；然后将图像J从'blood1.bmp'读取为矩阵；接着在图形窗口中创建一个2x2的子图布局，并在第一个子图中显示图像J；最后在第二个子图中显示图像J的直方图；
% J is set to values above 50;
% global threshold set at 50, the minimum value is approximately 50
J is assigned to values greater than 68;
% global threshold set at 68, the minimum value is approximately 68
J is assigned to values greater than 79;
% global threshold set at79, the minimum value is approximately79
K is assigned to values greater than89;
% global threshold set at89, the minimum value is approximately89
figure；subplot(131), imshow(i);
subplot(132),imshow(j); %全局阈值为120时的二值化图像
subplot(133),imshow(k); %全局阈值为130时的二值化图像
注：
在代码中, 第一行被注释掉, 未被利用. 接下来依次选取不同的全局阈值进行二值化处理, 生成了两个二值化图像j和k. 最后利用subplot函数在同一窗口展示三张图片: 原始图像以及基于120和130两个全局阈值得到的二值化图.
3、区域增长算法
disp(''); %清除命令行窗口
clear; %清除工作区
close all; %关闭前面所有窗口
I = imread('headCT.tif'); imshow(I); %导入头CT图像并显示
s = qtdecomp(I, 0.2); %进行四叉树分解操作
s2 = full(s); %生成完整的矩阵表示
注：
这段代码的核心功能是针对文件名 "headct.tif" 的二维灰度图像进行四叉树分割，并输出其分解结果。具体实现过程如下：
clc 用于清除 matlab 命令行窗口中已有的文本。
clear 用于清空工作区中定义的所有变量，以释放内存空间。
close all 用于关闭所有 matlab 图形界面窗口，以便于后续展示新的图像。
调用 imread 函数加载名为 'headct.tif' 的二维灰度阵列（tif 格式），并将该图像赋值给变量 i。
imshow(i) 将变量 i 中存储的二维灰度图像以弹出窗口形式展示出来。
qtdecomp(i, 0.2) 对变量i中的原始图像执行四叉树分解操作，并基于阈值0.2生成一个稀疏矩阵对象，并将其赋值给变量s
full(s) 实现了通过图像分解获得的稀疏矩阵转换为普通矩阵的过程，并将结果存储于变量s2中
%读取CT头部扫描图像
I = imread('headCT.tif');
%进行四叉树分解，使用0.26作为阈值来控制分裂
S = qtdecomp(I, .26);
%初始化每个块的值为0，大小与四叉树分解后相同
blocks = repmat(uint8(0), size(S));
%按从大到小的顺序（方便后面覆盖），遍历每个分裂的维度
for dim = [512 256 128 64 32 16 8 4 2 1]
%查找哪些块需要进一步分裂
numblocks = length(find(S == dim));
if (numblocks > 0)
%新建一个矩阵，宽高与dim相同，每个像素初始值为1
values = repmat(uint8(1), [dim dim numblocks]);
%将values中除第一行和第一列之外的所有元素赋值为0
values(2:dim, 2:dim, :) = 0;
%根据values对blocks进行分块处理
blocks = qtsetblk(blocks, S, dim, values);
end
end
%对最后一行和最后一列的块赋值为1，这是为了避免边缘问题
blocks(end, 1:end) = 1;
blocks(1:end, end) = 1;
%在两个子图中分别显示原始图像和四叉树分解后的图像
subplot(121);
imshow(I);
title('原始图像');
subplot(122);
imshow(blocks,[]);
title('块状四叉树分解图像');
注：
这段代码主要是执行了块状四叉树分解，并呈现了分解后的分块图像。具体步骤如下：
首先对读入的灰度图像i进行四叉树分解，使用qtdecomp函数，阈值为0.26。
随后设置blocks变量为一个具有与s相同维度的全零矩阵，用于表示最终的块状四叉树分解结果。
针对每个块大小dim（从最大尺寸512开始逐步减少），计算当前尺寸为dim的区块数量numblocks。若numblocks大于0，则将values变量设为1以指示该区块需进一步细分；若无则忽略。
最后使用qtsetblk将细分后的块写入blocks中。
将矩阵的末行和末列赋值为1，则表明所有子块已经完成计算，并且这将使得最终输出更加整洁。
最后将原始图像和块状四叉树分解图像以subplot的形式展示出来。
整体功能主要依赖于四叉树分解方法来实现图像压缩，并且该段代码还支持以块状形式组织信息以增强数据的局部性。
4、聚类分割算法
clear,clc;close all;I = imread('head80.png');I = rgb2gray(I);imagesc(I);注：
这段代码调用matlab对该彩色图像文件'head80.png'进行图像数据获取操作，并将其转码为灰度图像随后以图形界面的形式展示于图形界面窗口中
imread()：读取图像文件
rgb2gray()：将彩色图像转换为灰度图像
imagesc()：在图形窗口中显示图像
clear: 清除工作区中的变量
clc: 清空命令行窗口
close all: 关闭所有打开的图形窗口
% figure; % 注释掉的话，不会显示最初的图像
imshow(I); % 显示图像 I
[m,n] = size(I); % 获取图像 I 的尺寸
np = m * n; % np 表示图像中所有像素点的个数
data = reshape(double(I),np,1); % 将图像转化为一个列向量 data（每行一个像素）
k = 6; % 指定聚类簇数目 k
res = kmeans(data,k,'dist','sqEuclidean'); % 进行 K-Means 聚类，并输出聚类结果 res
result = reshape(res,m,n); % 实现了将聚类结果转译成图像形式，并将其恢复为与原始图像相同尺寸的矩阵
figure; % 创建新的窗口用于显示聚类后的图像
imagesc(result); % 显示聚类后的图像（类别之间用不同的颜色表示）
注：
该代码实现了K-means算法对图像数据进行聚类分割。kmeans()函数中的第三个参数'dist'用于表示所采用的距离度量方法来计算数据点之间的距离。其中'sqEuclidean'参数指定了使用欧氏平方距离作为度量标准。
四、讨论与结论
医学图像分割涉及从背景中分离出特定组织或病变区域的过程。它在医学影像学中发挥着重要作用，并且能够为医生提供更为精准和迅速的诊断支持。然而，在这一领域也面临诸多挑战与难点。
以下是医学图像分割过程中的重点：
数据预处理工作：医学图像在获取过程中往往伴随着噪声污染以及边缘模糊等问题。为了改善这些问题的影响，通常需要开展数据预处理工作，如去除噪声干扰、消除模糊边缘等常规步骤。
在图像处理领域中进行分类或分割任务时, 人们需要明确识别所需的关键属性, 并采用合适的算法来实现有效的数据处理. 特征提取过程不仅可以依赖传统的特定模式识别方法, 还可以通过结合复杂的人工智能模型来实现更加精确的结果.
分割算法：目前有很多分割算法被应用于医学图像领域中。这些方法主要包括基于传统技术（如阈值化、边缘检测以及区域增长）以及基于深度学习的模型（如卷积神经网络与UNet架构）。为此建议根据具体应用场景选择最适合的分割技术。
评估方法：为了准确地评价分割算法的表现效果，请采用科学合理的评价手段。其中常用的方法包括：Dice系数、SSIM指数以及灵敏度与特异性指标等。
难点如下：
图像中的噪声与低对比度问题是医学影像采集过程中的常见挑战，在这一背景下，在实际应用中往往会导致目标物体难以清晰明确地区分出来
在医学影像领域中，目标的形态与构造往往具有高度复杂性。例如，在CT或MRI等医学成像技术中可能会遇到肺部肿瘤、血管网状结构以及脑部组织等复杂的形态与构造。这些不同的构造之间彼此间的干扰因素包括穿透性问题以及分辨率限制等因素的存在都会显著影响分割效果。
多模态影像：临床医学影像通常涉及多种图像捕捉技术,如x射线断层扫描（CT）、磁共振成像（MRI）等,各自具有独特的优缺点,在实际应用中表现各异且各有不足之处.在捕捉图像时会呈现出不同视角,这无疑会对分割处理带来诸多困难.
在医疗领域中：由于患者隐私保护与安全监管等方面的限制，在获取、存储与应用医学影像数据方面面临着巨大挑战
人工标注难度：针对大规模数据集开展分割任务时必须投入大量的人力资源用于标注工作。该过程既耗时且复杂,同时会耗费精力并造成成本高昂,最终可能导致主观判断差异大以及结果的一致性难以保障
医学图像分割可以使用不同编程语言实现，以下是一些常见的代码语句：
Python
Python 是广泛应用于机器学习与深度学习的编程语言。在 Python 中提供了丰富的工具库用于执行图像分割任务，并包括如TensorFlow、Keras和PyTorch等。
加载图像文件
python
import cv2
image = cv2.imread('image.jpg')
预处理图像
python
import numpy as np
from skimage import filters
image_gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
sobel = filters.sobel(image_gray)
normalized_image = np.zeros((image.shape[0],image.shape[1]))
normalized_image = cv2.normalize(sobel, normalized_image, alpha=0,
beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
像素分类 / 分割
python
from sklearn.cluster import KMeans
kmeans被初始化为KMeans类的一个实例；然后对经过归一化处理的图像样本进行拟合
result = kmeans.labels_.reshape(image_gray.shape[0], image_gray.shape[1])
MATLAB
MATLAB已经成为医学图像处理与分割领域的重要工具之一，并拥有丰富全面的图像处理工具箱以及模块集合。基于上述原因，在进行医学图像分割时选择 MATLAB 是非常高效的。
加载图像文件
matlab
image = imread('image.jpg');
预处理图像
matlab
gray_image = rgb2gray(image);
sobel = edge(gray_image, 'Sobel');
normalized_image = imadjust(sobel, stretchlim(sobel),[0 1]);
像素分类 / 分割
matlab
num_of_classes = 2;
result = kmeans(normalized_image(:), num_of_classes);
result = reshape(result, size(normalized_image));
学生实验报告
| 学生姓名 | 景涛 | 学号 | B20040301002 | 组员： |
|---|---|---|---|---|
| 实验项目 | 医学图像配准与融合 | |||
| □必修 □选修 | □演示性实验 □验证性实验 □操作性实验 □综合性实验 | |||
| 实验地点 | E405、E503 | 实验日期 | 2023.2-2023.6 |
一、实验目的及要求
目的：
通过Matlab编程实现医学影像配准与融合技术的学习与实践，在实验过程中帮助学生全面了解这一领域的基本概念并熟练掌握以下几种方法：包括基于刚性变换的空间配准、基于平移变换的快速配准算法以及多种像素级融合法如极大极小法和加权平均法等基础理论与应用技能。此外还能深入学习区域能量最大化准则的小波变换域融合方法及其在医学图像处理中的实际应用。
内容及要求：
1、医学图像配准
2、图像像素极大极小融合法
3、图像像素加权法
4、区域能量最大法
5、小波变换下的融合
二、实验工具
| 仪器名称 | 规格/型号 | 数量 | 备注 |
|---|---|---|---|
| 计算机 | Windows10 | 1 | 有网络环境 |
| Matlab | 2016b及以上 | 1 | |
| Python | 3.6及以上 | 1 |
三、实验步骤及结果
1、医学图像配准
% 清除命令行窗口和工作区
clc;
clear;
close all;
% 读入图像 f2 和 f1，并进行旋转
f2 = imread('MRI.jpg');
f1 = imrotate(f2,30,'loose');
% 调用 cpselect 函数，手动选择配准点
cpselect(f1, f2);
% 在窗口中显示原始图像 f1 和 f2
figure;
subplot(2,2,1), imshow(f1);title('配准前的图1');
subplot(2,2,2), imshow(f2);title('配准前的图2');
% 利用非反射相似变换获取相似变换对象 tform
赋值给变量tform的几何变换函数接受移动点集和固定点集作为参数，并采用非反射相似性
% 进行仿射变换，获取处理后的图像 Iout
Iout = imwarp(f1,tform);
% 在窗口中显示处理后的图像
subplot(2,2,3),imshow(Iout);title('配准后的图1');
figure;imshow(Iout);title('配准后的图1');
注：
这段代码的主要功能是将两幅图像f1和f2进行配准操作。其中f1经过旋转后出现了一定的偏移现象。整个过程首先使用imread函数读取原始未变形的测试图片f2以及变形后的目标图片f1；接着，在图形界面中手动选择了关键控制点，并基于这些选中的关键控制点计算出相似变换矩阵tform；最后利用imwarp函数执行仿射变换得到处理后的变形校正图片Iout，并将其显示在图形界面中完成操作流程。
2、图像像素极大极小融合法
command window清空;
% command window清空;
工作区变量清除;
% work space variables cleared;
图形界面窗口关闭;
% graphical user interface windows closed;
读取mri图像到变量f1中;
% Assign the mri image to variable f1.
f2 = imread('ct.jpg'); %读取配准后的ct图像
[m1, n1] = size(f1); %获取mri图像的大小
[m2, n2] = size(f2); %获取配准后的ct图像的大小
if (m1=m2)||(n1=n2)
出现错误('图片尺寸不匹配'); %若输入图片尺寸不符则触发错误并终止程序
end
%像素灰度值极大法融合
for i = 1:m1
for j = 1:n1
if (abs(f1(i,j))>abs(f2(i,j)))
f3(i,j) = f1(i,j);
elseif (abs(f1(i,j))<abs(f2(i,j)))
f3(i,j) = f2(i,j);
end
end
end
%像素灰度值极小法融合
for i = 1:m1
for j = 1:n1
if (abs(f1(i,j))>abs(f2(i,j)))
f4(i,j) = f2(i,j);
elseif (abs(f1(i,j))<abs(f2(i,j)))
f4(i,j) = f1(i,j);
end
end
end
%显示图像
figure;
subplot(2,2,1);
imshow(f1),title('原始mri图像');
subplot(2,2,2);
imshow(f2),title('配准后的ct图像');
subplot(2,2,3);
imshow(f3),title('像素灰度值极大法融合后的图像');
subplot(2,2,4);
imshow(f4),title('像素灰度值极小法融合后的图像');
注：
这段代码执行了完成对两张医学影像图像的融合操作，并生成并展示了融合前后的四张图像。详细说明了具体的步骤。
读取两张图像（mri和ct），并获取它们的尺寸。
判断这两张图像的尺寸是否相同，如果不同就抛出一个错误提示信息。
扫描两张图片的所有像素点。对每个像素进行比较分析这两个像素值的大小关系。通过比较结果设置新图像该处对应的灰度级。
采用极大或极小的像素灰度值进行处理后，生成两张融合图像
调用matlab中的subplot、imshow和title函数以显示原始mri图像、对齐后的ct图像以及通过两种融合方法获得的结果图，并对比分析其效果
总体而言，这段代码不仅能够良好地融合两个医学图像，并且能提供直观的图像展示供医生或研究者进行进一步的研究和分析。
3、图像像素加权法
清空命令行窗口界面；
%clc;用于清空命令行窗口界面；
%clc;的作用是清除当前终端会话中的所有输出内容；
% 用于清空工作区内存内容；
clear;
%clear;的功能是删除工作区内所有的变量和定义；
%关闭当前及其之前的所有图形窗口；
close all;
%maintain the current figure window and close all previous ones;
% 读取两张图片；
f1 = imread('mri.jpg');
f2 = imread('CT.jpg');
% 将图像转换为double类型表示
f1 = im2double(f1);
f2 = im2double(f2);
% 对两幅图像做二维FFT变换，得到频域表示
y1 = fft2(f1);
y2 = fft2(f2);
% 通过不同的加权系数对频域表示进行加权平均
y3 = 0.3y1 + 0.7y2;
y4 = 0.5y1 + 0.5y2;
% 对两幅图像的频域表示做反变换，得到多通道融合图像
f3 = ifft2(y3);
f4 = ifft2(y4);
% 将多通道融合图像转换为无符号8位整数类型，方便显示
f3 = im2uint8(f3);
f4 = im2uint8(f4);
在窗口中分别呈现原始的MRI图像、经过配准的CT图像以及两幅经过加权融合的图像
subplot(2,2,1); imshow(f1), title('原始MRI图像');
subplot(2,2,2); imshow(f2), title('配准后的CT图像');
subplot(2,2,3); imshow(f3), title('加权系数0.3，0.7融合后的图像');
subplot(2,2,4); imshow(f4), title('加权系数0.5，0.5融合后的图像');
注：
这段代码用于图像融合和显示。这段代码中包括以下几个步骤：
调用imread()函数获取两张图像分别为mri.jpg和CT.jpg，并将其转换为双精度数据类型后存储于变量f1和f2中。这两张图片分别被保存在.mat文件中的变量f1和f2中。
完成对两幅图像的二维快速傅里叶（FFT）变换后, 分别计算出它们的频谱描述, 并将其存储到变量 y1 和 y2 中
采用加权平均的方法计算两个频域表示以获得变量y3和y4这些变量分别对应于第一种方案中权重分别为0.3与0.7以及第二种方案中权重均为0.5的情况
执行反变换操作于这两个经过加权处理的频域表示，从而获得变量f3和f4。这些变量分别对应两种不同的加权方案下的图像输出。
为了将两个结果图像的数据类型设置为无符号8位整数（uint8），以便正确显示
展示这些图像并利用 subplot()、imshow() 和 title() 函数分别呈现原始 MRI 图像、经配准的 CT 图像以及加权处理后的两幅图片于四个子图框中。
这段代码可作为展示图像融合技术效果的示例使用。你可以将此段代码复制到MATLAB命令窗口中执行它会自动读取图片并进行处理这将使程序自动读取图片并显示处理结果。另一种方法是将此段代码保存为.m文件并在MATLAB环境中通过执行该文件来实现功能
4、区域能量最大法
clc; %执行"clc"指令以清空命令行界面
clear; %执行"clear"指令以释放工作区空间
close all; %关闭当前及其所有相关窗口
首先读取mri.jpg文件并将其转换为灰度图像赋值给变量f1。
f2 = rgb2gray(imread('CT.jpg')); % 读取 CT.jpg 并转换为灰度图像
f3 = fenergy(f1,f2); % 对两张图像进行能量融合
subplot(131);imshow(uint8(f1));title('MRI'); % 在第一个子图中显示原始 mri 图像
subplot(132);imshow(uint8(f2));title('CT'); % 在第二个子图中显示原始 CT 图像
axes位置设置为1行3列第三块；使用imshow函数显示uint8类型的f3；设置标题为'区域能量最大值'；注：在第1行第3列的子图中展示融合后的图像。
function new = fenergy(a,b)
% 实现区域能量融合的函数
a = double(a); % 将 a 转换为双精度型数据
b = double(b); % 将 b 转换为双精度型数据
[m,n] = size(a); % 获取图像尺寸
temp_a = nlfilter(a,[3 3],@nengliang); % 使用 nlfilter 函数，对 a 进行邻域滑动操作
temp_b = nlfilter(b,[3 3],@nengliang); % 使用 nlfilter 函数，对 b 进行邻域滑动操作
for i=1:m % 遍历所有像素点
for j=1:n
if temp_a(i,j)>=temp_b(i,j) % 根据能量大小进行融合
new(i,j)=a(i,j); % 如果 a 的强度不低于 b，则会将 a 对应位置的像素值导入新图像中
else
new(i,j)=b(i,j); % 如果 b 的能量大于 a，则将 b 对应位置的像素值导入新图像中
end
end
end
% new=(new);
end
function c = nengliang(x)
% 计算给定区域的能量值的函数
A = [1,2,1]; % 定义加权矩阵
C = A*A'.x.(1/16); % 计算加权矩阵与当前滑动区域对应像素值的乘积
c = sum(sum(C)); % 对所有乘积进行求和，得到当前区域能量值
end
注：
在主程序运行中，在开始阶段调用 imread 函数加载两张彩色图片，并随后分别运用 rgb2gray 转换它们为灰度图象。随后调用 fenergy 函数对这两幅图象进行能量融合处理以获得融合后结果图象。最终运用 subplot 和 imshow 函数在同一图形界面中展示这三张图片，并附加图形标题。
在 fenergy 函数中，首先将输入的两张图像依次转换为双精度数据类型，并通过 nlfilter 函数对每张图像执行 3×3 邻域滑动滤波操作，在每个滤波窗口内计算对应的能量值。随后遍历所有像素点，并对每一点处的两张图像的能量值进行比较：当第一幅图像在该位置的能量值不小于第二幅时，则将该像素点处的第一幅图像素值赋入新图像；否则采用第二幅图对应位置的像素值。
Within the energy function, a weighting matrix A = [1,2,1] is defined. This weighting matrix is then multiplied by the corresponding pixel values in the current sliding window. The result of this multiplication is used to compute the current region's energy value. Finally, the sum of all products is calculated to obtain the current region's energy value. The function then returns this computed value.
5、小波变换下的融合
clc; %clc;清除命令行窗口 clear; %clear;清除工作区 close all; %关闭前面所有窗口 f1 = rgb2gray(imread('mri.jpg'));f2 = rgb2gray(imread('CT.jpg'));zt = 2; % 小波分解层数wtype = 'db1'; % 小波类型[c0,s0] = wavedec2(f1,zt,wtype); % 多尺度二维小波分解[c1,s1] = wavedec2(f2,zt,wtype);% 小波简单加权法c = (c0 + c1) * 0.5;% 高频部分系数选择绝对值极大法，低频部分系数采用二者求平均的方法KK = size(c1);Coef_Fusion1 = zeros(1,KK(2));% 低频系数的处理Coef_Fusion1(1 : s1(1,1)) = (c0(1 : s1(1,1))+c1(1 : s1(1,1)))/2;% 高频系数的处理MM1 = c0(s1(1,1) + 1 : KK(2));MM2 = c1(s1(1,1) + 1 : KK(2));mm = (abs(MM1)) > (abs(MM2));Y = (mm .* MM1) + ((~mm) .* MM2);Coef_Fusion1(s1(1,1)+1 : KK(2)) = Y;% 小波重构Y1 = waverec2(c,s0,wtype);Y2 = waverec2(Coef_Fusion1,s0,wtype);subplot(2,2,1);imshow(f1),title('原始MRI图像');subplot(2,2,2);imshow(f2),title('配准后的CT图像');subplot(2,2,3);imshow(Y1,[]),title('小波系数加权法融合后的图像');subplot(2,2,4);imshow(Y2,[]),title('低频加权融合，高频取绝对值极大法融合');
注：
这段代码实现了基于小波分解的图像融合操作。具体实现过程如下：
首先使用 imread 函数读取两张彩色图像，并将它们转换为灰度图像。
设定小波分解的层数与小波类型，并利用 wavedec2 函数对两张图像进行多级二维小波分解,获得相应的低频分量与高频分量
对高频图像的系数执行特定权重分配以达到预期效果；同时将低频区域的数据通过算术均值计算整合；从而生成一个新的组合向量
对每个位于高频系数向量中的元素进行遍历，在计算各个元素绝对值的基础上确定最大值，并将其结果作为新的高频系数向量 Y。
通过融合低频特征与最新的高频特征向量，并将其整合为一个整体序列,从而获得最终融合后的系数矩阵 Coef\_Fusion1
使用 waverec2 函数对两种系数进行小波重构，得到两张融合后的图像。
通过 subplot 和 matplotlib 的 imshow 函数分别呈现原始图像及其融合结果于同一绘图窗口，并在图形顶端添加适当的标题以明确展示内容
就目前而言，这段代码采用了基于小波分解的图像融合方法，并分别提供了两种不同的实现方案：主要采用了基于简单平均法的小波系数加权融合方案；另一种方案中低频分量采用了简单平均法、高频分量则采用了绝对值极大法来进行系数融合。同时，在代码实现过程中使用了abs函数来计算向量各元素的绝对值。
四、讨论与结论
医学影像配准与融合是指一种旨在精确对齐并整合多幅医学影像以生成完整合成影像的过程。该过程有助于医生更加高效地完成病情分析及诊断工作。其中配准阶段主要通过空间变换与形态重塑技术实现不同来源影像间的相互对应；而融合阶段则通过整合这些经过精确对齐的影像信息以获取更为全面的临床观察数据
在医学图像配准领域中，常用的配准手段主要包括特征点匹配技术、相似度量方法以及形变模型的应用等多类方法。其中一种常用的方法是通过提取关键点及其特征描述子来实现匹配过程；为了提高匹配质量，在实际应用中通常会结合RANSAC算法来剔除噪声匹配项；此外，在某些情况下还可以借助几何变换手段达到空间对齐的目的。例如，在仿射变换下图像的平行性及比例关系得以保持；而非线性变换则能够更好地适应复杂形状的变化需求；这些不同的变换策略可以根据具体的配准需求进行选择与应用
在医学图像融合领域中，广泛采用的方法主要包括基于像素级别的加权融合、基于小波分解的多分辨率处理以及基于深度学习的特征提取等技术手段。例如，在像素级别上进行操作时，默认情况下会对不同分辨率或质量的医学图像进行权重分配以达到平衡效果。此外，在深度学习框架下通过卷积神经网络（CNN）对多源医学影像进行统一建模与处理，在提升特征提取精度的同时实现了更为精确的图像重建效果。
医学图像配准的重点包括以下几个方面：
图像几何变换：在医学成像领域中常见的几何变化包括旋转、平移以及缩放操作。配准过程的核心在于对这些几何变化进行精确估算和有效还原。这一过程确保了对齐后的影像能够在三维空间中建立精确的一一对应关系。
开发精确的相似性度量模型：图像配准评估依靠相似性度量指标，包括互信息、相似性指数等。这些指标能够有效指导计算图像间的差异程度，并为配准过程提供可靠的评价依据。通过这些方法的应用，在计算图像之间的匹配程度时可获得关键的数据支持，在提升配准结果的信任度方面具有重要意义。
多模态医学影像配准：Within the context of multi-modal medical image registration, various imaging modalities produce images with distinct physical characteristics and signal distributions. Thus, the registration of medical images requires an understanding of the unique features and differences across various modalities.
医学图像融合的重点包括以下几个方面：
医学影像图象：在医疗领域中经常遇到噪声及其他干扰因素的影响,这些因素会影响不同医疗影像图象的融合效果.因此,在进行图象融合前必须对图象进行降噪以及提升处理,以便显著提高图片的质量与准确性.
在医学图像融合中进行模型选择时最关键的是要如何选取最适合的模型？在这一过程中使用的各种模型都必须要具备处理不同类型数据以及丰富图像特征的能力，并且这些模型不仅要能精确重构出清晰的图像还能为准确描述病灶位置提供可靠依据。
丰富多样的医学影像数据在现代医学影像融合技术中扮演着重要角色。无论是传统的CT与MRI所呈现的结构特征还是功能成像技术所揭示的功能信息，在这一过程中都发挥着不可或缺的作用。此外，在这一过程中还需要考虑放射性同位素显影技术以及显微镜下的组织样本分析等多种因素的影响。因此，在融合之前需要充分了解各类医学影像数据的特性及其应用价值
以下是基于Python的医学图像配准和融合的示例代码：
python
import SimpleITK as sitk
读取两幅医学图像
fixed_image = sitk.ReadImage('fixed_image.dcm')
moving_image = sitk.ReadImage('moving_image.dcm')
定义配准方法、转换及度量
registration_method = sitk.ImageRegistrationMethod()
initial_transform 被赋值为 sitk 中的 CenteredTransformInitializer 类函数的结果（基于 fixed_image 和 moving_image 的几何参数配置）
registration_method.SetInitialTransform(initial_transform)
registration_method.ConfigureMetricForMattesMutualInformation(numberOfHistogramBins=50)
registration_method->Assign a metric sampling strategy to the registration_method.
registration_method.SetMetricSamplingPercentage(0.01)
registration_method.SetInterpolator(sitk.sitkLinear)
optimizer is configured as gradient descent with a learning rate of 1.0, iterating for 100 steps, and estimating the learning rate per iteration.
registration_method.SetOptimizerScalesFromPhysicalShift()
registration_method.SetOptimizerStopConditionThreshold(0.001)
transform = registration_method.Execute(fixed_image, moving_image)
将移动图像进行配准
resampled_image := sitk.Resample(moving_image, fixed_image, transform, sitk.sitkLinear, 0.0, moving_image.GetPixelFormatIndex())
将两幅图像进行融合
alpha = 1.0
beta = 0.5
固定图像经过浮点转换后与重新采样的图像按权重相加生成融合图像
该代码段用于从输入图像生成一个可交互的可视化界面并对其结果进行可视化显示
该段代码的主要功能是：
首先将固定图像转换为无符号8位整数数据类型后与经重采样处理后的图像进行标签叠加显示
然后根据重采样后的图像生成对应的轮廓标签并将其与前面得到的复合标签图层叠加在一起以完成最终的可视化效果
这种方法能够有效地结合形态学操作和可视化展示过程
这段示例代码采用了SimpleITK库对医学图像进行了配准与融合操作。具体流程如下：首先获取了两幅医学图像数据，并设定好配准方法、空间变换模型以及评估指标等参数设置。随后通过调用 registration_method.Execute 函数完成了配准过程，并将移动图像中的像素值依据resampled_image进行加权平均计算后与固定图像结合生成最终融合图像fused_image。最后利用 sitk.Show 函数展示了配准效果及其边界信息图示。
学生实验报告
| 学生姓名 | 景涛 | 学号 | B20040301002 | 组员： |
|---|---|---|---|---|
| 实验项目 | 医学图像复原 | |||
| □必修 □选修 | □演示性实验 □验证性实验 □操作性实验 □综合性实验 | |||
| 实验地点 | E405、E503 | 实验日期 | 2023.2-2023.6 |
一、实验目的及要求
目的：
通过Matlab编程实现医学影像复原技术的教学实践，在实验过程中帮助学生深入理解并熟练掌握噪声对图像的影响机制以及多种图像复原方法的基本理论与应用。具体而言，在实验中将重点介绍逆滤波器的特点及其在去除高斯噪声中的作用，并结合维纳滤波器的原理探讨其在信号恢复中的有效性；同时深入探讨最小二乘法在图像恢复中的数学建模基础及其优化求解过程。
内容及要求：
Matlab导入医学影像数据后，在图像中分别加入高斯噪声分布、盐和 pepper 噬菌体分布、乘法性噪声分布以及泊松分布的 Poisson 噪声。通过分析不同参数大小的变化对图像质量的影响来评估算法效果
2、Matlab读入医学影像，使用运动模糊并叠加一定的噪声降质图像
3、Matlab读入医学影像，利用维纳滤波的三种形式复原图像
4、Matlab读入医学影像，利用有约束最小二乘复原图像
二、实验工具
| 仪器名称 | 规格/型号 | 数量 | 备注 |
|---|---|---|---|
| 计算机 | Windows10 | 1 | 有网络环境 |
| Matlab | 2016b及以上 | 1 | |
| Python | 3.6及以上 | 1 |
三、实验步骤及结果
使用Matlab导入医学影像，并施加高斯噪声、椒盐噪声等不同类型的人工干扰信号；研究不同参数强度变化对图像质量的具体影响。
lc；%clc；清空命令行窗口 clear；%clear；删除当前工作区中的所有变量或对象 close all；%关闭所有未保存的图形窗口 f = imread('MR.tif'); %从磁盘中加载MRI切片图像到当前工作区并存储为变量f % 生成不同参数组合下的高斯噪声图像 g6 = imnoise(f,'gaussian',0,0.01); g7 = imnoise(f,'gaussian',0.5, 2sqrt(2)pisqrt(epsilon)); g8 = imnoise(f,'gaussian', 2sqrt(2)pisqrt(epsilon), epsilon); subplot(131), imshow(g6), title('默认值参数下的高斯噪声污染图像'); subplot(132), imshow(g7), title('均值为{mu}={epsilon}、标准差为{sigma}=sqrt({epsilon})的高斯噪声污染图像');
subplot(133),imshow(g8),title('均值为0、方差为0.05的高斯噪声污染的图像');注：
这是使用matlab进行图像处理的代码片段，其中imnoise函数用于向输入图像f添加高斯噪声。具体地，g6、g7和g8分别表示添加了不同参数的高斯噪声后生成的新图像，其中：g6：意味着默认情况下生成的高斯噪声，在这里为均值0、方差0.01。g7：表示均值为0.5、方差为0.01的高斯噪声。g8：表示均值为0、方差为0.05的高斯噪声。subplot函数用于将这三个图像显示在单个图形窗口中的三个子区域中，并添加了相应的标题。这些参数可以用于研究图像处理技术在各种条件下对图像的影响，从而确定最优参数设置以达到所需的效果。
以下是符合要求的改写文本
注：
以下是对输入文本进行了同义改写的版本
% 乘性噪声g6 = imnoise(f,'speckle',0.01);g7 = imnoise(f,'speckle',0.04);g8 = imnoise(f,'speckle',0.08);figure;subplot(131),imshow(g6),title('参数为0.01乘性噪声污染的图像');subplot(132),imshow(g7),title('参数为0.04乘性噪声污染的图像');subplot(133),imshow(g8),title('参数为0.08乘性噪声污染的图像');创建了一个新的图形窗口，并将三个不同的带有不同强度的乘性噪声污染后的图像进行展示
该代码调用Matlab图像处理工具箱中的imnoise函数生成乘性噪声，并将其附加到图像变量f中。其中,g6、g7和g8分别代表使用不同参数生成的噪声图像序列。这些图像随后被组织成一个包含三个子图的图形窗口,每个子图均展示了不同参数下乘性噪声污染后的图像。通过subplot函数,在同一图形窗口中显示了多个子图,并通过imshow函数显示所选区域中的图片, title函数则在每个子图旁添加了相应的标题。这种类型的噪声会导致图像呈现粗粒状纹理并伴有灰度级的变化
% 泊松噪声g9 = imnoise(f,'poisson');figure;
imshow(g9),title('泊松噪声污染的图像');
注：
以下是一段用于在图像f上添加泊松噪声并显示结果的Matlab代码。具体而言,该代码通过调用imnoise函数向输入图像f添加泊松分布类型的随机噪声,其中默认的强度参数决定了噪声的具体分布情况（可以通过指定不同的参数值来调节这一强度）。随后利用imshow函数将增强处理后的高斯模糊滤波结果g9进行可视化展示,并在图中标注相应的处理效果指标数值。值得注意的是,由于该过程涉及的数学模型较为复杂,因此在实际应用中建议结合相应的理论知识进行深入理解与操作实践。
2、Matlab读入医学影像，使用运动模糊并叠加一定的噪声降质图像
f = checkerboard(8); % 生成一个由8×8个正方形组成的测试板
PSF = fspecial('motion',7,45); % 由于运动模糊特性, PSF被定义为空间滤波器
gb = imfilter(f,PSF,'circular'); % 应用运动模糊后,图像边缘采用圆周延伸处理
noise = imnoise(zeros(size(f)),'gaussian',0,0.001); % 在图像基础上叠加均值为0,方差为0.001的高斯噪声
g = gb + noise; % 添加高斯噪声至模糊图像以模拟真实观测结果
figure;
subplot(1,3,1);imshow(g);title('运动模糊+高斯噪声输入图像(g)');
subplot(1,3,2);imshow(fr1);title('基于全范围正则化的重建结果');
subplot(1,3,3);imshow(fr2);title('基于缩小范围正则化的重建结果');
注：
这段代码是一个MATLAB的图像处理任务，其步骤如下：（1）生成一个8*8的棋盘格测试板。 f = checkerboard(8);（2）产生一个运动模糊的点扩散函数(motion)，滤波器长度为7，方向为45度。 PSF = fspecial('motion',7,45);（3）使用循环卷积实现空间滤波器，并减少边界效应，得到产生退化的图像gb。 gb = imfilter(f,PSF,'circular');（4）通过imnoise()函数，给gb加入高斯噪声，用于构建退化的图像模型。 noise = imnoise(zeros(size(f)),'gaussian',0,0.001);g = gb + noise;（5）分别使用deconvreg()函数进行正则滤波，其中第三个参数为正则化程度，在这里设置为4和0.4。 fr1 = deconvreg(g,PSF,4); fr2 = deconvreg(g,PSF,0.4,[1e-7 1e7]);（6）最后使用subplot()函数展示原始图(g)、第一次正则滤波后的图(fr1)和第二次正则滤波后缩小范围的图(fr2)。 subplot(1,3,1);imshow(g);title('运动模糊+高斯噪声(g)'); subplot(1,3,2);imshow(fr1);title('正则滤波器'); subplot(1,3,3);imshow(fr2);title('正则滤波器(缩小范围)');
清空命令行窗口。
释放工作区内存。
关闭所有打开的图形窗口。
生成一个8x8的棋盘板图像是不必要的。
将MRI图像读取为双精度矩阵并赋值给变量f。
建立一个四子图中的第一个子图，并显示原始图像。
创建运动模糊滤波器。
对模糊滤波器进行卷积处理以生成退化图像。
将退化图像与高斯噪声叠加生成混合信号g。
在四子图中的第四个位置显示混合信号g并设置标题为"模糊加噪声后的信号"。
注：
这段代码实现了以下操作：
读入名为mr.tif的图像，并显示
创建一个大小为7x7，角度为45度的运动模糊滤波psf
将原始图像使用psf进行循环卷积，生成退化图像gb，并显示
创建一个与原始图像大小相同的高斯噪声，方差为0.002
将退化图像gb和高斯噪声相加，生成模糊加噪声后的图像g，并显示
其中，subplot函数用来在同一图形窗口中绘制多个子图；imshow函数用来呈现图片；title函数用来设定图表标题。
3、Matlab读入医学影像，利用维纳滤波的三种形式复原图像
clc; %clc;清除命令行窗口clear; %clear;清除工作区close all; %关闭前面所有窗口f = checkerboard(8); %产生棋盘板图像% f = im2double(imread('MR.tif'));PSF = fspecial('motion',7,45);gb = imfilter(f,PSF,'circular'); %产生退化图像减少边界效应noise = imnoise(zeros(size(f)),'gaussian',0,0.001);%产生高斯噪声g = gb + noise;subplot(221),imshow(g,[]);%模糊加噪声后的图像title('模糊加噪声后的图像');frest1 = deconvwnr(g,PSF); %维纳滤波,信噪比为0的逆滤波subplot(222),imshow(frest1,[]),title('维纳滤波后的结果');Sn = abs(fft2(noise)).^2; %噪声功率谱nA = sum(Sn(:))/numel(noise); %平均噪声功率谱Sf = abs(fft2(f)).^2; %信号功率谱fA = sum(Sf(:))/numel(f); %平均信号功率谱R = nA/fA; %求信噪比frest2 = deconvwnr(g,PSF,R);subplot(223),imshow(frest2,[]),title('使用常数比例的维纳滤波');NACORR = fftshift(real(ifft(Sn)));FACORR = fftshift(real(ifft(Sf)));frest3 = deconvwnr(g,PSF,NACORR,FACORR); %自相关后的维纳滤波subplot(224),imshow(frest3,[]),title('使用自相关函数的维纳滤波');
注：
该代码主要应用于数字信号处理领域的图象修复技术研究中。具体而言,该代码首先通过创建标准棋盘图案并结合运动模糊的点扩散函数（PSF）来生成原始图象数据,随后又在此基础之上叠加高斯噪声以模拟退化后的观测图象情况
在展示了一个退化图像之后，代码应用了维纳滤波器来进行退化图像的修复。同时比较了两种调整信噪比的方法：一种采用将信噪比设为0的逆滤波；另一种则采用了自相关函数来计算信噪比。实际上，在这种情况下，它本质上是从退化图像中提取原始信息的过程，在此过程中利用输入与噪声功率谱以及原始信号功率谱的估算值进行计算。
最终，在代码中使用 subplot 函数将结果显示于四个子图时，则能清晰地观察到不同维纳滤波方法之间的差异。
4、Matlab读入医学影像，利用有约束最小二乘复原图像
clc; %clc;清除命令行窗口clear; %clear;清除工作区close all; %关闭前面所有窗口f = checkerboard(8); %产生棋盘板图像% f = im2double(imread('MR.tif'));PSF = fspecial('motion',7,45);gb = imfilter(f,PSF,'circular'); %产生退化图像noise = imnoise(zeros(size(f)),'gaussian',0,0.001);%产生高斯噪声g = gb + noise;subplot(221),imshow(f,[]);subplot(222),imshow(g,[]);%模糊加噪声后的图像title('模糊加噪声后的图像');f1 = deconvwnr(g,PSF,0.02); %维纳滤波subplot(223),imshow(f1,[]),title('维纳滤波后的结果');% f = edgetaper(f,PSF);f2 = deconvreg(g,PSF,0.4,[1e-7,1e7]); %约束最小二乘滤波subplot(224),imshow(f2,[]),title('使用约束最小二乘滤波');
注：
这段代码执行了对棋盘板图像实施模糊化与去噪处理，并随后应用维纳滤波器与约束最小二乘滤波器来进行图像的复原。具体而言，则包括以下步骤：
清空命令行界面及其工作区，并关闭所有先前的窗口。
使用checkerboard函数生成棋盘板图案或从外部导入一张mr图片，并利用fspecial函数制造一个7×7的运动模糊卷积核以对原始图片进行卷积处理从而得到模糊图片。
生成方差设定为0.001的高斯噪声并将此噪音叠加至模糊图片以生成带有噪音干扰的新图片。
展示原始图片与其添加了高斯噪音后的版本。
调用matlab内置的deconvwnr函数执行维纳滤波以去除模糊影响并恢复出清晰图片并将结果展示在第三个子图区域。
调用约束最小二乘法滤波器（deconvreg）来处理带有噪音的影响并恢复出更清晰的版本并将结果显示在第四个子图区域。
四、讨论与结论
医学图像复原涉及利用计算机技术和数学方法来还原医学图像的详细信息，并从而有助于提升诊断和治疗的效果。以下是一些与这一领域相关的核心知识点：
图像去噪：消除噪声是医学图像复原的核心环节之一。常见的方法涉及中值滤波、高斯滤波以及小波变换等多种技术。
图像增强：采用提升对比度、锐化边缘等方式从而提高图像质量。常见手段包括直方图均衡化（HE）、梯度 sharpening 和基于模板的图像增强等多种技术。
图像恢复：将受退化影响的图像恢复到其原始状态，并通过消除散射光线和运动模糊的影响来实现这一目标。常用的手段包括逆滤波、维纳滤波以及最小均方误差等技术。
三维重建：在医学领域中的图像处理过程中通常会采用三维重建技术以便于观察病灶形态并辅助临床判断。常用的3D成像技术主要有CT断层扫描（CT Scan）、磁共振成像（MRI）、单光子发射断层扫描（SPECT）以及正电子发射断射线（PET）等方法
考虑到医学图像数据量较大，在传输与存储上都面临着诸多挑战。因此，在图像处理流程中实施压缩技术是不可或缺的关键步骤之一。常用的图像压缩方法有哪些？其中包括但不限于JPEG标准及其升级版JPEG 2000以及PNG格式等多种方案。
噪声种类：
- Gaussian noise arises from electronic circuit noise and is caused by low illumination or high temperatures affecting sensor noise.
- Salt-and-pepper noise is generated by imperfections in switching devices.
- The log-normal distribution noise is characterized by the sizes of silver grains in photographic film, which are modeled as random variables following a log-normal distribution.
- Rayleigh noise is due to the characteristics of band imaging technology.
（3）指数和厄兰（伽玛）噪声:在描逑激光成像中的噪声方面是很有用的
图像是可以通过恢复来修正的。在信息形成、记录以及处理和传输的过程中,由于成像系统、记录设备以及处理手段等存在一定的局限性,会导致所获得的信息质量有所降低,这一现象被我们定义为"图像退化"。其典型表现在于信息呈现模糊性、失真性和噪声干扰等特征。
常见的退化因素包括：成像系统的像差或受限的光孔直径；系统中存在的离焦问题；以及景物在运行过程中的相对位置配置等技术指标未满足要求的现象；此外还涉及显示器输出过程中出现的失真现象；遥感过程中的大气散射和扰动效应；还有系统内部各组成部分产生的噪声干扰问题等
就其实质而言，在医学图像复原这一领域中确实具有较高的综合性。就其实质而言，在这一领域中确实具有较高的综合性。为了实现有效的图像复原过程就必须依赖于多方面的学科知识，并且必须要综合运用计算机图像处理技术、信号处理技术和数学等技术手段才能达到理想的效果。
学生实验报告
| 学生姓名 | 景涛 | 学号 | B20040301002 | 组员： 傅云昊B20040301048 李佳一B20040301052 户景仁B20040301051 刘若谷B20040301053 |
|---|---|---|---|---|
| 实验项目 | 基于深度学习的图像分割 | |||
| □必修 □选修 | □演示性实验 □验证性实验 □操作性实验 □综合性实验 | |||
| 实验地点 | E405 | 实验日期 | 2023.5 |
一、实验目的及要求
1、目的
深入学习深度学习模型：包括全连接神经网络、卷积神经网络和U-网架构；全面理解其架构与计算原理等理论知识
2、内容及要求
1 利用U-net模型完成图像的分割
2 重点阐述模型的搭建与训练，可以分隔缺陷
3 数据集是有关 “钢铁表面缺陷”的图像
二、实验工具
| 仪器名称 | 规格/型号 | 数量 | 备注 |
|---|---|---|---|
| 计算机 | Windows10 | 1 | 有网络环境 |
| Matlab | 2016b及以上 | 1 | |
| Python | 3.6及以上 | 1 |
三、实验步骤及结果
1 导入库
import numpy as np 用于在Python中进行科学计算的基础库
import matplotlib.pyplot as plt 一个用于绘制数据图表的Python库
import pandas as pd 用于数据分析和处理的Python库
import os 提供了一些与操作系统交互的函数，可以用于文件和目录操作
from tqdm import tqdm_notebook 该库提供实时跟踪程序运行状态的功能
import cv2；OpenCV是一个功能强大的Python库；它提供了丰富的图像处理和计算机视觉工具
import Keras 也采用了Keras深度学习框架中的相关组件及其函数模块
from keras.layers.convolutional import Conv2DTranspose
from keras.layers import concatenate
#from keras.layers.merge import concatenate
从Keras的层导入上采样、卷积、激活函数等
from keras import Model
from keras import backend as K
from keras.layers.core import Lambda
注：
这段代码是一段Python代码，主要包括以下模块和库：
numpy：一个用于进行科学计算的基础库。
matplotlib.pyplot：一个用于绘制数据图表的Python库。
pandas：用于数据分析和处理的Python库。
os：提供了一些与操作系统交互的函数，可以用于文件和目录操作。
tqdm.tqdm_notebook：一个Python的进度条库，用于显示脚本运行时的进度。
cv2：OpenCV的Python接口，集成了许多图像处理和计算机视觉相关的函数。
keras：深度学习框架，用于构建神经网络模型。
该层在Keras中用于实现反卷积操作，在卷积神经网络中的应用非常广泛
keras中的连接层：该函数用于整合或连接多个输入张量
keras.layers.UpSampling2D：keras中的上采样层，用于增加图像的分辨率。
keras.layers.Conv2D：keras中的卷积层，用于提取特征。
keras.layers.Activation：在深度学习框架Keras中定义的激活函数层，在神经网络模型中作为关键组件用于对输入数据特征进行非线性转换处理
keras.layers.Input：keras中的输入层，用于定义网络的输入。
keras.layers.Dropout：keras中的随机失活层，用于减少过拟合。
keras.layers.MaxPooling2D：keras中的最大池化层，用于降低图像的尺寸。
keras.Model：keras中的模型类，用于定义神经网络的结构。
keras.backend：keras中的后端模块，用于执行计算图的操作。
keras.layers.core.Lambda：keras中的Lambda层，用于自定义函数层。
2 读取数据
tr = pd.read_csv('E:/imageprocess/Steel Defect Detection/train.csv')
读取csv文件并转换为pandas数据框
print(len(tr)) # 输出数据框的长度，即有多少行数据
tr.head() # 显示数据框的前5行数据
df_train = tr[tr['EncodedPixels'].notnull()].reset_index(drop=True)
筛选有标注信息的行，并重置索引
print(len(df_train)) # 输出筛选后的数据框长度
df_train.head() # 显示筛选后的数据框的前5行数据
注：
这段代码负责处理一个CSV文件（'E:/imageprocess/Steel Defect Detection/train.csv'），并将其解析为一个包含多行数据的pandas DataFrame对象。
初始代码行调用pandas库中的read_csv()函数以读取指定的CSV文件内容，并将获取的数据结果赋值给tr变量进行后续处理操作。其中，在pandas库中read_csv方法用于执行数据导入操作的具体功能与实现细节，请参考官方文档获取详细说明。
'E:/imageprocess/Steel Defect Detection/train.csv'：csv文件的文件路径；
delimiter：csv文件中的字段分隔符，默认是','；
header：指定csv文件的首行是否为列名称，默认值为0，表示第0行不是列名；
index_col：指定哪些列是索引列；
dtype：指定每列的数据类型。
第二行代码通过调用print()函数显示tr字段中数据框的总行数。
第三行代码涉及DataFrame的基本操作——根据特定条件筛选出符合条件的数据。该运算tr['EncodedPixels'].notnull()将返回一个布尔型的一维数组（Series），每个元素对应原始数据中的相应位置。然后利用该布尔数组作为筛选条件，在原始数据中提取出所有EncodedPixels值不为空的数据行。最后通过reset_index(drop=True)方法来重置df_train数据框的索引结构，默认生成新的递增整数序列并舍弃原有索引信息
第四行代码通过调用print()函数获取或显示df_train这个数据框中的行数，并说明该数据框共有多少条记录。
第五行代码采用DataFrame.head()函数来呈现df_train数据框的前五行数据，并有助于直观了解其数据布局。
3 数据处理
（1）图像处理
def rle2mask(rle, imgshape):
width = imgshape[0]
height= imgshape[1]
#根据输入图像的大小，创建一张全零的二维矩阵（掩膜图像）。
mask= np.zeros( width*height ).astype(np.uint8)
将RLE编码字符串转换为一个整数数组，在该数组中偶数索引对应的元素具体说明了像素值的起始位置信息，并且奇数索引对应的元素则详细说明了像素值的长度信息。
array = np.asarray([int(x) for x in rle.split()])
starts = array[0::2]
lengths = array[1::2]
#遍历整数数组中的每个元素，依次在掩膜图像上标记出像素值为“1”的区域。
current_position = 0
for index, start in enumerate(starts):
mask[int(start):int(start+lengths[index])] = 1
current_position += lengths[index]
#对生成的掩膜图像实施垂直翻折和平移逆时针旋转处理流程,实现预期效果的掩膜图像。
return np.flipud( np.rot90( mask.reshape(height, width), k=1 ) )
def keras_generator(batch_size):
while True: #函数使用一个无限循环语句while True来不断生成训练数据。
x_batch = []
y_batch = []
img_size = 256
for i in range(batch_size):
fn = df_train['ImageId'].iloc[i].split('_')[0]
#对于每个数据样本来说，在函数中会逐个读取其对应的图像文件，并将其转换为RGB格式。
img = cv2.imread( 'E:/imageprocess/Steel Defect Detection/train_images/'+fn )
img = cv2.cvtColor(img, cv2.COLOR_BGR2RGB)
mask = rle2mask(df_train['EncodedPixels'].iloc[i],img.shape)
img = cv2.resize(img, (img_size, img_size))
mask = cv2.resize(mask, (img_size, img_size))
x_batch += [img]
y_batch += [mask]
通过 numpy 将生成的输入图像序列以及对应的标签序列转为数组，并利用 yield 机制将它们以批处理的形式返回。
x_batch = np.array(x_batch)
y_batch = np.array(y_batch)
yield x_batch, np.expand_dims(y_batch, -1)
for x, y in keras_generator(4):
break
print(x.shape, y.shape)
注：
这段代码的主要组成部分包括两个核心功能模块：一个是基于Run-Length Encoding转置矩阵的实现（rle2mask），另一个是Keras数据生成器（keras_generator）。此外还包括用于测试数据生成器的一段辅助代码。
rle2mask()函数接受两个参数：rle和imgshape. 其中,rle是以字符串形式表示的Run Length Encoding编码结果,用于表示一张图像的掩膜; imgshape是以元组形式表示的图像尺寸,用于确定掩膜对应的图像大小.
该函数首先基于输入图像的尺寸生成一个全零的二维矩阵（称为掩膜图像）。接着将RLE编码字符串解析为一个整数数组，在此过程中偶数索引处存储的是像素值的起始位置信息，在奇数索引处存储的是像素值持续存在的长度信息。随后遍历该整数数组中的每一个元素，在对应的掩膜图像坐标位置上标记出相应的像素点区域设为1。最后对该掩膜图像执行垂直翻转和逆时针旋转操作以获得最终结果。
该函数接收名为batch_size的参数。
该函数通过不断循环调用while True机制完成数据生成工作。
每个数据样本依次读取与其对应的图像文件，并将其转换为RGB格式。
随后将该样本的掩膜利用rle2mask()函数转换为掩膜图像。
接着使用cv2.resize()函数对图像与掩膜进行尺寸调整（img_size = 256）。
最后将生成的一系列输入图像列表与标签列表转换为numpy数组，并通过yield关键字输出一批次的数据内容。
在测试过程中使用的代码调用keras_generator()函数，并从该生成器中获取一批大小为4的数据样本。随后将这些样本赋值给变量x和y，并从该生成器中获取一批大小为4的数据样本。随后打印变量x和y的形状信息以查看每个样本的具体结构信息。
（2）读图像
plt.imshow(x[3])
plt.imshow(np.squeeze(y[3]))
注：
这段代码调用Matplotlib库中的imshow()函数以呈现两张图像。第一个图片来自变量x中索引为3的位置；第二个图片来自变量y中索引为3的位置，并通过调用np.squeeze()函数去除掩膜图像中多余的维度之后的结果即是y数组中索引为3的部分。
需要注意的是，绘制图像时需要使用plt.show()函数将其显示出来。
3 模型搭建
（1）模型搭建采用
卷积层 + dropout层 + 卷积层 + 池化层(最大)
由 U-Net 结构构成的卷积神经网络模型具有识别输入彩色图像中每一个像素是否存在缺陷的功能。该模型由编码（downsampling）模块和解码（upsampling）模块组成，在逐层进行卷积和下采样操作后能够提取出图像的抽象特征，并通过上采样将其整合为精细的掩膜图像。
该模型接收一个大小设定为 (256×256×3) 的三维张量作为输入数据，在其中存储了一幅彩色图像的信息。经过 Lambda 层处理后将输入数据的标度范围调整至 (0,1) 区间内。随后通过五个相同的卷积-下采样模块进行处理，在每个模块中包含两个连续的卷积层、一个整流线性单元(ELU)作为激活函数以及一个从 0.1 到 0.3 范围内逐渐减小的失活层（Dropp layer）。每个卷积-下采样模块都负责提取输入图像中的抽象特征信息，并使输出特征图的空间尺寸减半。
然后，在解码过程中采用了四个相同的反卷积（上采样）模块组成了编码回路结构，并包含了一个转置卷积层以及一个跳跃连接层。这些组件的作用是将编码器阶段提取出的关键信息与解码器生成的特征图进行融合处理，在逐步放大生成图像尺寸的同时实现了对目标区域二进制掩膜信息的有效捕捉。其中，在解码器中嵌入了跳跃连接机制以实现跨尺度特征融合，在此基础之上结合多尺度上下文信息能够有效提升网络识别小型物体及显著边界的能力
最后，在经过卷积层处理后，神经网络借助sigmoid激活函数的作用使得输出特征图被限制在0到1之间范围之内。该模型采用二分类交叉熵损失函数并配合Adam优化器进行训练优化过程
inputs = Input((256, 256, 3)) # 输入层：256x256 大小的 3 通道图像
s = Lambda(lambda x: x / 255) (inputs) # 归一化，将像素值除以255
编码器，提取特征
c1 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (s)
c1 = Dropout(0.1) (c1)
c1 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c1)
p1 = MaxPooling2D((2, 2)) (c1)
c2 = Conv2D(32, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (p1)
c2 = Dropout(0.1) (c2)
c2 = Conv2D(32, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c2)
p2 = MaxPooling2D((2, 2)) (c2)
c3 = Conv2D(64, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (p2)
c3 = Dropout(0.2) (c3)
c3 = Conv2D(64, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c3)
p3 = MaxPooling2D((2, 2)) (c3)
c4 = Conv2D(128, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (p3)
c4 = Dropout(0.2) (c4)
c4 = Conv2D(128, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c4)
p4 = MaxPooling2D(pool_size=(2, 2)) (c4)
c5 = Conv2D(256, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (p4)
c5 = Dropout(0.3) (c5)
c5 = Conv2D(256, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c5)
解码器，将特征图还原为原始图像尺寸，并进行像素级别的分类
u6 = Conv2DTranspose(128, (2, 2), strides=(2, 2), padding='same') (c5)
u6 = concatenate([u6, c4])
c6 = Conv2D(128, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (u6)
c6 = Dropout(0.2) (c6)
c6 = Conv2D(128, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c6)
u7 = Conv2DTranspose(64, (2, 2), strides=(2, 2), padding='same') (c6)
u7 = concatenate([u7, c3])
c7 = Conv2D(64, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (u7)
c7 = Dropout(0.2) (c7)
c7 = Conv2D(64, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c7)
u8 = Conv2DTranspose(32, (2, 2), strides=(2, 2), padding='same') (c7)
u8 = concatenate([u8, c2])
c8 = Conv2D(32, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (u8)
c8 = Dropout(0.1) (c8)
c8 = Conv2D(32, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c8)
u9 = Conv2DTranspose(16, (2, 2), strides=(2, 2), padding='same') (c8)
u9 = concatenate([u9, c1], axis=3)
c9 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (u9)
c9 = Dropout(0.1) (c9)
c9 = Conv2D(16, (3, 3), activation='elu', kernel_initializer='he_normal', padding='same') (c9)
# 输出层，输出二进制分类结果
outputs = Conv2D(1, (1, 1), activation='sigmoid') (c9)
定义模型
model = Model(inputs=[inputs], outputs=[outputs])
编译模型，指定优化器和损失函数
model.compile(optimizer='adam', loss='binary_crossentropy')
注：
这段代码是定义了一个 UNet 模型用于图像分割，主要分为两部分：
（1）定义模型结构：
输入层：256*256 大小的 RGB 图像；
归一化：将像素值除以 255，使所有像素值在 [0, 1] 之间；
编码器：五个卷积层和池化层，用于提取输入图像的特征；
解码模块：该模块包含四个反卷积层与卷积层，并负责将特征图还原至与原始图像相同的空间维度，并完成像素级别的分类任务；
输出层：1 个卷积层，输出二进制分类结果。
（2）模型编译：
优化器：Adam；
损失函数：二元交叉熵。
（2） 模型训练
%%time
Fit model
这段代码使用给定的生成器和超参数训练模型
batch_size = 16 # 每个批次的大小为16
results = model.fit_generator(keras_generator(batch_size),
steps_per_epoch=100,
epochs=5)
使用keras_generator(batch_size)作为训练数据生成器
使用keras_generator(batch_size)作为训练数据生成器
每个epoch分为100步来迭代训练数据
训练5个epoch
注：
这段代码调用fit_generator函数来训练一个Keras模型，并且其中keras_generator函数负责生成用于Keras模型训练的数据集。
steps_per_epoch参数被用来设置每个epoch中安排多少次训练步骤。这实际上决定了每次从生成器中取样的数量。在这里我们将其设定为每16次取样一个批次共运行100次小步骤以完成一个完整的epoch循环
epochs 参数指定要训练多少个 epoch。在这里，我们训练了5个 epoch。
（3）模型预测
该代码运用Keras库中的深度学习模型预测函数对该输入数据集x进行预测运算，并展示第三个预测结果的可视化分析。
pred = model.predict(x) # 对输入数据进行预测，返回预测结果
plt.imshow(np.squeeze(pred[2])) # 可视化第三个预测结果
注：
此段代码采用Keras模型中的predict函数对输入数据x进行推断，并通过图像展示预测结果中的第三个样本。其中,predict函数能够处理输入数据并输出相应的预测结果。np.squeeze函数可用于去除数组中所有维度大小为1的那一部分。plt.imshow函数可实现图像的可视化展示。
4 测试
（1）测试数据预处理
该代码解析测试目录中的全部图片文件，并对每张图片进行尺寸调整和推断，最终返回推断结果的数量。
通过调用os模块的listdir方法获取指定目录下的所有文件列表到变量test_files中
获取测试图片目录下的所有文件名
len(testfiles) # 输出测试图片的数量
img_size = 256 # 缩放后的图片尺寸
test_img = [] # 定义一个列表存储缩放后的测试图片数据
for fn in tqdm_notebook(testfiles): # 遍历每个测试图片文件名
img = cv2.imread('E:/imageprocess/Steel Defect Detection/test_images/' + fn) 用于加载测试图像文件
img = cv2.resize(img,(img_size,img_size)) # 缩放测试图片为指定的大小
test_img.append(img) # 将缩放后的图片数据添加到列表中
Assign predict_img as the output of model.predict applied to the array representation of test_img.
print(len(predict_img)) # 输出预测结果的长度
注：
这段代码负责处理测试图片目录下的所有图片文件，并在获取这些文件后对其进行缩放。随后利用训练好的Keras模型对这些测试图片进行预测。最终输出预测结果的数量。
第一步，在代码运行时会调用os.listdir()函数来获取测试目录下的所有文件名。随后，在循环迭代每个文件名的过程中会执行以下操作：读取每一张测试图片并对其进行尺寸缩放处理后得到的图像数据会被添加到test_img列表中。接着，在循环结束后会利用模型中的predict函数来对所有处理后的图像执行预测操作，并将结果存储在predict_img变量中。为了方便后续的操作，在这一步骤之后会通过numpy.asarray()函数将test_img列表转换成数组格式以供后续处理。最后，在整个流程结束后代码会打印出预测结果的数量即为参与评估的测试图片数量。
这段代码根据预测结果生成 RLE 编码，并可视化原始测试图片。
pred_rle = [] # 定义一个列表存储每个预测结果的 RLE 编码
for img in tqdm_notebook(predict_img): # 遍历每个预测结果
img = cv2.resize(img, (1600, 256)) # 将预测结果缩放为原始图片大小
tmp = np.copy(img)
tmp[tmp<np.mean(img)] = 0 # 根据阈值转换为二值图像
tmp[tmp>0] = 1
pred_rle.append(mask2rle(tmp)) # 将预测结果的 RLE 编码添加到列表中
def mask2rle(img):
tmp = np.rot90( np.flipud( img ), k=3 ) # 对预测结果进行旋转和翻转操作
rle = []
lastColor = 0;
startpos = 0
endpos = 0
tmp = tmp.reshape(-1,1)
for i in range( len(tmp) ):
if (lastColor==0) and tmp[i]>0:
startpos = i
lastColor = 1
elif (lastColor==1)and(tmp[i]==0):
endpos = i-1
lastColor = 0
rle.append( str(startpos)+' '+str(endpos-startpos+1) )
将连续像素的起始位置和数量以字符串形式添加到 RLE 编码列表中
return " ".join(rle) # 将列表转换成字符串
读取变量img_t并使用OpenCV库中的imread函数加载指定的测试图像文件
读取测试图片
注：
这段代码将预测结果转换为 RLE 编码，并可视化原始测试图片。
首先，程序依次处理每一个预测结果。对每个预测结果按照原始图片的比例进行缩放处理，并通过设定阈值将其转为二值图像。
接着, 代码利用 mask2rle() 函数创建 RLE 编码, 该函数对预测结果执行旋转和平移操作, 并将连续像素的起始位置和数量以字符串形式附加到 RLE 编码列表中. 最后将所有预测结果的 RLE 编码附加到 pred_rle 列表中.
在处理流程的末尾段落中，最终应用 OpenCV 库中的 cv2 imread 函数来加载原始测试图像文件，并通过 matplotlib 库中的 imshow 方法来完成图像展示
mask_t = rle2mask(pred_rle[4], img.shape)
plt.imshow(mask_t)
注：
这段代码基于第五张测试图片的预测结果进行RLE编码，并对所得结果进行视觉化处理。
rle2mask() 函数利用 RLE 编码解析机制来生成二维掩膜矩阵。其中 pred_rle[4] 作为第五张测试图片预测输出的 RLE 编码标识符，在构建掩膜过程中起到了关键作用。 img.shape 则是原始测试图片的空间维度参数，在生成最终掩膜时需要提供这一信息以确定其尺寸大小。
代码最后使用 plt.imshow() 函数对生成的掩膜进行可视化。
sub = pd.read_csv(r'E:\imageprocess\Steel Defect Detection\sample_submission.csv')
sub.head()
注：
为了实现从指定路径 'E:/imageprocess/Steel Defect Detection/sample_submission.csv' 中读取数据的目标，并将其成功加载到变量 sub 中。随后可以通过调用 sub 的 head 方法来查看数据的前几行以确保读取过程顺利完成。
%%time
for fn, rle in zip(testfiles, pred_rle):
sub['EncodedPixels'][sub['ImageId'].apply(lambda x: x.split('_')[0]) == fn] = rle
img_s可以通过OpenCV库中的imread函数加载位于文件路径E:\imageprocess\Steel Defect Detection下的test_images文件夹中特定的图像文件
plt.imshow(img_s)
mask_s = rle2mask(sub['EncodedPixels'][16], (256, 1600))
plt.imshow(mask_s)
注：
在当前行中，调用 Jupyter Notebook 的魔术命令 %%time 来测量整个代码块的运行时长。
for循环依次遍历testfiles和pred_rle列表中的对应元素对（包含文件名及其预测结果），逐项对sub表格中对应位置的EncodedPixels字段进行数据更新操作；具体而言，在每次迭代过程中会将当前记录中EncodedPixels字段的值设置为相应的rle编码结果。这里所使用的lambda函数的功能在于通过对文件名进行分解处理，并将其与fn参数进行对比匹配，在此过程中自动识别出哪些行目的实现达到了预期效果并应予以保留
读取第 17 行图片文件，并使用 plt.imshow() 方法显示图片。
将 sub['EncodedPixels'][16] 的参数值传递给 rle2mask() 方法以生成相应的遮罩图像。其中的第一个参数是 EncodedPixels 类型的字符串；第二个参数指定的是目标遮罩图像的尺寸设置；具体来说，则是以高度 256 和宽度 1600 的形式进行定义；该方法生成的结果将被存储在变量 mask_s 中；最后一步使用 plt.imshow() 方法显示了所生成的遮罩图像。
四、讨论与结论
基于深度学习的图像分割相关的重点：
基本思想：将输入图像划分为若干具有特征的区域部分，并使每个区域内部像素特性具有显著一致性；同时确保不同区域之间像素特性存在明显差异性。
传统方法：采用基于预设规则或特征提取的方式进行图像分割操作；如经典的阈值分割方法、边缘检测算法以及基于区域增长模型的操作流程等；这些经典算法在处理复杂背景环境以及目标细节刻画等方面均存在一定的局限性。
深度学习方法：构建全连接网络模型进行图像分割操作；该类方法通过卷积神经网络实现端到端的学习过程；能够自适应地提取并表征图像特征信息；从而实现高质量分割效果。
数据集：深度学习算法体系需要建立高质量的数据训练集合；目前广泛使用的图像分割标准数据集包括PASCAL VOC、COCO以及ADE20K等多种类型。
损失函数：设计有效的损失函数计算机制成为模型训练的重要环节；此类函数用于评估预测结果与真实标签之间的差异程度；常用损失函数包括交叉熵损失函数以及Dice系数损失函数等。
数据增强：通过一系列图像变换操作手段来提升训练样本多样性；主要采用图像旋转操作、随机裁剪处理以及尺寸缩放等技术手段来实现数据增广目的。
U-Net模型的具体实现过程如下：第一步是利用编码器模块对输入图像进行特征提取，在这一阶段中提取出关键的特征信息；接着是经过解码器模块的处理后进行上采样操作；最后通过这一系列步骤生成与原始图像尺寸一致的输出图像。在这一阶段中采用拼接（Concat）技术，在channel维度上将各层特征进行融合，并最终形成厚度更大的复合特征¹³。
受限于我们的计算资源有限，在完成五轮实验后（约80%的成功率），我们的团队初步实现了图像分割功能。在深入研究U-Net模型时，我们认识到其高效的网络架构设计及其精确度令人印象深刻。传统的U-Net架构通过镜像对称编码器和解码器模块构建，并通过跳跃连接整合编码层与解码层的信息。针对输入的高分辨率图像，在处理过程中系统能够有效地提取并识别不同尺度的空间特征，并输出生成精确的空间分割结果。即使面对包含大量噪声干扰的数据集表现依然出色。
在实践中发现,U-Net模型具有广泛的应用潜力,能够在多个领域中发挥作用,包括医疗图像分割、语义分割以及半监督学习中的图像分割等场景。此外,该模型因其网络架构设计科学且简洁明了,使得在图像分割应用中能够快速部署
总体而言而言,U-Net model is a widely-used deep learning architecture for image segmentation tasks.The fundamental structure of U-Net model is inspired by autoencoder networks, typically comprising two key components: an encoder module and a decoder module.The encoder module systematically compresses the input image information to generate a feature vector, while the decoder module then reconstructs this feature vector into an output image that matches the spatial resolution of the original input data.This innovative design enables precise pixel-level annotations in medical imaging applications, among other domains, achieving state-of-the-art performance across various segmentation tasks as demonstrated in prior studies¹³