set_slice

 SET_SLICE  Set parameters for slice display.
   IMG = SET_SLICE(IMG, LEVEL, SPEC) sets the parameters on the image
   struct IMG the parameters required by MDL_SLICE_MAPPER.

   PARAMETERS:
       IMG     - an eidors image or fwd_model
       LEVEL   - any definition accepted by LEVEL_MODEL_SLICE for 
                 a single slice, or [] to leave current specification
                 unchanged.
       SPEC    - specification of interpolations points, in one of the
                 following forms
           R           - scalar resolution (points per unit)
                         Sets mdl_slice_mapper.resolution.
           [npx, npy]  - 2-element vector specifying horizontal and
                         vertical number of points
                         Sets mdl_slice_mapper.npx/npy.
           [np, 0]
           [0, np]     - 2-element vector, one element is zero. 
                         Sets img.mdl_slice_mapper.npoints,
                         the number of points across the larger.
           {xpts,ypts} - 2-elements cell array containing vectors of
                         points in horizontal and vertical directsion.
                         Sets mdl_slice_mapper.x_pts/y_pts.

   SET_SLICE will remove any alternative specifications of point locations
   from the mdl_slice_mapper field (fields npx/npy, xpts/ypts, resolution,
   npoints). If LEVEL is [], mdl_slice_mapper.level is not changed if
   present, and not added if absent. Otherwise mdl_slice_mapper.level will
   be replaced, and mdl_slice_mapper.center and .rotate (deprecated) will
   be removed.

   SET_SLICE does not modify the mdl_slice_mapper.model_2d field.

   SET_SLICE sets img.show_slices.axes_msm = true. See SHOW_SLICES.

 Example:
       img = set_slice(img, LEVEL, [npx, npy]);
   is equivalent to:
       img.fwd_model.mdl_slice_mapper.level = LEVEL;
       img.fwd_model.mdl_slice_mapper.npx = npx;
       img.fwd_model.mdl_slice_mapper.npy = npy;

 See also MDL_SLICE_MAPPER, LEVEL_MODEL_SLICE, SHOW_SLICES

