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