view util/sparco utils/genPDF.m @ 221:c1efdd5d6250 luisf_dev

added ompbox_fast as sparse approximation solver, which calls the functions in ompbox using the fastest implementation

function [pdf,val] = genPDF(imSize,p,pctg,distType,radius,disp)

%[pdf,val] = genPDF(imSize,p,pctg [,distType,radius,disp])
%
%	generates a pdf for a 1d or 2d random sampling pattern
%	with polynomial variable density sampling
%
%	Input:
%		imSize - size of matrix or vector
%		p - power of polynomial
%		pctg - partial sampling factor e.g. 0.5 for half
%		distType - 1 or 2 for L1 or L2 distance measure
%		radius - radius of fully sampled center
%		disp - display output
%
%	Output:
%		pdf - the pdf
%		val - min sampling density
%
% 
%	Example:
%	[pdf,val] = genPDF([128,128],2,0.5,2,0,1);
%
%	(c) Michael Lustig 2007

%   This file is used with the kind permission of Michael Lustig
%   ([e-mail]), and originally appeared in the
%   SparseMRI toolbox, http://www.stanford.edu/~mlustig/ .
%
% $Id: genPDF.m 1040 2008-06-26 20:29:02Z ewout78 $


if nargin < 4
	distType = 2;
end

if nargin < 5
	radius = 0;
end

if nargin < 6
	disp = 0;
end


minval=0;
maxval=1;
val = 0.5;

if length(imSize)==1
	imSize = [imSize,1];
end

sx = imSize(1);
sy = imSize(2);
PCTG = floor(pctg*sx*sy);


if sum(imSize==1)==0  % 2D
	[x,y] = meshgrid(linspace(-1,1,sy),linspace(-1,1,sx));
	switch distType
		case 1
			r = max(abs(x),abs(y));
		otherwise
			r = sqrt(x.^2+y.^2);
			r = r/max(abs(r(:)));			
	end

else %1d
	r = abs(linspace(-1,1,max(sx,sy)));
end




idx = find(r<radius);

pdf = (1-r).^p; pdf(idx) = 1;
if floor(sum(pdf(:))) > PCTG
	error('infeasible without undersampling dc, increase p');
end

% begin bisection
while(1)
	val = minval/2 + maxval/2;
	pdf = (1-r).^p + val; pdf(find(pdf>1)) = 1; pdf(idx)=1;
	N = floor(sum(pdf(:)));
	if N > PCTG % infeasible
		maxval=val;
	end
	if N < PCTG % feasible, but not optimal
		minval=val;
	end
	if N==PCTG % optimal
		break;
	end
end

if disp
	figure, 
	subplot(211), imshow(pdf)
	if sum(imSize==1)==0
		subplot(212), plot(pdf(end/2+1,:));
	else
		subplot(212), plot(pdf);
	end
end