function [map,I,Rg,By,hue,saturation,MAD] = skinfilt(red,green,blue)
%
%SKIN_FILT  takes the RGB input matrices of a color image and outputs the
%B&W matrix image with skin regions marked as white and background as black.
%usage: [map,I,Rg,By,hue,saturation,MAD]= skinfilt(red,green,blue)
%

%Get image dimensions
%
[row,column] = size(red);

red = red*255;
blue = blue*255;
green = green*255;

%Account for zero response of camera. Subtract the smallest
%pixel value >10 pixels from any edge from the rgb matrices
%
red_zero_resp = min(min(red(10:row-10,10:column-10)));
blue_zero_resp = min(min(blue(10:row-10,10:column-10)));
green_zero_resp = min(min(green(10:row-10,10:column-10)));

zero_resp = min([red_zero_resp, blue_zero_resp, green_zero_resp]);

red = red-zero_resp;
green = green-zero_resp;
blue = blue-zero_resp;

%Transform the RGB values into log-opponent values I, Rg, By
%The formula for conversion to log-opponent is L(x) = 105*log10(x+1+n)
%I = L(G)
%Rg = L(R)-L(G)
%By = L(B)-(L(G)+L(R))/2
%
I = (logopp(red)+logopp(blue)+logopp(green))/3;
Rg = logopp(red)-I;
By = logopp(blue)-(I+logopp(red))/2;

SCALE = (row+column)/320;
SCALE = round(SCALE);

if SCALE == 0
	SCALE = 1;
end

Rg = medfilt2(Rg,[2*2*SCALE 2*2*SCALE]);

By = medfilt2(By,[2*2*SCALE 2*2*SCALE]);

%Compute texture amplitude
%
I_filt = medfilt2(I,[2*4*SCALE 2*4*SCALE]);
MAD = I-I_filt;
MAD = abs(MAD);
MAD = medfilt2(MAD,[2*6*SCALE 2*6*SCALE]);

hue = (atan2(Rg,By)*(180/pi));
saturation = sqrt(Rg.^2 + By.^2);

map = zeros(size(red));

%Detect skin texture regions
%
for y=1:row,
	for x=1:column,

		if (MAD(y,x)<4.5 & 120<hue(y,x) & hue(y,x)<160 & 10<saturation(y,x) & saturation(y,x)<60),
			map(y,x) = 1;
		end
		
		if (MAD(y,x)<4.5 & 150<hue(y,x) & hue(y,x)<180 & saturation(y,x)>20 & saturation(y,x)<80),
			map(y,x) = 1;
		end

	end
end

%Expand skin regions
%
%se = ones(10*SCALE);
%map = dilate(map,se,5);
map = bwmorph(map,'dilate',5);


%Shrink regions
%
for y=1:row,
	for x=1:column,
		if (map(y,x)==1 & 110<=hue(y,x) & hue(y,x)<=180 & 0<=saturation(y,x) & saturation(y,x)<=130)
			map(y,x)=1;
		else
			map(y,x)=0;
		end
	end
end

end