comparison Problems/private/genPDF.m @ 61:42fcbcfca132

(none)
author idamnjanovic
date Tue, 15 Mar 2011 12:21:31 +0000
parents d5155eaa3f68
children
comparison
equal deleted inserted replaced
60:ad36f80e2ccf 61:42fcbcfca132
1 function [pdf,val] = genPDF(imSize,p,pctg,distType,radius,disp)
2
3 %[pdf,val] = genPDF(imSize,p,pctg [,distType,radius,disp])
4 %
5 % generates a pdf for a 1d or 2d random sampling pattern
6 % with polynomial variable density sampling
7 %
8 % Input:
9 % imSize - size of matrix or vector
10 % p - power of polynomial
11 % pctg - partial sampling factor e.g. 0.5 for half
12 % distType - 1 or 2 for L1 or L2 distance measure
13 % radius - radius of fully sampled center
14 % disp - display output
15 %
16 % Output:
17 % pdf - the pdf
18 % val - min sampling density
19 %
20 %
21 % Example:
22 % [pdf,val] = genPDF([128,128],2,0.5,2,0,1);
23 %
24 % (c) Michael Lustig 2007
25
26 % This file is used with the kind permission of Michael Lustig
27 % (mlustig@stanford.edu), and originally appeared in the
28 % SparseMRI toolbox, http://www.stanford.edu/~mlustig/ .
29 %
30 % $Id: genPDF.m 1040 2008-06-26 20:29:02Z ewout78 $
31
32
33 if nargin < 4
34 distType = 2;
35 end
36
37 if nargin < 5
38 radius = 0;
39 end
40
41 if nargin < 6
42 disp = 0;
43 end
44
45
46 minval=0;
47 maxval=1;
48 val = 0.5;
49
50 if length(imSize)==1
51 imSize = [imSize,1];
52 end
53
54 sx = imSize(1);
55 sy = imSize(2);
56 PCTG = floor(pctg*sx*sy);
57
58
59 if sum(imSize==1)==0 % 2D
60 [x,y] = meshgrid(linspace(-1,1,sy),linspace(-1,1,sx));
61 switch distType
62 case 1
63 r = max(abs(x),abs(y));
64 otherwise
65 r = sqrt(x.^2+y.^2);
66 r = r/max(abs(r(:)));
67 end
68
69 else %1d
70 r = abs(linspace(-1,1,max(sx,sy)));
71 end
72
73
74
75
76 idx = find(r<radius);
77
78 pdf = (1-r).^p; pdf(idx) = 1;
79 if floor(sum(pdf(:))) > PCTG
80 error('infeasible without undersampling dc, increase p');
81 end
82
83 % begin bisection
84 while(1)
85 val = minval/2 + maxval/2;
86 pdf = (1-r).^p + val; pdf(find(pdf>1)) = 1; pdf(idx)=1;
87 N = floor(sum(pdf(:)));
88 if N > PCTG % infeasible
89 maxval=val;
90 end
91 if N < PCTG % feasible, but not optimal
92 minval=val;
93 end
94 if N==PCTG % optimal
95 break;
96 end
97 end
98
99 if disp
100 figure,
101 subplot(211), imshow(pdf)
102 if sum(imSize==1)==0
103 subplot(212), plot(pdf(end/2+1,:));
104 else
105 subplot(212), plot(pdf);
106 end
107 end
108
109
110
111
112
113