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 %          - inv model (experimental), or
0013 %          - string specifying prepackaged models
0014 %
0015 %   radius - requested weighting matrix  (recommend 0.25 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 %     Nsim        - number of training points (default: 1000)
0022 %     distr       - distribution of training points:
0023 %         0 -> original (as per GREITv1, default)
0024 %         1 -> random, centre-heavy
0025 %         2 -> random, uniform
0026 %         3 -> fixed, uniform (debug)
0027 %     target_size - size of simulated targets as proportion of mesh radius
0028 %         (default: 0.02). Can be specified as [min_size max_size] for
0029 %         random variation
0030 %     target_plane - the (mean) height z at which simulation targets are
0031 %         placed. This controls the image plane. Default: mean electrode
0032 %         height
0033 %     target_offset - maximum allowed vertical displacement from the
0034 %         target_plane (default: 0). Can be specified as
0035 %         [down_offset up_offset].
0036 %     noise_figure - the noise figure (NF) to achieve. Overwrites weight
0037 %         which will be optimised to achieve the target NF.
0038 %     extra_noise - extra noise samples (such as electrode movement)
0039 %
0040 % NOTE
0041 %   currently extra_noise is not supported
0042 %   currently weighting matrix must be scalar
0043               
0044 % Examples
0045 %   imdl =  mk_GREIT_model( 'c=1;h=2;r=.08;ce=16;bg=1;st=1;me=1;nd', 0.25, 10);
0046 % OR
0047 %   fmdl = mk_library_model('adult_male_16el');
0048 %   fmdl.stimulation = stim;
0049 %   fmdl.normalize_measurements = 1;
0050 %   opt.noise_figure = 0.5;
0051 %   imdl = mk_GREIT_model(fmdl,0.25,5,opt);
0052 
0053 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0054 % $Id: mk_GREIT_model.html 2819 2011-09-07 16:43:11Z aadler $
0055 
0056 if isstr(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0057 
0058 if nargin < 4, options = [];end
0059 [imdl,fmdl,imgs] = parse_fmdl(fmdl);
0060 options = parse_options(options,fmdl,imdl);
0061 
0062 cache_obj= { fmdl, imdl, imgs, radius, weight, options};
0063 
0064 out= eidors_obj('get-cache', cache_obj, 'mk_GREIT_model');
0065 if ~isempty(out)
0066    eidors_msg('mk_GREIT_model: using cached value', 3);
0067    imdl= out{1}; weight=out{2}; return
0068 end
0069 [imdl, weight]= mk_GREIT_model_calc( fmdl, imdl, imgs, radius, weight, options);
0070 
0071 eidors_obj('set-cache', cache_obj, 'mk_GREIT_model', {imdl,weight});
0072 eidors_msg('mk_GREIT_model: setting cached value', 3);
0073 
0074 function [imdl, weight]= mk_GREIT_model_calc( fmdl, imdl, imgs, radius, weight, opt)
0075 
0076 Nsim = opt.Nsim;
0077 [vi,vh,xy,bound,elec_loc,opt]= stim_targets(imgs, Nsim, opt );
0078 
0079 
0080 xgrid = linspace(opt.minnode(1),opt.maxnode(1),opt.imgsz(1)+1);
0081 ygrid = linspace(opt.minnode(2),opt.maxnode(2),opt.imgsz(2)+1);
0082 rmdl = mk_grid_model([],xgrid,ygrid);
0083 x_avg = conv2(xgrid, [1,1]/2,'valid');
0084 y_avg = conv2(ygrid, [1,1]/2,'valid');
0085 [x,y] = ndgrid( x_avg, y_avg);
0086 
0087 %Calculate rec_model (if absent) and find the inside array
0088 if ~isfield(imdl,'rec_model');
0089  inside = inpolygon(x(:),y(:),bound(:,1),bound(:,2) );
0090  
0091  ff = find(~inside);
0092  rmdl.elems([2*ff, 2*ff-1],:)= [];
0093  rmdl.coarse2fine([2*ff, 2*ff-1],:)= [];
0094  rmdl.coarse2fine(:,ff)= [];
0095  
0096  imdl.rec_model = rmdl;
0097 else
0098  % this assumes the original grid model was created the same way
0099  inside = ismember(rmdl.elems,imdl.rec_model.elems,'rows');
0100  inside = inside(1:2:end);
0101 end
0102 
0103 imdl.solve = @solve_use_matrix;
0104 log_level = eidors_msg( 'log_level', 1);
0105 
0106 if ~isempty(opt.noise_figure)
0107     target = opt.noise_figure;
0108     if ~isempty(weight)
0109         eidors_msg('mk_GREIT_model: Using weight parameter as a guess, options.noise_figure is non-empty');
0110     else
0111         weight = target;
0112     end
0113 
0114     xyzr = mean(fmdl.nodes);
0115     xyzr(3) = opt.target_plane;
0116     xyzr(4) = opt.target_size;
0117     [jnk,vi_NF] = simulate_movement(imgs,xyzr');
0118     eidors_msg('mk_GREIT_model: Finding noise weighting for given Noise Figure',1);
0119     eidors_msg('mk_GREIT_model: This will take a while...',1);
0120     f = @(X) to_optimise(vh,vi,xy, radius, X, opt, inside, imdl, target, vi_NF);
0121     fms_opts.TolFun = 0.01*target; %don't need higher accuracy
0122     [weight, NF] = fminsearch(f, weight);
0123     eidors_msg(['mk_GREIT_model: Optimal solution gives NF=' ... 
0124         num2str(NF+target) ' with weight=' num2str(weight)],1);
0125 end
0126 eidors_msg( 'log_level', log_level);
0127 RM= calc_GREIT_RM(vh,vi, xy, radius, weight, opt );
0128 imdl.solve_use_matrix.RM = resize_if_reqd(RM,inside);
0129 %imdl.solve_use_matrix.map = inside;
0130 
0131 function out = to_optimise(vh,vi,xy,radius,weight, opt, inside, imdl, ...
0132     target,vi_NF)
0133 
0134    % calculate GREIT matrix as usual
0135    RM = calc_GREIT_RM(vh,vi,xy, radius, weight, opt);
0136    imdl.solve_use_matrix.RM = resize_if_reqd(RM,inside);
0137    NF = calc_noise_params(imdl,vh, vi_NF);
0138    eidors_msg(['NF = ', num2str(NF), ' weight = ', num2str(weight)],1);
0139    out = (NF - target)^2;
0140 %    out = (mean(NF) - target)^2 + std(NF);
0141    
0142 function  imgs = get_prepackaged_fmdls( fmdl );
0143   switch fmdl
0144     case 'c=1;h=2;r=.08;ce=16;bg=1;st=1;me=1;nd'
0145       fmdl = ng_mk_cyl_models([2,1,0.18],[16,1],[0.05]); 
0146       fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0147       fmdl.normalize_measurements = 1;
0148       imgs= mk_image( fmdl, 1);
0149     otherwise
0150       error('specified fmdl (%s) is not understood', fmdl);
0151   end
0152 
0153 function [vi,vh,xy,bound,elec_loc,opt]= stim_targets(imgs, Nsim, opt );
0154     fmdl = imgs.fwd_model;
0155    ctr =  mean(fmdl.nodes);  
0156    maxx = max(abs(fmdl.nodes(:,1) - ctr(1)));
0157    maxy = max(abs(fmdl.nodes(:,2) - ctr(2)));
0158 
0159    % Calculate the position of the electrodes
0160    Nelecs = length(imgs.fwd_model.electrode);
0161    for i=1:Nelecs
0162        enodesi =     imgs.fwd_model.electrode(i).nodes;
0163        elec_loc(i,:) = mean( imgs.fwd_model.nodes( enodesi,:),1 );
0164    end
0165    
0166    if opt.target_plane == 1i
0167        opt.target_plane = mean(elec_loc(:,3));
0168    end
0169    
0170    % calculate the boundary (for external use)
0171    F = fourier_fit(elec_loc(:,1:2));
0172    v = linspace(0,1,100+1); v(end)=[];
0173    bound = fourier_fit(F,v);
0174 
0175    
0176    switch opt.distr 
0177        case 0 % original
0178            r = linspace(0,0.9, Nsim);
0179            th = r*4321; % want object to jump around in radius
0180            xyzr = [maxx*r.*cos(th); maxy*r.*sin(th); 
0181                opt.target_plane*ones(1,Nsim);
0182                0.05/mean([maxx,maxy])*ones(1,Nsim)];
0183        
0184        case 1 %centre-heavy
0185            % Now, use elec_loc to figure out the shape. We can assume the obj is
0186            % extruded in z
0187            F = fourier_fit(elec_loc(:,1:2));
0188            v = linspace(0,1,Nsim*100+1); v(end)=[];
0189            pts = fourier_fit(F,v);
0190            idx_p = floor(rand(Nsim,1)*Nsim*100);
0191            xyzr = pts(idx_p,:)'.*repmat(rand(Nsim,1),[1 2])';
0192            xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0193            
0194            % TODO: What size is good here and how to figure it out?
0195            xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0196        case 2 %uniform
0197            F = fourier_fit(elec_loc(:,1:2));
0198            v = linspace(0,1,101); v(end)=[];
0199            pts = fourier_fit(F,v);
0200            % avoid edges
0201            pts = 0.9*( pts - repmat(ctr(1:2),length(pts),1) ) + repmat(ctr(1:2),length(pts),1);
0202            % using maxx and maxy below would in general not produce a
0203            % uniform distribution
0204            lim = max(maxx, maxy);
0205            x = ctr(1) + (rand(Nsim*10,1)-0.5)*2*lim;
0206            y = ctr(2) + (rand(Nsim*10,1)-0.5)*2*lim;
0207            IN = inpolygon(x,y,pts(:,1),pts(:,2));
0208            xyzr(1,:) = x(find(IN,Nsim));
0209            xyzr(2,:) = y(find(IN,Nsim));
0210            xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0211            % TODO: What size is good here and how to figure it out?
0212            xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0213        case 3 % uniform, non-random
0214            F = fourier_fit(elec_loc(:,1:2));
0215            v = linspace(0,1,101); v(end)=[];
0216            pts = fourier_fit(F,v);
0217            lim = max(maxx, maxy);
0218            frac = polyarea(pts(:,1),pts(:,2)) / (2*lim)^2;
0219            [x,y] = ndgrid( linspace(-lim,lim,ceil(sqrt(Nsim/frac))), ...
0220                            linspace(-lim,lim,ceil(sqrt(Nsim/frac))));
0221                       
0222            x = x+ctr(1); y = y + ctr(2);    
0223            IN = inpolygon(x,y,pts(:,1),pts(:,2));
0224            xyzr(1,:) = x(find(IN));
0225            xyzr(2,:) = y(find(IN));
0226            xyzr(3,:) = calc_offset(opt.target_plane,opt,size(xyzr,2));
0227            % TODO: What size is good here and how to figure it out?
0228            xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,size(xyzr,2));
0229            eidors_msg(['mk_GREIT_model: Using ' num2str(size(xyzr,2)) ' points']);
0230    end
0231    before = size(xyzr,2);
0232    [vh,vi,xyzr] = simulate_movement(imgs, xyzr);
0233    after = size(xyzr,2);
0234    if(after~=before)
0235        eidors_msg(['mk_GREIT_model: Now using ' num2str(after) ' points']);
0236    end
0237    xy = xyzr(1:2,:);
0238 
0239 function z = calc_offset(z0,opt,Nsim)
0240     if opt.random_offset
0241         l_bnd = opt.target_offset(1);
0242         width = sum(opt.target_offset(1:2));
0243         z = z0 - l_bnd + rand(Nsim,1)*width;
0244     else
0245         z = z0*ones(Nsim,1);
0246     end
0247 
0248 function r = calc_radius(R,opt,Nsim)
0249    if opt.random_size
0250        min_sz = opt.target_size(1);
0251        max_sz = opt.target_size(2);
0252        range = max_sz - min_sz;
0253        r = (min_sz + rand(Nsim,1)*range)*R;
0254    else
0255        r = opt.target_size(1)*ones(Nsim,1)*R;
0256    end
0257            
0258    
0259    
0260 function RM = resize_if_reqd(RM,inside);
0261    szRM = size(RM,1);
0262    if sum(inside) == szRM
0263       % RM is fine
0264    elseif size(inside,1) == szRM
0265       RM = RM(inside,:);
0266    else
0267       error('mismatch in size of provided RecMatrix');
0268    end
0269 
0270 
0271 function [imdl,fmdl,imgs] = parse_fmdl(fmdl);
0272    imdl = []; 
0273    if isstr(fmdl)
0274       imgs = get_prepackaged_fmdls( fmdl );
0275       fmdl = imgs.fwd_model;
0276    elseif isfield(fmdl,'type');
0277      switch fmdl.type
0278    %  if we get a fwd_model, assume uniform conductivity backgnd of 1
0279        case 'fwd_model'; imgs = mk_image( fmdl, 1);
0280    %  if we get an image, use it. It may have a non-uniform backgnd
0281        case 'image';     imgs = fmdl; % fmdl was an image
0282                          fmdl = imgs.fwd_model; % now it's a fmdl
0283        case 'inv_model'; imdl = fmdl;
0284                          fmdl = imdl.fwd_model;
0285                          imgs = mk_image( fmdl, 1);
0286        otherwise; error('unrecognized eidors object');
0287      end
0288    else
0289       error('specified parameter must be an object or a string');
0290    end
0291    % Prepare model
0292    if isempty(imdl)
0293       imdl = select_imdl( fmdl,{'Basic GN dif'});
0294    end
0295 
0296 function opt = parse_options(opt,fmdl,imdl);
0297     
0298     maxnode = max(fmdl.nodes); minnode = min(fmdl.nodes);
0299     opt.maxnode = maxnode;     opt.minnode = minnode; 
0300     
0301     if ~isfield(opt, 'imgsz'),     opt.imgsz = [32 32]; end
0302     % Allow imdl.rec_model to overwrite options.imgsz
0303     if isfield(imdl,'rec_model') && ~isempty(imdl.rec_model)
0304         % this assumes rec_model is a rectangular grid, as it should
0305         opt.imgsz(1) = numel(unique(imdl.rec_model.nodes(:,1)))-1;
0306         opt.imgsz(2) = numel(unique(imdl.rec_model.nodes(:,2)))-1;
0307     end  
0308     
0309     if ~isfield(opt, 'distr'),     opt.distr = 3;       end 
0310     if ~isfield(opt, 'Nsim' ),     opt.Nsim  = 1000;    end
0311     if ~isfield(opt, 'noise_figure'), opt.noise_figure = []; end
0312     if isfield(opt,'extra_noise')
0313       error('mk_GREIT_model: doesn''t currently support extra_noise');
0314     end
0315     if ~isfield(opt, 'target_size')
0316         opt.target_size = 0.02;
0317     end
0318     if sum(size(opt.target_size)) > 2
0319         if opt.target_size(1) == opt.target_size(2);
0320             opt.random_size = false;
0321         else
0322             opt.random_size = true;
0323         end
0324     end
0325     if sum(size(opt.target_size)) == 2
0326             opt.random_size = false;
0327     end
0328     
0329     if ~isfield(opt, 'target_plane')
0330         opt.target_plane = 1i;
0331     else
0332         t = opt.target_plane;
0333         if t<minnode(3) || t>maxnode(3)
0334             warning('options.target_plane is outside the model!');
0335             eidors_msg('mk_GREIT_model: Resorting to default target_plane');
0336             opt.target_plane = 1i;
0337         end
0338     end
0339     if ~isfield(opt, 'target_offset')
0340         opt.target_offset = 0;
0341     end
0342     if sum(size(opt.target_offset)) == 2
0343         if opt.target_offset < 0, opt.target_offset = 0; end
0344         opt.target_offset(2) = opt.target_offset(1);
0345     end
0346     if any(opt.target_offset > 0)
0347         opt.random_offset = true;
0348     else
0349         opt.random_offset = false;
0350     end
0351 
0352     try, opt.normalize = fmdl.normalize_measurements;
0353     catch, 
0354         opt.normalize = 0;
0355         eidors_msg('mk_GREIT_model: fmdl.normalize_measurements not specified, assuming 0');
0356     end
0357     opt.meshsz = [minnode(1) maxnode(1) minnode(2) maxnode(2)];
0358 
0359 function do_unit_test
0360 % Create a 3D elliptical cylinder with 16 circular electrodes
0361 fmdl_1= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1]); %show_fem(fmdl);
0362 % Put two balls into the elliptical cylinder
0363 extra={'ball','solid ball = sphere(0.5,0.5,0.5;0.1) or sphere(0.5,-0.5,0.5;0.1);'};
0364 [fmdl_2,mat_idx]= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1],extra); 
0365 % Set the model to use adjacent current patterns
0366 stim = mk_stim_patterns(16,1,[0,1],[0,1],{}); 
0367 fmdl_1.stimulation = stim;
0368 fmdl_2.stimulation = stim;
0369 % Simulate homogeneous voltages (background conductivity = 0.5);
0370 img = mk_image(fmdl_2, 0.5); vh = fwd_solve(img); %show_fem(img);
0371 % Simulate inhomogeneous voltages (ball conductivity = 1.0);
0372 img.elem_data(mat_idx{2})= 1.0; vi = fwd_solve(img); 
0373 show_fem(img);
0374 % Reconstruct the image using GREITv1
0375 imdl= mk_common_gridmdl('GREITc1'); 
0376 img= inv_solve(imdl,vh,vi);
0377 figure, show_slices(img)
0378 
0379 % Create a GREIT model for the ellipse
0380 opt.noise_figure = 0.5; opt.distr = 3; %other options are defaults
0381 fmdl_2.normalize_measurements = 0;
0382 % use the true model (inverse crime)
0383 imdl1 = mk_GREIT_model(mk_image(fmdl_2,0.5), 0.25, [], opt);
0384 img1= inv_solve(imdl1,vh,vi); 
0385 
0386 % use honogenous model
0387 fmdl_1.normalize_measurements = 0;
0388 imdl2 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0389 img2= inv_solve(imdl2,vh,vi); 
0390 
0391 %% repeat with normalized data
0392 fmdl_2.normalize_measurements = 1;
0393 % use the true model (inverse crime)
0394 imdl3 = mk_GREIT_model(mk_image(fmdl_2,0.5), 0.25, [], opt);
0395 img3= inv_solve(imdl3,vh,vi); 
0396 
0397 % use honogenous model
0398 fmdl_1.normalize_measurements = 1;
0399 imdl4 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0400 img4= inv_solve(imdl4,vh,vi); 
0401 
0402 figure
0403 show_slices([img1 img2 img3 img4])
0404 
0405 
0406 %% Use a prepackaged model
0407 fmdl = mk_library_model('adult_male_16el_lungs');
0408 fmdl.stimulation = stim;
0409 fmdl.normalize_measurements = 1;
0410 img = mk_image(fmdl,1);
0411 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.3;
0412 vh = fwd_solve(img);
0413 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.4;
0414 vi = fwd_solve(img);
0415 
0416 
0417 fmdl2 = mk_library_model('adult_male_16el');
0418 fmdl2.stimulation = stim;
0419 fmdl2.normalize_measurements = 1;
0420 imdl = mk_GREIT_model(fmdl2,0.25,3,opt);
0421 
0422 img = inv_solve(imdl,vh, vi);
0423 figure
0424 show_slices(img);
0425 
0426 
0427
mk_GREIT_model

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
