streamsliceのmethod指定について

streamsliceのドキュメンテーションにしたがって内挿法としてnearestを指定したのですが、グラフを見る限り流線は明らかに速度場0の領域を通ってしまっています。また、’linear’に変えたときも同じ結果が出てくるのでうまく内挿法が指定できていないのではないかと考えています。

どこか間違っているのでしょうか。

streamslice(X,Y,U1,V1,4,'nearest');

[Atsushi Ueno]

>どこか間違っているのでしょうか。

streamslice(X,Y,U1,V1,4,'nearest');

いいえ、文法はどこも間違っていません。引数method='nearest'は意図通り渡されています。(下記parseargs()参照)

しかし、 WとZの無い(2次元)引数構成のstreamslice関数では引数methodが使われていません。(下記コード参照)

• 3次元：ベクトルU,V,WをX,Y,Zに合わせinterp3関数(引数method)で補間してスライス面を生成
• 2次元：ベクトルU,Vをそのまま使用

type streamslice

function [ret1, ret2] = streamslice(varargin)
    %STREAMSLICE  Streamlines in slice planes.
    %   STREAMSLICE(X,Y,Z,U,V,W,Sx,Sy,Sz) draws well spaced streamlines
    %   (with direction arrows) from vector data U,V,W in axis aligned
    %   x,y,z planes at the points in the vectors Sx,Sy,Sz. The arrays
    %   X,Y,Z define the coordinates for U,V,W and must be monotonic
    %   and 3D plaid (as if produced by MESHGRID).  V must be an
    %   M-by-N-by-P volume array. It should not be assumed that the
    %   flow is parallel to the slice plane; for example, in a
    %   streamslice at a constant z, the z component of the vector
    %   field, W, is ignored when calculating the streamlines for that
    %   plane. Streamslices are useful for determining where to start
    %   streamlines, streamtubes, and streamribbons.
    %
    %   STREAMSLICE(U,V,W,Sx,Sy,Sz) assumes
    %       [X Y Z] = meshgrid(1:N, 1:M, 1:P) where [M,N,P]=SIZE(V).
    %
    %   STREAMSLICE(X,Y,U,V)  draws well spaced streamlines (with
    %   direction arrows) from vector data U,V. The arrays X,Y define
    %   the coordinates for U,V and must be monotonic and 2-D plaid (as
    %   if produced by MESHGRID).
    %
    %   STREAMSLICE(U,V) assumes
    %       [X Y] = meshgrid(1:N, 1:M) where [M,N]=SIZE(V).
    %
    %   STREAMSLICE(..., DENSITY) modifies the automatic spacing of the
    %   streamlines. DENSITY must be greater than 0. The default value
    %   is 1; higher values will produce more streamlines on each
    %   plane. For example, 2 will produce approximately twice as many
    %   streamlines while 0.5 will produce approximately half as many.
    %
    %   STREAMSLICE(...,'arrows'), (the default) draws direction
    %   arrows. STREAMSLICE(...,'noarrows') suppresses drawing the
    %   direction arrows.
    %
    %   STREAMSLICE(...,'method') specifies the interpolation method to
    %   use. 'method' can be 'linear', 'cubic', or 'nearest'.  'linear'
    %   is the default (see INTERP3).
    %
    %   STREAMSLICE(AX,...) plots into AX instead of GCA.
    %
    %   H = STREAMSLICE(...) returns a vector of handles to LINE
    %   objects.
    %
    %   [VERTICES ARROWVERTICES] = STREAMSLICE(...) returns 2 cell
    %   arrays of vertices for drawing the streamlines and the
    %   arrows.  These can be passed to any streamline drawing function
    %   (streamline) to draw the plot.
    %
    %   Example 1:
    %      load wind
    %      streamslice(x,y,z,u,v,w,[],[],[5]);
    %
    %   Example 2:
    %      load wind
    %      [verts averts] = streamslice(u,v,w,10,10,10);
    %      streamline([verts averts]);
    %      spd = sqrt(u.*u + v.*v + w.*w);
    %      hold on; slice(spd,10,10,10);
    %      colormap(hot)
    %      shading interp
    %      view(30,50); axis(volumebounds(spd));
    %      camlight; material([.5 1 0])
    %
    %   Example 3:
    %      z = peaks;
    %      surf(z); hold on
    %      shading interp;
    %      [c ch] = contour3(z,20); 
	%      ch.EdgeColor = 'b';
    %      [u v] = gradient(z);
    %      h = streamslice(-u,-v);  % downhill
    %      set(h,'Color','k');
    %      for i=1:length(h);
    %        zi = interp2(z,h(i).XData, h(i).YData);
    %        h(i).ZData = zi;
    %      end
    %      view(30,50); axis tight
    %
    %   See also STREAMLINE, SLICE, CONTOURSLICE.
    
    %   Copyright 1984-2017 The MathWorks, Inc.
    
    % not implemented:
    %   The colors of the streamlines
    %   indicate if the flow is in or out of the slice plane. 

    if nargin > 0
        [varargin{:}] = convertStringsToChars(varargin{:});
    end
    
    [cax,args,nargs] = axescheck(varargin{:});
    [x, y, z, u, v, w, sx, sy, sz, density, arrows, method] = parseargs(nargs,args);
    
    % Take this out when other data types are handled
    u = double(u);
    v = double(v);
    w = double(w);
    
    
    if isempty(method)
        method = 'linear';
    end
    
    if isempty(density)
        density = 1;
    else
        density = sqrt(density);
    end
    
    if isempty(arrows)
        arrows = 1;
    end
    
    
    vertsout = {};
    avertsout = {};
    
    if isempty(cax)
        f = get(0, 'currentfigure');
        if ~isempty(f)
            cax = get(f, 'currentaxes');
        end
    end
    
    if isempty(cax) && nargout < 2
        cax = newplot;
    end
        
    if isempty(w)  % 2D
        
        [msg, x, y] = xyuvcheck(x,y,u,v);  error(msg) 
        [ny,nx] = size(u);
        if isempty(x)
            [x, y] = meshgrid(1:nx, 1:ny);
        end
        
        [vertsout, avertsout] = nicestreams(x,y,u,v,density,arrows);
        
    else           % 3D
        % For each plane in each axis, generate stream slices.
        [msg, x, y, z] = xyzuvwcheck(x,y,z,u,v,w);  error(msg) 
        
        [ny,nx,nz] = size(u);
        if isempty(x)
            [x, y, z] = meshgrid(1:nx, 1:ny, 1:nz);
        end
        
        % X axis
        [xi,yi,zi] = meshgrid(sx,y(:,1,1),z(1,1,:));
        vi = interp3(x,y,z,v,xi,yi,zi,method);
        wi = interp3(x,y,z,w,xi,yi,zi,method);
        
        for j = 1:length(sx)
            if ~( anynan(wi(:,j,:)) || anynan(vi(:,j,:)) )
                [verts, averts] = nicestreams(...
                    reshape(zi(:,j,:),[ny nz]),...
                    reshape(yi(:,j,:),[ny nz]),...
                    reshape(wi(:,j,:),[ny nz]),...
                    reshape(vi(:,j,:),[ny nz]),density,arrows);
                for k = 1:length(verts)
                    vv = verts{k};
                    if ~isempty(vv)
                        verts{k} = [vv(:,1)*0+sx(j) vv(:,2) vv(:,1)];
                    end
                end
                for k = 1:length(averts)
                    vv = averts{k};
                    if ~isempty(vv)
                        averts{k} = [vv(:,1)*0+sx(j) vv(:,2) vv(:,1)];
                    end
                end
                vertsout = [vertsout verts];
                avertsout = [avertsout averts];
            end
        end
        
        
        % Y axis
        [xi,yi,zi] = meshgrid(x(1,:,1),sy,z(1,1,:));
        ui = interp3(x,y,z,u,xi,yi,zi,method);
        wi = interp3(x,y,z,w,xi,yi,zi,method);
        
        for j = 1:length(sy)
            if ~( anynan(wi(j,:,:)) || anynan(ui(j,:,:)) )
                [verts, averts] = nicestreams(...
                    reshape(zi(j,:,:),[nx nz]),...
                    reshape(xi(j,:,:),[nx nz]),...
                    reshape(wi(j,:,:),[nx nz]),...
                    reshape(ui(j,:,:),[nx nz]),density,arrows);
                for k = 1:length(verts)
                    vv = verts{k};
                    if ~isempty(vv)
                        verts{k} = [vv(:,2) vv(:,1)*0+sy(j) vv(:,1)];
                    end
                end
                for k = 1:length(averts)
                    vv = averts{k};
                    if ~isempty(vv)
                        averts{k} = [vv(:,2) vv(:,1)*0+sy(j) vv(:,1)];
                    end
                end
                vertsout = [vertsout verts];
                avertsout = [avertsout averts];
            end
        end
        
        
        % Z axis
        [xi,yi,zi] = meshgrid(x(1,:,1),y(:,1,1),sz);
        ui = interp3(x,y,z,u,xi,yi,zi,method);
        vi = interp3(x,y,z,v,xi,yi,zi,method);
        
        for j = 1:length(sz)
            if ~( anynan(ui(:,:,j)) || anynan(vi(:,:,j)) )
                [verts, averts] = nicestreams(xi(:,:,j),yi(:,:,j),ui(:,:,j),vi(:,:,j),density,arrows);
                for k = 1:length(verts)
                    vv = verts{k};
                    if ~isempty(vv)
                        verts{k} = [vv(:,1)  vv(:,2) vv(:,1)*0+sz(j)];
                    end
                end
                for k = 1:length(averts)
                    vv = averts{k};
                    if ~isempty(vv)
                        averts{k} = [vv(:,1)  vv(:,2) vv(:,1)*0+sz(j)];
                    end
                end
                vertsout = [vertsout verts];
                avertsout = [avertsout averts];
            end
        end
        
    end
    
    if nargout >= 2
        ret1 = vertsout;
        ret2 = avertsout;
    else
        h = streamline(cax,[vertsout avertsout]);
        if nargout==1
            ret1 = h;
        end
    end
