Description of mk_GREIT

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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
