/ase/se10.m
Objective C | 486 lines | 447 code | 39 blank | 0 comment | 52 complexity | 722410771c866ec17f243a2f8fa78698 MD5 | raw file
- function se10(parameter_name,parameter_values,varargin)
- %% run PARSE with several different values of a parameter
- % examples:
- % se10; % tests the code
- % se10('freq_off',-52:2:100,'dataset','fourtube_se');
- % se10('echo_time',1000:1000:11000,'dataset','fourtube_se');
- % se10('dataset','fourtube_se');
- % This is the currently the top level script for sePARSE code.
- % add your datasets and specific parameters in dataset_load.m
- % add new optional parameter try..catch blocks in se_opt.m
- % Varargin is a matlab keyword that is used to pull in optional name, value
- % pairs which are then parsed by the local function varparser
- % separse4 does most of the work
- % uses NIfTI toolbox from
- % http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8797&objectType=File
- %
- % msb 05 Feb '07
- % msb 19 Feb '07
- % msb 22 Feb '07
- %
- % dependencies:
- % se10
- % |- plotter
- % |- make_nii, from mathworks user site.
- % |- save_nii, from mathworks user site.
- % |- dataset_load
- % |- se_opts
- % |- varparser
- % |- load_procpar
- % '- query_procpar
- seversion=10;
- %% init some things...
- tic % start the process timer
- if exist('parameter_name','var')
- disp([parameter_name ' will be evaluated at ']);
- disp(parameter_values)
- else % if the code is called without parameters then execute test case
- disp('running sanity check.')
- parameter_name='test';
- parameter_values='0';
- s=['mv psis.mat psis' num2str(seversion) '.bak'];
- system(s);
- end
- try % flag to show progress or not with frequently updated figs
- plot_a.show=varparser('show',varargin);
- if plot_a.show
- disp('progress display is turned on.')
- else
- disp('progress display is turned off.')
- end
- catch
- plot_a.show=1;
- end
- if strcmp(parameter_name,'dataset')
- dataset=parameter_values;
- disp(['analyzing ' dataset])
- parameter_name=parameter_values;
- parameter_values=1;
- else
- try % choose the dataset using input value or default value...
- dataset=varparser('dataset',varargin);
- disp(['using dataset ' dataset])
- catch
- dataset='fourtube'; % an real example dataset to see if the code is working.
- disp(['using a default dataset of ' dataset])
- end
- end
- %% select dataset =======================================================
- % also contains some dataset specific defaults
- settings={}; % init with empty cell array that dataset_load fills in
- dataset_load % call the dataset and opt parameter loading script
- %% make a place to save the data
- savedir=['sweepresults_' num2str(seversion) '_' parameter_name]; % where should the results go
- try
- mkdirresult=mkdir( savedir );
- if mkdirresult==1
- disp(['running will write to directory ' savedir '. continuing in 5 sec.'])
- pause(5)
- end
- catch
- error('couldn not make dir for results')
- end
- %% run the sweep
- sweep_count=0; % parameter value trial counter
- figure(2);
- figx=size(parameter_values,2); % one little figure for each trial
- figy=4;
- yres=64;
- yrng=1:yres;
- for parameter_value=parameter_values % try each parameter value
- sweep_count=sweep_count+1;
- coefficientmaps=separse4(mr_data,plot_a,parameter_name,parameter_value,varargin{:},settings{:});
- paramvolume(:,:,sweep_count)=coefficientmaps;
- save([savedir '/' parameter_name num2str(parameter_value)],'coefficientmaps') % save current result
- set(0,'CurrentFigure',2);
- for idx=1:figy
- subplot('Position',[(sweep_count-1)/figx (idx-1)/figy (1/figx) (1/figy)]) % plot them all in a row.
- imagesc(coefficientmaps(yrng+((idx-1)*yres),:)) % display the current result
- axis off
- colormap gray;
- end
- drawnow
- elapsed_time=toc;
- disp([ 'elapsed time: ' num2str(elapsed_time)])
- end
- %% save all the results as a single 3D nifti format volume for inspection
- paramvolume(65:128,:,:)=paramvolume(65:128,:,:)*100000; % scaling factor because Mag image is small valued.
- save paramvolume paramvolume
- try
- nii=make_nii(paramvolume,[2 2 2]);
- save_nii(nii,[savedir '/sigvolume']);
- catch
- disp('could not make the nifti volumes. Do you have the nii tools installed?')
- end
- %% functions ==================================================
- function coefficientmaps = separse4(mr_data,plot_a,varargin)
- % example: cfm = separse2('fourtube',1)
- % mark bolding mbolding@uab.edu
- % 29 Nov '06 msb
- % 03 Jan '07 msb new trajectories
- % 19 Jan '07 msb converting to SE
- % 05 Feb '07 msb converting to SE second try, merging pre- and main code
- % 19 Feb '07 msb is "try catch" crazy
- % 22 Feb '07 msb is adding more parameters... SE third try rolling in
- % functions and making p into a structure
- % 27 Feb '07 msb SE works I think, now trying to search from echo and not
- % first data point.
- %
- % original:
- % D Twieg UAB Biomedical Engineering
- % uses hexagonal evaluation grid.
- %% constants ===========================================================
- % also see the tweakers
- plot_a.fig=1; % figure number for plot
- maxIter=200; % the maximum number of iterations
- sq2=sqrt(2);
- twopi=2*pi;
- % gam=4258.0;
- % gmax=3.643; % low-g
- % fov = 12.8; % field of view in cm calibration traj dependent
- % nres = 64; % nominal resolution for reconstruction
- %% tweakers <<<<<<<<<<<<<<<<<<<<
- trimoff=220; % how many initial data points to discard
- % trimoff=8280 % echo time for current spin echo data
- shift_fft=0; % calibration v. signal time offset
- dc_offset=21+1i*32;
- hammfilter=0; % hamming filter in k-space 1=on, 0=off
- %% look for varargs, optional inputs to the function =====================
- disp(['=== loading ' num2str(size(varargin,2)/2) ' optional settings ==='])
- se_opts % call the option loading script to parse the varargin
- disp('=== end of optional settings ===')
- %% try to grab some values from "procpar" ===============================
- % Varian stores all of the acquisition parameters in procpar.
- % it contains offsets, FOV, thk, powers, etc.
- try
- % using resto the resonance offset frequency as an example
- procpar_vars=load_procpar(mr_data);
- showprocvar('thk',procpar_vars)
- showprocvar('pss0',procpar_vars)
- catch
- disp('could not read values from the procpar file...')
- end
- %% bring all of the plots to the front at the start
- if(plot_a.show==1)
- figure(plot_a.fig)
- end
- %% load the calibration trajectory, kss =================================
- phi=[];
- ktraj=[];
- try
- load(calibrationdata); % brings in phi and ktraj
- catch
- error('Could not load the calibration data. Booo.');
- end
- disp(['processing map ' num2str(whichfid)]);
- % gmax=4.6833; % hi-g maximum read gradient used. Not the Varian
- % system maximum gradient!
- Ta=.068; % duration of acquisition in ms
- % Ta=.025; % post echo
- % Ta=0.045;
- % kf=nres/(2*fov);
- % delt=1.0/(gam*gmax*fov); % theoretical time step
- delt=1.0/198511.166253; % defined from sw in procpar
- numTimePts=floor(Ta/delt); % number of acquisition points
- kss=ktraj(shift_fft+trimoff+1:numTimePts+trimoff+shift_fft); % k trajectory with lead in and lead out trimmed off
- phi=real(phi)+imag(phi);
- phi=phi(shift_fft+trimoff+1:numTimePts+trimoff+shift_fft); % trim phi to the correct length and offset
- kss=kss*1.61; % empirical scaling
- %% build Hamming filter to correct for sample spacing ====================
- if hammfilter==1
- kr=abs(kss); % find the radius of each point
- kfe=max(kr); % kfe set as max radius
- kch=(0.54+0.46*cos(pi*kr/(sq2*kfe))); % circular Hamming, circumscribed
- kss=kss.*kch;
- end
- %% lengths and stopping tolerances for progressive conjugate =============
- % gradient search. For different rosettes, the lengths should be
- % integer multiples of a single "echo time" -- i.e., the number of
- % samples between successive crossings of the k-origin.
- % You can determine this length by plotting the magnitude of the signal
- % and counting the number of points between peaks. The first number, N1, in
- % NLIST should be an integer multiple of the echo time, and it should be
- % sufficiently long to give a complete circle of rosette
- % lobes, covering k-space more or less symmetrically. Check this by
- % plotting kss(1:N1).
- swoopsize=121.0; % one pass through the center of k space
- N1=floor(5*swoopsize); % one "slice"
- NLIST=[(N1*cumsum(ones(1,5)*4))';numTimePts]; % analyze a little data then a little more
- NLENGTHS = size(NLIST,1); % how many sections is it divided into
- % FLIST=[0.3 ; 0.12 ; 0.1 ; 0.1 ; 0.1]; % tighten the tolerance as more data is added
- FLIST=[0.2 ; 0.1 ; 0.05 ; 0.05 ; 0.05 ; 0.05]; % tighten the tolerance as more data is added
- whichplot='slice';plotter;fprintf('%s',whichplot)
- %% makehex1; % builds numVox, xx, yy, xa, ya, XX, YY hexagonal evaluation grid
- % makehex1.m
- % generate hexagonal sampling grid
- % D. Twieg
- fov = 12.8; % field of view in cm calibration traj dependent
- nres = 64; % nominal resolution for reconstruction
- delxnom=fov/nres;
- %% make hexagonal grid
- [xx,yy]=hexgrid80(delxnom);
- %% count hex grid points and define arrays
- RA=1.0*fov/2; % radius
- nl=1;
- xa=zeros(80*80,1);
- ya=xa;
- ia=xa;
- ja=xa;
- for nc=1:80 % number of columns
- for nr=1:80 % number of rows
- if sqrt(xx(nr,nc).^2+yy(nr,nc).^2) < RA % if it is in the radius
- xa(nl)=xx(nr,nc); % x coordinates
- ya(nl)=yy(nr,nc); % y coordinates
- ia(nl)=nr; % row index
- ja(nl)=nc; % col index
- nl=nl+1;
- end
- end
- end
- [XX,YY]=meshgrid(((1:64)-33)*fov/64); % standard rectilinear grid for data display
- %% initialize parameter arrays and structures
- numVox = nl-1; % numVox is the number of grid points
- Voxels = 1:numVox; % handy array of voxel indicies
- p.Mperp = zeros(1,numVox); % parameter estimate array. initial mag
- p.fRa = zeros(1,numVox); % parameter estimate array. freq and decay
- p.Rb = zeros(1,numVox); % parameter to account for 2*R2 prime
- p.echo_time = echo_time;
- ExponentA = zeros(size(xx)); % storage for hex grid representation
- ExponentB = zeros(size(xx));
- Magnitude = zeros(size(xx));
- CC = zeros(1,maxIter); % array to store current cost result
- %% load the data.=========================================================
- disp(['loading ' mr_data]);
- %try
- for signal_index=whichfid % load the set of wanted signals
- [the_data_RE,the_data_IM]=LOAD_FID(mr_data,signal_index);
- the_data(:,signal_index)=the_data_RE+1i*the_data_IM + dc_offset; % combine re and im parts and dc offset
- end
- % catch
- % error('could not load the data.')
- % end
- theSignal=sum(the_data(trimoff+1:numTimePts+trimoff,:),2)'; % trim the signal to the proper size
- theSignal=theSignal.*exp(-1i*phi)'; % phi correction
- theSignal=theSignal.*exp(-1i*twopi*freq_off*(1:numTimePts)*delt); % frequency offset
- theSignal=theSignal*sig_weight*numVox*delt/max(abs(theSignal)); % regularization
- clear the_data_RE the_data_IM
- %% Load psis, the basis functions representing the effect of k-traj on
- % signal. only need calculate psis once for a given traj and grid
- psis=psisGen(kss,numVox,numTimePts,xa,ya);
- %% prep done. main code.==================================================
- disp('setup is finished. beginning search...');
- old_timeLength = 0; % the time length from previous iteration
- iteration = 1; % record number of steps and estimations in parameter space
- CCold = 1e+15; % last "current cost"
- CCnew = 1e+14; % current "current cost"
- for nlen = 1:NLENGTHS % for each data length, this is the main code loop
- timeLength = NLIST(nlen);
- timeVec=1:timeLength;
- while (abs(CCold-CCnew) > FLIST(nlen)*CCold || timeLength ~= old_timeLength ) % while not converged or new data
- sigEstimate = estimateSignal(p,psis,timeVec); % estimate the signal
- sigdiff = (theSignal(timeVec)-sigEstimate); % estimate the error
-
- %---------------------parameter gradient evaluation----------------
- % define direction of 1-D search
- if (timeLength ~= old_timeLength) % more data, reset search direction
- g = calcGrad(psis,timeVec,sigdiff,p);
- h = g;
- whichplot='gradients';plotter;fprintf('%s',whichplot)
- else
- gold = g;
- g = calcGrad(psis,timeVec,sigdiff,p);
- ga = max(0,((g-gold)*g')/(gold*gold'));
- h = g + ga*h;
- whichplot='gradients';plotter;fprintf('%s',whichplot)
- end
- delpf = NewtonMeth(p,psis,timeVec,sigEstimate,theSignal,h); % how far do we step?
- p.Mperp = p.Mperp + delpf*(h(Voxels)); % take a step magnitude
- p.fRa = p.fRa + delpf*(h(Voxels+numVox)); % take a step freq and Ra
- p.Rb = p.Rb + real(delpf*(h(Voxels+numVox*2))); % take a step Rb
- for nl = Voxels % read parameter values into the hex grid locations
- Magnitude(ia(nl),ja(nl)) = p.Mperp(nl); % amplitude, complex magnitude image
- ExponentA(ia(nl),ja(nl)) = p.fRa(nl); % freq and part of R that may reverse after echo?
- ExponentB(ia(nl),ja(nl)) = p.Rb(nl); % R correction after echo, now we can find R2 and R2'
- end
- M0Rect =griddata(yy,xx,abs(Magnitude),XX,YY,'cubic'); % resample M on a rect grid
- frRect =griddata(yy,xx,imag(ExponentA)/(twopi*delt),XX,YY,'cubic'); % resample f on a rect grid
- R2Rect =griddata(yy,xx,real(-1*ExponentA)/delt,XX,YY,'cubic'); % resample R2 on a rect grid
- R2primeRect =griddata(yy,xx,real(-1*ExponentB)/delt,XX,YY,'cubic');
-
- %----------------------------------------
- sigEstimate = estimateSignal(p,psis,timeLength);
- sigdif = theSignal(1:timeLength)-sigEstimate;
- CC(iteration)= sigdif*sigdif';
- whichplot='progress';plotter;fprintf('%s',whichplot)
- old_timeLength = timeLength;
- CCold = CCnew;
- CCnew = CC(iteration);
- iteration = iteration + 1;
- end
- end
- %% all done, show results -----------------
- coefficientmaps=[frRect;M0Rect;R2Rect;R2primeRect];
- return % FIN!
- %% sub-functions ====================================
- %% Newtons method
- function delpf = NewtonMeth(p,psis,timeVec,signalEstimate,theSignal,h)
- numVox=size(p.Mperp,2);
- Voxels=1:numVox;
- delta = 1e-6; % step size for gradient estimation. <<<
- % delta = 1e-7;
- % delta = 1e-8;
- sigdiff = signalEstimate-theSignal(timeVec);
- f_0 = sigdiff*sigdiff';
- p0.echo_time=p.echo_time;
- p0.Mperp = p.Mperp + delta*(h(Voxels)); % M perp
- p0.fRa = p.fRa + delta*(h(Voxels+numVox)); % f and R
- p0.Rb = p.Rb + delta*(h(Voxels+numVox*2));% R prime
- sigdiff = estimateSignal(p0,psis,timeVec) - theSignal(timeVec);
- f_a_p = sigdiff*sigdiff';
- p0.Mperp = p.Mperp - delta*(h(Voxels)); % M perp
- p0.fRa = p.fRa - delta*(h(Voxels+numVox)); % f and R
- p0.Rb = p.Rb - delta*(h(Voxels+numVox*2));% R prime
- sigdiff = estimateSignal(p0,psis,timeVec) - theSignal(timeVec);
- f_a_n = sigdiff*sigdiff';
- f_1_0 = (f_a_p-f_a_n)/(2*delta); % First derivative @ alpha = 0
- f_1_a_p = (f_a_p-f_0)/delta; % First derivative @ alpha = delta/2
- f_1_a_n = (f_0-f_a_n)/delta; % First derivative @ alpha = -delta/2
- f_2_0 = (f_1_a_p-f_1_a_n)/delta; % Second derivative @ alpha = 0
- delpf = max(0,-f_1_0/f_2_0); % min point
- return
- %% signal estimation
- function sigEst = estimateSignal(p,psis,timeVec)
- sigEst=zeros(size(timeVec)); % array to hold the signal
- b=p.Mperp; % initial magnetization
- preechodecay =exp(p.fRa - real(p.Rb));
- postechodecay=exp(p.fRa + real(p.Rb)); % imag part p.fRa is freq, real parts are relaxation
- sigEst(1)=b*psis(:,1); % first time point
- for t=timeVec(2:end) % iterate over time points
- if t < p.echo_time
- b=b.*(preechodecay); % phase evolution and decay due to model
- else
- b=b.*(postechodecay); % after the echo some of the phase changes reverse
- end
- sigEst(t)=b*psis(:,t); % effect of motion through k-space
- end
- return
- %% gradient calculation
- function g = calcGrad(psis,timeVec,sigdiff,p)
- numVox=size(p.fRa,2);
- Voxels=1:numVox;
- CC_conj = conj(-sigdiff);
- fR_grad = zeros(numVox,1); % freq and decay
- Rp_grad = zeros(numVox,1); % R prime
- M = transpose(p.Mperp); % M perp
- wfpre = transpose(exp(p.fRa - real(p.Rb))); % before spin echo
- wfpost = transpose(exp(p.fRa + real(p.Rb))); % after spin echo
- wf_n = ones(numVox,1);
- M_grad = psis(:,1)*CC_conj(1);
- for t = timeVec
- if t < p.echo_time
- wf_n = wf_n.*wfpre;
- else
- wf_n = wf_n.*wfpost;
- end
- temp = wf_n.*psis(:,t)*CC_conj(t);
- M_grad = M_grad + temp;
- fR_grad = fR_grad + (t-1)*temp;
- if t > p.echo_time
- Rp_grad = Rp_grad + (t-p.echo_time)*temp;
- end
- end
- fR_grad = M.*fR_grad;
- Rp_grad = M.*Rp_grad;
- g(Voxels) = -conj(2*M_grad);
- g(Voxels+numVox) = -conj(2*fR_grad);
- g(Voxels+numVox*2) = -conj(2*Rp_grad);
- return
- function psis=psisGen(kss,numVox,numTimePts,xa,ya)
- try
- disp('looking for cached psis, delete the psis.mat file to recalculate');
- disp('psis depends on the k traj and evaluation grid');
- load psis
- catch
- disp('psis could not be loaded. regenerating psis...');
- kx=real(kss);
- ky=imag(kss);
- psis = zeros(numVox,numTimePts);
- % here psis is the array of ideal local responses to the encoding gradients for unit-magnitude on-resonance
- % magnetization vectors at each pixel (xa,ya) location. These are "basis"
- % functions of a sort.
- figure(3)
- for nloc = 1:numVox % for each location...
- psis(nloc,:) = exp(-2i*pi*(kx(1:numTimePts)*xa(nloc)+ky(1:numTimePts)*ya(nloc))); % effect at each Voxel due to k-traj
- if mod(nloc,500) == 0 % show calculation progress once every 500 steps
- imagesc(real(psis));
- drawnow;
- end
- end
- fprintf('%s','saving psis ')
- save('psis','psis')
- disp('done')
- end
- return
- function [xx,yy]=hexgrid80(delxnom)
- w=1.128/(2*delxnom);
- delx=1.0/(sqrt(3)*w); % circumscribed hexagon;
- % gives minimum bias error
- dely=delx*sqrt(3)/2;
- xx=zeros(80,80);
- yy=xx;
- for nc=1:80
- for nr=1:2:79
- yy(nr,nc)=(nr-41)*dely; % y locations
- yy(nr+1,nc)=(nr-40)*dely;
- xx(nr,nc)=(nc-33)*delx; % x locations
- xx(nr+1,nc)=((nc-33)+0.5)*delx;
- end
- end
- return
- %% display a value from procpar
- function showprocvar(procvar,procpar_vars)
- value=query_procpar(procvar,procpar_vars);
- disp([procvar ' was '])
- disp(value)
- return
- %% eq. 8.15, 8.18 from Haacke
- % { exp(-t*R2)*exp(-t*R2') | 0 < t < tau
- % Mtransverse(t) = Mtransverse(0) * { exp(-t*R2)*exp((t-TE)*R2') | tau < t < 2tau
- % { exp(-t*R2)*exp((TE-t)*R2') | 2tau < t