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-2024 Andy Adler and Bartek Grychtol.
0023 % License: GPL version 2 or version 3
0024 % $Id: calc_slices.m 6964 2024-09-27 13:30:40Z aadler $
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 end
0058 
0059 
0060 warnStruct = warning('query','EIDORS:calc_slices:level_conflict');
0061 rimg = [];
0062 for i=1:length(img)
0063    rimg = cat(3, rimg, calc_this_slice( img(i), levels, np) );
0064 end
0065 warning(warnStruct)
0066 
0067 function rimg = calc_this_slice( img, levels, np)
0068     % If scalar levels then we just create that many cut planes on z-dimension
0069     fwd_model = img.fwd_model;
0070     
0071     if ~isfield(fwd_model,'mdl_slice_mapper') || ...
0072           nnz(isfield(fwd_model.mdl_slice_mapper, {'x_pts', 'npoints', 'npx', 'resolution'})) == 0
0073 
0074         fwd_model.mdl_slice_mapper.npoints = np;
0075     end
0076     if isfield(fwd_model.mdl_slice_mapper, 'level')
0077        if ~isempty(levels)
0078           warning('EIDORS:calc_slices:level_conflict',...
0079              'Ignoring provided level definition. Using mdl_slice_mapper.level instead.');
0080           % only warn once
0081           warning('off', 'EIDORS:calc_slices:level_conflict'); 
0082        else
0083           eidors_msg('Using level definition from the mdl_slice_mapper.level field.', 3);
0084        end
0085        levels = fwd_model.mdl_slice_mapper.level;
0086     end
0087     if isempty(levels)
0088        z = max(img.fwd_model.nodes(:,3)) - min(img.fwd_model.nodes(:,3));
0089        levels = [Inf Inf z/2];
0090        eidors_msg('calc_slices: no levels specified, using mid-z xy plane',2);
0091     end
0092 
0093     
0094     nodes = {fwd_model.nodes};
0095     if size(fwd_model.nodes,2)==3 && size(fwd_model.elems,2) > 3 
0096         nodes = level_model_slice(fwd_model.nodes, levels, 'all');
0097         
0098         % we'll level by replacing nodes. Disable.
0099         lvl = struct('rotation_matrix',eye(3), 'centre', zeros(1,3));
0100         fwd_model.mdl_slice_mapper.level = lvl;
0101     end
0102     
0103     [data, n_images] = get_img_data(img);
0104     
0105     if ~any(isfield(img, {'elem_data', 'node_data'}))
0106        error('img does not have a data field');
0107     end
0108     n_levels = numel(nodes);
0109     for lev_no = fliplr(1:n_levels)
0110         fwd_model.nodes = nodes{lev_no};
0111         if isfield(img,'elem_data')
0112           rimg(:,:,:,lev_no) = calc_image_elems( data, fwd_model);
0113         elseif isfield(img,'node_data')
0114           rimg(:,:,:,lev_no) = calc_image_nodes( data, fwd_model);
0115         end
0116     end
0117     
0118     % FILTER IMAGE
0119     try   
0120         filt = img.calc_slices.filter; 
0121     catch
0122         filt = 1; 
0123     end
0124     try   
0125         scal = img.calc_slices.scale; 
0126     catch
0127         scal = 1; 
0128     end
0129 
0130     filt = scal * ( filt/sum(filt(:)) );
0131     rimg = filter_image(rimg, filt);
0132 
0133 % Calculate an image by mapping it onto the node_ptr matrix
0134 % This makes a blocky image to nearest node -> no longer used
0135 function rimg= calc_image_nearestnodes( node_data, fwd_model)
0136 
0137    node_ptr = mdl_slice_mapper( fwd_model, 'node' );
0138 
0139    backgnd= NaN;
0140    n_images= size(node_data,2);
0141    rval= [backgnd*ones(1,n_images); node_data];
0142    rimg= reshape( rval(node_ptr+1,:), [size(node_ptr), n_images]);
0143 
0144 % Calculate an image by interpolating it onto the elem_ptr matrix
0145 function rimg= calc_image_nodes( node_data, fwd_model)
0146 
0147    nd_interp= mdl_slice_mapper( fwd_model, 'nodeinterp' );
0148    elem_ptr = mdl_slice_mapper( fwd_model, 'elem' );
0149    [sx,sy]= size(elem_ptr);
0150 
0151    node_ptr = fwd_model.elems; node_ptr = [0*node_ptr(1,:);node_ptr];
0152    node_ptr = reshape( node_ptr( elem_ptr+1, :), sx, sy, []);
0153 
0154    n_images = size(node_data,2);
0155    rimg= zeros(sx, sy, n_images);
0156    backgnd= NaN;
0157    
0158    for ni = 1:n_images
0159      znd = [backgnd;node_data(:,ni)]; % add NaN for background
0160      rimg(:,:,ni) = sum( znd(node_ptr+1) .* nd_interp, 3); 
0161    end
0162 
0163 
0164 % Calculate an image by mapping it onto the elem_ptr matrix
0165 function rimg= calc_image_elems( elem_data, fwd_model)
0166 
0167    elem_ptr = mdl_slice_mapper( fwd_model, 'elem');
0168 
0169    backgnd= NaN;
0170    n_images= size(elem_data,2);
0171    rval= backgnd*ones(size(elem_data)+[1,0]);
0172    rval(2:end,:) = elem_data;
0173    rimg= reshape( rval(elem_ptr+1,:), [size(elem_ptr), n_images]);
0174 
0175 
0176 function  rimg = filter_image(rimg, filt);
0177     
0178    %%% Total MATLAB BS )(*&#$)(*#&@
0179    %%% the && operator used to short circuit. Now it doesn't
0180    %%% How the (*&)(*& can anyone take this language seriously
0181    % all(size(filt*scal) == [1,1]) && filt*scal == 1
0182    if all(size(filt)==1) if filt == 1; return; end ; end
0183 
0184    [sz1,sz2,sz3,sz4] = size(rimg);
0185    for j1 = 1:sz3; for j2 = 1:sz4; 
0186       rsl = rimg(:,:,j1,j2);
0187 
0188       rna = isnan(rsl);
0189       rsl(rna) = 0;
0190       rsl = conv2(rsl, filt, 'same');
0191       rsl(rna) = NaN;
0192 
0193       rimg(:,:,j1,j2) = rsl;
0194    end; end
0195 
0196 function do_unit_test
0197    img = mk_image( mk_common_model('a2c2',8));
0198    img.calc_colours.npoints = 8; 
0199 
0200    imc = calc_slices(img);
0201    imt = NaN*ones(8); imt(2:7,2:7) = 1; imt(3:6,:) = 1; imt(:,3:6) = 1; 
0202    unit_test_cmp('cs 2d 1', imc, imt);
0203 
0204    imn = rmfield(img,'elem_data');
0205    imn.node_data = ones(size(imn.fwd_model.nodes,1),1);
0206    imc = calc_slices(imn);
0207    unit_test_cmp('cs 2d 2', imc, imt, 1e-14);
0208 
0209    img.elem_data(1:4) = 2;
0210    imc = calc_slices(img);
0211    imt(4:5,4:5) = 2; 
0212    unit_test_cmp('cs 2d 3', imc, imt);
0213 
0214    imn.node_data(1:5) = 2;
0215    imc = calc_slices(imn);
0216    imt(3:6,3:6) = 1; imt(4:5,3:6) = 1.292893218813452;
0217    imt(3:6,4:5) = 1.292893218813452; imt(4:5,4:5) = 2;
0218    unit_test_cmp('cs 2d 4', imc, imt, 1e-14);
0219 
0220    imn.node_data(:) = 1; imn.node_data(1) = 4;
0221    imc = calc_slices(imn);
0222    imt(3:6,3:6) = 1; imt(4:5,4:5) = 1.878679656440358; 
0223    unit_test_cmp('cs 2d 5', imc, imt, 1e-14);
0224 
0225    imn.calc_colours.npoints = 7; 
0226    imc = calc_slices(imn);
0227    imt = NaN*ones(7); imt(3:5,:) = 1; imt(:,3:5)= 1; imt(2:6,2:6)= 1;imt(4,4) = 4; 
0228    unit_test_cmp('cs 2d 6', imc, imt, 1e-14);
0229 
0230 
0231 % 3D Tests
0232    img = calc_jacobian_bkgnd( mk_common_model('n3r2',[16,2]));
0233    xyz = interp_mesh(img);
0234    img.calc_colours.npoints = 8; 
0235    imn = calc_slices(img,[inf,inf,1]);
0236 
0237    imt = NaN*ones(8); imt(3:6,:) = 1; imt(:,3:6) = 1; imt(2:7,2:7) = 1; 
0238    unit_test_cmp('cs 3d 1', imn, imt);
0239 
0240    imgx= img;
0241    imgx.elem_data(xyz(:,3)>1.5) = 2;
0242    imn = calc_slices(imgx,[inf,inf,0.5;inf,inf,2.5]);
0243    unit_test_cmp('cs 3d 2', imn, cat(4,imt,2*imt));
0244 
0245    imgx.elem_data = imgx.elem_data*[1,2,3];
0246    imn = calc_slices(imgx,[inf,inf,0.5]);
0247    imtx = cat(3,imt,2*imt,3*imt);
0248    unit_test_cmp('cs 3d 3', imn, imtx);
0249 
0250    imn = calc_slices(imgx,[inf,inf,0.5;inf,inf,2.5]);
0251    unit_test_cmp('cs 3d 3', imn, cat(4,imtx,2*imtx));
0252 
0253    imgx.fwd_model.mdl_slice_mapper=struct('x_pts',(-1:1)/2,'y_pts',(-1:1)/2, ...
0254          'level',[inf,inf,0.5;inf,inf,2.5]);
0255    imn = calc_slices(imgx);
0256    imt = ones(3); imtx = cat(3,imt,2*imt,3*imt);
0257    unit_test_cmp('cs 3d 3', imn, cat(4,imtx,2*imtx));
0258 
0259 
0260    img.calc_colours.npoints = 12; 
0261    imn = calc_slices(img,[inf,0,inf]);
0262    imt = NaN*ones(12); imt(:,3:10) = 1; 
0263    unit_test_cmp('cs 3d 4', imn, imt);
0264 
0265    % Should have no effect
0266    img.fwd_model.nodes(:,3) = img.fwd_model.nodes(:,3)-1;
0267    imn = calc_slices(img,[inf,0,inf]); 
0268    unit_test_cmp('cs 3d 5', imn, imt);
0269    
0270 
0271    % Default level
0272    imt = calc_slices(img, [inf inf 0.5]);
0273    imn = calc_slices(img);
0274    unit_test_cmp('default', imn, imt);
0275 
0276     % Default level
0277    img.fwd_model.mdl_slice_mapper.level = [inf inf 0.5];
0278    imn = calc_slices(img);
0279    unit_test_cmp('no message', imn, imt);  
0280 
0281    % Warning test
0282    warning('off', 'EIDORS:calc_slices:level_conflict');  %don't need to print
0283    [old, old_id] = lastwarn;
0284    lastwarn('','');
0285    img.fwd_model.mdl_slice_mapper.level = [inf,inf, 1.5];
0286    imt = calc_slices(img,[inf,inf, .5]); 
0287    [~,id] = lastwarn;
0288    lastwarn(old, old_id);
0289    warning('on', 'EIDORS:calc_slices:level_conflict'); 
0290    unit_test_cmp('warning', id, 'EIDORS:calc_slices:level_conflict');
0291 
0292 
0293 % multi image struct
0294    img = mk_image( mk_common_model('a2c2',8));
0295    img.calc_colours.npoints = 8; 
0296    img(2) = img;
0297    imc = calc_slices(img); 
0298    imt = NaN*ones(8); imt(3:6,:) = 1; imt(:,3:6) = 1; imt(2:7,2:7) = 1; 
0299    unit_test_cmp('cs mult 1', imc, cat(3,imt,imt));
0300 
0301    imgb = mk_image( mk_common_model('b2c2',8));
0302    imgb.calc_colours.npoints = 8; 
0303    img(3)=imgb;
0304    imc = calc_slices(img); 
0305    unit_test_cmp('cs mult 2', imc, cat(3,imt,imt,imt));
0306 
0307    imdl = mk_common_model('c2t2',16);
0308    img = mk_image(imdl,1);
0309    imc = calc_slices(img);
0310    unit_test_cmp('size e  1', size(imc), [64,64]);
0311 
0312    img.calc_colours.npoints = 40;
0313    imc = calc_slices(img);
0314    unit_test_cmp('size e  2', size(imc), [40,40]);
0315 
0316    img.fwd_model.mdl_slice_mapper.npx = 22;
0317    img.fwd_model.mdl_slice_mapper.npy = 32;
0318    imc = calc_slices(img);
0319    unit_test_cmp('size e  3', size(imc), [32,22]);
0320 
0321    img.fwd_model = rmfield(img.fwd_model,'mdl_slice_mapper');
0322    img.fwd_model.mdl_slice_mapper.x_pts = linspace(-150,150,20);
0323    img.fwd_model.mdl_slice_mapper.y_pts =-linspace(-150,150,23);
0324    imc = calc_slices(img); 
0325    unit_test_cmp('size e  4', size(imc), [23,20]);
0326 
0327    img = rmfield(img,'elem_data');
0328    img.node_data =  ones(size(img.fwd_model.nodes,1),1);
0329    img.node_data =  (1:size(img.fwd_model.nodes,1))';
0330    img.fwd_model = rmfield(img.fwd_model,'mdl_slice_mapper');
0331    imc = calc_slices(img);
0332    unit_test_cmp('size n  1', size(imc), [40,40]);
0333 
0334    im2= img;
0335    im2.node_data =  ones(size(img.fwd_model.nodes,1),5);
0336    imc = calc_slices(im2);
0337    unit_test_cmp('size n x5', size(imc), [40,40,5]);
0338    im2.get_img_data.frame_select = 2:3;
0339    imc = calc_slices(im2);
0340    unit_test_cmp('size n x2', size(imc), [40,40,2]);
0341 
0342    img.fwd_model.mdl_slice_mapper.x_pts = linspace(-150,150,20);
0343    img.fwd_model.mdl_slice_mapper.y_pts =-linspace(-150,150,23);
0344    imc = calc_slices(img);
0345    unit_test_cmp('size n  2', size(imc), [23,20]);