数字图像处理之素描

　　这是我的随机信号作业，直接从word复制过来的，凑合凑合看吧。

　　涉及到的原理主要就是一个sobel边缘提取：

　　Sobel算法：

　　令 r，g，b 是 RGB 彩色空间沿 R，G，B 轴的单位向量，可定义向量为

　　数量定义为这些向量的点乘，如下所示：

　　利用该表示法，c(x,y)的最大变化率方向可以由角度给出：

　　（x，y）点在方向上变化率的幅值由下式给出：

　　过程与方法：

　　将Sobel算子的实现划分为五个步骤：

　　（1）计算Gx与Gy与模板每行的乘积。

　　（2）两个3x3矩阵的卷积即将每一行每一列对应相乘然后相加。

　　（3）求得3*3模板运算后的Gx、Gy。

　　（4）求Gx^2 + Gy^2的平方根或者直接对Gx和Gy取绝对值后求和。

　　（5）设置一个阈值，运算后的像素值大于该阈值输出为全1，小于该阈值输出为全0。

　　结果与分析

　　　　　　图3-1　　　　　　　　　　　　　　　　　　　　　图3-2

　　　　　　图3-3 T=180

　　　　　　　　图3-4 T=230

　　分析：

　　图3-1，通过sobel算法，提取出图像的边缘细节，边缘部分以白色显示，这是一个二值图像。

　　图3-2，3-1图像背景为黑色，不易辨别，也不美观，将背景色换为白色，边缘为黑色显示，如图3-3左上角，但是图片的边缘勾勒不是很清晰，

　　于是通过加强边缘的够了，将灰度在一个阈值范围内的灰度进行增强。

　　图3-3下面的两个直方图，反映了对应的上方图的灰度等级，左边的直方图由于左上角图偏白，边缘不够清晰，所以灰度值靠近255的居多；

　　右边的直方图相比之下，灰度值靠近0的有所增加，边缘更加明显。

　　图3-4与3-3相比，由于改变了T的值，提高了对比度，靠近黑色的灰度有所增加。

%白色为255，黑色0
close all;clear all;clc;
f = imread('manhua.jpg');
[VG,A,PPG] = colorgrad(f);
%把图像数据类型转换为无符号八位整型
ppg = im2uint8(PPG);
%由于边缘检测出来的边缘为白色，背景为黑色，
%所以如果想要背景成为白色，边缘成为黑色，需要用255-pgg
ppgf = 255 - ppg;
%[r,c]=size（A）
%将矩阵A的行数返回到第一个输出变量r，将矩阵的列数返回到第二个输出变量c
[M,N] = size(ppgf);T=230;
ppgf1 = zeros(M,N);
for ii = 1:M
    for jj = 1:N
        if ppgf(ii,jj)<T %边缘区域
            ppgf1(ii,jj)=0;
        else%平滑区域
            ppgf1(ii,jj)=235/(255-T)*(ppgf(ii,jj)-T);
        end
    end
end
% unit8 255
ppgf1 = uint8(ppgf1);
figure;
subplot(221);imshow(ppgf);
subplot(222);imshow(ppgf1);
subplot(223);imhist(ppgf);
subplot(224);imhist(ppgf1);
figure;imshow(ppgf1);
figure;
imshow(ppg);

　　colorgrad.m

%梯度边缘提取算法，计算彩色图像的梯度
%VG是RGB向量梯度Fθ(x, y)
%A是以弧度计的角度θ(x, y)
%PPG是由单独彩色平面的2D梯度之和形成的梯度图像
%输出VG和PPG通过clorgrad被归一化到[0, 1]范围内，并且把它们作为阈值进行处理
%ndims是matlab中求一个数组维数的函数，size(f,3)返回f的行数
function [VG, A, PPG] = colorgrad(f, ~)
if (ndims(f)~=3) || (size(f,3)~=3)              
    error('Input image must be RGB');
end
sh = fspecial('sobel');                          %sobel用于边缘提取
sv = sh';                                        %转置
%r,g,b值 x,y方向, 3个色度卷
%线性空间滤波函数；replicate表示图像大小通过复制外边界的值来扩展，
Rx = imfilter(double(f(:,:,1)), sh, 'replicate');
Ry = imfilter(double(f(:,:,1)), sv, 'replicate');
Gx = imfilter(double(f(:,:,2)), sh, 'replicate');
Gy = imfilter(double(f(:,:,2)), sv, 'replicate');
Bx = imfilter(double(f(:,:,3)), sh, 'replicate');
By = imfilter(double(f(:,:,3)), sv, 'replicate');
 
gxx = Rx.^2 + Gx.^2 + Bx.^2;
gyy = Ry.^2 + Gy.^2 + By.^2;
gxy = Rx.*Ry + Gx.*Gy + Bx.*By;
A = 0.5*(atan(2*gxy./(gxx-gyy+eps)));            
%atan计算反正切函数                
%./点除                                          
%如果a、b是矩阵，a./b就是a、b中对应的每个元素相除，得到一个新的矩阵；
%如果a、b是两个数，那么a./b就是普通的除法
%%eps是一个函数，可以返回某一个数N的最小浮点数精度
G1 = 0.5*((gxx+gyy) + (gxx-gyy).*cos(2*A) + 2*gxy.*sin(2*A));  
A = A + pi/2;
G2 = 0.5*((gxx+gyy) + (gxx-gyy).*cos(2*A) + 2*gxy.*sin(2*A));
G1 = G1.^0.5;
G2 = G2.^0.5;
%mat2gray实现图像矩阵的归一化操作
VG = mat2gray(max(G1, G2));            
RG = sqrt(Rx.^2 + Ry.^2);
GG = sqrt(Gx.^2 + Gy.^2);
BG = sqrt(Bx.^2 + By.^2);
PPG = mat2gray(RG + GG + BG);