mk_GREIT_model

 MK_GREIT_MODEL: make EIDORS inverse models using the GREIT approach
   [imdl, weight]= mk_GREIT_model( mdl, radius, weight, options )

 Output: 
   imdl   - GREIT inverse model
   weight - value of the weight paramater chosed to satisfy the prescribed
            noise figure (NF). See options.noise_figure below.

 Parameters:
   mdl    - fwd model on which to do simulations, or
          - inv model (experimental), or
          - string specifying prepackaged models

   radius - requested weighting matrix  (recommend 0.2 for 16 electrodes)
   weight - weighting matrix (weighting of noise vs signal). Can be empty
            options.noise_figure is specified
   options- structure with fields:
     imgsz         - [xsz ysz] reconstructed image size in pixels 
                     (default: [32 32])
     square_pixels - forces square pixels if 1 (default: 0)
     Nsim          - number of training points (default: 1000)
     distr         - distribution of training points:
         0 -> original (as per GREITv1)
         1 -> random, centre-heavy 
         2 -> random, uniform (default)
         3 -> fixed, uniform (debug)
     target_size - size of simulated targets as proportion of mesh radius
         (default: 0.02). Can be specified as [min_size max_size] for 
         random variation
     target_plane - the (mean) height z at which simulation targets are
         placed. This controls the image plane. Default: mean electrode
         height
     target_offset - maximum allowed vertical displacement from the
         target_plane (default: 0). Can be specified as
         [down_offset up_offset].
     noise_figure - the noise figure (NF) to achieve. Overwrites weight 
         which will be optimised to achieve the target NF.     
     noise_figure_targets - circular target(s) to use for NF calculation
         as an array of coordinates and radius xyzr [4xN] (default: single
         target at the center at average electrode height with radius of
         opt.target_size. Note that multiple targets are simultaneously
         simulated in a single measurement, meaning they should not
         overlap.
     image_SNR - an alternative method (apart from the NF) to specify the 
         noise performance of the algorithm to achieve.    
     image_SNR_targets - circular targets used for image SNR calculation 
         see xyzr_targets in calc_image_SNR for more information
     extra_noise - extra noise samples (such as electrode movement)
     desired_solution_fn - specify a function to calculate the desired 
         image. It must have the signature:
         D = my_function( xyc, radius, options); 
         See CALC_GREIT_RM for details.
     keep_intermediate_results - if true, stores additional data of 
         reconstruction matrix computation to be used later on, 
         e.g. for faulty electrode compensation
     show_NF_chosen - Show the NF value chosen

 NOTE
   currently extra_noise is not supported
   currently weighting matrix must be scalar
               
 Examples
   fmdl = mk_library_model('adult_male_16el');
   fmdl.stimulation = mk_stim_patterns(16,1,'{ad}','{ad}',{},1);
   fmdl.normalize_measurements = 1;
   imdl = mk_GREIT_model(fmdl,0.25,5); % uses weight 5
   OR
   opt.noise_figure = 0.5; 
   imdl = mk_GREIT_model(fmdl,0.25,[],opt); % optimises weight for NF=0.5

 CITATION_REQUEST:
 AUTHOR: A Adler et al.
 TITLE: GREIT: a unified approach to 2D linear EIT reconstruction of lung
 images
 JOURNAL: Phys Meas
 YEAR: 2009
 VOL: 30
 NUM: 6
 PAGE: S35-55
 LINK: http://iopscience.iop.org/0967-3334/30/6/S03

 See also CALC_GREIT_RM

0001 function [imdl, weight]= mk_GREIT_model( fmdl, radius, weight, options )
0002 % MK_GREIT_MODEL: make EIDORS inverse models using the GREIT approach
0003 %   [imdl, weight]= mk_GREIT_model( mdl, radius, weight, options )
0004 %
0005 % Output:
0006 %   imdl   - GREIT inverse model
0007 %   weight - value of the weight paramater chosed to satisfy the prescribed
0008 %            noise figure (NF). See options.noise_figure below.
0009 %
0010 % Parameters:
0011 %   mdl    - fwd model on which to do simulations, or
0012 %          - inv model (experimental), or
0013 %          - string specifying prepackaged models
0014 %
0015 %   radius - requested weighting matrix  (recommend 0.2 for 16 electrodes)
0016 %   weight - weighting matrix (weighting of noise vs signal). Can be empty
0017 %            options.noise_figure is specified
0018 %   options- structure with fields:
0019 %     imgsz         - [xsz ysz] reconstructed image size in pixels
0020 %                     (default: [32 32])
0021 %     square_pixels - forces square pixels if 1 (default: 0)
0022 %     Nsim          - number of training points (default: 1000)
0023 %     distr         - distribution of training points:
0024 %         0 -> original (as per GREITv1)
0025 %         1 -> random, centre-heavy
0026 %         2 -> random, uniform (default)
0027 %         3 -> fixed, uniform (debug)
0028 %     target_size - size of simulated targets as proportion of mesh radius
0029 %         (default: 0.02). Can be specified as [min_size max_size] for
0030 %         random variation
0031 %     target_plane - the (mean) height z at which simulation targets are
0032 %         placed. This controls the image plane. Default: mean electrode
0033 %         height
0034 %     target_offset - maximum allowed vertical displacement from the
0035 %         target_plane (default: 0). Can be specified as
0036 %         [down_offset up_offset].
0037 %     noise_figure - the noise figure (NF) to achieve. Overwrites weight
0038 %         which will be optimised to achieve the target NF.
0039 %     noise_figure_targets - circular target(s) to use for NF calculation
0040 %         as an array of coordinates and radius xyzr [4xN] (default: single
0041 %         target at the center at average electrode height with radius of
0042 %         opt.target_size. Note that multiple targets are simultaneously
0043 %         simulated in a single measurement, meaning they should not
0044 %         overlap.
0045 %     image_SNR - an alternative method (apart from the NF) to specify the
0046 %         noise performance of the algorithm to achieve.
0047 %     image_SNR_targets - circular targets used for image SNR calculation
0048 %         see xyzr_targets in calc_image_SNR for more information
0049 %     extra_noise - extra noise samples (such as electrode movement)
0050 %     desired_solution_fn - specify a function to calculate the desired
0051 %         image. It must have the signature:
0052 %         D = my_function( xyc, radius, options);
0053 %         See CALC_GREIT_RM for details.
0054 %     keep_intermediate_results - if true, stores additional data of
0055 %         reconstruction matrix computation to be used later on,
0056 %         e.g. for faulty electrode compensation
0057 %     show_NF_chosen - Show the NF value chosen
0058 %
0059 % NOTE
0060 %   currently extra_noise is not supported
0061 %   currently weighting matrix must be scalar
0062 %
0063 % Examples
0064 %   fmdl = mk_library_model('adult_male_16el');
0065 %   fmdl.stimulation = mk_stim_patterns(16,1,'{ad}','{ad}',{},1);
0066 %   fmdl.normalize_measurements = 1;
0067 %   imdl = mk_GREIT_model(fmdl,0.25,5); % uses weight 5
0068 %   OR
0069 %   opt.noise_figure = 0.5;
0070 %   imdl = mk_GREIT_model(fmdl,0.25,[],opt); % optimises weight for NF=0.5
0071 %
0072 % CITATION_REQUEST:
0073 % AUTHOR: A Adler et al.
0074 % TITLE: GREIT: a unified approach to 2D linear EIT reconstruction of lung
0075 % images
0076 % JOURNAL: Phys Meas
0077 % YEAR: 2009
0078 % VOL: 30
0079 % NUM: 6
0080 % PAGE: S35-55
0081 % LINK: http://iopscience.iop.org/0967-3334/30/6/S03
0082 %
0083 % See also CALC_GREIT_RM
0084 
0085 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0086 % $Id: mk_GREIT_model.m 5870 2018-12-16 15:00:04Z aadler $
0087 
0088 citeme(mfilename);
0089 
0090 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0091 
0092 if nargin < 4, options = [];end
0093 [imdl,fmdl,imgs] = parse_fmdl(fmdl);
0094 options = parse_options(options,fmdl,imdl, weight);
0095 
0096 copt.cache_obj= { fmdl, imdl, imgs, radius, weight, options};
0097 copt.fstr = 'mk_GREIT_model';
0098 params = {fmdl, imdl, imgs, radius, weight, options};
0099 
0100 [imdl, weight] = eidors_cache(@mk_GREIT_model_calc, params, copt);
0101 
0102 
0103 function [imdl, weight]= mk_GREIT_model_calc( fmdl, imdl, imgs, radius, weight, opt)
0104 
0105 Nsim = opt.Nsim;
0106 [vi,vh,xyz,opt]= stim_targets(imgs, Nsim, opt );
0107 
0108 %Calculate rec_model (if absent)
0109 if ~isfield(imdl,'rec_model');
0110    opt.do_coarse2fine = 0;
0111    [imdl.rec_model imdl.fwd_model] = mk_pixel_slice(fmdl,opt.target_plane,opt);
0112    imdl.rec_model.nodes(:,3) = []; % the third dimension complicated display
0113    % medical orientation
0114    imdl.rec_model.mdl_slice_mapper.y_pts = fliplr(imdl.rec_model.mdl_slice_mapper.y_pts);
0115 end
0116 
0117 opt.rec_model = imdl.rec_model; % for desired image calculation
0118 
0119 imdl.solve = @solve_use_matrix;
0120 %
0121 
0122 if ~isempty(opt.noise_figure) || ~isempty(opt.image_SNR)
0123     if ~isempty(opt.noise_figure)
0124         % we'll optimise the weight for a given noise figure (NF)
0125         target = opt.noise_figure;
0126         NoisPerfName = 'Noise Figure';
0127     elseif ~isempty(opt.image_SNR)
0128         % we'll optimise the weight for a given image SNR
0129         NoisPerfName = 'Image SNR';
0130         target = opt.image_SNR;        
0131         if isfield(opt, 'SigmaN')
0132             imdl.hyperparameter.SigmaN = opt.SigmaN;
0133         end
0134         if isfield(opt, 'image_SNR_targets')
0135             imdl.hyperparameter.xyzr_targets = opt.image_SNR_targets;
0136         end
0137     else
0138         error('internal bug: shouldn''t get here');
0139     end
0140     
0141     if ~isempty(weight)
0142         eidors_msg('mk_GREIT_model: Using weight parameter as a guess, options.noise_figure or opt.image_SNR is non-empty');
0143     else
0144         if ~isempty(opt.noise_figure)
0145             weight = target;
0146         elseif ~isempty(opt.image_SNR)
0147             weight = 1/target;   % the inverse, as image SNR \propto 1/NF
0148         else
0149             error('internal bug: shouldn''t get here');
0150         end
0151     end
0152     
0153     xyzr = opt.noise_figure_targets;
0154     [jnk,vi_NF] = simulate_movement(imgs,xyzr');
0155     vi_NF = sum(vi_NF,2); % sum the targets
0156     eidors_msg(['mk_GREIT_model: Finding noise weighting for given ', NoisPerfName],1);
0157     eidors_msg('mk_GREIT_model: This will take a while...',1);
0158     f = @(X) to_optimise(vh,vi,xyz, radius, X, opt, imdl, target, vi_NF);
0159     fms_opts.TolFun = 0.01*target; %don't need higher accuracy
0160     % The first call can take a long time. Take it out of the loop to
0161     % allow progress messages.
0162     imdl.solve_use_matrix.RM = calc_GREIT_RM(vh,vi,xyz, radius, weight, opt);
0163     log_level = eidors_msg('log_level');
0164     if  log_level > 1
0165        log_level = eidors_msg( 'log_level', 1); % suppress messages
0166     end
0167     [weight, NF] = fminsearch(f, weight,fms_opts);
0168     eidors_msg(['mk_GREIT_model: Optimal solution gives ', NoisPerfName, '=' ... 
0169         num2str(NF+target) ' with weight=' num2str(weight)],1);
0170     assert((sqrt(NF) / target) < 0.01, ...
0171             ['Cannot find an accurate enough match for desired ', NoisPerfName]');
0172      eidors_msg( 'log_level', log_level); % restore
0173 end
0174 %
0175 [RM, PJt, M] = calc_GREIT_RM(vh,vi, xyz, radius, weight, opt );
0176 imdl.solve_use_matrix.RM = RM;
0177 if opt.keep_intermediate_results
0178    % store additional data to be used for faulty electrode compensation
0179    imdl.solve_use_matrix.PJt = PJt;
0180    imdl.solve_use_matrix.M = M;
0181    imdl.solve_use_matrix.X = inv(M); % Not sure why we are keeping X - perhaps remove(AA)
0182 end
0183 % imdl.solve_use_matrix.RM = resize_if_reqd(RM,inside,imdl.rec_model);
0184 imdl.jacobian_bkgnd = imgs;
0185 %imdl.solve_use_matrix.map = inside;
0186 imdl.hyperparameter.value = weight;     % store the applied weight as "hyperparameter" value
0187 if isfield(opt,'show_NF_chosen') && opt.show_NF_chosen
0188    xyzr = opt.noise_figure_targets;
0189    [jnk,vi_NF] = simulate_movement(imgs,xyzr');
0190    NF = calc_noise_figure(imdl,vh, vi_NF);
0191    imdl.show_NF_chosen = NF; 
0192    eidors_msg(['NF = ', num2str(NF), ' weight = ', num2str(weight)],1);
0193 end
0194 
0195 function out = to_optimise(vh,vi,xy,radius,weight, opt, imdl, ...
0196     target,vi_NF)
0197 
0198    % calculate GREIT matrix as usual
0199    imdl.solve_use_matrix.RM = calc_GREIT_RM(vh,vi,xy, radius, weight, opt);
0200    imdl.hyperparameter.value = weight;
0201 %    imdl.solve_use_matrix.RM = resize_if_reqd(RM,inside,imdl.rec_model);
0202    if ~isempty(opt.noise_figure)
0203       NF = calc_noise_figure(imdl,vh, vi_NF);
0204       eidors_msg(['NF = ', num2str(NF), ' weight = ', num2str(weight)],1);
0205    elseif ~isempty(opt.image_SNR)
0206       NF = calc_image_SNR(imdl);
0207       eidors_msg(['SNR = ', num2str(NF), ' weight = ', num2str(weight)],1);
0208    else
0209       error('internal bug: shouldn''t get here');       
0210    end
0211    out = (NF - target)^2;
0212 %    out = (mean(NF) - target)^2 + std(NF);
0213 
0214 
0215 function  imgs = get_prepackaged_fmdls( fmdl );
0216   switch fmdl
0217     case 'c=1;h=2;r=.08;ce=16;bg=1;st=1;me=1;nd'
0218       fmdl = ng_mk_cyl_models([2,1,0.18],[16,1],[0.05]); 
0219       fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0220       fmdl = mdl_normalize(fmdl,1);
0221       imgs= mk_image( fmdl, 1);
0222     otherwise
0223       error('specified fmdl (%s) is not understood', fmdl);
0224   end
0225 
0226 function [vi,vh,xyz,opt]= stim_targets(imgs, Nsim, opt );
0227     fmdl = imgs.fwd_model;
0228    ctr =  mean(fmdl.nodes);  
0229    maxx = max(abs(fmdl.nodes(:,1) - ctr(1)));
0230    maxy = max(abs(fmdl.nodes(:,2) - ctr(2)));
0231    if numel(opt.distr) > 1
0232       xyzr = opt.distr;
0233       xyzr(4,:) = calc_radius(mean([maxx,maxy]),opt,size(opt.distr,2));
0234    else
0235        switch opt.distr
0236            case 0 % original
0237                r = linspace(0,0.9, Nsim);
0238                th = r*4321; % want object to jump around in radius
0239                xyzr = [maxx*r.*cos(th); maxy*r.*sin(th);
0240                    opt.target_plane*ones(1,Nsim);
0241                    0.05*mean([maxx,maxy])*ones(1,Nsim)];
0242                
0243            case 1 %centre-heavy
0244                F = fourier_fit(opt.contour_boundary(:,1:2));
0245                v = linspace(0,1,Nsim*100+1); v(end)=[];
0246                pts = fourier_fit(F,v);
0247                idx_p = floor(rand(Nsim,1)*Nsim*100);
0248                xyzr = pts(idx_p,:)'.*repmat(rand(Nsim,1),[1 2])';
0249                xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0250                
0251                % TODO: What size is good here and how to figure it out?
0252                xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0253            case 2 %uniform
0254                %            F = fourier_fit(opt.contour_boundary(:,1:2));
0255                %            v = linspace(0,1,101); v(end)=[];
0256                %            pts = fourier_fit(F,v);
0257                pts = opt.contour_boundary(:,1:2);
0258                % avoid edges
0259                pts = 0.9*( pts - repmat(ctr(1:2),length(pts),1) ) + repmat(ctr(1:2),length(pts),1);
0260                % using maxx and maxy below would in general not produce a
0261                % uniform distribution
0262                lim = max(maxx, maxy);
0263                x = ctr(1) + (rand(Nsim*10,1)-0.5)*2*lim;
0264                y = ctr(2) + (rand(Nsim*10,1)-0.5)*2*lim;
0265                IN = inpolygon(x,y,pts(:,1),pts(:,2));
0266                xyzr(1,:) = x(find(IN,Nsim));
0267                xyzr(2,:) = y(find(IN,Nsim));
0268                xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0269                % TODO: What size is good here and how to figure it out?
0270                xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0271            case 3 % uniform, non-random
0272                %            F = fourier_fit(opt.elec_loc(:,1:2));
0273                %            v = linspace(0,1,101); v(end)=[];
0274                %            pts = fourier_fit(F,v);
0275                pts = opt.contour_boundary(:,1:2);
0276                lim = max(maxx, maxy);
0277                frac = polyarea(pts(:,1),pts(:,2)) / (2*lim)^2;
0278                [x,y] = ndgrid( linspace(-lim,lim,ceil(sqrt(Nsim/frac))), ...
0279                    linspace(-lim,lim,ceil(sqrt(Nsim/frac))));
0280                
0281                x = x+ctr(1); y = y + ctr(2);
0282                IN = inpolygon(x,y,pts(:,1),pts(:,2));
0283                xyzr(1,:) = x(find(IN));
0284                xyzr(2,:) = y(find(IN));
0285                xyzr(3,:) = calc_offset(opt.target_plane,opt,size(xyzr,2));
0286                % TODO: What size is good here and how to figure it out?
0287                xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,size(xyzr,2));
0288                eidors_msg(['mk_GREIT_model: Using ' num2str(size(xyzr,2)) ' points']);
0289            otherwise; error('GREIT: opt.distr no such case=%d',opt.distr);
0290        end
0291    end
0292    before = size(xyzr,2);
0293    [vh,vi,xyzr] = simulate_movement(imgs, xyzr);
0294    after = size(xyzr,2);
0295    if(after~=before)
0296        eidors_msg(['mk_GREIT_model: Now using ' num2str(after) ' points']);
0297    end
0298    xyz = xyzr(1:3,:);
0299 
0300 function z = calc_offset(z0,opt,Nsim)
0301     if opt.random_offset
0302         l_bnd = opt.target_offset(1);
0303         width = sum(opt.target_offset(1:2));
0304         z = z0 - l_bnd + rand(Nsim,1)*width;
0305     else
0306         z = z0*ones(Nsim,1);
0307     end
0308 
0309 function r = calc_radius(R,opt,Nsim)
0310    if opt.random_size
0311        min_sz = opt.target_size(1);
0312        max_sz = opt.target_size(2);
0313        range = max_sz - min_sz;
0314        r = (min_sz + rand(Nsim,1)*range)*R;
0315    else
0316        r = opt.target_size(1)*ones(Nsim,1)*R;
0317    end
0318            
0319    
0320    
0321 function RM = resize_if_reqd(RM,inside,rmdl)
0322    szRM = size(RM,1);
0323    if sum(inside) == szRM || ...
0324         szRM == size(rmdl.elems,1) || ...
0325         (isfield(rmdl,'coarse2fine') && szRM == size(rmdl.coarse2fine,2))
0326       % RM is fine
0327    elseif any(size(inside)==szRM) && any(size(inside) == 1)
0328       RM = RM(inside,:);
0329    else
0330       error('mismatch in size of provided RecMatrix');
0331    end
0332 
0333 
0334 function [imdl,fmdl,imgs] = parse_fmdl(fmdl);
0335    imdl = []; 
0336    if ischar(fmdl)
0337       imgs = get_prepackaged_fmdls( fmdl );
0338       fmdl = imgs.fwd_model;
0339    elseif isfield(fmdl,'type');
0340      switch fmdl.type
0341    %  if we get a fwd_model, assume uniform conductivity backgnd of 1
0342        case 'fwd_model'; imgs = mk_image( fmdl, 1);
0343    %  if we get an image, use it. It may have a non-uniform backgnd
0344        case 'image';     imgs = fmdl; % fmdl was an image
0345                          fmdl = imgs.fwd_model; % now it's a fmdl
0346        case 'inv_model'; imdl = fmdl;
0347                          fmdl = imdl.fwd_model;
0348                          imgs = calc_jacobian_bkgnd(imdl);
0349        otherwise; error('unrecognized eidors object');
0350      end
0351    else
0352       error('specified parameter must be an object or a string');
0353    end
0354    % Prepare model
0355    if isempty(imdl)
0356       imdl = select_imdl( fmdl,{'Basic GN dif'});
0357    end
0358    
0359    
0360     function opt = parse_options(opt,fmdl,imdl, weight)
0361 
0362     if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end
0363     if ~isfield(opt, 'square_pixels')
0364         opt.square_pixels = 0;
0365     end
0366     % Allow imdl.rec_model to overwrite options.imgsz
0367     if isfield(imdl,'rec_model') && ~isempty(imdl.rec_model)
0368         % this assumes rec_model is a rectangular grid, as it should
0369         opt.imgsz(1) = numel(unique(imdl.rec_model.nodes(:,1)))-1;
0370         opt.imgsz(2) = numel(unique(imdl.rec_model.nodes(:,2)))-1;
0371         try
0372             opt.imgsz(3) = numel(unique(imdl.rec_model.nodes(:,3)))-1;
0373         end
0374     end  
0375     
0376     if ~isfield(opt, 'distr'),     opt.distr = 3;       end 
0377     if ~isfield(opt, 'Nsim' ),     opt.Nsim  = 1000;    end
0378     if ~isfield(opt, 'noise_figure'), opt.noise_figure = []; end
0379     if ~isfield(opt, 'image_SNR'), opt.image_SNR = []; end
0380     if isempty(opt.noise_figure) && isempty(opt.image_SNR) && isempty(weight)
0381         error('EIDORS:WrongInput', ...
0382             'The weight parameter must be specified if opt.noise_figure or opt.image_SNR are empty or absent');
0383     end
0384     if isfield(opt,'extra_noise')
0385       error('mk_GREIT_model: doesn''t currently support extra_noise');
0386     end
0387     if ~isfield(opt, 'target_size')
0388         opt.target_size = 0.05;
0389     end
0390     if sum(size(opt.target_size)) > 2
0391         if opt.target_size(1) == opt.target_size(2);
0392             opt.random_size = false;
0393         else
0394             opt.random_size = true;
0395         end
0396     end
0397     if sum(size(opt.target_size)) == 2
0398             opt.random_size = false;
0399     end
0400     
0401     % Calculate the position of the electrodes
0402     Nelecs = length(fmdl.electrode);
0403     for i=1:Nelecs
0404        enodesi = fmdl.electrode(i).nodes;
0405        elec_loc(i,:) = mean( fmdl.nodes( enodesi,:),1 );
0406     end
0407     opt.elec_loc = elec_loc;
0408     
0409     if ~isfield(opt, 'target_plane')
0410           opt.target_plane = mean(elec_loc(:,3));
0411     else
0412         t = opt.target_plane;
0413         minnode = min(fmdl.nodes);
0414         maxnode = max(fmdl.nodes);
0415         if t<minnode(3) || t>maxnode(3)
0416             warning('options.target_plane is outside the model!');
0417             eidors_msg('mk_GREIT_model: Resorting to default target_plane');
0418             opt.target_plane = mean(elec_loc(:,3));
0419         end
0420     end
0421     if ~isfield(opt, 'target_offset')
0422         opt.target_offset = 0;
0423     end
0424     if sum(size(opt.target_offset)) == 2
0425         if opt.target_offset < 0, opt.target_offset = 0; end
0426         opt.target_offset(2) = opt.target_offset(1);
0427     end
0428     if any(opt.target_offset > 0)
0429         opt.random_offset = true;
0430     else
0431         opt.random_offset = false;
0432     end
0433 
0434     if ~isfield(opt,'noise_figure_targets');
0435        R = max(max(fmdl.nodes(:,1:2)) - min(fmdl.nodes(:,1:2)));
0436        xyzr = mean(fmdl.nodes);
0437        xyzr(3) = opt.target_plane;
0438        xyzr(4) = mean(opt.target_size)*0.5*R;
0439        opt.noise_figure_targets = xyzr;
0440     end
0441 
0442        
0443     if ~isfield(opt,'keep_intermediate_results');
0444        opt.keep_intermediate_results = false;
0445     end
0446     
0447     
0448     try, opt.normalize = fmdl.normalize_measurements;
0449     catch, 
0450         opt.normalize = 0;
0451         eidors_msg('mk_GREIT_model: fmdl.normalize_measurements not specified, assuming 0');
0452     end
0453     
0454     % find the boundary at target level (needed in many places)
0455     slc = mdl_slice_mesher(fmdl,[inf inf opt.target_plane]);
0456     bnd = find_boundary(slc.fwd_model);
0457     opt.contour_boundary = order_loop(slc.fwd_model.nodes(unique(bnd),:));
0458     
0459 
0460 function do_unit_test
0461 
0462 sidx= 1; subplot(4,4,sidx);
0463  do_performance_test; 
0464 
0465 % Create a 3D elliptical cylinder with 16 circular electrodes
0466 fmdl_1= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1]); %show_fem(fmdl);
0467 % Put two balls into the elliptical cylinder
0468 extra={'ball','solid ball = sphere(0.5,0.5,0.5;0.1);'};
0469 [fmdl_2,mat_idx]= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1],extra); 
0470 % Set the model to use adjacent current patterns
0471 stim = mk_stim_patterns(16,1,[0,1],[0,1],{}); 
0472 fmdl_1.stimulation = stim;
0473 fmdl_2.stimulation = stim;
0474 % Simulate homogeneous voltages (background conductivity = 0.5);
0475 img = mk_image(fmdl_2, 0.5); vh = fwd_solve(img); %show_fem(img);
0476 % Simulate inhomogeneous voltages (ball conductivity = 1.0);
0477 img.elem_data(mat_idx{2})= 1.0; vi = fwd_solve(img); 
0478 sidx= sidx+1; subplot(4,4,sidx);
0479 show_fem(img);
0480 % Reconstruct the image using GREITv1
0481 imdl= mk_common_gridmdl('GREITc1'); 
0482 img= inv_solve(imdl,vh,vi);
0483 sidx= sidx+1; subplot(4,4,sidx);
0484 show_slices(img)
0485 
0486 % Create a GREIT model for the ellipse
0487 opt.noise_figure = 0.5; opt.distr = 3;opt.square_pixels = 1; %other options are defaults
0488 fmdl_2 = mdl_normalize(fmdl_2,0);
0489 % use the true model (inverse crime)
0490 img_2 = mk_image(fmdl_2,0.5);
0491 imdl1 = mk_GREIT_model(img_2, 0.25, [], opt);
0492 img1= inv_solve(imdl1,vh,vi);  
0493 sidx= sidx+1; subplot(4,4,sidx);
0494 show_slices(img1);
0495 
0496 % now do the same but using image SNR and not NF
0497 opt = rmfield(opt,'noise_figure');
0498 opt.image_SNR = 1e-3; 
0499 weight = 90; % need to choose a weight that works with SNR
0500 imdl1 = mk_GREIT_model(img_2, 0.25, weight, opt);
0501 img1= inv_solve(imdl1,vh,vi);  
0502 sidx= sidx+1; subplot(4,4,sidx); show_slices(img1);
0503 
0504 unit_test_cmp('Expect no PJT or X', ~isfield(imdl1.solve_use_matrix, 'PJt') & ...
0505                                     ~isfield(imdl1.solve_use_matrix, 'X'), true);
0506 
0507 weight = [];
0508 opt.keep_intermediate_results = true;
0509 imdl1 = mk_GREIT_model(img_2, 0.25, weight, opt);
0510 img1= inv_solve(imdl1,vh,vi);  
0511 sidx= sidx+1; subplot(4,4,sidx); show_slices(img1);
0512 
0513 unit_test_cmp('Expect PJT and X', isfield(imdl1.solve_use_matrix, 'PJt') & ...
0514                                   isfield(imdl1.solve_use_matrix, 'X'), true);
0515 
0516 opt = rmfield(opt,{'image_SNR', 'keep_intermediate_results'}); opt.noise_figure = 0.5;
0517 
0518 % use honogenous model
0519 fmdl_1 = mdl_normalize(fmdl_1,0);
0520 imdl2 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0521 img2= inv_solve(imdl2,vh,vi); 
0522 sidx= sidx+1; subplot(4,4,sidx); show_slices(img2);
0523 
0524 
0525 % specify targets for NF calc
0526 opt.noise_figure_targets = [-.5 0 .5 .2;.5 0 .5 .2;];
0527 imdl3 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0528 img3= inv_solve(imdl3,vh,vi); 
0529 sidx= sidx+1; subplot(4,4,sidx); show_slices(img3);
0530 % cleanup
0531 opt = rmfield(opt,'noise_figure_targets');
0532 
0533 
0534 %% repeat with normalized data
0535 fmdl_2 = mdl_normalize(fmdl_2,1);
0536 % use the true model (inverse crime)
0537 imdl3 = mk_GREIT_model(mk_image(fmdl_2,0.5), 0.25, [], opt);
0538 img3= inv_solve(imdl3,vh,vi); 
0539 
0540 % use honogenous model
0541 fmdl_1 = mdl_normalize(fmdl_1,1);
0542 imdl4 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0543 img4= inv_solve(imdl4,vh,vi); 
0544 
0545 sidx= sidx+1; subplot(4,4,sidx);
0546 show_slices([img1 img2 img3 img4])
0547 
0548 
0549 %% Use a prepackaged model
0550 fmdl = mk_library_model('adult_male_16el_lungs');
0551 fmdl.stimulation = stim;
0552 fmdl = mdl_normalize(fmdl,1);
0553 img = mk_image(fmdl,1);
0554 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.3;
0555 vh = fwd_solve(img);
0556 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.4;
0557 vi = fwd_solve(img);
0558 
0559 
0560 fmdl2 = mk_library_model('adult_male_16el');
0561 fmdl2.stimulation = stim;
0562 fmdl2 = mdl_normalize(fmdl2,1);
0563 
0564 opt.imgsz = [50 30];
0565 opt.square_pixels = 1;
0566 imdl = mk_GREIT_model(fmdl2,0.25,3,opt);
0567 
0568 img = inv_solve(imdl,vh, vi);
0569 sidx= sidx+1; subplot(4,4,sidx);
0570 show_slices(img);
0571 
0572 
0573 % do the same again with image SNR and not NF
0574 opt = rmfield(opt,'noise_figure');
0575 opt.image_SNR = 1e-4; 
0576 weight = 0.5; % need to choose a weight that works with SNR
0577 imdl = mk_GREIT_model(fmdl2, 0.25, weight, opt);
0578 img = inv_solve(imdl,vh, vi);
0579 sidx= sidx+1; subplot(4,4,sidx); show_slices(img);
0580 
0581 opt.image_SNR_targets = [0.3 0.3  0.5 0.05;  0.3 -0.3 0.5 0.05; ...
0582                          0.3 -0.3 0.5 0.05; -0.3 -0.3 0.5 0.05; ...
0583                          0.3 0    0.5 0.05; -0.3  0   0.5 0.05; ...
0584                          0   0.3  0.5 0.05;  0   -0.3 0.5 0.05]';
0585 opt.image_SNR = 3e-4; 
0586 weight = 1E-2; % need to choose a weight that works with SNR
0587 imdl = mk_GREIT_model(fmdl2, 0.25, weight, opt);
0588 img = inv_solve(imdl,vh, vi);
0589 sidx= sidx+1; subplot(4,4,sidx); show_slices(img);
0590 
0591 
0592 function do_performance_test
0593 % Reconstruct GREIT Images
0594 imdl_v1 = mk_common_gridmdl('GREITc1');
0595 imdl_v1.inv_solve.calc_solution_error = false;
0596 
0597 % Reconstruct backprojection Images
0598 imdl_bp = mk_common_gridmdl('backproj');
0599 
0600 % Recosntruct with new GREIT
0601 % fmdl = ng_mk_cyl_models([2,1,0.05],[16,1],[0.05]);
0602 fmdl = mk_library_model('cylinder_16x1el_fine');
0603 fmdl.nodes = fmdl.nodes/15; % make radius 1;
0604 fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
0605 opt.noise_figure = 0.88;
0606 opt.target_size = 0.1;
0607 opt.distr = 0;
0608 imdl_gr = mk_GREIT_model(fmdl, 0.2, [], opt);
0609 
0610 opt = struct();
0611 opt.noise_figure = 0.5; % current recommendation
0612 imdl_def = mk_GREIT_model(fmdl,0.2,[],opt);
0613 
0614 opt.desired_solution_fn = 'GREIT_desired_img_original';
0615 imdl_org = mk_GREIT_model(fmdl,0.2,[],opt);
0616 
0617 test_performance( { imdl_v1, imdl_gr, imdl_def, imdl_org},fmdl );