view Problems/Cardiac_MRI_problem.m @ 77:62f20b91d870

add routines from sparco problems privite folder to {root}\util
some changes to ksvd vs rlsdla image denoising example

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