数字图像处理-图像的点运算

1.点运算概述

　　数字图像处理系统处理流程：预处理--->特征提取--->识别分析.

　　点运算是对图像中的每个像素依次进行同样的灰度变换运算。设输入图像f(x,y)和输出图像g(x,y)在任一点(x,y)的灰度值，点运算表示为：g(x,y) = T(f(x,y))。其中T为采用的点运算算子，表示原始图像和输出图像之间的某种灰度级映射关系。

　　点运算常用于改变图像的灰度范围及分布，点运算因其作用的性质有时也被称为对比度增强、对比度拉伸或灰度变换。(对比度：灰度图像中最大亮度与最小亮度的比值)

2.灰度直方图

　　灰度直方图描述了一幅图像的灰度级统计信息，主要用于图像分割和图像灰度变换等处理过程。

　　从数学较读来说，图像直方图描述的是图像中各个灰度级统计特性，它是关于图像灰度值的函数，用于统计一幅图像中各个灰度级出现的次数或概率(即各个灰度级像素出现的次数或占比)。归一化直方图可以直接反映不同灰度级出现的比率。

　　从图形上来说，灰度直方图是一个二维图，横坐标为图像中各个像素点的灰度级别，纵坐标为各个灰度级别的像素出现的次数或概率。

 　　Matlab生成图像的灰度直方图：

I = imread('pout.tif');
figure;                                          % 创建一个窗口 
imshow(I); title('Source');              % 显示原图，设置标题
figure;
imhist(I);title('Graph');                  % 生成灰度直方图，默认是256个灰度值，区间也是256个(0~255)，此时统计的是在256个区间上像素出现的次数
figure;
>> imhist(I, 64); title('Graph64');  % 将灰度级划分成64个区间  

 　　绘制归一化直方图：

I = imread('pout.tif');
figure;
[M, N] = size(I);
[counts,x] = imhist(I, 32); % 32个区间，返回落在各个区间的像素个数counts, x是各个区间
counts = counts / (M * N); % 归一化
stem(x, counts); % 绘制归一化直方图

 

 3.灰度的线性变换

　　Matlab示例：

I = imread('coins.png');
I = im2double(I); % 转成浮点表示
[M,N] = size(I);

figure(1);
imshow(I);
title('Source Image');

figure(2);
[H,x] = imhist(I, 64);
stem(x, (H / (M *  N)), '.');
title('Source Stem'); % 归一化直方图

% 增加对比度
Fa = 2; Fb = -55;
O =  Fa .* I + Fb / 255;
figure(3);
subplot(2,2,1);
imshow(O);
title('Fa = 2 Fb = -55 Add Contrast');
figure(4);
subplot(2,2,1);
[H,x] = imhist(O, 64);
stem(x, (H/M/N), '.');
title('Fa = 2 Fb = -55 Add Contrast');

% 减少对比度
Fa = 0.5; Fb = -55;
O =  Fa .* I + Fb / 255;
figure(3);
subplot(2,2,2);
imshow(O);
title('Fa = 2 Fb = -55 Sub Contrast');
figure(4);
subplot(2,2,2);
[H,x] = imhist(O, 64);
stem(x, (H/M/N), '.');
title('Fa = 2 Fb = -55 Sub Contrast');

% 增加亮度
Fa = 1; Fb = 55;
O =  Fa .* I + Fb / 255;
figure(3);
subplot(2,2,3);
imshow(O);
title('Fa = 1 Fb = 55 Add Brightness');
figure(4);
subplot(2,2,3);
[H,x] = imhist(O, 64);
stem(x, (H/M/N), '.');
title('Fa = 1 Fb = 55 Add Brightness');

% 反向显示
Fa = -1; Fb = 255;
O =  Fa .* I + Fb / 255;
figure(3);
subplot(2,2,4);
imshow(O);
title('Fa = -1 Fb = 255 Reverse');
figure(4);
subplot(2,2,4);
[H,x] = imhist(O, 64);
stem(x, (H/M/N), '.');
title('Fa = -1 Fb = 255 Sub Reverse');

 

 

 　　从上面的示例可以看出改变图像的对比度是对直方图的缩放与平移，改变图像的亮度则只是平移直方图在横轴上的位置，而反向就是将直方图水平镜像。单纯的线性灰度变换只能在一定程度上解决视觉上整体对比度问题，但是对于图像细节部分的增强较为有限，一般还需要结合非线性变换来解决。

4.灰度对数变换

