clear;

X = imread(‘testimages/sample1.png‘); % sample1.png is a gray-scale CG generated face image

[H, W] = size(X);
S1 = cat(3,X,X,X); % used for displaying final result (Geometric mean)
S2 = cat(3,X,X,X); % used for displaying final result (Arithmetic mean)
X = double(X); % convert data type to double
I1 = cvtIntegralImage(X);       % calculate integral image
P1 = cvtIntegralImage(X.^2);    % calculate integral image of squared pixel value
I2 = cvtIntegralImage45(X);     % calculate 45 degrees integral image
P2 = cvtIntegralImage45(X.^2);  % calculate 45 degrees integral image of squared pixel value

nR = 3; % filter size parameter
nTH = 0.55; % threshold for finding local peaks

P = zeros(H,W,4);  % variable to store separability map
P(:,:,1:2) = cvtCombSimpRectFilter(I1,P1,nR);   % apply vertical and horizontal rectangular filters
P(:,:,3:4) = cvtCombSimpRectFilter45(I2,P2,nR); % apply diagonal left and right filters
P(P<0) = 0;
finalMap1 = prod(P(:,:,:),3).^(1/4.0);
finalMap2 = mean(P(:,:,:),3);

figure(10);clf;

for i=1:6
    subplot(2,4,i);
    if (i < 5)
        imagesc(P(:,:,i));
        axis equal tight;
        title([‘separability map #‘ num2str(i)]);
    elseif (i==5)
        imagesc(finalMap1);
        axis equal tight;
        title(‘Geometric mean‘);
    elseif (i==6)
        imagesc(finalMap2);
        axis equal tight;
        title(‘Arithmetic mean‘);
    end
end

% find local peaks (Geometric mean)
PL1 = cvtFindLocalPeakX(finalMap1,1,nTH);
% draw circle and cross at each local peak with radius of the filter (nR)
for H=1:size(PL1,2)
    S1 = cvtDrawCircle(S1, PL1(2,H),PL1(1,H),nR,[255,0,0],20);
    S1 = cvtDrawCross(S1,PL1(2,H),PL1(1,H),nR,[255,255,255]);
end
subplot(2,4,7); 
image(S1); % display original
title({[‘Local peaks > ‘ num2str(nTH)]; ‘Geometric mean‘});
axis equal tight;

% find local peaks (Arithmetic mean)
PL2 = cvtFindLocalPeakX(finalMap2,1,nTH);
% draw circle and cross at each local peak with radius of the filter (nR)
for H=1:size(PL2,2)
    S2 = cvtDrawCircle(S2, PL2(2,H),PL2(1,H),nR,[255,0,0],20);
    S2 = cvtDrawCross(S2,PL2(2,H),PL2(1,H),nR,[255,255,255]);
end
subplot(2,4,8);
image(S2); % display original with marks for the local peak
axis equal tight;
title({[‘Local peaks > ‘ num2str(nTH)]; ‘Arithmetic mean‘});