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

 rimg= np x np x I x L where np is 128 by default
 np can be adjusted by calc_colours('npoints')

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 %
0019 % rimg= np x np x I x L where np is 128 by default
0020 % np can be adjusted by calc_colours('npoints')
0021 
0022 % (C) 2006 Andy Adler. License: GPL version 2 or version 3
0023 % $Id: calc_slices.html 2819 2011-09-07 16:43:11Z aadler $
0024 
0025 if isstr(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0026 
0027 np = calc_colours('npoints');
0028 try   np = img(1).calc_colours.npoints;
0029 end
0030 
0031 if nargin < 2
0032     try 
0033         levels = img.calc_slices.levels; 
0034     catch
0035         levels = [];
0036     end
0037 end
0038 
0039 % Assume all fwd_models are same dimension (all 3D or 2D no mixed dims)
0040 if mdl_dim(img(1))==2 
0041    if nargin>1 && ~isempty(levels);
0042        if ~all(levels(1,:) == [inf,inf,0])
0043           warning('specified levels ignored for 2D FEM');
0044        end
0045    end
0046    levels= [Inf,Inf,0];
0047 elseif mdl_dim(img(1))==3 && isempty(levels)
0048    levels = [Inf Inf mean(img.fwd_model.nodes(:,3))];
0049    eidors_msg('calc_slices: no levels specified, assuming an xy plane',2);
0050 end
0051 
0052 
0053 rimg = [];
0054 for i=1:length(img)
0055    rimg = cat(3, rimg, calc_this_slice( img(i), levels, np) );
0056 end
0057 
0058 function rimg = calc_this_slice( img, levels, np)
0059     % If scalar levels then we just create that many cut planes on z-dimension
0060     fwd_model = img.fwd_model;
0061     if size(levels)== [1,1]
0062        zmax= max(fwd_model.nodes(:,3));
0063        zmin= min(fwd_model.nodes(:,3));
0064        levels = linspace(zmin,zmax, levels+2);
0065        levels = levels(2:end-1)'*[Inf,Inf,1];
0066     end
0067     num_levs= size(levels,1);
0068     if isfield(img,'elem_data')
0069        [elem_data, n_images] = get_img_data(img);
0070        rimg=zeros(np,np,n_images,num_levs);
0071 
0072        for lev_no = 1:num_levs
0073           level= levels( lev_no, 1:3 );
0074 
0075           rimg(:,:,:,lev_no) = calc_image_elems( elem_data, level, fwd_model, np);
0076        end
0077     elseif isfield(img,'node_data')
0078        node_data= [img.node_data];
0079        if size(node_data,1)==1; node_data=node_data';end
0080        n_images= size(node_data,2);
0081        rimg=zeros(np,np,n_images,num_levs);
0082 
0083        for lev_no = 1:num_levs
0084           level= levels( lev_no, 1:3 );
0085 
0086           rimg(:,:,:,lev_no) = calc_image_nodes( node_data, level, fwd_model, np);
0087        end
0088     else
0089        error('img does not have a data field');
0090     end
0091 
0092     % FILTER IMAGE
0093     try   filt = img.calc_slices.filter; 
0094     catch filt = 1; end
0095     try   scal = img.calc_slices.scale; 
0096     catch scal = 1; end
0097 
0098     
0099     if filt*scal ~= 1
0100        filt = scal * ( filt/sum(filt(:)) );
0101        rimg = filter_image(rimg, filt);
0102     end
0103 
0104 % Calculate an image by mapping it onto the node_ptr matrix
0105 % This makes a blocky image to nearest node -> no longer used
0106 function rimg= calc_image_nearestnodes( node_data, level, fwd_model, np);
0107    fwd_model.mdl_slice_mapper.npx  = np;
0108    fwd_model.mdl_slice_mapper.npy  = np;
0109    fwd_model.mdl_slice_mapper.level= level;
0110    node_ptr = mdl_slice_mapper( fwd_model, 'node' );
0111 
0112    backgnd= NaN;
0113    n_images= size(node_data,2);
0114    rval= [backgnd*ones(1,n_images); node_data];
0115    rimg= reshape( rval(node_ptr+1,:), np,np, n_images );
0116 
0117 % Calculate an image by interpolating it onto the elem_ptr matrix
0118 function rimg= calc_image_nodes( node_data, level, fwd_model, np)
0119 
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 
0124    nd_interp= mdl_slice_mapper( fwd_model, 'nodeinterp' );
0125    elem_ptr = mdl_slice_mapper( fwd_model, 'elem' );
0126    [sx,sy]= size(elem_ptr);
0127 
0128    node_ptr = fwd_model.elems; node_ptr = [0*node_ptr(1,:);node_ptr];
0129    node_ptr = reshape( node_ptr( elem_ptr+1, :), sx, sy, []);
0130 
0131    n_images = size(node_data,2);
0132    rimg= zeros(sx, sy, n_images);
0133    backgnd= NaN;
0134    
0135    for ni = 1:n_images
0136      znd = [backgnd;node_data(:,ni)]; % add NaN for background
0137      rimg(:,:,ni) = sum( znd(node_ptr+1) .* nd_interp, 3); 
0138    end
0139 
0140 
0141 % Calculate an image by mapping it onto the elem_ptr matrix
0142 function rimg= calc_image_elems( elem_data, level, fwd_model, np)
0143 
0144    fwd_model.mdl_slice_mapper.npx  = np;
0145    fwd_model.mdl_slice_mapper.npy  = np;
0146    fwd_model.mdl_slice_mapper.level= level;
0147    elem_ptr = mdl_slice_mapper( fwd_model, 'elem' );
0148 
0149    backgnd= NaN;
0150    n_images= size(elem_data,2);
0151    rval= backgnd*ones(size(elem_data)+[1,0]);
0152    rval(2:end,:) = elem_data;
0153    rimg= reshape( rval(elem_ptr+1,:), np,np, n_images );
0154 
0155 
0156 function  rimg = filter_image(rimg, filt);
0157    [sz1,sz2,sz3,sz4] = size(rimg);
0158    for j1 = 1:sz3; for j2 = 1:sz4; 
0159       rsl = rimg(:,:,j1,j2);
0160 
0161       rna = isnan(rsl);
0162       rsl(rna) = 0;
0163       rsl = conv2(rsl, filt, 'same');
0164       rsl(rna) = NaN;
0165 
0166       rimg(:,:,j1,j2) = rsl;
0167    end; end
0168 
0169 function do_unit_test
0170    img = mk_image( mk_common_model('a2c2',8));
0171    img.calc_colours.npoints = 8; 
0172 
0173    imc = calc_slices(img);
0174    imt = NaN*ones(8); imt(3:6,2:7) = 1; imt(2:7,3:6) = 1; 
0175    unit_test_cmp('cs 2d 1', imc, imt);
0176 
0177    imn = rmfield(img,'elem_data');
0178    imn.node_data = ones(size(imn.fwd_model.nodes,1),1);
0179    imc = calc_slices(imn);
0180    unit_test_cmp('cs 2d 2', imc, imt, 1e-14);
0181 
0182    img.elem_data(1:4) = 2;
0183    imc = calc_slices(img);
0184    imt(3:6,3:6) = 1; imt(4:5,4:5) = 2; 
0185    unit_test_cmp('cs 2d 3', imc, imt);
0186 
0187    imn.node_data(1:5) = 2;
0188    imc = calc_slices(imn);
0189    imt(3:6,3:6) = 1; imt(4:5,3:6) = 1.049020821501088;
0190    imt(3:6,4:5) = 1.049020821501088; imt(4:5,4:5) = 2;
0191    unit_test_cmp('cs 2d 4', imc, imt, 1e-14);
0192 
0193    imn.node_data(:) = 1; imn.node_data(1) = 4;
0194    imc = calc_slices(imn);
0195    imt(3:6,3:6) = 1; imt(4:5,4:5) = 1.575633893074693; 
0196    unit_test_cmp('cs 2d 5', imc, imt, 1e-14);
0197 
0198    imn.calc_colours.npoints = 7; 
0199    imc = calc_slices(imn);
0200    imt = NaN*ones(7); imt(2:6,2:6) = 1; imt(4,1:7)= 1; imt(1:7,4)= 1;imt(4,4) = 4; 
0201    unit_test_cmp('cs 2d 6', imc, imt, 1e-14);
0202 
0203 
0204 % 3D Tests
0205    img = calc_jacobian_bkgnd( mk_common_model('n3r2'));
0206    img.calc_colours.npoints = 8; 
0207    imn = calc_slices(img,[inf,inf,1]);
0208 
0209    imt = NaN*ones(8); imt(3:6,2:7) = 1; imt(2:7,3:6) = 1; 
0210    unit_test_cmp('cs 3d 1', imn, imt);
0211 
0212    imn = calc_slices(img,[inf,0,inf]);
0213    imt = NaN*ones(8); imt(1:8,3:6) = 1; 
0214    unit_test_cmp('cs 3d 2', imn, imt);
0215 
0216    % Should have no effect
0217    img.fwd_model.nodes(:,3) = img.fwd_model.nodes(:,3)-1;
0218    imn = calc_slices(img,[inf,0,inf]); 
0219    imt = NaN*ones(8); imt(1:8,3:6) = 1; 
0220    unit_test_cmp('cs 3d 3', imn, imt);
0221 
0222 
0223 % multi image struct
0224    img = mk_image( mk_common_model('a2c2',8));
0225    img.calc_colours.npoints = 8; 
0226    img(2) = img;
0227    imc = calc_slices(img); 
0228    imt = NaN*ones(8); imt(3:6,2:7) = 1; imt(2:7,3:6) = 1; 
0229    unit_test_cmp('cs mult 1', imc, cat(3,imt,imt));
0230 
0231    imgb = mk_image( mk_common_model('b2c2',8));
0232    imgb.calc_colours.npoints = 8; 
0233    img(3)=imgb;
0234    imc = calc_slices(img); 
0235    unit_test_cmp('cs mult 2', imc, cat(3,imt,imt,imt));
0236