function [D,xm,ym,imin,imax,jmin,jmax] = f(im,S,U,Xmean,n,imin,imax,jmin,jmax)
% Distance-from-Feature-Space Map
%
% Syntax: [D,rowm,colm,imin,imax,jmin,jmax] = dffs_map(IM,Size,U,Xmean,n,
%						       imin,imax,jmin,jmax)
%
% input:
% IM		input image matrix
% Size  	2-by-1 vector denoting the dimensions of the eigenvectors in U
% U		matrix of eigenvectors (one per column)
% Xmean 	column vector of mean
% n     	index vector denoting the eigenvectors in the feature-space 
% imin,imax	(optional) arguments denoting image region to process
% jmin,jmax
%
% output:
% D     	output DFFS map
% rowm,colm	row,col indices of global minimum of D
% imin,imax	region overwhich DFFS values have been computed 
% jmin,jmax

% ----------------------------------------------------------------------
% Author: baback@media.mit.edu
% Date  : 5/9/93

% determine dimensions of input image and eigenvectors
Nrow = S(1); Ncol = S(2);
if (Nrow*Ncol) ~= length(Xmean)
	disp('ERROR: dimensions in Size do not agree with U,Xmean')
	return;
end
[M,N] = size(im);

% definedefault processing window region
rowl = fix(Nrow/2);
rowr = Nrow - rowl - 1;
coll = fix(Ncol/2);
colr = Ncol - coll - 1;

% define processing window limits
if nargin<6
  imin = rowl + 1; 
  imax = M - rowr;
  jmin = coll + 1;
  jmax = N - colr;
end

% setup output map and trim the U matrix
D = zeros(size(im));
U = U(:,n);

% compute DFFS at each pixel within processing window
for i=imin:imax
for j=jmin:jmax
 X = im(i-rowl:i+rowr,j-coll:j+colr);
 Xn = X(:) - Xmean;
 C = U'*Xn;
 E = Xn'*Xn;
 D(i,j) = E - C'*C;
end;
end;

% determine coordinates of global minimum
dmin = min(min(D(imin:imax,jmin:jmax)));
[xm,ym] = find(D==dmin);