Matlab图像处理 参考： https://www.cnblogs.com/tiandsp/p/14337176.html
1、图像处理基本函数 MATLAB图像处理函数：读取、写入、显示和修改图像
https://ww2.mathworks.cn/help/matlab/images_btfntr_-1.html
显示图像
读取、写入和修改图像
转换图像类型
修改图像颜色
属性 Image Processing Toolbox 本工具箱可执行图像处理、可视化和分析 https://ww2.mathworks.cn/help/images/index.html
Image Processing Toolbox™ 为图像处理、分析、可视化和算法开发提供了一套全面的参考标准算法和工作流 App。您可以使用深度学习和传统图像处理技术执行图像分割、图像增强、去噪、几何变换和图像配准。工具箱支持处理二维、三维和任意大的图像。
Image Processing Toolbox App 可让您自动完成常见的图像处理工作流。您可以交互方式分割图像数据、比较图像配准技术，以及对大型数据集进行批量处理。利用可视化函数和 App，您可以探查图像、三维体和视频；调整对比度；创建直方图；以及对关注区域 (ROI) 执行操作。
matlab影像处理 一、Introduction to digital image（介绍数字影像）
（一）Digital Image and Its Acquisition（数字图像及其采集）
1、An image is an artifact that depicts or records visual perception（图像是描述或记录视觉感知的人工制品）
2、Typically acquired by using charge-coupled device（CCD）or complementary metal-oxide-semiconductor（CMOS）device（通常通过电荷耦合器件（CCD）或互补金属氧化物半导体（CMOS）器件获得）
（二）Where Does the Color Come from？
1、Three kinds of light-sensitive photoreceptor cells in the human eye（i.e.，cone cells）respond most to red，green and blue（人眼中三种感光细胞（即视锥细胞）对红、绿、蓝的反应最为强烈）
（三）Typical RGB Image
（四） 图像在MATLAB中的存储格式Type of Digital Image
1、Binary（黑白图）：Each pixel is just black or white（每个像素都是黑白的）
2、Grayscale（灰度图）：Each pixel is a shade of gray，normal from 0（black）to 255（white）（每个像素都是灰色，数值一般在0到255之间）
3、True color or RGB（彩图）：Each pixel has a particular color described by the amount of red，green and blue in it（每个像素都有一种特定的颜色，用其中的红、绿和蓝的数量来描述）
（五）Elements of Images（图像元素）
1、Binary（黑白图）：
2、Grayscale（灰度图）：
3、True color or RGB（彩图）：
二、Read and show write image（图片的读取、显示、写入）
（一）读取与写入图像指令
1、读取图片：imread（）
2、显示图片：imshow（）
示例代码：
clc;
close all;
I = imread('Pout.tif');%read
imshow(I);%show 输出结果：
3、写入图片：imwrite（） 将图像写入图形文件
语法 imwrite(A,filename)
说明
imwrite(A,filename) 将图像数据 A 写入 filename 指定的文件，并从扩展名推断出文件格式。imwrite 在当前文件夹中创建新文件。输出图像的位深取决于 A 的数据类型和文件格式。对于大多数格式来说：
如果 A 属于数据类型 uint8，则 imwrite 输出 8 位值。
如果 A 属于数据类型 uint16 且输出文件格式支持 16 位数据（JPEG、PNG 和 TIFF），则 imwrite 将输出 16 位的值。如果输出文件格式不支持 16 位数据，则 imwrite 返回错误。
如果 A 是灰度图像或者属于数据类型 double 或 single 的 RGB 彩色图像，则imwrite 假设动态范围是 [0,1]，并在将其作为 8 位值写入文件之前自动按 255 缩放数据。如果 A 中的数据是 single，则在将其写入 GIF 或 TIFF 文件之前将 A 转换为 double。
如果 A 属于 logical 数据类型，则 imwrite 会假定数据为二值图像并将数据写入位深为 1 的文件（如果格式允许）。BMP、PNG 或 TIFF 格式以输入数组形式接受二值图像。
如果 A 包含索引图像数据，则应另外指定 map 输入参数。
如：将一个 50×50 的灰度值数组写入当前文件夹中的 PNG 文件。
A = rand(50);
（二）Image Variable in Workspace（工作区中的图像变量）
1、示例代码：
whos 输出结果：
2、示例代码：
clc; 输出结果：
左边为原图，右边为偶数位置置0后的图片（这些位置为黑色）
（三）Image Info：显示影像信息
示例代码：
imageinfo('pout.tif') 输出结果：
（四）Image Viewer（图像查看器）
Get pixel information in image viewer（在图像查看器中获取像素信息）
示例代码：
imtool('pout.tif') 输出结果：
（五）Image Processing（图像处理）
1、Any form of signal processing for which the input is an image（输入为图像的任何形式的信号处理）
注意：
noisy lena（透镜噪声）
Gaussian filter（高斯滤波）
median filter（中值滤波器）
Wiener filter（维纳滤波器）
三、Image arithmetic（影像四则运算）
（一）影像四则运算函数
（二）Image Multiplication：图像乘法 i mmultiply（）
示例代码：
I = imread('rice.png'); 输出结果：
Excise：How to reduce the brightness of the image？
答案代码：
I = imread('rice.png'); 输出结果：
（三）Image Addition：图像加法 imadd（）
示例代码：
I = imread('rice.png'); 输出结果：
Practice
1、Adjust the “brightness”and“contrast”of rice.png and display it on the screen（调整rice.png 的“亮度”和“对比度”，并将其显示在屏幕上）
答案代码：
I=imread('rice.png'); 输出结果：
（四）Image Histogram：图像直方图 imhist（）
示例代码：
I = imread('Pout.tif');%read 输出结果：
======================================
直方图
直方图就是用来统计一幅图像各个亮度的像素的个数，并在一个2维图像上显示。横向（x轴方向）是灰度值，最左边像素值为0（黑），最右边像素值为255（白），纵向（y轴方向）是各像素值在图像中出现的个数。一张照片的明暗可以通过直方图看出来，也就是说我们可以通过直方图知道图像的欠曝和过曝。对于欠曝，形象的可以理解图像太暗了，过曝可理解为图像太亮了。
需要了解的是，相机在记录信息的时候会同时产生噪点。在较亮区域，因为记录下的光线多，所以噪点并不明显。而在较暗区域，相机记录的光线信号很少，噪点就会一下子突显出来。至于为什么在较暗的区域噪点会更明显的原因主要是，因为信噪比（芯片收到的光线信号与芯片本身噪点数量的比值）太低。所以在这种情况下，如果后期再对暗部提亮，噪点就会更加明显。
直方图如下：
那具体如何从直方图中得到这些信息？就需要简单地聊一聊。 直方图的x轴方向还可以被进一步地划分为5个区域，分别是：黑色，阴影，中间调，高光和白色。
一张理想的图像应该是直方图堆积在中部，最左侧和最右侧都没有被切断。切断或者溢出，是指在直方图左右两个边缘，有很高的柱子堆积，也就是存在很多像素点在左（过暗）右（过亮）两边。 一幅 正常 的图像的直方图：
一幅 欠曝 图像的直方图：在视觉上看整体偏暗。而在直方图中看，会看到最左边的像素很多。
一幅 过曝 图像的直方图：在视觉上看整体很亮。而在直方图中看，最右边的像素（255附近）很多，使得那一块的柱形很高。
灰度直方图的作用 任何一幅图像都有灰度直方图，但相同的灰度直方图可能对应不同的图像，因为在形成灰度直方图的时候丢失了空间信息，即你并不能知道该灰度对应于图像中的位置。
直方图应用：利用直方图进行图像中边界阈值的选择
======================================
Practice
1、Plot the histograms of the images before and after the “brightness”and“contrast”adjustment for rice.png （在rice.png 的“亮度”和“对比度”调整前后绘制图像的直方图）
答案代码：
I=imread('rice.png'); 输出结果：
（五）Histogram Equalization（ 直方图均 衡） ：histeq（）
如果图像的对比度太差，常用的方法就是灰度直方图均衡化。在matlab中，能达到这个目的的函 数就是histeq.那么直方图均衡化的原理是什么呢？下面就主要讲解一下：
histeq的原理：
Matlab中的histeq函数——图像灰度直方图均衡化
首先来看一下没有进行灰度直方图均衡化的图片和经过灰度直方图均衡化后的图片的对比：
左为原图，右为经过灰度直方图均衡化处理后的图片，很显然，经过处理之后的图片在对比度方面有了显著的提升，下面附上两张图片的灰度直方图：
左图为原图像的灰度直方图，右图为处理后的图片的灰度直方图。从灰度直方图上来看，二者有很大的不同，特别是在灰度值处于100-255范围内的直方图，但细细观察我们会发现，右图相较于左图，其实可以看做是左图直方图的扩展，即原图像的直方图范围可能是0-80，经过均衡化后，将其范围扩展到0-255，因此，图像灰度直方图均衡化处理的关键点就在于扩充，即找到原图中某个灰度级的点在处理后图片中的位置并映射过去，从而达到扩充的目的。下面附上matlab源码：
function dst = my_histeq(src)
%histeq函数实现图片灰度直方图均衡化的原理
[m, n] = size(src);
dst = ones(m, n);
% 创建一个和原图像同等大小的矩阵
h = imhist(src);
I = length(h);
PDF = h/numel(src);
%概率密度PDF和概率分布函数CDF
CDF = cumsum(PDF);
j = CDF.*256;
% 取整扩展，结果是均衡化之后的灰度直方图
J = round(j);
for y = 1:m
for x = 1:n
src_y_x = src(y, x);
% 原图中该点的灰度值
dst(y, x) = J(src_y_x + 1);
% 找到原图中该点的灰度值在概率分布函数中的值
end
end
1、Enhances the contrast of the image（增强图像的对比度）
示例代码：
I=imread('pout.tif'); 输出结果：
2、Practice：Write your own equalization function，try it on pout.tif ，and display it on the screen（编写自己的均衡功能，在pout.tif 上尝试，并显示在屏幕上）
（1）Find the new gray level for each pixel（找到每个像素的新灰度）
（2）Write a loop to update the gray level of all the pixels（写一个循环来更新所有像素的灰度）
答案代码：（经典！）
f = imread('pout.tif');ld = find(J == i); %找出扩展后的级数对应的扩展前的级数 %chunchunlei_code20210422 255* (f-P_MIN)/(P_MAX-P_MIN))); 输出结果：
注意：这个不使用histeq（） 重在原理的理解，本人参考以下内容
MATLAB中histeq的原理以及自写的具体实现函数_matlab_yutong5818的博客-CSDN博客
直方图匹配
https://blog.csdn.net/cjsh_123456/article/details/79218179
在上一篇中，我们使用matlab完成了直方图的均衡化，在这一篇中，我们要实现图像的直方图匹配。与直方图均衡不同的一点是，直方图匹配要得到的是具有特定形状的直方图的图像。
当图像中接近0的像素过多，进行直方图均衡的时候，会导致统计概率变大，直接映射到高灰度，直接导致图像整体变亮，
直方图匹配的计算过程：
matlab代码如下：
1. 直方图匹配函数如下：
该匹配函数为histmatching, 输入需要进行直方图匹配的图像，模板直方图，输出为直方图匹配之后的图像，变换向量。
当给出模板图像时，需要在调用该函数之前使用imhist(I)作为函数的模板直方图。该函数使用灰度级别为256的灰度图像。
function [image_out, T3] = histmatching(image_in, hist)
% 直方图匹配（规定化）函数
% 输入为需要进行直方图匹配的灰度图像，模板直方图
% 输出为直方图匹配后的灰度图像，进行变换的向量。
% Level为灰度级别
% T1, T2分别为输入图像，模板直方图的均衡化用到的变换向量
% T3为输入图像匹配模板直方图用到的变换向量
Level = 256;
[m, n] = size(image_in);
image_hist = imhist(image_in);
image_out = image_in;
% 求解T1
ac1 = zeros(Level, 1);
T1 = zeros(Level, 1, 'uint8');
ac1(1) = image_hist(1);
for i = 2 : Level
ac1(i) = ac1(i - 1) + image_hist(i);
end
ac1 = ac1 * (Level - 1);
for i = 1 : 256
T1(i) = uint8(round((ac1(i)) / (m * n)));
end
% 求解T2
ac2 = zeros(Level, 1);
T2 = zeros(Level, 1, 'uint8');
ac2(1) = hist(1);
for i = 2 : Level
ac2(i) = ac2(i - 1) + hist(i);
end
ac2 = ac2 * (Level - 1);
hist_sum = sum(hist);
for i = 1 : 256
T2(i) = uint8(round((ac2(i)) / hist_sum));
end
% 求解T3
% T1映射到T2^(-1)时，若有多个值，选取最小的那个值。
% 产生0 到 255 之间的256个点，即产生0,1,2,...,255的大小为256的数组
temp = zeros(Level, 1, 'uint8');
T3 = T1;
for i = 1 : 256
for j = 1 : 256
temp(j) = abs(T1(i) - T2(j));
end
[~, B] = min(temp);
T3(i) = B - 1;
end
% 根据T3转换输入图像的值
for i = 1 : m
for j = 1 : n
image_out(i, j) = T3(uint32(image_in(i, j)) + 1);
end
end
end
测试代码如下： clear all;
clc;
image2 = imread('2.png');
image3 = imread('3.png');
hist3 = imhist(image3);
match1 = histeq(image2, hist3);
[match2, T] = histmatching(image2, hist3);
figure;
subplot(2, 4, 1), imshow(image2), title('原图像');
subplot(2, 4, 2), imshow(image3), title('模板图像');
subplot(2, 4, 3), imshow(match2), title('匹配后得到的图像');
subplot(2, 4, 4), imshow(match1), title('调用histeq的结果');
subplot(2, 4, 5), imhist(image2), title('原图像的直方图');
subplot(2, 4, 6), imhist(image3), title('模板图像的直方图');
subplot(2, 4, 7), imhist(match2), title('匹配后的直方图');
subplot(2, 4, 8), imhist(match1), title('调用histeq得到图像对应的直方图');
关于用到的matlab函数的一些说明：
1. [A, B] = min(V); V为一个向量，函数返回值A表示V中最小的值，B表示该最小值的位置，当有多个最小值时，返回第一个出现的位置。
2. matlab的函数调用中，对于有多个返回值的函数，例如有2个返回值[A, B]，若在之后的过程中不需要使用到所有的值，可以在接受时将相应位置的参数用~代替。例如代码中对于min(V)函数，我们不需要用到其最小值，只需要用到最小值所在的位置，可以用[~,B] = min(V)。若只需要用到前几个返回值，可以省略后面的。例如：[A] = min(V); 返回的是V的最小值（A = min(V) 与 [A] = min(V)等价）。
3. histeq函数(可以在matlab 通过 help histeq 查看该函数的相关说明)
J = histeq(I, hist); I为输入图像，hist为匹配直方图，可以是用不同方式表示的灰度图像的直方图，例如doule/single类型（横坐标范围0~1），uint8类型（横坐标范围0~255），uint16和int16
J = histeq(I, n); 直方图均衡函数，I为输入图像，n为均衡后的灰度级别，缺省值为64。
[J, T] = histeq(I); T为将灰度图像I中灰度级映射到J中灰度级的灰度变换，值范围为[0, 1]。
另外，还有3个针对索引图像调色板的函数如下：
newmap = histeq(X,map,hgram)
newmap = histeq(X,map)
[newmap,T] = histeq(X,...)
（六）Geometric Transformation（几何变换）：
1、Moving the coordinate（Not the gray-levels）of the pixels in an image（移动图像中像素的坐标（不是灰度））
2、Geometic Transformation Matrices（2D）
https://www.mathworks.com/help/images/performing-general-2-d-spatial-transformations.html
（七）Image Rotation：图像旋转 imrotate（）
示例代码：
I = imread('rice.png'); 输出结果：
Practice：In two dimensions，rotation of a point（x，y）for an angle θ“counter-clockwise”can be written as：（在二维中，一个点（x，y）的旋转角度θ“逆时针”可以写成：）
答案代码：
jiaodu=45; %要旋转的角度，旋转方向为顺时针 输出结果：
注意：一般有 最邻近差值 与 双线性差值 两种方法
最邻近差值：
matlab练习程序（图像旋转，最邻近插值） - Dsp Tian - 博客园
双线性差值：
matlab练习程序（图像旋转，双线性插值） - Dsp Tian - 博客园
（八）Write Image：imwrite（）
1、Format support：
'bmp' 示例代码：
imwrite(I,'pout2.png'); %上面示例的格式均可，不一定png
=============================================================
图像类和类型间的转换
im2uint8 将输入中所有小于0的设置为0,而将输入中所有大于1的设置为255 其他的所有乘以255
im2uint16 将输入中所有小于0的设置为0,而将输入中所有大于1的设置为65535
mat2gray 把一个double类的任意数组转换成值范围在[0,1]的归一化double类数组
im2double 将输入转换为double类.若输入是uint8类 uint16 类 logical类则函数将其转换为范围[0,1]之间的类. 若输入是double类,则函数im2double将返回一个与输入相等的数组.
g=im2bw(f,T) 将一副亮度图像f转换成一副二值图像g ,输出二值图像g中值为0的像素,对应于输入亮度图像f中值小于T的的像素点,输出二值图像g中的1对应于输入亮度图像中大于T的像素点..不管输入是何种数据类图像,T的取值必须在[0,1]内.
亮度变换函数 imadjust 是对灰度图像进行亮度转换的基本IPT工具
此函数将图像f中的亮度值映像到g中的新值,即将low_in至high_in之间的值的映射到low_out和high_out之间的值.low_in以下high_in以上的值被剪切了. 参数gamma给出曲线的形状.该曲线用来映射f的亮度值,以便生成图像g.若gamma小于1,则映射被加权至更高的输出值.
>>A = uint8([235 200 89 20])
mat2gray
功能：将矩阵转化为灰度图像。
>> A = randint(5, 5, [0 255])
mat2gray函数原理分析:
matlab里面很多函数都非常好用，但是当我们需要转化成C的时候就必须理解它的原理，mat2gray这个函数的原理经过单步调式发现其中有个函数imlincombc并没有提供源代码，这里我通过测试几组有意义的数据，基本弄明白了它的原理,具体算法步骤如下： 1、mat2gray是将输入数据F归一化为0-1之间的double型数据
==============================================================
影像处理二 graythresh()
T=graythresh() ,从灰度图像 I 计算一个全局阈值T，使用 Otsu的方法。这个方法是选择一个相对于黑像素和白像素最小的一个方差。全局阈值T可以用来将灰度图像转换为二值图像，使用 imbinarize. 阈值T肯定是属于 [0 1] 区间。
BW= imbinarize( I )
从2D 或 3D灰度图像 I 创建一个二值化图像，通过将所有数值用全局阈值转换为0 或者是1。
im2bw()
BW=im2bw(I ,level) , 将灰度图转换为二值化图像。通过将所有的像素用亮于 level 的变为1， 低于 level 的变为0。
I =imread('rice.png');
level=graythresh(I); %计算最佳阈值
bw=im2bw(I,level); %将灰度图转换为二值化图
subplot(1,2,1); imshow(I);
subplot(1,2,2); imshow(bw);
一、Image thresholding（图像阈值化/二值化）
（一）Problem Setup
1、Count the rice grains and identify their sizes in this image（数一数米粒并确定它们的大小）
2、What are your strategies？
（二）Image Thresholding（图像阈值化/二值化）
1、A gray-level image can be turned into a binary image by using a threshold（使用阈值可以将灰度图像转换为二值图像）
示例代码：
I = imread('rice.png'); 输出结果：
（三）graythresh（）and im2bw（）
1、graythresh（） computes an optimal threshold level（计算最佳阈值级别）
2、im2bw（） converts an images into binary image（将图像转换为二进制图像）
示例代码：
I = imread('rice.png'); 输出结果：
（四）Practice：
1、Write a program to convert the image rice.png into a binary image using a threshold（编写程序，使用阈值将图像rice.png 转换为二进制图像）
2、Do NOT use im2bw（） （不要用函数im2bw ）
3、Try different threshold values to see if you program works（尝试不同的阈值以查看程序的执行情况）
答案代码：
I = imread('rice.png'); 输出结果：
level = 100时的效果
（五）Issue of the Binary Image（二值图像问题）
1、The binary image is not “perfect”（二值化不完美）
（1）Existing sparkles in the background（背景中现有的噪点）
（2）Some grains are missing（一些米粒消失了）
2、What is the cause of this issue？（什么原因导致这个问题）
3、How do we solve it？（如何解决）
二、Background estimation（（在二值化过程中）背景估算）
（一）Background Estimation（预测背景）
1、Estimation for the gray level of the background：
示例代码：
I = imread('rice.png');
BG = imopen(I,strel('disk',15));%strel()与imopen()见注意事项
imshow(BG); 输出结果：
注意：
·strel('disk',15)。创建一个半径为15的圆盘形结构元素
如：
Morphological structuring element（形态学上的结构元素）
https://ww2.mathworks.cn/help/images/ref/strel.html
·imopen(I,strel('disk',15))，
以形态学的方式打开图片
Morphologically open image - MATLAB imopen - MathWorks 中国
（二）Background Subtraction（背景相减）
示例代码：
I = imread('rice.png');
subplot(131);imshow(I);
BG = imopen(I,strel('disk',15));
subplot(132);imshow(BG);
I2 = imsubtract(I,BG);
subplot(133);imshow(I2); 输出结果：
（三）Thresholding on Background Removed Image（将背景移除图像）
示例代码：
I = imread('rice.png');
%% 直接使用graythresh给出的预估值进行二值化
level = graythresh(I);
bw = im2bw(I,level);
subplot(121);imshow(bw);
%% 先对图像进行除背景操作，并将背景与原图进行相减
BG = imopen(I,strel('disk',15));
I2 = imsubtract(I,BG);
%% 使用相减后的图像，然后通过graythresh函数进行二值化
level = graythresh(I2);
bw2 = im2bw(I2,level);
subplot(122);imshow(bw2); 输出结果：
（四）Now What？
1、Want to identify how many grains are there in the image（想知道图像中有多少谷粒）
How？
三、Connected-component labeling（（图像处理过程中）连接点标签）
（一）A procedure for assigning a unique label to each object（为每个对象指定唯一标签的过程）
（二）Connected-component Labeling（Cont'd）
1、Finish labeling of a component（完成组件的标记）
2、Iterative process until all the pixel are checked（迭代过程，直到检查所有像素）
（三）Connected-component Labeling ：bwlabel（）
1、Built-in connected-component labeling algorithm（连通分量标注算法插件）
示例代码：
I = imread('rice.png');
BG = imopen(I,strel('disk',15));
I2 = imsubtract(I,BG);
level = graythresh(I2);
BW = im2bw(I2,level);
[labeled,numObjects] = bwlabel(BW,8);%8联通区域 输出结果：
注意：
·bwlabel(BW,8)函数，用于找联通区域
2、练习1：What is the size of the largest grain？
练习2：What is mean size of the grains？
答案代码：
输出结果：
（四）Color-coding Objects：label2rgb（）
1、Converts a label matrix into an RGB color image
2、Visualize the labeled regions
示例代码：
输出结果：
注意：
·label2rgb(labeled)将标签矩阵转换为RGB图像
（五）Practice
1、Plot the histogram of grain size（绘制米粒大小直方图）
2、Identify all the grains in the image by painting them in red（将图像中的所有颗粒涂成红色，以识别它们）
答案1:
I=imread('rice.png');
BG=imopen(I, strel('disk', 15));
I2=imsubtract(I, BG); level=graythresh(I2);
BW=im2bw(I2, level);
[labeled, numObjects]=bwlabel(BW, 8);
%% 直方图
num_of_lab = zeros(1,numObjects); %预为哈希表先分配空间
for i=1:size(labeled, 1)
for j=1:size(labeled, 2)
index = labeled(i,j);
if index ~= 0 %如果该点不为0，即被标记了，对应该标记的米粒大小加一
num_of_lab(index)=num_of_lab(index)+1; %思想类似哈希表
end
end
end
histogram(num_of_lab)
%% 上色
r = BW * 255; %把白色部分的red分量置为全红
g = zeros(size(BW)); %green分量和blue分量为全0
b = zeros(size(BW));
red_rice = cat(3,r,g,b); %将RGB三个分量拼接为一个三维的彩图矩阵
figure; imshow(red_rice)
答案2:
I=imread('rice.png'); 输出结果：
（六）Object Properties：regionprops（）
1、Provides a set of properties for each connected component（为每个连接的组件提供一组属性）
示例代码：
I=imread('rice.png'); 输出结果：
注意：
·Centroid 中心点
·BoundingBox 边界框
（七）Interactive Selection：bwselect（）
1、Lets you select object using the mouse
示例代码：
I=imread('rice.png'); level=graythresh(I); 输出结果：
======================================================
======================================================
4.使用MATLAB分析图像:目标计数 我们想要通过MATLAB分析rice.png图片中米粒的个数.
要分析图像中的米粒个数,我们需要对图像进行两步预处理:
去除图像的背景:
I = imread('rice.png');
subplot(1,3,1); imshow(I);
BG = imopen(I, strel('disk', 15));
subplot(1,3,2); imshow(BG);
I2 = imsubtract(I, BG);
subplot(1,3,3); imshow(I2);
对图像进行二值化:
I2 = imsubtract(I, BG); level=graythresh(I2);
bw2 = im2bw(I2, level); 下面代码展示了是否去除背景对图像二值化结果的影响:
% 直接对图像进行二值化
I = imread('rice.png');
level=graythresh(I); bw = im2bw(I, level);
subplot (1,2,1); imshow(bw); title('直接进行二值化');
% 去除背景后对图像进行二值化
BG = imopen(I, strel('disk', 15)); I2 = imsubtract(I, BG);
level=graythresh(I2); bw2 = im2bw(I2, level);
subplot(1,2,2); imshow(bw2); title('去除背景后进行二值化');
识别米粒个数的关键在于识别连通区域.
在这里,我们使用MATLAB自带的bwlabel()函数计算连通区域,该函数使用了连通区域标记算法,将每个连通区域内的像素点赋值为同一个值.
I=imread('rice.png');
BG=imopen(I, strel('disk', 15));
I2=imsubtract(I, BG); level=graythresh(I2);
BW=im2bw(I2, level);
[labeled, numObjects]=bwlabel(BW, 8); 得到labeled为标记好的矩阵,其尺寸与原图片相同,每个连通区域都被赋值为一个相同的整数,其他区域被赋值为0.numObjects为计算出的连通区域个数,为99.
使用label2rgb()函数可以将标记结果以彩色图片的形式展示
RGB_label=label2rgb(labeled);
imshow(RGB_label);
使用regionprops()函数可以将检测结果封装成结构体数组.
graindata = regionprops(labeled, 'basic');
graindata(51) Area: 155
Centroid: [112.4258 245.8645]
BoundingBox: [108.5000 234.5000 8 22] 使用bwselect()函数可以交互式选择连通区域
ObjI = bwselect(BW);
imshow(ObjI);
==================================
==================================
Matlab中image、imagesc和imshow函数用法解析
1、显示RGB图像
相同点：这三个函数都是把m*n*3的矩阵中的数值当做RGB值来显示的。
区别：imshow将图像以原始尺寸显示，image和imagesc则会对图像进行适当的缩放（显示出来的尺寸大小）。
2、显示灰度图像
说明：先搞明白什么是索引图像？（灰度图像也是索引图像的一种）
当用Matlab中的imread函数将图像读入并存入矩阵时，我们知道如果是RGB图像，得到是m*n*3的矩阵，但如果是索引图像，得到就是m*n的矩阵，这个矩阵的每个元素只是1个数值，那么怎么确定它的RGB值来显示图像呢？这就需要colormap了，colormap是一个m*3的矩阵，每一行有3列元素构成RGB组，也就是一种颜色，一个m*3的colormap中有m中颜色，而索引图像存储的数值和colormap中的行号对应起来就可以像RGB那样显示图片了，至于对应方法，可以直接对应（比如1对应1,2对应2）也可以是线性映射对应（比如[-128,128]映射到[1,256]）。还有一点要说明的是，默认情况下每一个figure都有且仅有一个colormap,而且默认的是 jet(64)，可在figure窗口通过，edit->colormap...查看，另外在弹出的窗口colormap editor中，可通过Tools->Standard colormap来修改当前figure的colormap，这里是Matlab已经做好的一些colormap。
（1）当灰度图像转化成矩阵后，矩阵中的元素都介于[0,255]，下面我们结合具体实例来看看这三个函数的调用效果，并解释原因。代码：
clear all;clc;close all;
img = imread('lenna.bmp');
% my picture is named lenna.bmp while yours may be not
I = rgb2gray(img);
% Attention: we use the axis off to get rid of the axis.
figure(1),image(I); %equals to imagesc(I,[1 64]);you can try it.
colorbar,title('show by image in figure1');axis off;
figure(2),imagesc(I);
%equals to imagesc(I,[min(I(:)) max(I(:))]);you can try it.
colorbar,title('show by imagesc in figure2');axis off;
%colormap(gray) %use this statement you can get a gray image.
figure(3),imshow(I),colorbar,title('show by imshow in figure3');
显示效果：
我们看到现象是image 和imagesc 显示出来是彩色的，只有imshow 显示出来是灰度图像，为什么会出现这种情况呢？还记得前面所说的吗，索引图像是矩阵和colormap 配合起来显示的，而每个figure 默认使用的colormap 是jet(64) ，而不是gray(gray 和gray(64) 是一样的) ，这个jet(64) 就使得figure1 和figure2 中显示出来时是彩色的，当然你也可以修改当前figure 的colormap 使用colormap(gray) （使用64 个等级的灰度色图），或者colormap(gray(256)) （使用256 个等级的灰度色图，这就是调用imshow 函数时使用的colormap ，后面有讲解）。而figure3 为什么会是灰度图像呢，这是因为当调用imshow 来显示索引图像时，这个函数就会把当前的figure 的colormap 设置成gray(256), 这下明白为什么会出现这种情况了吧。我们再仔细观察一下figure1 和figure2 会发现，figure2 中人物的轮廓显示的还算可以，而figure1 中则出现了大面积的红色的区域，人物的轮廓被抹掉了很多。
colormap 中的索引（ index ）的对应关系了。
image： 这个函数，直接把矩阵中的数当做索引值（我称为直接映射），例如colormap 中索引为 1 的是颜色 RGB1 ，索引为 2 的是颜色 RGB2 ，……，索引为 64 的是颜色 RGB64 。那么矩阵中为 1 的数就显示成颜色 RGB1 ，矩阵中为 2 的数就显示成颜色 RGB2 ，……，矩阵中为 64 的数就显示成颜色 RGB64 ，值得注意的是当矩阵中的数小于 1 时，此时该数也将被显示成颜色 RGB1 ，同样，而矩阵中大于 64 的数将被显示成颜色 RGB64 （类似于信号处理里面的限幅，也可以认为是削顶或者削底了），这下我们就能明白为什么 figure1 中会出现大面积的红色区域，这说明这些地方的数值都大于等于 64 。
imagesc: 在figure2 中我们用 imagesc 来显示图像与 figure1 相比能较好的显示出来，同样我们也得搞明白调用 imagesc 时矩阵的数和 colormap 中索引的对应关系，与 image 不同的是 imagesc 采用的不是直接映射而是线性映射，至于什么是线性映射，我粗略的说一下，比如把区间 A = [0,a] 映射到区间 B = [0,b] 我们对 A 中的元素做 A/a*b 就可以了，矩阵的数到 colormap 索引的线性映射大概就是这样， Matlab 会自动获取矩阵中数的最小值和最大值，并把区间 [Cmin,Cmax] 映射到 colormap 的 [ 最小索引，最大索引 ] 比如 [1,64], 然后再根据这个对应关系把图像显示出来，具体的算法细节是 Matlab 确定的，当然也可以自己指定显示范围，比如一副索引图像 I 范围为 [27,218] ，而我只想显示 [1 64 ], 使用命令 imagesc(I,[1 64]) 就可以了， 如果你把上面程序中的imagesc(I) 换成 imagesc(I,[1,64]) ，那么 figure2 中的效果就和 figure1 中一样了， 因为只是把[1,64] 这个范围映射到色图，超过的都被认为是 64 。关于映射，我截图 Matlab 中 imagesc 的 help 页给大家看看，这里要自己慢慢体会哦，使用 imagesc(I) 这种线性映射就可以用到整个色图从而将图像较好的显示出来，这就是 figure2 中的显示效果比 figure1 中好的原因。
imshow： 调用这个函数会把当前figure 的 colormap 设置成 gray （ 256 ） ，这个前面也有提到，我们先讨论矩阵元素是uint8 型（范围： 0~255 ，整数，一般使用 imread 和 rgb2gray 返回的都是 uint8 型的），同样我们也要搞明白矩阵中的数和 colormap 中颜色索引的对应关系， imshow 的功能是比较全的，它即可使用像 image 那样的直接映射，也可使用像 imagesc 那样的线性映射，当我们使用 imshow(I) ，即只有一个矩阵作为参数，这时采用的是直接映射，比如矩阵中元素0 就显示成 colormap 中索引为 1 的颜色也就是黑色，矩阵中元素 255 就显示成 colormap 中索引为 256 的颜色也就是白色， （注意：uint8 范围是 0~255 ，而 gray （ 256 ）的索引是 1:256 ，当然这些我们只要了解就可以了，编程并不会用到，因为这些对应的细节 Matlab 已经帮我们做好了） 如果这样调用imshow(I,[ ])， 此时矩阵中的数和颜色表就是线性映射，为什么会这样，我解释一下，我们看这种调用方式和 imagesc(I,[1 64]) 很相似，其实原理是一样的，第二个参数是一个向量， 这个向量指定了矩阵中映射到颜色表的数的范围，也就是显示范围（ Matlab 里叫做 display range ） 前面已经介绍了， Matlab 中 imshow 的 help 中说如果采用 imshow(I,[low high])调用 imshow 的话而且你用 [ ] 代替 [low high]那么 imshow 会使用 [min(I(:)) max(I(:))]作为显示范围，也就说 I 中的最小值会显示成黑色，最大值会显示成白色，这其实就是整个范围的线性映射（没有削顶也没有削底），此时的 imshow(I,[ ]) 函数就相当于 imagesc(I);
imshow 不仅具有 image 的功能也具有 imagesc 的功能，同时体会一下直接映射和线性映射的区别，我们来写一段小程序来测试一下， 程序如下：
clear all;clc;close all;
img = imread('lenna.bmp');
I = rgb2gray(img);
figure(1),image(I); colormap(gray(256));
colorbar,title('show by image in figure1');axis off;
figure(2),imagesc(I);colormap(gray(256));
% We can see that the image showed in figure(2) is a little bright
% than in figure(1).
colorbar,title('show by imagesc in figure2');axis off;
% When we use the imshow, we do not need to set the colormap,
% because the imshow set the colormap to gray(256) automatically.
figure(3),imshow(I),colorbar,title('show by imshow(I) in figure3');
% The effectiveness in figure(3) is the same as in figure(1).
figure(4),imshow(I,[]),colorbar,title('show by imshow(I,[]) in figure3');
% The effectiveness in figure(4) is the same as in figure(2).
显示效果：
figure2 中的图像比 figure1 中的图像要亮一些，而且， figure3 中的显示效果和 figure1 中是一样的， figure4 中的显示效果和 figure2 中是一样的，为什么会这样呢？这是因为 image(I) 和 imshow(I) 是将 I 中的值直接作为 colormap(gray(256)) 中的索引，也就是我所说的直接映射，我这里读到的索引图像矩阵也就是 I 中的数值的范围是 [27,218], 也就是说直接映射显示 I, 只用到的色图 (colormap) 上 [27,218] 范围的颜色，（比如表示白色的索引 255 就没有用到），从右边的colorbar也可以看出来。 但线性映射就不一样了，imagesc(I), 和 imshow(I,[ ]) 采用的就是线性映射，线性映射把 [27,218] 按照线性算法（ Matlab 写的）映射到色图索引 [1,256] 然后再显示出来，这样整个色图的颜色都被用到了，这里也可以认为把 [27,218] 放大到 [1,256] ，这就是 figure2 中显示效果比 figure1 中亮的原因。
小结： 直接映射和线性映射的区别在于映射到色图的数值范围，如果是[min(I(:)) max(I(:))]就是线性映射，如果是0-255 或者 1-64 或者 0-1 就是直接映射。这个数值范围就叫做显示范围（ display range ）。
3、最后再说说 image ， imagesc,imshow 在显示 double 型数据时的用法，
我们做图像处理就会对图像进行运算，使用uint8 型数据精度不高，因为当运算结果超过 255 时会被认为是 255 ，而负数就会被认为是 0 （ 注意在Matlab 中数据默认采用 double 型（ 64 位）进行存储和运算 ）所以，我们读到灰度图像后一般都会将图像转换成double 型（ I = double(I) ）然后再参与运算，运算的结果有正有负，也有小数，正的还可能超过 255 ，比如我经过运算后的得到图像矩阵 I ，假如 I 的范围是 [-187,152] ，当然你也可以用 max(I(:)) 和 min(I(:)) 去获取，这时怎么显示图像呢？ image,imagesc,imshow 都可以用来显示 double 型数据的图像矩阵，主要区别如下：
image： 将double 型数据取整（正数取整就是把小数部分舍掉）然后使用直接映射的方法按照颜色表显示。
imagesc： 这个函数很好，会对数据进行缩放再显示，也就是把显示范围自动设置成[min(I(:)) max(I(:))]，也就是线性
imshow： 这个函数调用方式不同，显示效果也不同，如下：
imshow(I)： 直接调用，因为当图像为double 型时 imshow 函数会把显示范围设置成 [0 , 1], 这样小于 0 的就变成黑色了， 1 的就变成白色了，所以处理不当就会出现全白的情况。
imshow( I/(max(I(:)))： 针对直接调用imshow 函数出现的问题，用 max(I(:) ) 对图像矩阵进行归一化再显示，这样负数
imshow(uint8(I))： 这种方式把I 转化成 uint8 ，负数会被归零，超过 255 的被置为 255 ，而且小数也会被 round( 四舍五 ) ，当参数为 uint8 型时， imshow 函数把显示范围设置成 [0,255] ，这样图像虽然也能显示出来，但与原始数据
imshow(I,[ ])： 这种方式就是把imshow 的显示范围设置成 [min(I(:)) max(I(:))]，也就是线性映射，相当于magesc(I),colormap(gray(256)) 可以将整幅图像的信息显示出来。
imshow与image的用法小节：
相同点：
imshow与image都会产生一个image对象。
区别：
1. imshow的用法：
1）imshow（路径与文件名字符串）
2）imshow（图像矩阵）
3）若当前figure存在坐标轴，imshow会将产生的image对象（即图像对象）显示在当前坐标轴内；
若当前figure不存在坐标轴，imshow会产生一个隐藏的坐标轴，并将产生的 image 对象显示于其中；
4） imshow不会扩展填充图像数据 ， 即不会拉伸图像使其铺满坐标轴 ， 而是改变坐标轴宽高比使其适应图像数据；
2. image的用法
image不会改变坐标轴的大小尺寸 ， 而是扩展填充图像矩阵，使其铺满坐标轴区域 。
image是用来显示附标图像，即显示的图像上有x，y坐标轴的显示，可以看到图像的像素大小。
3. imagesc的用法
imagese函数具有image的功能，所不同的是imagesc函数还自动将输入数据比例化，以全色图的方式显示。
示例
close all; clear; clc size(a) class(a) % displays array a as image % makes the pixels square % set standard 8-bit grayscale colormap
MATLAB中image和imshow都可以用来图像显示，image函数的语法调用常有以下几种格式：image（A）；image（x,y,A）;
其中，image（A）是将矩阵A 作为一个图像显示，A中的每一个元素都被指定一种颜色；image（x,y,A）其中的x,y分别表示显示图像左上角的坐标，其它与image（A）含义相同。当然image还有别的调用格式，这里先不一一介绍了。 对于imshow函数，它的调用格式常见的有：imshow（A,n）；imshow（A,[low high]）;imshow(BW)等。 其中，imshow（A,n）表示利用n个灰度等级来显示一幅灰度图像A，当忽略n时，对于24位显示系统来说n的默认值是256，对于其他系统n默认值是64. imshow（A,[low high]）表示显示灰度图像A，并且指定A的数据范围。A中的数据小于或者等于low的数值被显示为黑的，大于或者等于high的数值被显示为白的，属于区间[low high]的数值自动按照灰度等级进行显示。如果使用空矩阵‘[]’来代替[low high]的话，imshow函数此时自动设置为[min(A) max(A)]就是说，A中最小的值显示为黑色，最大值显示为白色。 imshow（BW）用于显示二进制图像BW，BW中数值为0的像素显示为黑色，数值为1的像素显示为白色。 所以，你说的这三种格式的函数处理图像的时候效果是不一样的，而且要注意处理的图像格式也有区别。
MATLAB图像处理——傅里叶变换
pic=imread(' F:\hua01.jpg '); %读取自己的图片。
pic2=rgb2gray(pic); %RGB彩图转灰度
F=fft2(pic2); %二维fft变换
F=fftshift(F);
F=abs(F);
F =mapminmax(F ,0,1);
T=log(F+1);
T=mapminmax(T,0,1);
subplot(2,2,1); imshow(pic
subplot(2,2,2); imshow(pic2 ); %灰度图
subplot(2,2,3); imshow(F %fft之后图像的频谱分布
subplot(2,2,4); imshow(T %fft之后，数据经过log处理之后，图像的频谱分布
mapminmax介绍：
在利用matlab构建神经网络时，通常需要对数据预处理，一种方式就是将待处理数据归一化到[-1 +1]，如没有指定YMIN和YMAX，默认情形下，将输入矩阵按行进行归一化映射，将输入矩阵中行元素的最小值映射为-1，行元素的最大值映射为+1。
默认：
（1）将矩阵的每行分别进行归一化；
（2）每行的最大值最小值作为每行归一化的xmin和xmax；
（3）将数据归一化到[-1,1].
若要将数据归一化到0到1之间，即y∈[0,1]，使用
b = mapminmax(a,0,1);
————————————————
fftshift介绍：
在matlab中，经过fft变换后，数据的频率范围是从[0,fs]排列的。而一般，我们在画图或者讨论的时候，是从[-fs/2,fs/2]的范围进行分析。因此，需要将经过fft变换后的图像的[fs/2,fs]部分移动到[-fs/2,0]这个范围内。
而fftshift就是完成这个功能。通常，如果想得到所见的中间是0频的图像，经过fft变换后，都要再经过fftshift这个过程。
Matlab fftshift 详解
一.实信号情况
因为实信号以fs为采样速率的信号在 fs/2处混叠，所以实信号fft的结果中前半部分对应[0, fs/2],后半部分对应[ -fs/2, 0]
1）实信号fft的结果前半部分对应[0, fs/2]是正频率的结果,后半部分对应[ -fs/2, 0]是负频率的结果。大于fs/2的部分的频谱实际上是实信号的 负频率加fs的结果 。故要得到正确的结果，只需将 视在频率减去fs即可得到频谱对应的真实负频率
2）如果要让实信号fft的结果与[-fs/2, fs/2]对应，则要fft后fftshift一下即可，fftshift的操作是将fft结果以fs/2为中心左右互换
3）如果实信号fft的绘图频率f从[-fs/2, fs/2]，并且没有fftshift，则fft正频谱对应f在[0, fs/2]的结果将混叠到(f - fs/2)的位置;
fft负频谱对应f在[-fs/2, 0]的结果混叠到 f + fs - fs/2 的位置，注意这里f为负值，也就是说此种情况下fft负频谱对应的 视在频率减去fs/2即可得到频谱对应的真实负频率
二.复信号情况
1）复信号没有负频率，以fs为采样速率的信号，fft的频谱结果是从[0,fs]的。
2）在 f> fs/2时，对复信号的fft结果进行fftshift会产生频率混叠(将下面的示例2中的频率从f=15改为f=85可以验证f=85的谱线在fftshift后跑到 f= -15 = 85 - fs = 85 - 100的位置了)，所以复信号也一般要求 f <= fs/2
3）在对雷达的慢时间维(复信号)进行fft后，由于要用doppler= ((0:LFFT-1)/LFFT - 0.5)*PRF; 计算多普勒频率，所以对该慢时间信号fft后要fftshift下，以便和正确的频率单元相对应。注意多普勒频率fd < = PRF/2 时才测的准！
fftshift
作用：将零频点移到频谱的中间
用法：
Y=fftshift(X)
描述：fftshift移动零频点到频谱中间，重新排列fft,fft2和fftn的输出结果。将零频点放到频谱的中间对于观察傅立叶变换是有用的。
示例1 -实信号的情况：
clf;
fs=100;N=256; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序列 x=0.5*sin(2*pi*15*t)+2*sin(2*pi*40*t); %信号
y1=fft(x,N); %对信号进行快速Fourier变换 y2=fftshift(y1);
mag1=abs(y1); %求得Fourier变换后的振幅 mag2=abs(y2);
f1=n*fs/N; %频率序列 f2=n*fs/N-fs/2;
subplot(3,1,1),plot(f1,mag1,'r'); %绘出随频率变化的振幅 xlabel('频率/Hz');
subplot(3,1,2),plot(f2,mag1,'b'); %绘出随频率变化的振幅 xlabel('频率/Hz');
subplot(3,1,3),plot(f2 , mag2 ,'c'); %绘出随频率变化的振幅 xlabel('频率/Hz');
结论：
1）如果期望绘制的幅频图的频率范围为0~fs，则无需运行fftshift变换，正频率对应在[0, fs/2],
大于fs/2的频谱的频率值为 对应[-fs/2 , 0 ]负频率f + fs，注意f是负频率，是个负数。如图1。
2）如果期望绘制的幅频图的频率范围为-fs/2~fs/2，则需要运行fftshift变换，如图3；
如果不变换，图示的响应频点会发生变换，如图2，分析见顶端。
示例2 -复信号的情况：
close all;
clear;
clf;
fs=100;N=256; %采样频率和数据点数 n=0:N-1;t=n/fs; %时间序列 x=0.5*exp(j*2*pi*15*t)+2*exp(j*2*pi*40*t); %信号
y1=fft(x,N); %对信号进行快速Fourier变换 y2=fftshift(y1);
mag1=abs(y1); %求得Fourier变换后的振幅 mag2=abs(y2);
f1=n*fs/N; %频率序列 f2=n*fs/N-fs/2;
subplot(3,1,1),plot(f1,mag1,'r'); %绘出随频率变化的振幅 xlabel('频率/Hz');
subplot(3,1,2),plot(f2,mag1,'b'); %绘出随频率变化的振幅 xlabel('频率/Hz');
subplot(3,1,3),plot(f2,mag2,'c'); %绘出随频率变化的振幅 xlabel('频率/Hz');
我们知道Fourier分析是信号处理里很重要的技术，matlab提供了强大的信号处理能力，但是有一些细节部分需要我们注意。
记信号f(t)的起始时间为t_start, 终止时间为t_end, 采样周期为t_s, 可以计算信号的持续时间Duration为 t_end – t_start, 信号离散化造成的采样点数 N = Duration/t_s + 1;
根据Fourier分析的相关结论，我们知道时域的采样将会造成频域的周期化，该周期为采样频率f_s（著名的香农采样定理基于此）.
于是， 经过matlab的fft函数处理后，得到数据的横坐标为0:f_s/(N-1):f_s。相关代码如下所示：
%matlab fft 测试代码
t_s = 0.01; %把 t_s 改更小试试？t_s就是采样时间间隔
t_start = 0.5; t_end = 5;
subplot(3,1,1);
Duration = t_end - t_start;
f_s = 1/t_s;
subplot(3,1,2);
subplot(3,1,3);
也就是说，如果我们不使用fftshift，其变换后的横坐标为0:f_s/(N-1):f_s， 如果使用fftshift命令，0频率分量将会移到坐标中心，这也正是matlab中帮助中心给出的意思：对fft的坐标进行了处理。实际上由于频谱的周期性，我们这样做是合理的，可以接受的。
请读者特别要注意横坐标的差别。另外，根据函数的特性，频谱应当只有在15Hz，40Hz出现峰值，但是fft变换后在60Hz，及85Hz处同样出现了峰值，应当可以从fft的计算过程中得到相应的解释（采样率不够高引起的）。
事实上，如果我们用15Hz，60Hz来测试fft变换，也即是 y = 0.5*sin(2*pi*15*t)+2*sin(2*pi*60
t_s=0.001时的结果：
这种现象提醒我们，频率在f_s以内，即 0<f<f_s，f 以及 f_s – f 都有可能是测试信号的频率谱，这就给我们带来了歧义。并且从第三个子图也可以看出，这时候的fftshift会给我们带来错误的引导，也就是说，如果我们试图采样fft或者fftshift来分析信号的频率谱显得不那么靠谱了，matlab的fft谱线与信号的实际频率并不是一一对应的映射关系。这当然不是我们所期望看到的结果，所以实际分析信号时，有关这个问题需要额外的注意。
实际上，这也就间接地证明了Nyquist采样定理的合理性：采样频率要高于截止频率的两倍，上面的处理中我们所使用的采样频率为100Hz，于是当截止频率超过50Hz时，就会出现混叠效应，特殊情况就如上图所示：完全一样。于是，这也就告诉我们若要正确的显示频谱，需要仔细地考量采样频率与截止频率的关系，若太小，则有可能出现混叠，若太大，则计算代价过高。
