Mercurial > hg > smallbox

changeset 79:ee2a4d0f0c4c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Mon, 28 Mar 2011 11:19:41 +0100
parents f69ae88b8be5 (current diff) 62f20b91d870 (diff)
children 16df822019f1
files
diffstat 9 files changed, 485 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/sparco utils/completeOps.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,78 @@
+function data = completeOps(data)
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: completeOps.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+operators = {};
+flagM = 0; if isfield(data,'M'), flagM = 1; end;
+flagB = 0; if isfield(data,'B'), flagB = 1; end;
+
+if (~flagM) && (~flagB)
+   error('At least one of the operators M or B has be to given.');
+end
+
+% Define measurement matrix if needed
+if ~flagM
+  info   = data.B([],0);
+  data.M = opDirac(info{1});
+else
+  operators{end+1} = data.M;
+end
+
+% Define sparsity basis if needed
+if ~flagB
+  info   = data.M([],0);
+  data.B = opDirac(info{2});
+else
+  operators{end+1} = data.B;
+end
+
+% Define operator A if needed
+if ~isfield(data,'A')
+  if (length(operators) > 1)
+     data.A = opFoG(operators{:});
+  else
+     data.A = operators{1};
+  end
+end
+
+% Define empty solution if needed
+if ~isfield(data,'x0')
+  data.x0 = [];
+end
+
+% Define the operator size and string
+opInfo       = data.A([],0);
+data.sizeA   = [opInfo{1},opInfo{2}];
+opInfo       = data.B([],0);
+data.sizeB   = [opInfo{1},opInfo{2}];
+opInfo       = data.M([],0);
+data.sizeM   = [opInfo{1},opInfo{2}];
+data.op.strA = opToString(data.A);
+data.op.strB = opToString(data.B);
+data.op.strM = opToString(data.M);
+
+% Get the size of the desired signal
+if ~isfield(data,'signalSize')
+   if ~isfield(data,'signal')
+     error(['At least one of the fields signal ', ...
+            'or signalSize must be given.']);
+   end
+   data.signalSize = size(data.signal);
+end
+
+% Reconstruct signal from sparse coefficients
+if ~isfield(data,'reconstruct')
+   data.reconstruct = @(x) reshape(data.B(x,1),data.signalSize);
+end
+
+% Reorder the fields (sort alphabetically)
+m = fieldnames(data);
+m = sort(m);
+dataReorder = struct();
+for i=1:length(m)
+  eval(sprintf('dataReorder.%s = data.%s;',m{i},m{i}));
+end
+
+data = dataReorder;
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/sparco utils/genPDF.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,113 @@
+function [pdf,val] = genPDF(imSize,p,pctg,distType,radius,disp)
+
+%[pdf,val] = genPDF(imSize,p,pctg [,distType,radius,disp])
+%
+%	generates a pdf for a 1d or 2d random sampling pattern
+%	with polynomial variable density sampling
+%
+%	Input:
+%		imSize - size of matrix or vector
+%		p - power of polynomial
+%		pctg - partial sampling factor e.g. 0.5 for half
+%		distType - 1 or 2 for L1 or L2 distance measure
+%		radius - radius of fully sampled center
+%		disp - display output
+%
+%	Output:
+%		pdf - the pdf
+%		val - min sampling density
+%
+% 
+%	Example:
+%	[pdf,val] = genPDF([128,128],2,0.5,2,0,1);
+%
+%	(c) Michael Lustig 2007
+
+%   This file is used with the kind permission of Michael Lustig
+%   (mlustig@stanford.edu), and originally appeared in the
+%   SparseMRI toolbox, http://www.stanford.edu/~mlustig/ .
+%
+% $Id: genPDF.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+
+if nargin < 4
+	distType = 2;
+end
+
+if nargin < 5
+	radius = 0;
+end
+
+if nargin < 6
+	disp = 0;
+end
+
+
+minval=0;
+maxval=1;
+val = 0.5;
+
+if length(imSize)==1
+	imSize = [imSize,1];
+end
+
+sx = imSize(1);
+sy = imSize(2);
+PCTG = floor(pctg*sx*sy);
+
+
+if sum(imSize==1)==0  % 2D
+	[x,y] = meshgrid(linspace(-1,1,sy),linspace(-1,1,sx));
+	switch distType
+		case 1
+			r = max(abs(x),abs(y));
+		otherwise
+			r = sqrt(x.^2+y.^2);
+			r = r/max(abs(r(:)));			
+	end
+
+else %1d
+	r = abs(linspace(-1,1,max(sx,sy)));
+end
+
+
+
+
+idx = find(r<radius);
+
+pdf = (1-r).^p; pdf(idx) = 1;
+if floor(sum(pdf(:))) > PCTG
+	error('infeasible without undersampling dc, increase p');
+end
+
+% begin bisection
+while(1)
+	val = minval/2 + maxval/2;
+	pdf = (1-r).^p + val; pdf(find(pdf>1)) = 1; pdf(idx)=1;
+	N = floor(sum(pdf(:)));
+	if N > PCTG % infeasible
+		maxval=val;
+	end
+	if N < PCTG % feasible, but not optimal
+		minval=val;
+	end
+	if N==PCTG % optimal
+		break;
+	end
+end
+
+if disp
+	figure, 
+	subplot(211), imshow(pdf)
+	if sum(imSize==1)==0
+		subplot(212), plot(pdf(end/2+1,:));
+	else
+		subplot(212), plot(pdf);
+	end
+end
+
+
+
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/sparco utils/genSampling.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,53 @@
+function [minIntrVec,stat,actpctg] = genSampling(pdf,iter,tol)
+
+%[mask,stat,N] = genSampling(pdf,iter,tol)
+%
+% a monte-carlo algorithm to generate a sampling pattern with 
+% minimum peak interference. The number of samples will be
+% sum(pdf) +- tol
+%
+%	pdf - probability density function to choose samples from
+%	iter - number of tries
+%	tol  - the deviation from the desired number of samples in samples
+%
+% returns:
+%	mask - sampling pattern
+%	stat - vector of min interferences measured each try
+%	actpctg    - actual undersampling factor
+%
+%	(c) Michael Lustig 2007
+
+%   This file is used with the kind permission of Michael Lustig
+%   (mlustig@stanford.edu), and originally appeared in the
+%   SparseMRI toolbox, http://www.stanford.edu/~mlustig/ .
+%
+% $Id: genSampling.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+% h = waitbar(0);
+
+pdf(find(pdf>1)) = 1;
+K = sum(pdf(:));
+
+minIntr = 1e99;
+minIntrVec = zeros(size(pdf));
+
+for n=1:iter
+	tmp = zeros(size(pdf));
+	while abs(sum(tmp(:)) - K) > tol
+		tmp = rand(size(pdf))<pdf;
+	end
+	
+	TMP = ifft2(tmp./pdf);
+	if max(abs(TMP(2:end))) < minIntr
+		minIntr = max(abs(TMP(2:end)));
+		minIntrVec = tmp;
+	end
+	stat(n) = max(abs(TMP(2:end)));
+      % waitbar(n/iter,h);
+end
+
+actpctg = sum(minIntrVec(:))/prod(size(minIntrVec));
+
+% close(h);
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/sparco utils/pwsmoothfield.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,32 @@
+function f = pwsmoothfield(n,var,alpha)
+% PWSMOOTHFIELD(N,VAR,ALPHA)
+%   Generate an image of piecewise smooth filtered white noise
+%
+%   N = sidelength of the field in pixels
+%   VAR = variance of original white nose
+%   ALPHA = fraction of FFT coefficents to keep
+%
+%   Returns an N-by-N array
+%
+%   This file is used with the kind permission of Stephen J. Wright
+%   (swright@cs.wisc.edu), and was originally included in the GPSR
+%   v1.0 distribution: http://www.lx.it.pt/~mtf/GPSR .
+
+% $Id: pwsmoothfield.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+f = sqrt(var)*randn(n);
+F = fft2(f);
+a = ceil(n*alpha/2);
+b = fix(n*(1-alpha));
+F(a+1:a+b,:) = 0;
+F(:,a+1:a+b) = 0;
+f = real(ifft2(F));
+
+for i = 1:n
+    for j = 1:n
+        if (j/n >= 15*(i/n - 0.5)^3 + 0.4)
+            f(i,j) = f(i,j) + sqrt(var);
+        end
+    end
+end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/sparco utils/scalarToRGB.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,18 @@
+function s = scalarToRGB(x,colors)
+% input values are assumed to lie between 0 and 1
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: scalarToRGB.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+l  = size(colors,1);
+m  = size(x,1);
+n  = size(x,2);
+s  = zeros(m,n,3);
+
+for i=1:m
+  for j=1:n
+     idx = max(1,min(l,1+floor((l-1) * x(i,j))));
+     s(i,j,:) = colors(idx,:);
+  end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/sparco utils/thumbFromOp.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,42 @@
+function P = thumbFromOp(op,m,n,sm,sn,grayscale)
+% Output matrix P, of size m x n, sampled at
+% sm x sn upper top
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: thumbFromOp.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+if nargin < 6, grayscale = 0; end
+    
+info = op([],0);
+
+sm = min(info{1},sm);
+sn = min(info{2},sn);
+M  = zeros(sm,sn);
+
+for i=1:sn
+   v = zeros(info{2},1); v(i) = 1;
+   w = real(op(v,1));
+
+   M(:,i) = w(1:sm);
+end
+
+mn = min(min(M));
+mx = max(max(M));
+M  = (M - mn) / (mx-mn);
+
+idxm = floor(linspace(1,sm+1,m+1)); idxm = idxm(1:end-1);
+idxn = floor(linspace(1,sn+1,n+1)); idxn = idxn(1:end-1);
+
+if grayscale
+   P = 1-M(idxm,idxn);
+else
+   clrmap = hsv;
+   M = 1 + round(M * (length(clrmap)-1));
+   P = zeros(m,n,3);
+   for j1=1:m
+       for j2=1:n
+           P(j1,j2,:) = clrmap(M(idxm(j1),idxn(j2)),:);
+       end
+   end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/sparco utils/thumbPlot.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,38 @@
+function P = thumbPlot(P,x,y,color)
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: thumbPlot.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+m = size(P,1);
+n = size(P,2);
+if (size(P,3) == 0) & (length(color) == 3)
+  % Convert to gray-scale
+  color = 0.30*color(1) + 0.59*color(2) + 0.11*color(3);
+end
+
+mnx = min(x);  % Minimum x
+mxx = max(x);  % Maximum x
+mny = min(y);  % Minimum y
+mxy = max(y);  % Maximum y
+dy  = (mxy - mny) * 0.1;   % Offset on vertical axis
+sx  = (mxx - mnx) * 1.0;   % Scale of horizontal axis
+sy  = (mxy - mny) * 1.2;   % Scale of vertical axis
+
+if (sx < 1e-6), sx = 1; end
+if (sy < 1e-6), sy = 1; end
+
+for i=1:length(x)-1
+   x0 = floor(1 + (n-1) * (x(i  ) - mnx) / sx);
+   x1 = floor(1 + (n-1) * (x(i+1) - mnx) / sx);
+   y0 = floor(    (n-1) * (y(i  ) - mny + dy) / sy);
+   y1 = floor(    (n-1) * (y(i+1) - mny + dy) / sy);
+   
+   samples = 1+2*max(abs(x1-x0)+1,abs(y1-y0)+1);
+   c       = linspace(0,1,samples);
+   idx     = round((1-c)*x0 + c*x1);
+   idy     = n - round((1-c)*y0 + c*y1);
+   for j=1:samples
+      P(idy(j),idx(j),:) = color;
+   end
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/sparco utils/thumbwrite.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,18 @@
+function thumbwrite(data,name,opts)
+
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: thumbwrite.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+%
+% data       Thumbnail data in range 0-1
+% name       Name of file (no extension)
+% opts
+%   .thumbtype  Type of image (png, eps, ps, ...)
+%   .thumbdir   Output directory
+%
+
+[type,ext] = getFigureExt(opts.thumbtype);
+data = round(data * 255) / 255;
+imwrite(data,[opts.thumbpath,name,'.',ext],type);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/sparco utils/updateFigure.m	Mon Mar 28 11:19:41 2011 +0100
@@ -0,0 +1,93 @@
+function updateFigure(opts, figTitle, figFilename)
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: updateFigure.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+% Ensure default values are available
+opts.linewidth  = getOption(opts,'linewidth', []);
+opts.fontsize   = getOption(opts,'fontsize',  []);
+opts.markersize = getOption(opts,'markersize',[]);
+
+% Output the plots
+if opts.update
+  % Set the line width, font size and marker size
+  chld = [gca; get(gca,'Children')];
+  lnwd = ones(length(chld),1) * NaN;
+  fnts = ones(length(chld),1) * NaN;
+  mrks = ones(length(chld),1) * NaN;
+  for i=1:length(chld)
+    conf = get(chld(i));
+    if ~isempty(opts.linewidth) && isfield(conf,'LineWidth')
+      lnwd(i) = get(chld(i),'LineWidth');
+      if (lnwd(i) == 0.5) % Default
+        set(chld(i),'Linewidth',opts.linewidth);
+      end
+    end
+    if ~isempty(opts.fontsize) && isfield(conf,'FontSize')
+      fnts(i) = get(chld(i),'FontSize');
+      if (fnts(i) == 10) % Default
+        set(chld(i),'FontSize',opts.fontsize);
+      end
+    end
+    if ~isempty(opts.markersize) && isfield(conf,'MarkerSize')
+      mrks(i) = get(chld(i),'MarkerSize');
+      if (mrks(i) == 6) % Default
+        set(chld(i),'MarkerSize',opts.markersize);
+      end
+    end
+  end
+  
+  for i=1:length(opts.figtype)
+     updateFigureType(opts.update, 0, opts.figtype{i}, ...
+                      opts.figpath, figTitle, figFilename);
+  end
+
+  % Restore the line-widths, font size
+  for i=1:length(chld)
+    if ~isnan(lnwd(i))
+      set(chld(i),'LineWidth',lnwd(i));
+    end
+    if ~isnan(fnts(i))
+      set(chld(i),'FontSize',fnts(i));
+    end
+    if ~isnan(mrks(i))
+      set(chld(i),'MarkerSize',mrks(i));
+    end
+  end
+    
+end
+
+% Show the plot
+if opts.show
+  updateFigureType(0,opts.show,'','',figTitle,'');
+end
+
+
+
+function updateFigureType(update,show,figtype,figpath,figTitle,figFilename)
+filename = [figpath,figFilename];
+
+switch lower(figtype)
+ case {'pdf'}
+   cmdPostprocess = sprintf('!pdfcrop %s.pdf %s.pdf >& /dev/null', ...
+                            filename, filename);
+ otherwise
+   cmdPostprocess = [];
+end
+
+[figtype,figext] = getFigureExt(figtype);
+
+% Print the figure for output (without title)
+if update
+  evalc(sprintf('print -d%s %s.%s;', figtype, filename, figext));
+  
+  if ~isempty(cmdPostprocess)
+    eval(cmdPostprocess);
+  end
+end
+
+% Add title if needed
+if show
+  title(figTitle);
+end