mk_GREIT_model

PURPOSE ^

MK_GREIT_MODEL: make EIDORS inverse models using the GREIT approach

SYNOPSIS ^

function [imdl, weight]= mk_GREIT_model( fmdl, radius, weight, options )

DESCRIPTION ^

 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, default)
         1 -> random, centre-heavy 
         2 -> random, uniform
         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

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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, default)
0025 %         1 -> random, centre-heavy
0026 %         2 -> random, uniform
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 %
0058 % NOTE
0059 %   currently extra_noise is not supported
0060 %   currently weighting matrix must be scalar
0061 %
0062 % Examples
0063 %   fmdl = mk_library_model('adult_male_16el');
0064 %   fmdl.stimulation = mk_stim_patterns(16,1,'{ad}','{ad}',{},1);
0065 %   fmdl.normalize_measurements = 1;
0066 %   imdl = mk_GREIT_model(fmdl,0.25,5); % uses weight 5
0067 %   OR
0068 %   opt.noise_figure = 0.5;
0069 %   imdl = mk_GREIT_model(fmdl,0.25,[],opt); % optimises weight for NF=0.5
0070 %
0071 % CITATION_REQUEST:
0072 % AUTHOR: A Adler et al.
0073 % TITLE: GREIT: a unified approach to 2D linear EIT reconstruction of lung
0074 % images
0075 % JOURNAL: Phys Meas
0076 % YEAR: 2009
0077 % VOL: 30
0078 % NUM: 6
0079 % PAGE: S35-55
0080 % LINK: http://iopscience.iop.org/0967-3334/30/6/S03
0081 %
0082 % See also CALC_GREIT_RM
0083 
0084 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0085 % $Id: mk_GREIT_model.m 5419 2017-04-25 13:17:11Z aadler $
0086 
0087 citeme(mfilename);
0088 
0089 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0090 
0091 if nargin < 4, options = [];end
0092 [imdl,fmdl,imgs] = parse_fmdl(fmdl);
0093 options = parse_options(options,fmdl,imdl, weight);
0094 
0095 copt.cache_obj= { fmdl, imdl, imgs, radius, weight, options};
0096 copt.fstr = 'mk_GREIT_model';
0097 params = {fmdl, imdl, imgs, radius, weight, options};
0098 
0099 [imdl, weight] = eidors_cache(@mk_GREIT_model_calc, params, copt);
0100 
0101 
0102 function [imdl, weight]= mk_GREIT_model_calc( fmdl, imdl, imgs, radius, weight, opt)
0103 
0104 Nsim = opt.Nsim;
0105 [vi,vh,xyz,opt]= stim_targets(imgs, Nsim, opt );
0106 
0107 %Calculate rec_model (if absent)
0108 if ~isfield(imdl,'rec_model');
0109    opt.do_coarse2fine = 0;
0110    [imdl.rec_model imdl.fwd_model] = mk_pixel_slice(fmdl,opt.target_plane,opt);
0111    imdl.rec_model.nodes(:,3) = []; % the third dimension complicated display
0112    % medical orientation
0113    imdl.rec_model.mdl_slice_mapper.y_pts = fliplr(imdl.rec_model.mdl_slice_mapper.y_pts);
0114 end
0115 
0116 opt.rec_model = imdl.rec_model; % for desired image calculation
0117 
0118 imdl.solve = @solve_use_matrix;
0119 %
0120 
0121 if ~isempty(opt.noise_figure) || ~isempty(opt.image_SNR)
0122     if ~isempty(opt.noise_figure)
0123         % we'll optimise the weight for a given noise figure (NF)
0124         target = opt.noise_figure;
0125         NoisPerfName = 'Noise Figure';
0126     elseif ~isempty(opt.image_SNR)
0127         % we'll optimise the weight for a given image SNR
0128         NoisPerfName = 'Image SNR';
0129         target = opt.image_SNR;        
0130         if isfield(opt, 'SigmaN')
0131             imdl.hyperparameter.SigmaN = opt.SigmaN;
0132         end
0133         if isfield(opt, 'image_SNR_targets')
0134             imdl.hyperparameter.xyzr_targets = opt.image_SNR_targets;
0135         end
0136     else
0137         error('internal bug: shouldn''t get here');
0138     end
0139     
0140     if ~isempty(weight)
0141         eidors_msg('mk_GREIT_model: Using weight parameter as a guess, options.noise_figure or opt.image_SNR is non-empty');
0142     else
0143         if ~isempty(opt.noise_figure)
0144             weight = target;
0145         elseif ~isempty(opt.image_SNR)
0146             weight = 1/target;   % the inverse, as image SNR \propto 1/NF
0147         else
0148             error('internal bug: shouldn''t get here');
0149         end
0150     end
0151     
0152     xyzr = opt.noise_figure_targets;
0153     [jnk,vi_NF] = simulate_movement(imgs,xyzr');
0154     vi_NF = sum(vi_NF,2); % sum the targets
0155     eidors_msg(['mk_GREIT_model: Finding noise weighting for given ', NoisPerfName],1);
0156     eidors_msg('mk_GREIT_model: This will take a while...',1);
0157     f = @(X) to_optimise(vh,vi,xyz, radius, X, opt, imdl, target, vi_NF);
0158     fms_opts.TolFun = 0.01*target; %don't need higher accuracy
0159     % The first call can take a long time. Take it out of the loop to
0160     % allow progress messages.
0161     imdl.solve_use_matrix.RM = calc_GREIT_RM(vh,vi,xyz, radius, weight, opt);
0162     log_level = eidors_msg('log_level');
0163     if  log_level > 1
0164        log_level = eidors_msg( 'log_level', 1); % suppress messages
0165     end
0166     [weight, NF] = fminsearch(f, weight,fms_opts);
0167     eidors_msg(['mk_GREIT_model: Optimal solution gives ', NoisPerfName, '=' ... 
0168         num2str(NF+target) ' with weight=' num2str(weight)],1);
0169     assert((sqrt(NF) / target) < 0.01, ...
0170             ['Cannot find an accurate enough match for desired ', NoisPerfName]');
0171      eidors_msg( 'log_level', log_level); % restore
0172 end
0173 %
0174 [RM, PJt, M] = calc_GREIT_RM(vh,vi, xyz, radius, weight, opt );
0175 imdl.solve_use_matrix.RM = RM;
0176 if opt.keep_intermediate_results
0177    % store additional data to be used for faulty electrode compensation
0178    imdl.solve_use_matrix.PJt = PJt;
0179    imdl.solve_use_matrix.X = inv(M);
0180 end
0181 % imdl.solve_use_matrix.RM = resize_if_reqd(RM,inside,imdl.rec_model);
0182 imdl.jacobian_bkgnd = imgs;
0183 %imdl.solve_use_matrix.map = inside;
0184 imdl.hyperparameter.value = weight;     % store the applied weight as "hyperparameter" value
0185 
0186 function out = to_optimise(vh,vi,xy,radius,weight, opt, imdl, ...
0187     target,vi_NF)
0188 
0189    % calculate GREIT matrix as usual
0190    imdl.solve_use_matrix.RM = calc_GREIT_RM(vh,vi,xy, radius, weight, opt);
0191    imdl.hyperparameter.value = weight;
0192 %    imdl.solve_use_matrix.RM = resize_if_reqd(RM,inside,imdl.rec_model);
0193    if ~isempty(opt.noise_figure)
0194       NF = calc_noise_figure(imdl,vh, vi_NF);
0195       eidors_msg(['NF = ', num2str(NF), ' weight = ', num2str(weight)],1);
0196    elseif ~isempty(opt.image_SNR)
0197       NF = calc_image_SNR(imdl);
0198       eidors_msg(['SNR = ', num2str(NF), ' weight = ', num2str(weight)],1);
0199    else
0200       error('internal bug: shouldn''t get here');       
0201    end
0202    out = (NF - target)^2;
0203 %    out = (mean(NF) - target)^2 + std(NF);
0204 
0205 
0206 function  imgs = get_prepackaged_fmdls( fmdl );
0207   switch fmdl
0208     case 'c=1;h=2;r=.08;ce=16;bg=1;st=1;me=1;nd'
0209       fmdl = ng_mk_cyl_models([2,1,0.18],[16,1],[0.05]); 
0210       fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0211       fmdl = mdl_normalize(fmdl,1);
0212       imgs= mk_image( fmdl, 1);
0213     otherwise
0214       error('specified fmdl (%s) is not understood', fmdl);
0215   end
0216 
0217 function [vi,vh,xyz,opt]= stim_targets(imgs, Nsim, opt );
0218     fmdl = imgs.fwd_model;
0219    ctr =  mean(fmdl.nodes);  
0220    maxx = max(abs(fmdl.nodes(:,1) - ctr(1)));
0221    maxy = max(abs(fmdl.nodes(:,2) - ctr(2)));
0222    if numel(opt.distr) > 1
0223       xyzr = opt.distr;
0224       xyzr(4,:) = calc_radius(mean([maxx,maxy]),opt,size(opt.distr,2));
0225    else
0226        switch opt.distr
0227            case 0 % original
0228                r = linspace(0,0.9, Nsim);
0229                th = r*4321; % want object to jump around in radius
0230                xyzr = [maxx*r.*cos(th); maxy*r.*sin(th);
0231                    opt.target_plane*ones(1,Nsim);
0232                    0.05*mean([maxx,maxy])*ones(1,Nsim)];
0233                
0234            case 1 %centre-heavy
0235                F = fourier_fit(opt.contour_boundary(:,1:2));
0236                v = linspace(0,1,Nsim*100+1); v(end)=[];
0237                pts = fourier_fit(F,v);
0238                idx_p = floor(rand(Nsim,1)*Nsim*100);
0239                xyzr = pts(idx_p,:)'.*repmat(rand(Nsim,1),[1 2])';
0240                xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0241                
0242                % TODO: What size is good here and how to figure it out?
0243                xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0244            case 2 %uniform
0245                %            F = fourier_fit(opt.contour_boundary(:,1:2));
0246                %            v = linspace(0,1,101); v(end)=[];
0247                %            pts = fourier_fit(F,v);
0248                pts = opt.contour_boundary(:,1:2);
0249                % avoid edges
0250                pts = 0.9*( pts - repmat(ctr(1:2),length(pts),1) ) + repmat(ctr(1:2),length(pts),1);
0251                % using maxx and maxy below would in general not produce a
0252                % uniform distribution
0253                lim = max(maxx, maxy);
0254                x = ctr(1) + (rand(Nsim*10,1)-0.5)*2*lim;
0255                y = ctr(2) + (rand(Nsim*10,1)-0.5)*2*lim;
0256                IN = inpolygon(x,y,pts(:,1),pts(:,2));
0257                xyzr(1,:) = x(find(IN,Nsim));
0258                xyzr(2,:) = y(find(IN,Nsim));
0259                xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0260                % TODO: What size is good here and how to figure it out?
0261                xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0262            case 3 % uniform, non-random
0263                %            F = fourier_fit(opt.elec_loc(:,1:2));
0264                %            v = linspace(0,1,101); v(end)=[];
0265                %            pts = fourier_fit(F,v);
0266                pts = opt.contour_boundary(:,1:2);
0267                lim = max(maxx, maxy);
0268                frac = polyarea(pts(:,1),pts(:,2)) / (2*lim)^2;
0269                [x,y] = ndgrid( linspace(-lim,lim,ceil(sqrt(Nsim/frac))), ...
0270                    linspace(-lim,lim,ceil(sqrt(Nsim/frac))));
0271                
0272                x = x+ctr(1); y = y + ctr(2);
0273                IN = inpolygon(x,y,pts(:,1),pts(:,2));
0274                xyzr(1,:) = x(find(IN));
0275                xyzr(2,:) = y(find(IN));
0276                xyzr(3,:) = calc_offset(opt.target_plane,opt,size(xyzr,2));
0277                % TODO: What size is good here and how to figure it out?
0278                xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,size(xyzr,2));
0279                eidors_msg(['mk_GREIT_model: Using ' num2str(size(xyzr,2)) ' points']);
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 = select_imdl( fmdl,{'Basic GN dif'});
0347    end
0348    
0349    
0350     function opt = parse_options(opt,fmdl,imdl, weight)
0351 
0352     if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end
0353     if ~isfield(opt, 'square_pixels')
0354         opt.square_pixels = 0;
0355     end
0356     % Allow imdl.rec_model to overwrite options.imgsz
0357     if isfield(imdl,'rec_model') && ~isempty(imdl.rec_model)
0358         % this assumes rec_model is a rectangular grid, as it should
0359         opt.imgsz(1) = numel(unique(imdl.rec_model.nodes(:,1)))-1;
0360         opt.imgsz(2) = numel(unique(imdl.rec_model.nodes(:,2)))-1;
0361         try
0362             opt.imgsz(3) = numel(unique(imdl.rec_model.nodes(:,3)))-1;
0363         end
0364     end  
0365     
0366     if ~isfield(opt, 'distr'),     opt.distr = 3;       end 
0367     if ~isfield(opt, 'Nsim' ),     opt.Nsim  = 1000;    end
0368     if ~isfield(opt, 'noise_figure'), opt.noise_figure = []; end
0369     if ~isfield(opt, 'image_SNR'), opt.image_SNR = []; end
0370     if isempty(opt.noise_figure) && isempty(opt.image_SNR) && isempty(weight)
0371         error('EIDORS:WrongInput', ...
0372             'The weight parameter must be specified if opt.noise_figure or opt.image_SNR are empty or absent');
0373     end
0374     if isfield(opt,'extra_noise')
0375       error('mk_GREIT_model: doesn''t currently support extra_noise');
0376     end
0377     if ~isfield(opt, 'target_size')
0378         opt.target_size = 0.05;
0379     end
0380     if sum(size(opt.target_size)) > 2
0381         if opt.target_size(1) == opt.target_size(2);
0382             opt.random_size = false;
0383         else
0384             opt.random_size = true;
0385         end
0386     end
0387     if sum(size(opt.target_size)) == 2
0388             opt.random_size = false;
0389     end
0390     
0391     % Calculate the position of the electrodes
0392     Nelecs = length(fmdl.electrode);
0393     for i=1:Nelecs
0394        enodesi = fmdl.electrode(i).nodes;
0395        elec_loc(i,:) = mean( fmdl.nodes( enodesi,:),1 );
0396     end
0397     opt.elec_loc = elec_loc;
0398     
0399     if ~isfield(opt, 'target_plane')
0400           opt.target_plane = mean(elec_loc(:,3));
0401     else
0402         t = opt.target_plane;
0403         minnode = min(fmdl.nodes);
0404         maxnode = max(fmdl.nodes);
0405         if t<minnode(3) || t>maxnode(3)
0406             warning('options.target_plane is outside the model!');
0407             eidors_msg('mk_GREIT_model: Resorting to default target_plane');
0408             opt.target_plane = mean(elec_loc(:,3));
0409         end
0410     end
0411     if ~isfield(opt, 'target_offset')
0412         opt.target_offset = 0;
0413     end
0414     if sum(size(opt.target_offset)) == 2
0415         if opt.target_offset < 0, opt.target_offset = 0; end
0416         opt.target_offset(2) = opt.target_offset(1);
0417     end
0418     if any(opt.target_offset > 0)
0419         opt.random_offset = true;
0420     else
0421         opt.random_offset = false;
0422     end
0423 
0424     if ~isfield(opt,'noise_figure_targets');
0425        R = max(max(fmdl.nodes(:,1:2)) - min(fmdl.nodes(:,1:2)));
0426        xyzr = mean(fmdl.nodes);
0427        xyzr(3) = opt.target_plane;
0428        xyzr(4) = mean(opt.target_size)*0.5*R;
0429        opt.noise_figure_targets = xyzr;
0430     end
0431 
0432        
0433     if ~isfield(opt,'keep_intermediate_results');
0434        opt.keep_intermediate_results = false;
0435     end
0436     
0437     
0438     try, opt.normalize = fmdl.normalize_measurements;
0439     catch, 
0440         opt.normalize = 0;
0441         eidors_msg('mk_GREIT_model: fmdl.normalize_measurements not specified, assuming 0');
0442     end
0443     
0444     % find the boundary at target level (needed in many places)
0445     slc = mdl_slice_mesher(fmdl,[inf inf opt.target_plane]);
0446     bnd = find_boundary(slc.fwd_model);
0447     opt.contour_boundary = order_loop(slc.fwd_model.nodes(unique(bnd),:));
0448     
0449 
0450 function do_unit_test
0451 
0452 sidx= 1; subplot(4,4,sidx);
0453  do_performance_test; 
0454 
0455 % Create a 3D elliptical cylinder with 16 circular electrodes
0456 fmdl_1= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1]); %show_fem(fmdl);
0457 % Put two balls into the elliptical cylinder
0458 extra={'ball','solid ball = sphere(0.5,0.5,0.5;0.1);'};
0459 [fmdl_2,mat_idx]= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1],extra); 
0460 % Set the model to use adjacent current patterns
0461 stim = mk_stim_patterns(16,1,[0,1],[0,1],{}); 
0462 fmdl_1.stimulation = stim;
0463 fmdl_2.stimulation = stim;
0464 % Simulate homogeneous voltages (background conductivity = 0.5);
0465 img = mk_image(fmdl_2, 0.5); vh = fwd_solve(img); %show_fem(img);
0466 % Simulate inhomogeneous voltages (ball conductivity = 1.0);
0467 img.elem_data(mat_idx{2})= 1.0; vi = fwd_solve(img); 
0468 sidx= sidx+1; subplot(4,4,sidx);
0469 show_fem(img);
0470 % Reconstruct the image using GREITv1
0471 imdl= mk_common_gridmdl('GREITc1'); 
0472 img= inv_solve(imdl,vh,vi);
0473 sidx= sidx+1; subplot(4,4,sidx);
0474 show_slices(img)
0475 
0476 % Create a GREIT model for the ellipse
0477 opt.noise_figure = 0.5; opt.distr = 3;opt.square_pixels = 1; %other options are defaults
0478 fmdl_2 = mdl_normalize(fmdl_2,0);
0479 % use the true model (inverse crime)
0480 img_2 = mk_image(fmdl_2,0.5);
0481 imdl1 = mk_GREIT_model(img_2, 0.25, [], opt);
0482 img1= inv_solve(imdl1,vh,vi);  
0483 sidx= sidx+1; subplot(4,4,sidx);
0484 show_slices(img1);
0485 
0486 % now do the same but using image SNR and not NF
0487 opt = rmfield(opt,'noise_figure');
0488 opt.image_SNR = 1e-3; 
0489 weight = 90; % need to choose a weight that works with SNR
0490 imdl1 = mk_GREIT_model(img_2, 0.25, weight, opt);
0491 img1= inv_solve(imdl1,vh,vi);  
0492 sidx= sidx+1; subplot(4,4,sidx); show_slices(img1);
0493 
0494 unit_test_cmp('Expect no PJT or X', ~isfield(imdl1.solve_use_matrix, 'PJt') & ...
0495                                     ~isfield(imdl1.solve_use_matrix, 'X'), true);
0496 
0497 weight = [];
0498 opt.keep_intermediate_results = true;
0499 imdl1 = mk_GREIT_model(img_2, 0.25, weight, opt);
0500 img1= inv_solve(imdl1,vh,vi);  
0501 sidx= sidx+1; subplot(4,4,sidx); show_slices(img1);
0502 
0503 unit_test_cmp('Expect PJT and X', isfield(imdl1.solve_use_matrix, 'PJt') & ...
0504                                   isfield(imdl1.solve_use_matrix, 'X'), true);
0505 
0506 opt = rmfield(opt,{'image_SNR', 'keep_intermediate_results'}); opt.noise_figure = 0.5;
0507 
0508 % use honogenous model
0509 fmdl_1 = mdl_normalize(fmdl_1,0);
0510 imdl2 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0511 img2= inv_solve(imdl2,vh,vi); 
0512 sidx= sidx+1; subplot(4,4,sidx); show_slices(img2);
0513 
0514 
0515 % specify targets for NF calc
0516 opt.noise_figure_targets = [-.5 0 .5 .2;.5 0 .5 .2;];
0517 imdl3 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0518 img3= inv_solve(imdl3,vh,vi); 
0519 sidx= sidx+1; subplot(4,4,sidx); show_slices(img3);
0520 % cleanup
0521 opt = rmfield(opt,'noise_figure_targets');
0522 
0523 
0524 %% repeat with normalized data
0525 fmdl_2 = mdl_normalize(fmdl_2,1);
0526 % use the true model (inverse crime)
0527 imdl3 = mk_GREIT_model(mk_image(fmdl_2,0.5), 0.25, [], opt);
0528 img3= inv_solve(imdl3,vh,vi); 
0529 
0530 % use honogenous model
0531 fmdl_1 = mdl_normalize(fmdl_1,1);
0532 imdl4 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0533 img4= inv_solve(imdl4,vh,vi); 
0534 
0535 sidx= sidx+1; subplot(4,4,sidx);
0536 show_slices([img1 img2 img3 img4])
0537 
0538 
0539 %% Use a prepackaged model
0540 fmdl = mk_library_model('adult_male_16el_lungs');
0541 fmdl.stimulation = stim;
0542 fmdl = mdl_normalize(fmdl,1);
0543 img = mk_image(fmdl,1);
0544 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.3;
0545 vh = fwd_solve(img);
0546 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.4;
0547 vi = fwd_solve(img);
0548 
0549 
0550 fmdl2 = mk_library_model('adult_male_16el');
0551 fmdl2.stimulation = stim;
0552 fmdl2 = mdl_normalize(fmdl2,1);
0553 
0554 opt.imgsz = [50 30];
0555 opt.square_pixels = 1;
0556 imdl = mk_GREIT_model(fmdl2,0.25,3,opt);
0557 
0558 img = inv_solve(imdl,vh, vi);
0559 sidx= sidx+1; subplot(4,4,sidx);
0560 show_slices(img);
0561 
0562 
0563 % do the same again with image SNR and not NF
0564 opt = rmfield(opt,'noise_figure');
0565 opt.image_SNR = 1e-4; 
0566 weight = 0.5; % need to choose a weight that works with SNR
0567 imdl = mk_GREIT_model(fmdl2, 0.25, weight, opt);
0568 img = inv_solve(imdl,vh, vi);
0569 sidx= sidx+1; subplot(4,4,sidx); show_slices(img);
0570 
0571 opt.image_SNR_targets = [0.3 0.3  0.5 0.05;  0.3 -0.3 0.5 0.05; ...
0572                          0.3 -0.3 0.5 0.05; -0.3 -0.3 0.5 0.05; ...
0573                          0.3 0    0.5 0.05; -0.3  0   0.5 0.05; ...
0574                          0   0.3  0.5 0.05;  0   -0.3 0.5 0.05]';
0575 opt.image_SNR = 3e-4; 
0576 weight = 1E-2; % 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 
0582 function do_performance_test
0583 % Reconstruct GREIT Images
0584 imdl_v1 = mk_common_gridmdl('GREITc1');
0585 imdl_v1.inv_solve.calc_solution_error = false;
0586 
0587 % Reconstruct backprojection Images
0588 imdl_bp = mk_common_gridmdl('backproj');
0589 
0590 % Recosntruct with new GREIT
0591 % fmdl = ng_mk_cyl_models([2,1,0.05],[16,1],[0.05]);
0592 fmdl = mk_library_model('cylinder_16x1el_fine');
0593 fmdl.nodes = fmdl.nodes/15; % make radius 1;
0594 fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
0595 opt.noise_figure = 0.88;
0596 opt.target_size = 0.1;
0597 opt.distr = 0;
0598 imdl_gr = mk_GREIT_model(fmdl, 0.2, [], opt);
0599 
0600 opt = struct();
0601 opt.noise_figure = 0.5; % current recommendation
0602 imdl_def = mk_GREIT_model(fmdl,0.2,[],opt);
0603 
0604 opt.desired_solution_fn = 'GREIT_desired_img_original';
0605 imdl_org = mk_GREIT_model(fmdl,0.2,[],opt);
0606 
0607 test_performance( { imdl_v1, imdl_gr, imdl_def, imdl_org},fmdl );

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005