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.
     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.

 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 %     extra_noise - extra noise samples (such as electrode movement)
0046 %     desired_solution_fn - specify a function to calculate the desired
0047 %         image. It must have the signature:
0048 %         D = my_function( xyc, radius, options);
0049 %         See CALC_GREIT_RM for details.
0050 %
0051 % NOTE
0052 %   currently extra_noise is not supported
0053 %   currently weighting matrix must be scalar
0054 %
0055 % Examples
0056 %   fmdl = mk_library_model('adult_male_16el');
0057 %   fmdl.stimulation = mk_stim_patterns(16,1,'{ad}','{ad}',{},1);
0058 %   fmdl.normalize_measurements = 1;
0059 %   imdl = mk_GREIT_model(fmdl,0.25,5); % uses weight 5
0060 %   OR
0061 %   opt.noise_figure = 0.5;
0062 %   imdl = mk_GREIT_model(fmdl,0.25,[],opt); % optimises weight for NF=0.5
0063 %
0064 % CITATION_REQUEST:
0065 % AUTHOR: A Adler et al.
0066 % TITLE: GREIT: a unified approach to 2D linear EIT reconstruction of lung
0067 % images
0068 % JOURNAL: Phys Meas
0069 % YEAR: 2009
0070 % VOL: 30
0071 % NUM: 6
0072 % PAGE: S35-55
0073 % LINK: http://iopscience.iop.org/0967-3334/30/6/S03
0074 %
0075 % See also CALC_GREIT_RM
0076 
0077 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0078 % $Id: mk_GREIT_model.m 4886 2015-04-17 08:12:43Z bgrychtol-ipa $
0079 
0080 citeme(mfilename);
0081 
0082 if isstr(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0083 
0084 if nargin < 4, options = [];end
0085 [imdl,fmdl,imgs] = parse_fmdl(fmdl);
0086 options = parse_options(options,fmdl,imdl, weight);
0087 
0088 copt.cache_obj= { fmdl, imdl, imgs, radius, weight, options};
0089 copt.fstr = 'mk_GREIT_model';
0090 params = {fmdl, imdl, imgs, radius, weight, options};
0091 
0092 [imdl, weight] = eidors_cache(@mk_GREIT_model_calc, params, copt);
0093 
0094 
0095 function [imdl, weight]= mk_GREIT_model_calc( fmdl, imdl, imgs, radius, weight, opt)
0096 
0097 Nsim = opt.Nsim;
0098 [vi,vh,xyz,opt]= stim_targets(imgs, Nsim, opt );
0099 
0100 %Calculate rec_model (if absent)
0101 if ~isfield(imdl,'rec_model');
0102 %    opt.do_coarse2fine = 0;
0103    [imdl.rec_model imdl.fwd_model] = mk_pixel_slice(fmdl,opt.target_plane,opt);
0104    imdl.rec_model.nodes(:,3) = []; % the third dimension complicated display
0105    % medical orientation
0106    imdl.rec_model.mdl_slice_mapper.y_pts = fliplr(imdl.rec_model.mdl_slice_mapper.y_pts);
0107 end
0108 
0109 opt.rec_model = imdl.rec_model; % for desired image calculation
0110 
0111 imdl.solve = @solve_use_matrix;
0112 %
0113 
0114 if ~isempty(opt.noise_figure)
0115     target = opt.noise_figure;
0116     if ~isempty(weight)
0117         eidors_msg('mk_GREIT_model: Using weight parameter as a guess, options.noise_figure is non-empty');
0118     else
0119         weight = target;
0120     end
0121     
0122     xyzr = opt.noise_figure_targets;
0123     [jnk,vi_NF] = simulate_movement(imgs,xyzr');
0124     vi_NF = sum(vi_NF,2); % sum the targets
0125     eidors_msg('mk_GREIT_model: Finding noise weighting for given Noise Figure',1);
0126     eidors_msg('mk_GREIT_model: This will take a while...',1);
0127     f = @(X) to_optimise(vh,vi,xyz, radius, X, opt, imdl, target, vi_NF);
0128     fms_opts.TolFun = 0.01*target; %don't need higher accuracy
0129     % The first call can take a long time. Take it out of the loop to
0130     % allow progress messages.
0131     imdl.solve_use_matrix.RM = calc_GREIT_RM(vh,vi,xyz, radius, weight, opt);
0132     log_level = eidors_msg( 'log_level', 1); % suppress messages
0133     if exist('OCTAVE_VERSION')
0134        % octave doesn't currently (2013 Apr) include an fminsearch function
0135        [weight, NF] = fminsearch_octave(f, weight,fms_opts);
0136     else
0137        [weight, NF] = fminsearch(f, weight,fms_opts);
0138     end
0139     eidors_msg(['mk_GREIT_model: Optimal solution gives NF=' ... 
0140         num2str(NF+target) ' with weight=' num2str(weight)],1);
0141      eidors_msg( 'log_level', log_level); % restore
0142 end
0143 %
0144 imdl.solve_use_matrix.RM= calc_GREIT_RM(vh,vi, xyz, radius, weight, opt );
0145 % imdl.solve_use_matrix.RM = resize_if_reqd(RM,inside,imdl.rec_model);
0146 imdl.jacobian_bkgnd = imgs;
0147 %imdl.solve_use_matrix.map = inside;
0148 
0149 function out = to_optimise(vh,vi,xy,radius,weight, opt, imdl, ...
0150     target,vi_NF)
0151 
0152    % calculate GREIT matrix as usual
0153    imdl.solve_use_matrix.RM = calc_GREIT_RM(vh,vi,xy, radius, weight, opt);
0154 %    imdl.solve_use_matrix.RM = resize_if_reqd(RM,inside,imdl.rec_model);
0155    NF = calc_noise_figure(imdl,vh, vi_NF);
0156    eidors_msg(['NF = ', num2str(NF), ' weight = ', num2str(weight)],1);
0157    out = (NF - target)^2;
0158 %    out = (mean(NF) - target)^2 + std(NF);
0159 
0160 
0161 function  imgs = get_prepackaged_fmdls( fmdl );
0162   switch fmdl
0163     case 'c=1;h=2;r=.08;ce=16;bg=1;st=1;me=1;nd'
0164       fmdl = ng_mk_cyl_models([2,1,0.18],[16,1],[0.05]); 
0165       fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0166       fmdl = mdl_normalize(fmdl,1);
0167       imgs= mk_image( fmdl, 1);
0168     otherwise
0169       error('specified fmdl (%s) is not understood', fmdl);
0170   end
0171 
0172 function [vi,vh,xyz,opt]= stim_targets(imgs, Nsim, opt );
0173     fmdl = imgs.fwd_model;
0174    ctr =  mean(fmdl.nodes);  
0175    maxx = max(abs(fmdl.nodes(:,1) - ctr(1)));
0176    maxy = max(abs(fmdl.nodes(:,2) - ctr(2)));
0177    if numel(opt.distr) > 1
0178       xyzr = opt.distr;
0179       xyzr(4,:) = calc_radius(mean([maxx,maxy]),opt,size(opt.distr,2));
0180    else
0181        switch opt.distr
0182            case 0 % original
0183                r = linspace(0,0.9, Nsim);
0184                th = r*4321; % want object to jump around in radius
0185                xyzr = [maxx*r.*cos(th); maxy*r.*sin(th);
0186                    opt.target_plane*ones(1,Nsim);
0187                    0.05*mean([maxx,maxy])*ones(1,Nsim)];
0188                
0189            case 1 %centre-heavy
0190                F = fourier_fit(opt.contour_boundary(:,1:2));
0191                v = linspace(0,1,Nsim*100+1); v(end)=[];
0192                pts = fourier_fit(F,v);
0193                idx_p = floor(rand(Nsim,1)*Nsim*100);
0194                xyzr = pts(idx_p,:)'.*repmat(rand(Nsim,1),[1 2])';
0195                xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0196                
0197                % TODO: What size is good here and how to figure it out?
0198                xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0199            case 2 %uniform
0200                %            F = fourier_fit(opt.contour_boundary(:,1:2));
0201                %            v = linspace(0,1,101); v(end)=[];
0202                %            pts = fourier_fit(F,v);
0203                pts = opt.contour_boundary(:,1:2);
0204                % avoid edges
0205                pts = 0.9*( pts - repmat(ctr(1:2),length(pts),1) ) + repmat(ctr(1:2),length(pts),1);
0206                % using maxx and maxy below would in general not produce a
0207                % uniform distribution
0208                lim = max(maxx, maxy);
0209                x = ctr(1) + (rand(Nsim*10,1)-0.5)*2*lim;
0210                y = ctr(2) + (rand(Nsim*10,1)-0.5)*2*lim;
0211                IN = inpolygon(x,y,pts(:,1),pts(:,2));
0212                xyzr(1,:) = x(find(IN,Nsim));
0213                xyzr(2,:) = y(find(IN,Nsim));
0214                xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0215                % TODO: What size is good here and how to figure it out?
0216                xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0217            case 3 % uniform, non-random
0218                %            F = fourier_fit(opt.elec_loc(:,1:2));
0219                %            v = linspace(0,1,101); v(end)=[];
0220                %            pts = fourier_fit(F,v);
0221                pts = opt.contour_boundary(:,1:2);
0222                lim = max(maxx, maxy);
0223                frac = polyarea(pts(:,1),pts(:,2)) / (2*lim)^2;
0224                [x,y] = ndgrid( linspace(-lim,lim,ceil(sqrt(Nsim/frac))), ...
0225                    linspace(-lim,lim,ceil(sqrt(Nsim/frac))));
0226                
0227                x = x+ctr(1); y = y + ctr(2);
0228                IN = inpolygon(x,y,pts(:,1),pts(:,2));
0229                xyzr(1,:) = x(find(IN));
0230                xyzr(2,:) = y(find(IN));
0231                xyzr(3,:) = calc_offset(opt.target_plane,opt,size(xyzr,2));
0232                % TODO: What size is good here and how to figure it out?
0233                xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,size(xyzr,2));
0234                eidors_msg(['mk_GREIT_model: Using ' num2str(size(xyzr,2)) ' points']);
0235        end
0236    end
0237    before = size(xyzr,2);
0238    [vh,vi,xyzr] = simulate_movement(imgs, xyzr);
0239    after = size(xyzr,2);
0240    if(after~=before)
0241        eidors_msg(['mk_GREIT_model: Now using ' num2str(after) ' points']);
0242    end
0243    xyz = xyzr(1:3,:);
0244 
0245 function z = calc_offset(z0,opt,Nsim)
0246     if opt.random_offset
0247         l_bnd = opt.target_offset(1);
0248         width = sum(opt.target_offset(1:2));
0249         z = z0 - l_bnd + rand(Nsim,1)*width;
0250     else
0251         z = z0*ones(Nsim,1);
0252     end
0253 
0254 function r = calc_radius(R,opt,Nsim)
0255    if opt.random_size
0256        min_sz = opt.target_size(1);
0257        max_sz = opt.target_size(2);
0258        range = max_sz - min_sz;
0259        r = (min_sz + rand(Nsim,1)*range)*R;
0260    else
0261        r = opt.target_size(1)*ones(Nsim,1)*R;
0262    end
0263            
0264    
0265    
0266 function RM = resize_if_reqd(RM,inside,rmdl)
0267    szRM = size(RM,1);
0268    if sum(inside) == szRM || ...
0269         szRM == size(rmdl.elems,1) || ...
0270         (isfield(rmdl,'coarse2fine') && szRM == size(rmdl.coarse2fine,2))
0271       % RM is fine
0272    elseif any(size(inside)==szRM) && any(size(inside) == 1)
0273       RM = RM(inside,:);
0274    else
0275       error('mismatch in size of provided RecMatrix');
0276    end
0277 
0278 
0279 function [imdl,fmdl,imgs] = parse_fmdl(fmdl);
0280    imdl = []; 
0281    if isstr(fmdl)
0282       imgs = get_prepackaged_fmdls( fmdl );
0283       fmdl = imgs.fwd_model;
0284    elseif isfield(fmdl,'type');
0285      switch fmdl.type
0286    %  if we get a fwd_model, assume uniform conductivity backgnd of 1
0287        case 'fwd_model'; imgs = mk_image( fmdl, 1);
0288    %  if we get an image, use it. It may have a non-uniform backgnd
0289        case 'image';     imgs = fmdl; % fmdl was an image
0290                          fmdl = imgs.fwd_model; % now it's a fmdl
0291        case 'inv_model'; imdl = fmdl;
0292                          fmdl = imdl.fwd_model;
0293                          imgs = calc_jacobian_bkgnd(imdl);
0294        otherwise; error('unrecognized eidors object');
0295      end
0296    else
0297       error('specified parameter must be an object or a string');
0298    end
0299    % Prepare model
0300    if isempty(imdl)
0301       imdl = select_imdl( fmdl,{'Basic GN dif'});
0302    end
0303    
0304    
0305     function opt = parse_options(opt,fmdl,imdl, weight)
0306 
0307     if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end
0308     if ~isfield(opt, 'square_pixels')
0309         opt.square_pixels = 0;
0310     end
0311     % Allow imdl.rec_model to overwrite options.imgsz
0312     if isfield(imdl,'rec_model') && ~isempty(imdl.rec_model)
0313         % this assumes rec_model is a rectangular grid, as it should
0314         opt.imgsz(1) = numel(unique(imdl.rec_model.nodes(:,1)))-1;
0315         opt.imgsz(2) = numel(unique(imdl.rec_model.nodes(:,2)))-1;
0316         try
0317             opt.imgsz(3) = numel(unique(imdl.rec_model.nodes(:,3)))-1;
0318         end
0319     end  
0320     
0321     if ~isfield(opt, 'distr'),     opt.distr = 3;       end 
0322     if ~isfield(opt, 'Nsim' ),     opt.Nsim  = 1000;    end
0323     if ~isfield(opt, 'noise_figure'), opt.noise_figure = []; end
0324     if isempty(opt.noise_figure) && isempty(weight)
0325         error('EIDORS:WrongInput', ...
0326             'The weight parameter must be specified if opt.noise_figure is empty or absent');
0327     end
0328     if isfield(opt,'extra_noise')
0329       error('mk_GREIT_model: doesn''t currently support extra_noise');
0330     end
0331     if ~isfield(opt, 'target_size')
0332         opt.target_size = 0.05;
0333     end
0334     if sum(size(opt.target_size)) > 2
0335         if opt.target_size(1) == opt.target_size(2);
0336             opt.random_size = false;
0337         else
0338             opt.random_size = true;
0339         end
0340     end
0341     if sum(size(opt.target_size)) == 2
0342             opt.random_size = false;
0343     end
0344     
0345     % Calculate the position of the electrodes
0346     Nelecs = length(fmdl.electrode);
0347     for i=1:Nelecs
0348        enodesi = fmdl.electrode(i).nodes;
0349        elec_loc(i,:) = mean( fmdl.nodes( enodesi,:),1 );
0350     end
0351     opt.elec_loc = elec_loc;
0352     
0353     if ~isfield(opt, 'target_plane')
0354           opt.target_plane = mean(elec_loc(:,3));
0355     else
0356         t = opt.target_plane;
0357         minnode = min(fmdl.nodes);
0358         maxnode = max(fmdl.nodes);
0359         if t<minnode(3) || t>maxnode(3)
0360             warning('options.target_plane is outside the model!');
0361             eidors_msg('mk_GREIT_model: Resorting to default target_plane');
0362             opt.target_plane = mean(elec_loc(:,3));
0363         end
0364     end
0365     if ~isfield(opt, 'target_offset')
0366         opt.target_offset = 0;
0367     end
0368     if sum(size(opt.target_offset)) == 2
0369         if opt.target_offset < 0, opt.target_offset = 0; end
0370         opt.target_offset(2) = opt.target_offset(1);
0371     end
0372     if any(opt.target_offset > 0)
0373         opt.random_offset = true;
0374     else
0375         opt.random_offset = false;
0376     end
0377 
0378     if ~isfield(opt,'noise_figure_targets');
0379        R = max(max(fmdl.nodes(:,1:2)) - min(fmdl.nodes(:,1:2)));
0380        xyzr = mean(fmdl.nodes);
0381        xyzr(3) = opt.target_plane;
0382        xyzr(4) = mean(opt.target_size)*0.5*R;
0383        opt.noise_figure_targets = xyzr;
0384     end
0385 
0386        
0387     
0388     
0389     try, opt.normalize = fmdl.normalize_measurements;
0390     catch, 
0391         opt.normalize = 0;
0392         eidors_msg('mk_GREIT_model: fmdl.normalize_measurements not specified, assuming 0');
0393     end
0394     
0395     % find the boundary at target level (needed in many places)
0396     slc = mdl_slice_mesher(fmdl,[inf inf opt.target_plane]);
0397     bnd = find_boundary(slc.fwd_model);
0398     opt.contour_boundary = order_loop(slc.fwd_model.nodes(unique(bnd),:));
0399     
0400 
0401 function do_unit_test
0402 
0403 do_performance_test; 
0404 % return;
0405 figure
0406 % Create a 3D elliptical cylinder with 16 circular electrodes
0407 fmdl_1= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1]); %show_fem(fmdl);
0408 % Put two balls into the elliptical cylinder
0409 extra={'ball','solid ball = sphere(0.5,0.5,0.5;0.1);'};
0410 [fmdl_2,mat_idx]= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1],extra); 
0411 % Set the model to use adjacent current patterns
0412 stim = mk_stim_patterns(16,1,[0,1],[0,1],{}); 
0413 fmdl_1.stimulation = stim;
0414 fmdl_2.stimulation = stim;
0415 % Simulate homogeneous voltages (background conductivity = 0.5);
0416 img = mk_image(fmdl_2, 0.5); vh = fwd_solve(img); %show_fem(img);
0417 % Simulate inhomogeneous voltages (ball conductivity = 1.0);
0418 img.elem_data(mat_idx{2})= 1.0; vi = fwd_solve(img); 
0419 show_fem(img);
0420 % Reconstruct the image using GREITv1
0421 imdl= mk_common_gridmdl('GREITc1'); 
0422 img= inv_solve(imdl,vh,vi);
0423 figure, subplot(2,2,1);
0424 show_slices(img)
0425 
0426 % Create a GREIT model for the ellipse
0427 opt.noise_figure = 0.5; opt.distr = 3;opt.square_pixels = 1; %other options are defaults
0428 fmdl_2 = mdl_normalize(fmdl_2,0);
0429 % use the true model (inverse crime)
0430 imdl1 = mk_GREIT_model(mk_image(fmdl_2,0.5), 0.25, [], opt);
0431 img1= inv_solve(imdl1,vh,vi);  
0432 subplot(2,2,2);show_slices(img1);
0433 
0434 % use honogenous model
0435 fmdl_1 = mdl_normalize(fmdl_1,0);
0436 imdl2 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0437 img2= inv_solve(imdl2,vh,vi); 
0438 subplot(2,2,3); show_slices(img2);
0439 
0440 
0441 % specify targets for NF calc
0442 opt.noise_figure_targets = [-.5 0 .5 .2;.5 0 .5 .2;];
0443 imdl3 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0444 img3= inv_solve(imdl3,vh,vi); 
0445 subplot(2,2,4); show_slices(img3);
0446 % cleanup
0447 opt = rmfield(opt,'noise_figure_targets');
0448 
0449 
0450 %% repeat with normalized data
0451 fmdl_2 = mdl_normalize(fmdl_2,1);
0452 % use the true model (inverse crime)
0453 imdl3 = mk_GREIT_model(mk_image(fmdl_2,0.5), 0.25, [], opt);
0454 img3= inv_solve(imdl3,vh,vi); 
0455 
0456 % use honogenous model
0457 fmdl_1 = mdl_normalize(fmdl_1,1);
0458 imdl4 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0459 img4= inv_solve(imdl4,vh,vi); 
0460 
0461 figure
0462 show_slices([img1 img2 img3 img4])
0463 
0464 
0465 %% Use a prepackaged model
0466 fmdl = mk_library_model('adult_male_16el_lungs');
0467 fmdl.stimulation = stim;
0468 fmdl = mdl_normalize(fmdl,1);
0469 img = mk_image(fmdl,1);
0470 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.3;
0471 vh = fwd_solve(img);
0472 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.4;
0473 vi = fwd_solve(img);
0474 
0475 
0476 fmdl2 = mk_library_model('adult_male_16el');
0477 fmdl2.stimulation = stim;
0478 fmdl2 = mdl_normalize(fmdl2,1);
0479 
0480 opt.imgsz = [50 30];
0481 opt.square_pixels = 1;
0482 imdl = mk_GREIT_model(fmdl2,0.25,3,opt);
0483 
0484 img = inv_solve(imdl,vh, vi);
0485 figure
0486 show_slices(img);
0487 
0488 
0489 function do_performance_test
0490 % Reconstruct GREIT Images
0491 imdl_v1 = mk_common_gridmdl('GREITc1');
0492 imdl_v1.inv_solve.calc_solution_error = false;
0493 
0494 % Reconstruct backprojection Images
0495 imdl_bp = mk_common_gridmdl('backproj');
0496 
0497 % Recosntruct with new GREIT
0498 % fmdl = ng_mk_cyl_models([2,1,0.05],[16,1],[0.05]);
0499 fmdl = mk_library_model('cylinder_16x1el_fine');
0500 fmdl.nodes = fmdl.nodes/15; % make radius 1;
0501 fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
0502 opt.noise_figure = 0.88;
0503 opt.target_size = 0.1;
0504 opt.distr = 0;
0505 imdl_gr = mk_GREIT_model(fmdl, 0.2, [], opt);
0506 
0507 opt = struct();
0508 opt.noise_figure = 0.5; % current recommendation
0509 imdl_def = mk_GREIT_model(fmdl,0.2,[],opt);
0510 
0511 opt.desired_solution_fn = 'GREIT_desired_img_original';
0512 imdl_org = mk_GREIT_model(fmdl,0.2,[],opt);
0513 
0514 test_performance( { imdl_v1, imdl_gr, imdl_def, imdl_org},fmdl );

Generated on Tue 12-May-2015 16:00:42 by m2html © 2005