# HG changeset patch # User Ivan Damnjanovic lnx # Date 1301307581 -3600 # Node ID ee2a4d0f0c4cde7e7b49817d29b3ff0572e4c2bb # Parent f69ae88b8be5f011705fe9c5a40f793d121e123c# Parent 62f20b91d870e344e747b5615cd9654f866daa29 Merge diff -r f69ae88b8be5 -r ee2a4d0f0c4c util/sparco utils/completeOps.m --- /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; diff -r f69ae88b8be5 -r ee2a4d0f0c4c util/sparco utils/genPDF.m --- /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 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 + + + + + + diff -r f69ae88b8be5 -r ee2a4d0f0c4c util/sparco utils/genSampling.m --- /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))= 15*(i/n - 0.5)^3 + 0.4) + f(i,j) = f(i,j) + sqrt(var); + end + end +end + diff -r f69ae88b8be5 -r ee2a4d0f0c4c util/sparco utils/scalarToRGB.m --- /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 diff -r f69ae88b8be5 -r ee2a4d0f0c4c util/sparco utils/thumbFromOp.m --- /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 diff -r f69ae88b8be5 -r ee2a4d0f0c4c util/sparco utils/thumbPlot.m --- /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 diff -r f69ae88b8be5 -r ee2a4d0f0c4c util/sparco utils/thumbwrite.m --- /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); diff -r f69ae88b8be5 -r ee2a4d0f0c4c util/sparco utils/updateFigure.m --- /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