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 = Matrix [Lx3] of L image levels
          each row of the matrix specifies the intercepts
          of the slice on the x, y, z axis. To specify a z=2 plane
          parallel to the x,y: use levels= [inf,inf,2]

 if levels is scalar, then make levels equispaced horizontal
          cuts through the object

 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.calc_slices.levels % an alternative way to specify levels
   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

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