mk_common_gridmdl

 MK_COMMON_MODEL: make EIT on reconstruction grids (GREIT)

 Usage
      inv_mdl = mk_common_gridmdl( mdl_string, RecMtx )

 2D models
   mk_common_gridmdl('b2c', RM)  - 32x32 circular shape, 16 elec
   mk_common_gridmdl('b2d', RM)  - 32x32 diamond shape, 16 elec
   mk_common_gridmdl('b2t?', RM) - 32x32 thorax shape, 16 elec
             - thorax levels 1-5 are provided

 Note that the electrodes added to the model are just to 
   indicate location, it does not necessarily correspond to the
   Reconstruction Matrix RM provided. 

 GREIT MODELS - calculated for library models
  mk_common_gridmdl('GREIT', PARAMS TO MK_LIBRARY_MODEL)
    (run mk_library_model('list') to get list)
  eg. mk_common_gridmdl('GREIT','cylinder_16x1el_coarse')
  eg. mk_common_gridmdl('GREIT','adult_male_32el')

 GREIT V1 (GREIT matrix for reconstruction to circle, 2009)
   mk_common_gridmdl('GREITc1') - 32x32 with Circle shape

 Sheffield Backprojection
   mk_common_gridmdl('b2d','backproj') - 32x32 with Diamond shape
   mk_common_gridmdl('b2c','backproj') - 32x32 with Circular shape
   mk_common_gridmdl('backproj') - 32x32 with Circular shape

 COPYRIGHT NOTICE FOR BACKPROJECTION MATRIX:
   This matrix is copyright DC Barber and BH Brown at
   University of Sheffield. It may be used free of
   charge for research and non-commercial purposes.
   Commercial applications require a licence from the
   University of Sheffield.

