view general/arrutils/zipmap.m @ 61:eff6bddf82e3 tip

Finally implemented perceptual brightness thing.

function varargout=zipmap(spec,f,varargin)
% zipmap - Map and zip function over multidimensional arrays
%
% ZIPMAP allows a certain simple functions to be applied to multidimensional
% array arguments, manipulating the dimensions to combine elementwise and
% outer-product-like behaviour. (The name comes from functional programming
% terminology: functions/operators like .* and + which apply elementwise to
% vector arguments are said to zip their arguments, while outer products are
% like mapping a function over a list, or a list of functions over another list.)
%
% Let f be a vectorised function which knows how to zip array arguments. Eg,
% f(x,y)=x+y. If x and y are vectors, they must be the same size, and
% the function is applied to corresponding pairs of elements, resulting in
% a vector of the same size. With ZIPMAP, you can build functions, that, eg,
% add EVERY pair of elements (not just corresponding pairs) resulting in a
% 2D array: zipmap(0,@(x,y)x+y,[1:4]',[1:6]') will do this
%
% The domain of a multidimensional array can be described a list containing the
% size of each dimension, as returned by the SIZE function (or SIZE1 which strips
% trailing 1s). The arguments to f must have a certain domain structure:
%
% 	dom arg1 = [p1 zipdom]
% 	dom arg2 = [p2 zipdom]
%	         :
%
% where p1, p2, etc are the domains of the elementary function of which f
% is a vectorisation, eg the elementary domain of .* can be scalars and
% hence p1=p2=[], but the two elementary domains of CROSS are both [3], because it 
% it operates on pairs of 3-D vectors. The zipdom part must be the same for
% all arguments, otherwise, we cannot extract matching argument tuples.
% If this is the case, the f must return a result:
%
% 	dom f(...) : [rdom zipdom]
%
% where rdom is the domain of the elementary function, eg [3] for the vector
% valued CROSS product, but [] for the scalar valued DOT product. Note that
% these domains can be vectors of any length from zero up, eg the determinant
% of a square matrix has a domain of [n n].
% 
% ZIPDIM allows the domains of the arguments and result to be generalised as follows:
%
% 	dom arg1 : [p1 m1 zipdom]
% 	dom arg2 : [p2 m2 zipdom]
%	         :
%
%	dom zipdom(...) : [rdom m1 m2 ... zipdom ]
%
% All the mapped dimensions m1, m2 etc go into a sort of outer product where 
% every combination of values is evaluated. In addition, the order of the
% m1, m2 etc can be permuted in the result at essentially no extra computational cost.
%
% Usage: r=zipmap(spec,f,<varargs>)
% 
% spec defines how the domain of each of the arguments is to be interpreted, ie
% how is is split into elementary, mappable, and zippable domains.
%
% If spec is a number z, then the last z dimensions of each argument are marked for
% zipping, while the rest are for mapping. The elementary domains are assumed to
% be scalar and hence p1, p2 etc =[]. So, zipmap(0,...) does the complete outer
% product type of evaluation.
%
% If spec is structure, the field 'zipdims' determines how many dimensions are zipped.
% The field 'pdims' is a the vector [length(p1) length(p2) ..], ie the dimensionality
% of each elementary domain: 0 for scalars, 1 for vectors, 2 for matrices etc.
% The optional field 'mapord' determines the reordering of the mapped domains, eg
% if mapord=[2 1 3], then the result domain is [m2 m1 m3] instead of [m1 m2 m3].
%
% If spec is a cell array, it is interpreted as {zipdims pdims mapord}.
%
% Extra options
%
% If spec is a structure, it can contain further options: if spec.shift=1,
% then the domains of all arguments are pre-shifted to remove singleton
% dimensions before anything else is done. Amongst other effects, this means
% that row vectors can be used in place of column vectors.

	% unpack args from spec
	if iscell(spec)
		[zipdims,pdims,mapord]=getargs(spec,{0 [] []});
	elseif isstruct(spec)
		zipdims=getparam(spec,'zipdims',0);
		pdims=getparam(spec,'pdims',[]);
		mapord=getparam(spec,'mapord',[]);
		if getparam(spec,'shift',0),
			varargin=cellmap(@shiftdim,varargin);
		end
	elseif isscalar(spec)
		zipdims=spec; pdims=[]; mapord=[];
	end


	argdoms=cellmap(@size1,varargin);
	if isempty(pdims), pdims=zeros(size(varargin)); end
	for i=1:length(varargin)
		[pdoms{i},mapdoms{i},zipdoms{i}]=split(argdoms{i},pdims(i),zipdims);
	end


	% check all zipdoms match and use first
	zipdom=zipdoms{1};
	for i=2:length(zipdoms)
		if ~all(zipdoms{i}==zipdom), error('Zip domains do not match'); end
	end


	if isempty(mapord), mapord=1:length(varargin); end
	outerdom=cat(2,mapdoms{mapord});
	outer1=cellmap(@to1s,mapdoms);
	for i=1:length(varargin)
		% need to reshape each arg and then repmat up to correct size
		newsizes=outer1; newsizes{i}=mapdoms{i};
		newsize=[pdoms{i} cat(2,newsizes{mapord}) zipdom];
		varargin{i}=reshape(varargin{i},tosize([pdoms{i} cat(2,newsizes{mapord}) zipdom]));
		newsizes=mapdoms; newsizes{i}=outer1{i};
		newsize=[to1s(pdoms{i}) cat(2,newsizes{mapord}) to1s(zipdom)];
		varargin{i}=repmat(varargin{i},tosize(newsize));
	%	size(varargin{i})
	end

	[varargout{1:max(1,nargout)}]=feval(f,varargin{:});
end

function [a,b,c]=split(x,n,m)
	a=x(1:n);
	b=x(n+1:end-m);
	c=x(end-m+1:end);
end

% getargs - get values of optional parameters with defaults
%
% getargs :: {I:[N]->(A(I) | nil)}, {I:[N]->A(I)} -> A(1), ..., A(N).
function varargout=getargs(args,defs)
	nout=max(length(args),length(defs));
	varargout=defs;
	for i=1:length(args)
		if ~isempty(args{i}), varargout(i)=args(i); end
	end
end