annotate Problems/private/genSampling.m @ 51:217a33ac374e

(none)
author idamnjanovic
date Mon, 14 Mar 2011 16:52:27 +0000
parents
children
rev   line source
idamnjanovic@51 1 function [minIntrVec,stat,actpctg] = genSampling(pdf,iter,tol)
idamnjanovic@51 2
idamnjanovic@51 3 %[mask,stat,N] = genSampling(pdf,iter,tol)
idamnjanovic@51 4 %
idamnjanovic@51 5 % a monte-carlo algorithm to generate a sampling pattern with
idamnjanovic@51 6 % minimum peak interference. The number of samples will be
idamnjanovic@51 7 % sum(pdf) +- tol
idamnjanovic@51 8 %
idamnjanovic@51 9 % pdf - probability density function to choose samples from
idamnjanovic@51 10 % iter - number of tries
idamnjanovic@51 11 % tol - the deviation from the desired number of samples in samples
idamnjanovic@51 12 %
idamnjanovic@51 13 % returns:
idamnjanovic@51 14 % mask - sampling pattern
idamnjanovic@51 15 % stat - vector of min interferences measured each try
idamnjanovic@51 16 % actpctg - actual undersampling factor
idamnjanovic@51 17 %
idamnjanovic@51 18 % (c) Michael Lustig 2007
idamnjanovic@51 19
idamnjanovic@51 20 % This file is used with the kind permission of Michael Lustig
idamnjanovic@51 21 % (mlustig@stanford.edu), and originally appeared in the
idamnjanovic@51 22 % SparseMRI toolbox, http://www.stanford.edu/~mlustig/ .
idamnjanovic@51 23 %
idamnjanovic@51 24 % $Id: genSampling.m 1040 2008-06-26 20:29:02Z ewout78 $
idamnjanovic@51 25
idamnjanovic@51 26 % h = waitbar(0);
idamnjanovic@51 27
idamnjanovic@51 28 pdf(find(pdf>1)) = 1;
idamnjanovic@51 29 K = sum(pdf(:));
idamnjanovic@51 30
idamnjanovic@51 31 minIntr = 1e99;
idamnjanovic@51 32 minIntrVec = zeros(size(pdf));
idamnjanovic@51 33
idamnjanovic@51 34 for n=1:iter
idamnjanovic@51 35 tmp = zeros(size(pdf));
idamnjanovic@51 36 while abs(sum(tmp(:)) - K) > tol
idamnjanovic@51 37 tmp = rand(size(pdf))<pdf;
idamnjanovic@51 38 end
idamnjanovic@51 39
idamnjanovic@51 40 TMP = ifft2(tmp./pdf);
idamnjanovic@51 41 if max(abs(TMP(2:end))) < minIntr
idamnjanovic@51 42 minIntr = max(abs(TMP(2:end)));
idamnjanovic@51 43 minIntrVec = tmp;
idamnjanovic@51 44 end
idamnjanovic@51 45 stat(n) = max(abs(TMP(2:end)));
idamnjanovic@51 46 % waitbar(n/iter,h);
idamnjanovic@51 47 end
idamnjanovic@51 48
idamnjanovic@51 49 actpctg = sum(minIntrVec(:))/prod(size(minIntrVec));
idamnjanovic@51 50
idamnjanovic@51 51 % close(h);
idamnjanovic@51 52
idamnjanovic@51 53