end

function [vertsout, avertsout] = nicestreams(x,y,u,v,density,arrows)
    
    stepsize = min(.1, (min(size(v))-1)/100);
    streamoptions = [stepsize min(10000,sum(size(v))*2/stepsize)];
    
    vertsout = {};
    avertsout = {};
    
    num = 20;
    
    nrstart = ceil(num*density);
    ncstart = ceil(num*density);
    nrend   = ceil(num*density*4);
    ncend   = ceil(num*density*4);
    nrarrow = ceil(num*density);
    ncarrow = ceil(num*density);
    
    xmin = min(x(:));
    xmax = max(x(:));
    ymin = min(y(:));
    ymax = max(y(:));
    xrange = xmax-xmin;
    yrange = ymax-ymin;
    
    incstartx = xrange/ncstart;
    inrstarty = yrange/nrstart;
    ixrangecs = ncstart/xrange*(1-eps);
    iyrangers = nrstart/yrange*(1-eps);
    ixrangece = ncend/xrange*(1-eps);
    iyrangere = nrend/yrange*(1-eps);
    
    % startgrid and endgrid are used to keep track of the location/density
    % of streamlines. As new streamlines are created, the values in these
    % matrices will indicate whether a streamline has passed through each
    % quadrant of the data space. startgrid is coarser grid, while endgrid
    % is a finer grid. startgrid is used to decide whether to start a new
    % streamline. If an existing streamline has already passed through a
    % quadrant, we won't start a new streamline. endgrid is used to limit
    % the density of the final streamlines. New streamlines will stop when
    % the reach a quandrant that is already occupied by an existing
    % streamline.
    startgrid = zeros(nrstart,ncstart);
    endgrid = zeros(nrend,ncend);
    
    if arrows
        arrowgrid = ones(nrarrow,ncarrow);
        arrowgrid(2:3:end,2:3:end) = 0;
        arrowscale = .01*(xrange+yrange)/density;
        ixrangeca = ncarrow/xrange*(1-eps);
        iyrangera = nrarrow/yrange*(1-eps);
    end
    
    [r, c] = meshgrid(1:nrstart, 1:ncstart);
    rc = [r(:) c(:)];
    %rc = rc(randperm(size(rc,1)),:);
    
    for k = 1:size(rc,1)
        %if mod(k,100)==0, disp([num2str(k) '/' num2str(size(rc,1))]), end
        r = rc(k,1);
        c = rc(k,2);
        if startgrid(r,c)==0
            
            startgrid(r,c)=1;
            
            xstart = xmin+(c-.5)*incstartx;
            ystart = ymin+(r-.5)*inrstarty;
            vertsf = stream2(x,y, u, v,xstart,ystart,streamoptions);
            vertsb = stream2(x,y,-u,-v,xstart,ystart,streamoptions);
            verts = [vertsf vertsb];
            vo = cell(1,2);
            for q = 1:2
                vv = verts{q};
                
                nanpos = find(isnan(vv));
                if ~isempty(nanpos)
                    if nanpos(1)==1
                        vv = [];
                    else
                        vv = vv(1:nanpos(1)-1,:);
                    end
                end
                
                if ~isempty(vv)
                    tcc = floor( (vv(1,1)-xmin)*ixrangece )+1;
                    trr = floor( (vv(1,2)-ymin)*iyrangere )+1;
                end
                
                for j=1:size(vv,1)
                    xc = vv(j,1);
                    yc = vv(j,2);
                    
                    % Calculate indices into startgrid (rr,cc), based on
                    % the coordinates of this particular data point on this
                    % streamline. As a streamline passes through
                    % coordinates, mark them off so that we do not start
                    % new streamlines in those coordinates.
                    cc = floor( (xc-xmin)*ixrangecs )+1;
                    rr = floor( (yc-ymin)*iyrangers )+1;
                    if cc > 0 && cc <= ncstart && rr > 0 && rr <= nrstart
                        startgrid(rr,cc)=1;
                    end
                    
                    % Now calculate rr and cc using a finer mesh so that
                    % they are indices into endgrid. As a streamline passes
                    % through coordinates, mark them off so that we do not
                    % start new streamlines in those coordinates. If a new
                    % streamline hits an existing streamline, then the new
                    % streamline will be truncated.
                    cc = floor( (xc-xmin)*ixrangece )+1;
                    rr = floor( (yc-ymin)*iyrangere )+1;
                    if cc <= 0 || cc > ncend || rr <= 0 || rr > nrend
                        break;
                    elseif endgrid(rr,cc)==1
                        if ~(any(cc==tcc) && any(rr==trr))
                            break;
                        end
                    else
                        tcc=cc; trr=rr;
                        endgrid(rr,cc)=1;
                    end
                    if arrows
                        cc = floor( (xc-xmin)*ixrangeca )+1;
                        rr = floor( (yc-ymin)*iyrangera )+1;
                        if j>1 && cc > 0 && cc <= ncarrow && rr > 0 && rr <= nrarrow && arrowgrid(rr,cc)==0
                            arrowgrid(rr,cc)=1;
                            d = vv(j,:)-vv(j-1,:);
                            d = d/norm(d);
                            if q==2, d = -d; end
                            
                            arrowlen = 1;
                            arrowwidth = .6;
                            p2 = [xc yc] - arrowlen*d .* arrowscale;
                            
                            av(1,:) = p2 - arrowwidth*[-d(2) d(1)].*arrowscale;
                            av(2,:) = [xc yc];
                            av(3,:) = p2 + arrowwidth*[-d(2) d(1)].*arrowscale;
                            
                            avertsout{end+1} = av;
                        end
                    end
                end
                vo{q} = vv(1:j-1,:);
                
            end
            vertsout{end+1} = [flipud(vo{2}); vo{1}(2:end,:)];
        end
    end
