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 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 V1 (GREIT matrix for reconstruction to circle, 2009)
0018 %   mk_common_gridmdl('GREITc1') - 32x32 with Circle shape
0019 %
0020 % Sheffield Backprojection
0021 %   mk_common_gridmdl('b2d','backproj') - 32x32 with Diamond shape
0022 %   mk_common_gridmdl('b2c','backproj') - 32x32 with Circular shape
0023 %   mk_common_gridmdl('backproj') - 32x32 with Circular shape
0024 %
0025 % COPYRIGHT NOTICE FOR BACKPROJECTION MATRIX:
0026 %   This matrix is copyright DC Barber and BH Brown at
0027 %   University of Sheffield. It may be used free of
0028 %   charge for research and non-commercial purposes.
0029 %   Commercial applications require a licence from the
0030 %   University of Sheffield.
0031 %
0032 
0033 % (C) 2008 Andy Adler. License: GPL version 2 or version 3
0034 % $Id: mk_common_gridmdl.html 2819 2011-09-07 16:43:11Z aadler $
0035 
0036 name = str;
0037 if strcmp(str,'backproj')
0038    str= 'b2c';
0039    RM= 'backproj';
0040 elseif strcmp(str,'GREITc1');
0041    str= 'b2c';
0042    RM= 'GREITc1';
0043 else
0044    name = str;
0045 end
0046 
0047 n_elec= 16;
0048 
0049       space = linspace( -1, 1, 32+1 );
0050       fmdl= mk_grid_model( [], space, space);
0051       space_avg = conv2(space, [1,1]/2,'valid'); % average of each box
0052       [x,y] = ndgrid( space_avg, space_avg);
0053 switch str
0054    case 'b2c'
0055       inside=  (x(:).^2 + y(:).^2)<1.10 ;
0056 
0057    case 'b2d'
0058       inside=  (abs(x(:)) + abs(y(:)))<1.55 ;
0059 
0060    case 'b2t1'; inside = inside_thorax(x,y,1);
0061    case 'b2t2'; inside = inside_thorax(x,y,2);
0062    case 'b2t3'; inside = inside_thorax(x,y,3);
0063    case 'b2t4'; inside = inside_thorax(x,y,4);
0064    case 'b2t5'; inside = inside_thorax(x,y,5);
0065 
0066 
0067    otherwise
0068       error(['mdl_string ',str,' not understood']);
0069 end
0070 
0071 if isstr(RM)
0072    switch RM
0073      case 'backproj'; RM= get_Sheffield_Backproj;
0074      case 'GREITc1';  RM= get_GREIT_c1;
0075      otherwise;       error('RM string not understood');
0076    end 
0077 end
0078       ff = find(~inside);
0079       fmdl.elems([2*ff, 2*ff-1],:)= [];
0080       fmdl.coarse2fine([2*ff, 2*ff-1],:)= [];
0081       fmdl.coarse2fine(:,ff)= [];
0082 
0083 fmdl.electrode = mk_electrode_locns( fmdl.nodes, n_elec );
0084 
0085 inv_mdl = eidors_obj('inv_model',['mk_common_gridmdl: ',name]);
0086 inv_mdl.reconst_type= 'difference';
0087 inv_mdl.fwd_model= fmdl;
0088 inv_mdl.solve_use_matrix.RM = resize_if_reqd(RM,inside);
0089 % solve_use_matrix has a c2f mapping field map, which
0090 % it may be useful to populate in this case
0091 %  inv_mdl.solve_use_matrix.map = map;
0092 inv_mdl.solve = @solve_use_matrix;
0093 inv_mdl.fwd_model.normalize_measurements= 1;
0094 [st, els]= mk_stim_patterns(16, 1, '{ad}','{ad}', {}, 10);
0095 inv_mdl.fwd_model.meas_select= els;
0096 
0097 function inside = inside_thorax(x,y,level);
0098    [x_bdy, y_bdy ] = thorax_geometry(level,1); % normalized
0099    inside = inpolygon(x(:), y(:), x_bdy, y_bdy); 
0100 
0101 function RM = resize_if_reqd(RM,inside);
0102    szRM = size(RM,1);
0103    if sum(inside) == szRM
0104       % RM is fine
0105    elseif size(inside,1) == szRM
0106       RM = RM(inside,:);
0107    else
0108       error('mismatch in size of provided RecMatrix');
0109    end
0110 
0111 function RM = get_Sheffield_Backproj
0112    load Sheffield_Backproj_Matrix.mat
0113    RM = unpack_matrix(Sheffield_Backproj_Matrix);
0114 
0115 function RM= get_GREIT_c1;
0116    load GREIT_v10_Circ_Matrix.mat
0117    RM = unpack_matrix(GREIT_v10_Circ_Matrix);
0118 
0119 
0120 function RM= unpack_matrix(MM)
0121    RM = unpack_reconst_matrix(MM, 16, 32);
0122 
0123 
0124 function elec = mk_electrode_locns( nodes, n_elec );
0125    phi = linspace(0, 2*pi, n_elec+1); 
0126    phi(end) = [];
0127    rad = 1;
0128 
0129    for i= 1:n_elec
0130       posn_x = rad*sin(phi(i));
0131       posn_y = rad*cos(phi(i));
0132       dist = (nodes(:,1)-posn_x).^2 + (nodes(:,2)-posn_y).^2;
0133       [jnk,e_node] = min(dist);
0134 
0135       elec(i).z_contact = 0.001;
0136       elec(i).nodes     = e_node;
0137    end
0138     
0139