　　对数变换的一般表达式：

 　　如下图的对数函数曲线图，函数自变量低值时曲线斜率很高，自变量高值时曲线斜率变小。由对数函数曲线可知，这种变换可以增强一幅图像中较暗部分(相当于低值自变量)的细节，从而可以用来扩展被压缩的高值图像中的较暗像素，因此对数变换被广泛运用于频谱图像的显示中。一个典型的运用是傅里叶频谱，其动态范围可能宽达0~10^6.直接显示频谱时，图像显示设备的动态范围往往不能满足要求，从而会丢失大量的暗部细节，而在使用对数变换之后，图像的动态范围被合理的非线性压缩，从而可以清晰的显示。

 　　使用如下数学函数实现对图像I的对数变换：

 　　

　　下面使用Matlab实现对傅里叶频谱图像进行对数变换，从

% 读取图像
I = imread('coins.png');

% 计算图像的傅里叶频谱
F = fft2(im2double(I));
F = fftshift(F);
F = abs(F);

% 对数变换
T = log(F + 1);

subplot(1,2,1);
imshow(F, []);
title('未经变换的傅里叶频谱');

subplot(1,2,2);
imshow(T, []);
title('经对数变换后的频谱图');

 5.伽玛变换

　　伽玛变换又叫指数变换或幂次变换，时另一种常用的灰度非线性变换。

 　　伽马变换的映射关系如下图，在进行变换时通常要将0~255的灰度范围先先变换到0~1的范围，执行完伽马变换后再映射到原来的范围。

　　如下程序实现了gamma系数不同取值时的效果：

I = imread('pout.tif');

% Gamma = 0.75
subplot(1,3,1);
imshow(imadjust(I, [], [], 0.75));
title('Gamma 0.75');

% Gamma = 1.0
subplot(1,3,2);
imshow(imadjust(I, [], [], 1.0));
title('Gamma 1.0');

% Gamma = 1.5
subplot(1,3,3);
imshow(imadjust(I, [], [], 1.5));
title('Gamma 1.5');

 从效果可以看出伽马因子给图像的整体明暗程度带来的变化，以及图像暗部和亮部细节清晰度的影响。当伽马银子取1时，图像没有任何改变。下面生成3幅图像的灰度直方图：

subplot(1,3,1);
imhist(imadjust(I, [], [], 0.75));
title('Gamma=0.75');

subplot(1,3,2);
imhist(imadjust(I, [], [], 1));
title('Gamma=1');

subplot(1,3,3);
imhist(imadjust(I, [], [], 1.5));
title('Gamma=1.5');

 　　观察直方图非零值区间位置的变化，以及这些变化给图像带来的影响。由于伽马变换并不是线性变换，所以它不仅可以改变图像的对比度，还能增强细节，从而带来整体效果的改善。

6.灰度阈值变换

　　灰度阈值变换可以将一幅灰度图像转换成黑白的二值图像，这种过程也叫二值化或阈值化。用户指定一个起到分解限作用的灰度值，如果图像中某像素的灰度值小于该灰度值，则将该像素的灰度值设置为0，否则设置为255，这个起到分解先作用的灰度值称为阈值。通过阈值化可以直接将图像内容划分为我们关心和不关心的2个部分，是图像分割的重要手段之一。

 

 　　Matlab与阈值相关的函数主要2个：im2bw和graythresh

BW = im2bw(I, level)
将灰度图像I转换为二值图像BW，其中level是阈值等级它是一个[0,1]之间的double类型，可以映射到对应的灰度值，未指定时默认是0.5。
例如输入图像I为灰度范围在0~255之间的uint8图像，指定level=0.5则对应的分割阈值为128.

BW = im2bw(X, map, level)
将带有颜色映射表的索引图像X转换成二值图像

BW = im2bw(RGB, level)
将真彩色图像RGB转换成二值图像

如果输入图像不是灰度图像，则im2bw将输入图像转换为灰度图像，然后通过阈值将此灰度图像转换为二进制图像。



level = graythresh(I)
使用Otsu算法(大津算法)计算一个出全局阈值，该阈值可用于将灰度图像转换为二值图像。I为需要计算阈值的输入图像。

具体说明可以查看Matlab帮助文档，使用doc+函数名 命令即可查看。


Matlab2017B版本

为了使图片转换为二值图，利用 imbinarize 函数,其中的参数必须为二维的灰度图：

bw=imbinarize( g ); 或 bw=im2bw( g );

在 matlab2018 中建议用 imbinarize 来将图片转换为二值图，其参数必须为灰度图。

在 matlab2016 中，只有 im2bw 函数，其参数可以是灰度图或 rgb 图。

即在 matlab2016 中可以用 im2bw 直接将 rgb 图像转换为二值图；但是在 matlab2018 中只能先将 rgb 图片转换为灰度图（用 rgb2gray() 函数），再用 imbinarize 转换为二值图。

两个函数都是将图片转换为灰度图，除了参数的类型不同，两个函数都可再传入一个 [0,1] 的阈值，即：

