http://www.mathworks.com/matlabcentral/newsreader/view_thread/255623
MATLAB Central Newsreader  Spline ,Histopolation or Histospline  Matlab
Feed for thread: Spline ,Histopolation or Histospline  Matlab
enus
©19942015 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Wed, 08 Jul 2009 15:53:01 +0000
Spline ,Histopolation or Histospline  Matlab
http://www.mathworks.com/matlabcentral/newsreader/view_thread/255623#663625
Paul Meissner
Hello,<br>
<br>
I have a very tricky problem, I want to create an interpolation of a StepFunction. The problem is that the interpolation should keep the condition of areapreserving for every single step. This sounds like a histopolation, does it? Using Google (search: histopolation, histospline) had shown many theoretical abstracts but no Matlabcode. I found only a codeexample for Mathematica but at all this is not working for me. (The spline is to "nervous"(rippled) for making a good interpolation) <br>
<br>
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.<br>
<br>
Best regards,<br>
Paul M.

Wed, 08 Jul 2009 16:11:24 +0000
Re: Spline ,Histopolation or Histospline  Matlab
http://www.mathworks.com/matlabcentral/newsreader/view_thread/255623#663631
John D'Errico
"Paul Meissner" <magnesit@gmx.at> wrote in message <h32fct$aut$1@fred.mathworks.com>...<br>
> Hello,<br>
> <br>
> I have a very tricky problem, I want to create an interpolation of a StepFunction. The problem is that the interpolation should keep the condition of areapreserving for every single step. This sounds like a histopolation, does it? Using Google (search: histopolation, histospline) had shown many theoretical abstracts but no Matlabcode. I found only a codeexample for Mathematica but at all this is not working for me. (The spline is to "nervous"(rippled) for making a good interpolation) <br>
> <br>
> 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.<br>
> <br>
> Best regards,<br>
> Paul M.<br>
<br>
So instead of using the spline to interpolate specific<br>
data points, you wish to find a piecewise cubic (C1?<br>
C2?) that reproduces the area in each interval?<br>
<br>
Should be doable. The spline is linear in its parameters.<br>
There will be a minimum energy configuration that<br>
satisfies the areal constraints. I already have such a<br>
constraint in my slm tools, but only for the overall<br>
area of the curve.<br>
<br>
John

