Path: news.mathworks.com!not-for-mail
From: "John D'Errico" <woodchips@rochester.rr.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Spline ,Histopolation or Histospline - Matlab
Date: Wed, 8 Jul 2009 19:08:02 +0000 (UTC)
Organization: John D'Errico (1-3LEW5R)
Lines: 172
Message-ID: <h32qqi$l71$1@fred.mathworks.com>
References: <h32fct$aut$1@fred.mathworks.com> <h32gfc$d2q$1@fred.mathworks.com>
Reply-To: "John D'Errico" <woodchips@rochester.rr.com>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1247080082 21729 172.30.248.35 (8 Jul 2009 19:08:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Wed, 8 Jul 2009 19:08:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 869215
Xref: news.mathworks.com comp.soft-sys.matlab:553872


"John D'Errico" <woodchips@rochester.rr.com> wrote in message <h32gfc$d2q$1@fred.mathworks.com>...
> "Paul Meissner" <magnesit@gmx.at> wrote in message <h32fct$aut$1@fred.mathworks.com>...
> > Hello,
> > 
> > I have a very tricky problem, I want to create an interpolation of a Step-Function. The problem is that the interpolation should keep the condition of area-preserving for every single step. This sounds like a histopolation, does it? Using Google (search: histopolation, histospline) had shown many theoretical abstracts but no Matlab-code. I found only a code-example for Mathematica but at all this is not working for me. (The spline is to "nervous"(rippled) for making a good interpolation)  
> > 
> > I hope you can give me an advice, how to make such an interpolation in Matlab. I would also be happy if you know a software which is able to do such an interpolation.
> > 
> > Best regards,
> > Paul M.
> 
> So instead of using the spline to interpolate specific
> data points, you wish to find a piecewise cubic (C1?
> C2?) that reproduces the area in each interval?
> 
> Should be doable. The spline is linear in its parameters.
> There will be a minimum energy configuration that
> satisfies the areal constraints. I already have such a
> constraint in my slm tools, but only for the overall
> area of the curve.
> 
> John

Very much hacked together, this uses my lse from the
file exchange. The result is in the form of a SLM
spline. If you wish a pp form, my slm tools include a
converter slm2pp.

A quick test is all I have time for. In a few days I'll
have time to clean up the code, but not now. Sorry.
I've not even written the help, or cleaned up any
error checks.


bins = 0:.1:1;
areas  = sqrt(bins(2:end));
areas = areas / sum(areas);

slm = histospline(bins,areas)
slm = 
      form: 'slm'
    degree: 3
     knots: [11x1 double]
      coef: [11x2 double]

slmpar(slm,'integral')
ans =
         0.999999999999997

plotslm(slm)

All seems reasonable. 

====================================
function slm = histospline(bins,areas,histstyle,continuity)
% generates a cubic spline that integrates to the given set of areas in each bin

if (nargin < 3) || isempty(histstyle)
  % set the default
  histstyle = 'area';
else
  % check to see that either 'height' or 'area'
  % was provided here.


end
if (nargin < 4) || isempty(continuity)
  % set the default
  continuity = 'C2';
else
  % 'C1' or 'C2' are the options.


  
end

x = bins(:);
areas = areas(:);

nk = length(x);
nc = 2*nk;
if (nk-1) ~= length(areas)
  error('HISTOSPLINE:incompatibledata', ...
    'There must be exactly one more bin boundary than area supplied')
end

% the knots are just the bin boundaries
dx = diff(x);

if strcmp(histstyle,'height')
  % if the heights were specified, then scale to form an area
  areas = areas .* dx;
end

% Regularizer
% We are integrating the piecewise linear f''(x), as
% a quadratic form in terms of the (unknown) second
% derivatives at the knots.
Mreg=zeros(nk,nk);
Mreg(1,1:2)=[dx(1)/3 , dx(1)/6];
Mreg(nk,nk+[-1 0])=[dx(end)/6 , dx(end)/3];
for i=2:(nk-1)
  Mreg(i,i+[-1 0 1])=[dx(i-1)/6 , (dx(i-1)+dx(i))/3 , dx(i)/6];
end
% do a matrix square root. cholesky is simplest & ok since regmat is
% positive definite. this way we can write the quadratic form as:
%    s'*r'*r*s,
% where s is the vector of second derivatives at the knots.
Mreg=chol(Mreg);

% next, write the second derivatives as a function of the
% function values and first derivatives.   s = [sf,sd]*[f;d]
sf = zeros(nk,nk);
sd = zeros(nk,nk);
for i = 1:(nk-1)
  sf(i,i+[0 1]) = [-1 1].*(6/(dx(i)^2));
  sd(i,i+[0 1]) = [-4 -2]./dx(i);
end
sf(nk,nk+[-1 0]) = [1 -1].*(6/(dx(end)^2));
sd(nk,nk+[-1 0]) = [2 4]./dx(end);
Mreg = Mreg*[sf,sd];
rhsreg = zeros(nk,1);

if strcmp(continuity,'C2')
  % C2 continuity across knots
  Meq = zeros(nk-2,nc);
  rhseq = zeros(nk-2,1);
  
  for i = 1:(nk-2)
    Meq(i,[i,i+1]) = [6 -6]./(dx(i).^2);
    Meq(i,[i+1,i+2]) = Meq(i,[i+1,i+2]) + [6 -6]./(dx(i+1).^2);
    
    Meq(i,nk+[i,i+1]) = [2 4]./dx(i);
    Meq(i,nk+[i+1,i+2]) = Meq(i,nk+[i+1,i+2]) + [4 2]./dx(i+1);
  end
else
  % allow the spline to be C1
  Meq = zeros(0,nc);
  rhseq = zeros(0,1);
end

% set up the integral equality constraints. We can
% be very lazy here. Simpson's rule over the support
% will be exact, evaluating at each midpoint. Don't
% forget that the effective "stepsize" is actually
% dx(i)/2
Mint = zeros(nk-1,nc);
for i = 1:(nk-1)
  Mint(i,i+[0 1]) = dx(i)/6;
  % the midpoint
  Mint(i,i+[0 1 nk nk+1]) = Mint(i,i+[0 1 nk nk+1]) + ...
    [1/2 , 1/2 , (1/8).*dx(i) , (-1/8).*dx(i)]*4*dx(i)/6;
end
Meq = [Meq;Mint];
rhseq = [rhseq;areas];

% Solve for the spline coefficients. I could have
% set this up in a quadprog form, but I hacked the
% pieces for this out of slmengine. So lse is
% perfect for the task, since it does not require
% the optimization toolbox and there are no
% inequality constraints. Use the pinv solver in lse.
solverflag = 'pinv';
coef = lse(Mreg,rhsreg,Meq,rhseq,solverflag);

% output in a slm form
slm.form = 'slm';
slm.degree = 3;
slm.knots = bins';
slm.coef = reshape(coef,nk,2);

=================================