end

function [x, y, z, u, v, w, sx, sy, sz, density, arrows, method] = parseargs(nin, vargin)
    
    [x, y, z, w, sx, sy, sz, density, arrows, method] = deal([]);
    
    for j=1:2
        if nin>0
            lastarg = vargin{nin};
            if ischar(lastarg) % streamslice(...,'method'),  streamslice(...,'noarrows')
                if ~isempty(lastarg)
                    lastarg = lower(lastarg);
                    if lastarg(1)=='n' && length(lastarg)>=2 && lastarg(2)=='o'
                        arrows = 0;
                    elseif lastarg(1)=='a' && length(lastarg)>=2 && lastarg(2)=='r'
                        arrows = 1;
                    else
                        method = lastarg;
                    end
                end
                nin = nin - 1;
            end
        end
    end
    
    if nin==2 || nin==3       % streamslice(u,v) or streamslice(u,v,density)
        u = vargin{1};
        v = vargin{2};
        if nin==3, density = vargin{3}; end
    elseif nin==4 || nin==5   % streamslice(x,y,u,v) or streamslice(x,y,u,v,density)
        [x, y, u, v] = deal(vargin{1:4});
        if nin==5, density = vargin{5}; end
    elseif nin==6 || nin==7   % streamslice(u,v,w,sx,sy,sz) or streamslice(u,v,w,sx,sy,sz,density)
        [u, v, w, sx, sy, sz] = deal(vargin{1:6});
        if nin==7, density = vargin{7}; end
    elseif nin==9 || nin==10   % streamslice(x,y,z,u,v,w,sx,sy,sz) or streamslice(x,y,z,u,v,w,sx,sy,sz,density)
        [x, y, z, u, v, w, sx, sy, sz] = deal(vargin{1:9});
        if nin==10, density = vargin{10}; end
    else
        error(message('MATLAB:streamslice:WrongNumberOfInputs'));
    end
    
    sx = sx(:);
    sy = sy(:);
    sz = sz(:);
end

function out = anynan(x)
    out = any(isnan(x(:)));
end

[saito ishiue]

ありがとうございます。納得できました。