function [VG, G_ver, G_hor, A, PPG]= colorgrad(f, T)
%COLORGRAD Computes the vector gradient of an RGB image.
% [VG, VA, PPG] = COLORGRAD(F, T) computes the vector gradient, VG,
% and corresponding angle array, VA, (in radians) of RGB image
% F. It also computes PPG, the per-plane composite gradient
% obtained by summing the 2-D gradients of the individual color
% planes. Input T is a threshold in the range [0, 1]. If it is
% included in the argument list, the values of VG and PPG are
% thresholded by letting VG(x,y) = 0 for values <= T and VG(x,y) =
% VG(x,y) otherwise. Similar comments apply to PPG. If T is not
% included in the argument list then T is set to 0. Both output
% gradients are scaled to the range [0, 1].
% Copyright 2002-2004 R. C. Gonzalez, R. E. Woods, & S. L. Eddins
% Digital Image Processing Using MATLAB, Prentice-Hall, 2004
% $Revision: 1.6 $ $Date: 2003/11/21 14:27:21 $
if (ndims(f) ~= 3) | (size(f, 3) ~= 3)
error('Input image must be RGB.');
end
% Compute the x and y derivatives of the three component images
% using Sobel operators.
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');
% Compute the parameters of the vector gradient.
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));
ind1 = find(G1<0);
G1(ind1) = 0;
% Now repeat for angle + pi/2. Then select the maximum at each point.
A = A + pi/2;
G2 = 0.5*((gxx + gyy) + (gxx - gyy).*cos(2*A) + 2*gxy.*sin(2*A));
ind2 = find(G2<0);
G2(ind2) = 0;
G1 = G1.^0.5;
G2 = G2.^0.5;
% Form VG by picking the maximum at each (x,y) and then scale
% to the range [0, 1].
G_ver = mat2gray(G1);
G_hor = mat2gray(G2);
VG = mat2gray(max(G1, G2));
% Compute the per-plane gradients.
RG = sqrt(Rx.^2 + Ry.^2);
GG = sqrt(Gx.^2 + Gy.^2);
BG = sqrt(Bx.^2 + By.^2);
% Form the composite by adding the individual results and
% scale to [0, 1].
PPG = mat2gray(RG + GG + BG);
% Threshold the result.
if nargin == 2
VG = (VG > T).*VG;
PPG = (PPG > T).*PPG;
end