mk_GN_model

 MK_GN_MODEL: make EIDORS inverse models using the GREIT approach
   imdl = mk_GN_model(img,  opt, lambda)

 Output: 
   imdl      - GREIT inverse model

 Parameters:
   img       - image of the forward model including stimulation pattern, etc.
   options   - structure with fields:
     imgsz         - [xsz ysz] reconstructed image size in pixels 
                     (default: [32 32])
     square_pixels - forces square pixels if 1 (default: 0)
     RtRprior      - desired regularization prior:
                     'laplace' (default), 'noser' or 'tikhonov'
     noise_figure - the noise figure (NF) to achieve. 
     meas_icov    - the inverse of the noise covariance matrix.
   lambda    - set a fixed hyperparameter, if empty (default value) 
               the hyperparameter is chosen acording to the desired noise figure

 See also MK_GREIT_MODEL (from which this script was adapted)

0001 function imdl = mk_GN_model(img, opt, lambda)
0002 % MK_GN_MODEL: make EIDORS inverse models using the GREIT approach
0003 %   imdl = mk_GN_model(img,  opt, lambda)
0004 %
0005 % Output:
0006 %   imdl      - GREIT inverse model
0007 %
0008 % Parameters:
0009 %   img       - image of the forward model including stimulation pattern, etc.
0010 %   options   - structure with fields:
0011 %     imgsz         - [xsz ysz] reconstructed image size in pixels
0012 %                     (default: [32 32])
0013 %     square_pixels - forces square pixels if 1 (default: 0)
0014 %     RtRprior      - desired regularization prior:
0015 %                     'laplace' (default), 'noser' or 'tikhonov'
0016 %     noise_figure - the noise figure (NF) to achieve.
0017 %     meas_icov    - the inverse of the noise covariance matrix.
0018 %   lambda    - set a fixed hyperparameter, if empty (default value)
0019 %               the hyperparameter is chosen acording to the desired noise figure
0020 %
0021 % See also MK_GREIT_MODEL (from which this script was adapted)
0022 %
0023 
0024 % (C) 2016 Fabian Braun. License: GPL version 2 or version 3
0025 % $Id: mk_GN_model.m 5424 2017-04-25 17:45:19Z aadler $
0026 
0027     %% do unit testing?
0028     if ischar(img) && strcmpi(img, 'unit_test')
0029         do_unit_test();
0030         return;
0031     end    
0032 
0033     %% parse and prepare inputs
0034     if ~exist('lambda', 'var')
0035       lambda = [];
0036     end
0037     if strcmp(img.type, 'fwd_model')
0038         img = mk_image(img, 1);
0039     end
0040     fmdl = img.fwd_model;
0041     imdl = select_imdl( fmdl,{'Basic GN dif'});
0042    
0043     %% parse and default options
0044     opt = parse_options(opt,fmdl,imdl);
0045 
0046     %% Calculate rec_model (if absent)
0047     if ~isfield(img,'rec_model');
0048         opt.do_coarse2fine = 0;  
0049         [imdl.rec_model, imdl.fwd_model] = mk_pixel_slice(fmdl, opt.target_plane, opt);
0050         imdl.rec_model.nodes(:,3) = []; % the third dimension complicated display
0051         % medical orientation: NO, DO flip the x-axis of the model beforehand
0052         imdl.rec_model.mdl_slice_mapper.y_pts = fliplr(imdl.rec_model.mdl_slice_mapper.y_pts);
0053     else
0054         imdl.rec_model = img.rec_model;
0055     end
0056 
0057 
0058     %% create GaussNewton reconstruction matrix
0059     imdl = select_imdl(imdl, {'Basic GN dif'});
0060     
0061     imdl.fwd_model.coarse2fine = mk_coarse_fine_mapping(imdl.fwd_model, imdl.rec_model);
0062     
0063     imdl.inv_solve.calc_solution_error = 0 ;
0064     
0065     if isfield(img, 'elem_data') 
0066       % non-homogenous
0067       imdl.jacobian_bkgnd.value = img.elem_data;        
0068     else
0069       imdl.jacobian_bkgnd.value = 1;
0070     end
0071     
0072     % assign the desired prior
0073     if strcmpi(opt.RtRprior, 'laplace');
0074         imdl.RtR_prior = @prior_laplace;  
0075     elseif strcmpi(opt.RtRprior, 'noser');
0076         imdl.RtR_prior = @prior_noser;  
0077         imdl.prior_use_fwd_not_rec = true;    % else we have issues with c2f mapping
0078     elseif strcmpi(opt.RtRprior, 'tikhonov');
0079         imdl.RtR_prior = @prior_tikhonov;  
0080     else
0081         error(['undefined prior: ', opt.RtRprior]);
0082     end 
0083     
0084     % assign inverse noise covariance matrix
0085     if ~isempty(opt.meas_icov)
0086         imdl.meas_icov = opt.meas_icov;
0087     end
0088     
0089     %% determine hyperparameter (either via noise figure or set manually)
0090     if ~isempty(opt.noise_figure)        
0091         %% take the four elements closest to the center (taken from select_imdl)
0092         % BUT: take the elements at height of the BELT!
0093         xyz_elems = interp_mesh( imdl.fwd_model );
0094         ctr_elems = mean(xyz_elems, 1);
0095         ctr_elems(3) = opt.target_plane;  % at belt height!
0096         xyz_elems = xyz_elems - ones(size(xyz_elems,1),1)*ctr_elems;
0097         d_elems   = sqrt( sum( xyz_elems.^2, 2 ));
0098         [~, e_idx] = sort(d_elems);
0099 
0100         imdl.hyperparameter.tgt_elems = e_idx(1:4);
0101         imdl.hyperparameter.noise_figure = opt.noise_figure;
0102 
0103         sv_log = eidors_msg('log_level'); eidors_msg('log_level', 2);
0104         imdl.hyperparameter.value = choose_noise_figure( imdl );
0105         eidors_msg('log_level', sv_log);
0106     elseif ~isempty(opt.image_SNR)       
0107         imdl.hyperparameter.image_SNR = opt.image_SNR;
0108         imdl.hyperparameter.value = choose_image_SNR(imdl);
0109     else
0110         imdl.hyperparameter.value = lambda;
0111     end
0112 
0113 end
0114 
0115 
0116 function opt = parse_options(opt,fmdl,imdl)
0117 
0118     if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end
0119     if ~isfield(opt, 'square_pixels')
0120         opt.square_pixels = 0;
0121     end
0122     % Allow imdl.rec_model to overwrite options.imgsz
0123     if isfield(imdl,'rec_model') && ~isempty(imdl.rec_model)
0124         % this assumes rec_model is a rectangular grid, as it should
0125         opt.imgsz(1) = numel(unique(imdl.rec_model.nodes(:,1)))-1;
0126         opt.imgsz(2) = numel(unique(imdl.rec_model.nodes(:,2)))-1;
0127     end
0128     
0129     if ~isfield(opt, 'noise_figure'), opt.noise_figure = []; end
0130     
0131     if ~isfield(opt, 'meas_icov'), opt.meas_icov = []; end
0132     
0133     if ~isfield(opt, 'image_SNR'), opt.image_SNR = []; end
0134     
0135     % prior type for regularization
0136     if ~isfield(opt, 'RtRprior')
0137         opt.RtRprior = 'laplace'; 
0138     end
0139         
0140     % Calculate the position of the electrodes
0141     Nelecs = length(fmdl.electrode);
0142     for i=1:Nelecs
0143        enodesi = fmdl.electrode(i).nodes;
0144        elec_loc(i,:) = mean( fmdl.nodes( enodesi,:),1 );
0145     end
0146     opt.elec_loc = elec_loc;
0147     
0148     try
0149         opt.normalize = fmdl.normalize_measurements;
0150     catch 
0151         opt.normalize = 0;
0152         eidors_msg('mk_GN_model: fmdl.normalize_measurements not specified, assuming 0');
0153     end
0154     
0155     % electrode target planes
0156     if ~isfield(opt, 'target_plane')
0157           opt.target_plane = mean(elec_loc(:,3));
0158     else
0159         t = opt.target_plane;
0160         minnode = min(fmdl.nodes);
0161         maxnode = max(fmdl.nodes);
0162         if t<minnode(3) || t>maxnode(3)
0163             warning('options.target_plane is outside the model!');
0164             eidors_msg('mk_GN_model: Resorting to default target_plane');
0165             opt.target_plane = mean(elec_loc(:,3));
0166         end
0167     end
0168     
0169 end
0170 
0171 
0172 function do_unit_test()
0173 
0174     % Recosntruct with GREIT
0175     fmdl = mk_library_model('cylinder_16x1el_fine');
0176     fmdl.nodes = fmdl.nodes/15; % make radius 1;
0177     fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
0178     opt.noise_figure = 1.0;
0179     imdl_gr = mk_GREIT_model(fmdl, 0.2, [], opt);
0180 
0181     opt = struct();
0182     opt.noise_figure = 1.0; 
0183     imdl_gn_lap = mk_GN_model(fmdl, opt);
0184 
0185     opt = struct();
0186     opt.noise_figure = 1.0; 
0187     opt.RtRprior = 'tikhonov';
0188     imdl_gn_tik = mk_GN_model(fmdl, opt);
0189 
0190     test_performance( { imdl_gr, imdl_gn_lap, imdl_gn_tik}, fmdl );
0191 
0192 end