Wed, 08 Jul 2009 19:08:02 +0000
Re: Spline ,Histopolation or Histospline  Matlab
http://www.mathworks.com/matlabcentral/newsreader/view_thread/255623#663691
John D'Errico
"John D'Errico" <woodchips@rochester.rr.com> wrote in message <h32gfc$d2q$1@fred.mathworks.com>...<br>
> "Paul Meissner" <magnesit@gmx.at> wrote in message <h32fct$aut$1@fred.mathworks.com>...<br>
> > Hello,<br>
> > <br>
> > I have a very tricky problem, I want to create an interpolation of a StepFunction. The problem is that the interpolation should keep the condition of areapreserving for every single step. This sounds like a histopolation, does it? Using Google (search: histopolation, histospline) had shown many theoretical abstracts but no Matlabcode. I found only a codeexample for Mathematica but at all this is not working for me. (The spline is to "nervous"(rippled) for making a good interpolation) <br>
> > <br>
> > 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.<br>
> > <br>
> > Best regards,<br>
> > Paul M.<br>
> <br>
> So instead of using the spline to interpolate specific<br>
> data points, you wish to find a piecewise cubic (C1?<br>
> C2?) that reproduces the area in each interval?<br>
> <br>
> Should be doable. The spline is linear in its parameters.<br>
> There will be a minimum energy configuration that<br>
> satisfies the areal constraints. I already have such a<br>
> constraint in my slm tools, but only for the overall<br>
> area of the curve.<br>
> <br>
> John<br>
<br>
Very much hacked together, this uses my lse from the<br>
file exchange. The result is in the form of a SLM<br>
spline. If you wish a pp form, my slm tools include a<br>
converter slm2pp.<br>
<br>
A quick test is all I have time for. In a few days I'll<br>
have time to clean up the code, but not now. Sorry.<br>
I've not even written the help, or cleaned up any<br>
error checks.<br>
<br>
<br>
bins = 0:.1:1;<br>
areas = sqrt(bins(2:end));<br>
areas = areas / sum(areas);<br>
<br>
slm = histospline(bins,areas)<br>
slm = <br>
form: 'slm'<br>
degree: 3<br>
knots: [11x1 double]<br>
coef: [11x2 double]<br>
<br>
slmpar(slm,'integral')<br>
ans =<br>
0.999999999999997<br>
<br>
plotslm(slm)<br>
<br>
All seems reasonable. <br>
<br>
====================================<br>
function slm = histospline(bins,areas,histstyle,continuity)<br>
% generates a cubic spline that integrates to the given set of areas in each bin<br>
<br>
if (nargin < 3)  isempty(histstyle)<br>
% set the default<br>
histstyle = 'area';<br>
else<br>
% check to see that either 'height' or 'area'<br>
% was provided here.<br>
<br>
<br>
end<br>
if (nargin < 4)  isempty(continuity)<br>
% set the default<br>
continuity = 'C2';<br>
else<br>
% 'C1' or 'C2' are the options.<br>
<br>
<br>
<br>
end<br>
<br>
x = bins(:);<br>
areas = areas(:);<br>
<br>
nk = length(x);<br>
nc = 2*nk;<br>
if (nk1) ~= length(areas)<br>
error('HISTOSPLINE:incompatibledata', ...<br>
'There must be exactly one more bin boundary than area supplied')<br>
end<br>
<br>
% the knots are just the bin boundaries<br>
dx = diff(x);<br>
<br>
if strcmp(histstyle,'height')<br>
% if the heights were specified, then scale to form an area<br>
areas = areas .* dx;<br>
end<br>
<br>
% Regularizer<br>
% We are integrating the piecewise linear f''(x), as<br>
% a quadratic form in terms of the (unknown) second<br>
% derivatives at the knots.<br>
Mreg=zeros(nk,nk);<br>
Mreg(1,1:2)=[dx(1)/3 , dx(1)/6];<br>
Mreg(nk,nk+[1 0])=[dx(end)/6 , dx(end)/3];<br>
for i=2:(nk1)<br>
Mreg(i,i+[1 0 1])=[dx(i1)/6 , (dx(i1)+dx(i))/3 , dx(i)/6];<br>
end<br>
% do a matrix square root. cholesky is simplest & ok since regmat is<br>
% positive definite. this way we can write the quadratic form as:<br>
% s'*r'*r*s,<br>
% where s is the vector of second derivatives at the knots.<br>
Mreg=chol(Mreg);<br>
<br>
% next, write the second derivatives as a function of the<br>
% function values and first derivatives. s = [sf,sd]*[f;d]<br>
sf = zeros(nk,nk);<br>
sd = zeros(nk,nk);<br>
for i = 1:(nk1)<br>
sf(i,i+[0 1]) = [1 1].*(6/(dx(i)^2));<br>
sd(i,i+[0 1]) = [4 2]./dx(i);<br>
end<br>
sf(nk,nk+[1 0]) = [1 1].*(6/(dx(end)^2));<br>
sd(nk,nk+[1 0]) = [2 4]./dx(end);<br>
Mreg = Mreg*[sf,sd];<br>
rhsreg = zeros(nk,1);<br>
<br>
if strcmp(continuity,'C2')<br>
% C2 continuity across knots<br>
Meq = zeros(nk2,nc);<br>
rhseq = zeros(nk2,1);<br>
<br>
for i = 1:(nk2)<br>
Meq(i,[i,i+1]) = [6 6]./(dx(i).^2);<br>
Meq(i,[i+1,i+2]) = Meq(i,[i+1,i+2]) + [6 6]./(dx(i+1).^2);<br>
<br>
Meq(i,nk+[i,i+1]) = [2 4]./dx(i);<br>
Meq(i,nk+[i+1,i+2]) = Meq(i,nk+[i+1,i+2]) + [4 2]./dx(i+1);<br>
end<br>
else<br>
% allow the spline to be C1<br>
Meq = zeros(0,nc);<br>
rhseq = zeros(0,1);<br>
end<br>
<br>
% set up the integral equality constraints. We can<br>
% be very lazy here. Simpson's rule over the support<br>
% will be exact, evaluating at each midpoint. Don't<br>
% forget that the effective "stepsize" is actually<br>
% dx(i)/2<br>
Mint = zeros(nk1,nc);<br>
for i = 1:(nk1)<br>
Mint(i,i+[0 1]) = dx(i)/6;<br>
% the midpoint<br>
Mint(i,i+[0 1 nk nk+1]) = Mint(i,i+[0 1 nk nk+1]) + ...<br>
[1/2 , 1/2 , (1/8).*dx(i) , (1/8).*dx(i)]*4*dx(i)/6;<br>
end<br>
Meq = [Meq;Mint];<br>
rhseq = [rhseq;areas];<br>
<br>
% Solve for the spline coefficients. I could have<br>
% set this up in a quadprog form, but I hacked the<br>
% pieces for this out of slmengine. So lse is<br>
% perfect for the task, since it does not require<br>
% the optimization toolbox and there are no<br>
% inequality constraints. Use the pinv solver in lse.<br>
solverflag = 'pinv';<br>
coef = lse(Mreg,rhsreg,Meq,rhseq,solverflag);<br>
<br>
% output in a slm form<br>
slm.form = 'slm';<br>
slm.degree = 3;<br>
slm.knots = bins';<br>
slm.coef = reshape(coef,nk,2);<br>
<br>
=================================