bw = im2bw( g , level ); 或 bw = imbinarize( g , level );

阈值（level）的作用为：用值 1 (白色)替换输入图像中亮度大于 level 比例的所有像素，用值 0 (黑色)替换所有其他像素。

当 im2bw 的 level 省略时，默认为 0.5 。

下面是一个示例：

I = imread('rice.png');
thresh = graythresh(I); % 使用大津算法计算最优阈值
bw1 = im2bw(I, thresh); % 阈值化
bw2 = imbinarize(I, 130/255); % 使用130阈值做阈值化
subplot(1,3,1);imshow(I);title('Source Img');
subplot(1,3,2);imshow(bw1);title('Threshed By Otsu');
subplot(1,3,3);imshow(bw2);title('Threshed By 130');

 7.分段线性变换

　　分段线性变换又很多种，包括灰度拉伸、灰度窗口变换等，本节只分析最常用的灰度拉伸。

 

 

 以下是示例，先定义了2个自定义函数：

% parse_inputs.m
% 分析输入参数个数和有效性
% A         输入图像，RGB图(3D)，灰度图(2D)，索引图(X)
% map       索引图的调色板(颜色表) (:,3)
% [x1, x2]  参数组1，分段线性曲线中2个转折点的横坐标
% [y1, y2]  参数组2，分段线性曲线中2个转折点的纵坐标
function [A, map, x1, x2, y1, y2] = parse_inputs(varargin)

map = [];

% iptchecknargin(low, high, num_inputs, func_name)函数
% 检查nargin中包含的输入遍历个数是否在low和high指定的范围，不在范围则报错。
% 使用方法参考帮助文档
iptchecknargin(3,4,nargin,mfilename);

% iptcheckinput(A, classes, attributes, func_name, var_name, arg_pos)函数
% 检查矩阵A中元素是否属于给的的类型列表，如果有元素不属于给定类型则报错
iptcheckinput(varargin{1},...
            {'uint8', 'uint16', 'int16', 'double'},...
            {'real', 'nonsparse'},mfilename,'I, X or RGB',1);
            
% 根据参数个数返回不同的值
switch nargin
    case 3 % 参数可能是(I, [x1,x2],[y1,y2])或(RGB, [x1,x2],[y1,y2])
    A = varargin{1};
    x1 = varargin{2}(1);
    x2 = varargin{2}(2);
    y1 = varargin{3}(1);
    y2 = varargin{3}(2);
    case 4 % (I, map,[x1,x2],[y1,y2])
    A = varargin{1};
    map = varargin{2};
    x1 = varargin{3}(1);
    x2 = varargin{3}(2);
    y1 = varargin{4}(1);
    y2 = varargin{4}(2);
end

% 检查输入参数的有效性
% 检查RGB数组
if(ndims(A) == 3) && (size(A,3) ~= 3)
    msg = sprintf('%s:真彩色图像应使用M-N-3维度数组', mfilename);
    eid = sprintf('Images:%s:RGB Color Image Must Be M-N-3', mfilename);
    error(eid, '%s', msg);
end

% 检查调色板
if ~isempty(map)
    if (size(map,2) ~= 3) || ndims(map)>2
        msg1 = sprintf('%s:调色板必须是一个矩阵，有3列(R G B)', mfilename);
        error('Error', '%s', msg1);
    elseif (min(map(:))<0) || (max(map(:)>1))
        msg1 = sprintf('%s:调色板中各个分量的强度应在[0,1]范围', mfilename);
        error('Error', '%s', msg1);
    end
end

% int16转uint16
if isa(A, 'int16')
    A = int16touint16(A);
end 

% imgrayscaling.m
% 执行灰度拉伸
% 语法：
%     out = imgrayscaling(I, [x1,x2],[y1,y2]);
%     out = imgrayscaling(X, map, [x1,x2],[y1,y2]);
%     out = imgrayscaling(RGB, [x1,x2],[y1,y2]);
% 执行灰度拉伸，输入图像应为灰度图像，但是如果不是灰度图像函数会自动将其转换为灰度图x1,x2,y1,y2应使用double类型
function out = imgrayscaling(varargin)
% 解析参数
[A, map, x1, x2, y1, y2] = parse_inputs(varargin{:});

% 计算输入图像A中数据类型对应的取值范围
range = getrangefromclass(A);
range = range(2);

% 如果不是灰度图则转换
if ndims(A)==3 % A矩阵为3维，RGB图像
    A = rgb2gray(A);
elseif ~isempty(map) % 索引图像
    A = ind2gray(A, map);
end

% 读取原始图像大小，并初始化输出图像
[M,N] = size(A);
I = im2double(A); % 转double
out = zeros(M,N);

