/matlab/plf_preimages.m

  1% libsrvf 
  2% =======
  3%
  4% A shape analysis library using the square root velocity framework.
  5% 
  6% Copyright (C) 2012   FSU Statistical Shape Analysis and Modeling Group
  7% 
  8% This program is free software: you can redistribute it and/or modify
  9% it under the terms of the GNU General Public License as published by
 10% the Free Software Foundation, either version 3 of the License, or
 11% (at your option) any later version.
 12% 
 13% This program is distributed in the hope that it will be useful,
 14% but WITHOUT ANY WARRANTY; without even the implied warranty of
 15% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 16% GNU General Public License for more details.
 17% 
 18% You should have received a copy of the GNU General Public License
 19% along with this program.  If not, see <http://www.gnu.org/licenses/>
 20% --------------------------------------------------------------------
 21
 22
 23% Computes preimages of the values in xv under the non-decreasing 
 24% piecewise-linear function F,T.  The function must be a 1-D function, 
 25% and xv must be non-decreasing.
 26%
 27% If F is not a one-to-one function (i.e. F(i)==F(i+1) for some i), 
 28% then there may be more than one preimage for a given value xv(j).
 29% In this case, there is at the moment no convention determining 
 30% which preimage is returned for such a value.  The only guarantee 
 31% is that if you do this
 32%
 33% xvi=plf_preimages(F,T,xv);
 34% Fxvi=plf_evaluate(F,T,xvi);
 35%
 36% then Fxvi will be equal (or at least close) to xv.
 37% TODO: is there a sensible convention for choosing preimages?
 38%
 39% Inputs
 40%  F,T:  the function.  F and T must be non-decreasing.
 41%  xv:   the sample points.  Must be non-decreasing.
 42%
 43% Outputs
 44%  xvi:  preimages of xv
 45% --------------------------------------------------------------------------
 46function xvi = plf_preimages( F, T, xv )
 47  assert( size(F,1) == 1 );
 48  assert( min(diff(F)) >= 0 );
 49  assert( min(diff(T)) >= 0 );
 50
 51  xvi = zeros(1,length(xv));
 52  Fidx = 1;
 53  for xvidx = 1:length(xv)
 54    while( Fidx < length(F)-1 && xv(xvidx) > F(Fidx+1) )
 55      Fidx = Fidx+1;
 56    end
 57
 58    dF = F(Fidx+1) - F(Fidx);
 59    if ( dF > 1e-4 )
 60      w1 = (F(Fidx+1) - xv(xvidx)) / dF;
 61      w2 = (xv(xvidx) - F(Fidx)) / dF;
 62      xvi(xvidx) = w1 * T(Fidx) + w2 * T(Fidx+1);
 63    else
 64      xvi(xvidx) = T(Fidx);
 65    end
 66  end
 67end
 68
 69
 70%!function v=_random_increasing_vector(v0,n)
 71%!  v=zeros(1,n);
 72%!  v(1)=v0;
 73%!  for i=2:n
 74%!    v(i)=v(i-1)+rand();
 75%!  end
 76%! end
 77%!
 78%!test
 79%! F=linspace(0,1,5);
 80%! T=linspace(0,1,5);
 81%! xv=linspace(0,1,100);
 82%! xvi=plf_preimages(F,T,xv);
 83%! assert(xvi,xv,1e-5);
 84%!
 85%!test
 86%! F=[0 0 1/2 1/2 1 1];
 87%! T=linspace(0,1,6);
 88%! tv=linspace(0,1,100);
 89%! xv=[0 0.2 0.499 0.5 0.501 0.999 1];
 90%! xvie=[0 0.28 0.3996 0.4 0.6004 0.7996 0.8];
 91%! xvi=plf_preimages(F,T,xv);
 92%! assert(xvi,xvie,1e-4);
 93%!
 94%!test
 95%! F=_random_increasing_vector(-100*rand(),50);
 96%! T=_random_increasing_vector(-5*rand(),50);
 97%! tv=linspace(T(1),T(end),100);
 98%! xv=interp1(T,F,tv);
 99%! xvi=plf_preimages(F,T,xv);
100%! assert(xvi,tv,1e-4)
101%!
102%!test
103%! F=_random_increasing_vector(0,1000);
104%! F=F/F(end);
105%! T=linspace(0,1,length(F));
106%! tv=linspace(T(1),T(end),2000);
107%! xv=interp1(T,F,tv);
108%! xvi=plf_preimages(F,T,xv);
109%! Fxvi=plf_evaluate(F,T,xvi);
110%! assert(Fxvi,xv,1e-4)