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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