Thu, 09 Jul 2009 11:13:02 +0000
Re: Spline ,Histopolation or Histospline  Matlab
http://www.mathworks.com/matlabcentral/newsreader/view_thread/255623#663897
Paul Meissner
"John D'Errico" <woodchips@rochester.rr.com> wrote in message <h32qqi$l71$1@fred.mathworks.com>...<br>
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <h32gfc$d2q$1@fred.mathworks.com>...<br>
> > "Paul Meissner" <magnesit@gmx.at> wrote in message <h32fct$aut$1@fred.mathworks.com>...<br>
> > > Hello,<br>
> > > <br>
> > > I have a very tricky problem, I want to create an interpolation of a StepFunction. The problem is that the interpolation should keep the condition of areapreserving for every single step. This sounds like a histopolation, does it? Using Google (search: histopolation, histospline) had shown many theoretical abstracts but no Matlabcode. I found only a codeexample for Mathematica but at all this is not working for me. (The spline is to "nervous"(rippled) for making a good interpolation) <br>
> > > <br>
> > > 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.<br>
> > > <br>
> > > Best regards,<br>
> > > Paul M.<br>
> > <br>
> > So instead of using the spline to interpolate specific<br>
> > data points, you wish to find a piecewise cubic (C1?<br>
> > C2?) that reproduces the area in each interval?<br>
> > <br>
> > Should be doable. The spline is linear in its parameters.<br>
> > There will be a minimum energy configuration that<br>
> > satisfies the areal constraints. I already have such a<br>
> > constraint in my slm tools, but only for the overall<br>
> > area of the curve.<br>
> > <br>
> > John<br>
> <br>
> Very much hacked together, this uses my lse from the<br>
> file exchange. The result is in the form of a SLM<br>
> spline. If you wish a pp form, my slm tools include a<br>
> converter slm2pp.<br>
> <br>
> A quick test is all I have time for. In a few days I'll<br>
> have time to clean up the code, but not now. Sorry.<br>
> I've not even written the help, or cleaned up any<br>
> error checks.<br>
> <br>
> <br>
> bins = 0:.1:1;<br>
> areas = sqrt(bins(2:end));<br>
> areas = areas / sum(areas);<br>
> <br>
> slm = histospline(bins,areas)<br>
> slm = <br>
> form: 'slm'<br>
> degree: 3<br>
> knots: [11x1 double]<br>
> coef: [11x2 double]<br>
> <br>
> slmpar(slm,'integral')<br>
> ans =<br>
> 0.999999999999997<br>
> <br>
> plotslm(slm)<br>
> <br>
> All seems reasonable. <br>
> <br>
> ====================================<br>
> function slm = histospline(bins,areas,histstyle,continuity)<br>
> % generates a cubic spline that integrates to the given set of areas in each bin<br>
> <br>
> if (nargin < 3)  isempty(histstyle)<br>
> % set the default<br>
> histstyle = 'area';<br>
> else<br>
> % check to see that either 'height' or 'area'<br>
> % was provided here.<br>
> <br>
> <br>
> end<br>
> if (nargin < 4)  isempty(continuity)<br>
> % set the default<br>
> continuity = 'C2';<br>
> else<br>
> % 'C1' or 'C2' are the options.<br>
> <br>
> <br>
> <br>
> end<br>
> <br>
> x = bins(:);<br>
> areas = areas(:);<br>
> <br>
> nk = length(x);<br>
> nc = 2*nk;<br>
> if (nk1) ~= length(areas)<br>
> error('HISTOSPLINE:incompatibledata', ...<br>
> 'There must be exactly one more bin boundary than area supplied')<br>
> end<br>
> <br>
> % the knots are just the bin boundaries<br>
> dx = diff(x);<br>
> <br>
> if strcmp(histstyle,'height')<br>
> % if the heights were specified, then scale to form an area<br>
> areas = areas .* dx;<br>
> end<br>
> <br>
> % Regularizer<br>
> % We are integrating the piecewise linear f''(x), as<br>
> % a quadratic form in terms of the (unknown) second<br>
> % derivatives at the knots.<br>
> Mreg=zeros(nk,nk);<br>
> Mreg(1,1:2)=[dx(1)/3 , dx(1)/6];<br>
> Mreg(nk,nk+[1 0])=[dx(end)/6 , dx(end)/3];<br>
> for i=2:(nk1)<br>
> Mreg(i,i+[1 0 1])=[dx(i1)/6 , (dx(i1)+dx(i))/3 , dx(i)/6];<br>
> end<br>
> % do a matrix square root. cholesky is simplest & ok since regmat is<br>
> % positive definite. this way we can write the quadratic form as:<br>
> % s'*r'*r*s,<br>
> % where s is the vector of second derivatives at the knots.<br>
> Mreg=chol(Mreg);<br>
> <br>
> % next, write the second derivatives as a function of the<br>
> % function values and first derivatives. s = [sf,sd]*[f;d]<br>
> sf = zeros(nk,nk);<br>
> sd = zeros(nk,nk);<br>
> for i = 1:(nk1)<br>
> sf(i,i+[0 1]) = [1 1].*(6/(dx(i)^2));<br>
> sd(i,i+[0 1]) = [4 2]./dx(i);<br>
> end<br>
> sf(nk,nk+[1 0]) = [1 1].*(6/(dx(end)^2));<br>
> sd(nk,nk+[1 0]) = [2 4]./dx(end);<br>
> Mreg = Mreg*[sf,sd];<br>
> rhsreg = zeros(nk,1);<br>
> <br>
> if strcmp(continuity,'C2')<br>
> % C2 continuity across knots<br>
> Meq = zeros(nk2,nc);<br>
> rhseq = zeros(nk2,1);<br>
> <br>
> for i = 1:(nk2)<br>
> Meq(i,[i,i+1]) = [6 6]./(dx(i).^2);<br>
> Meq(i,[i+1,i+2]) = Meq(i,[i+1,i+2]) + [6 6]./(dx(i+1).^2);<br>
> <br>
> Meq(i,nk+[i,i+1]) = [2 4]./dx(i);<br>
> Meq(i,nk+[i+1,i+2]) = Meq(i,nk+[i+1,i+2]) + [4 2]./dx(i+1);<br>
> end<br>
> else<br>
> % allow the spline to be C1<br>
> Meq = zeros(0,nc);<br>
> rhseq = zeros(0,1);<br>
> end<br>
> <br>
> % set up the integral equality constraints. We can<br>
> % be very lazy here. Simpson's rule over the support<br>
> % will be exact, evaluating at each midpoint. Don't<br>
> % forget that the effective "stepsize" is actually<br>
> % dx(i)/2<br>
> Mint = zeros(nk1,nc);<br>
> for i = 1:(nk1)<br>
> Mint(i,i+[0 1]) = dx(i)/6;<br>
> % the midpoint<br>
> Mint(i,i+[0 1 nk nk+1]) = Mint(i,i+[0 1 nk nk+1]) + ...<br>
> [1/2 , 1/2 , (1/8).*dx(i) , (1/8).*dx(i)]*4*dx(i)/6;<br>
> end<br>
> Meq = [Meq;Mint];<br>
> rhseq = [rhseq;areas];<br>
> <br>
> % Solve for the spline coefficients. I could have<br>
> % set this up in a quadprog form, but I hacked the<br>
> % pieces for this out of slmengine. So lse is<br>
> % perfect for the task, since it does not require<br>
> % the optimization toolbox and there are no<br>
> % inequality constraints. Use the pinv solver in lse.<br>
> solverflag = 'pinv';<br>
> coef = lse(Mreg,rhsreg,Meq,rhseq,solverflag);<br>
> <br>
> % output in a slm form<br>
> slm.form = 'slm';<br>
> slm.degree = 3;<br>
> slm.knots = bins';<br>
> slm.coef = reshape(coef,nk,2);<br>
> <br>
> =================================<br>
<br>
Hallo John,<br>
<br>
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! <br>
<br>
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 codeexample, works fine, but is not following the shape. The second picture shows a spline which is handdrawn by me ;=) (sorry it is not exact and the last section spline is really bad, but it should be enough for illustrating the problem). <br>
<br>
Maybe you are interested why I need such a spline? I study MineralProcessing and therefore I have to analyse some data for my masterthesis 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 stepfunction with a spline (areapreserving). In practise it is necessary to draw the spline by hand, there is no automated standardconform solution available. Everyone is drawing the splines in an other way and so you cannot reproduce the analysis again.<br>
<br>
Thanking you in anticipation,<br>
<br>
Paul M.<br>

