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 6926 2024-05-31 15:34:13Z bgrychtol $
0041 
0042 if strcmp(str,'UNIT_TEST'), do_unit_test; return, end
0043 
0044 name = str;
0045 if strcmp(str,'backproj')
0046    str= 'b2c';
0047    RM= 'backproj';
0048 elseif strcmp(str,'GREITc1');
0049    str= 'b2c';
0050    RM= 'GREITc1';
0051 else
0052    name = str;
0053 end
0054 
0055 n_elec= 16;
0056 
0057       space = linspace( -1, 1, 32+1 );
0058       fmdl= mk_grid_model( [], space, space);
0059       space_avg = conv2(space, [1,1]/2,'valid'); % average of each box
0060       [x,y] = ndgrid( space_avg, space_avg);
0061 switch str
0062    case 'b2c'
0063       inside=  (x(:).^2 + y(:).^2)<1.10 ;
0064 
0065    case 'b2d'
0066       inside=  (abs(x(:)) + abs(y(:)))<1.55 ;
0067 
0068    case 'b2t1'; inside = inside_thorax(x,y,1);
0069    case 'b2t2'; inside = inside_thorax(x,y,2);
0070    case 'b2t3'; inside = inside_thorax(x,y,3);
0071    case 'b2t4'; inside = inside_thorax(x,y,4);
0072    case 'b2t5'; inside = inside_thorax(x,y,5);
0073    
0074    case 'GREIT'; inv_mdl = greit_library_model(RM); return
0075 
0076    otherwise
0077       error(['mdl_string ',str,' not understood']);
0078 end
0079 
0080 if ischar(RM)
0081    switch RM
0082      case 'backproj'; RM= get_Sheffield_Backproj;
0083      case 'GREITc1';  RM= get_GREIT_c1;
0084      otherwise;       error('RM string not understood');
0085    end 
0086 end
0087       ff = find(~inside);
0088       fmdl.elems([2*ff, 2*ff-1],:)= [];
0089       fmdl.coarse2fine([2*ff, 2*ff-1],:)= [];
0090       fmdl.coarse2fine(:,ff)= [];
0091 
0092 fmdl.electrode = mk_electrode_locns( fmdl.nodes, n_elec );
0093 
0094 inv_mdl = eidors_obj('inv_model',['mk_common_gridmdl: ',name]);
0095 inv_mdl.reconst_type= 'difference';
0096 inv_mdl.fwd_model= fmdl;
0097 inv_mdl.solve_use_matrix.RM = resize_if_reqd(RM,inside);
0098 % solve_use_matrix has a c2f mapping field map, which
0099 % it may be useful to populate in this case
0100 %  inv_mdl.solve_use_matrix.map = map;
0101 inv_mdl.solve = @solve_use_matrix;
0102 inv_mdl.fwd_model = mdl_normalize(inv_mdl.fwd_model, 1);
0103 [st, els]= mk_stim_patterns(16, 1, '{ad}','{ad}', {}, 10);
0104 inv_mdl.fwd_model.meas_select= els;
0105 
0106 % standard field order
0107 inv_mdl = eidors_obj('set', inv_mdl);
0108 
0109 function inside = inside_thorax(x,y,level);
0110    [x_bdy, y_bdy ] = thorax_geometry(level,1); % normalized
0111    inside = inpolygon(x(:), y(:), x_bdy, y_bdy); 
0112 
0113 function RM = resize_if_reqd(RM,inside);
0114    szRM = size(RM,1);
0115    if sum(inside) == szRM
0116       % RM is fine
0117    elseif size(inside,1) == szRM
0118       RM = RM(inside,:);
0119    else
0120       error('mismatch in size of provided RecMatrix');
0121    end
0122 
0123 function RM = get_Sheffield_Backproj
0124    load Sheffield_Backproj_Matrix.mat
0125    RM = unpack_matrix(Sheffield_Backproj_Matrix);
0126 
0127 function RM= get_GREIT_c1;
0128    load GREIT_v10_Circ_Matrix.mat
0129    RM = unpack_matrix(GREIT_v10_Circ_Matrix);
0130 
0131 
0132 function RM= unpack_matrix(MM)
0133    RM = unpack_reconst_matrix(MM, 16, 32);
0134 
0135 
0136 function elec = mk_electrode_locns( nodes, n_elec );
0137    phi = linspace(0, 2*pi, n_elec+1); 
0138    phi(end) = [];
0139    rad = 1;
0140 
0141    for i= 1:n_elec
0142       posn_x = rad*sin(phi(i));
0143       posn_y = rad*cos(phi(i));
0144       dist = (nodes(:,1)-posn_x).^2 + (nodes(:,2)-posn_y).^2;
0145       [jnk,e_node] = min(dist);
0146 
0147       elec(i).z_contact = 0.001;
0148       elec(i).nodes     = e_node;
0149    end
0150    
0151 function inv_mdl = greit_library_model(str)
0152 fmdl = mk_library_model(str);
0153 fmdl = mdl_normalize(fmdl,1);
0154 nelec = numel(fmdl.electrode);
0155 fmdl.stimulation = ...
0156    mk_stim_patterns(nelec,1,[0,1],[0,1],{'no_meas_current'}, 1);
0157 
0158 opt.noise_figure = 1;
0159 opt.target_size = 0.1;
0160 opt.square_pixels = 1;
0161 if strcmp(str(1:8),'cylinder');
0162    opt.distr = 0; % seems best, but only works for cylinders
0163 else
0164    opt.distr = 4; % best bet for the time being
0165 end
0166 
0167 inv_mdl = mk_GREIT_model(fmdl,0.25,5,opt);
0168 
0169 function do_unit_test
0170 % GREIT MODELS - calculated for library models
0171 mk_common_gridmdl('GREIT','cylinder_16x1el_coarse');
0172 
0173 % GREIT V1 (GREIT matrix for reconstruction to circle, 2009)
0174 mk_common_gridmdl('GREITc1'); % 32x32 with Circle shape
0175 
0176 % Sheffield Backprojection
0177 mk_common_gridmdl('b2d','backproj'); % 32x32 with Diamond shape
0178 mk_common_gridmdl('b2c','backproj'); % 32x32 with Circular shape
0179 mk_common_gridmdl('backproj'); % 32x32 with Circular shape