Code covered by the BSD License  

Highlights from
Regular Control Point Interpolation Matrix with Boundary Conditions

image thumbnail

Regular Control Point Interpolation Matrix with Boundary Conditions

by

 

06 Jan 2010 (Updated )

Creates Toeplitz-like matrices representing interpolation operations with edge conditions.

T=interpMatrix(kernel, origin, numCtrlPoints, CtrlPointSep, extraprule)
function T=interpMatrix(kernel, origin, numCtrlPoints, CtrlPointSep, extraprule)
% INTERPMATRIX - This function creates a sparse Toeplitz-like matrix representing 
% a regularly-spaced interpolation operation between a set of control points. The 
% user can specify the interpolation kernel, the number of control points, the
% spacing between the control points, and certain boundary conditions governing  
% the behavior at the first and last control point. 
%
% The tool has obvious applications to interpolation, curve fitting, and signal 
% reconstruction. More generally, the ability to represent interpolation as a 
% matrix is useful for minimizing cost functions involving interpolation operations.
% For such functions, the interpolation matrix and its transpose inevitably arise
% in the gradient.
% 
% Although the matrix generated by interpMatrix() is for 1D interpolation, it can be
% generalized to n-dimensional tensorial interpolation using kron(). However, a more
% efficient alternative to kron() is this tool,
% 
% http://www.mathworks.com/matlabcentral/fileexchange/25969-efficient-object-oriented-kronecker-product-manipulation
% 
%USAGE:
%
%
%  T=interpMatrix(kernel, origin, numCtrlPoints, CtrlPointSep, extraprule)
%
%out:
%
%     T: sparse output matrix. The columns of T are copies of a common interpolation
%        kernel (with adjustments for boundary conditions), but  shifted in
%        increments of the CtrlPointSep parameter (see below) to different control
%        point locations. The result is that if x is a vector of coefficients,         
%        then T*x is the interpolation of these coefficients on the interval 
%        enclosed by the control points.
% 
%
%in:
%
%  kernel: vector containing samples of an interpolation function, shifted copies 
%          of which will be used to create the columns of T. This vector never
%          needs to be zero-padded. Zero-padding is derived automatically from  
%          the other input arguments below.
%
%  origin: Index i such that kernel(i) is located at the first control point.
%          It is also possible to specify the origin using the following
%          string options:
%
%            'max': origin i will be selected where kernel(i) is maximized.
%            'ctr': origin i will be selected as ceil((length(kernel)+1)/2)
%
%  numCtrlPoints: number of control points in the system.
%
%  CtrlPointSep: a stricly positive integer indicating the number of samples between
%                control points.
%
%  extraprule: String argument. Initially, the shifted copies of "kernel" form 
%            the columns of T. The columns are then modified to satisfy edge conditions
%            indicated by the "extraprules" parameter (options for this parameter
%            are discussed below).
%
%
%Choices of extraprule
%-------------------------
%
%The following 4 choices of extraprule result in a matrix T of size MxN where
%N=numCtrlPoints and M=(numCtrlPoints-1)*CtrlPointSep+1.
% 
%'zero' (default) - Extrapolation with zeroes. The columns of T are initially
%                   linear shifted copies of the kernel samples. Rows
%                   from the top and bottom of T are then discarded
%                   so that T(1,:)*x is the interpolated value at the 1st
%                   control point, for some coefficient vector x, while T(end,,:)*x 
%                   is the interpolated value at the last control point. 
%                   When CtrlPointSep=1, the resulting matrix T is a Toeplitz matrix.
%                  
%                   Equivalently, this means that T*x gives the same result as if there
%                   were infinite control points and if the coefficient vector x
%                   were extrapolated with zeros, i.e.,  x(i)=0 for i<=1 and for i>=N.
%
%
%'circ'           - Similar to 'zero' except the columns are related by a
%                   circulant shift of the kernel samples. When CtrlPointSep=1, 
%                   T is a circulant matrix.
%
%
%'rep'            - Extrapolation by replication. This option is the same as 'zero' 
%                   except that the first and final columns of T are modified.
%                   The modification is such that T*x gives the same result
%                   as if there were infinite control points and under the assumption                   
%                   that the coefficient vector x is extrapolated so 
%                   that x(i)=x(1) for i<=1 and x(i)=x(N) for i>=N.
%
%
%'mirror'         - Extrapolation by mirroring. Similar to 'rep' except the 
%                   extrapolation of x(i) uses mirroring at the array
%                   boundaries instead of replication. 
%
%
%The final option for the extraprule parameter is 'allcontrib', which is different from the
%other options in that the output matrix T has more columns:
%
%
%'allcontrib'       - This option adds control points to the system,
%                     with corresponding additional columns in T, 
%                     such that all control points that can affect the interval,                     
%                     Q, enclosed by the initial set of control points are now
%                     included. 
%
%                     As before, however, T(1,:)*x is the interpolated value at the 1st
%                     sample point in Q, while T(end,,:)*x is the interpolated value
%                     at the last sample point in Q. 
%
%
% 
%Copyright 2005, Matt Jacobson, The University of Michigan 