Thu, 09 Jul 2009 12:44:02 +0000
Re: Spline ,Histopolation or Histospline  Matlab
http://www.mathworks.com/matlabcentral/newsreader/view_thread/255623#663924
John D'Errico
"Paul Meissner" <magnesit@gmx.at> wrote in message <h34jbu$27q$1@fred.mathworks.com>...<br>
<br>
> Hallo John,<br>
> <br>
> 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! <br>
> <br>
> 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 codeexample, works fine, but is not following the shape. The second picture shows a spline which is handdrawn by me ;=) (sorry it is not exact and the last section spline is really bad, but it should be enough for illustrating the problem). <br>
> <br>
> Maybe you are interested why I need such a spline? I study MineralProcessing and therefore I have to analyse some data for my masterthesis 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 stepfunction with a spline (areapreserving). In practise it is necessary to draw the spline by hand, there is no automated standardconform solution available. Everyone is drawing the splines in an other way and so you cannot reproduce the analysis again.<br>
> <br>
<br>
Hi Paul,<br>
<br>
To be honest, I'd suggest that the best solution is not<br>
to use the code I posted in this thread. It was an<br>
interesting exercise, but not of tremendous value, as<br>
it is not that intelligent. Better is to use a tool that can<br>
employ your knowledge of the system  my SLM tools.<br>
I don't have an explicit area preserving constraint built<br>
in there, although this is an idea that I'll follow up on if<br>
I can think of the proper way to provide that information<br>
to slmengine.<br>
<br>
Using the slm tools, I might specify the midpoints of the<br>
bars from your histogram as your data points. Then use<br>
a monotonicity constraint on the curve, and force it to<br>
be everywhere positive.<br>
<br>
It appears from your figures that you are working with a<br>
cumulative histogram. So the call with slmengine might<br>
look like this:<br>
<br>
S = slmengine(x,y,'plot','on','knots',10,'incr','on', ...<br>
'leftvalue',0,'rightvalue',100);<br>
<br>
 Thus, I've told it to generate a plot with the data and<br>
