show_3d_slices

 show_3d_slices(img, z_cuts, x_cuts, y_cuts)
 Show a 3d view of an object with many slices through it
  z_cuts = planes in z to do a cut
  x_cuts = planes in x to do a cut
  y_cuts = planes in y to do a cut
 Default show 2 z_cuts and 1 x and 1 y cut

0001 function show_3d_slices(img, varargin);
0002 % show_3d_slices(img, z_cuts, x_cuts, y_cuts)
0003 % Show a 3d view of an object with many slices through it
0004 %  z_cuts = planes in z to do a cut
0005 %  x_cuts = planes in x to do a cut
0006 %  y_cuts = planes in y to do a cut
0007 % Default show 2 z_cuts and 1 x and 1 y cut
0008 
0009 % (C) 2007 Andy Adler. License: GPL version 2 or version 3
0010 % $Id: show_3d_slices.html 2819 2011-09-07 16:43:11Z aadler $
0011 
0012 if isstr(img) && strcmp(img,'UNIT_TEST'); do_unit_test; return; end
0013 
0014 cla;
0015 hold on
0016 
0017 [xyz_max, xyz_min, rimg, cimg, ...
0018           x_cuts, y_cuts, z_cuts] = get_slices(img, varargin{:});
0019 
0020 for zi= 1:length(z_cuts)
0021    M_trans= [1,0;0,1;0,0];
0022    M_add= [0,0,z_cuts(zi)];
0023    idx= zi;
0024    surf_slice(rimg(:,:,idx), cimg(:,:,idx), xyz_min, xyz_max, ...
0025               M_trans, M_add, 1);
0026 end
0027 
0028 for xi= 1:length(x_cuts)
0029    M_trans= [0,0;1,0;0,1];
0030    M_add= [x_cuts(xi),0,0];
0031    idx= length(z_cuts) + xi;
0032    surf_slice(rimg(:,:,idx), cimg(:,:,idx), xyz_min, xyz_max, ...
0033               M_trans, M_add, 1);
0034 end
0035 
0036 for yi= 1:length(y_cuts)
0037    M_trans= [1,0;0,0;0,1];
0038    M_add= [0,y_cuts(yi),0];
0039    idx= length(z_cuts) + length(x_cuts) + yi;
0040    surf_slice(rimg(:,:,idx), cimg(:,:,idx), xyz_min, xyz_max, ...
0041               M_trans, M_add, 1);
0042 end
0043 
0044 hold off
0045 
0046 function [xyz_max, xyz_min, rimg, cimg, ...
0047           x_cuts, y_cuts, z_cuts] = get_slices(img, varargin);
0048    xyz_max= max(img.fwd_model.nodes);
0049    xyz_min= min(img.fwd_model.nodes);
0050    if nargin==1;
0051       % Default show 2 z_cuts and 1 x and 1 y cut
0052        x_cuts= linspace(xyz_min(1), xyz_max(1), 3); x_cuts([1,3])=[];
0053        y_cuts= linspace(xyz_min(2), xyz_max(2), 3); y_cuts([1,3])=[];
0054        z_cuts= linspace(xyz_min(3), xyz_max(3), 4); z_cuts([1,4])=[];
0055    elseif nargin==2;
0056        z_cuts= varargin{1};
0057        x_cuts= [];
0058        y_cuts= [];
0059    elseif nargin==3;
0060        z_cuts= varargin{1};
0061        x_cuts= varargin{2};
0062        y_cuts= [];
0063    elseif nargin==4;
0064        z_cuts= varargin{1};
0065        x_cuts= varargin{2};
0066        y_cuts= varargin{3};
0067    else 
0068        error('too many inputs');
0069    end
0070 
0071    limts= [ z_cuts(:)*[inf,inf,1  ]; ...
0072             x_cuts(:)*[1  ,inf,inf]; ...
0073             y_cuts(:)*[inf,1  ,inf]];
0074 
0075    rimg0= calc_slices( img, limts);
0076 % SURF DOESN'T SHOW THE BLOODY OUTER BOUNDARY, WE NEED TO ADD 4 POINTS
0077    rimg= NaN*ones(size(rimg0,1)+4,size(rimg0,2)+4,size(limts,1));
0078    rimg(3:end-2,3:end-2,:)= rimg0;
0079    cimg = calc_colours( rimg, img);
0080 
0081 function xyz= linspace_plus4(lim_min, lim_max, np);
0082    ooo= ones(length(lim_min),1);
0083    delta = lim_max(:)-lim_min(:);
0084    delta = max(delta)*ones(size(delta)); 
0085    middle=(lim_max(:)+lim_min(:))/2;
0086    xyz= middle(:)*ones(1,np+4) + ...
0087         delta/(np-1)*[-(np+3)/2:(np+3)/2]; % for plus4
0088 %       delta/(np-1)*[-(np+1)/2:(np+1)/2]; % for plus2
0089 %       delta/(np-1)*[-(np-1)/2:(np-1)/2]; % for plus0
0090 
0091 function surf_slice(rimg, cimg, xyz_min, xyz_max, M_trans, M_add, show_surf);
0092    np= size(rimg,1)-4;
0093    lim_min= xyz_min*M_trans;
0094    lim_max= xyz_max*M_trans;
0095    xyz= linspace_plus4( lim_min, lim_max, np);
0096    [x,y]= meshgrid(xyz(1,:),xyz(2,:));
0097    xyz= reshape( [x(:),y(:)]*M_trans', np+4,np+4,3);
0098 
0099    ff=isnan(rimg);
0100    bdr= (conv2(double(~ff),ones(3),'same')>0) & ff;
0101    outbdr = ff & ~bdr;
0102 
0103    ver = eidors_obj('interpreter_version');
0104    if ver.isoctave ==0 && ver.is64bit
0105       eidors_msg(['Poor you. You''re using a 64 bit version of matlab. ' ...
0106                   'Unfortunately, these versions have serious graphics ' ...
0107                   'bugs and we can''t show a nice image. Sorry. ' ...
0108                   'Please bug Mathworks for a fix (good luck!).'],0);
0109    else
0110       cimg(outbdr)= NaN;
0111    end
0112 
0113    if show_surf
0114       hh=surf(xyz(:,:,1)+M_add(1), ...
0115               xyz(:,:,2)+M_add(2), ...
0116               xyz(:,:,3)+M_add(3), flipud(cimg));
0117 
0118       set(hh,'EdgeAlpha',0); % Remove background grid
0119 
0120       % WHY WOULD CDataMapping BE ANYTHING ELSE  - STUPID MATLAB !!!
0121       % In 64 bit matlab, doing CDataMapping direct on a matrix with
0122       %   NaN's will crash matlab. Damn.
0123       set(hh,'CDataMapping','direct');
0124    end
0125 
0126 %  draw_line_around(cimg, rimg, x,y, M_trans, M_add);
0127 
0128 
0129 function draw_line_around(cimg, rimg, x,y, M_trans, M_add);
0130 % The MATLAB contour functions are inflexible crap. We
0131 % Need to completely break them to get it to work
0132 % [x_ax;y_ax;z_ax] = M_trans * [x_ax_matlab;y_ax_matlab] + M_add
0133 % For x,z plane at y level 4 we have
0134 % M_trans = [1 0;0 0;0 1]; M_add= [0;4;0];
0135 
0136    cimg(isnan(rimg))= -1e50;
0137    [jnk,hh]=contour(x,y,cimg, [-1e49,-1e49]);
0138    Contour_paths = M_trans*[ get(hh,'Xdata'), ...
0139                              get(hh,'Ydata') ]';
0140    set(hh,'EdgeColor',[0,0,0], ...
0141           'Xdata', Contour_paths(1,:) + M_add(1), ...
0142           'Ydata', Contour_paths(2,:) + M_add(2), ...
0143           'Zdata', Contour_paths(3,:) + M_add(3));
0144 
0145 function do_unit_test
0146    img = mk_image(mk_common_model('n3r2',16),1);
0147    show_3d_slices(img,[1,2],0,0.5);
0148    view(10,18);