数字图像处理之素描
这是我的随机信号作业，直接从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);
