function [VG, A, PPG] = colorgrad(f, T)
if (ndims(f)~=3) || (size(f,3)~=3)
    error(\'Input image must be RGB\');
end
sh = fspecial(\'sobel\');
sv = sh\';
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)));
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;
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);
 
if nargin ==2
    VG = (VG>T).*VG;
    PPG = (PPG>T).*PPG;
end

close all;clear all;clc;
f = imread(\'wo.jpg\');
[VG,A,PPG] = colorgrad(f);
ppg = im2uint8(PPG);
ppgf = 255 - ppg;
[M,N] = size(ppgf);
T=150;
ppgf1 = zeros(M,N);
for ii = 1:M
    for jj = 1:N
        if ppgf(ii,jj)<T %边缘区域
%             ppgf1(ii,jj)=0;
            ppgf1(ii,jj)=(T-ppgf1(ii,jj))/5;
        else %平滑区域
            ppgf1(ii,jj)=235/(255-T)*(ppgf(ii,jj)-T);
        end
    end
end
ppgf1 = uint8(ppgf1);
figure;
subplot(221);imshow(ppgf);
subplot(222);imshow(ppgf1);
subplot(223);imhist(ppgf);
subplot(224);imhist(ppgf1);
 
figure;imshow(ppgf1);