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