sliceomatic.m

matlab源代码,适用于开发研究,带来很好的学习效果.

    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')
    d=getappdata(gcf,'sliceomatic');
    data=getappdata(isosurface,'isosurfacedata');
    if isnan(d.xmesh)==1
      caps=isocaps(data,getappdata(isosurface,'isosurfacevalue'));
    else
      caps=isocaps(d.xmesh,d.ymesh,d.zmesh,data,getappdata(isosurface,'isosurfacevalue'));
    end
  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');

  d=getappdata(gcf,'sliceomatic');
  
  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
    try
      set(oldiso,fv,'facecolor',cmap(idx,:));
    catch
      set(oldiso,fv,'facecolor','none');
    end
    p=oldiso;
    cap=getappdata(p,'isosurfacecap');
    if ~isempty(cap)
      localisocaps(p,cap);
    end
  else
    if isnan(d.xmesh)==1
      p=patch(fv,'edgecolor','none','facecolor',cmap(idx,:),...
              'tag', 'sliceomaticisosurface');
    else
      p=patch(fv,'edgecolor','none','facecolor',cmap(idx,:),...
              'tag', 'sliceomaticisosurface');
    end
    % 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 isnan(d.xmesh) == 1
      if xi > 0 && xi <= ds(2)
        cdata=reshape(data(:,xi,:),ds(1),ds(3));
        [xdata ydata zdata]=meshgrid(xi,1:ds(1),1:ds(3));
        st = 'X';
      else
        return
      end
    else
      if isequal(d.xdir,'reverse')==1
        locate_xi=histc(xi,flipdim(d.xmesh,2));
        slice_number=find(locate_xi);
        slice_number=length(d.xmesh)-slice_number+1;
      else 
        locate_xi=histc(xi,d.xmesh);
        slice_number=find(locate_xi);
      end
      if ~isempty(slice_number) && slice_number > 0 && slice_number <= ds(2)
        cdata=reshape(data(:,slice_number,:),ds(1),ds(3));
        [xdata ydata zdata]=meshgrid(X,d.ymesh,d.zmesh);
        st = 'X';
      else
        return
      end
    end
    
  elseif ~isempty(Y)
    yi=round(Y);
    if isnan(d.ymesh) == 1
      if yi > 0 && yi <= ds(1)
        cdata=reshape(data(yi,:,:),ds(2),ds(3));
        [xdata ydata zdata]=meshgrid(1:ds(2),yi,1:ds(3));
        st = 'Y';
      else
        return    
      end
    else
      if isequal(d.ydir,'reverse')==1
        locate_yi=histc(yi,flipdim(d.ymesh,2));
        slice_number=find(locate_yi);
        slice_number=length(d.ymesh)-slice_number+1;
      else 
        locate_yi=histc(yi,d.ymesh);
        slice_number=find(locate_yi);
      end
      if ~isempty(slice_number) && slice_number > 0 && slice_number <= ds(1)
        cdata=reshape(data(slice_number,:,:),ds(2),ds(3));
        [xdata ydata zdata]=meshgrid(d.xmesh,Y,d.zmesh);
        st = 'Y';
      else
        return
      end
    end
    
  elseif ~isempty(Z)
    zi=round(Z);
    if isnan(d.zmesh) == 1
      if zi > 0 && zi <= ds(3)
        cdata=reshape(data(:,:,zi),ds(1),ds(2));
        [xdata ydata zdata]=meshgrid(1:ds(2),1:ds(1),zi);
        st = 'Z';
      else
        return
      end
    else
      if isequal(d.zdir,'reverse')==1
        locate_zi=histc(zi,flipdim(d.zmesh,2));
        slice_number=find(locate_zi);
        slice_number=length(d.zmesh)-slice_number+1;
      else 
        locate_zi=histc(zi,d.zmesh);
        slice_number=find(locate_zi);
      end
      if ~isempty(slice_number) && slice_number > 0 && slice_number <= ds(3)
        cdata=reshape(data(:,:,slice_number),ds(1),ds(2));
        [xdata ydata zdata]=meshgrid(d.xmesh,d.ymesh,Z);
        st = 'Z';
      else
        return
      end
    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
    setappdata(s,'slicetype',st);
  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;
    setappdata(s,'slicetype',st);
    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    
  end
  
  contour = getappdata(s,'contour');
  if ~isempty(contour)
    try 
      levels = getappdata(s, 'contourlevels');
      if isempty(levels)~=1
        localcontour(s,contour,levels);
      else
        localcontour(s, contour);
      end
    catch
      localcontour(s, contour);
    end
  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,levels)
% 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.
  if nargin < 3
    if isnan(d.zmesh)==1
      c = contourc(cdata);
    else
      switch st
       case 'X'
        c = contours(d.zmesh,d.ymesh,cdata);
       case 'Y'
        c = contours(d.zmesh,d.xmesh,cdata);
       case'Z'
        c = contours(d.xmesh,d.ymesh,cdata);
      end
    end
  else
    if isnan(d.zmesh)==1
      c = contourc(cdata,levels);
    else
      switch st
       case 'X'
        c = contours(d.zmesh,d.ymesh,cdata,levels);
       case 'Y'
        c = contours(d.zmesh,d.xmesh,cdata,levels);
       case 'Z'
        c = contours(d.xmesh,d.ymesh,cdata,levels);
      end
    end
  end

  newvertices = [];
  newfaces = {};
  longest = 1;
  cdata = [];
  
  limit = size(c,2);
  i = 1;
  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

  % Glom a NAN on the end for loop-breaking
  newvertices = [ newvertices ; nan nan nan ];
  cdata = [ cdata ; nan ];
  
  vertmax = size(newvertices,1);
  
  % 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))*vertmax vertmax ];
  end
  
  if isempty(oldcontour)
    oldcontour = patch('facecolor','none', 'edgecolor',d.defcontourcolor,...
                       'linewidth',d.defcontourlinewidth);
    try
      set(oldcontour,'linesmoothing',d.defcontoursmooth);
    catch
    end
    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);