Path: news.mathworks.com!not-for-mail
From: "Paul Meissner" <magnesit@gmx.at>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Spline ,Histopolation or Histospline - Matlab
Date: Thu, 9 Jul 2009 11:13:02 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 186
Message-ID: <h34jbu$27q$1@fred.mathworks.com>
References: <h32fct$aut$1@fred.mathworks.com> <h32gfc$d2q$1@fred.mathworks.com> <h32qqi$l71$1@fred.mathworks.com>
Reply-To: "Paul Meissner" <magnesit@gmx.at>
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 1247137982 2298 172.30.248.35 (9 Jul 2009 11:13:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Thu, 9 Jul 2009 11:13:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1900131
Xref: news.mathworks.com comp.soft-sys.matlab:554078


"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.