function F = features(image)
% compute a table of feature images
%   1 -> x gradient, no prefiltering
%   2 -> y gradient, no prefiltering
%   3 -> x gradient, gaussian prefiltered
%   4 -> y gradient, gaussian prefiltered
%   5 -> gradient length
%   6 -> edge detector 
%   7 -> basic texture detector

F = [-1 1];
avg = replicateconv(image, [1 1;1 1]/4);
dx = replicateconv(image, F);
dy = replicateconv(image, F');
fimage = replicateConv(image, fspecial('gaussian', 10, 2));
fdx = replicateconv(fimage, F);
fdy = replicateconv(fimage, F');
%figure;
%imshow(fdx/100+0.5);
%figure;
%imshow(dx/256+0.5);

grad = dx.^2 + dy.^2;
%figure
%imshow(grad / 2500);


f = fspecial('gaussian', 10, 1);
A = replicateConv(fdx.^2, f);
B = replicateConv(fdy.^2, f);
C = replicateConv(fdx.*fdy, f);

edges = A.*B-C.^2 + 0.04 * (A + B);
%figure;
%imshow(edges/1000);

himage = ((image)./(fimage+eps)-1).^2;

himage = filter2(f, himage);
himage = filter2(f, himage);
%himage = min(himage, 0.05);
%figure;
%imshow(himage*5);

F = zeros([size(image) 7]);

F(:,:, 1) = dx/256+0.5;
F(:,:, 2) = dy/256+0.5;
F(:,:, 3) = smoothstep(8,fdx/256+0.5);
F(:,:, 4) = smoothstep(4,fdy/256+0.5);
F(:,:, 5) = grad/256/256;
F(:,:, 6) = edges/1000;
F(:,:, 7) = (himage*3).^(0.5);

F = max(0, min(1, F));