Description of mk_common

0001 function inv_mdl= mk_common_model( str, n_elec, varargin )
0002 % MK_COMMON_MODEL: make common EIT models
0003 %
0004 % Utility function to create common EIT FEM models,
0005 % so that users do not need to re-write common code
0006 %
0007 % Usage:
0008 %      inv_mdl = mk_common_model( mdl_string, [n_elec/plane, n_planes])
0009 %
0010 % 2D Models using distmesh (D = show distmesh graphics, d= no graphics)
0011 %   mk_common_model('a2d0c',16)  - 2D circ model using distmesh
0012 %   mk_common_model('b2d1c',16)  - 2D circ model using distmesh ~ 1300 elems
0013 %   mk_common_model('d2d4c',16)  - 2D circ model using distmesh ~ 3200 elems
0014 %      a-j => mesh density
0015 %      2d  => 2d Distmesh model
0016 %      0-4 => element refinement
0017 %      c   => circular mesh
0018 %   mk_common_model('b2d1t2',16)  - 2D circ model using distmesh ~ 1300 elems
0019 %       deformed to T2 thorax shape
0020 %
0021 % 2D Models using distmesh using fixed point electrodes (faster but worse refinement)
0022 %   mk_common_model('a2d0d',16)  - 2D circ model using distmesh
0023 %
0024 % 2D Models circular models:
0025 %   mk_common_model('a2C',16)   - 2D circ model (64 elems) with 16 elecs
0026 %   mk_common_model('b2C',16)   - 2D circ model (256 elems)
0027 %   mk_common_model('c2C',16)   - 2D circ model (576 elems)
0028 %   mk_common_model('d2C',16)   - 2D circ model (1024 elems)
0029 %   mk_common_model('e2C',16)   - 2D circ model (1600 elems)
0030 %   mk_common_model('f2C',16)   - 2D circ model (2304 elems)
0031 %   mk_common_model('g2C',16)   - 2D circ model (3136 elems)
0032 %   mk_common_model('h2C',16)   - 2D circ model (4096 elems)
0033 %   mk_common_model('i2C',16)   - 2D circ model (5184 elems)
0034 %   mk_common_model('j2C',16)   - 2D circ model (6400 elems)
0035 %
0036 %   models with 'c' are point electrode models,
0037 %   models with 'C' use the complete electrode model (with 2 nodes/elec)
0038 %
0039 %   models ??c or ??c0 are rotated by zero.
0040 %   models ??c1, ??c2, ??c3 are rotated by 22.5, 45, 67.5 degrees
0041 %
0042 % 2D Thorax models (levels 1 - 5 from shoulders to abdomen)
0043 %   mk_common_model('b2t2',16)  - 2D Thorax#2 (chest) (256 elems)
0044 %   mk_common_model('c2t4',16)  - 2D Thorax#3 (upper abdomen) (576 elems)
0045 %   - all t1-t5 are available for each a-f models
0046 %
0047 % 2D square models:
0048 %   mk_common_model('a2s',8)   - 2D square model (4x4x2 elems) (max 8 elecs)
0049 %   mk_common_model('b2s',16)  - 2D square model (8x8x2 elems) (16 elecs)
0050 %   mk_common_model('c2s',16)  - 2D square model (16x16x2 elems)
0051 %   mk_common_model('d2s',16)  - 2D square model (24x24x2 elems)
0052 %   mk_common_model('e2s',16)  - 2D square model (32x32x2 elems)
0053 %   mk_common_model('f2s',16)  - 2D square model (40x40x2 elems)
0054 %
0055 %   models ??c or ??c0 are rotated by zero.
0056 %   models ??c1, ??c2, ??c3 are rotated by 22.5, 45, 67.5 degrees
0057 %
0058 % 3D Models:
0059 %   mk_common_model('n3r2',16)  - NP's 3D model with 2 ring electrodes
0060 %
0061 %   mk_common_model('b3cr',[16,3])  - cylinder with 3 rings of 16 elecs
0062 %   mk_common_model('b3t2r',[16,1]) - t2 thorax shape with 1 ring of 16 elecs
0063 %   mk_common_model('b3cz2',[16,1]) - cylinder with 2 rows of 8
0064 %           zigzag pattern elecs. Stimulation treats this as 16x1 pattern
0065 %   mk_common_model('b3cp2',16)      - cylinder with 2 rows of 8
0066 %           elecs in 'planar' pattern. Stim treats this as 16x1 pattern
0067 %
0068 %   mk_common_model('a3cr',16)      - 64 elems * 4 planes
0069 %   mk_common_model('b3cr',16)      - 256 elems * 10 planes
0070 %   mk_common_model('c3cr',16)      - 576 elems * 20 planes
0071 %   mk_common_model('d3cr',16)      - 1024 elems * 40 planes
0072 %
0073 
0074 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0075 % $Id: mk_common_model.html 2819 2011-09-07 16:43:11Z aadler $
0076 
0077 if isstr(str) && strcmp(str,'UNIT_TEST'); do_unit_test; return; end
0078 
0079 
0080 options = {'no_meas_current','no_rotate_meas'};
0081 % n_elec is number of [elec/ring n_rings]
0082 if nargin<2
0083     n_elec= [16,1]; % default
0084 end
0085 if length(n_elec)==1
0086     n_elec= [n_elec,1]; % default
0087 end
0088     
0089 if length(str)<3
0090    error('format specified not recognized')
0091 end
0092 
0093 if strcmp(str(2:3),'2c') || strcmp(str(2:3),'2C')
0094 % 2D circular models
0095    if     str(1)=='a'; layers=  4;
0096    elseif str(1)=='b'; layers=  8;
0097    elseif str(1)=='c'; layers= 12;
0098    elseif str(1)=='d'; layers= 16;
0099    elseif str(1)=='e'; layers= 20;
0100    elseif str(1)=='f'; layers= 24;
0101    elseif str(1)=='g'; layers= 28;
0102    elseif str(1)=='h'; layers= 32;
0103    elseif str(1)=='i'; layers= 36;
0104    elseif str(1)=='j'; layers= 40;
0105    else;  error(['don`t know what to do with option=%s',str]);
0106    end
0107 
0108    inv_mdl = mk_2c_model( n_elec, layers, options );   
0109 
0110    if length(str)==3; str= [str,'0'];end
0111       
0112    inv_mdl = rotate_model( inv_mdl, str2num(str(4)));
0113 
0114    if str(3)=='C' % complete electrode model
0115       inv_mdl = mk_complete_elec_mdl( inv_mdl, layers);
0116    end
0117 
0118 elseif lower(str(2:3))=='2d'
0119    global distmesh_do_graphics;
0120    if str(3)=='d'; distmesh_do_graphics= 0;
0121    else          ; distmesh_do_graphics= 1;
0122    end
0123 
0124    switch str(5)
0125       case 'd' % Deprecated circle functions
0126          inv_mdl= distmesh_2d_model_depr(str, n_elec, options);
0127       case 'c'
0128          inv_mdl= distmesh_2d_model(str, n_elec, options);
0129       case 't'
0130          inv_mdl= distmesh_2d_model(str, n_elec, options);
0131          inv_mdl.fwd_model = thorax_geometry( inv_mdl.fwd_model, str2num(str(6)));
0132       otherwise;
0133          error(['can''t parse command string:', str]);
0134    end
0135 elseif str(2:3)=='2s'
0136 % 2D square models
0137    if     str(1)=='a'; layers=  4;
0138    elseif str(1)=='b'; layers=  8;
0139    elseif str(1)=='c'; layers= 16;
0140    elseif str(1)=='d'; layers= 24;
0141    elseif str(1)=='e'; layers= 32;
0142    elseif str(1)=='f'; layers= 40;
0143    else;  error('don`t know what to do with option=%s',str);
0144    end
0145 
0146    if rem( layers, n_elec(1)/2)~=0; 
0147       error('the %s model can`t support %d electrodes',str,n_elec(1));
0148    end
0149    inv_mdl = mk_2r_model( n_elec, layers, options);
0150 
0151 elseif ( strcmp(str(2:3),'2t') || strcmp(str(2:3),'2T')) && length(str)==4
0152    if     str(1)=='a'; layers=  4;
0153    elseif str(1)=='b'; layers=  8;
0154    elseif str(1)=='c'; layers= 12;
0155    elseif str(1)=='d'; layers= 16;
0156    elseif str(1)=='e'; layers= 20;
0157    elseif str(1)=='f'; layers= 24;
0158    elseif str(1)=='g'; layers= 28;
0159    elseif str(1)=='h'; layers= 32;
0160    elseif str(1)=='i'; layers= 36;
0161    elseif str(1)=='j'; layers= 40;
0162    else;  error(['don`t know what to do with option(1)=',str]);
0163    end
0164 
0165    inv_mdl = mk_2c_model( n_elec, layers, options );   
0166    inv_mdl = rotate_model( inv_mdl, 2); % 45 degrees
0167 
0168    if length(str)==0; str= [str,' '];end
0169 
0170    if str(3)=='T' % complete electrode model
0171       inv_mdl = mk_complete_elec_mdl( inv_mdl, layers);
0172    end
0173       
0174    inv_mdl.fwd_model = thorax_geometry( inv_mdl.fwd_model, str2num(str(4)));
0175 
0176 elseif strcmp( str, 'n3r2')
0177     inv_mdl = mk_n3r2_model( n_elec, options );
0178 elseif strcmp( str, 'n3z') || strcmp(str, 'n3z2')
0179     inv_mdl = mk_n3z_model( n_elec, options );
0180 elseif strcmp( str(2:3), '3c') || strcmp(str(2:3),'3t')
0181    if     str(1)=='a'; xy_layers=  4; z_layers= linspace(-.5,.5,5);
0182    elseif str(1)=='b'; xy_layers=  8; z_layers= linspace(-.7,.7,11);
0183    elseif str(1)=='c'; xy_layers= 12; z_layers= linspace(-.9,.9,21);
0184    elseif str(1)=='d'; xy_layers= 16; z_layers= linspace(-1,1,41);
0185    else;  error(['don`t know what to do with option(1)=',str]);
0186    end
0187 
0188    if str(3)== 'c'
0189       elec_cfg_str= str(4:end);
0190    elseif str(3) == 't'
0191       elec_cfg_str= str(5:end);
0192    else
0193       error(['don`t know what to do with option(3)=',str]);
0194    end
0195 
0196    elec_per_plane = n_elec(1);
0197    spacing=.5;
0198    if     elec_cfg_str=='r';
0199       elec_conf= 'planes';
0200       ne_1  = (n_elec(2)-1)/2;
0201       elec_space= [ne_1:-1:-ne_1]*spacing;
0202    elseif elec_cfg_str=='z2';
0203       elec_conf= 'zigzag';
0204       elec_space= [1,-1]*spacing/2;
0205    elseif elec_cfg_str=='p2';
0206       elec_conf= 'planes';
0207       elec_space= [1,-1]*spacing/2;
0208       elec_per_plane = n_elec(1)/2;
0209    else;
0210       error(['don`t know what to do with option(4)=',str]);
0211    end
0212 
0213    inv_mdl = mk_3c_model( n_elec, xy_layers, z_layers, elec_space, ...
0214                               elec_per_plane, elec_conf, options );
0215 
0216    if str(3) == 't' % thorax models
0217       inv_mdl = rotate_model( inv_mdl, 2); % 45 degrees
0218       inv_mdl.fwd_model = thorax_geometry( inv_mdl.fwd_model, str2num(str(4)));
0219    end
0220 else
0221     error(['Don`t know what to do with option=',str]);
0222 end
0223 
0224 inv_mdl.name= ['EIDORS common_model_',str]; 
0225 inv_mdl= eidors_obj('inv_model', inv_mdl);
0226     
0227 function inv_mdl = distmesh_2d_model(str, n_elec, options);
0228 % This function is an interface to distmesh_2d_model for some common values
0229 % fmdl = dm_2d_circ_pt_elecs( elec_pts, pfix, spacing);
0230 %  params.base_spacing = spacing(1);
0231 %  params.refine_ratio = spacing(2);
0232 %  params.gradient     = spacing(3);
0233    Elec_width= 4; % 2 degrees - electrode width
0234    switch [str(1),str(4)]
0235       case 'j0'; params = [ 55,0,100]./[1000,1,100];
0236       case 'i0'; params = [ 67,0,100]./[1000,1,100];
0237       case 'h0'; params = [ 77,0,100]./[1000,1,100];
0238       case 'g0'; params = [ 87,0,100]./[1000,1,100];
0239       case 'f0'; params = [100,0,100]./[1000,1,100];
0240       case 'e0'; params = [120,0,100]./[1000,1,100];
0241       case 'd0'; params = [150,0,100]./[1000,1,100];
0242       case 'c0'; params = [200,0,100]./[1000,1,100];
0243       case 'b0'; params = [270,0,100]./[1000,1,100];
0244       case 'a0'; params = [500,0,100]./[1000,1,100];
0245 
0246       case 'j1'; params = [ 21,3,5]./[1000,1,100];
0247       case 'i1'; params = [ 23,3,5]./[1000,1,100];
0248       case 'h1'; params = [ 26,3,5]./[1000,1,100];
0249       case 'g1'; params = [ 30,3,5]./[1000,1,100];
0250       case 'f1'; params = [ 35,3,5]./[1000,1,100];
0251       case 'e1'; params = [ 39,3,5]./[1000,1,100];
0252       case 'd1'; params = [ 54,3,5]./[1000,1,100];
0253       case 'c1'; params = [100,3,5]./[1000,1,100];
0254       case 'b1'; params = [180,3,5]./[1000,1,100];
0255       case 'a1'; params = [400,3,5]./[1000,1,100];
0256 
0257       case 'j2'; params = [ 12,5,3]./[1000,1,100];
0258       case 'i2'; params = [ 14,5,3]./[1000,1,100];
0259       case 'h2'; params = [ 16,5,3]./[1000,1,100];
0260       case 'g2'; params = [ 19,5,3]./[1000,1,100];
0261       case 'f2'; params = [ 21,5,3]./[1000,1,100];
0262       case 'e2'; params = [ 39,5,3]./[1000,1,100];
0263       case 'd2'; params = [ 50,5,3]./[1000,1,100];
0264       case 'c2'; params = [100,5,3]./[1000,1,100];
0265       case 'b2'; params = [200,5,3]./[1000,1,100];
0266       case 'a2'; params = [500,5,3]./[1000,1,100];
0267 
0268       case 'j3'; params = [  6,10,2]./[1000,1,100];
0269       case 'i3'; params = [  7,10,2]./[1000,1,100];
0270       case 'h3'; params = [  8,10,2]./[1000,1,100];
0271       case 'g3'; params = [  9,10,2]./[1000,1,100];
0272       case 'f3'; params = [ 10,10,2]./[1000,1,100];
0273       case 'e3'; params = [ 13,10,2]./[1000,1,100];
0274       case 'd3'; params = [ 20,10,2]./[1000,1,100];
0275       case 'c3'; params = [ 70,10,2]./[1000,1,100];
0276       case 'b3'; params = [150,10,2]./[1000,1,100];
0277       case 'a3'; params = [250,10,2]./[1000,1,100];
0278 
0279 % We set refinement 4==3. This is OK for this mdl density
0280       case 'j4'; params = [  6,10,2]./[1000,1,100];
0281       case 'i4'; params = [  7,10,2]./[1000,1,100];
0282       case 'h4'; params = [  8,10,2]./[1000,1,100];
0283       case 'g4'; params = [  9,10,2]./[1000,1,100];
0284       case 'f4'; params = [ 10,10,2]./[1000,1,100];
0285       case 'e4'; params = [ 13,10,2]./[1000,1,100];
0286       case 'd4'; params = [ 20,10,2]./[1000,1,100];
0287       case 'c4'; params = [ 70,10,2]./[1000,1,100];
0288       case 'b4'; params = [150,10,2]./[1000,1,100];
0289       case 'a4'; params = [250,10,2]./[1000,1,100];
0290 
0291       otherwise; error('don`t know what to do with option=%s',str);
0292    end
0293    ea = Elec_width/2 *(2*pi/360);
0294    for i=1:n_elec(1); 
0295      ai = (i-1)/n_elec(1) * 2*pi;
0296      elec_pts{i} = [sin(ai+ea),cos(ai+ea);sin(ai-ea),cos(ai-ea)];
0297    end
0298    fwd_mdl= dm_2d_circ_pt_elecs( elec_pts, [], params);
0299    inv_mdl= add_params_2d_mdl( fwd_mdl, n_elec(1), options);
0300 
0301 % THIS FUNCTION IS DEPRECATED (from EIDORS 3.3)
0302 function inv2d = distmesh_2d_model_depr(str, n_elec, options);
0303    switch str(1)
0304       case 'a'; n_nodes=  50;
0305       case 'b'; n_nodes= 100;
0306       case 'c'; n_nodes= 200;
0307       case 'd'; n_nodes= 400;
0308       case 'e'; n_nodes= 800;
0309       case 'f'; n_nodes=1200;
0310       case 'g'; n_nodes=1800;
0311       case 'h'; n_nodes=2400;
0312       case 'i'; n_nodes=3000;
0313       case 'j'; n_nodes=4000;
0314       otherwise; error('don`t know what to do with option=%s',str);
0315    end
0316  
0317    refine_level= abs(str(4))-'0';
0318 
0319    elec_width= .1;
0320    th=linspace(0,2*pi,n_elec(1)+1)';th(end)=[];
0321    elec_posn= [sin(th),cos(th)];
0322    [elec_nodes, refine_nodes] = dm_mk_elec_nodes( elec_posn, ...
0323           elec_width, refine_level);
0324    fd=inline('sqrt(sum(p.^2,2))-1','p');
0325    bbox = [-1,-1;1,1];
0326    z_contact = 0.01;
0327    fwd_mdl= dm_mk_fwd_model( fd, [], n_nodes, bbox, ...
0328                              elec_nodes, refine_nodes, z_contact);
0329 
0330    inv2d= add_params_2d_mdl( fwd_mdl, n_elec(1), options);
0331 
0332 
0333 function inv2d= mk_2c_model( n_elec, n_circles, options )
0334 
0335     n_elec= n_elec(1);
0336     params= mk_circ_tank(n_circles, [], n_elec); 
0337     inv2d= add_params_2d_mdl( params, n_elec, options);
0338    
0339 
0340 function inv2d= mk_2r_model( n_elec, xy_size, options)
0341     if length(xy_size)==1; xy_size= xy_size*[1,1]; end
0342     xy_size= xy_size+1;
0343 
0344 % TODO: To keep elements square, we should scale the 1's
0345     [x,y]= meshgrid( linspace(-1,1,xy_size(1)), ...
0346                      linspace(-1,1,xy_size(2)));
0347     x=x';y=y';
0348     fmdl.nodes= [x(:),y(:)];
0349     k= 1:xy_size(1)-1;
0350     elem_frac = [ k;k+1;k+xy_size(2); ...
0351                   k+1;k+xy_size(2);k+xy_size(2)+1];
0352     elem_frac= reshape(elem_frac, 3,[])';
0353     fmdl.elems=  [];
0354     for j=0:xy_size(2)-2
0355        fmdl.elems=  [fmdl.elems; elem_frac + xy_size(2)*j];
0356     end
0357 
0358     fmdl.boundary = find_boundary( fmdl.elems);
0359 
0360     % put 1/4 of elecs on each side
0361     tb_elecs= linspace(1, xy_size(1), 1+2*n_elec(1)/4); 
0362     tb_elecs= tb_elecs(2:2:end);
0363     sd_elecs= linspace(1, xy_size(2), 1+2*n_elec(1)/4);
0364     sd_elecs= sd_elecs(2:2:end);
0365     
0366     el_nodes= [];
0367     % Top nodes -left to right
0368     bdy_nodes= (1:xy_size(1)) + xy_size(1)*(xy_size(2)-1); 
0369     el_nodes= [el_nodes, bdy_nodes(tb_elecs)];
0370     % Right nodes - top to bottom
0371     bdy_nodes= (1:xy_size(2))*xy_size(1); 
0372     el_nodes= [el_nodes, bdy_nodes(fliplr(sd_elecs))];
0373     % Bottom nodes - right to left
0374     bdy_nodes= 1:xy_size(1); 
0375     el_nodes= [el_nodes, bdy_nodes(fliplr(tb_elecs))];
0376     % Left nodes - bottom to top
0377     bdy_nodes= (0:xy_size(2)-1)*xy_size(1)+1; 
0378     el_nodes= [el_nodes, bdy_nodes(sd_elecs)];
0379 
0380 %   trimesh(fmdl.elems,fmdl.nodes(:,1), fmdl.nodes(:,2));
0381     for i=1:n_elec(1)
0382        n= el_nodes(i);
0383        fmdl.electrode(i).nodes= n;
0384        fmdl.electrode(i).z_contact= .001; % choose a low value
0385 %      plot(fmdl.nodes(n,1),fmdl.nodes(n,2),'*'); pause;
0386     end
0387     inv2d= add_params_2d_mdl( fmdl, n_elec(1), options);
0388 
0389 % params is the part of the fwd_model
0390 function inv2d= add_params_2d_mdl( params, n_elec, options);
0391     n_rings= 1;
0392     [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0393     params.stimulation= st;
0394     params.meas_select= els;
0395     params.solve=      'aa_fwd_solve';
0396     params.system_mat= 'aa_calc_system_mat';
0397     params.jacobian=   'aa_calc_jacobian';
0398     params.normalize_measurements= 0;
0399     params.np_fwd_solve.perm_sym= '{n}';
0400     mdl_2d   = eidors_obj('fwd_model', params);
0401 
0402     inv2d.solve=       'aa_inv_solve';
0403     inv2d.hyperparameter.value = 3e-2;
0404     %inv2d.hyperparameter.func = 'choose_noise_figure';
0405     %inv2d.hyperparameter.noise_figure= 1;
0406     %inv2d.hyperparameter.tgt_elems= 1:4;
0407      inv2d.RtR_prior= 'laplace_image_prior';
0408     %inv2d.RtR_prior= 'gaussian_HPF_prior';
0409     inv2d.jacobian_bkgnd.value= 1;
0410     inv2d.reconst_type= 'difference';
0411     inv2d.fwd_model= mdl_2d;
0412 
0413 function inv3d = mk_3c_model( n_elec, xy_layers, z_layers, elec_space, ...
0414                               elec_per_plane, elec_conf, options );
0415 
0416     e_layers=[];
0417     for es= elec_space;
0418        ff= abs(z_layers  -es);
0419        ff= find(ff==min(ff));
0420        e_layers= [e_layers,ff(1)];
0421     end
0422 
0423     params= mk_circ_tank( xy_layers, z_layers, ...
0424            { elec_conf, elec_per_plane, e_layers} );
0425 
0426     [st, els]= mk_stim_patterns(n_elec(1), n_elec(2), '{ad}','{ad}', options, 10);
0427 
0428     params.stimulation= st;
0429     params.meas_select= els;
0430     params.solve=      'aa_fwd_solve';
0431     params.system_mat= 'aa_calc_system_mat';
0432     params.jacobian=   'aa_calc_jacobian';
0433     params.normalize_measurements= 0;
0434     params.np_fwd_solve.perm_sym= '{n}';
0435     fm3d = eidors_obj('fwd_model', params);
0436 
0437     inv3d.name=  'EIT inverse: 3D';
0438     inv3d.solve= @time_prior_solve;
0439     inv3d.time_prior_solve.time_steps=   0;
0440     inv3d.time_smooth_prior.space_prior = @noser_image_prior;
0441     inv3d.time_smooth_prior.time_weight = 0;
0442     inv3d.time_prior_solve.time_steps   = 0;
0443 
0444     inv3d.hyperparameter.value = 3e-2;
0445     inv3d.reconst_type= 'difference';
0446     inv3d.jacobian_bkgnd.value= 1;
0447     inv3d.fwd_model= fm3d;
0448 
0449 function inv_mdl = mk_n3z_model( n_elec, options );
0450    inv_mdl= mk_n3r2_model( n_elec, options);
0451    fwd_mdl= inv_mdl.fwd_model;
0452    renumber= [1:2:15; 18:2:32];
0453    fwd_mdl.electrode= fwd_mdl.electrode(renumber(:));
0454    n_rings= 1;
0455    [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0456    fwd_mdl.stimulation= st;
0457    fwd_model.meas_select= els;
0458    inv_mdl.fwd_model= fwd_mdl;
0459    inv_mdl.name= 'NP 3D model with zigzag electrodes';
0460 
0461 function inv_mdl = mk_n3r2_model( n_elec, options );
0462    if ~exist('OCTAVE_VERSION');
0463       load( 'datareal.mat' );
0464    else
0465       load( file_in_loadpath( 'datareal.mat' ));
0466    end
0467    fmdl.nodes= vtx;
0468    fmdl.elems= simp;
0469    fmdl.boundary= find_boundary( simp );
0470 
0471    fmdl.solve=      'np_fwd_solve';
0472    fmdl.jacobian=   'np_calc_jacobian';
0473    fmdl.system_mat= 'np_calc_system_mat';
0474 
0475    for i=1:length(zc)
0476        electrodes(i).z_contact= zc(i);
0477        electrodes(i).nodes=     unique( elec(i,:) );
0478    end
0479 
0480    fmdl.gnd_node=           gnd_ind;
0481    fmdl.electrode =         electrodes;
0482    fmdl.np_fwd_solve.perm_sym =     '{n}';
0483 
0484    [I,Ib] = set_3d_currents(protocol, elec, ...
0485                fmdl.nodes, fmdl.gnd_node, no_pl);
0486 
0487    % get the measurement patterns, only indH is used in this model
0488    %   here we only want to get the meas pattern from 'get_3d_meas',
0489    %   not the voltages, so we enter zeros
0490    [jnk,jnk,indH,indV,jnk] = get_3d_meas( elec, ...
0491             fmdl.nodes, zeros(size(I)), Ib, no_pl );
0492    n_elec= size(elec,1);
0493    n_meas= size(indH,1) / size(Ib,2);
0494    for i=1:size(Ib,2)
0495       fmdl.stimulation(i).stimulation= 'mA';
0496       fmdl.stimulation(i).stim_pattern= Ib(:,i);
0497 
0498       idx= ( 1+ (i-1)*n_meas ):( i*n_meas );
0499       fmdl.stimulation(i).meas_pattern= sparse ( ...
0500               (1:n_meas)'*[1,1], ...
0501               indH( idx, : ), ...
0502               ones(n_meas,2)*[1,0;0,-1], ...
0503               n_meas, n_elec );
0504    end
0505    fmdl= eidors_obj('fwd_model', fmdl);
0506 
0507    inv_mdl.name=         'Nick Polydorides EIT inverse';
0508    inv_mdl.solve=       'np_inv_solve';
0509    inv_mdl.hyperparameter.value = 1e-2;
0510    inv_mdl.RtR_prior= 'np_calc_image_prior';
0511    inv_mdl.np_calc_image_prior.parameters= [3 1]; % see iso_f_smooth: deg=1, w=1
0512    inv_mdl.reconst_type= 'difference';
0513    inv_mdl.jacobian_bkgnd.value= 1;
0514    inv_mdl.fwd_model= fmdl;
0515 
0516 
0517 function inv3d= mk_b3r1_model( n_elec, options )
0518     n_rings= 1;
0519     levels= [-.5:.1:.5]; 
0520     e_levels= 6; 
0521     nr= 8;
0522 
0523     params= mk_circ_tank( nr, levels, { 'planes', n_elec, e_levels } );
0524     [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0525 
0526     params.stimulation= st;
0527     params.meas_select= els;
0528     params.solve=      'aa_fwd_solve';
0529     params.system_mat= 'aa_calc_system_mat';
0530     params.jacobian=   'aa_calc_jacobian';
0531     params.normalize_measurements= 0;
0532     params.np_fwd_solve.perm_sym= '{n}';
0533     mdl_3d = eidors_obj('fwd_model', params);
0534 
0535     inv3d.name = 'EIT inverse: 3D';
0536     inv3d.solve=       'aa_inv_solve';
0537     %inv3d.solve=       'aa_inv_conj_grad';
0538     inv3d.hyperparameter.value = 1e-5;
0539     inv3d.RtR_prior= 'laplace_image_prior';
0540     %inv3d.RtR_prior= 'gaussian_HPF_prior';
0541     inv3d.jacobian_bkgnd.value= 1;
0542     inv3d.reconst_type= 'difference';
0543     inv3d.fwd_model= mdl_3d;
0544 
0545 function inv3d= mk_b3r2_model( n_elec, nr, options )
0546     n_rings= 2;
0547     z_axis = [0:.1:1];
0548     e_levels= [4,8]; 
0549     nr= 4;
0550     n_elec = 8;
0551    
0552     params= mk_circ_tank( nr, z_axis, { 'planes', n_elec, e_levels } );
0553     [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0554 
0555     params.stimulation= st;
0556     params.meas_select= els;
0557     params.solve=      'aa_fwd_solve';
0558     params.system_mat= 'aa_calc_system_mat';
0559     params.jacobian=   'aa_calc_jacobian';
0560     params.normalize_measurements= 0;
0561     params.np_fwd_solve.perm_sym= '{n}';
0562     mdl_3d = eidors_obj('fwd_model', params);
0563     
0564     % Specify number of levels in mesh for imaging slices
0565     num_levs = length(e_levels);
0566     levels = inf*ones(num_levs,3);
0567     levels(:,3) = e_levels / (length(z_axis)-1);
0568     levels(:,4) = ones(num_levs,1);
0569     levels(:,5) = (1:num_levs)';    
0570     mdl_3d.levels = levels;
0571     
0572     inv3d.name = 'EIT inverse: 3D';
0573     inv3d.solve=       'aa_inv_solve';
0574     %inv3d.solve=       'aa_inv_conj_grad';
0575     inv3d.hyperparameter.value = 1e-5;
0576     inv3d.RtR_prior= 'laplace_image_prior';
0577     %inv3d.RtR_prior= 'gaussian_HPF_prior';
0578     inv3d.jacobian_bkgnd.value= 1;
0579     inv3d.reconst_type= 'difference';
0580     inv3d.fwd_model= mdl_3d;
0581 
0582 % rotate model rotate_model/16 times around
0583 function inv_mdl = rotate_model( inv_mdl, rotate_mdl);
0584     inv_mdl = turn_model( inv_mdl, pi/8*rotate_mdl );
0585 
0586     n_elec= length( inv_mdl.fwd_model.electrode );
0587     renum = rem(  rotate_mdl*n_elec/16 + (0:n_elec-1),n_elec)+1;
0588     renum = floor(renum); % round is not quite right for all model
0589                           % cases, but no errors.
0590     inv_mdl.fwd_model.electrode = ...
0591        inv_mdl.fwd_model.electrode( renum);
0592 
0593 % rotate model rotate_model/16 times around
0594 function inv_mdl = turn_model( inv_mdl, angle );
0595     nodes= inv_mdl.fwd_model.nodes;
0596     cos_rot = cos( angle );
0597     sin_rot = sin( angle );
0598     nodes(:,1)= inv_mdl.fwd_model.nodes(:,1:2)*[ cos_rot;-sin_rot];
0599     nodes(:,2)= inv_mdl.fwd_model.nodes(:,1:2)*[ sin_rot; cos_rot];
0600     inv_mdl.fwd_model.nodes= nodes;
0601 
0602 function inv_mdl = mk_complete_elec_mdl( inv_mdl, layers);
0603       inv_mdl = turn_model( inv_mdl, 2*pi/4/layers/2 );
0604 
0605       bdy= inv_mdl.fwd_model.boundary;
0606       for i=1:length(inv_mdl.fwd_model.electrode);
0607          enode= inv_mdl.fwd_model.electrode(i).nodes;
0608          ff= find( enode== bdy(:,1) );
0609          inv_mdl.fwd_model.electrode(i).nodes = bdy(ff,:);
0610       end
0611 
0612 %%%%%%%%%%%%%%%%%%%%%% TESTS %%%%%%%%%%%%%%%%%%%%%%%
0613 function do_unit_test
0614 
0615 % 2D Circular Models
0616 for j=('a'+0):('j'+0)
0617     mk_common_model(sprintf('%c2C2',j),16);
0618     mk_common_model(sprintf('%c2c0',j),16);
0619     mk_common_model(sprintf('%c2t3',j),16);
0620     mk_common_model(sprintf('%c2T4',j),16);
0621 end;
0622 
0623 for j=('a'+0):('f'+0)
0624     mk_common_model(sprintf('%c2s',j),8);
0625 end;
0626 
0627 % 3D Models:
0628     mk_common_model('n3r2',16);
0629  %  mk_common_model('n3z',16);
0630  
0631     mk_common_model('b3cr',[16,3]);
0632     mk_common_model('b3t2r',[16,1]);
0633     mk_common_model('b3cz2',[16,1]);
0634     mk_common_model('b3cp2',16);
0635  
0636     mk_common_model('a3cr',16);
0637     mk_common_model('b3cr',16);
0638     mk_common_model('c3cr',16);
0639     mk_common_model('d3cr',16);
0640 
0641 % Distmesh models
0642 for i=0:4; for j=('a'+0):('j'+0)
0643     mk_common_model(sprintf('%c2d%dd',j,i),16);
0644 end; end
0645
mk_common_model

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
