# HG changeset patch
# User idamnjanovic
# Date 1300117404 0
# Node ID 2953097411d43b1cacefe2f9cf3b59cb8e6ed4db
# Parent  6a37442514c5947935164a60744fd63ade5292ba


diff -r 6a37442514c5 -r 2953097411d4 Problems/Cardiac_MRI_problem.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Problems/Cardiac_MRI_problem.m	Mon Mar 14 15:43:24 2011 +0000
@@ -0,0 +1,155 @@
+function data = Cardiac_MRI_Problem(varargin)
+% CHANGE!!!!PROB503  Shepp-Logan phantom, partial Fourier with sample mask,
+%         complex domain, total variation.
+%
+%   PROB503 creates a problem structure.  The generated signal will
+%   consist of a N = 256 by N Shepp-Logan phantom. The signal is
+%   sampled at random locations in frequency domain generated
+%   according to a probability density function.
+%
+%   The following optional arguments are supported:
+%
+%   PROB503('n',N,flags) is the same as above, but with a
+%   phantom of size N by N. The 'noseed' flag can be specified to
+%   suppress initialization of the random number generators. Both
+%   the parameter pair and flags can be omitted.
+%
+%   Examples:
+%   P = prob503;  % Creates the default 503 problem.
+%
+%   References:
+%
+%   [LustDonoPaul:2007] M. Lustig, D.L. Donoho and J.M. Pauly,
+%     Sparse MRI: The application of compressed sensing for rapid MR
+%     imaging, Submitted to Magnetic Resonance in Medicine, 2007.
+%
+%   [sparsemri] M. Lustig, SparseMRI,
+%     http://www.stanford.edu/~mlustig/SparseMRI.html
+%
+%   See also GENERATEPROBLEM.
+%
+%MATLAB SPARCO Toolbox.
+
+%   Copyright 2008, Ewout van den Berg and Michael P. Friedlander
+%   http://www.cs.ubc.ca/labs/scl/sparco
+%   $Id: prob503.m 1040 2008-06-26 20:29:02Z ewout78 $
+
+% Parse parameters and set problem name
+
+[opts,varg] = parseDefaultOpts(varargin{:});
+[parm,varg] = parseOptions(varg,{'noseed'},{'n','fold','sigma','slice'});
+n           = getOption(parm,'n',256);
+info.name   = 'Cardiac_MRI';
+opts.show   =   1;
+
+
+fold        = getOption(parm,'fold', 6);  % undersampling level
+sigma       = getOption(parm,'sigma', 0.05);;  % noise level 
+z           = getOption(parm,'slice', 5);;  % slice number (1-10) 
+szt         = 20; % number of time samples
+
+% Return problem name if requested
+if opts.getname, data = info.name; return; end;
+
+% Initialize random number generators
+if (~parm.noseed), randn('state',0); rand('twister',2000); end;
+
+% Set up the data
+% if allowed use variable density
+%pdf  = genPDF([n,n],5,0.1,2,0.1,0);
+
+
+
+%load heart images
+FS=filesep;
+TMPpath=pwd;
+    [pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
+    cd([pathstr1,FS,'data',FS,'images',FS,'Cardiac_MRI_dataset',FS,'Images']);
+    [filename,pathname] = uigetfile({'*.mat;'},'Select a patient MRI image set');
+    [pathstr, name, ext, versn] = fileparts(filename);
+load(filename);
+data.name=name;
+cd(TMPpath);
+
+% Set up the problem
+
+% Get 3D matrix of heart images (size 256x256, 20 frames) and stack them to
+% 2D matrix (256 x 256*20)
+data.signal          = reshape(sol_yxzt(:,:,z,:), [n n*szt]);
+
+% make a noise matrix 
+
+noise_var=sqrt(sigma*var(reshape(data.signal, [n*n*szt 1])));
+data.noise           = randn(n,n*szt)*noise_var + sqrt(-1)*randn(n,n*szt)*noise_var;
+
+% make a mask of random lines in phase encode and time domain random - vector
+% of 0 and 1 of size n*szt multiplied with vector of 1 of size n
+
+mask = rand(n*szt,1);
+mask(mask>(1-1/fold))=1;
+mask(mask<=(1-1/fold))=0;
+mask=(mask*ones(1,n))'; 
+data.op.mask         = opMask(mask);
+data.op.padding      = opPadding([n,n*szt],[n,n*szt]);
+
+% make an fft 2D dictionary. It will do 2D fft on evry image in the stack
+data.op.fft2d        = opKron(opDiag(szt,1), opFFT2C(n,n));
+
+% make measurement operator mask*padding*fft2d
+data.M               = opFoG(data.op.mask, data.op.padding, ...
+                             data.op.fft2d);
+
+% make a mesurement vector b = M* (signal + noise) where s+n is stack to 1d vector                     
+data.b               = data.M(reshape(data.signal + data.noise,[n*n*szt,1]),1);
+
+
+data                 = completeOps(data);
+
+% Additional information
+info.title           = 'Cardiac-MRI';
+info.thumb           = 'figcardiacProblem';
+info.citations       = {'LustDonoPaul:2007','sparsemri'};
+info.fig{1}.title    = 'Cardiac MRI';
+% info.fig{1}.filename = 'figProblemCardiac';
+% info.fig{2}.title    = 'Probability density function';
+% info.fig{2}.filename = 'figProblem503PDF';
+% info.fig{3}.title    = 'Sampling mask';
+% info.fig{3}.filename = 'figProblem503Mask';
+
+% Set the info field in data
+data.info = info;
+opts.figinc=1;
+% Plot figures
+if opts.update || opts.show
+  
+    %figure(opts.figno); opts.figno = opts.figno + opts.figinc;
+    
+    mov=reshape(data.signal/500, [n n szt]);
+  
+    implay(mov);
+    clear mov;
+  
+    %updateFigure(opts, info.fig{1}.title, info.fig{1}.filename);
+  
+    movMeas=reshape(abs(data.A(data.b,2))/500, [n n szt]);
+    implay(movMeas);
+    clear movMeas;
+%   figure(opts.figno); opts.figno = opts.figno + opts.figinc;
+%   imagesc(pdf), colormap gray;
+%   updateFigure(opts, info.fig{2}.title, info.fig{2}.filename)
+
+    implay(reshape(mask, [n n szt]));
+
+%   figure(opts.figno); opts.figno = opts.figno + opts.figinc;
+%   imagesc(mask), colormap gray
+%   updateFigure(opts, info.fig{3}.title, info.fig{3}.filename)
+%   
+%   if opts.update
+%      mn = min(min(data.signal + real(data.noise)));
+%      mx = max(max(data.signal + real(data.noise)));
+%      P = (data.signal + real(data.noise) - mn) / (mx - mn);
+%      P = scaleImage(P,128,128);
+%      P = P(1:2:end,1:2:end,:);
+%      thumbwrite(P, info.thumb, opts);
+%   end
+end