"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 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)
> > >
> > > 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 (nk1) ~= 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:(nk1)
> Mreg(i,i+[1 0 1])=[dx(i1)/6 , (dx(i1)+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:(nk1)
> 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(nk2,nc);
> rhseq = zeros(nk2,1);
>
> for i = 1:(nk2)
> 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(nk1,nc);
> for i = 1:(nk1)
> 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 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).
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.
Thanking you in anticipation,
Paul M.
