fourD_prior_solve

 fourD_prior_solve-- inverse solver to account for temporal
 and 3D spatial correlation
 img= fourD_prior_solve( inv_model, data1, data2)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time

 both data1 and data2 may be matrices (MxT) each of
  M measurements at T times
 if either data1 or data2 is a vector, then it is expanded
  to be the same size matrix

0001 function img= fourD_prior_solve( inv_model, data1, data2)
0002 % fourD_prior_solve-- inverse solver to account for temporal
0003 % and 3D spatial correlation
0004 % img= fourD_prior_solve( inv_model, data1, data2)
0005 % img        => output image (or vector of images)
0006 % inv_model  => inverse model struct
0007 % data1      => differential data at earlier time
0008 % data2      => differential data at later time
0009 %
0010 % both data1 and data2 may be matrices (MxT) each of
0011 %  M measurements at T times
0012 % if either data1 or data2 is a vector, then it is expanded
0013 %  to be the same size matrix
0014 
0015 % (C) 2007, Tao Dai and Andy Adler. Licenced under the GPL Version 2
0016 % $Id: fourD_prior_solve.html 2819 2011-09-07 16:43:11Z aadler $
0017 
0018 fwd_model= inv_model.fwd_model;
0019 time_steps = inv_model.fourD_prior.time_steps;
0020 
0021 % The one_step reconstruction matrix is cached
0022 one_step_inv = eidors_obj('get-cache', inv_model, 'fourD_prior_solve');
0023 if ~isempty(one_step_inv)
0024     eidors_msg('fourD_prior_solve: using cached value', 2);
0025 else
0026     img_bkgnd = calc_jacobian_bkgnd( inv_model );
0027     J = calc_jacobian( fwd_model, img_bkgnd);
0028 
0029     %%%%%Calc temporal correlation%%%%%%%%%%%%%%%%%%%%%%%%%%
0030     [gama,Gama_x] = calc_Gama( inv_model, data1,data2);
0031     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0032     inv_model.fourD_prior.temporal_corr_mat = Gama_x;
0033     inv_model.fourD_prior.temporal_weight = gama;
0034 
0035     one_step_inv= data_form( inv_model, J );
0036 
0037     eidors_obj('set-cache', inv_model, 'fourD_prior_solve', one_step_inv);
0038     eidors_msg('fourD_prior_solve: setting cached value', 2);
0039 end
0040 
0041 if fwd_model.normalize_measurements
0042     dva= data2 ./ data1 - 1;
0043 else
0044     dva= data2 - data1;
0045 end
0046 l_dva = size( dva, 2);
0047 
0048 idx= [-time_steps:time_steps]'*ones(1,l_dva) + ...
0049     ones(time_steps*2+1,1)*(1:l_dva);
0050 % replicate first and last measurements
0051 idx(idx<1) = 1;
0052 idx(idx>l_dva) = l_dva;
0053 
0054 dvat= reshape(dva(:,idx),[],l_dva);
0055 %%%%%% calc averaged measurement  %%%%%%%%%%%%%%%
0056 % ts_vec = -time_steps:time_steps;
0057 % weight = (gama.^abs(ts_vec));
0058 % y_avg = dvat * weight(:) / sum(weight);
0059 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0060 sol = one_step_inv *dvat(:,1+time_steps); %y_avg;%
0061 
0062 % create a data structure to return
0063 img.name= 'solved by fourD_prior_solve';
0064 img.elem_data = sol;
0065 img.inv_model= inv_model;
0066 img.fwd_model= fwd_model;
0067 
0068 % calculate the one_step_inverse using the data form
0069 % CovX * J' * inv(J*CovX*J' + CovZ)
0070 %   iRtR*Jt/(Ji*RtR*Jt +  hp^2*iW);
0071 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0072 function one_step_inv= data_form( inv_model, J)
0073 space_prior= inv_model.fourD_prior.space_prior;
0074 time_weight= inv_model.fourD_prior.temporal_weight;
0075 ts         = inv_model.fourD_prior.time_steps;
0076 
0077 space_Reg= feval(space_prior, inv_model);
0078 
0079 % %%%% calc space_Reg with inter-layer(z direction) coef considered
0080 
0081 P = calc_P_3d(inv_model,inv(space_Reg));
0082 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0083 
0084 iRtRJt_frac=  (P*J');
0085 JiRtRJt_frac= J*iRtRJt_frac;
0086 
0087 % JiRtRJt_mult accounts for different parts of the
0088 % frame being taken at different times
0089 delta_vec= calc_delta( inv_model, J);
0090 delta_vec1= delta_vec*ones(1,length(delta_vec));
0091 JiRtRJt_mult = time_weight.^abs(delta_vec1 - delta_vec1');
0092 
0093 [x,y]= meshgrid(-ts:ts,  -ts:ts);
0094 %     time_w_mat= time_weight.^abs(x-y)
0095 time_w_mat= inv_model.fourD_prior.temporal_corr_mat;
0096 
0097 JiRtRJt= kron( time_w_mat, JiRtRJt_frac .* JiRtRJt_mult );
0098 
0099 %FIXME: do we multiply be JiRtRJt_mult here?
0100 iRtRJt=  kron( time_w_mat(ts+1,:), iRtRJt_frac );
0101 
0102 iW   = kron( speye(1+2*ts), inv( ...
0103     calc_meas_icov( inv_model ) ));
0104 hp   = calc_hyperparameter( inv_model );
0105 
0106 one_step_inv= iRtRJt/(JiRtRJt +  hp^2*iW);
0107 
0108 % if measurements are taken at different times,
0109 % then calculate a delta of each wrt the centre time
0110 function delta_vec= calc_delta( inv_model, J)
0111 stimulation= inv_model.fwd_model.stimulation;
0112 n_N= size(J,1);
0113 
0114 if isfield(stimulation(1),'delta_time')
0115     delta_time= [stimulation(:).delta_time];
0116     if diff(delta_time) ~= 0;
0117         error('All time steps must be same for kalman filter');
0118     end
0119 else
0120     delta_time=0;
0121 end
0122 
0123 % sequence is a vector location of each stimulation in the frame
0124 if delta_time == 0
0125     seq= size(J,1);
0126 else
0127     for i=1:length(stimulation)
0128         seq(i) = size(stimulation(i).meas_pattern,1);
0129     end
0130     seq= cumsum( seq );
0131 end
0132 
0133 delta_time= cumsum(delta_time);
0134 
0135 delta_vec= zeros(size(J,1),1);
0136 seq= [0;seq(:)];
0137 for i=1:length(seq)-1
0138     delta_vec( (seq(i)+1):seq(i+1) )= delta_time(i);
0139 end
0140 
0141 % normalize so middle time is centre, and max time is 1
0142 delta_vec= (delta_vec - mean(delta_vec)) / ...
0143     (sum(delta_time) + eps );
0144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0145 % Calculate correlation matrix of images from measurements
0146 function [gama,Gama] = calc_Gama( inv_model, data1,data2)
0147 cov_n = inv_model.fourD_prior.cov_noise;%cov matrix of noise is
0148 % calculated in advance in real meas,
0149 % it can be obtained from baseline test.
0150 ts = inv_model.fourD_prior.time_steps;
0151 
0152 data = data2-data1;
0153 Gama_y = corrcoef(data);%corr matrix of data
0154 cov_y = cov(data');%cov matrix of data
0155 
0156 JcovxJt = (cov_y-cov_n);%cov_y=JcovxJt+cov_n
0157 
0158 % find out the best r to describe correlation between
0159 % adjacent images.
0160 r=.01:.01:.99;
0161 [k,l]= meshgrid(-ts:ts,-ts:ts);
0162 I = eye(ts*2+1);
0163 % cov_n_telda = kron(I,cov_n);
0164 % cov_y_telda = kron(Gama_y,cov_y);%
0165 MP = Gama_y*norm(cov_y,'fro') -I*norm( cov_n,'fro');%matrix precalculated
0166 for i =1:length(r)
0167     G = r(i).^abs(k-l);
0168     e(i) = norm(MP - G*norm(JcovxJt,'fro'),'fro'); %this is equivalent to above. but save mem.
0169 end
0170 [e_min,ind] = min(e);
0171 gama = r(ind)%correlation coefficient of adjacent img
0172 Gama = gama.^abs(k-l);%correlation matrix of image set
0173 
0174 function P_3d = calc_P_3d(inv_model,P_2d)
0175 
0176 % P_C = exponential_covar_prior( inv_model );
0177 P_C = calc_covar_prior( inv_model );%a simplified calculation of "exponential_covar_prior"
0178 
0179 P_N = sqrt(P_2d);
0180 P_3d = P_N*P_C*P_N;