sliceomatic_mia.m

医学图像处理matlab工具箱

  popset(gcf,'windowbuttonupfcn');
  popset(gcf,'windowbuttonmotionfcn');

  popset(d.tip,'visible');
  
  % Make sure whatever buttonupfcn on the figure is run now to "turn
  % off" whatever was going on before we got our callback on the
  % arrow.

  buf = get(gcf,'windowbuttonupfcn');
  if ~strcmp(buf,'')
    eval(buf);
  end

function movetipforarrow(tip, ax, value, position, va, ha)
% Setup the current data tip for a slice arrow, and show it's
% control value
  set(tip,'parent',ax, ...
          'string',sprintf('Value: %1.3f',value),...
          ... 'string','o', ...
          'units','data', ...
          'position', position, ...
          'verticalalignment', va,...
          'horizontalalignment', ha);
  set(tip,'units','pixels');
  % Put it onto d.axisiso so that
  % it always appears on top.
  %set(t,'parent',d.axiso);
  
function p=arrow(parent,dir,pos)

%   21012    21012      12345     12345
% 5  *-*   5   *     2   *     2   *  
% 4  | |   4  / \    1 *-*\    1  /*-*
% 3 ** **  3 ** **   0 |   *   0 *   |
% 2  \ /   2  | |   -1 *-*/   -1  \*-*
% 1   *    1  *-*   -2   *    -2   *  

  switch dir
   case 'down'
    pts=[ 0 1; -2 3; -1 3; -1 5; 1 5; 1 3; 2 3 ];
    mp = 'SOM leftright';
   case 'up'
    pts=[ 0 5; 2 3; 1 3; 1 1; -1 1; -1 3; -2 3; ];
    mp = 'SOM leftright';
   case 'right'
    pts=[ 5 0; 3 -2; 3 -1; 1 -1; 1 1; 3 1; 3 2 ];
    mp = 'SOM topbottom';
   case 'left'
    pts=[ 1 0; 3 2; 3 1; 5 1; 5 -1; 3 -1; 3 -2 ];
    mp = 'SOM topbottom';
  end

  f=[1 2 7; 3 4 5; 3 5 6 ];

  % Modify the arrows to look good no matter what
  % the data aspect ratio may be.
  if pos(1)
    lim=get(parent,'xlim');
    fivep=abs(lim(1)-lim(2))/15/5;
    pts(:,1)=pts(:,1)*fivep+pos(1);
  elseif pos(2)
    lim=get(parent,'ylim');
    fivep=abs(lim(1)-lim(2))/15/5;
    pts(:,2)=pts(:,2)*fivep+pos(2);
  end

  % Create the patches, and add app data to them to remember what
  % They are associated with.
  p(1)=patch('vertices',pts,'faces',1:size(pts,1),'facec','n','edgec','k',...
             'linewidth',2,'hittest','off',...
             'parent',parent);
  p(2)=patch('vertices',pts,'faces',f,'facec','g','facea',.5,'edgec','n',...
             'parent',parent,'tag','sliceomaticarrow');
  setappdata(p(2),'arrowcenter',pos);
  setappdata(p(2),'arrowedge',p(1));
  setappdata(p(2),'motionpointer',mp);

  
function p=localisocaps(isosurface,isocap)
% Isocap management
  
% Get relevant info from the isosurface.
  if nargin<2 | ~strcmp(get(isocap,'visible'),'off')
    data=getappdata(isosurface,'isosurfacedata');
    caps=isocaps(data,getappdata(isosurface,'isosurfacevalue'));
  end

  if nargin==2
    if ~strcmp(get(isocap,'visible'),'off')
      set(isocap,caps);
    end
    p=isocap;
  else
    p=patch(caps,'edgecolor','none','facecolor','flat',...
            'facelighting','none',...
            'tag','sliceomaticisocap');

    setappdata(p,'isosurface',isosurface);
    setappdata(isosurface,'isosurfacecap',p);
    
    d=getappdata(gcf,'sliceomatic');
    
    switch d.defcolor
     case 'faceted'
      set(p,'facec','flat','edgec','black');
     case 'flat'
      set(p,'facec','flat','edgec','none');
     case 'interp'
      set(p,'facec','interp','edgec','none');
     case 'texture'
      set(p,'facec','flat','edgec','none');
     case 'none'
      set(p,'facec','none','edgec','none');
    end
    switch d.defalpha
     case 'none'
      set(p,'facea',1);
     case 'flat'
      set(p,'facea','flat');
     case 'interp'
      set(p,'facea','interp');
     case 'texture'
      set(p,'facea','flat');
    end    
  end


