Home > atmlab > math > pointinterp.m

pointinterp

PURPOSE ^

INTERPD Point interpolation of 1D-5D data

SYNOPSIS ^

function B = pointinterp(agrids,A,pos,varargin)

DESCRIPTION ^

 INTERPD  Point interpolation of 1D-5D data

    Atmlabs function for interpolation of data to a set of points. This is
    basically an interface to to the standard interp1 and interpn functions.

    For "grid interpolation", see accompanying function *gridinterp*.

    The grids of the data to be interpolated (*A*) are packed into the vector
    array *agrids*. The interpolation points are given as a matrix where each
    column corresponds to a dimension of *A*. Accordingly, each row in *pos*
    gives the coordinate for an interpolation point.

    The Atmlab convention of that 1D objects (including grids) are column
    vectors is here strictly followed.

    The *extrap* optional argument works here as for *gridinterp*.

 FORMAT   B = pointinterp(agrids,A,pos,iopt)
        
 OUT   B        Result of interpolation. A column vector.
 IN    agrids   Grids of A, as an array of vectors.
       A        Data to be interpolated. 
       pos      Matrix where each row is a position.
 OPT   iopt     Interpolation option. See *interpn*. Default is 'linear'.
       extrap   Special treatment of extrapolation. See above. 
                Default is false.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

pointinterp.m

SOURCE CODE ^