0001 function img = set_slice(img, level, spec)
0002 % SET_SLICE  Set parameters for slice display.
0003 %   IMG = SET_SLICE(IMG, LEVEL, SPEC) sets the parameters on the image
0004 %   struct IMG the parameters required by MDL_SLICE_MAPPER.
0005 %
0006 %   PARAMETERS:
0007 %       IMG     - an eidors image or fwd_model
0008 %       LEVEL   - any definition accepted by LEVEL_MODEL_SLICE for
0009 %                 a single slice, or [] to leave current specification
0010 %                 unchanged.
0011 %       SPEC    - specification of interpolations points, in one of the
0012 %                 following forms
0013 %           R           - scalar resolution (points per unit)
0014 %                         Sets mdl_slice_mapper.resolution.
0015 %           [npx, npy]  - 2-element vector specifying horizontal and
0016 %                         vertical number of points
0017 %                         Sets mdl_slice_mapper.npx/npy.
0018 %           [np, 0]
0019 %           [0, np]     - 2-element vector, one element is zero.
0020 %                         Sets img.mdl_slice_mapper.npoints,
0021 %                         the number of points across the larger.
0022 %           {xpts,ypts} - 2-elements cell array containing vectors of
0023 %                         points in horizontal and vertical directsion.
0024 %                         Sets mdl_slice_mapper.x_pts/y_pts.
0025 %
0026 %   SET_SLICE will remove any alternative specifications of point locations
0027 %   from the mdl_slice_mapper field (fields npx/npy, xpts/ypts, resolution,
0028 %   npoints). If LEVEL is [], mdl_slice_mapper.level is not changed if
0029 %   present, and not added if absent. Otherwise mdl_slice_mapper.level will
0030 %   be replaced, and mdl_slice_mapper.center and .rotate (deprecated) will
0031 %   be removed.
0032 %
0033 %   SET_SLICE does not modify the mdl_slice_mapper.model_2d field.
0034 %
0035 %   SET_SLICE sets img.show_slices.axes_msm = true. See SHOW_SLICES.
0036 %
0037 % Example:
0038 %       img = set_slice(img, LEVEL, [npx, npy]);
0039 %   is equivalent to:
0040 %       img.fwd_model.mdl_slice_mapper.level = LEVEL;
0041 %       img.fwd_model.mdl_slice_mapper.npx = npx;
0042 %       img.fwd_model.mdl_slice_mapper.npy = npy;
0043 %
0044 % See also MDL_SLICE_MAPPER, LEVEL_MODEL_SLICE, SHOW_SLICES
0045 
0046 % (C) 2024 Bartek Grychtol.
0047 % License: GPL version 2 or version 3
0048 % $Id: set_slice.m 6799 2024-04-20 19:44:10Z aadler $
0049 
0050 if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0051 
0052 switch img.type
0053     case 'image'
0054         fmdl = img.fwd_model;
0055         isimage = true;
0056     case 'fwd_model'
0057         fmdl = img;
0058         isimage = false;
0059     otherwise
0060         error('EIDORS:WrongInput', 'Expected eidors image or fwd_model.')
0061 end
0062 
0063 fmdl = do_set_slice(fmdl, level, spec);
0064 
0065 if isimage
0066     img.fwd_model = fmdl;
0067     img.show_slices.axes_msm = true;
0068 else
0069     img = fmdl;
0070 end
0071 
0072 
0073 function fmdl = do_set_slice(fmdl, level, spec)
0074 
0075 if ~isempty(level)
0076     fmdl.mdl_slice_mapper.level = level;
0077     try fmdl.mdl_slice_mapper = rmfield(fmdl.mdl_slice_mapper, {'centre', 'rotate'}); end
0078 end
0079 
0080 if ~isempty(spec)
0081     try fmdl.mdl_slice_mapper = rmfield(fmdl.mdl_slice_mapper, {'x_pts', 'y_pts'}); end
0082     try fmdl.mdl_slice_mapper = rmfield(fmdl.mdl_slice_mapper, {'npx', 'npy'}); end
0083     try fmdl.mdl_slice_mapper = rmfield(fmdl.mdl_slice_mapper, 'npoints'); end
0084     try fmdl.mdl_slice_mapper = rmfield(fmdl.mdl_slice_mapper, 'resolution'); end
0085 end
0086 
0087 if iscell(spec) % {xpts, ypts}
0088     fmdl.mdl_slice_mapper.x_pts = spec{1};
0089     fmdl.mdl_slice_mapper.y_pts = spec{2};
0090 elseif isnumeric(spec)
0091     switch numel(spec)
0092         case 1 % resolution
0093             fmdl.mdl_slice_mapper.resolution = spec;
0094         case 2 % [npx, npy] or [np, 0]
0095             if any(spec == 0)
0096                 fmdl.mdl_slice_mapper.npoints = max(spec);
0097             else
0098                 fmdl.mdl_slice_mapper.npx = spec(1);
0099                 fmdl.mdl_slice_mapper.npy = spec(2);
0100             end
0101         otherwise
0102             error('EIDORS:WrongInput', 'Expected 1- or 2- element SPEC');
0103             
0104     end      
0105 else
0106     error('EIDORS:WrongInput','Point definition SPEC not understood.');
0107 end
0108 
0109 
0110 
0111 function do_unit_test()
0112     figure(1)
0113     show_slices('UNIT_TEST');
0114     figure(2)
0115     do_unit_test_show_slices()
0116     figure(3);
0117     mdl_slice_mapper('UNIT_TEST');
0118     figure(4);
0119     do_unit_test_mdl_slice_mapper()
0120     
0121 
0122 function do_unit_test_show_slices()
0123    clf; sp=0;
0124    %1
0125    img=calc_jacobian_bkgnd(mk_common_model('a2t3',8)); 
0126    img.elem_data=toeplitz(1:size(img.fwd_model.elems,1),1);
0127    sp=sp+1;subplot(4,5,sp); show_slices(img) 
0128    %2
0129    %img.calc_colours.npoints= 128;
0130    img = set_slice(img, [], [128, 0]);
0131    sp=sp+1;subplot(4,5,sp); show_slices(img) 
0132    %3
0133 %    img.calc_colours.npoints= 32;
0134    img = set_slice(img, [], [0 32]);
0135    img.elem_data=toeplitz(1:size(img.fwd_model.elems,1),1:3);
0136    sp=sp+1;subplot(4,5,sp); show_slices(img) 
0137 
0138    %4
0139    img.show_slices.img_cols= 1;
0140    sp=sp+1;subplot(4,5,sp); show_slices(img) 
0141 
0142    %5
0143    imgn = rmfield(img,'elem_data');
0144    imgn.node_data=toeplitz(1:size(img.fwd_model.nodes,1),1);
0145 
0146    img.elem_data = img.elem_data(:,1);
0147 %    img.fwd_model.mdl_slice_mapper.npx = 10;
0148 %    img.fwd_model.mdl_slice_mapper.npy = 20;
0149 %    img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0150    img = set_slice(img, [inf,inf,0], [10, 20]);
0151    sp=sp+1;subplot(4,5,sp); show_slices(img);
0152 
0153    %6
0154    img.elem_data = img.elem_data(:,1);
0155 %    img.fwd_model.mdl_slice_mapper = rmfield(img.fwd_model.mdl_slice_mapper, {'npx','npy'});
0156 %    img.fwd_model.mdl_slice_mapper.x_pts = linspace(-100,100,20);
0157 %    img.fwd_model.mdl_slice_mapper.y_pts = linspace(-150,150,30);
0158 %    img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0159    img = set_slice(img,[inf,inf,0], {linspace(-100,100,20),linspace(-150,150,30)}); 
0160    sp=sp+1;subplot(4,5,sp); show_slices(img);
0161 
0162    %7
0163    sp=sp+1;subplot(4,5,sp); show_slices(imgn) 
0164 
0165    %8
0166 %    imgn.fwd_model.mdl_slice_mapper.x_pts = linspace(-100,100,20);
0167 %    imgn.fwd_model.mdl_slice_mapper.y_pts = linspace(-150,150,30);
0168 %    imgn.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0169    imgn = set_slice(imgn,[inf,inf,0], {linspace(-100,100,20),linspace(-150,150,30)}); 
0170    sp=sp+1;subplot(4,5,sp); show_slices(imgn) 
0171 
0172 
0173 % 3D images
0174    %9
0175    img=calc_jacobian_bkgnd(mk_common_model('n3r2',[16,2])); 
0176 %    img.calc_colours.npoints= 16;
0177    img = set_slice(img, [], [0 16]);
0178    img.elem_data=toeplitz(1:size(img.fwd_model.elems,1),1);
0179    sp=sp+1;subplot(4,5,sp); show_slices(img,2) 
0180 
0181    %10
0182    img.elem_data=img.elem_data*[1:3];
0183    sp=sp+1;subplot(4,5,sp); show_slices(img,2) 
0184 
0185    %11
0186    img.elem_data=img.elem_data(:,1:2);
0187    sp=sp+1;subplot(4,5,sp); show_slices(img,[inf,inf,1;0,inf,inf;0,1,inf]);
0188 
0189    %12
0190    img.show_slices.sep = 5;
0191 %    img.fwd_model.mdl_slice_mapper.x_pts = linspace(-1,1,20);
0192 %    img.fwd_model.mdl_slice_mapper.y_pts = linspace(-1,1,30);
0193 %    img.fwd_model.mdl_slice_mapper.level = [inf,inf,0];
0194    img = set_slice(img, [inf,inf,0], {linspace(-1,1,20), linspace(-1,1,30)});
0195    sp=sp+1;subplot(4,5,sp); show_slices(img,2) 
0196 
0197    % tests 13--19 of calc_slices are not relevant
0198 
0199 function do_unit_test_mdl_slice_mapper()
0200 % 2D NUMBER OF POINTS
0201    imdl = mk_common_model('a2c2',8); fmdl = imdl.fwd_model;
0202    fmdl.nodes = 1e-15 * round(1e15*fmdl.nodes);
0203 %    fmdl.mdl_slice_mapper.level = [inf,inf,0];
0204 %    fmdl.mdl_slice_mapper.npx = 5;
0205 %    fmdl.mdl_slice_mapper.npy = 5;
0206    fmdl = set_slice(fmdl, [inf,inf,0], [5, 5]);
0207    eptr = mdl_slice_mapper(fmdl,'elem');
0208    
0209    si = @(x,y) sub2ind([5,5],x, y);
0210    % points lie on nodes and edges, we allow any associated element
0211    els = {si(2,2), [25,34];
0212           si(2,4), [27,30];
0213           si(3,3), [1,2,3,4];
0214           si(4,2), [33,36];
0215           si(4,4), [28,31];
0216           };
0217    
0218    tst = eptr == [    0   37   38   39    0
0219                      46   25    5   27   42
0220                      47    8    1    6   41
0221                      48   33    7   28   40
0222                       0   45   44   43    0];
0223                       
0224    for i = 1:size(els,1)
0225     tst(els{i,1}) = ismember(eptr(els{i,1}), els{i,2});
0226    end
0227    unit_test_cmp('eptr01',tst,true(5,5));
0228 
0229    nptr = mdl_slice_mapper(fmdl,'node');
0230    test = [   0   27   28   29    0
0231              41    6    7    8   31
0232              40   13    1    9   32
0233              39   12   11   10   33
0234               0   37   36   35    0];
0235    unit_test_cmp('nptr01',nptr,test);
0236 
0237    nint = mdl_slice_mapper(fmdl,'nodeinterp');
0238    unit_test_cmp('nint01a',nint(2:4,2:4,1),[ 0.2627, 0.6909, 0.2627; 
0239                                              0.6906, 1,      0.6906;
0240                                              0.2627, 0.6909, 0.2627], 1e-3);
0241 
0242 %    fmdl.mdl_slice_mapper.npx = 6;
0243 %    fmdl.mdl_slice_mapper.npy = 4;
0244    fmdl = set_slice(fmdl, [], [6 4]);
0245    eptr = mdl_slice_mapper(fmdl,'elem');
0246    res = [  0   49   38   38   52    0
0247            62   23    9   10   20   55
0248            63   24   14   13   19   54
0249             0   60   44   44   57    0];
0250    unit_test_cmp('eptr02',eptr,res);
0251 
0252    nptr = mdl_slice_mapper(fmdl,'node');
0253    res = [  0   27   15   16   29    0
0254            25    6    2    3    8   18
0255            24   12    5    4   10   19
0256             0   37   22   21   35    0];
0257    unit_test_cmp('nptr02',nptr,res);
0258 
0259 % DIRECT POINT TESTS
0260    imdl = mk_common_model('a2c2',8); fmdl = imdl.fwd_model;
0261 %    fmdl.mdl_slice_mapper.level = [inf,inf,0];
0262 %    fmdl.mdl_slice_mapper.x_pts = linspace(-.95,.95,4);
0263 %    fmdl.mdl_slice_mapper.y_pts = [0,0.5];
0264    fmdl = set_slice(fmdl, [inf,inf,0], {linspace(-.95,.95,4), [0,0.5]});
0265    eptr = mdl_slice_mapper(fmdl,'elem');
0266    unit_test_cmp('eptr03',eptr,[ 47 8 6 41; 0 25 27 0]);
0267 
0268    nptr = mdl_slice_mapper(fmdl,'node');
0269    unit_test_cmp('nptr03',nptr,[ 40 13 9 32; 0 6 8 0]);
0270 
0271 % 3D NPOINTS
0272    imdl = mk_common_model('n3r2',[16,2]); fmdl = imdl.fwd_model;
0273 %    fmdl.mdl_slice_mapper.level = [inf,inf,1]; % co-planar with faces
0274 %    fmdl.mdl_slice_mapper.npx = 4;
0275 %    fmdl.mdl_slice_mapper.npy = 4;
0276    fmdl = set_slice(fmdl,[inf,inf,1],[4, 4]);
0277    eptr = mdl_slice_mapper(fmdl,'elem');
0278    si = @(x,y) sub2ind([4,4],x, y);
0279    % points lie on nodes and edges, we allow any associated element
0280    els = {si(1,2), [ 42,317];
0281           si(1,3), [ 30,305];
0282           si(2,1), [ 66,341];
0283           si(2,2), [243,518];
0284           si(2,3), [231,506];
0285           si(2,4), [  6,281];
0286           si(3,1), [ 78,353];
0287           si(3,2), [252,527];
0288           si(3,3), [264,539];
0289           si(3,4), [138,413];
0290           si(4,2), [102,377];
0291           si(4,3), [114,389];};
0292    tst = eptr == zeros(4); 
0293    for i = 1:size(els,1)
0294     tst(els{i,1}) = ismember(eptr(els{i,1}), els{i,2});
0295    end
0296    unit_test_cmp('eptr04',tst, true(4));
0297    nptr = mdl_slice_mapper(fmdl,'node');
0298    test = [   0   101    99     0
0299             103   116   113    97
0300             105   118   121   111
0301               0   107   109     0];
0302    unit_test_cmp('nptr04',nptr, test);
0303 
0304    fmdl.mdl_slice_mapper.level = [inf,0.01,inf];
0305    eptr = mdl_slice_mapper(fmdl,'elem'); 
0306    test = [621   792   777   555
0307            345   516   501   279
0308            343   515   499   277
0309             69   240   225     3];
0310    unit_test_cmp('eptr05',eptr,test);
0311 
0312    nptr = mdl_slice_mapper(fmdl,'node');
0313    test = [230   250   248   222
0314            167   187   185   159
0315            104   124   122    96
0316             41    61    59    33];
0317    unit_test_cmp('nptr05',nptr,test);
0318 
0319    nint = mdl_slice_mapper(fmdl,'nodeinterp');
0320    unit_test_cmp('nint05a',nint(2:3,2:3,4),[.1250,.1250;.8225,.8225],1e-3);
0321    
0322 % Centre and Rotate
0323    fmdl.mdl_slice_mapper = rmfield(fmdl.mdl_slice_mapper,'level');
0324    fmdl.mdl_slice_mapper.centre = [0,0,0.9];
0325    fmdl.mdl_slice_mapper.rotate = eye(3);
0326 %    fmdl.mdl_slice_mapper.npx = 4;
0327 %    fmdl.mdl_slice_mapper.npy = 4;
0328    fmdl = set_slice(fmdl, [], [4,4]);
0329    eptr = mdl_slice_mapper(fmdl,'elem');
0330    test = [  0    42    30     0
0331             66   243   229     6
0332             78   250   264   138
0333              0   102   114     0];
0334    unit_test_cmp('eptr06',eptr, test);
0335 
0336 % SLOW
0337    imdl = mk_common_model('d3cr',[16,3]); fmdl = imdl.fwd_model;
0338    fmdl.nodes = 1e-15*round(1e15*fmdl.nodes);
0339 %    fmdl.mdl_slice_mapper.level = [inf,inf,1];
0340 %    fmdl.mdl_slice_mapper.npx = 64;
0341 %    fmdl.mdl_slice_mapper.npy = 64;
0342    fmdl = set_slice(fmdl, [inf,inf,1], [64, 64]);
0343    t = cputime;
0344    eptr = mdl_slice_mapper(fmdl,'elem');
0345    txt = sprintf('eptr10 (t=%5.3fs)',cputime - t);
0346 % Note that triangulation gives different
0347 % results near the edge
0348    test = [122872   122872   122748
0349            122809   122749   122689
0350            122749   122749   122689];
0351    unit_test_cmp(txt,eptr(10:12,11:13),test);  
0352    
0353 
0354 % CHECK ORIENTATION
0355    imdl=mk_common_model('a2c0',16);
0356    img= mk_image(imdl,1); img.elem_data(26)=1.2;
0357    subplot(231);show_fem(img);
0358    subplot(232);show_slices(img);
0359 %    img.fwd_model.mdl_slice_mapper.npx= 20;
0360 %    img.fwd_model.mdl_slice_mapper.npy= 30;
0361 %    img.fwd_model.mdl_slice_mapper.level= [inf,inf,0];
0362    img = set_slice(img, [inf,inf,0], [20, 30]);
0363    
0364 % Remaining tests in mdl_slice_mapper are not relevant to set_slice