function p=localisosurface(volume, data, datanormals, value, oldiso)
% Isosurface management
  
  pushset(gcf, 'pointer','watch');

  fv = isosurface(volume{:},data, value);
  clim=get(gca,'clim');
  cmap=get(gcf,'colormap');
  clen=clim(2)-clim(1);
  idx=floor((value-clim(1))*length(cmap)/clen);

  if nargin==5
    set(oldiso,fv,'facecolor',cmap(idx,:));
    p=oldiso;
    cap=getappdata(p,'isosurfacecap');
    if ~isempty(cap)
      localisocaps(p,cap);
    end
  else
    p=patch(fv,'edgecolor','none','facecolor',cmap(idx,:),...
            'tag', 'sliceomaticisosurface');
    d=getappdata(gcf,'sliceomatic');
    switch d.deflight
     case 'flat'
      set(p,'facelighting','flat');
     case 'smooth'
      set(p,'facelighting','phong');
    end
    setappdata(p,'isosurfacecap',[]);
  end

  setappdata(p,'isosurfacevalue',value);
  setappdata(p,'isosurfacedata',data);

  reducepatch(p,10000);
  isonormals(volume{:},datanormals,p);

  popset(gcf,'pointer');

function s=localslice(data, X, Y, Z, oldslice)
% Slice Management.  Uses specialized slicomatic slices, not slices
% created with the SLICE command.

  s=[];
  d=getappdata(gcf,'sliceomatic');

  ds=size(data);

  if ~isempty(X)
    xi=round(X);
    if xi > 0 & xi < ds(2)
      cdata=reshape(data(:,xi,:),ds(1),ds(3));
      [xdata ydata zdata]=meshgrid(X,1:ds(1),1:ds(3));
      st = 'X';
    else
      return
    end
  elseif ~isempty(Y)
    yi=round(Y);
    if yi > 0 & yi < ds(1)
      cdata=reshape(data(yi,:,:),ds(2),ds(3));
      [xdata ydata zdata]=meshgrid(1:ds(2),Y,1:ds(3));
      st = 'Y';
    else
      return    
    end
  elseif ~isempty(Z)
    zi=round(Z);
    if zi > 0 & zi < ds(3)
      cdata=reshape(data(:,:,zi),ds(1),ds(2));
      [xdata ydata zdata]=meshgrid(1:ds(2),1:ds(1),Z);
      st = 'Z';
    else
      return
    end
  else
    error('Nothing was passed into LOCALSLICE.');
  end

  cdata=squeeze(cdata);
  xdata=squeeze(xdata);
  ydata=squeeze(ydata);
  zdata=squeeze(zdata);

  if nargin == 5
    % Recycle the old slice
    set(oldslice,'cdata',cdata,'alphadata',cdata, 'xdata',xdata, ...
                 'ydata',ydata, 'zdata',zdata);
    s=oldslice;
    %delete(news);
    if propcheck(s,'facec','texturemap')
      textureizeslice(s,'on');
    end
    
  else
    % setup the alphadata
    news=surface('cdata',cdata,'alphadata',cdata, 'xdata',xdata, ...
                 'ydata',ydata, 'zdata',zdata);
    set(news,'alphadata',cdata,'alphadatamapping','scaled','tag','sliceomaticslice',...
             'facelighting','none',...
             'uicontextmenu',d.uic);
    s=news;
    switch d.defcolor
     case 'faceted'
      set(s,'facec','flat','edgec','k');
     case 'flat'
      set(s,'facec','flat','edgec','n');
     case 'interp'
      set(s,'facec','interp','edgec','n');
     case 'texture'
      set(s,'facec','texture','edgec','n');
    end
    switch d.defalpha
     case 'none'
      set(s,'facea',1);
     case 'flat'
      set(s,'facea','flat');
     case 'interp'
      set(s,'facea','interp');
     case 'texture'
      set(s,'facea','texture');
    end    

    setappdata(s,'slicetype',st);
    
    if strcmp(d.defcolor,'texture')
      textureizeslice(s,'on');
    end
  end

  contour = getappdata(s,'contour');
  if ~isempty(contour)
    localcontour(s, contour);
  end
  

