MATLAB Central Newsreader - Spline ,Histopolation or Histospline - Matlab

Spline ,Histopolation or Histospline - Matlab

Paul Meissner — Wed, 08 Jul 2009 15:53:01 +0000

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.

Re: Spline ,Histopolation or Histospline - Matlab

John D'Errico — Wed, 08 Jul 2009 16:11:24 +0000

"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

Re: Spline ,Histopolation or Histospline - Matlab

John D'Errico — Wed, 08 Jul 2009 19:08:02 +0000

"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);

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

Re: Spline ,Histopolation or Histospline - Matlab

Paul Meissner — Thu, 09 Jul 2009 11:13:02 +0000

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <h32qqi$l71$1@fred.mathworks.com>...
> "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);
>
> =================================

Hallo John,

I am very happy that you answered so quick! Thank you for your code example! I had downloaded your SLM Toolpak und your LSE function. Nice work!

I have uploaded two examples (h**p://yfrog.com/09ex1ktljx), both based on the same data. By viewing the examples you will see my problem, the first picture is made by your code-example, works fine, but is not following the shape. The second picture shows a spline which is hand-drawn by me ;=) (sorry it is not exact and the last section spline is really bad, but it should be enough for illustrating the problem).

Maybe you are interested why I need such a spline? I study Mineral-Processing and therefore I have to analyse some data for my master-thesis This data represents the quality of separating mineral grains by there characteristics. An old DIN standard exists of creating such a diagram by smoothing the step-function with a spline (area-preserving). In practise it is necessary to draw the spline by hand, there is no automated standard-conform solution available. Everyone is drawing the splines in an other way and so you cannot reproduce the analysis again.

Thanking you in anticipation,

Paul M.

Re: Spline ,Histopolation or Histospline - Matlab

John D'Errico — Thu, 09 Jul 2009 12:44:02 +0000

"Paul Meissner" <magnesit@gmx.at> wrote in message <h34jbu$27q$1@fred.mathworks.com>...

> Hallo John,
>
> I am very happy that you answered so quick! Thank you for your code example! I had downloaded your SLM Toolpak und your LSE function. Nice work!
>
> I have uploaded two examples (h**p://yfrog.com/09ex1ktljx), both based on the same data. By viewing the examples you will see my problem, the first picture is made by your code-example, works fine, but is not following the shape. The second picture shows a spline which is hand-drawn by me ;=) (sorry it is not exact and the last section spline is really bad, but it should be enough for illustrating the problem).
>
> Maybe you are interested why I need such a spline? I study Mineral-Processing and therefore I have to analyse some data for my master-thesis This data represents the quality of separating mineral grains by there characteristics. An old DIN standard exists of creating such a diagram by smoothing the step-function with a spline (area-preserving). In practise it is necessary to draw the spline by hand, there is no automated standard-conform solution available. Everyone is drawing the splines in an other way and so you cannot reproduce the analysis again.
>

Hi Paul,

To be honest, I'd suggest that the best solution is not
to use the code I posted in this thread. It was an
interesting exercise, but not of tremendous value, as
it is not that intelligent. Better is to use a tool that can
employ your knowledge of the system - my SLM tools.
I don't have an explicit area preserving constraint built
in there, although this is an idea that I'll follow up on if
I can think of the proper way to provide that information
to slmengine.

Using the slm tools, I might specify the midpoints of the
bars from your histogram as your data points. Then use
a monotonicity constraint on the curve, and force it to
be everywhere positive.

It appears from your figures that you are working with a
cumulative histogram. So the call with slmengine might
look like this:

S = slmengine(x,y,'plot','on','knots',10,'incr','on', ...
'leftvalue',0,'rightvalue',100);

- Thus, I've told it to generate a plot with the data and
the resulting curve overlaid on top.

- I set the number of knots to be 10, equally spaced along
the curve by default. You may choose more or fewer knots
as you find necessary, but the curve shape itself will often
be more controlled by an intelligent set of constraints
applied by the user.

- The curve will be a monotone increasing function.

- It will be everywhere positive, since it is an increasing
function that is zero at the left hand end point of the
domain of the spline.

- Since this is apparently a cumulative distribution, the
right end point is set to 100%.

You should find that the curve drawn, while it is not
forced to have the exact areal constraints for each bin
of the histogram, will be quite reasonable. If not, then
only a moderate amount of tweaking should allow a
reasonable fit. For example, you might even choose a
few more knots, then add an extra property to the list,
like this:

S = slmengine(x,y,'plot','on','knots',15,'incr','on', ...
'leftvalue',0,'rightvalue',100,'regularization','cross');

As I said above, if I can think of the proper format
to input data oints to be interpreted as a histogram, I'll
see if I might implement explicit areal constraints in a
future release. Until then, I hope the prescriptive rules
above might be acceptable for your task.

John

Re: Spline ,Histopolation or Histospline - Matlab

Irma Hernandez — Fri, 30 Oct 2009 00:35:05 +0000

Dear John,

I am not sure if my question is related to the one posted by Paul. I am interested in the case where we think in 2d and data is grouped in regions (bins) I want to fit a 2d spline where the conservation of volume is happening in each region. It is sometimes called "thin plate histospline".

Do you know if there is a matlab code related to this problem? Thanks

Irma