calc_slices

 calc_slices (img, levels, clim  ) show slices at levels of an
             using a fast rendering algorithm
 img    = EIDORS image struct, or array of I images
 levels = any level definition accepted by level_model_slice

 PARAMETERS:
   img.calc_slices.filter % Filter to be applied to images
      Example:    img.calc_slices.filter = ones(3)/9
   img.calc_slices.scale  % Scaling to apply to images
   img.get_img_data.frame_select = which frames of image to display

 rimg= np x np x I x L where np is 128 by default

   np can be adjusted by calc_colours('npoints')
     or by setting
   img.fwd_model.mdl_slice_mapper.{npx, npy}
        see help of mdl_slice_mapper for more options

 See also: LEVEL_MODEL_SLICE

0001 function rimg = calc_slices( img, levels );
0002 % calc_slices (img, levels, clim  ) show slices at levels of an
0003 %             using a fast rendering algorithm
0004 % img    = EIDORS image struct, or array of I images
0005 % levels = any level definition accepted by level_model_slice
0006 %
0007 % PARAMETERS:
0008 %   img.calc_slices.filter % Filter to be applied to images
0009 %      Example:    img.calc_slices.filter = ones(3)/9
0010 %   img.calc_slices.scale  % Scaling to apply to images
0011 %   img.get_img_data.frame_select = which frames of image to display
0012 %
0013 % rimg= np x np x I x L where np is 128 by default
0014 %
0015 %   np can be adjusted by calc_colours('npoints')
0016 %     or by setting
0017 %   img.fwd_model.mdl_slice_mapper.{npx, npy}
0018 %        see help of mdl_slice_mapper for more options
0019 %
0020 % See also: LEVEL_MODEL_SLICE
0021 
0022 % (C) 2006-2022 Andy Adler and Bartek Grychtol.
0023 % License: GPL version 2 or version 3
0024 % $Id: calc_slices.m 6364 2022-05-05 12:55:54Z bgrychtol $
0025 
0026 if ischar(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0027 
0028 np = calc_colours('npoints');
0029 try np = img(1).calc_colours.npoints; end
0030 
0031 if nargin < 2, levels = []; end
0032 
0033 if isfield(img, 'calc_slices') && isfield(img.calc_slices, 'levels')
0034     warning('EIDORS:CALC_SLICES:DeprecatedInterface', ...
0035         ['Ingoring deprecated img.calc_slices.levels definition. '...
0036         'Use img.fwd_model.mdl_slice_mapper.level instead. '...
0037         'See help level_model_slice for valid inputs.'])
0038     img.calc_slices = rmfield(img.calc_slices, 'levels');
0039 end
0040 
0041 
0042 
0043 img = data_mapper(img);
0044 
0045 % Assume all fwd_models are same dimension (all 3D or 2D no mixed dims)
0046 if mdl_dim(img(1))==2 
0047     if nargin>1
0048         eidors_msg('Specified levels ignored for 2D FEM', 4);
0049     end
0050     if isfield(img(1).fwd_model,'mdl_slice_mapper') && isfield(img(1).fwd_model.mdl_slice_mapper, 'level')
0051         eidors_msg('mdl_slice_mapper.level definition ignored for 2D FEM',4); 
0052         for i = 1:numel(img)
0053             img(i).fwd_model.mdl_slice_mapper = rmfield(img(i).fwd_model.mdl_slice_mapper,'level');
0054         end
0055     end
0056     levels= [Inf,Inf,0];
0057 elseif mdl_dim(img(1))==3 && isempty(levels)
0058    levels = [Inf Inf mean(img.fwd_model.nodes(:,3))];
0059    eidors_msg('calc_slices: no levels specified, assuming an xy plane',2);
0060 end
0061 
0062 
0063 rimg = [];
0064 for i=1:length(img)
0065    rimg = cat(3, rimg, calc_this_slice( img(i), levels, np) );
0066 end
0067 
0068 function rimg = calc_this_slice( img, levels, np)
0069     % If scalar levels then we just create that many cut planes on z-dimension
0070     fwd_model = img.fwd_model;
0071     
0072     if ~isfield(fwd_model,'mdl_slice_mapper')
0073         fwd_model.mdl_slice_mapper.npoints = np;
0074     end
0075     if isfield(fwd_model.mdl_slice_mapper, 'level')
0076       eidors_msg(3, 'Using level definition from the mdl_slice_mapper.level field.');
0077       levels = fwd_model.mdl_slice_mapper.level; 
0078     end
0079     
0080     nodes = {fwd_model.nodes};
0081     if size(fwd_model.nodes,2)==3 && size(fwd_model.elems,2) > 3 
0082         nodes = level_model_slice(fwd_model.nodes, levels, 'all');
0083         
0084         % we'll level by replacing nodes. Disable.
0085         lvl = struct('rotation_matrix',eye(3), 'centre', zeros(1,3));
0086         fwd_model.mdl_slice_mapper.level = lvl;
0087     end
0088     
0089     [data, n_images] = get_img_data(img);
0090     
0091     if ~any(isfield(img, {'elem_data', 'node_data'}))
0092        error('img does not have a data field');
0093     end
0094     n_levels = numel(nodes);
0095     for lev_no = fliplr(1:n_levels)
0096         fwd_model.nodes = nodes{lev_no};
0097         if isfield(img,'elem_data')
0098           rimg(:,:,:,lev_no) = calc_image_elems( data, fwd_model);
0099         elseif isfield(img,'node_data')
0100           rimg(:,:,:,lev_no) = calc_image_nodes( data, fwd_model);
0101         end
0102     end
0103     
0104     % FILTER IMAGE
0105     try   
0106         filt = img.calc_slices.filter; 
0107     catch
0108         filt = 1; 
0109     end
0110     try   
0111         scal = img.calc_slices.scale; 
0112     catch
0113         scal = 1; 
0114     end
0115 
0116     filt = scal * ( filt/sum(filt(:)) );
0117     rimg = filter_image(rimg, filt);
0118 
0119 % Calculate an image by mapping it onto the node_ptr matrix
0120 % This makes a blocky image to nearest node -> no longer used
0121 function rimg= calc_image_nearestnodes( node_data, fwd_model)
0122 
0123    node_ptr = mdl_slice_mapper( fwd_model, 'node' );
0124 
0125    backgnd= NaN;
0126    n_images= size(node_data,2);
0127    rval= [backgnd*ones(1,n_images); node_data];
0128    rimg= reshape( rval(node_ptr+1,:), [size(node_ptr), n_images]);
0129 
0130 % Calculate an image by interpolating it onto the elem_ptr matrix
0131 function rimg= calc_image_nodes( node_data, fwd_model)
0132 
0133    nd_interp= mdl_slice_mapper( fwd_model, 'nodeinterp' );
0134    elem_ptr = mdl_slice_mapper( fwd_model, 'elem' );
0135    [sx,sy]= size(elem_ptr);
0136 
0137    node_ptr = fwd_model.elems; node_ptr = [0*node_ptr(1,:);node_ptr];
0138    node_ptr = reshape( node_ptr( elem_ptr+1, :), sx, sy, []);
0139 
0140    n_images = size(node_data,2);
0141    rimg= zeros(sx, sy, n_images);
0142    backgnd= NaN;
0143    
0144    for ni = 1:n_images
0145      znd = [backgnd;node_data(:,ni)]; % add NaN for background
0146      rimg(:,:,ni) = sum( znd(node_ptr+1) .* nd_interp, 3); 
0147    end
0148 
0149 
0150 % Calculate an image by mapping it onto the elem_ptr matrix
0151 function rimg= calc_image_elems( elem_data, fwd_model)
0152 
0153    elem_ptr = mdl_slice_mapper( fwd_model, 'elem');
0154 
0155    backgnd= NaN;
0156    n_images= size(elem_data,2);
0157    rval= backgnd*ones(size(elem_data)+[1,0]);
0158    rval(2:end,:) = elem_data;
0159    rimg= reshape( rval(elem_ptr+1,:), [size(elem_ptr), n_images]);
0160 
0161 
0162 function  rimg = filter_image(rimg, filt);
0163     
0164    %%% Total MATLAB BS )(*&#$)(*#&@
0165    %%% the && operator used to short circuit. Now it doesn't
0166    %%% How the (*&)(*& can anyone take this language seriously
0167    % all(size(filt*scal) == [1,1]) && filt*scal == 1
0168    if all(size(filt)==1) if filt == 1; return; end ; end
0169 
0170    [sz1,sz2,sz3,sz4] = size(rimg);
0171    for j1 = 1:sz3; for j2 = 1:sz4; 
0172       rsl = rimg(:,:,j1,j2);
0173 
0174       rna = isnan(rsl);
0175       rsl(rna) = 0;
0176       rsl = conv2(rsl, filt, 'same');
0177       rsl(rna) = NaN;
0178 
0179       rimg(:,:,j1,j2) = rsl;
0180    end; end
0181 
0182 function do_unit_test
0183    img = mk_image( mk_common_model('a2c2',8));
0184    img.calc_colours.npoints = 8; 
0185 
0186    imc = calc_slices(img);
0187    imt = NaN*ones(8); imt(2:7,2:7) = 1; imt(3:6,:) = 1; imt(:,3:6) = 1; 
0188    unit_test_cmp('cs 2d 1', imc, imt);
0189 
0190    imn = rmfield(img,'elem_data');
0191    imn.node_data = ones(size(imn.fwd_model.nodes,1),1);
0192    imc = calc_slices(imn);
0193    unit_test_cmp('cs 2d 2', imc, imt, 1e-14);
0194 
0195    img.elem_data(1:4) = 2;
0196    imc = calc_slices(img);
0197    imt(4:5,4:5) = 2; 
0198    unit_test_cmp('cs 2d 3', imc, imt);
0199 
0200    imn.node_data(1:5) = 2;
0201    imc = calc_slices(imn);
0202    imt(3:6,3:6) = 1; imt(4:5,3:6) = 1.292893218813452;
0203    imt(3:6,4:5) = 1.292893218813452; imt(4:5,4:5) = 2;
0204    unit_test_cmp('cs 2d 4', imc, imt, 1e-14);
0205 
0206    imn.node_data(:) = 1; imn.node_data(1) = 4;
0207    imc = calc_slices(imn);
0208    imt(3:6,3:6) = 1; imt(4:5,4:5) = 1.878679656440358; 
0209    unit_test_cmp('cs 2d 5', imc, imt, 1e-14);
0210 
0211    imn.calc_colours.npoints = 7; 
0212    imc = calc_slices(imn);
0213    imt = NaN*ones(7); imt(3:5,:) = 1; imt(:,3:5)= 1; imt(2:6,2:6)= 1;imt(4,4) = 4; 
0214    unit_test_cmp('cs 2d 6', imc, imt, 1e-14);
0215 
0216 
0217 % 3D Tests
0218    img = calc_jacobian_bkgnd( mk_common_model('n3r2',[16,2]));
0219    img.calc_colours.npoints = 8; 
0220    imn = calc_slices(img,[inf,inf,1]);
0221 
0222    imt = NaN*ones(8); imt(3:6,:) = 1; imt(:,3:6) = 1; imt(2:7,2:7) = 1; 
0223    unit_test_cmp('cs 3d 1', imn, imt);
0224 
0225    img.calc_colours.npoints = 12; 
0226    imn = calc_slices(img,[inf,0,inf]);
0227    imt = NaN*ones(12); imt(:,3:10) = 1; 
0228    unit_test_cmp('cs 3d 2', imn, imt);
0229 
0230    % Should have no effect
0231    img.fwd_model.nodes(:,3) = img.fwd_model.nodes(:,3)-1;
0232    imn = calc_slices(img,[inf,0,inf]); 
0233    unit_test_cmp('cs 3d 3', imn, imt);
0234 
0235 
0236 % multi image struct
0237    img = mk_image( mk_common_model('a2c2',8));
0238    img.calc_colours.npoints = 8; 
0239    img(2) = img;
0240    imc = calc_slices(img); 
0241    imt = NaN*ones(8); imt(3:6,:) = 1; imt(:,3:6) = 1; imt(2:7,2:7) = 1; 
0242    unit_test_cmp('cs mult 1', imc, cat(3,imt,imt));
0243 
0244    imgb = mk_image( mk_common_model('b2c2',8));
0245    imgb.calc_colours.npoints = 8; 
0246    img(3)=imgb;
0247    imc = calc_slices(img); 
0248    unit_test_cmp('cs mult 2', imc, cat(3,imt,imt,imt));
0249 
0250    imdl = mk_common_model('c2t2',16);
0251    img = mk_image(imdl,1);
0252    imc = calc_slices(img);
0253    unit_test_cmp('size e  1', size(imc), [64,64]);
0254 
0255    img.calc_colours.npoints = 40;
0256    imc = calc_slices(img);
0257    unit_test_cmp('size e  2', size(imc), [40,40]);
0258 
0259    img.fwd_model.mdl_slice_mapper.npx = 22;
0260    img.fwd_model.mdl_slice_mapper.npy = 32;
0261    imc = calc_slices(img);
0262    unit_test_cmp('size e  3', size(imc), [32,22]);
0263 
0264    img.fwd_model = rmfield(img.fwd_model,'mdl_slice_mapper');
0265    img.fwd_model.mdl_slice_mapper.x_pts = linspace(-150,150,20);
0266    img.fwd_model.mdl_slice_mapper.y_pts =-linspace(-150,150,23);
0267    imc = calc_slices(img);
0268    unit_test_cmp('size e  4', size(imc), [23,20]);
0269 
0270    img = rmfield(img,'elem_data');
0271    img.node_data =  ones(size(img.fwd_model.nodes,1),1);
0272    img.node_data =  (1:size(img.fwd_model.nodes,1))';
0273    img.fwd_model = rmfield(img.fwd_model,'mdl_slice_mapper');
0274    imc = calc_slices(img);
0275    unit_test_cmp('size n  1', size(imc), [40,40]);
0276 
0277    im2= img;
0278    im2.node_data =  ones(size(img.fwd_model.nodes,1),5);
0279    imc = calc_slices(im2);
0280    unit_test_cmp('size n x5', size(imc), [40,40,5]);
0281    im2.get_img_data.frame_select = 2:3;
0282    imc = calc_slices(im2);
0283    unit_test_cmp('size n x2', size(imc), [40,40,2]);
0284 
0285    img.fwd_model.mdl_slice_mapper.x_pts = linspace(-150,150,20);
0286    img.fwd_model.mdl_slice_mapper.y_pts =-linspace(-150,150,23);
0287    imc = calc_slices(img);
0288    unit_test_cmp('size n  2', size(imc), [23,20]);