function textureizeslice(slice,onoff)
% Convert a regular slice into a texture map slice, or a texture
% slice into a regular slice.
  
  for k=1:prod(size(slice))

    d=getappdata(slice(k),'textureoptimizeations');

    switch onoff
     case 'on'
      d.xdata=get(slice(k),'xdata');
      d.ydata=get(slice(k),'ydata');
      d.zdata=get(slice(k),'zdata');
      setappdata(slice(k),'textureoptimizeations',d);
      if max(size(d.xdata)==1)
        nx=[d.xdata(1) d.xdata(end)];
      else
        nx=[d.xdata(1,1)   d.xdata(1,end);
            d.xdata(end,1) d.xdata(end,end)];
      end
      if max(size(d.ydata)==1)
        ny=[d.ydata(1) d.ydata(end)];
      else
        ny=[d.ydata(1,1)   d.ydata(1,end);
            d.ydata(end,1) d.ydata(end,end)];
      end
      if max(size(d.zdata)==1)
        nz=[d.zdata(1) d.zdata(end)];
      else
        nz=[d.zdata(1,1)   d.zdata(1,end);
            d.zdata(end,1) d.zdata(end,end)];
      end
      set(slice(k),'xdata',nx, 'ydata', ny, 'zdata', nz,...
                   'facec','texturemap');
      if ischar(get(slice(k),'facea'))
        set(slice(k),'facea','texturemap');
      end
      if ischar(get(slice(k),'facec'))
        set(slice(k),'facec','texturemap');
      end
     case 'off'
      if ~isempty(d)
        set(slice(k),'xdata',d.xdata,'ydata',d.ydata,'zdata',d.zdata);
        setappdata(slice(k),'textureoptimizeations',[]);
      end
      if ischar(get(slice(k),'facea')) & strcmp(get(slice(k),'facea'),'texturemap')
        set(slice(k),'facea','flat');
      end
      if ischar(get(slice(k),'facec')) & strcmp(get(slice(k),'facec'),'texturemap')
        set(slice(k),'facec','flat');
      end
    end
  end

function localcontour(slice,oldcontour)
% Create a contour on SLICE
% When OLDCONTROUR, recycle that contour patch.
% This does not use the CONTOURSLICE command, but instead uses a
% specialized slice created for sliceomantic.

  d=getappdata(gcf,'sliceomatic');
  
  cdata = get(slice,'cdata');
  st = getappdata(slice,'slicetype');

  % Calculate the new contour for CDATA's values.
  c = contourc(cdata);

  newvertices = [];
  newfaces = {};
  longest = 1;
  cdata = [];
  
  limit = size(c,2);
  i = 1;
  h = [];
  color_h = [];
  while(i < limit)
    z_level = c(1,i);
    npoints = c(2,i);
    nexti = i+npoints+1;

    xdata = c(1,i+1:i+npoints);
    ydata = c(2,i+1:i+npoints);

    switch st
     case 'X'
      xv = get(slice,'xdata');
      lzdata = xv(1,1) + 0*xdata;
      vertices = [[lzdata].', [ydata].', [xdata].'];
     case 'Y'
      yv = get(slice,'ydata');
      lzdata = yv(1,1) + 0*xdata;
      vertices = [[ydata].', [lzdata].', [xdata].'];
     case 'Z'
      zv = get(slice,'zdata');
      lzdata = zv(1,1) + 0*xdata;
      vertices = [[xdata].', [ydata].', [lzdata].'];
    end    
    
    faces = 1:length(vertices);
    faces = faces + size(newvertices,1);

    longest=max(longest,size(faces,2));
    
    newvertices = [ newvertices ; vertices ];
    newfaces{end+1} = faces;
    
    tcdata =  (z_level + 0*xdata).';
    
    cdata = [ cdata; tcdata ]; % need to be same size as faces
    
    i = nexti;
  end

  % Fix up FACES, which is a cell array.
  faces = [];
  for i = 1:size(newfaces,2)
    faces = [ faces;
              newfaces{i} ones(1,longest-size(newfaces{i},2))*nan nan ];
    % Nans don't work in patches in OpenGL with trailing NaNs, but
    % if I fake it out, edges with breaks in them don't work.
    % Bummer!  I tried filling in with the last data point, but
    % then unbroken segmets would fill in again.    
  end
  
  if isempty(oldcontour)
    oldcontour = patch('facecolor','none', 'edgecolor',d.defcontourcolor,...
                       'linewidth',d.defcontourlinewidth);
    setappdata(slice,'contour',oldcontour);
  end

  set(oldcontour,'vertices',newvertices,...
                 'faces',faces,...
                 'facevertexcdata',cdata);
  
function ss=allSlices
  ss=findobj(gcf,'type','surface','tag','sliceomaticslice');

function ss=allIsos
  ss=findobj(gcf,'type','patch','tag','sliceomaticisosurface');

function ss=allCaps
  ss=findobj(gcf,'type','patch','tag','sliceomaticisocap');

function working(onoff)

  ax=getappdata(gcf,'workingaxis');

  if isempty(ax)
    ax=axes('units','norm','pos',[.3 .4 .4 .2],...
            'box','on','ytick',[],'xtick',[],...
            'xlim',[-1 1],'ylim',[-1 1],...
            'color','none','handlevis','off');
    text('parent',ax,'string','Working...','fontsize',64,...
         'pos',[0 0], ...
         'horizontalalignment','center',...
         'verticalalignment','middle',...
         'erasemode','xor');
    setappdata(gcf,'workingaxis',ax);
  end

  disp(['Working...' onoff]);
  set([ax get(ax,'children')],'vis',onoff);