0001 % INTERPD  Point interpolation of 1D-5D data
0002 %
0003 %    Atmlabs function for interpolation of data to a set of points. This is
0004 %    basically an interface to to the standard interp1 and interpn functions.
0005 %
0006 %    For "grid interpolation", see accompanying function *gridinterp*.
0007 %
0008 %    The grids of the data to be interpolated (*A*) are packed into the vector
0009 %    array *agrids*. The interpolation points are given as a matrix where each
0010 %    column corresponds to a dimension of *A*. Accordingly, each row in *pos*
0011 %    gives the coordinate for an interpolation point.
0012 %
0013 %    The Atmlab convention of that 1D objects (including grids) are column
0014 %    vectors is here strictly followed.
0015 %
0016 %    The *extrap* optional argument works here as for *gridinterp*.
0017 %
0018 % FORMAT   B = pointinterp(agrids,A,pos,iopt)
0019 %
0020 % OUT   B        Result of interpolation. A column vector.
0021 % IN    agrids   Grids of A, as an array of vectors.
0022 %       A        Data to be interpolated.
0023 %       pos      Matrix where each row is a position.
0024 % OPT   iopt     Interpolation option. See *interpn*. Default is 'linear'.
0025 %       extrap   Special treatment of extrapolation. See above.
0026 %                Default is false.
0027 
0028 % 2007-05-22   Created by Patrick Eriksson.
0029 
0030 
0031 function B = pointinterp(agrids,A,pos,varargin)
0032 %
0033 [iopt,extrap] = optargs( varargin, { 'linear', false } );
0034                                                                           %&%
0035 %- Check input                                                            %&%
0036 %                                                                         %&%
0037 rqre_nargin( 3, nargin )                                                  %&%
0038 rqre_datatype( agrids, @iscell );                                         %&%
0039 rqre_datatype( A, @isnumeric );                                           %&%
0040 rqre_datatype( pos, @istensor2 );                                         %&%
0041 indim  = length(agrids);
0042 outdim = size(pos,2);
0043 if dimens(A) > length(agrids)                                             %&%
0044   error( 'Dimensionality of A is higher than number of given grids.' );   %&%
0045 end                                                                       %&%
0046 if outdim < indim                                                         %&%
0047   error('There can not be fewer columns in *pos* than grids in *agrids*.');%&%
0048 end                                                                       %&%
0049 for d = 1 : indim                                                         %&%
0050   rqre_datatype( agrids{d}, @istensor1, 'Grids in *agrids*' );            %&%
0051   rqre_gridmatch( A, d, agrids{d} );                                      %&%
0052 end                                                                       %&%
0053 rqre_datatype( extrap, @isboolean );                                      %&%
0054 
0055 
0056 %- Check a or determine a's effective dimensions
0057 %
0058 if extrap    % Ignore singleton dimensions
0059   dims = find( size( A ) > 1 );
0060   A    = getdims( A, dims );
0061 else
0062   if outdim ~= indim                                                      %&%
0063     error( ['With *extrap* false, the number of columns in *pos* ',...    %&%
0064                 'and the number of grids in *agrids* must be equal.'] );  %&%
0065   end                                                                     %&%
0066   asize = size( A );                                                      %&%
0067   if length(asize) < outdim  ||  any( asize(1:outdim)==1 )                %&%
0068     error( ['With *extrap* false there can not be any singleton ',...     %&%
0069                                               'dimmensions in *A*.'] );   %&%
0070   end                                                                     %&%
0071   dims = 1:outdim;
0072 end
0073 
0074 dim = length( dims );
0075 
0076 if dim == 0
0077   %
0078   B = repmat( A, size(pos,1), 1 );
0079 
0080 elseif dim == 1
0081   %
0082   if extrap
0083     xi = handle_expand( agrids{dims(1)}, pos(:,dims(1)) );
0084   else
0085     xi = pos(:,dims(1));
0086   end
0087   %
0088   B = interp1( agrids{dims(1)}, A, xi, iopt );
0089 
0090 elseif dim == 2
0091   %
0092   if extrap
0093     xi = handle_expand( agrids{dims(1)}, pos(:,dims(1)) );
0094     yi = handle_expand( agrids{dims(2)}, pos(:,dims(2)) );
0095   else
0096     xi = pos(:,dims(1));
0097     yi = pos(:,dims(2));
0098   end
0099   %
0100   B = interpn( agrids{dims(1)}, agrids{dims(2)}, A, xi, yi, iopt ); 
0101 
0102 elseif dim == 3
0103   %
0104   if extrap
0105     xi = handle_expand( agrids{dims(1)}, pos(:,dims(1)) );
0106     yi = handle_expand( agrids{dims(2)}, pos(:,dims(2)) );
0107     zi = handle_expand( agrids{dims(3)}, pos(:,dims(3)) );
0108   else
0109     xi = pos(:,dims(1));
0110     yi = pos(:,dims(2));
0111     zi = pos(:,dims(3));
0112   end
0113   %
0114   B = interpn( agrids{dims(1)}, agrids{dims(2)}, agrids{dims(3)}, ...
0115                                                       A, xi, yi, zi, iopt ); 
0116 
0117 elseif dim == 4
0118   %
0119   if extrap
0120     xi = handle_expand( agrids{dims(1)}, pos(:,dims(1)) );
0121     yi = handle_expand( agrids{dims(2)}, pos(:,dims(2)) );
0122     zi = handle_expand( agrids{dims(3)}, pos(:,dims(3)) );
0123     ui = handle_expand( agrids{dims(4)}, pos(:,dims(4)) );
0124   else
0125     xi = pos(:,dims(1));
0126     yi = pos(:,dims(2));
0127     zi = pos(:,dims(3));
0128     ui = pos(:,dims(4));
0129   end
0130   %
0131   [xi,yi,zi,ui] = ndgrid( xi, yi, zi, ui );
0132   %
0133   B = interpn( agrids{dims(1)}, agrids{dims(2)}, ...
0134                agrids{dims(3)}, agrids{dims(4)}, A, xi, yi, zi, ui, iopt ); 
0135 
0136 elseif dim == 5
0137   %
0138   if extrap
0139     xi = handle_expand( agrids{dims(1)}, pos(:,dims(1)) );
0140     yi = handle_expand( agrids{dims(2)}, pos(:,dims(2)) );
0141     zi = handle_expand( agrids{dims(3)}, pos(:,dims(3)) );
0142     ui = handle_expand( agrids{dims(4)}, pos(:,dims(4)) );
0143     vi = handle_expand( agrids{dims(5)}, pos(:,dims(5)) );
0144   else
0145     xi = pos(:,dims(1));
0146     yi = pos(:,dims(2));
0147     zi = pos(:,dims(3));
0148     ui = pos(:,dims(4));
0149     vi = pos(:,dims(5));
0150   end
0151   %
0152   B = interpn( agrids{dims(1)}, agrids{dims(2)}, agrids{dims(3)}, ...
0153                agrids{dims(4)}, agrids{dims(5)}, A, xi, yi, zi, ui, vi, iopt ); 
0154 
0155 else                                                             %&%
0156   error( 'Interpolation above 5D is not handled.' );             %&%
0157 end
0158 
0159 return
0160 
0161 
0162 %---
0163 function xi = handle_expand(x,xi)
0164   %
0165   v1 = min( x );
0166   i1 = find( xi<v1 );
0167   %
0168   v2 = max( x );
0169   i2 = find( xi>v2 );
0170   %
0171   if ~isempty(v1) 
0172     xi(i1) = v1;
0173   end
0174   %
0175   if ~isempty(v2) 
0176     xi(i2) = v2;
0177   end
0178   %
0179 return

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005