MATLAB Newsgroup

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.

"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

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

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

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

"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

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

You can think of your watch list as threads that you have bookmarked.

You can add tags, authors, threads, and even search results to your watch list. This way you can easily keep track of topics that you're interested in. To view your watch list, click on the "My Newsreader" link.

To add items to your watch list, click the "add to watch list" link at the bottom of any page.

To add search criteria to your watch list, search for the desired term in the search box. Click on the "Add this search to my watch list" link on the search results page.

You can also add a tag to your watch list by searching for the tag with the directive "tag:tag_name" where tag_name is the name of the tag you would like to watch.

To add an author to your watch list, go to the author's profile page and click on the "Add this author to my watch list" link at the top of the page. You can also add an author to your watch list by going to a thread that the author has posted to and clicking on the "Add this author to my watch list" link. You will be notified whenever the author makes a post.

To add a thread to your watch list, go to the thread page and click the "Add this thread to my watch list" link at the top of the page.

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Got questions?

Get answers.

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi test

Learn moreDiscover what MATLAB ^{®} can do for your career.

Opportunities for recent engineering grads.

Apply TodayThe newsgroups are a worldwide forum that is open to everyone. Newsgroups are used to discuss a huge range of topics, make announcements, and trade files.

Discussions are threaded, or grouped in a way that allows you to read a posted message and all of its replies in chronological order. This makes it easy to follow the thread of the conversation, and to see what’s already been said before you post your own reply or make a new posting.

Newsgroup content is distributed by servers hosted by various organizations on the Internet. Messages are exchanged and managed using open-standard protocols. No single entity “owns” the newsgroups.

There are thousands of newsgroups, each addressing a single topic or area of interest. The MATLAB Central Newsreader posts and displays messages in the comp.soft-sys.matlab newsgroup.

**MATLAB Central**

You can use the integrated newsreader at the MATLAB Central website to read and post messages in this newsgroup. MATLAB Central is hosted by MathWorks.

Messages posted through the MATLAB Central Newsreader are seen by everyone using the newsgroups, regardless of how they access the newsgroups. There are several advantages to using MATLAB Central.

**One Account**

Your MATLAB Central account is tied to your MathWorks Account for easy access.

**Use the Email Address of Your Choice**

The MATLAB Central Newsreader allows you to define an alternative email address as your posting address, avoiding clutter in your primary mailbox and reducing spam.

**Spam Control**

Most newsgroup spam is filtered out by the MATLAB Central Newsreader.

**Tagging**

Messages can be tagged with a relevant label by any signed-in user. Tags can be used as keywords to find particular files of interest, or as a way to categorize your bookmarked postings. You may choose to allow others to view your tags, and you can view or search others’ tags as well as those of the community at large. Tagging provides a way to see both the big trends and the smaller, more obscure ideas and applications.

**Watch lists**

Setting up watch lists allows you to be notified of updates made to postings selected by author, thread, or any search variable. Your watch list notifications can be sent by email (daily digest or immediate), displayed in My Newsreader, or sent via RSS feed.

- Use a newsreader through your school, employer, or internet service provider
- Pay for newsgroup access from a commercial provider
- Use Google Groups
- Mathforum.org provides a newsreader with access to the comp.soft sys.matlab newsgroup
- Run your own server. For typical instructions, see: http://www.slyck.com/ng.php?page=2