% 灰度拉伸
for m=1:M
    for n=1:N
        if I(m,n)<x1
            out(m,n)=y1*I(m,n)/x1;
        elseif I(m,n)>x2
            out(m,n)=(I(m,n)-x2)*(range-y2)/(range-x2)+y2;
        else
            out(m,n)=(I(m,n)-x1)*(y2-y1)/(x2-x1)+y1;
        end
    end
end

% 将输出图像的格式转换成与输入图像相同
if isa(A,'uint8') 
    out = im2uint8(out);
elseif isa(A,'uint16')
    out=im2uint16(out);
% 其他情况输出double
end

% 输出处理
if nargout==0 % 外面没有接收返回值
    imshow(out);
    return;
end

I = imread('coins.png');
J1=imgrayscaling(I,[0.3 0.7],[0.15 0.85]);
figure,imshow(J1,[]);title('[x1 x2]=[0.3 0.7], [y1 y2]=[0.15 0.85]');
J2=imgrayscaling(I,[0.15,0.85],[0.3,0.7]);
figure,imshow(J2,[]);title('[x1 x2]=[0.15,0.85], [y1 y2]=[0.3,0.7]');
figure(3);
subplot(1,2,1);
imhist(J1);
subplot(1,2,2);
imhist(J2);

 

 8.直方图均衡化

　　直方图均衡化也叫灰度均衡化，是指通过灰度映射使输入图像转换为在每一灰度级上都有近似相同的像素点数的输出图像(即输出的直方图是均匀的)。经过处理的图像中像素将占有尽可能多的灰度级并且分布均匀，因此具有较高的对比度和动态范围。

　　首先考虑灰度级范围为[0,1]连续的情况，此时图像的归一化直方图即为概率密度函数(在数学中，连续型随机变量的概率密度函数是一个描述这个随机变量的输出值在某个确定的取值点附近的可能性的函数；随机变量的取值落在某个区域之内的概率则为概率密度函数在这个区域上的积分。当概率密度函数存在的时候，累积分布函数是概率密度函数的积分。百度百科：https://baike.baidu.com/item/%E6%A6%82%E7%8E%87%E5%AF%86%E5%BA%A6%E5%87%BD%E6%95%B0/5021996?fr=aladdin)：

I=imread('pout.tif');
I=im2double(I);

% 对于对比度较大的图像
I1=2*I-55/255;
subplot(4,4,1);
imshow(I1);
subplot(4,4,2);
imhist(I1);
subplot(4,4,3);
% [J,T]=histeq(I)函数用于直方图均衡化，输入图像I，输出图像J和变换矩阵T
imshow(histeq(I));
subplot(4,4,4);
imhist(histeq(I1));

% 对于对比度较小的图像
I2=0.5*I+55/255;
subplot(4,4,5);
imshow(I2);
subplot(4,4,6);
imhist(I2);
subplot(4,4,7);
imshow(histeq(I2));
subplot(4,4,8);
imhist(histeq(I2));

% 对于线性增加亮度的图像
I3=I+55/255;
subplot(4,4,9);
imshow(I3);
subplot(4,4,10);
imhist(I3);
subplot(4,4,11);
imshow(histeq(I3));
subplot(4,4,12);
imhist(histeq(I3));

% 对于线性减小亮度的图像
I4=I-55/255;
subplot(4,4,13);
imshow(I4);
subplot(4,4,14);
imhist(I4);
subplot(4,4,15);
imshow(histeq(I4));
subplot(4,4,16);
imhist(histeq(I4));

 

 9.直方图规定话(匹配)

 

 　　Matlab示例：

% 实现从图像I分别到图像I1和I2的直方图匹配
I=imread('pout.tif');          % 原图
I1=imread('coins.png');        % 要匹配直方图的图像
I2=imread('circuit.tif');      % 要匹配直方图的图像

% 计算直方图
[hgram1,x]=imhist(I1);
[hgram2,x]=imhist(I2);

% 执行直方图均衡化
J1=histeq(I,hgram1);
J2=histeq(I,hgram2);

% 绘图
subplot(2,3,1);
imshow(I);title('原图');
subplot(2,3,2);
imshow(I1);title('标准图1');
subplot(2,3,3);
imshow(I2);title('标准图2');
subplot(2,3,5);
imshow(J1);title('匹配到1');
subplot(2,3,6);
imshow(J2);title('匹配到2');

% 直方图
figure;
subplot(2,3,1);
imhist(I);title('原图');
subplot(2,3,2);
imhist(I1);title('标准图1');
subplot(2,3,3);
imhist(I2);title('标准图2');
subplot(2,3,5);
imhist(J1);title('规定化到1');
subplot(2,3,6);
imhist(J2);title('规定化到2');