ivan@77: function [pdf,val] = genPDF(imSize,p,pctg,distType,radius,disp)
ivan@77: 
ivan@77: %[pdf,val] = genPDF(imSize,p,pctg [,distType,radius,disp])
ivan@77: %
ivan@77: %	generates a pdf for a 1d or 2d random sampling pattern
ivan@77: %	with polynomial variable density sampling
ivan@77: %
ivan@77: %	Input:
ivan@77: %		imSize - size of matrix or vector
ivan@77: %		p - power of polynomial
ivan@77: %		pctg - partial sampling factor e.g. 0.5 for half
ivan@77: %		distType - 1 or 2 for L1 or L2 distance measure
ivan@77: %		radius - radius of fully sampled center
ivan@77: %		disp - display output
ivan@77: %
ivan@77: %	Output:
ivan@77: %		pdf - the pdf
ivan@77: %		val - min sampling density
ivan@77: %
ivan@77: % 
ivan@77: %	Example:
ivan@77: %	[pdf,val] = genPDF([128,128],2,0.5,2,0,1);
ivan@77: %
ivan@77: %	(c) Michael Lustig 2007
ivan@77: 
ivan@77: %   This file is used with the kind permission of Michael Lustig
ivan@77: %   (mlustig@stanford.edu), and originally appeared in the
ivan@77: %   SparseMRI toolbox, http://www.stanford.edu/~mlustig/ .
ivan@77: %
ivan@77: % $Id: genPDF.m 1040 2008-06-26 20:29:02Z ewout78 $
ivan@77: 
ivan@77: 
ivan@77: if nargin < 4
ivan@77: 	distType = 2;
ivan@77: end
ivan@77: 
ivan@77: if nargin < 5
ivan@77: 	radius = 0;
ivan@77: end
ivan@77: 
ivan@77: if nargin < 6
ivan@77: 	disp = 0;
ivan@77: end
ivan@77: 
ivan@77: 
ivan@77: minval=0;
ivan@77: maxval=1;
ivan@77: val = 0.5;
ivan@77: 
ivan@77: if length(imSize)==1
ivan@77: 	imSize = [imSize,1];
ivan@77: end
ivan@77: 
ivan@77: sx = imSize(1);
ivan@77: sy = imSize(2);
ivan@77: PCTG = floor(pctg*sx*sy);
ivan@77: 
ivan@77: 
ivan@77: if sum(imSize==1)==0  % 2D
ivan@77: 	[x,y] = meshgrid(linspace(-1,1,sy),linspace(-1,1,sx));
ivan@77: 	switch distType
ivan@77: 		case 1
ivan@77: 			r = max(abs(x),abs(y));
ivan@77: 		otherwise
ivan@77: 			r = sqrt(x.^2+y.^2);
ivan@77: 			r = r/max(abs(r(:)));			
ivan@77: 	end
ivan@77: 
ivan@77: else %1d
ivan@77: 	r = abs(linspace(-1,1,max(sx,sy)));
ivan@77: end
ivan@77: 
ivan@77: 
ivan@77: 
ivan@77: 
ivan@77: idx = find(r<radius);
ivan@77: 
ivan@77: pdf = (1-r).^p; pdf(idx) = 1;
ivan@77: if floor(sum(pdf(:))) > PCTG
ivan@77: 	error('infeasible without undersampling dc, increase p');
ivan@77: end
ivan@77: 
ivan@77: % begin bisection
ivan@77: while(1)
ivan@77: 	val = minval/2 + maxval/2;
ivan@77: 	pdf = (1-r).^p + val; pdf(find(pdf>1)) = 1; pdf(idx)=1;
ivan@77: 	N = floor(sum(pdf(:)));
ivan@77: 	if N > PCTG % infeasible
ivan@77: 		maxval=val;
ivan@77: 	end
ivan@77: 	if N < PCTG % feasible, but not optimal
ivan@77: 		minval=val;
ivan@77: 	end
ivan@77: 	if N==PCTG % optimal
ivan@77: 		break;
ivan@77: 	end
ivan@77: end
ivan@77: 
ivan@77: if disp
ivan@77: 	figure, 
ivan@77: 	subplot(211), imshow(pdf)
ivan@77: 	if sum(imSize==1)==0
ivan@77: 		subplot(212), plot(pdf(end/2+1,:));
ivan@77: 	else
ivan@77: 		subplot(212), plot(pdf);
ivan@77: 	end
ivan@77: end
ivan@77: 
ivan@77: 
ivan@77: 
ivan@77: 
ivan@77: 
ivan@77: