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 %   mk_common_model('g2s',16)  - 2D square model (56x56x2 elems)
0055 %   mk_common_model('h2s',16)  - 2D square model (80x80x2 elems)
0056 %   mk_common_model('i2s',16)  - 2D square model (96x96x2 elems)
0057 %   mk_common_model('j2s',16)  - 2D square model (144x144x2 elems)
0058 %   mk_common_model('k2s',16)  - 2D square model (200x200x2 elems)
0059 %
0060 %   models ??c or ??c0 are rotated by zero.
0061 %   models ??c1, ??c2, ??c3 are rotated by 22.5, 45, 67.5 degrees
0062 %
0063 % 3D Models:
0064 %   mk_common_model('n3r2',[16,2])  - NP's 3D model with 2 rings of 16 elecs
0065 %
0066 %   mk_common_model('b3cr',[16,3])  - cylinder with 3 rings of 16 elecs
0067 %   mk_common_model('b3t2r',[16,1]) - t2 thorax shape with 1 ring of 16 elecs
0068 %   mk_common_model('b3cz2',[16,1]) - cylinder with 2 rows of 8
0069 %           zigzag pattern elecs. Stimulation treats this as 16x1 pattern
0070 %   mk_common_model('b3cp2',16)      - cylinder with 2 rows of 8
0071 %           elecs in 'planar' pattern. Stim treats this as 16x1 pattern
0072 %
0073 %   mk_common_model('a3cr',16)      - 64 elems * 4 planes
0074 %   mk_common_model('b3cr',16)      - 256 elems * 10 planes
0075 %   mk_common_model('c3cr',16)      - 576 elems * 20 planes
0076 %   mk_common_model('d3cr',16)      - 1024 elems * 40 planes
0077 %   mk_common_model('e3cr',16)      - 1600 elems * 60 planes
0078 %   mk_common_model('f3cr',16)      - 2304 elems * 80 planes
0079 
0080 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0081 % $Id: mk_common_model.m 5112 2015-06-14 13:00:41Z aadler $
0082 
0083 if ischar(str) && strcmp(str,'UNIT_TEST'); do_unit_test; return; end
0084 
0085 
0086 options = {'no_meas_current','no_rotate_meas'};
0087 % n_elec is number of [elec/ring n_rings]
0088 if nargin<2
0089     n_elec= [16,1]; % default
0090 end
0091 if length(n_elec)==1
0092     n_elec= [n_elec,1]; % default
0093 end
0094     
0095 if length(str)<3
0096    error('format specified not recognized')
0097 end
0098 
0099 if strcmp(str(2:3),'2c') || strcmp(str(2:3),'2C')
0100 % 2D circular models
0101    switch str(1)
0102       case 'a'; layers=  4;
0103       case 'b'; layers=  8;
0104       case 'c'; layers= 12;
0105       case 'd'; layers= 16;
0106       case 'e'; layers= 20;
0107       case 'f'; layers= 24;
0108       case 'g'; layers= 28;
0109       case 'h'; layers= 32;
0110       case 'i'; layers= 36;
0111       case 'j'; layers= 40;
0112       case 'k'; layers= 44;
0113       case 'l'; layers= 48;
0114       case 'm'; layers= 52;
0115       otherwise; error(['don`t know what to do with option=%s',str]);
0116    end
0117 
0118    inv_mdl = mk_2c_model( n_elec, layers, options );   
0119 
0120    if length(str)==3; str= [str,'0'];end
0121       
0122    inv_mdl = rotate_model( inv_mdl, str2num(str(4)));
0123 
0124    if str(3)=='C' % complete electrode model
0125       inv_mdl = mk_complete_elec_mdl( inv_mdl, layers);
0126    end
0127 
0128 elseif lower(str(2:3))=='2d'
0129    global distmesh_do_graphics;
0130    if str(3)=='d'; distmesh_do_graphics= 0;
0131    else          ; distmesh_do_graphics= 1;
0132    end
0133 
0134    switch str(5)
0135       case 'd' % Deprecated circle functions
0136          inv_mdl= distmesh_2d_model_depr(str, n_elec, options);
0137       case 'c'
0138          inv_mdl= distmesh_2d_model(str, n_elec, options);
0139       case 't'
0140          inv_mdl= distmesh_2d_model(str, n_elec, options);
0141          inv_mdl.fwd_model = thorax_geometry( inv_mdl.fwd_model, str2num(str(6)));
0142       otherwise;
0143          error(['can''t parse command string:', str]);
0144    end
0145 elseif str(2:3)=='2s'
0146 % 2D square models
0147    if     str(1)=='a'; layers=  4;
0148    elseif str(1)=='b'; layers=  8;
0149    elseif str(1)=='c'; layers= 16;
0150    elseif str(1)=='d'; layers= 24;
0151    elseif str(1)=='e'; layers= 32;
0152    elseif str(1)=='f'; layers= 40;
0153    elseif str(1)=='g'; layers= 56;
0154    elseif str(1)=='h'; layers= 80;
0155    elseif str(1)=='i'; layers= 96;
0156    elseif str(1)=='j'; layers= 144;
0157    elseif str(1)=='k'; layers= 200;
0158    else;  error('don`t know what to do with option=%s',str);
0159    end
0160 
0161    if rem( layers, n_elec(1)/2)~=0; 
0162       error('the %s model can`t support %d electrodes',str,n_elec(1));
0163    end
0164    inv_mdl = mk_2r_model( n_elec, layers, options);
0165 
0166 elseif ( strcmp(str(2:3),'2t') || strcmp(str(2:3),'2T')) && length(str)==4
0167    if     str(1)=='a'; layers=  4;
0168    elseif str(1)=='b'; layers=  8;
0169    elseif str(1)=='c'; layers= 12;
0170    elseif str(1)=='d'; layers= 16;
0171    elseif str(1)=='e'; layers= 20;
0172    elseif str(1)=='f'; layers= 24;
0173    elseif str(1)=='g'; layers= 28;
0174    elseif str(1)=='h'; layers= 32;
0175    elseif str(1)=='i'; layers= 36;
0176    elseif str(1)=='j'; layers= 40;
0177    else;  error(['don`t know what to do with option(1)=',str]);
0178    end
0179 
0180    inv_mdl = mk_2c_model( n_elec, layers, options );   
0181    inv_mdl = rotate_model( inv_mdl, 2); % 45 degrees
0182 
0183    if length(str)==0; str= [str,' '];end
0184 
0185    if str(3)=='T' % complete electrode model
0186       inv_mdl = mk_complete_elec_mdl( inv_mdl, layers);
0187    end
0188       
0189    inv_mdl.fwd_model = thorax_geometry( inv_mdl.fwd_model, str2num(str(4)));
0190 
0191 elseif strcmp( str, 'n3r2')
0192     inv_mdl = mk_n3r2_model( n_elec, options );
0193 elseif strcmp( str, 'n3z') || strcmp(str, 'n3z2')
0194     inv_mdl = mk_n3z_model( n_elec, options );
0195 elseif strcmp( str(2:3), '3c') || strcmp(str(2:3),'3t')
0196    if     str(1)=='a'; xy_layers=  4; z_layers= linspace(-.5,.5,5);
0197    elseif str(1)=='b'; xy_layers=  8; z_layers= linspace(-.7,.7,11);
0198    elseif str(1)=='c'; xy_layers= 12; z_layers= linspace(-.9,.9,21);
0199    elseif str(1)=='d'; xy_layers= 16; z_layers= linspace(-1,1,41);
0200    elseif str(1)=='e'; xy_layers= 20; z_layers= linspace(-1.3,1.3,61);
0201    elseif str(1)=='f'; xy_layers= 24; z_layers= linspace(-1.6,1.6,81);
0202    else;  error(['don`t know what to do with option(1)=',str]);
0203    end
0204 
0205    if str(3)== 'c'
0206       elec_cfg_str= str(4:end);
0207    elseif str(3) == 't'
0208       elec_cfg_str= str(5:end);
0209    else
0210       error(['don`t know what to do with option(3)=',str]);
0211    end
0212 
0213    elec_per_plane = n_elec(1);
0214    spacing=.5;
0215    if     elec_cfg_str=='r';
0216       elec_conf= 'planes';
0217       ne_1  = (n_elec(2)-1)/2;
0218       elec_space= [ne_1:-1:-ne_1]*spacing;
0219    elseif elec_cfg_str=='z2';
0220       elec_conf= 'zigzag';
0221       elec_space= [1,-1]*spacing/2;
0222    elseif elec_cfg_str=='p2';
0223       elec_conf= 'planes';
0224       elec_space= [1,-1]*spacing/2;
0225       elec_per_plane = n_elec(1)/2;
0226    else;
0227       error(['don`t know what to do with option(4)=',str]);
0228    end
0229 
0230    inv_mdl = mk_3c_model( n_elec, xy_layers, z_layers, elec_space, ...
0231                               elec_per_plane, elec_conf, options );
0232 
0233    if str(3) == 't' % thorax models
0234       inv_mdl = rotate_model( inv_mdl, 2); % 45 degrees
0235       inv_mdl.fwd_model = thorax_geometry( inv_mdl.fwd_model, str2num(str(4)));
0236    end
0237 else
0238     error(['Don`t know what to do with option=',str]);
0239 end
0240 
0241 inv_mdl.name= ['EIDORS common_model_',str]; 
0242 inv_mdl= eidors_obj('inv_model', inv_mdl);
0243 
0244 % check we are giving back a good specimen
0245 valid_inv_model(inv_mdl);
0246     
0247 function inv_mdl = distmesh_2d_model(str, n_elec, options);
0248 % This function is an interface to distmesh_2d_model for some common values
0249 % fmdl = dm_2d_circ_pt_elecs( elec_pts, pfix, spacing);
0250 %  params.base_spacing = spacing(1);
0251 %  params.refine_ratio = spacing(2);
0252 %  params.gradient     = spacing(3);
0253    Elec_width= 4; % 2 degrees - electrode width
0254    switch [str(1),str(4)]
0255       case 'j0'; params = [ 55,0,100]./[1000,1,100];
0256       case 'i0'; params = [ 67,0,100]./[1000,1,100];
0257       case 'h0'; params = [ 77,0,100]./[1000,1,100];
0258       case 'g0'; params = [ 87,0,100]./[1000,1,100];
0259       case 'f0'; params = [100,0,100]./[1000,1,100];
0260       case 'e0'; params = [120,0,100]./[1000,1,100];
0261       case 'd0'; params = [150,0,100]./[1000,1,100];
0262       case 'c0'; params = [200,0,100]./[1000,1,100];
0263       case 'b0'; params = [270,0,100]./[1000,1,100];
0264       case 'a0'; params = [500,0,100]./[1000,1,100];
0265 
0266       case 'j1'; params = [ 21,3,5]./[1000,1,100];
0267       case 'i1'; params = [ 23,3,5]./[1000,1,100];
0268       case 'h1'; params = [ 26,3,5]./[1000,1,100];
0269       case 'g1'; params = [ 30,3,5]./[1000,1,100];
0270       case 'f1'; params = [ 35,3,5]./[1000,1,100];
0271       case 'e1'; params = [ 39,3,5]./[1000,1,100];
0272       case 'd1'; params = [ 54,3,5]./[1000,1,100];
0273       case 'c1'; params = [100,3,5]./[1000,1,100];
0274       case 'b1'; params = [180,3,5]./[1000,1,100];
0275       case 'a1'; params = [400,3,5]./[1000,1,100];
0276 
0277       case 'j2'; params = [ 12,5,3]./[1000,1,100];
0278       case 'i2'; params = [ 14,5,3]./[1000,1,100];
0279       case 'h2'; params = [ 16,5,3]./[1000,1,100];
0280       case 'g2'; params = [ 19,5,3]./[1000,1,100];
0281       case 'f2'; params = [ 21,5,3]./[1000,1,100];
0282       case 'e2'; params = [ 39,5,3]./[1000,1,100];
0283       case 'd2'; params = [ 50,5,3]./[1000,1,100];
0284       case 'c2'; params = [100,5,3]./[1000,1,100];
0285       case 'b2'; params = [200,5,3]./[1000,1,100];
0286       case 'a2'; params = [500,5,3]./[1000,1,100];
0287 
0288       case 'j3'; params = [  6,10,2]./[1000,1,100];
0289       case 'i3'; params = [  7,10,2]./[1000,1,100];
0290       case 'h3'; params = [  8,10,2]./[1000,1,100];
0291       case 'g3'; params = [  9,10,2]./[1000,1,100];
0292       case 'f3'; params = [ 10,10,2]./[1000,1,100];
0293       case 'e3'; params = [ 13,10,2]./[1000,1,100];
0294       case 'd3'; params = [ 20,10,2]./[1000,1,100];
0295       case 'c3'; params = [ 70,10,2]./[1000,1,100];
0296       case 'b3'; params = [150,10,2]./[1000,1,100];
0297       case 'a3'; params = [250,10,2]./[1000,1,100];
0298 
0299 % We set refinement 4==3. This is OK for this mdl density
0300       case 'j4'; params = [  6,10,2]./[1000,1,100];
0301       case 'i4'; params = [  7,10,2]./[1000,1,100];
0302       case 'h4'; params = [  8,10,2]./[1000,1,100];
0303       case 'g4'; params = [  9,10,2]./[1000,1,100];
0304       case 'f4'; params = [ 10,10,2]./[1000,1,100];
0305       case 'e4'; params = [ 13,10,2]./[1000,1,100];
0306       case 'd4'; params = [ 20,10,2]./[1000,1,100];
0307       case 'c4'; params = [ 70,10,2]./[1000,1,100];
0308       case 'b4'; params = [150,10,2]./[1000,1,100];
0309       case 'a4'; params = [250,10,2]./[1000,1,100];
0310 
0311       otherwise; error('don`t know what to do with option=%s',str);
0312    end
0313    ea = Elec_width/2 *(2*pi/360);
0314    for i=1:n_elec(1); 
0315      ai = (i-1)/n_elec(1) * 2*pi;
0316      elec_pts{i} = [sin(ai+ea),cos(ai+ea);sin(ai-ea),cos(ai-ea)];
0317    end
0318    fwd_mdl= dm_2d_circ_pt_elecs( elec_pts, [], params);
0319    inv_mdl= add_params_2d_mdl( fwd_mdl, n_elec(1), options);
0320 
0321 % THIS FUNCTION IS DEPRECATED (from EIDORS 3.3)
0322 function inv2d = distmesh_2d_model_depr(str, n_elec, options);
0323    switch str(1)
0324       case 'a'; n_nodes=  50;
0325       case 'b'; n_nodes= 100;
0326       case 'c'; n_nodes= 200;
0327       case 'd'; n_nodes= 400;
0328       case 'e'; n_nodes= 800;
0329       case 'f'; n_nodes=1200;
0330       case 'g'; n_nodes=1800;
0331       case 'h'; n_nodes=2400;
0332       case 'i'; n_nodes=3000;
0333       case 'j'; n_nodes=4000;
0334       otherwise; error('don`t know what to do with option=%s',str);
0335    end
0336  
0337    refine_level= abs(str(4))-'0';
0338 
0339    elec_width= .1;
0340    th=linspace(0,2*pi,n_elec(1)+1)';th(end)=[];
0341    elec_posn= [sin(th),cos(th)];
0342    [elec_nodes, refine_nodes] = dm_mk_elec_nodes( elec_posn, ...
0343           elec_width, refine_level);
0344    fd=inline('sqrt(sum(p.^2,2))-1','p');
0345    bbox = [-1,-1;1,1];
0346    z_contact = 0.01;
0347    fwd_mdl= dm_mk_fwd_model( fd, [], n_nodes, bbox, ...
0348                              elec_nodes, refine_nodes, z_contact);
0349 
0350    inv2d= add_params_2d_mdl( fwd_mdl, n_elec(1), options);
0351 
0352 
0353 function inv2d= mk_2c_model( n_elec, n_circles, options )
0354 
0355     n_elec= n_elec(1);
0356     params= mk_circ_tank(n_circles, [], n_elec); 
0357     inv2d= add_params_2d_mdl( params, n_elec, options);
0358    
0359 
0360 function inv2d= mk_2r_model( n_elec, xy_size, options)
0361     if length(xy_size)==1; xy_size= xy_size*[1,1]; end
0362     xy_size= xy_size+1;
0363 
0364     xvec = linspace(-1,1,xy_size(1));
0365     yvec = linspace(-1,1,xy_size(2));
0366     fmdl = mk_grid_model([],xvec,yvec);
0367 
0368     % put 1/4 of elecs on each side
0369     tb_elecs= linspace(1, xy_size(1), 1+2*n_elec(1)/4); 
0370     tb_elecs= tb_elecs(2:2:end);
0371     sd_elecs= linspace(1, xy_size(2), 1+2*n_elec(1)/4);
0372     sd_elecs= sd_elecs(2:2:end);
0373     
0374     el_nodes= [];
0375     % Top nodes -left to right
0376     bdy_nodes= (1:xy_size(1)) + xy_size(1)*(xy_size(2)-1); 
0377     el_nodes= [el_nodes, bdy_nodes(tb_elecs)];
0378     % Right nodes - top to bottom
0379     bdy_nodes= (1:xy_size(2))*xy_size(1); 
0380     el_nodes= [el_nodes, bdy_nodes(fliplr(sd_elecs))];
0381     % Bottom nodes - right to left
0382     bdy_nodes= 1:xy_size(1); 
0383     el_nodes= [el_nodes, bdy_nodes(fliplr(tb_elecs))];
0384     % Left nodes - bottom to top
0385     bdy_nodes= (0:xy_size(2)-1)*xy_size(1)+1; 
0386     el_nodes= [el_nodes, bdy_nodes(sd_elecs)];
0387 
0388 %   trimesh(fmdl.elems,fmdl.nodes(:,1), fmdl.nodes(:,2));
0389     for i=1:n_elec(1)
0390        n= el_nodes(i);
0391        fmdl.electrode(i).nodes= n;
0392        fmdl.electrode(i).z_contact= .001; % choose a low value
0393 %      plot(fmdl.nodes(n,1),fmdl.nodes(n,2),'*'); pause;
0394     end
0395     inv2d= add_params_2d_mdl( fmdl, n_elec(1), options);
0396 
0397 % params is the part of the fwd_model
0398 function inv2d= add_params_2d_mdl( params, n_elec, options);
0399     n_rings= 1;
0400     [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0401     params.stimulation= st;
0402     params.meas_select= els;
0403     params.solve=      'eidors_default';
0404     params.system_mat= 'eidors_default';
0405     params.jacobian=   'eidors_default';
0406     params.normalize_measurements= 0;
0407     mdl_2d   = eidors_obj('fwd_model', params);
0408 
0409     inv2d.solve=       'eidors_default';
0410     inv2d.hyperparameter.value = 3e-2;
0411     %inv2d.hyperparameter.func = 'choose_noise_figure';
0412     %inv2d.hyperparameter.noise_figure= 1;
0413     %inv2d.hyperparameter.tgt_elems= 1:4;
0414      inv2d.RtR_prior= 'eidors_default';
0415     %inv2d.RtR_prior= 'prior_gaussian_HPF';
0416     inv2d.jacobian_bkgnd.value= 1;
0417     inv2d.reconst_type= 'difference';
0418     inv2d.fwd_model= mdl_2d;
0419 
0420 function inv3d = mk_3c_model( n_elec, xy_layers, z_layers, elec_space, ...
0421                               elec_per_plane, elec_conf, options );
0422 
0423     e_layers=[];
0424     for es= elec_space;
0425        ff= abs(z_layers  -es);
0426        ff= find(ff==min(ff));
0427        e_layers= [e_layers,ff(1)];
0428     end
0429 
0430     params= mk_circ_tank( xy_layers, z_layers, ...
0431            { elec_conf, elec_per_plane, e_layers} );
0432 
0433     [st, els]= mk_stim_patterns(n_elec(1), n_elec(2), '{ad}','{ad}', options, 10);
0434 
0435     params.stimulation= st;
0436     params.meas_select= els;
0437     params.solve=      'eidors_default';
0438     params.system_mat= 'eidors_default';
0439     params.jacobian=   'eidors_default';
0440     params.normalize_measurements= 0;
0441     fm3d = eidors_obj('fwd_model', params);
0442 
0443     inv3d.name=  'EIT inverse: 3D';
0444     inv3d.solve= 'eidors_default';
0445     inv3d.RtR_prior= 'eidors_default';
0446 %     inv3d.inv_solve_time_prior.time_steps=   0;
0447 %     inv3d.prior_time_smooth.space_prior = @prior_noser;
0448 %     inv3d.prior_time_smooth.time_weight = 0;
0449 %     inv3d.inv_solve_time_prior.time_steps   = 0;
0450 
0451     inv3d.hyperparameter.value = 3e-2;
0452     inv3d.reconst_type= 'difference';
0453     inv3d.jacobian_bkgnd.value= 1;
0454     inv3d.fwd_model= fm3d;
0455 
0456 function inv_mdl = mk_n3z_model( n_elec, options );
0457    inv_mdl= mk_n3r2_model( n_elec, options);
0458    fwd_mdl= inv_mdl.fwd_model;
0459    renumber= [1:2:15; 18:2:32];
0460    fwd_mdl.electrode= fwd_mdl.electrode(renumber(:));
0461    n_rings= 1;
0462    [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0463    fwd_mdl.stimulation= st;
0464    fwd_model.meas_select= els;
0465    inv_mdl.fwd_model= fwd_mdl;
0466    inv_mdl.name= 'NP 3D model with zigzag electrodes';
0467 
0468 function inv_mdl = mk_n3r2_model( n_elec, options );
0469    if ~isempty(n_elec) if  ~all(n_elec == [16,2]);
0470       if length(n_elec)~=2; n_elec= [n_elec(1),1]; end
0471       warning(sprintf(['You have requested an "n3" model with [%d,%d] electrodes.' ...
0472          'Note that these models always have 32 electrodes (16x2)'], n_elec));
0473    end; end
0474    if ~exist('OCTAVE_VERSION');
0475       load( 'datareal.mat' );
0476    else
0477       load( file_in_loadpath( 'datareal.mat' ));
0478    end
0479    fmdl.nodes= vtx;
0480    fmdl.elems= simp;
0481    fmdl.boundary= find_boundary( simp );
0482 
0483    fmdl.solve=      @fwd_solve_1st_order;
0484    fmdl.jacobian=   @jacobian_adjoint;
0485    fmdl.system_mat= @system_mat_1st_order;
0486    fmdl.normalize_measurements = 0;
0487 
0488    for i=1:length(zc)
0489        electrodes(i).z_contact= zc(i);
0490        electrodes(i).nodes=     unique( elec(i,:) );
0491    end
0492 
0493    fmdl.gnd_node=           gnd_ind;
0494    fmdl.electrode =         electrodes;
0495 
0496    fmdl.stimulation = mk_stim_patterns(16,2,[0,1],[0,1], ...
0497              {'no_meas_current','no_rotate_meas'},-1);
0498 
0499    fmdl= eidors_obj('fwd_model', fmdl);
0500 
0501    inv_mdl.name=         'Nick Polydorides EIT inverse';
0502    inv_mdl.solve=       @inv_solve_diff_GN_one_step;
0503    inv_mdl.hyperparameter.value = 1e-2;
0504    inv_mdl.RtR_prior= @prior_laplace;
0505    inv_mdl.reconst_type= 'difference';
0506    inv_mdl.jacobian_bkgnd.value= 1;
0507    inv_mdl.fwd_model= fmdl;
0508 
0509 
0510 function inv3d= mk_b3r1_model( n_elec, options )
0511     n_rings= 1;
0512     levels= [-.5:.1:.5]; 
0513     e_levels= 6; 
0514     nr= 8;
0515 
0516     params= mk_circ_tank( nr, levels, { 'planes', n_elec, e_levels } );
0517     [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0518 
0519     params.stimulation= st;
0520     params.meas_select= els;
0521     params.solve=      'fwd_solve_1st_order';
0522     params.system_mat= 'system_mat_1st_order';
0523     params.jacobian=   'jacobian_adjoint';
0524     params.normalize_measurements= 0;
0525     mdl_3d = eidors_obj('fwd_model', params);
0526 
0527     inv3d.name = 'EIT inverse: 3D';
0528     inv3d.solve=       'inv_solve_diff_GN_one_step';
0529     %inv3d.solve=       'aa_inv_conj_grad';
0530     inv3d.hyperparameter.value = 1e-5;
0531     inv3d.RtR_prior= 'prior_laplace';
0532     %inv3d.RtR_prior= 'prior_gaussian_HPF';
0533     inv3d.jacobian_bkgnd.value= 1;
0534     inv3d.reconst_type= 'difference';
0535     inv3d.fwd_model= mdl_3d;
0536 
0537 function inv3d= mk_b3r2_model( n_elec, nr, options )
0538     n_rings= 2;
0539     z_axis = [0:.1:1];
0540     e_levels= [4,8]; 
0541     nr= 4;
0542     n_elec = 8;
0543    
0544     params= mk_circ_tank( nr, z_axis, { 'planes', n_elec, e_levels } );
0545     [st, els]= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', options, 10);
0546 
0547     params.stimulation= st;
0548     params.meas_select= els;
0549     params.solve=      'fwd_solve_1st_order';
0550     params.system_mat= 'system_mat_1st_order';
0551     params.jacobian=   'jacobian_adjoint';
0552     params.normalize_measurements= 0;
0553     mdl_3d = eidors_obj('fwd_model', params);
0554     
0555     % Specify number of levels in mesh for imaging slices
0556     num_levs = length(e_levels);
0557     levels = inf*ones(num_levs,3);
0558     levels(:,3) = e_levels / (length(z_axis)-1);
0559     levels(:,4) = ones(num_levs,1);
0560     levels(:,5) = (1:num_levs)';    
0561     mdl_3d.levels = levels;
0562     
0563     inv3d.name = 'EIT inverse: 3D';
0564     inv3d.solve=       'inv_solve_diff_GN_one_step';
0565     %inv3d.solve=       'aa_inv_conj_grad';
0566     inv3d.hyperparameter.value = 1e-5;
0567     inv3d.RtR_prior= 'prior_laplace';
0568     %inv3d.RtR_prior= 'prior_gaussian_HPF';
0569     inv3d.jacobian_bkgnd.value= 1;
0570     inv3d.reconst_type= 'difference';
0571     inv3d.fwd_model= mdl_3d;
0572 
0573 % rotate model rotate_model/16 times around
0574 function inv_mdl = rotate_model( inv_mdl, rotate_mdl);
0575     inv_mdl = turn_model( inv_mdl, pi/8*rotate_mdl );
0576 
0577     n_elec= length( inv_mdl.fwd_model.electrode );
0578     renum = rem(  rotate_mdl*n_elec/16 + (0:n_elec-1),n_elec)+1;
0579     renum = floor(renum); % round is not quite right for all model
0580                           % cases, but no errors.
0581     inv_mdl.fwd_model.electrode = ...
0582        inv_mdl.fwd_model.electrode( renum);
0583 
0584 % rotate model rotate_model/16 times around
0585 function inv_mdl = turn_model( inv_mdl, angle );
0586     nodes= inv_mdl.fwd_model.nodes;
0587     cos_rot = cos( angle );
0588     sin_rot = sin( angle );
0589     nodes(:,1)= inv_mdl.fwd_model.nodes(:,1:2)*[ cos_rot;-sin_rot];
0590     nodes(:,2)= inv_mdl.fwd_model.nodes(:,1:2)*[ sin_rot; cos_rot];
0591     inv_mdl.fwd_model.nodes= nodes;
0592 
0593 function inv_mdl = mk_complete_elec_mdl( inv_mdl, layers);
0594       inv_mdl = turn_model( inv_mdl, 2*pi/4/layers/2 );
0595 
0596       bdy= inv_mdl.fwd_model.boundary;
0597       for i=1:length(inv_mdl.fwd_model.electrode);
0598          enode= inv_mdl.fwd_model.electrode(i).nodes;
0599          ff= find( enode== bdy(:,1) );
0600          inv_mdl.fwd_model.electrode(i).nodes = bdy(ff,:);
0601       end
0602 
0603 %%%%%%%%%%%%%%%%%%%%%% TESTS %%%%%%%%%%%%%%%%%%%%%%%
0604 function do_unit_test
0605 
0606 test_circ_models
0607 test_3d_models
0608 test_np_get_3d_meas
0609 test_distmesh_models
0610 
0611 function test_np_get_3d_meas
0612    imdl=mk_common_model('n3r2',[16,2]); st1=imdl.fwd_model.stimulation;
0613    fmdl = imdl.fwd_model;
0614    
0615    % get the measurement patterns, only indH is used in this model
0616    %   here we only want to get the meas pattern from 'get_3d_meas',
0617    %   not the voltages, so we enter zeros
0618    load( 'datareal.mat' );
0619    [I,Ib] = set_3d_currents(protocol, elec, ...
0620                fmdl.nodes, fmdl.gnd_node, no_pl);
0621    [jnk,jnk,indH,indV,jnk] = get_3d_meas( elec, ...
0622             fmdl.nodes, zeros(size(I)), Ib, no_pl );
0623    n_elec= size(elec,1);
0624    n_meas= size(indH,1) / size(Ib,2);
0625    for i=1:size(Ib,2)
0626       fmdl.stimulation(i).stimulation= 'Amp';
0627       fmdl.stimulation(i).stim_pattern= Ib(:,i);
0628 
0629       idx= ( 1+ (i-1)*n_meas ):( i*n_meas );
0630       fmdl.stimulation(i).meas_pattern= sparse ( ...
0631               (1:n_meas)'*[1,1], ...
0632               indH( idx, : ), ...
0633               ones(n_meas,2)*[1,0;0,-1], ...
0634               n_meas, n_elec );
0635    end
0636    st2 = fmdl.stimulation;
0637    
0638 
0639    sp1=[]; sp2=[]; mp1=[]; mp2=[]; 
0640    for i=1:32;
0641       sp1= [sp1, st1(i).stim_pattern]; mp1= [mp1, st1(i).meas_pattern];
0642       sp2= [sp2, st2(i).stim_pattern]; mp2= [mp2, st2(i).meas_pattern];
0643    end
0644    unit_test_cmp('STIM_PAT:',sp1,sp2);
0645    unit_test_cmp('MEAS_PAT:',mp1,mp2);
0646      
0647       
0648 
0649 function test_circ_models
0650 % 2D Circular Models
0651 for j=('a'+0):('j'+0)
0652     mk_common_model(sprintf('%c2C2',j),16);
0653     mk_common_model(sprintf('%c2c0',j),16);
0654     mk_common_model(sprintf('%c2t3',j),16);
0655     mk_common_model(sprintf('%c2T4',j),16);
0656 end;
0657 
0658 for j=('a'+0):('f'+0)
0659     mk_common_model(sprintf('%c2s',j),8);
0660 end;
0661 
0662 function test_3d_models
0663 % 3D Models:
0664     mk_common_model('n3r2',[16,2]);
0665  %  mk_common_model('n3z',[16,2]);
0666  
0667     mk_common_model('b3cr',[16,3]);
0668     mk_common_model('b3t2r',[16,1]);
0669     mk_common_model('b3cz2',[16,1]);
0670     mk_common_model('b3cp2',16);
0671  
0672     mk_common_model('a3cr',16);
0673     mk_common_model('b3cr',16);
0674     mk_common_model('c3cr',16);
0675     mk_common_model('d3cr',16);
0676 
0677 function test_distmesh_models
0678 % Distmesh models
0679 for i=0:4; for j=('a'+0):('j'+0)
0680     mk_common_model(sprintf('%c2d%dd',j,i),16);
0681 end; end
0682
mk_common_model

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