if nargin<5
 extraprule='zero';
else
 extraprule=strrep(lower(extraprule),' ',''); %shouldn't contain spaces
end

if ~any(ismember(extraprule,{'zero','rep','mirror','circ','allcontrib'}))
   error 'Unrecognized extraprule selection.' 
end

if nargin<4,
 CtrlPointSep=1;
end

if ischar(origin)
    
    switch origin
       
        case 'max'
        
            [~,origin]=max(kernel);
            
        case 'ctr'
            
            origin = ceil( (length(kernel)+1)/2 );
            
        otherwise
            error 'Unrecognized option for "origin"'
    end
    
end

%%


kernel=kernel(:);
causal=kernel(origin:end);
anticausal=kernel(1:origin-1);

lenkernel=length(kernel);
lencausal=length(causal);
lenanticausal=length(anticausal);



if ~isequal(extraprule,'circ')
    


  extraleft=floor((lencausal-1)/CtrlPointSep);
  extraright=floor(lenanticausal/CtrlPointSep);

  numCtrlPoints=numCtrlPoints+extraleft+extraright;

  T=cell(1,numCtrlPoints);

  len=CtrlPointSep*(numCtrlPoints-1)+lenkernel;
  prototype=col0s(len);
  prototype(1:lenkernel)=kernel;


  T=rawMatrix(prototype, CtrlPointSep,numCtrlPoints);
  
  
  
    aa=extraleft*CtrlPointSep+origin;
    bb=aa+(numCtrlPoints-extraleft-extraright-1)*CtrlPointSep;
  
   if isequal(extraprule,'rep')


     xx=col1s(numCtrlPoints); 
     xx(extraleft+2:end)=0;


     yy=col1s(numCtrlPoints); 
     yy(1:end-(extraright)-1)=0;


     xx=T*xx;     yy=T*yy;
     T=T(:,extraleft+1 : end-extraright);
    
     T(:,1)=xx; T(:,end)=yy;    


   elseif isequal(extraprule,'zero') %ZERO-PADDED (TEOPLITZ) END CONDTIONS

      T=T(:,extraleft+1 : end-extraright);


   elseif isequal(extraprule,'mirror') %MIRROR END CONDTIONS

       
     Tc=T(:,extraleft+1 : end-extraright);  
       
     mod1b=@(i,N) mod(i-1,N)+1; %1-based index modulus
     IndexMap=@(i,N) round( (1-abs(  (mod1b(i,2*N) -(N+.5)  )/(N+.5) ) )*(N+.5));
       
     colAxis=(1:numCtrlPoints)-extraleft;
     
     nn=size(Tc,2);
     mappedCols=IndexMap(colAxis,nn);
    
     V=bsxfun(@eq,  mappedCols.' , (1:nn) );
    
     T=T*V;


   end


     T=T(aa:bb , :);



else %CIRCULANT END CONDITIONS

  len=CtrlPointSep*(numCtrlPoints-1)+1;
  
  prototype=col0s(len); 
     p2=prototype;
     
  prototype(1:lencausal)=causal;
    p2(end+1-lenanticausal:end)=anticausal;

  prototype=prototype+p2;  
    


  T=rawMatrix(prototype, CtrlPointSep,numCtrlPoints);
  

  
end

T=sparse(T);



function c=col0s(nn)

    c=sparse(nn,1);
    %c=zeros(nn,1);


function c=col1s(nn)

    c=ones(nn,1);
    
    
function T=rawMatrix(prototype, CtrlPointSep,numCtrlPoints)

    N=length(prototype);
    mod1b=@(t) mod(t-1,N)+1;
  
    [ii,jj,ss]=find(prototype(:));
    
    nn=length(ii);
    
    II=repmat(ii,1,numCtrlPoints);
       II=mod1b(bsxfun(@plus,II,(0:numCtrlPoints-1)*CtrlPointSep));
    JJ=repmat(1:numCtrlPoints,nn,1);   
    SS=repmat(ss,1,numCtrlPoints);
    
 
    T=sparse(II(:),JJ(:),SS(:));

Contact us