the resulting curve overlaid on top.<br>
<br>
 I set the number of knots to be 10, equally spaced along<br>
the curve by default. You may choose more or fewer knots<br>
as you find necessary, but the curve shape itself will often<br>
be more controlled by an intelligent set of constraints<br>
applied by the user.<br>
<br>
 The curve will be a monotone increasing function.<br>
<br>
 It will be everywhere positive, since it is an increasing<br>
function that is zero at the left hand end point of the<br>
domain of the spline.<br>
<br>
 Since this is apparently a cumulative distribution, the<br>
right end point is set to 100%.<br>
<br>
You should find that the curve drawn, while it is not<br>
forced to have the exact areal constraints for each bin<br>
of the histogram, will be quite reasonable. If not, then<br>
only a moderate amount of tweaking should allow a<br>
reasonable fit. For example, you might even choose a<br>
few more knots, then add an extra property to the list,<br>
like this:<br>
<br>
S = slmengine(x,y,'plot','on','knots',15,'incr','on', ...<br>
'leftvalue',0,'rightvalue',100,'regularization','cross');<br>
<br>
As I said above, if I can think of the proper format<br>
to input data oints to be interpreted as a histogram, I'll<br>
see if I might implement explicit areal constraints in a<br>
future release. Until then, I hope the prescriptive rules<br>
above might be acceptable for your task.<br>
<br>
John

Fri, 30 Oct 2009 00:35:05 +0000
Re: Spline ,Histopolation or Histospline  Matlab
http://www.mathworks.com/matlabcentral/newsreader/view_thread/255623#690860
Irma Hernandez
Dear John,<br>
<br>
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".<br>
<br>
Do you know if there is a matlab code related to this problem? Thanks<br>
<br>
Irma