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