0001 function inv_mdl= mk_common_gridmdl( str, RM)
0002 % MK_COMMON_MODEL: make EIT on reconstruction grids (GREIT)
0003 %
0004 % Usage
0005 %      inv_mdl = mk_common_gridmdl( mdl_string, RecMtx )
0006 %
0007 % 2D models
0008 %   mk_common_gridmdl('b2c', RM)  - 32x32 circular shape, 16 elec
0009 %   mk_common_gridmdl('b2d', RM)  - 32x32 diamond shape, 16 elec
0010 %   mk_common_gridmdl('b2t?', RM) - 32x32 thorax shape, 16 elec
0011 %             - thorax levels 1-5 are provided
0012 %
0013 % Note that the electrodes added to the model are just to
0014 %   indicate location, it does not necessarily correspond to the
0015 %   Reconstruction Matrix RM provided.
0016 %
0017 % GREIT MODELS - calculated for library models
0018 %  mk_common_gridmdl('GREIT', PARAMS TO MK_LIBRARY_MODEL)
0019 %    (run mk_library_model('list') to get list)
0020 %  eg. mk_common_gridmdl('GREIT','cylinder_16x1el_coarse')
0021 %  eg. mk_common_gridmdl('GREIT','adult_male_32el')
0022 %
0023 % GREIT V1 (GREIT matrix for reconstruction to circle, 2009)
0024 %   mk_common_gridmdl('GREITc1') - 32x32 with Circle shape
0025 %
0026 % Sheffield Backprojection
0027 %   mk_common_gridmdl('b2d','backproj') - 32x32 with Diamond shape
0028 %   mk_common_gridmdl('b2c','backproj') - 32x32 with Circular shape
0029 %   mk_common_gridmdl('backproj') - 32x32 with Circular shape
0030 %
0031 % COPYRIGHT NOTICE FOR BACKPROJECTION MATRIX:
0032 %   This matrix is copyright DC Barber and BH Brown at
0033 %   University of Sheffield. It may be used free of
0034 %   charge for research and non-commercial purposes.
0035 %   Commercial applications require a licence from the
0036 %   University of Sheffield.
0037 %
0038 
0039 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0040 % $Id: mk_common_gridmdl.m 5112 2015-06-14 13:00:41Z aadler $
0041 
0042 name = str;
0043 if strcmp(str,'backproj')
0044    str= 'b2c';
0045    RM= 'backproj';
0046 elseif strcmp(str,'GREITc1');
0047    str= 'b2c';
0048    RM= 'GREITc1';
0049 else
0050    name = str;
0051 end
0052 
0053 n_elec= 16;
0054 
0055       space = linspace( -1, 1, 32+1 );
0056       fmdl= mk_grid_model( [], space, space);
0057       space_avg = conv2(space, [1,1]/2,'valid'); % average of each box
0058       [x,y] = ndgrid( space_avg, space_avg);
0059 switch str
0060    case 'b2c'
0061       inside=  (x(:).^2 + y(:).^2)<1.10 ;
0062 
0063    case 'b2d'
0064       inside=  (abs(x(:)) + abs(y(:)))<1.55 ;
0065 
0066    case 'b2t1'; inside = inside_thorax(x,y,1);
0067    case 'b2t2'; inside = inside_thorax(x,y,2);
0068    case 'b2t3'; inside = inside_thorax(x,y,3);
0069    case 'b2t4'; inside = inside_thorax(x,y,4);
0070    case 'b2t5'; inside = inside_thorax(x,y,5);
0071    
0072    case 'GREIT'; inv_mdl = greit_library_model(RM); return
0073 
0074    otherwise
0075       error(['mdl_string ',str,' not understood']);
0076 end
0077 
0078 if ischar(RM)
0079    switch RM
0080      case 'backproj'; RM= get_Sheffield_Backproj;
0081      case 'GREITc1';  RM= get_GREIT_c1;
0082      otherwise;       error('RM string not understood');
0083    end 
0084 end
0085       ff = find(~inside);
0086       fmdl.elems([2*ff, 2*ff-1],:)= [];
0087       fmdl.coarse2fine([2*ff, 2*ff-1],:)= [];
0088       fmdl.coarse2fine(:,ff)= [];
0089 
0090 fmdl.electrode = mk_electrode_locns( fmdl.nodes, n_elec );
0091 
0092 inv_mdl = eidors_obj('inv_model',['mk_common_gridmdl: ',name]);
0093 inv_mdl.reconst_type= 'difference';
0094 inv_mdl.fwd_model= fmdl;
0095 inv_mdl.solve_use_matrix.RM = resize_if_reqd(RM,inside);
0096 % solve_use_matrix has a c2f mapping field map, which
0097 % it may be useful to populate in this case
0098 %  inv_mdl.solve_use_matrix.map = map;
0099 inv_mdl.solve = @solve_use_matrix;
0100 inv_mdl.fwd_model = mdl_normalize(inv_mdl.fwd_model, 1);
0101 [st, els]= mk_stim_patterns(16, 1, '{ad}','{ad}', {}, 10);
0102 inv_mdl.fwd_model.meas_select= els;
0103 
0104 function inside = inside_thorax(x,y,level);
0105    [x_bdy, y_bdy ] = thorax_geometry(level,1); % normalized
0106    inside = inpolygon(x(:), y(:), x_bdy, y_bdy); 
0107 
0108 function RM = resize_if_reqd(RM,inside);
0109    szRM = size(RM,1);
0110    if sum(inside) == szRM
0111       % RM is fine
0112    elseif size(inside,1) == szRM
0113       RM = RM(inside,:);
0114    else
0115       error('mismatch in size of provided RecMatrix');
0116    end
0117 
0118 function RM = get_Sheffield_Backproj
0119    load Sheffield_Backproj_Matrix.mat
0120    RM = unpack_matrix(Sheffield_Backproj_Matrix);
0121 
0122 function RM= get_GREIT_c1;
0123    load GREIT_v10_Circ_Matrix.mat
0124    RM = unpack_matrix(GREIT_v10_Circ_Matrix);
0125 
0126 
0127 function RM= unpack_matrix(MM)
0128    RM = unpack_reconst_matrix(MM, 16, 32);
0129 
0130 
0131 function elec = mk_electrode_locns( nodes, n_elec );
0132    phi = linspace(0, 2*pi, n_elec+1); 
0133    phi(end) = [];
0134    rad = 1;
0135 
0136    for i= 1:n_elec
0137       posn_x = rad*sin(phi(i));
0138       posn_y = rad*cos(phi(i));
0139       dist = (nodes(:,1)-posn_x).^2 + (nodes(:,2)-posn_y).^2;
0140       [jnk,e_node] = min(dist);
0141 
0142       elec(i).z_contact = 0.001;
0143       elec(i).nodes     = e_node;
0144    end
0145    
0146 function inv_mdl = greit_library_model(str)
0147 fmdl = mk_library_model(str);
0148 fmdl = mdl_normalize(fmdl,1);
0149 nelec = numel(fmdl.electrode);
0150 fmdl.stimulation = ...
0151    mk_stim_patterns(nelec,1,[0,1],[0,1],{'no_meas_current'}, 1);
0152 
0153 opt.noise_figure = 1;
0154 opt.target_size = 0.1;
0155 opt.square_pixels = 1;
0156 if strcmp(str(1:8),'cylinder');
0157    opt.distr = 0; % seems best, but only works for cylinders
0158 else
0159    opt.distr = 4; % best bet for the time being
0160 end
0161 
0162 inv_mdl = mk_GREIT_model(fmdl,0.25,5,opt);