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

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
